\def\ie{i.~e.}
\def\domega{\partial \Omega}
\def\domegai{\partial \Omega_i}
\def\domegao{\partial \Omega_o}
\def\scal#1{\left< #1 \right>}
\def\ts{t_S}
\def\td{t_D}
\def\ex#1#2#3#4{\chi(#1,#2,#3,#4)}
\def\dep#1{\left(#1\right)}
\def\depb#1{\left[#1\right]}
\def\depc#1{\left\{#1\right\}}
\def\dem{1/2}
\def\eq#1{Eq.~\ref{eq:#1}}
\def\sec#1{Section~\ref{sec:#1}}
\def\dq{{\dep{\delta c}}}
\def\grad{{\mbox{\bf grad}}}
\def\div#1{{\mbox{div}}\dep{#1}}
\def\kdiff#1{\frac{1}{\rho}\frac{\partial}{\partial z}
\dep{\rho K_z \frac{\partial #1}{\partial z}}}
\def\deriv#1#2{\frac{\partial #1}{\partial #2}}


\def\x{{\bf x}}
\def\V{{\bf v}}
\def\n{{\bf n}}
\def\dx{{\bf dx}}
\def\dxdt{\dx dt}

\def\dt#1{\frac{\partial #1}{\partial t}}



\def\local{retro}
\def\retrodir{\local/FIGURES}

\chapter{Inversion du transport atmosphérique\label{ch:retro}}

\section{Contexte}

Le troisième volet de ce travail concerne l'inversion du transport
atmosphérique. Il a été initié suite à une demande du CEA visant
à évaluer l'efficacité du réseau global de mesure de la radioactivité
atmosphérique déployé dans le cadre
du Traité d'Interdiction Complète des Essais nucléaires (TICE).
La surveillance des essais reposera sur quatre réseaux complémentaires.
Les trois premiers mesureront les ondes sismiques, hydroacoustiques et
infra-sonores et devraient permettre de localiser relativement bien
les explosions
pour des essais effectués respectivement sous terre, dans la mer ou dans
l'atmosphère. En parallèle de ces trois technologies, 
un réseau global de 80 stations, en cours de déploiement,
mesurera en permanence la concentration en
radio-éléments dans l'atmosphère permettant de mieux caractériser le type
d'évènement détecté. Toutes ces stations détecterons les aérosols
radioactifs. Un sous-réseau de 40 stations détectera également les isotopes
du xénon dont on pense qu'il sont relâchés en quantité
significative dans l'atmosphère même lors d'essais souterrains ou sous-marins.

La question qui nous était posée par le CEA était l'évaluation de la capacité de
détection de ce réseau.
A l'époque (été 1997), nous disposions d'une première version de LMDZT qui 
permettait d'effectuer des calcul de dispersion de polluants ou de
radio-éléments atmosphériques. Pour répondre à la question, une approche
directe aurait consisté à simuler des essais nucléaires
en injectant un traceur en chaque point d'un maillage de la sphère (en
chaque point de grille du modèle de transport par exemple) et de
comptabiliser les fois où la concentration simulée aux différentes
stations du réseau excédait le seuil de détection. On voit vite qu'une
telle approche conduit à des coûts informatiques prohibitifs.
Avec Robert Sadourny (LMD),
nous nous sommes convaincus que ce problème
pouvait être abordé en inversant le sens du temps dans notre modèle de
transport Eulérien,
le calcul du transport étant effectué alors en remontant à rebours
dans le temps le long des trajectoires atmosphériques.


L'idée d'utiliser des {\em modèles de récepteurs} pour ce type
de problème n'est pas nouvelle.
De façon générale, l'identification des sources pour un traceur atmosphérique
(localisation spatiale et temporelle des sources ainsi que la quantification
de la quantité de traceur émise) est une question très importante pour
beaucoup d'aspects des sciences de l'environnement.
La caractérisation  des puits et sources naturels et anthropiques
de CO$_2$ \cite[]{Kami:99b,Rayn:99,Bous:00,Gurn:02}
est par exemple un sujet très sensible dans la perspective du contrôle
des émissions de gaz à effet de serre.
La surveillance de possibles rejets accidentels de pollution par des
centrales nucléaires ou des installations chimiques est un autre bon
exemple. Les inventaires d'émissions d'espèces chimiques qui alimentent
les modèles de prévision de la pollution atmosphérique doivent aussi souvent
être en partie construits en utilisant des méthodes inverses, à partir de la
comparaison des prévisions du modèle et des observations 
\cite[voir e.~g.][]{Menu:00}.
Autre exemple encore, l'interprétation fine en termes de paléoclimats
des carottages effectués sur les calottes de glace en Antarctique ou au
Gro\"enland, nécessite qu'on soit capable de remonter à l'origine de l'eau
qui a précipité sur les calottes ou des molécules gazeuses et aérosols
qui s'y sont sont accumulés.

Nous nous intéressons donc dans ce chapitre à 
{\em  l'inversion du transport atmosphérique}, c'est à dire aux approches
permettant de remonter
à l'identification des sources à partir de mesures de concentration
dans l'atmosphère.
A noter que l'inversion du transport, au-delà des nombreuses applications
mentionnées ci-dessus, est un outil très intéressant pour analyser la 
dynamique et la chimie de l'atmosphère.

Une première approche, fréquemment utilisée pour interpréter des mesures
de concentration dans l'atmosphère, consiste à remonter à rebours
dans le temps le long des trajectoires des masses d'air.
On parle alors de {\em rétro-transport Lagrangien} ou de
{\em rétro-trajectoires}.
Cette technique permet d'obtenir facilement
une description qualitative de l'origine de la masse d'air
échantillonnée lors d'une mesure de concentration
\cite[]{Hess:96,Merr:94,Ramo:96,Chia:97,Vero:92}.
Comme dans les modèles lagrangiens directs, la composante turbulente
du transport peut être prise en compte au moyen de perturbations
aléatoires des trajectoires \cite[]{Fles:95,Vaut:01,Sieb:03}.
Cependant, les calculs de rétro-trajectoires sont souvent limités à
la composante grande échelle.

Curieusement, alors que de nombreux modèles directs sont basés sur une
description eulérienne du transport, le rétro-transport est généralement
réservé à une description lagrangienne.
On montre ici qu'on peut de façon générale, comme dans le monde direct, 
définir le {\em rétro-transport} le long des trajectoires
d'air parcourues à rebours dans le temps, et lui appliquer 
les outils et approches utilisées habituellement pour le transport direct~:
formalismes Lagrangien ou Eulérien, décomposition de Reynolds entre transport
grande échelle et turbulent, schémas numériques sophistiqués pour l'advection
(monotones, conservatifs, peu diffusifs, etc.).
L'équation du rétro-transport se déduit de l'équation directe par des
transformations simples. La diffusion turbulente, 
basée sur l'idée de mélange par des fluctuations symétriques du champs de vent,
produit par exemple la même diffusion pour le transport
direct et à rebours.
Pour les paramétrisations en flux de masse de la convection, il faut
en revanche inverser la direction des flux de masse dans les différents
compartiments.
C'est sur ces idées physiques que nous avons
développé une version rétro-transport du modèle de transport 
LMDZT \cite[]{Hour:99CRAS}  dans sa version débranchée (il s'agit d'inverser le transport de
traceurs à météo connue).

Le rétro-transport, ainsi défini,
peut être utilisé pour des inversions quantitatives
en tirant avantage d'une {\em symétrie temporelle du transport}.
En absence d'autres sources ou puits, la concentration 
d'un traceur mesurée par un détecteur à un instant $\td$, faisant
suite à l'injection d'une certaine quantité de traceur à un point
source à un instant antérieur $\ts$, peut être calculée également comme
la concentration au niveau de la source à l'instant $\ts$ 
d'un rétro-traceur émis en même quantité au lieu de détection à l'instant
$\td$.
Ce rétro-traceur correspond (à un facteur multiplicatif près) à la distribution
de l'air qui va être échantillonné ultérieurement au niveau du détecteur.
A noter qu'on ne viole pas ici l'irréversibilité du transport atmosphérique. 
Une plus forte dispersion ou diffusion dans le monde direct, correspondant
à un gain d'entropie, se traduit également par une origine de l'air
échantillonnée plus diffuse et donc également par un gain d'entropie pour 
le rétro-transport.
Le gain d'entropie ou la perte d'information est en fait identique dans les
deux cas.

Une seconde approche utilisée classiquement pour l'inversion du transport
atmosphérique consiste à partir
d'un modèle existant du transport direct, aussi complexe soit-il,
et à lui appliquer  les techniques classiques de l'assimilation
des observations telles que des variantes de l'interpolation optimale,
ou l'assimilation variationnelle \cite[voir par exemple][]{Petr:02},
les filtres ou lissages de Kalman \cite[]{Haas:96,Zhan:99}.
En particulier, les 
techniques adjointes, qui remontent aux travaux de 
Lions et Marchuk \cite[se reporter par exemple à][]{Lion:71,Marc:74,Marc:82},
fournissent une méthode systématique et puissante pour déterminer les sensibilités
d'un modèle par rapport à ses variables d'état ou à des paramètres de
contrôle.
Ces méthodes sont généralement introduite sur des bases purement 
mathématiques sans référence à une quelconque signification physique.

Les techniques
adjointes sont utilisées pour de nombreuses applications tant
météorologiques qu'océanographiques, en particulier pour l'assimilation
variationnelle d'observations \cite[]{Pene:76,LeDi:86,Tala:87,Cour:87}.
Les techniques adjointes ont également été utilisées dans de nombreuses
études relatives à l'inversion du transport atmosphérique
\cite[]{Ulia:91,Pudy:98,Kami:99a,Kami:99b,Vuki:00,Robe:91,Houw:99}.
Si l'observable considérée est une mesure de concentration
à un instant donné et à une station particulière,
l'équation de transport adjointe est en fait un modèle "orienté récepteur"
de la mesure, restreint
aux processus qui vont effectivement
influencer cette observation particulière. Le calcul adjoint détermine
la sensibilité ou encore la {\em fonction d'influence}, qui, une fois
combinée avec les sources et la concentration initiale du traceur, permet
de calculer effectivement l'équivalent de la mesure.

A partir de la propriété de symétrie du transport,
on montre facilement  que l'équation 
de rétro-transport, définie plus haut sur des considérations physiques,
est en fait l'adjoint de l'équation du transport direct pour un 
produit scalaire particulier, le {\em produit scalaire pondéré par
l'air}
\begin{equation}\label{eq:scal1}
\scal{\phi,\psi}=\int \rho \phi \psi \dxdt
\end{equation}
où $\phi$ et $\psi$ sont des concentrations massiques de traceurs.
Dans ce cas particulier, l'adjoint d'un modèle peut
donc être obtenu sur des considérations purement physiques, sans faire appel
au techniques algébriques utilisées classiquement, et qui consistent
en particulier à enchaîner des intégrations par partie.
De plus, pour ce produit scalaire particulier, le modèle adjoint
se déduit du modèle direct en changeant simplement le signe
de certains termes ce qui évite en pratique certaines difficultés relatives au
développement et à la maintenance des codes adjoints.

Si on utilise un produit scalaire plus classique de la forme
$[\phi,\psi]=\int\phi\psi \dxdt$, la symétrie est mise à mal et les équations
directes et adjointes prennent des formes différentes
 \cite[cf. e.~g.][]{Pudy:98,Vuki:00}.
\cite{Ulia:91} étaient en fait déjà arrivés sur un jeu d'équations symétriques
en utilisant pourtant un produit scalaire non pondéré, mais ils utilisaient
une approximation de Boussinesq pour le fluide, ce qui rend en fait
les deux produits scalaires équivalents. \cite{Ulia:91} avaient également
déjà remarqué l'essentiel à savoir
que ``la fonction d'influence peut se calculer à partir
de rétro-trajectoires particulaires quand on utilise un modèle Lagrangien.
Les modèles Eulériens, gouvernés par des équations aux dérivées partielles,
sont formulés dans un cadre variationnel et, dans ce cas,
la fonction d'influence peut
s'obtenir comme solution de l'équation adjointe à rebours dans le temps,
en utilisant le récepteur comme source."
Aux vues des considérations ci-dessus, on est tenté de récrire ce paragraphe
sous une forme un peu plus symétrique~: la fonction d'influence peut
être obtenue soit en suivant à rebours dans le temps la masse d'air
échantillonnée au niveau du détecteur (rétro-transport) soit comme solution 
de l'équation adjointe, et ce indépendamment du cadre (Eulérien ou Lagrangien)
choisi pour représenter le transport atmosphérique.

La symétrie du transport à la base de cette équivalence peut cependant
être perdue dans le monde numérique, c'est à dire que
le modèle numérique intégré à rebours dans le temps peut être différent
de l'adjoint du code numérique direct.
En particulier, seuls les opérateurs
linéaires par rapport à la concentration du traceur ont un adjoint et peuvent
donc respecter la symétrie temporelle.
Les schémas sophistiqués utilisés aujourd'hui pour l'advection des traceurs
dans la plupart des modèles introduisent des non linéarités pour garantir
un meilleur comportement physique (conservation, positivité, monotonie,
faible diffusion).
Ces schémas ne sont donc pas symétriques.
C'est le cas en particulier du schéma de Van Leer utilisé
dans LMDZ.
On montre en revanche que les codes numériques utilisés dans LMDZ pour calculer
le transport turbulent ou convectif respectent bien l'équivalence entre
transport rétro et adjoint.

Dans le cas où cette symétrie n'est pas respectée, la question se pose
de savoir s'il est préférable d'utiliser le code adjoint du code direct
ou d'intégrer à rebours dans le temps le modèle direct.
Dans le cas des schémas d'advection, l'utilisation du rétro-transport
garantira la positivité des rétro-panaches alors que le calcul adjoint
pourra faire apparaître des valeurs négatives non physiques (une mesure
de concentration pouvant augmenter à la suite d'une émission moindre).
Dans certains algorithmes d'inversion, la positivité de la source peut
être une contrainte importante, auquel cas on aura tendance à privilégier
le rétro-transport. Pour des algorithmes de minimisation linéaires, utilisant
des descentes de gradient, on peut penser a contrario que le code
adjoint, qui fournit un calcul exact du gradient, doit être privilégié.
On montre ici, sur des cas académiques, que cette conclusion
est un peu trop rapide et que, même dans ce cas, l'utilisation d'un
gradient approché mais garantissant certaines propriétés physiques peut
permettre une minimisation moins poussée mathématiquement mais plus
robuste.


Dans ce chapitre, nous revenons en détail sur l'ensemble des aspects introduits
ci-dessus.
Nous commençons par introduire (\sec{transport}) le rétro-transport et
la symétrie du transport atmosphérique
pour des sources et détecteurs ponctuels et
montrons comment le rétro-transport ainsi défini
s'étend aux processus diffusifs.
Nous expliquons dans la même section le lien avec l'équation adjointe et étendons
la théorie à des sources diffuses.
Nous explicitons en particulier
ce lien entre rétro-transport et équations adjointes
sur le cas des schémas de convection en flux de masse (\sec{param}).
La discussion, développée d'abord dans le monde analytique, est ensuite
étendue au monde numérique en présentant tout d'abord une illustration
numérique du rétro-transport dans le cas de la campagne ETEX (\sec{numeric}).
On démontre la symétrie des algorithmes utilisés dans LMDZ pour le transport
turbulent et convectif.
On étudie ensuite le cas des schémas non symétriques en se focalisant
sur le schéma de Van Leer utilisé dans LMDZ (\sec{inversion}).
Le chapitre se termine par quelques illustrations relatives à des applications
menées autour de l'utilisation inverse du modèle LMDZ, et notamment
des résultats relatifs à la surveillance des essais nucléaires.

Ce travail doit beaucoup aux discussions initiales avec Robert Sadourny.
La sollicitation et le financement des études par Jean-Pierre Issartel puis
Philippe Heinrich (CEA/DAM) ont permis de mener le travail au-delà du simple cadre
académique.
Ce travail a aussi bénéficié de nombreuses et fructueuses discussions
avec Bertrand Cabrit. Les tests et illustrations ETEX ont
été réalisés par Abderrahmane Idelkadi pendant sa thèse.
Le travail sur le TICE a vu  passer plusieurs stagiaires
comme Alexandre Maes et Elie Anselin.
Enfin Olivier Talagrand a été d'une aide précieuse pour mener
à bien la mise en forme de la partie adjointe et
l'écriture d'un article en deux parties sur
le lien entre transport rétro et adjoint \cite[]{Hour:05I,Hour:05II},
article sur lequel repose largement le présent chapitre.

