Introduction

Avant propos

Deux objectifs, en partie contradictoires, m'ont guidé dans la rédaction de ce document. Le premier est d'ordre statutaire : pour soutenir une Habilitation à Diriger des Recherches, il faut rédiger une synthèse de ses travaux antérieurs, souvent déjà publiés en anglais dans des revues à comité de lecture. Le second objectif que je me suis fixé est de profiter de cet exercice imposé pour essayer d'écrire un texte qui puisse servir plus ou moins directement de support de travail pour des collègues, voir de support pédagogique sur différents sujets auxquels je me suis intéressé au cours des dix dernières années. Faute d'avoir pu y consacrer suffisamment de temps, ces objectifs ne sont certainement que partiellement atteints. Les rapporteurs et examinateurs ne manqueront pas de trouver le document trop long. Les collègues ou étudiants voulant se renseigner sur un des sujets abordés le trouveront incomplet ou elliptique ; déséquilibré entre des parties trop techniques et d'autres trop vagues. J'espère que l'exercice restera malgré tout utile. Parmi les efforts pour accroître la lisibilité, j'ai essayé de réduire au maximum le nombre des acronymes utilisés et de répertorier en fin de document ceux qui sont utilisés.

Le texte est composé de quatre parties plus ou moins indépendantes, avec à chaque fois une introduction assez complète et des conclusions et perspectives. L'introduction et la conclusion du document sont donc plutôt là pour brosser le cadre général du travail pour l'une et tracer des perspectives générales pour l'autre. Le fil conducteur de ces différentes parties est le transport atmosphérique et la modélisation de ce transport dans les modèles globaux de climat.

Couplages entre composition et transport atmosphérique

Les couplages entre composition et transport atmosphérique occupent une place grandissante dans l'étude des atmosphères planétaires.

Ces couplages sont tout d'abord au centre d'une partie des questions relatives au changement climatique. En effet, une part des incertitudes relatives au réchauffement global du climat terrestre provient des incertitudes sur l'évolution de la composition même de l'atmosphère (CO$_2$, méthane, ozone troposphérique, aérosols). Or l'évolution de cette composition est étroitement liée au transport atmosphérique et au climat. Pour le CO$_2$, l'augmentation des concentrations atmosphériques sous l'effet des émissions anthropiques conduit à une augmentation du stockage dans les océans et les écosystèmes. Le puits biosphérique est lui-même sensible à l'évolution du climat. Les estimations actuelles prévoient une réduction du puits biosphérique consécutive au changement climatique qui pourrait correspondre à une rétroaction positive de 15$\%$ sur la teneur en CO$_2$ de l'atmosphère (, ). La modification de l'ozone troposphérique est également étroitement couplée à l'évolution de la température et de l'humidité. Les changements des vents en surface peuvent également modifier le soulèvement des poussières désertiques ou les émissions de DMS (précurseurs des aérosols soufrés) par les océans (, ).

Au-delà des modifications de la composition, la sensibilité du climat à un changement de concentration imposé est également souvent conditionnée par des processus de transport. La fameuse rétroaction vapeur d'eau (l'atmosphère plus chaude se charge en vapeur d'eau, ce qui augmente en retour l'effet de serre) peut par exemple être fortement modulée par les processus de transport. Si on suppose par exemple que le réchauffement résulte également en une augmentation de l'altitude de la pénétration de la convection (convection nuageuse ou grands systèmes de Hadley-Walker), l'air subsident autour de ces zones convectives sera au contraire plus sec. Même si les processus de transport sont plus complexes (cf. par exemple , ), il n'en reste pas moins que l'accroissement d'humidité lors d'un réchauffement climatique peut être modulé par les changements d'advection. Les effets indirects des aérosols (une augmentation du nombre de noyaux de condensation résultant en des nuages formés de plus petites gouttes, donc plus brillants et moins précipitants) sont également une source importante d'incertitude sur l'amplitude du changement climatique, étroitement liée au transport et à la microphysique.

Concernant les couplages entre composition et dynamique atmosphérique sur les autres planètes, on mentionnera d'abord le cas des poussières sur Mars. En dehors des bandes d'absorption du CO$_2$, constituant majoritaire de la fine atmosphère martienne, la poussière est le principal constituant actif radiativement. Cette poussière est en permanence soulevée du sol du grand désert martien par des rafales de vent, des tornades ou de petites tempêtes locales. Régulièrement, des tempêtes plus importantes se déclenchent, soulevant la poussière sur des milliers de kilomètres. A certaines périodes de l'année, ces tempêtes peuvent finalement dégénérer en évènements globaux spectaculaires, au cours desquels la surface de Mars est entièrement voilée pour l'il d'un observateur extérieur. La circulation est alors profondément modifiée (, ). Notons enfin les couplages entre la photochimie, la microphysique des brumes et la circulation stratosphérique sur Titan qui ont largement retenu notre attention au cours des années passées et qui font l'objet d'un chapitre du présent document.

Modélisation du climat

La modélisation numérique globale est devenue un outil de base pour aborder ces systèmes complexes. Les modèles de circulation générale atmosphérique, développés au début des années 70 pour les besoins de la prévision météorologique, se sont petit à petit enrichis tant sur le plan physique (représentation des nuages, modèles thermodynamiques du sol, paramétrisation de la convection) que par la prise en compte du couplage avec les autres composantes du système climatique. On pense en particulier pour la Terre au couplage avec l'océan, la végétation, la chimie ou les aérosols. Ces développements ont abouti dans les années récentes au concept de ``modèles intégrés du climat" (les soi-disant ``Earth system models") utilisés notamment pour essayer de prévoir les évolutions futures du climat. Un modèle de ce type est actuellement développé et utilisé à l'IPSL. Il comprend, couplé au modèle de circulation générale atmosphérique LMDZ, le modèle de circulation générale océanique ORCALIM, le modèle des surfaces continentales ORCHIDEE et le module aérosols-chimie INCA. Ce modèle est actuellement impliqué dans la réalisation de ``scénarios climatiques" pour le prochain rapport du GIEC.

Une des originalités de la recherche menée au LMD est d'étudier de front et avec le même outil, LMDZ, le climat de la Terre et celui d'autres planètes du système solaire, en particulier Mars et Titan. En parallèle du développement du modèle intégré terrestre, et souvent même avant, les versions planétaires ont connu des évolutions similaires. Sur Mars, ce sont les couplages avec le cycle du carbone (un quart de l'atmosphère de CO$_2$ se condense saisonnièrement dans les calottes polaires), des poussières (on en a parlé plus haut) et de l'eau (avec l'enjeu de déterminer les réservoirs d'eau sous la surface et de comprendre les évolutions passées du climat de la planète rouge) qui ont été inclus dans les modèles de circulation existants (, ,,,). Pour Titan, les couplages avec la photochimie et la brume (, ,,) ou le méthane (, ) ont également été inclus dans les modèles existants.

Représentation du transport des traceurs dans les modèles globaux

Pour la modélisation des couplages entre composition et climat, une étape essentielle du travail de développement consiste à introduire, dans le modèle de circulation générale atmosphérique, les algorithmes permettant de représenter le transport des espèces traces. Il faut traiter à la fois le transport par l'écoulement explicitement représenté dans le modèle de circulation (c'est à dire pour des échelles supérieures à quelques centaines ou dizaines de kilomètres) et le transport par les écoulements non résolus, turbulents ou convectifs. C'est en fait l'introduction, dans le modèle LMDZ, de schémas permettant de représenter le transport à grande échelle qui a donné le coup d'envoi à l'ensemble des études présentées ici. Nous avons plus particulièrement mis en uvre dans le modèle LMDZ des schémas en volumes finis développés à l'origine par (). Les versions terrestre et planétaires étant développées de front, ce travail d'introduction du transport a été effectué de fait à la fois pour la Terre, Mars et Titan. Cette version avec transport des traceurs du modèle LMDZ est à l'origine d'un grand nombre de développements et applications concernant ces trois planètes. Ce travail préliminaire est décrit dans le Chapitre 2.

Je me suis ensuite plus particulièrement intéressé au transport vertical par les structures méso-échelles de la couche limite convective. Ces structures (rouleaux, cellules thermiques), bien connues des amateurs de vol libre (deltaplanes, planeurs, parapentes), ne sont en général pas considérées de façon spécifique dans les modèles de circulation globaux, qui privilégient, à un bout, une vision en diffusion du transport turbulent dans la couche limite et, à l'autre bout, des schémas de convection profonde, contrôlés pour une bonne part par les changements de phase de l'eau. Les modèles doivent du coup inclure des traitements adhoc (``contre-gradients", ``ajustement convectif") pour pallier l'absence de paramétrisation des structures convectives de couche limite ; ces structures convectives qui, sur des régions désertiques ou sur une planète-désert comme Mars, peuvent dominer le transport vertical jusqu'à plusieurs kilomètres au-dessus de la surface. Cette nouvelle paramétrisation, le ``modèle du thermique", basée sur une formulation dite ``en flux de masse", est décrite en détail dans le Chapitre . Y sont également présentées différentes validations par rapport à des simulations numériques des grands tourbillons ou à des observations.

Calculs de dispersion et surveillance de l'environnement

En marge des questions relevant directement des couplages entre composition et climat, les outils développés pour le transport des espèces traces peuvent être utilisés pour étudier la dispersion de polluants atmosphériques à écoulement atmosphérique connu. On présente dans les chapitres 2 et de tels calculs de dispersion, réalisés ici à des fins de validation des algorithmes de transport. Dans ces simulations numériques, on force le modèle de circulation atmosphérique à suivre au plus près la situation synoptique afin de comparer, au jour le jour, les concentrations observées de certains constituants avec des données de terrain. La technique employée est baptisée ``nudging" en Anglais ce qu'on traduira ici par ``guidage". Cette technique consiste à relaxer en permanence les champs météorologiques du modèle vers des données d'``analyses" ou ``réanalyses" produites par les grands centres de prévisions météorologiques. Dans ce cas, le modèle de circulation joue le rôle d'une espèce d'interpolateur physique sur le maillage choisi et permet de recalculer un jeu cohérent de grandeurs physiques nécessaires à la représentation du transport grande échelle et sous-maille. Les grandeurs nécessaires pour les algorithmes de transport des traceurs peuvent être soit passées directement aux algorithmes concernés, soit stockées dans des fichiers puis relues pour le seul calcul du transport des espèces traces. On parlera de modes ``branché" et ``débranché".

Il s'avère que le modèle guidé et débranché, développé à l'origine pour des besoins de validation dans le cadre des études couplées chimie-climat, est un modèle global de dispersion atmosphérique parfaitement adapté à certaines questions relatives à la surveillance de l'environnement. Un travail spécifique a été entrepris dans ce domaine suite à une demande du CEA relative à la surveillance des essais nucléaires à partir de la mesure de la concentration en éléments radioactifs. Il s'agit d'un cas classique de problème inverse dans lequel on veut obtenir des contraintes sur les sources à partir de mesures de concentration. Avec Robert Sadourny, nous nous sommes convaincus, à l'époque de la demande, qu'il était légitime, pour répondre à cette question, d'inverser la direction du temps dans le modèle de transport Eulérien (débranché). En émettant un traceur au niveau des détecteurs, le modèle calcule alors la distribution d'origine de l'air échantillonné à la station.

L'utilisation du transport inverse et de ``modèles de détecteurs" pour aborder ce type de question n'est pas nouvelle. Mais le fait qu'on puisse utiliser directement les codes Eulériens à rebours dans le temps (ce que nous appellerons le ``rétro-transport Eulérien'') ne semblait pas vraiment acquis. Ce travail sur la détection des essais nucléaires a donc été l'occasion de clarifier la théorie sous-jacente. Le rétro-transport Eulérien peut être présenté d'un point de vue physique comme une formulation Eulérienne de l'approche des rétro-trajectoires Lagrangiennes, largement utilisée dans la communauté des chimistes de l'atmosphère pour interpréter des mesures de composition ponctuelles. Le rétro-transport Eulérien peut également être présenté d'un point de vue mathématique comme l'adjoint du transport direct pour un produit scalaire particulier, pondéré par la masse de l'air sous-tendant le transport. Les visions à la base des rétro-trajectoires et de l'approche adjointe sont cependant suffisamment différentes pour que les outils développés le soient aussi, avec des conséquences importantes sur l'efficacité des algorithmes d'inversion. Le Chapitre  présente de façon détaillée à la fois les aspects théoriques, des illustrations numériques et des exemples d'application pour la surveillance des essais nucléaires. Les outils développés à cette occasion sont en cours d'intégration dans une chaîne opérationnelle au CEA.

Couplages entre composition et dynamique sur Titan

La dernière partie de ce document concerne Titan, le plus gros satellite de Saturne. Titan fait partie de ces objets fascinants du système solaire révélés par l'épopée Voyager. En 1981, les responsables des missions Voyager choisissent de privilégier pour la sonde Voyager 1 un survol de Titan plutôt que de poursuivre la course vers Uranus et Neptune (périple magistralement réussi ensuite par Voyager 2). On sait en effet à l'époque que Titan est, avec la Terre, le seul corps tellurique du système solaire entouré d'une atmosphère dense d'azote (1,5 bar à la surface). Les photos renvoyées vers la Terre sont très décevantes. Une épaisse couche de brume orangée voile entièrement la surface. Tout juste peut-on distinguer un léger contraste entre les deux hémisphères, signe probable d'un effet saisonnier. Les mesures spectroscopiques permettent en revanche d'identifier un grand nombre de composés chimiques, hydrocarbures et nitriles. Ces espèces chimiques, créées dans la très haute atmosphère à partir de la photo-dissociation de l'azote moléculaire et du méthane (second constituant atmosphérique) sont ensuite transportées vers le bas dans la stratosphère où on pense qu'elles polymérisent pour donner naissance à la brume orange. L'analyse des contrastes latitudinaux de température dans la stratosphère suggère également que l'atmosphère tourne beaucoup plus vite que le satellite, lui-même en phase bloquée autour de Saturne, avec une durée du jour de 16 jours terrestres environ. Si la direction de la rotation de l'atmosphère ne peut être obtenue à partir des observations de la température, l'analogie avec Vénus et des arguments théoriques suggèrent que l'atmosphère est en régime de ``superrotation", l'atmosphère vers 200 km tournant une dizaine de fois plus vite que la surface et dans la même direction.

Suite au passage des sondes Voyager, une mission est programmée vers le système de Titan, sous l'impulsion de Toby Owen et Daniel Gautier. La sonde américaine Cassini se consacrera au système de Saturne et emmènera à son bord la sonde européenne Huygens qui plongera dans l'atmosphère de Titan. Mission parfaitement remplie. Le 14 janvier 2005, après 7 années de voyage dans le système solaire, la sonde Huygens a dévoilé sous l'épaisse couche de brume des paysages familiers où chacun reconnaît qui sa Côte d'Azur, qui son lac de montagne ; en tout cas des images qui évoquent un rivage.

Au début des années 90s, sous l'impulsion de Daniel Gautier (LESIA) et Christopher P. McKay (NASA/Ames), différents travaux de modélisation sont entrepris pour interpréter les résultats des missions Voyager et préparer la mission Cassini-Huygens. Le modèle de circulation du LMD est adapté aux conditions de Titan (, ) et prédit effectivement une forte superrotation sur Titan, superrotation confirmée depuis par des mesures Doppler. En parallèle, des modèles unidimensionnels sont développés pour la photochimie (, ) et la microphysique des brumes (, ,). Il apparaît cependant rapidement que les différentes composantes de ce système sont fortement couplées. Les brumes sont formées par la polymérisation des constituants chimiques et peuvent servir également de noyaux de condensation à ces dernières au niveau de la troposphère glaciale de Titan (70 K environ). Brumes et espèces chimiques sont évidemment transportées par les vents. En retour, les contrastes latitudinaux de la composition jouent les premiers rôles dans le forçage de la circulation. Vers 1995, nous décidons, avec Michel Cabane et Dominique Toublanc, de réunir les différents efforts de modélisation pour s'attaquer à ce système climatique complet.

En 1998, lors du colloque quadriénal du Programme National de Planétologie, le programme est déjà clair (, ) : ``L'arrivée sur Titan de la mission Cassini-Huygens est sans doute une des dernières occasions avant des décennies d'explorer un système physique analogue à la Terre mais encore très mal connu. Pour l'atmosphère et le climat en particulier, c'est une occasion unique avant longtemps de mettre à l'épreuve pour une planète tellurique les théories et modèles développés dans le contexte terrestre. Cette perspective ainsi que la préparation de la mission (étude en amont et préparation de l'analyse des résultats) ont motivé le développement d'un modèle de circulation générale de l'atmosphère de Titan au Laboratoire de Météorologie Dynamique du CNRS, sous l'impulsion de Daniel Gautier et en collaboration avec Christopher P. McKay (NASA/Ames) et Régis Courtin (LESIA/Obs. Paris Meudon). [...] Pour Titan, le modèle prédit une stratosphère tournant environ 10 fois plus vite que la planète solide avec des vents zonaux (d'ouest) de l'ordre de 100 m s$^{-1}$. En plus de ce phénomène dynamique spectaculaire, les résultats du modèle ont contribué à mettre en évidence l'importance des couplages entre dynamique atmosphérique, microphysique des aérosols, et photochimie. Ceci nous a conduit à bâtir et à proposer au PNP pour les années à venir, un projet de modélisation du climat de Titan intégrant ces différentes composantes. L'enjeu est d'importance et la tâche ardue quand on connaît les problèmes rencontrés dans la modélisation de ces problèmes sur Terre. Mais la perspective de la confrontation du modèle aux observations de la mission Cassini-Huygens en 2005 en font [un objectif] scientifique de tout premier plan.''

Après une mise en route souvent ardue, le modèle couplé a bien été développé. Il a permis de mettre en évidence le couplage très fort entre la brume et les vents (, ). Le modèle a également permis d'expliquer les contrastes latitudinaux observés dans la composition chimique (, ,). Sur la base de simulations numériques effectuées avec cet outil, une base de données a été constituée et mise à disposition de la communauté sur la toile avant l'arrivée de la mission (, ). Les données sont là. Les premières photos ont réservé leur lot de surprises. Que nous réservent les dépouillements en cours des enregistrements des spectromètres et "imageurs spectraux" ou des mesures in-situ de la composition par Huygens ?

C'est cette histoire qui clôt ce document (Chapitre ) avant quelques conclusions générales.

Représentation du transport des espèces traces dans un modèle de climat global

Le présent chapitre est consacré à la représentation du transport atmosphérique d'espèces traces dans les modèles de circulation atmosphérique de grande échelle ainsi qu'à la présentation d'un outil particulier : la version ``traceurs" du modèle de circulation générale LMDZ.

Le LMD développe et exploite depuis le début des années 70 un modèle de circulation générale atmosphérique. Comme beaucoup d'autres, ce modèle s'est petit à petit enrichi pour devenir un véritable modèle climatique, avec par exemple la prise en compte des couverts végétaux pour prédire le comportement thermodynamique des surfaces continentales ou le couplage avec l'océan. Le transport d'espèces traces est pour sa part introduit une première fois par Sylvie () dans une version précédente du modèle de climat du LMD pour étudier le cycle des poussières désertiques. Christophe Genthon introduit ensuite le Radon ($^{222}$Rn) et le plomb ($^{210}$Pb) dans le modèle en utilisant, pour représenter l'advection, le schéma des Pentes du NASA/GISS (, ,).

Au milieu des années 80, Robert Sadourny et Phu LeVan entreprennent la réécriture du noyau hydrodynamique du modèle, afin de le rendre plus modulaire, lisible et efficace (l'ancien modèle avait été écrit sur cartes perforées par Phu LeVan pour des machines ne pouvant pas contenir un champ entier en mémoire) et de généraliser l'idée de grille étirable qui avait été une première fois testée pour étudier un cyclone en baie du Benghal. C'est cette possibilité de raffinement de la grille qui donnera plus tard le "Z" (pour zoom) de LMDZ. Du fait de l'inertie inhérente à la modélisation du climat terrestre (on regarde des choses très précises avec un modèle robuste dont on connaît bien le fonctionnement et dont on a ``réglé" la climatologie), ce nouveau noyau dynamique est dans un premier temps utilisé sur Mars (, ,) et Titan (, ,). Pour ces deux planètes, un des prolongements possibles des études développées au LMD consistait à s'intéresser au cycle d'espèces transportées : les poussières sur Mars, avec les spectaculaires tempêtes globales qui peuvent voiler la surface de la planète pendant plusieurs dizaines de jours, et les brumes et espèces chimiques sur Titan, dont on montre dans le dernier chapitre de ce document qu'elles sont fortement couplées à la circulation atmosphérique.

Les traceurs sont introduits dans LMDZ en 1996, en utilisant, pour l'advection de grande échelle, des schémas en volumes finis^2.1. Les schémas en volumes finis sont basés sur une partition du domaine considéré en volumes de contrôle, aux frontières desquels on évalue les flux entrants ou sortants de traceurs. Ce sont des schémas basés sur une formulation intégrale de l'équation de transport. Nous avons plus précisément codé et testé une série de schémas en volumes finis proposés à l'origine par (). Ces schémas conduisent facilement à une mise en uvre tridimensionnelle et satisfont des propriétés physiques fondamentales du transport: localité, conservation, monotonie, positivité (plus généralement pas de création d'extrema numériques) et invariance par addition d'une constante au champ de traceur. () avait en fait proposé une hiérarchie de schémas dont les plus sophistiqués ont été introduits ultérieurement et indépendamment dans la communauté météorologique par () (schéma des pentes du NASA/GISS) et (). Avec Alexandre Armengaud alors en thèse au LGGE sous la direction de Christophe Genthon, nous avons testé dans LMDZ plusieurs de ces schémas (, ). Nous avons pu montrer que leurs performances étaient en fait assez semblables dès lors qu'on les comparait non pas à résolution spatiale fixée mais à coût numérique équivalent (un schéma plus précis mais plus coûteux se comporte comme un schéma moins précis mais utilisé sur une grille plus fine). Nous avons retenu pour le modèle du LMD le schéma le plus simple (le schéma I dans l'article original de Van Leer souvent appelé MUSCL ou MINMOD).

Ce travail a donné naissance à la version traceurs du modèle, baptisée un moment LMDZT. En parallèle des schémas d'advection, il a fallu inclure, dans le modèle, le transport associé aux paramétrisations des mouvements non résolus, turbulents ou convectifs. Cette composante est souvent essentielle pour contrôler le transport vertical des espèces. Pour ces schémas, on suit généralement ce qui est fait pour transporter l'humidité ou la température potentielle dans les paramétrisations d'origine. L'introduction des traceurs dans les paramétrisations de la turbulence de couche limite et de la convection nuageuse a été initialement réalisée par Olivier Boucher et Alexandre Armengaud. Cette composante traceurs fait depuis partie intégrante du modèle LMDZ qui permet de transporter un nombre arbitraire de traceurs. De nombreuses applications ont été développées à partir de cette version, à la fois pour la Terre, Mars ou Titan, en lui adjoignant des modules de chimie ou de microphysique des aérosols.

LMDZT a été conçu de façon très modulaire de façon à pouvoir y rajouter facilement des codes de chimie ou des modules d'aérosols. Il est utilisable soit en mode ``branché'' (on-line en anglais), soit en mode débranché, en relisant des fichiers météorologiques issus d'une simulation numérique précédente. Les simulations météorologiques elles-mêmes peuvent être effectuées soit en mode climatique, en laissant le modèle évoluer librement à partir d'une condition initiale unique, soit en le ``guidant'' par des analyses météorologiques. Marie-Angèle Filiberti et Abderrahmane Idelkadi ont largement contribué à fiabiliser et valider l'ensemble des outils développés et présentés ici.

Ce chapitre relativement technique revient à la fois sur la présentation des schémas de transport et sur la description de LMDZT. Nous présentons à la fin quelques exemples de validations ou utilisations du modèle.

Le modèle LMDZ

Avant d'introduire la composante traceurs, on présente rapidement le modèle LMDZ et notamment la discrétisation des équations de grande échelle sur laquelle s'appuie le transport des traceurs. On en profite pour donner un bref aperçu du contenu des paramétrisations physiques du modèle et des évolutions apportées récemment. ^2.2

LMDZ correspond à la seconde génération d'un modèle de climat développé depuis une trentaine d'années au LMD et décrit initialement par (). Ce modèle est plus modulaire et flexible^2.3que son prédécesseur. Le ``Z" de son nom se réfère à la capacité de raffinement de la grille (Zoom) obtenue grâce à une écriture généralisée de la formulation numérique avec un maillage dont les facteurs d'élongation dans les deux directions horizontales peuvent être choisis arbitrairement.

Description générale

Comme la plupart des modèles de climat, le modèle LMDZ intègre sur la sphère et dans le temps les ``équations primitives de la météorologie". Ces équations sont une version des équations de Navier Stokes simplifiées en supposant que l'atmosphère est à tout moment en équilibre hydrostatique sur la verticale et en négligeant les variations verticales de la géométrie horizontale (hypothèse de couche mince). Le moment cinétique par rapport à l'axe des pôles est par exemple calculé en utilisant comme distance à l'axe $a \cos\phi$ - où $a$ est le rayon de la planète et $\phi$ la latitude - plutôt que $(a+z) \cos\phi$, qui tiendrait compte de l'altitude $z$ d'une particule d'air au dessus de la surface.

Le modèle du LMD est bâti sur une discrétisation en différences finies de ces équations avec certaines propriétés intéressantes du schéma numérique comme la conservation de la masse, la conservation du moment cinétique par la composante axi-symétrique de l'écoulement et de la vorticité barotrope. On revient sur ces aspects en toute fin de chapitre.

Ces équations ne peuvent pas être intégrées jusqu'à l'échelle visqueuse. En pratique, dans les modèles globaux, on utilise des mailles de quelques dizaines à quelques centaines de kilomètres suivant les applications. L'impact des échelles sous-mailles sur la grande échelle doit donc être représenté au travers de paramétrisations. Il faut également représenter des processus fondamentaux comme le transfert de rayonnement visible et infrarouge dans l'atmosphère, les processus nuageux ou les interactions avec la surface. Dans le jargon des modélisateurs, ces paramétrisations sont regroupées dans la ``physique" du modèle, par opposition avec le code représentant la ``dynamique" de grande échelle. Dans le cas de l'utilisation du modèle de circulation pour d'autres planètes, c'est cette partie physique qui doit être largement modifiée, et notamment le calcul du transfert radiatif.

Discrétisation des équations primitives

Figure 2.1: Disposition des variables pour la grille du modèle LMDZ.

Dans le modèle LMDZ, les équations primitives sont discrétisées horizontalement sur une grille-C dans la classification d'Arakawa (cf. e. g. , ). On note $X$ et $Y$ les coordonnées horizontales : $X$ (resp. $Y$) est une fonction biunivoque de la longitude $\lambda$ (resp. de la latitude $\phi$). Les variables scalaires (la température potentielle $\theta$, le géopotentiel $\Phi$ et la pression de surface $p_s$) sont évaluées aux points correspondant à des couples de valeurs entières $(X,Y)=(i,j)$. Le vent zonal est calculé aux points $(X,Y)=(i+1/2,j)$ et le vent méridien aux points $(X,Y)=(i,j+1/2)$. La disposition des variables sur la grille est illustrée sur la Fig. 2.1. De façon à pouvoir modifier la distribution des longitudes et latitudes de la grille, on utilise en fait les composantes covariantes ($\tilde{u}$ et $\tilde{v}$) et contravariantes ( $\tilde{\tilde{u}}$ et $\tilde{\tilde{v}}$) du vent définies par

$$\begin{array}{llllllllll} \tilde{u}= c_uu & \mbox{et} & \til... ...c_v& \mbox{avec} & c_v= a \left( d\phi / dY \right) \end{array}$$ (2.1)

où $u$ et $v$ sont les composantes physiques du vecteur vent horizontal.

Comme beaucoup de modèles globaux, le modèle du LMD avait initialement été codé avec une coordonnée verticale de type $\sigma$. La coordonnée $\sigma=p/p_s$ (où $p$ est la pression et $p_s$ la pression à la surface au point considéré) est pratique parce qu'elle varie de 1 à la surface à 0 au sommet de l'atmosphère quelque soit le relief sous-jacent. Cette coordonnée devient cependant problématique plus haut dans l'atmosphère où on préfère travailler sur des isobares ou, mieux, sur des isentropes. Suivant là aussi beaucoup de modèles de climat, nous avons avec Phu LeVan récrit le modèle en coordonnée hybride $\sigma-p$. La coordonnée hybride est définie de façon implicite comme donnant la pression dans la couche $l$ du modèle sous la forme

$$p_l=A_l + B_lp_s$$ (2.2)

On choisit près de la surface $A\sim 0$ et $B\sim 1$. Vers le haut du modèle, on a $A\sim 0$ et $B\sim 0$ avec $A\gg B p_s$. On retrouve ainsi la coordonnée $\sigma$ près de la surface et des coordonnées pression plus haut dans l'atmosphère.

Les niveaux de pression définis par la relation 2.2 correspondent aux interfaces entre les couches du modèle avec $p_1=p_s$ et $p_{N+1}=0$, $N$ étant le nombre de couches dans le modèle. La masse d'air contenue dans une maille du modèle comprise entre les niveaux $p_k$ et $p_{k+1}$ est donnée par

$$m_k={\cal A}\frac{p_k-p_{k+1}}{g}$$ (2.3)

( ${\cal A}=c_uc_v$ est l'aire de la maille) ou encore

$$m_k= \frac{{\cal A}}{g} \left[A_k-A_{k+1} +\left(B_k-B_{k+1}\right) p_s\right]$$ (2.4)

On introduit alors les trois composantes du flux de masse :

$$U={\overline{m}}^{ X } \tilde{\tilde{u}}, V= {\overline{m}}^{ Y } \tilde{\tilde{v}} \mbox{et} \ W$$ (2.5)

où la notation ${\overline{a}}^{ X }$ signifie qu'on prend la moyenne arithmétique de la quantité $a$ suivant la direction $X$ et où le flux de masse vertical $W$ est défini à partir de l'équation de continuité

$$\frac{\partial m}{\partial t} + \delta_x U+ \delta_y V +\delta_z W=0$$ (2.6)

la notation $\delta_X$ correspondant à la différence entre deux points consécutifs suivant la direction $X$.^2.4

On ne présente pas ici la discrétisation des équations du mouvements qui n'a pas changé par rapport à la version $\sigma$ du modèle et dont ont discute en toute fin de chapitre. En revanche, l'introduction de la coordonnée hybride à entraîné la modification du schéma de calcul de la fonction d'Exner $\Pi =p^\kappa$.^2.5L'équilibre hydrostatique est intégré verticalement suivant le schéma :

$$\delta_z \Phi=- {\overline{\theta }}^{ z } \delta_z \Pi$$ (2.11)

où $\Phi$ est le géopotentiel, avec ${\overline{\theta }}^{ z }_1=\theta _1$ et ${\delta_z \Pi }_1=p_s^\kappa-\Pi _1$. Au lieu de calculer $\Pi _k$ au milieu de la couche $k$ à partir d'une moyenne arithmétique des pressions $p_k$ et $p_{k+1}$ (ce qui est entre autres coûteux numériquement), on utilise une relation numérique différente entre $p$ et $\Pi$ :

$$\frac{{\cal A}}{g} {\overline{p\delta_z \Pi }}^{ z }=-\kappa \Pi m$$ (2.12)

Cette relation garantit un équivalent numérique de la relation

$$\int_0^\infty \Phi \rho dz =\int_0^\infty RT \rho dz$$ (2.13)

entre énergie potentielle et énergie interne, à savoir ^2.6

$$\sum_{l=1}^N\Phi_lm_l -\Phi_s = \sum_{l=1}^N RT_lm_l = \sum_{l=1}^N \kappa \theta _l \Pi _l m_l$$ (2.21)

Les paramétrisations physiques du modèle de climat terrestre

La première utilisation de ce nouveau code hydrodynamique est martienne. Le modèle du LMD est le premier à simuler un cycle saisonnier entier sur la planète rouge (, ). Pendant ce temps le modèle climatique terrestre du LMD poursuit une vie trépidante. La version 4ter, qui donnait des résultats excellents sur la Mousson Indienne (présentez une simulation du modèle à un chercheur authentique du LMD, il commencera par regarder la ``précip'' de juillet sur l'Inde), est progressivement remplacée par les versions 5 et 6. Ces versions sont notamment couplées à SECHIBA pour la végétation et à ORCALIM pour l'océan.

C'est lors d'un séminaire interne du LMD, tenu en la royale abbaye de Fontevraud en 1990, qu'est prise collectivement la décision de faire de ce nouveau code dynamique l'ossature du futur modèle de climat du LMD. Ce choix est motivé par la volonté de disposer d'un outil souple et modulaire, permettant facilement l'échange et le test de modules et de procédures. Cette décision s'inscrit dans l'idée de modélisation communautaire alors dans l'air du temps.

Il faut donc adjoindre au code dynamique un jeu de paramétrisations physiques. Dans un premier temps, la décision est prise de partir de la physique du Centre Européen pour les Prévisions Météorologiques à Moyen Terme. Ce modèle contient en effet le code radiatif Fouquart/Morcrette, développé à l'origine pour le modèle du LMD, et les schémas de nuages d'Hervé LeTreut. Cette décision ne sera finalement pas tenue. Si Laurent Li s'est bien battu pour effectuer les portages nécessaires, le groupe décidera quelques années plus tard de privilégier la continuité avec les versions précédentes.^2.7C'est la combinaison de ce nouveau code dynamique et de la réécriture des paramétrisations physiques qui donnera naissance au modèle LMDZ d'aujourd'hui.

On décrit ci-dessous dans ses grandes lignes la version du modèle LMDZ3.3 qui a été utilisée pour les applications "traceurs" décrites dans ce document ainsi que les évolutions apportées au court de la mise au point du modèle couplé IPSLCM4 et qui définissent la version LMDZ4 (, ).

L'effet de la turbulence de petite échelle dans la couche limite est pris en compte au moyen d'une fermeture en super-viscosité avec un flux turbulent proportionnel au gradient vertical de la quantité transportée. La viscosité turbulente dépend du cisaillement vertical du vent et du nombre de Richardson, selon les formules présentées par (). Un "contre-gradient" est introduit sur la température potentielle pour permettre un transport de chaleur "en remontant le gradient" dans les couches limites marginalement stables. On revient très largement sur ces aspects relatifs à la couche limite dans le Chapitre .

Le code radiatif est celui du modèle du Centre Européen pour les Prévisions Météorologiques à Moyen Terme (, ). La partie concernant l'absorption et la diffusion du rayonnement solaire est une version raffinée du modèle de (). La partie infrarouge thermique a été développée par (). La condensation est paramétrée de façon différente pour la convection profonde et pour les autres nuages. Pour la convection nuageuse, on utilise la paramétrisation développée par (). Cette paramétrisation est basée sur une représentation en flux de masse d'une colonne convective nuageuse idéalisée. La colonne atmosphérique est divisée en trois sous-colonnes : une ascendance concentrée, une subsidence rapide associée à la pluie et enfin l'environnement, généralement en subsidence lente. Pour la partie non convective des nuages, on diagnostique une fraction nuageuse et un contenu en eau des nuages en se donnant a priori une fonction de distribution sous-maille de l'eau (, ). Les nuages ainsi diagnostiqués sont utilisés à la fois pour le calcul radiatif et pour calculer un taux de chauffage et un taux de précipitation. La pluie calculée ainsi peut se réévaporer dans les couches inférieures du modèle. L'eau vapeur et l'eau condensée sont ensuite advectées indépendamment dans la partie dynamique^2.8. Le transport de grande échelle de la vapeur d'eau et de l'eau liquide, traité à l'origine avec un schéma amont très diffusif sur l'horizontale et un schéma centré non positif sur la verticale, est maintenant calculé avec les schémas en volumes finis décrits un peu plus loin.

Dans les expériences présentées dans ce document, la température de surface de la mer est imposée. La conduction thermique dans les sols continentaux est quant à elle calculée au moyen d'un schéma multi-couches d'un sol homogène développé à l'origine pour la version martienne du modèle (, ). Si on note $\xi$ la conduction thermique du sol et $C$ la capacité thermique volumique, on peut montrer facilement que le flux conductif à la surface peut s'écrire sous la forme $F_c=-I\partial T/\partial z'$ et la conduction dans le sol $\partial T/\partial t=\partial^2 T/\partial {z'}^2$. $I=\sqrt{\xi C}$ est appelé inertie thermique et $z'$ est une pseudo profondeur définie par $z'=z\sqrt{C/\xi}$.

Suivant (), l'évaporation à la surface est calculée comme $E=1.35 \beta \rho {u_*}^2 [q_s(T_s)-q]$ où $u_*$ est la tension de vent en surface, $q$ est l'humidité spécifique de la première couche du modèle atmosphérique et $q_s(T_s)$ est l'humidité spécifique à saturation pour la température de surface $T_s$. Le coefficient d'aridité $\beta$ peut être soit imposé soit dépendre directement d'un contenu en eau du sol, calculé au fil du temps en fonction du bilan entre précipitation et évaporation, $d q_s/dt = P-E$. Le modèle utilisé par défaut, dit modèle du saut d'eau (bucket en anglais), fait l'hypothèse que $q_s$ peut croître jusqu'à une valeur de 150 mm d'eau. Entre 75 et 150 mm, le coefficient d'aridité $\beta$ est égal à 1 et l'évaporation est égale à l'évaporation potentielle. Ce coefficient varie linéairement entre 0 et 1 pour des contenus en eau du sol variant de 0 à 75 mm.

La partie physique du modèle est en fait en évolution constante. La description qui en est faite ci-dessus correspond à une vue instantanée de la version qui a servi à développer les aspects liés au transport des traceurs. Sauf mention explicite, c'est cette version de LMDZ3 qui est utilisée dans les simulations présentées dans ce document.

Figure 2.2: Climatologie des précipitations en janvier dans deux versions du modèle LMDZ. En haut : la version LMDZ3 avec le schéma de convection de Tiedtke et le modèle de seau d'eau en surface. Cette version a été beaucoup utilisée pour le développement des aspects chimie et traceurs. Au milieu : le modèle LMDZ4 avec le schéma de convection d'Emanuel, des nuages couplés à la convection et le schéma de surface SECHIBA du modèle ORCHIDEE. C'est la version utilisée pour le modèle couplé de l'IPSL. En bas : Observations (Climatologie de Xie-Arkin). Les simulations sont effectuées en prescrivant les températures de surface de l'océan et les cartes montrent des moyennes sur 5 années consécutives.

Figure 2.3: Climatologie des précipitations en juillet dans deux versions du modèle LMDZ. En haut : la version LMDZ3 avec le schéma de convection de Tiedtke et le modèle de seau d'eau en surface. Cette version a été beaucoup utilisée pour le développement des aspects chimie et traceurs. Au milieu : le modèle LMDZ4 avec le schéma de convection d'Emanuel, des nuages couplés à la convection et le schéma de surface SECHIBA du modèle ORCHIDEE. C'est la version utilisée pour le modèle couplé de l'IPSL. En bas : Observations (Climatologie de Xie-Arkin). Les simulations sont effectuées en prescrivant les températures de surface de l'océan et les cartes montrent des moyennes sur 5 années consécutives.

Depuis, notamment en vue du développement du modèle couplé IPSLCM4, la physique du modèle a évolué de façon significative. Les améliorations principales concernent :

1. le remplacement du schéma de convection nuageuse de Tiedtke par le schéma d'(). Comme celui de Tiedtke, le schéma d'Emanuel est basé sur une séparation de la colonne atmosphérique en trois compartiments. Les différences entre les deux schémas concernent principalement la fermeture (paramétrisation donnant par exemple le flux de masse à la base de l'ascendance adiabatique), l'échange d'air entre l'ascendance et l'environnement (décrit de façon plus sophistiquée dans le schéma d'Emanuel) ainsi que les descentes précipitantes (beaucoup plus actives dans le schéma d'Emanuel, cf. , ,). Un travail important est réalisé au LMD sur cette paramétrisation d'Emanuel (, ).
2. le schéma statistique de nuages qui a été couplé à la paramétrisation de la convection suivant l'approche de (). Un effort important a été consacré au réglage de la composante nuageuse.
3. le modèle ORCHIDEE qui est maintenant utilisé pour représenter la thermodynamique des surfaces continentales (, ,). Ce modèle distingue pour l'eau du sol un réservoir de surface et un réservoir profond. Les paramètres thermodynamiques de la surface sont prescrits en fonction de 10 types de végétation. Le routage de l'eau par les grands bassins fluviaux, depuis les zones de ruissellement jusqu'aux océans, est représenté. Différents degrés de complexité peuvent ensuite être activés, afin d'introduire par exemple une représentation explicite de la croissance des feuilles voir, à terme, un calcul évolutif de la végétation.
4. la prise en compte de 4 sous-surfaces avec des mailles fractionnaires pour les terres, les glaciers, les océans et la banquise. Ces 4 sous-surfaces sont en fait couplées à 4 sous-colonnes d'atmosphère sur lesquelles on effectue indépendamment 4 calculs de diffusion turbulente dans la couche limite.

Le modèle ainsi modifié - version LMDZ4 - couplé au modèle global d'océan et de banquise ORCALIM, est utilisé notamment pour effectuer des prévisions du changement climatique.

On ne s'appesantit pas dans ce document sur la climatologie du modèle d'atmosphère. On montre cependant une comparaison de la version précédente LMDZ3 du modèle et du nouveau modèle LMDZ4 en ce qui concerne la précipitation (Fig. 2.2 et Fig. 2.3). Les différences principales entre les deux versions proviennent en fait du changement de schéma de convection. L'ancien modèle avait tendance à exagérer l'intensité des précipitations tropicales dans les zones de convergences, notamment à l'Est de Madagascar ainsi que sur l'océan Pacifique. Ce biais est nettement diminué dans le nouveau modèle. La diminution (amélioration) des précipitations sur les continents nord en été (Sibérie, Canada), est elle davantage due au schéma de surface. La détérioration la plus notable avec le nouveau modèle est sans doute la distribution des précipitations de mousson autour de l'Inde en juillet. Lorsqu'on interprète des résultats relatifs à la simulation des concentrations atmosphériques d'aérosols ou d'espèces chimiques, il faut bien sûr garder à l'esprit ces différences quant à la capacité du modèle à simuler correctement le climat.

Organisation informatique - version unidimensionnelle - mode guidé

Figure 2.4: Séparation entre les parties dynamique 3D et physique 1D du code LMDZ.

Modularité

L'organisation du modèle LMDZ est très étroitement calquée sur cette distinction entre partie "dynamique" - la seule partie où soient pris en compte des échanges horizontaux entre des mailles du modèle - et la partie ``physique" qui peut être vue comme une juxtaposition de ``colonnes" d'atmosphère.

Cette spécificité de la partie physique est exploitée en ce sens que le codage de toutes les paramétrisations est fait avec un indice interne muet qui représente la grille horizontale. Cette organisation est illustrée sur la Fig. 2.4. Chacune des parties, physique et dynamique, a ses propres fichiers de conditions initiales et de sorties.

Cette écriture se prête à la fois à la vectorisation et à la parallélisation des codes. Cette approche permet également de disposer de façon transparente d'une version unidimensionnelle du modèle de circulation. Pour disposer d'un modèle unidimensionnel, il suffit d'écrire un programme dans lequel on initialise des profils météorologiques sur un point particulier du globe et dans lequel on appelle ensuite en boucle le moniteur physique.

Autre atout important, cette conception modulaire permet de gérer en parallèle des ``physiques" différentes interfacées avec le même code dynamique. Ce point est essentiel pour les études menées au LMD sur Mars et sur Titan. Enfin, on peut noter qu'il existe aussi une version bidimensionnelle, latitude-altitude, du noyau dynamique du modèle qui a été abondamment utilisée sur Titan comme on le voit à la fin de ce document.

Mode climatique et mode guidé

Le modèle LMDZ est avant tout développé pour des études climatiques. On effectue alors des intégrations longues, dans lesquelles l'état initial est vite ``oublié" et les résultats ne s'interprètent qu'en termes statistiques. Cependant, notamment pour des aspects de validation, il peut s'avérer utile de contraindre le modèle à suivre une situation météorologique observée. Dans ce mode ``guidé" (on parle souvent de ``nudging" en anglais), les champs du modèle sont rappelés avec un terme linéaire vers les champs des analyses ou des réanalyses

$$\frac{\partial X}{\partial t}=M(X) + \frac{X_a-X}{\tau}$$ (2.22)

où $X$ est une variable météorologique ($u$, $v$, $\theta$, $p_s$) et $X_a$ est la variable correspondante issue des analyses et interpolée sur la grille spatio-temporelle du modèle. $M$ représente le calcul des dérivées temporelles par le modèle. Les temps de relaxation $\tau$ peuvent être différents pour les différentes variables ou dépendre de la région considérée. Cette approche peut être vue comme une alternative légère à l'assimilation directe de données météorologiques dans le modèle telle qu'elle est mise en uvre dans les centres de prévision opérationnelle (, ).

Une approche intermédiaire a été mise en uvre dans LMDZ (, ). Elle consiste à minimiser une fonction mesurant la distance entre l'état du modèle et les analyses. Cette approche, qui s'apparente encore davantage à l'assimilation météorologique opérationnelle, permet d'introduire, en plus des analyses, des observations supplémentaires ou des contraintes comme des pénalités sur les modes de gravité excités par la procédure d'assimilation. Cette dernière méthode n'est pas utilisée ici.

L'équation de transport : séparation d'échelles

Séparation des processus

Dans les modèles de chimie-climat comme dans les modèles débranchés transport-chimie, les processus comme la microphysique des aérosols ou les réactions chimiques d'une part et le transport de l'autre sont généralement traités alternativement et séquentiellement. On parle en anglais d'``operator splitting". La partie transport proprement dite peut alors être traitée de façon systématique, indépendamment de l'espèce trace considérée, en assimilant cette espèce à un traceur conservé^2.9.

Dans cette partie, on s'intéresse donc à la modélisation du transport d'un traceur conservatif (une espèce trace suivant exactement l'air) et passif (n'affectant pas en retour la météorologie). L'équation du transport pour un traceur de ce type est simplement

$$\frac{d c}{d t}=0$$ (2.23)

ou

$$\frac{\partial c}{\partial t}+{\bf v}{\mbox{\bf grad}}c=0$$ (2.24)

Dans la présentation qui suit, on suppose le champ de vent ${\bf v}$ connu.

Dans les applications atmosphériques classiques, l'équation du transport ne peut être résolue jusqu'à l'échelle de la diffusion moléculaire. Cette constatation s'applique d'ailleurs aussi bien à l'observation qu'à la modélisation. Dans les deux cas, on travaille explicitement jusqu'à une échelle donnée mais on traite la petite échelle de façon statistique. Dans le cas d'un modèle numérique, la grande échelle est définie en pratique par le maillage (ou par la troncature pour les modèles spectraux). L'effet des grandeurs sous-mailles sur les variables de grande échelle ne peut être représenté que de façon statistique, au travers de paramétrisations.

Séparation d'échelles

Pour séparer échelle explicite et échelle turbulente, on introduit la notion de moyenne d'ensemble. La turbulence est considérée comme un processus aléatoire. Un élément du processus correspond à une réalisation complète de l'écoulement atmosphérique. La moyenne d'ensemble d'une grandeur $X$, qu'on notera $\overline{X}$, est simplement l'espérance mathématique de cette variable.^2.10

Pour introduire proprement le découpage pour un fluide compressible, il faut introduire en plus une moyenne pondérée par l'air $\tilde{X}=\overline{\rho X}/\overline{\rho}$. La fluctuation turbulente par rapport à cette moyenne $X'=X-\tilde{X}$ obéit à l'identité $\overline{\rho X'}=\tilde{X'}\overline{\rho}=0$.

Si on note ${\bf v}$ le champ de vent et $c$ la concentration massique d'un traceur conservatif ($dc/dt=0$), l'équation de transport non visqueux peut s'écrire soit

$$\frac{\partial c}{\partial t}+{\bf v}{\mbox{\bf grad}}{c}=0$$ (2.25)

sous sa forme advective, soit, en introduisant l'équation de continuité pour le fluide atmosphérique,

$$\frac{\partial \rho}{\partial t}+{\mbox{div}}\left(\rho{\bf v}\right)=0,$$ (2.26)

sous sa forme conservative ou flux

$$\frac{\partial \rho c}{\partial t}+{\mbox{div}}\left(\rho{\bf v}c\right)=0$$ (2.27)

En prenant la moyenne d'ensemble de la forme conservative et en remarquant que

$$\overline{\rho {\bf v}c}=\overline{\rho}\tilde{{\bf v}}\tilde{c}+\overline{\rho {\bf v}'c'},$$ (2.28)

on obtient l'équation

$$\frac{\partial \overline{\rho}\tilde{c}}{\partial t} +{\mbox... ...ight) + {\mbox{div}}\left(\overline{\rho{\bf v}'c'}\right) = 0$$ (2.29)

En remarquant que l'équation de continuité pour l'air est inchangée par la prise de moyenne, on peut repasser à la forme advective. En adoptant les notations $\rho$, $u$, $v$ et $c$ pour les variables grande échelle, on obtient finalement

$$\frac{\partial c}{\partial t}+{{\bf v}}.{\mbox{\bf grad}} {... ...{1}{\rho}{\mbox{div}}\left(\overline{\rho{\bf v}'c'}\right) =0$$ (2.30)

On voit donc qu'on retrouve, pour les moyennes pondérées par l'air, les équations initiales avec des termes supplémentaires liés aux corrélations entre les fluctuations turbulentes de ${\bf v}$ et de $c$.

Ci-dessous, on présente en détail les schémas introduits pour traiter du transport de grande échelle ainsi que les paramétrisations des termes turbulents calquées sur les paramétrisations d'origine du code LMDZ. Dans le chapitre suivant, on présente en détail un travail spécifique mené sur le transport turbulent dans la couche limite convective.

Le transport grande échelle

Les différentes approches

De nombreuses méthodes ont été développées au cours des dernières décennies pour représenter l'advection, c'est à dire les parties non turbulentes des Eqs 2.30 ou 2.29 formellement équivalentes aux Eqs 2.25 et 2.27.

On distingue notamment :

• Les méthodes lagrangiennes: elles consistent à ``ensemencer" le domaine d'étude (où la région source) avec des particules et à les suivre individuellement. Ces méthodes peuvent être très précises mais il faut généralement utiliser un grand nombre de particules, notamment pour éviter que des trous se créent dans la distribution spatiale de ces particules.
• Les méthodes semi-Lagrangiennes: (ou méthodes des caractéristiques dans le monde des sciences de l'ingénieur) on suppose qu'à un pas de temps $t$ donné, $q$ est connu sur une grille. On prédit la valeur sur le même maillage au temps $t+\delta t$ en remontant à $t$ ``le long du vent": $q({\bf x},t+\delta t)=q({\bf x}-\delta t {\bf v}, t)$. Le problème est que, si le point ${\bf x}$ se trouve sur le maillage, le point ${\bf x}-\delta t{\bf v}$ a peu de chances de s'y trouver. On interpole donc cette valeur à partir des valeurs sur le maillage au temps $t$. Cette méthode évite le problème principal des méthodes lagrangiennes. Sa précision est déterminée par la précision des routines d'interpolation. Les méthodes semi-Lagrangiennes ont été très largement utilisées dans la communauté atmosphérique. Elles présentent un défaut important: il est difficile de garantir la conservation de la quantité totale de traceur.
• L'advection des contours: elle consiste à advecter des particules uniquement sur des ``iso-niveaux" du champs considéré. Cette dernière méthode, potentiellement très précise, se couple difficilement à d'autres processus (il faut alors calculer l'effet des autres processus sur la topologie des ``iso-niveaux").
• Les différences finies : il existe également tout un tas de méthodes en différences finies (cf. par exemple , ). Les différences finies permettent de dériver des schémas à des ordres élevés mais ne garantissent pas à priori la conservation du traceur ou la positivité.
• Les volumes finis : ils sont de plus en plus utilisés dans la communauté. Basés sur l'écriture de bilans sur des volumes de contrôle, ils sont a priori appropriés pour la modélisation de processus physiques couplés et complexes pour laquelle le respect des lois de conservation s'avère souvent essentiel. C'est cette approche qui a été suivie pour introduire les traceurs dans LMDZ.
• Les éléments finis : basés sur une décomposition sur des fonctions de bases analytiques, ils sont encore relativement peu utilisés dans la communauté.

Les schémas en volumes finis

^2.11 L'intégration de l'équation de transport (2.27) et de l'équation de conservation de la masse (2.26) sur un polyèdre (volume de contrôle) à $N$ faces conduit aux formulations:

$$\frac{\partial m}{\partial t}=\sum_{n=1}^N {{\cal U}}_n$$ (2.31)

et

$$\frac{\partial cm}{\partial t}=\sum_{n=1}^N {{\cal F}}_n$$ (2.32)

où $m$ est la masse totale d'air dans le polyèdre, $c$ est la concentration massique moyenne du traceur dans le volume, ${{\cal U}}_n$ et ${{\cal F}}_n$ sont les flux de masse et flux de traceur vers l'intérieur du polyèdre pour la face $n$.

Les formulations en ``volumes finis" présentent l'avantage très important d'être conservatives par nature, dès lors que le même flux est utilisé pour les volumes en amont et en aval de l'interface considérée.

Pour le plus connu de ces schémas, introduit par (), ${{\cal F}}$ est simplement estimé comme le produit de ${\cal U}$ par la valeur de $c$ dans le volume amont (dans la maille d'où l'air provient). Ce schéma simple, souvent appelé schéma amont, garantit

1. la conservation de la quantité totale de traceur,
2. la positivité (un traceur positif partout reste positif),
3. la monotonie (pour l'advection en dimension 1, une distribution monotone reste monotone),
4. la ``décroissance de la variation totale" du champ de traceur, c'est à dire la décroissance de la somme des écarts absolus entre deux concentrations successives. Cette propriété garantit la stabilité numérique du schéma.
5. l'invariance par addition d'une constante au champ de traceur (le modèle de Godunov est en fait linéaire par rapport au champ de traceurs, ce qui garantit en plus la symétrie temporelle du transport comme on l'explique au Chapitre ).

Cependant, ces propriétés physiques fondamentales ne sont obtenues dans le schéma amont qu'au prix d'une très grande diffusion numérique.

() a proposé d'utiliser, pour la valeur amont de $c$, non pas la valeur moyenne dans la maille amont mais une valeur extrapolée à la frontière de celle-ci, en utilisant une approximation polynomiale de la distribution sous-maille du traceur dans la maille amont. Deux des schémas proposés par Van Leer en 1977 ont été introduits apparemment indépendamment dans la littérature météorologique par () - le schéma III dans la classification de Van Leer - et () - le schéma VI -. Alors que dans son article original, Van Leer conclut que la complexité du schéma VI est trop grande par rapport au gain en précision, ce schéma est devenu une référence dans la communauté météorologique.

Les principes de dérivation de ces différents schémas sont exposés ci-dessous.

Description des schémas en dimension un

En dimension 1 et après intégration sur un pas de temps, les Eq. 2.31 et 2.32 s'écrivent simplement

$$\delta_t (m)_i=U_{i-1/2}-U_{i+1/2}$$ (2.33)

et

$$\delta_t (cm)_i=F_{i-1/2}-F_{i+1/2}$$ (2.34)

où les indices correspondent à la position sur l'axe des $x$, $\delta_t$ représente la différence finie en temps et $U_{i+1/2}$ et $F_{i+1/2}$ sont les transferts de masse d'air et de traceur à travers l'interface $i+1/2$ durant un pas de temps.

Si on connaît complètement la distribution spatio-temporelle du vent, de la densité de l'air et de la concentration de traceur, on a

$$m_i= \int_{\Delta x} \rho(x,t) dx$$ (2.35)

$$c_i= \int_{\Delta x} c\rho(x,t) dx/\int_{\Delta x} \rho(x,t) dx$$ (2.36)

$$U_{i+1/2}= \int_{t}^{t+\Delta t} \rho(x_{i+1/2},t) u(x_{i+1/2},t) dt$$ (2.37)

et

$$F_{i+1/2}= \int_{t}^{t+\Delta t} \rho(x_{i+1/2},t) c(x_{i+1/2},t) u(x_{i+1/2},t) dt$$ (2.38)

Dans la présentation faite ici, on considère que les transferts de masse $U_{i+1/2}$ sont connus à chaque instant et sur tout le maillage. L'évolution temporelle de la masse d'air $m_i$ est donc également complètement déterminée. La seule chose qui reste à estimer est le flux de traceur.

De façon générale, ce flux $F_{i+1/2}$ peut être écrit comme le produit du transfert de masse $U_{i+1/2}$ par la valeur moyenne de la concentration du traceur ${\breve{c}}_{i+1/2}$ dans l'air qui traverse l'interface au cours du pas de temps.

Figure 2.5: Principe du schéma de Van Leer et notations. On montre le cas où la distribution sous-maille est représentée au moyen d'un polynôme du premier degré. L'axe vertical correspond à la concentration massique du traceur. Les deux axes horizontaux correspondent à l'indiçage des variables et à la masse d'air comptée à partir du centre de la maille $i$ et normalisée par la masse totale de cette même maille, $m_i$. La surface grisée correspond à la quantité de traceur qui est transféré au travers de l'interface durant un pas de temps.

La méthode proposée par Van Leer consiste à approximer la distribution sous-maille par un polynôme pour lequel le calcul de ${\breve{c}}_{i+1/2}$ - et donc du transfert de traceur - peut être fait exactement. Le principe et les notations sont illustrées sur la Fig. 2.5 pour le cas d'un polynôme du premier degré. Avec ces notations, les valeurs moyennes de $c$ et $m$ après un pas de temps d'advection s'écrivent simplement

$$c^*_i=\left(c_im_i+U_{i-1/2}{\breve{c}}_{i-1/2}-U_{i+1/2}{\breve{c}}_{i+1/2}\right)/m^*_i$$ (2.39)

et

$$m^*_i=m_i+U_{i-1/2}-U_{i+1/2}$$ (2.40)

Pour des valeurs positives de $U_{i-1/2}$ et $U_{i+1/2}$, l'Eq. 2.39 peut se récrire

$$c^*_i=\left[U_{i-1/2}{\breve{c}}_{i-1/2}+\left(m_i-U_{i+1/2}\right){\breve{c}}_i\right]/m^*_i$$ (2.41)

où ${\breve{c}}_i$ est la concentration massique moyenne dans l'air qui reste dans le volume $i$ pendant le pas de temps (cf. Fig. 2.5). Avec ces notations, le schéma de Van Leer peut être décrit comme la séquence d'opérations suivante : on commence tout d'abord par définir une nouvelle distribution ${\breve{c}}$ (les croix sur la Fig. 2.5) à partir de la distribution initiale $c$ (les points) en utilisant une approximation polynomiale : on advecte ensuite cette nouvelle distribution approchée à travers les mailles.

Un des points essentiels de l'article original de () est de remarquer qu'on peut garantir la monotonie du schéma d'advection - en dimension 1, une distribution monotone reste monotone après advection - en imposant que ${\breve{c}}$ prenne une valeur intermédiaire entre les deux valeurs voisines de $c$ dans les régions où $c$ est monotone.

Supposons pour fixer les idées que $c_{i-1}\leq c_i \leq c_{i+1}$ et $U\geq 0$. La condition ci-dessus s'écrit simplement

$$c_j \leq {\breve{c}}_{j+1/2} \leq c_{j+1} \mbox{ pour } j=i-1 \mbox{ et } j=i$$ (2.42)

et

$$c_{i-1} \leq {\breve{c}}_i \leq c_i.$$ (2.43)

En introduisant la condition 2.42 dans l'Eq. 2.39, on assure que $c^*_i \leq c_{i}$ tandis que l'Eq. 2.41 et les conditions 2.42 et 2.43 garantissent que $c_{i-1}\leq c^*_i$. On voit finalement que la condition suffisante énoncée ci-dessus garantit que $c^*_i$ est compris entre $c_{i-1}$ et $c_i$ (resp. $c_i$ et $c_{i+1}$) pour $U>0$ (resp. $U<0$). Or, comme le remarque (), cette proposition implique la monotonie.

Si on rajoute en plus la condition que ${\breve{c}}_i=c_i$ (distribution sous-maille constante) quand $c_i$ est un extremum, on interdit la croissance des extrema. Ceci implique de fait la positivité du schéma et interdit la création d'oscillations provenant du caractère dispersif du schéma numérique (provenant de l'advection avec des vitesses différentes des différentes composantes de fourier de la distribution, cf. e. g. , ).

Noter aussi qu'avec les définitions ci-dessus, une distribution uniforme de traceur sera inchangée par l'advection, même avec des champs de vents divergents (on s'en convainc en remplaçant ${\breve{c}}_{i-1/2}$ et ${\breve{c}}_{i+1/2}$ par une valeur constante dans l'Eq. 2.39).

Dans les dérivations ci-dessus, on a supposé implicitement que le pas de temps était suffisamment petit pour éviter de transférer plus que la maille d'une cellule en un pas de temps (nombre de Courant $U/m$ inférieur à 1). Si le flux de masse est exactement égal à la masse d'air dans le volume amont, on calculera la distribution de traceur exactement, quelque soit l'approximation choisie pour la distribution sous-maille.

Le schéma de Godunov

La première approximation consiste à supposer que le traceur est constant dans chaque maille (polynôme de degré zéro). ${\breve{c}}$ est alors simplement la valeur de $c$ dans la maille amont ( ${\breve{c}}_{i+1/2}=c_i$ si $U_{i+1/2}>0$ et $c_{i+1}$ sinon). Ce schéma proposé à l'origine par () est bon marché, positif et monotone mais au prix d'une diffusion numérique très forte. La diffusivité du schéma peut d'ailleurs être quantifiée. Dans le cas d'un champ de vent non divergent ( $u=U\delta x/(m\delta t)=$cste à une dimension), le schéma 2.39 s'écrit

$$c^*_i=c_i +\frac{U}{m}c_{i-1}-\frac{U}{m}c_{i}$$ (2.44)

et peut se récrire

$$\frac{c^*_i-c_i}{\delta t} = u\frac{c_{i-1}-c_{i}}{\delta x}$$ (2.45)

$$= u \frac{c_{i-1}-c_{i+1}}{2\delta x}+\frac{u\delta x}{2} \frac{c_{i-1}-2c_i+c_{i+1}}{{\delta x}^2}$$ (2.46)

comme la somme d'un schéma numérique centré d'ordre 2 et d'un terme de diffusion numérique avec une diffusivité $u\delta x/2$. ^2.12

Schémas du second ordre

Pour passer à un ordre supérieur, on suppose que la distribution sous-maille est linéaire avec une pente ${\left(\delta c\right)}_i$ donnant les valeurs aux bords gauche et droit de la maille $i$ comme $c_i^\pm=c_i\pm{\left(\delta c\right)}_i/2$ (voir l'illustration de la Fig. 2.5).

Dans ce cas, ${\breve{c}}_{i+1/2}$ est donné par

$${\breve{c}}_{i+1/2} = c_i + \frac{1}{2}\left(1-\frac{U_{i+1/2}}{m_{i}}\right) {\left(\delta c\right)}_i \mbox{ si } U_{i+1/2}>0$$ (2.47)

$$= c_{i+1} - \frac{1}{2}\left(1+\frac{U_{i+1/2}}{m_{i+1}}\right) {\left(\delta c\right)}_{i+1} \mbox{ , sinon.}$$ (2.48)

Différents schémas ont été proposés (, ,) correspondant à différentes estimations de ${\left(\delta c\right)}$, dont deux schémas particulièrement intéressants, l'un pour son faible coût et l'autre pour sa précision.

Figure 2.6: Schémas du 2nd ordre. Illustration de l'estimation de la pente par différences finies pour le schéma I de Van Leer (a) et par calcul par moindres carrés à partir de la distribution en ligne brisée résultant de l'advection au pas précédent (b). Cette seconde estimation correspond au schéma des pentes de () ou au schéma III de Van Leer (se reporter au texte pour plus de détails). c : Exemples de calculs d'advection unidimensionnelle sur un domaine périodique de 70 points (axe horizontal), avec une vitesse constante $u$ et trois distributions initiales. Les concentrations de traceur ($c$, unité arbitraire) sont montrées après une révolution complète, au pas de temps 350 pour un nombre de Courant $U/m=0.2$.

Dans le schéma I de Van Leer, la pente est simplement estimée à chaque pas de temps par différences finies ( $\delta c_i=\left(c_{i+1}-c_{i-1}\right)/2$) comme l'illustre le graphique a de la Fig. 2.6. Ce schéma peut être considéré comme une version volumes-finis du schéma de ().

Dans le second schéma, la pente au pas de temps $t+\delta t$ dans une maille donnée est calculée à partir de la distribution en ligne brisée résultant de l'advection au temps $t$ (illustration sur le graphique b de la Fig. 2.6). La nouvelle distribution sous-maille minimise la distance quadratique par rapport à cette distribution. Ce second schéma (le schéma III dans l'article original de Van Leer) a en fait été redécouvert quelques années plus tard par ()^2.13 dans un contexte météorologique tridimensionnel. Ce schéma a été popularisé dans le modèle de circulation générale du NASA/GISS. On l'appelle souvent schéma des pentes (``slopes scheme" en anglais). On retiendra ce nom dans ce qui suit.

On présente, sur le graphique c de la Fig. 2.6, un exemple typique de calcul d'advection unidimensionnel avec le schéma de Godunov, Van Leer I et le schéma des pentes. Le schéma des pentes est bien sûr le plus précis mais la différence principale se situe entre le schéma de Godunov, extrêmement diffusif, et les deux autres. On voit également apparaître pour les deux schémas d'ordre 2 des oscillations provenant de la non monotonie du schéma et la création de valeurs négatives.

Limiteurs de pentes

L'un des intérêts principaux du travail de () réside dans la possibilité qu'il offre d'assurer facilement la monotonie des schémas. Une première façon consiste à appliquer directement la contrainte explicitée à la fin de la Section 2.3.3.

Pour les schémas du second ordre, () propose aussi une condition suffisante à la fois plus brutale et plus simple. Il suffit, pour que le schéma soit monotone, d'imposer que la distribution dans une maille $i$ soit entièrement comprise entre les valeurs moyennes des deux mailles adjacentes et que la pente soit du même signe que dans ces deux mailles (à noter que ce dernier critère est automatiquement satisfait par le schéma I). Cette condition suffisante peut s'exprimer facilement comme un limiteur de pente. Pour le schéma I, la formulation complète du calcul de la pente avec limiteur s'écrit simplement

$${\left(\delta c\right)}_i = \mbox{sign}\left(c_{i+1}-c_i\right)\times$$ (2.49)

$$\mbox{min}\left(\frac{\vert c_{i+1}-c_{i-1}\vert}{2},2\vert c_{i+1}-c_i\vert,2\vert c_i-c_{i-1}\vert\right)$$

si $c_i$ est compris entre $c_{i-1}$ et $c_{i+1}$ et 0 sinon.

Figure 2.7: a) Illustration de l'application d'un limiteur de pente. A droite, impact des limiteurs de pentes sur le schéma des pentes (b) et sur le schéma I de Van Leer (c). Les cas sans limiteurs (les mêmes que sur la Fig. 2.6) et les limiteurs faibles et forts (se reporter au texte pour plus de détails) sont montrés pour les distributions carrée et gaussienne.

L'effet de l'application d'un limiteur fort (décrit ici) ou d'un limiteur faible (application directe des Eq. 2.42 et Eq. 2.43), qui correspondent respectivement aux Eq. (66) et (74) données par () est illustré sur la Fig. 2.7a. Un test numérique de l'impact de ces limiteurs est montré sur les Fig. 2.7c et d. Le limiteur fort dégrade de façon significative la précision du schéma des pentes tandis que le limiteur faible corrige ce schéma de façon sélective au niveau des oscillations, sans émousser par exemple le sommet de la gaussienne. Pour le schéma I, la différence entre les deux limiteurs est nettement plus marginale.

Il existe une alternative élégante à l'Eq. 2.49 qui consiste à utiliser la moyenne géométrique des deux pentes voisines

$${\left(\delta c\right)}_i = \frac{2\left(c_{i+1}-c_i\right)\left(c_i-c_{i-1}\right)}{c_{i+1}-c_{i-1}}$$ (2.50)

quand $\left(c_{i+1}-c_i\right)\left(c_i-c_{i-1}\right)>0$ et 0 sinon. Cette formulation satisfait automatiquement la limitation forte. Ce schéma est seulement légèrement plus diffusif que le schéma I mais le gain en temps est également relativement faible, le calcul de la pente ne représentant qu'une partie du coût du schéma.

Schémas du troisième ordre

La description de la distribution sous-maille du traceur peut être encore améliorée en utilisant un polynôme du second ordre. L'équivalent du schéma I de Van Leer (le schéma IV dans l'article original) consiste à évaluer les coefficients de ce polynôme par différences finies à partir des valeurs moyennes dans les mailles. Van Leer a attaché peu d'importance à ce schéma qui présente par exemple comme défaut de n'être pas plus précis que le schéma I pour un nombre de Courant de 0.5. Cependant, pour des valeurs plus petites, il peut devenir nettement plus précis. Comme le schéma I, ce schéma présente l'intérêt de n'utiliser qu'un traceur par point de grille. () ont développé une alternative au schéma IV en utilisant un filtre non-linéaire assez élaboré permettant de renforcer les pentes en cas de choc (se reporter également à , ,). La méthode qui en résulte a été baptisée Piecewise Parabolic Method. Elle a été utilisée pour les applications atmosphériques par exemple par (), () ou (). Dans ce manuscrit on ne présente pas les résultats avec PPM. Mais des tests ont été réalisés récemment avec ce schéma qui se comporte de façon assez similaire au schéma des pentes pour un coût numérique intermédiaire entre les schémas I et III (pentes) en termes de temps de calcul. Le coût en stockage est le même que pour le schéma I.

L'équivalent à l'ordre 3 du schéma III de Van Leer est le schéma VI connu dans la littérature météorologique sous le nom de schéma de (). Le schéma de Prather est évidemment beaucoup plus précis que tous les schémas présentés jusque là comme le montrent les illustrations de la section suivante. Cependant, il nécessite la conservation de 10 traceurs indépendants : la moyenne, les 3 pentes dans les 3 directions d'espace $c_x$, $c_y$, $c_z$ et les 6 moments du second ordre $c_{xx}$, $c_{yy}$, $c_{zz}$, $c_{xy}$, $c_{yz}$ et $c_{xz}$. Il est également nécessaire de lui adjoindre des algorithmes du type limiteurs de pente pour éviter complètement les oscillations numériques ou les valeurs négatives.

Introduction des schémas dans LMDZ

L'introduction des schémas en volumes finis est relativement facile dans le modèle de circulation générale LMDZ car les flux de masse sont déjà définis sur une grille décalée. Ces flux de masse sont ceux utilisés dans la partie météorologique pour intégrer l'équation de continuité pour l'air, à savoir l'Eq. 2.33 ou l'Eq. 2.88 donnée plus loin.

Flux alternés

Figure 2.8: Advection uniforme d'un pic gaussien de concentration le long de la diagonale d'un maillage bidimensionnel régulier. En haut, le calcul de la divergence des flux est effectué à partir de flux estimés indépendamment en $x$ et en $y$. La ligne du bas montre un calcul alterné avec d'abord une advection en $x$ puis une advection en $y$. Ces illustrations numériques ont été réalisées avec le schéma I de Van Leer, sur un maillage de 60 points dans chaque direction horizontale.

Les présentations ci-dessus ont été faites avec une seule dimension d'espace. Une solution naturelle pour passer en 3 dimensions consiste à calculer d'abord les flux dans chaque dimension puis à calculer la divergence des flux.

Une alternative classique consiste à faire trois calculs d'advection, successivement dans les trois directions. On parlera de flux alternés ou de directions alternées.

Cette seconde méthode est curieusement plus précise dans le cas d'une advection en diagonale par rapport au maillage. Ceci est illustré sur la Fig. 2.8 pour une configuration bidimensionnelle. Si on considère un pas de temps d'advection, le calcul direct ne tient pas compte de l'air passant directement dans la maille située en diagonale par rapport à la maille du milieu (le petit carré sombre sur les figures). En revanche, cet air est compté deux fois, une fois dans l'advection en $x$ et une fois dans l'advection en $y$. Ceci produit une très forte diffusion latéral en réduisant la propagation dans la direction du mouvement. En partant d'un panache gaussien, la distribution obtenue s'étire dans la direction transverse à l'écoulement (en haut à droite sur la figure).

Au contraire, pour l'advection en flux alternés, le carré sombre commence par passer dans la maille située à droite de la maille d'origine puis en haut et arrive bien finalement dans la maille située en diagonale par rapport à la maille d'origine. Si on pouvait à chaque instant connaître exactement la distribution sous-maille du traceur, on aurait donc un calcul exact.

Dans la méthode des flux alternés, on intègre successivement non seulement l'équation d'advection du traceur (Eq. 2.39) mais aussi l'advection de transport de l'air (Eq. 2.40)^2.14. Avec cette approche, une distribution uniforme de traceur est inchangée par l'advection indépendamment du caractère divergent ou non du champ de vent (, , montrent que ce n'est pas garanti par toutes les formulations en flux alternés).

Cette approche en flux alternés est utilisée assez systématiquement dans les calculs d'advection en volumes finis (, ,,,).

On utilise ici la séquence proposée par () :

Direction	pas de temps
X	$\delta t$/2
Y	$\delta t$/2
Z	$\delta t$
Y	$\delta t$/2
X	$\delta t$/2

Sauter des mailles

Si on veut ne travailler qu'avec des nombres de Courant ( $\delta t u/\delta x$ ou $U/m$) inférieurs à 1, le raffinement de la discrétisation près des pôles, inhérent aux grilles globales longitudes-latitudes, nécessite l'utilisation de pas de temps extrêmement petits dans la direction longitudinale. Une possibilité pour résoudre ce problème consiste à découper encore davantage le pas de temps pour l'advection longitudinale.

Il existe une alternative utilisant le fait que le calcul est exact pour un nombre de courant de 1. Dans ce cas en effet, le flux est simplement le produit du flux de masse par la concentration moyenne du traceur dans la maille amont. Si le nombre de courant est plus grand que 1, (pour fixer les idées, si on prend $U_{i+1/2}>0$ et $m_i<U_{i+1/2}< m_i+m_{i-1}$) on peut calculer le transfert de traceur comme la somme de la quantité de traceur dans la maille amont ($m_i c_i$) et de la quantité de traceur transférée depuis la maille $i-1$ avec un flux de masse $U_{i+1/2}-m_i$.

De façon générale, pour

$$\sum_{j=i-n+1}^i m_j < U_{i+1/2}\leq \sum_{j=i-n}^i m_j$$ (2.51)

on prend

$$F_{i+1/2}=\sum_{j=i-n+1}^i m_j c_j + \left(U_{i+1/2}-\sum_{j=i-n+1}^i m_j\right) {\breve{c}}_{i-n+1/2}$$ (2.52)

avec

$${\breve{c}}_{i-n+1/2}=c_{i-n} +\frac{1}{2}\left(1-\frac{U_{i... ...sum_{j=i-n+1}^i m_j}{m_{i-n}}\right) {\left(\delta c\right)}_i$$ (2.53)

Cette façon de voir les schémas en volumes finis présente de fortes analogies avec les approches semi-Lagrangiennes pour lesquelles ont prédit la valeur ponctuelle en un point du maillage en remontant en amont la trajectoire de la particule. Cette analogie a conduit () à introduire la notion de flux-form Semi-Lagrangian method. A noter qu'en parallèle, les adeptes des schémas semi-Lagrangiens classiques modifient leurs formulations pour garantir la conservation et les font ressembler de plus en plus à des schémas en volumes finis (e. g. , ).

On pourrait de façon générale et pour un pas de temps aussi grand qu'on veut, estimer les contours du maillage transportés à rebours sur un pas de temps. A partir d'une estimation polynomiale des distributions sous mailles, comme celles proposées par Van Leer, on pourrait calculer alors la distribution de traceur sur ce maillage déformé par le transport à rebours. Mais on paie vite en complexité du schéma et en coût informatique (introduction de branchements supplémentaires et gestion de la mémoire, particulièrement pénalisant sur des ordinateurs vectoriels ou parallèles) ce qu'on gagne avec l'utilisation d'un pas de temps plus long. De plus, on ne tient pas forcément à utiliser des pas de temps plus longs que les constantes caractéristiques des autres processus (transport turbulent dans la couche limite, chimie, etc...). On restreint donc ici le traitement des nombres de Courant plus grands que 1 à la direction longitudinale pour traiter spécifiquement le problème de raffinement du maillage en longitude à l'approche des pôles. Pour une résolution typique du modèle LMDZT, avec une grille d'environ 100 000 points et un pas de temps de 15 minutes, la condition $U/m>1$ n'est rencontrée typiquement que quelques dizaines de fois par pas de temps ce qui rend le coût marginal raisonnable.

Le point du pôle

Dans le modèle du LMD, les pôles correspondent à des centres de mailles triangulaires. La dimension de ces mailles est 2 fois plus petite dans la direction méridienne que celle des mailles normales. Toutes ces mailles ne sont en fait pas indépendantes et on impose que les valeurs scalaires (pression de surface, température, concentration de traceurs) soient toutes identiques. L'évolution de la masse d'air ou de la quantité d'un traceur au pôle est estimée à partir de la convergence totale des flux méridiens. En fait, on peut considérer le pôle comme un volume de contrôle consistant en un polygone avec autant d'arêtes que le nombre de points longitudinaux de la grille.

Il semble impossible de garantir strictement la monotonie du schéma en conservant des pentes non nulles pour ces mailles polaires. On retient donc aux pôles un schéma de Godunov.

Advection transpolaire

La Fig. 2.9 montre un exemple d'advection transpolaire d'une distribution sinusoïdale avec un écoulement en rotation solide le long du méridien de Greenwich. Les figures du bas montrent le résultat de l'advection après une révolution complète (la solution exacte, à gauche, est identique à la distribution initiale). La distribution obtenue avec le schéma I (cas a) est allongée dans la direction méridienne à cause de la diffusion numérique, dans la direction de l'écoulement. Pour la même résolution, le schéma de Prather est beaucoup moins diffusif. La forme des iso-lignes est cependant légèrement altérée par le passage au pôle. Pour ces calculs, on s'est arrangé pour que le nombre de Courant soit toujours inférieur à 1. Le dernier exemple sur la droite (cas b) correspond à un calcul avec le schéma I utilisant un pas de temps beaucoup plus grand, avec un nombre de Courant qui atteint 8 dans la direction longitudinale près du pôle. On remarque que les résultats sont plutôt meilleurs que ceux du cas a à cause du plus petit nombre de pas de temps nécessaire (160 au lieu de 16000).

Figure: Test numérique d'advection transpolaire. Suivant () et (), la distribution initiale de traceur est donnée par $c(\lambda,\phi)=(1+\cos[\mbox{min}(r[\lambda,\phi]/R,1)])/2$ avec $r=\arccos(\cos\lambda\cos\phi)$ et $R=7\times (2\pi )/128$. La grille utilisée comprend 128 points en longitude et 64 en latitude. On montre de gauche à droite, la solution exacte, un test du schéma I de Van Leer avec 16000 pas de temps (pour avoir un nombre de Courant plus petit que 1) pour une révolution complète, un test du schéma de Prather avec le même pas de temps et enfin une simulation avec le schéma I mais seulement 160 pas de temps. Les graphiques du haut montrent la distribution de traceur juste après le premier passage par un pôle. Les graphiques du bas montrent le résultat obtenu après une révolution complète autour du globe.

Tests bidimensionnels

Dans cette section, nous présentons des tests bidimensionnels des schémas d'advection tels qu'ils sont codés dans le modèle LMDZT. On teste le schéma de Godunov, le schéma I de Van Leer, le schéma des pentes et le schéma de Prather. Pour le schéma I, on utilise le limiteur fort. Pour les schémas de pentes et de Prather, on assure seulement la positivité en utilisant un limiteur de flux (suivant , , on se contente d'imposer de ne pas sortir plus de traceur d'une maille que ce qu'elle contient). Les schémas de Pentes et de Prather que nous utilisons ont été optimisés et interfacés avec le modèle ancien du LMD par Pascal Simon et Christophe Genthon.

On insiste dans les tests présentés ci-dessous sur le rapport entre précision et coût numérique. On montre en effet que l'arbitrage entre une résolution plus fine, qui rend tous les schémas plus précis, et l'utilisation d'un schéma intrinsèquement plus précis, comme Prather, n'est pas évident et peut dépendre du type de machine utilisé ou d'autres considérations relatives à d'autres composantes du modèle.

Coût numérique des différents schémas

De façon générale, on peut gagner en précision en utilisant une résolution spatiale plus fine, ce qui se fait évidemment au prix d'un coût numérique plus grand. Une question pratique importante en termes d'efficacité des schémas est de savoir ce qu'on gagne en changeant de schéma d'advection à coût numérique inchangé. Puisque notre but est l'advection tridimensionnelle dans un modèle de circulation générale, il faut d'abord se faire une idée du coût relatif des différents schémas dans une telle configuration.

La comparaison la plus facile est celle de l'occupation en mémoire. Les schémas de Godunov, Van Leer I et PPM ne nécessitent de conserver qu'une variable indépendante par maille et par champ de traceur.

Les schémas des pentes et de Prather sont nettement plus coûteux avec respectivement 4 et 10 variables d'état pour décrire un traceur physique.

La comparaison en termes de rapidité est moins évidente et peut dépendre du type de machine utilisé. Des test sur machines scalaires (stations SUN) et vectorielles montrent qu'il y a typiquement un facteur 2 entre le schéma I et les pentes et entre les pentes et Prather.

Les machines vectorielles favorisent de façon générale l'utilisation de schémas plus grossier sur des grilles plus fines. Pour les test présentés ci-dessous, effectués sur un CRAY-90, les pentes ont exactement le même coût numérique pour une grille horizontale de 60 par 43 points que le schéma I de Van Leer pour une grille de 120 par 85 points.

On voit que le fait de doubler la résolution horizontale dans chaque direction rend le schéma I équivalent à celui des pentes à la fois en termes de stockage et en termes de rapidité sur une machine vectorielle. De même, les pentes en résolution double sont comparables à Prather sur une grille deux fois plus grossière.

De façon générale, les machines scalaires sont plus favorables aux schémas plus précis.

Test avec solution exacte connue

Nous présentons ici des tests numériques bidimensionnels effectués avec un champ de vent analytique présentant une rotation différentielle et pour lequel l'advection peut être calculée exactement. Ce test est effectué sur la sphère en utilisant la discrétisation du modèle de circulation, pour permettre de valider directement les codes utilisés dans le modèle de circulation générale.

Pour un champ de vent horizontal non divergent dans le plan longitude-latitude ($\lambda$, $\phi$), les composantes zonale ($u$) et méridienne ($v$) du vent satisfont la relation suivante :

$$\frac{\partial u}{\partial \lambda}+\frac{\partial v\cos\phi}{\partial \phi}=0$$ (2.54)

Le champ de vent que nous utilisons dérive du potentiel suivant :

$$\Psi= a U_0\cos^2\frac{\lambda}{2}\cos^2\phi$$ (2.55)

(où $U_0$ est une vitesse caractéristique et $a$ est le rayon planétaire) de sorte qu'on a

$$u = \frac{1}{a} \frac{\partial \Psi}{\partial \phi} =- 2 U_0\cos\phi\sin\phi \cos^2\frac{\lambda}{2}$$ (2.56)

$$v= \frac{1}{a\cos\phi}\frac{\partial \Psi}{\partial \lambda} = -U_0\cos\phi\cos\frac{\lambda}{2}\sin\frac{\lambda}{2}$$ (2.57)

En suivant une trajectoire (valeur constante de $\Psi$), la vitesse méridienne $v=a \mbox{d}\phi/\mbox{d}t$ peut être récrite en combinant les équations 2.55 et 2.57 :

$$a \mbox{d}\phi/\mbox{d}t= -\sqrt{\Psi U_0/a}\mbox{ sign}\lambda\sqrt{1-\cos^2(\lambda/2)}$$ (2.58)

Si on introduit $\alpha=\sqrt{\Psi U_0/a^3}$ et qu'on remplace $\cos^2(\lambda/2)$ selon l'Eq. 2.55, cette équation devient

$$-\alpha\int \mbox{ sign}{\lambda} \mbox{d}t = \int \frac{\co... ...rcsin\left(\frac{\sin\phi}{\sqrt{1-\Psi/(aU_0)}}\right)\right]$$ (2.59)

Les trajectoires, iso-valeurs de $\Psi$, sont parcourues avec la loi horaire donnée ci-dessus. Le temps mis par une particule pour revenir à sa position de départ est donc $2\pi/\alpha=2\pi a/(U_0\cos\phi\cos\lambda/2)$. Ce temps est plus court au centre du domaine et infini sur les bords.

Résultats numériques

Figure 2.10: Tests d'advection bidimensionnelle avec un champ de vent analytique. La distribution initiale de traceur est une fonction gaussienne de la longitude. Les graphiques du haut montrent l'état initial (à gauche) et la solution exacte aux instants $T/2$ et $T$. Le temps $T$ correspond à une révolution complète au centre du domaine. En dessous, on montre les résultats numériques au temps $T$ pour différents schémas d'advection et 3 résolutions : résolution pleine, 1/2 et 1/3, correspondant respectivement à des grilles longitude-latitude de $120\times 60$, $60 \times 30$ et $40\times 20$ points.

Sur le graphique en haut à gauche de la Fig. 2.10, on montre à la fois la distribution initiale du traceur, une gaussienne ne dépendant que de la longitude, et le champ de vent (dans une unité arbitraire). Toujours en haut, les figures à droite montrent l'évolution exacte du traceur sous l'effet de l'advection. La rotation est plus rapide au centre qu'aux bords du domaine, ce qui produit cette forme en spirale et une filamentation.

Pour les tests numériques on utilise trois résolutions : $120\times 60$ (résolution pleine), $60 \times 30$ (résolution 1/2) et $40\times 20$ (résolution 1/3). Le pas de temps est choisi suffisamment petit pour que le nombre de Courant reste toujours plus petit que 1 (le traitement spécial en longitude n'étant codé que pour le schéma I de Van Leer). L'état final montré sur les figures correspond à une révolution complète au centre du maillage ce qui nécessite un nombre de pas de temps de 4000, 6000 et 12000 pour les différentes résolutions horizontales testées.

La reconstruction très fine des filaments (avec des valeurs maximum dépassant 0.9 pour une valeur initiale de 1) avec le schéma de Prather en pleine résolution est très impressionnante. Cependant, le schéma des pentes fait plutôt mieux à cette résolution que Prather dans une résolution deux fois plus grossière dans chaque direction horizontale. On observe la même chose entre le schéma I de Van Leer et le schéma des pentes.

Donc, si on change la résolution dans seulement deux directions, le schéma I se comporte plutôt mieux que les pentes en terme de rapport qualité/coût.

Si on change de résolution dans les trois directions, la comparaison doit être faite entre la pleine résolution et la résolution intermédiaire avec un avantage significatif pour les schémas précis.

Il est enfin intéressant de noter la différence flagrante de performances entre les schémas de Godunov et le schéma I en dépit d'un coût équivalent en termes de stockage.

Remarques pour conclure

Les schémas en volumes finis proposés par () conduisent facilement à des mises en uvre tridimensionnelles qui satisfont des propriétés essentielles du transport comme

• la conservation de la quantité totale de traceur,
• la non modification d'une distribution constante de traceurs,
• l'invariance par addition d'une constante au champ de traceur,
• la monotonie,
• la non création d'extrema d'origine numérique (d'où la positivité),
• la décroissance de la variation totale (TVD en anglais) - somme des écarts absolus entre deux points consécutifs - au cours du temps qui assure la stabilité du schéma. On peut montrer que cette propriété (valable uniquement en dimension 1) est respectée à la fois par le schéma I de Van Leer et par le schéma PPM (cf. , ).

En pratique, on constate que les schémas plus sophistiqués se comportent mieux, mais au pris d'un coût numérique additionnel du même ordre que celui qu'aurait entraîné l'utilisation d'une grille plus fine. A noter qu'il se peut que les schémas deviennent à partir d'un certain stade moins diffusifs que l'atmosphère elle-même. C'est particulièrement vrai pour la basse troposphère, dans laquelle la turbulence de couche limite induit une très forte diffusion verticale, qui, couplée à des cisaillements importants du vent horizontal, conduit aussi à une forte dispersion horizontale effective. Les tests présentés plus loin illustrent ce point.

Pour finir, il faut noter que nous avons présenté ici les schémas dérivés à l'origine par Van Leer. Le schéma PPM qui a également été testé dans LMDZ (résultats non présentés) semble supérieur à coût numérique égal (en tous cas en termes de stockage) à ces schémas d'origine.

L'ensemble des schémas décrits ci-dessus a été introduit dans le modèle LMDZ. Le schéma I de Van Leer a été retenu en standard pour sa robustesse et sa simplicité mais ce choix ne doit pas être considéré comme définitif et doit être reconsidéré en fonction du problème abordé.

Le transport sous-maille

Pour le terme turbulent ${\mbox{div}}\left(\overline{\rho{\bf v}'c'}\right)$, on distingue en fait trois contributions décrites ci-dessous.

Turbulence de couche limite

Dans la version standard de LMDZ, la turbulence de couche limite est traitée comme une super-viscosité ou viscosité turbulente. Dans ces formulations, comme pour la viscosité moléculaire, le flux d'une quantité transportée est proportionnel (avec un coefficient négatif) au gradient local de la quantité en question. Dans la couche limite planétaire, et si on s'intéresse à l'écoulement à grande échelle, le terme vertical domine de loin ce flux qui s'écrit alors

$$\overline{\rho w'c'}=-K_z\frac{\partial c}{\partial z}$$ (2.60)

Dans la version standard du modèle LMDZ, ce coefficient $K_z$ dépend du cisaillement vertical de vent et d'un nombre de Richardson. La formulation utilisée est donnée par (). Les limitations de cette paramétrisation ainsi que des approches alternatives sont discutées en détail dans le chapitre suivant.

Convection nuageuse

Figure 2.11: Notations pour le transport en flux de masse des traceurs par la convection.

De nombreux développements ont été consacrés ces dernières décennies à la paramétrisation de la convection nuageuse (profonde ou peu profonde), notamment dans le cadre de la modélisation du climat. Les paramétrisations à la mode sont basées sur des approches dites en flux de masse (, ,,). Elles ont en commun d'expliciter des ascendances concentrées, sensées représenter le cur des nuages convectifs. Dans ces ascendances, l'air monte rapidement sous l'effet de sa propre flottabilité, renforcée dans le nuage par le dégagement de chaleur latente.

Certaines de ces paramétrisations considèrent un spectre complet de panaches ascendants. Dans les développements présentés ici, on utilise la paramétrisation de () qui sépare la colonne atmosphérique en trois sous-colonnes : une pour les ascendances, une pour les descentes précipitantes et un troisième compartiment pour l'environnement dans lequel se produit une subsidence compensatoire plus lente.

L'ascendance est caractérisée par un flux de masse $\hat{f}(z)$ qui échange de l'air avec l'environnement. Cet échange est prescrit au travers d'un entraînement $\hat{e}$ et d'un détraînement $\check{d}$. Pour les descentes précipitantes, ont définit de même un flux de masse $\check{f}$, un entraînement $\check{e}$ et un détraînement $\check{d}$.

La colonne convective est supposée stationnaire de sorte que la conservation de la masse d'air entre les différents compartiments s'écrit

$$\frac{\partial \hat{f}}{\partial z}=\hat{e}-\hat{d}$$ (2.61)

et

$$-\frac{\partial \check{f}}{\partial z}=\check{e}-\check{d}$$ (2.62)

Le flux dans l'ascendance et les descentes précipitantes est compensé par un flux, en général plus lent, dans l'environnement, $f_e=-\hat{f}-\check{f}$.

Pour l'inclusion de la composante traceur, on fait les approximations suivantes en suivant la philosophie du schéma d'origine : on suppose que le traceur est dans un régime stationnaire à la fois dans l'ascendance et dans les descentes précipitantes. On suppose de plus que la fraction de la maille couverte par ces deux compartiments est suffisamment faible pour qu'on puisse confondre la concentration dans l'environnement $c_e$ avec la concentration moyenne dans la maille ($c_e=\tilde{c}$ ou $c$).

Sous ces hypothèses, la concentration dans l'ascendance $\hat{c}$ est donnée par

$$\frac{\partial \hat{f}\hat{c}}{\partial z}=\hat{e}{c} -\hat{d}\hat{c}$$ (2.63)

avec une équation similaire pour les descentes

$$-\frac{\partial \check{f}\check{c}}{\partial z}=\check{e}{c} -\check{d}\check{c}$$ (2.64)

(on pourra se reporter à la Fig. 2.11). Enfin, le flux de masse turbulent est donné par

$$\overline{\rho w'c'}=\hat{f}\hat{c}+\check{f}\check{c}-(\hat{f}+\check{f}) {c}$$ (2.65)

Afin d'assurer la stabilité numérique de ce schéma, les différents termes de transport (de la forme $fc$) sont traités avec un schéma amont du premier ordre. ^2.15La diffusion numérique n'est pas un problème ici puisque le processus physique lui-même est très diffusif. Les erreurs numériques associées sont certainement plus faibles que les incertitudes sur l'intensité et la description des échanges d'air dans la colonne convective.

Diffusion latérale

Les termes turbulents associés à du mélange vertical sont souvent nettement plus important que les termes horizontaux. Par exemple dans la basse troposphère, la combinaison d'un cisaillement de vent et d'un mélange vertical turbulent produit une dispersion horizontale des espèces traces extrêmement efficace. Tant que les mailles horizontales sont assez grossières, il est probable de plus que la diffusion numérique soit supérieure à la diffusion latérale réelle de l'atmosphère. Enfin, il faut noter que la théorie physique qui permet d'estimer la diffusivité latérale effective est loin d'être établie.

Cependant, il est probable que, notamment pour une grille zoomée très fine, il commence à être nécessaire d'inclure une paramétrisation de cette diffusion latérale. Ici, cette diffusion est plutôt introduite pour des tests de sensibilité et on retiendra une approche simple en longueur de mélange : comme pour la super-viscosité verticale, le flux horizontal de traceur est relié au gradient local de la quantité. L'effet de cette diffusion latérale sur le transport des traceurs s'écrit alors sous la forme d'un laplacien

$$\frac{\partial c}{\partial t}=\frac{{\delta x}^2}{\nu} \Delta c$$ (2.74)

Architecture

D'un point de vue informatique, l'introduction de la composante traceurs suit l'organisation du modèle LMDZ avec séparation entre dynamique et physique. Les schémas d'advection grande échelle sont interfacés avec le code dynamique. Les flux de masse aux bords des mailles sont en général cumulés dans le temps puisque les schémas d'advection admettent un pas de temps plus long que le code dynamique.

Les parties turbulentes et convectives sont gérées par un moniteur de la ``physique traceurs" interfacé avec la ``physique". Ici, on utilise en général le même pas de temps d'une demi-heure que pour la physique. C'est au niveau de ce moniteur que l'on branche les codes chimiques comme INCA développé par ().

L'organisation de cet outil est résumée sur la Fig. 2.12.

Figure 2.12: Organigramme du modèle LMDZT

Ce modèle est destiné avant tout a des études climatiques couplées dans lesquelles les distributions d'espèces chimiques ou d'aérosols rétroagissent sur les variables météorologiques. Cependant, notamment pour le développement et la validation, il est intéressant de disposer de versions ``débranchées" du modèle et de pouvoir forcer la situation météorologique a suivre au plus près les analyses.

Le modèle a donc été conçu de façon à pouvoir débrancher la météorologie. Les interfaces entre les parties météorologiques et traceurs ont été clairement identifiées. En mode ``branché", on passe à chaque pas de temps les flux de masse pour le transport grande échelle ou les coefficients de mélange turbulent pour la partie physique. On peut également stoker ces variables d'interface sur des fichiers (en pratique, on est obligé de les cumuler sur quelques heures) qui peuvent alors être relus pour effectuer à moindre coût des simulations de transport débranchées.

Sous réserve de cumuler proprement les champs et d'effectuer un découpage propre dans les simulations débranchées, on peut utiliser des pas de temps de quelques heures pour le stockage. Ce point a été documenté en détail par ().

Le mode débranché permet de travailler de façon un peu plus souple quand on s'intéresse à des développements spécifiques à la chimie par exemple, ou de faire des tests de sensibilité en utilisant les même champs météorologiques. Il devient surtout essentiel quand on intègre la dispersion à rebours dans le temps suivant l'approche détaillée dans le Chapitre .

LMDZT, en mode débranché et guidé par les analyses, s'apparente finalement à un modèle de type ``transport-chimie" (Chemistry Transport Models en anglais). Mais, alors que dans les modèles de transport-chimie la météorologie est en général directement issue des réanalyses, on effectue ici une première simulation météorologique guidée. Cette approche offre l'avantage de pouvoir extraire des paramétrisations physiques du modèle tous les paramètres jugés nécessaires pour le transport des espèces traces.

Remarquons enfin que cette version débranchée permet éventuellement de calculer le transport sur une grille plus fine que la grille météorologique. Cette approche est mise en uvre dans LMDZT en redécoupant par exemple chaque maille en 4 sous-mailles horizontalement.

Eléments de validation de la composante traceurs de LMDZ

On montre ici deux exemples de validation du modèle LMDZT en version guidée. Dans la première application, on s'intéresse à l'isotope 222 du Radon, un radio-élément émis par les surface continentales. Ce traceur a été beaucoup utilisé pour valider et inter-comparer les codes de grande échelle. Le second exemple concerne un cas de dispersion d'une source ponctuelle, dans le cadre de la campagne européenne ETEX. On utilise une grille zoomée pour ce cas.

Simulations Radon

Figure 2.13: Mesure (points) et simulation (trait) du $^{222}$Rn (mBq/m$^3$) à l'île d'Amsterdam pour l'année 2000. La simulation est réalisée avec le modèle LMDZT guidé par les analyses ECMWF en utilisant une grille régulière avec une résolution de $2^o\times 2^o$.

$\includegraphics[angle=-90,width=12cm]{lmdzt/FIGURES/AMS1.eps}$

Janvier, Février

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/AMS2.eps}$

Mars, Avril

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/AMS3.eps}$

Mai, Juin

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/AMS4.eps}$

Juillet, Août

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/AMS5.eps}$

Septembre, Octobre

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/AMS6.eps}$

Novembre, Décembre

Figure 2.14: Mesure (points) et simulation (trait) du $^{222}$Rn (mBq/m$^3$) à Mace Head pour l'année 2000. La simulation est réalisée avec le modèle LMDZT guidé par les analyses ECMWF en utilisant une grille régulière avec une résolution de $2^o\times 2^o$.

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/MHD1.epsi}$

Janvier, Février

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/MHD2.epsi}$

Mars, Avril

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/MHD3.epsi}$

Mai, Juin

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/MHD4.epsi}$

Juillet, Août

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/MHD5.epsi}$

Septembre, Octobre

$\includegraphics[height=12cm,angle=-90]{lmdzt/FIGURES/MHD6.epsi}$

Novembre, Décembre

L'isotope 222 du Radon ($^{222}$Rn) est un gaz trace radioactif particulièrement adapté à la validation des modèles de transport aux échelles temporelles de quelques heures à quelques semaines (, ). Il est formé à partir de la décomposition radioactive de l'isotope 238 de l'uranium, présent dans les sols. Les sources océaniques sont de ce fait de deux à trois ordres de grandeur plus faibles que les sources continentales. On peut donc voir le $^{222}$Rn comme un traceur de l'air continental. Le seul puits significatif du $^{222}$Rn est sa décomposition en Polonium-218 avec une période de décroissance de 5,5 jours.

Le $^{222}$Rn a été largement utilisé pour des validations et inter-comparaisons de modèles de transport en supposant une source uniforme sur les continents. Ici, on utilise une version un peu plus sophistiquée, développée par Christophe Genthon (LGGE) et Alexandre Armengaud (, ) dans laquelle le $^{222}$Rn diffuse vers la surface à travers une couche de sol poreuse. C'est alors le flux en bas de cette couche qui est prescrit.

Nous présentons une simulation guidée du $^{222}$Rn pour l'année 2000 avec une résolution horizontale de 2 degrés par 2 degrés. On compare les valeurs observées et simulées à deux stations d'observation : l'île d'Amsterdam dans l'océan indien austral (Fig. 2.13) et Mace Head à l'extrémité occidentale de l'Irlande (Fig. 2.14). Les deux jeux de données nous ont été aimablement communiqués par Michel Ramonet et Philippe Ciais (LSCE).

Ces deux stations sont sous forte influence océanique et font donc apparaître des contrastes marqués entre un fond relativement faible et des bouffées de Radon associées à des arrivées d'air continental, soit en situation de vent d'Est sur l'Europe pour Mace Head soit quand la circulation amène sur l'île d'Amsterdam des masses d'air en provenance d'Afrique du Sud.

On voit que le modèle reproduit assez bien à la fois les valeurs faibles dans les périodes sous influence océanique et une bonne partie des pics. Des simulations pour les années 1991 et 1992 (non montrées) ont été comparées aux résultats publiés par (). Cette comparaison montre que le modèle se comporte raisonnablement et qu'une partie des pics non simulés provient d'erreurs sur les champs de vents analysés. Noter que la forte dépendance de la source au contenu en eau des sols est négligée ici. Elle pourrait expliquer les surestimations systématiques observées à Mace Head du 15 mai au 10 juin ou fin septembre.

Evaluation sur la campagne ETEX

Nous présentons dans cette section quelques résultats d'un travail effectué par Abderrahmane Idelkadi durant sa thèse (, ) concernant l'évaluation du modèle LMDZT en se basant notamment sur la comparaison à une série de modèles de dispersion ayant participé à la campagne ETEX (, ). La comparaison est faite par rapport à des résultats publiés il y a quelques années de cela et il est donc possible que certains de ces modèles possèdent aujourd'hui des versions plus récentes et plus performantes.

Pour cette évaluation, nous utilisons les résultats de la campagne ETEX (European Tracer EXperiment) organisée en 1994 par l'organisation mondiale de la météorologie, la Commission Européenne et l'Agence Internationale de l'Energie Atomique. Une quantité totale de 340 kg d'un gaz insoluble, le PMCH (Perfluoro-Methyl-Cyclo-Hexane), a été émise le 23 octobre 1994, à partir de 16h00 UTC (T0) et durant 12 heures, à travers une cheminée de 8 m située à Monterfil près de Rennes. Les conditions météorologiques ont été choisies de manière à ce que le traceur soit vu par le plus grand nombre possible des 168 stations d'observation réparties sur l'Europe. L'expérience a duré trois jours et les mesures de concentration du polluant ont été effectuées toute les trois heures à chaque station.

Figure 2.15: Grille retenue pour les simulations ETEX-1 et localisation des 11 stations retenues pour les diagnostics d'intercomparaison.

$\includegraphics[height=17cm,angle=-90]{lmdzt/FIGURES/grille.eps}$

La Fig. 2.15 montre la grille LMDZ choisie pour la simulation avec un zoom sur l'Europe. On montre aussi en incrustation l'emplacement de 11 stations privilégiées par les organisateurs pour certains diagnostics lors des études d'intercomparaison.

Figure 2.16: Panaches de PMCH (ng m$^{-3}$) observés et simulés pour la campagne ETEX. Le panache est reconstitué à partir des mesures aux stations en appliquant un filtre de () utilisable depuis le logiciel graphique du domaine publique GrADS utilisé pour ces figures. Les concentrations sont en ng par kg d'air. Les heures sont comptées par rapport au temps T0 de relâchage du PMCH. Pour la simulation, plutôt que de tracer le panache directement, on commence par extraire pour chacune des 168 stations la séquence des concentrations simulées (en prenant la valeur au point de grille le plus proche) puis on reconstitue le panache avec la même méthode que pour les observations.

Observations	Simulation

La Fig. 2.16 montre les panaches observés et simulés à différents instants. On voit que le panache est globalement bien reproduit, à la fois en répartition spatiale et en intensité. Pour quantifier cet accord, on utilise l'un des critères statistiques retenus pour l'intercomparaison des modèles dans le cadre d'ETEX. Il s'agit du FMT pour Figure of Merit in Time. Le FMT se définit comme le rapport (traduit en pourcentage)

$$FMT_i=100\times \frac{\int_0^T MIN[c^{\mbox{\footnotesize ob... ...footnotesize obs}}_i(t),c^{\mbox{\footnotesize sim}}_i(t)] dt}$$ (2.75)

où $MIN(x,y)$ (resp. $MAX(x,y)$) est le minimum (resp. maximum) de $x$ et $y$; $c^{\mbox{\footnotesize obs}}_i(t)$ et $c^{\mbox{\footnotesize sim}}_i(t)$ sont les concentrations observées et simulées à la station $i$ à l'instant $t$. Le FMT vaut au maximum 100$\%$ quand observation et simulation coïncident. Il est sensible à la fois aux erreurs d'intensité et aux déphasages temporels.

Figure 2.17: Schéma illustrant le critère FMT (à gauche) et FMTs calculés pour les 11 stations de la Fig. 2.15 et pour un grand nombre de modèles utilisés dans le cadre d'ETEX et pour LMDZT (à droite).

La Fig. 2.17 montre les FMTs obtenus pour les 11 stations montrées sur la Fig. 2.15. Ces FMTs sont comparés à ceux obtenus avec une trentaine d'autres modèles. Les conclusions des inter-comparaisons ETEX avaient identifié une poignée de modèles qui se plaçaient au dessus du lot. Les performances relatives des différents modèles semblaient d'ailleurs peu liées à leurs caractéristiques, que ce soit au type d'approche (lagrangienne ou eulérienne par exemple) ou aux résolutions spatiales retenues. Il est fort probable d'ailleurs que les meilleurs modèles étaient en grande partie limités par le caractère imparfait des analyses météorologiques. Le modèle LMDZT n'a rien d'exceptionnel mais fait partie de ces bons modèles.

Figure 2.18: Sensibilité des simulations ETEX à la résolution horizontale et à la dissipation. La simulation de contrôle a une maille de 120 km de côté. On teste des grilles deux fois plus fines dans les deux directions. La simulation ``Haute résolution" correspond à un cas où on recalcule la météorologie et le transport sur une grille plus fine. Dans la simulation avec ``Redécoupage horizontal", on utilise les archives de la simulation de contrôle mais on redécoupe la maille horizontalement pour le transport suivant le schéma de droite. Cette solution permet de diminuer la diffusivité sans avoir à recalculer la météorologie. Ses résultats sont très proches de ceux de la simulation ``Haute résolution". On teste également l'impact d'introduire une dissipation horizontale avec pour coefficient de dissipation $K_h$, 10$^5$ ou 10$^6$ m$^2$ s$^{-1}$.

Les cartes suggèrent que le modèle est un peu trop diffusif. Les panaches à 24 et 48 heures sont en particulier plus étendus que les panaches observés. Des tests de sensibilité à la résolution (Fig. 2.18) indiquent que, quand on raffine la grille d'un facteur 2 dans les deux directions horizontales, on diminue la diffusion horizontale, mais en réduisant plutôt l'accord avec les mesures. A partir d'une maille de quelques dizaines de kilomètres, le modèle semble donc sous-estimer la diffusion horizontale. Ces résultats sont présentés plus en détail par ().

Discussion

On voit que le code LMDZT peut être utilisé comme un code de dispersion atmosphérique. Mais LMDZT se distingue par de nombreux aspects des codes développés spécifiquement pour les calculs de dispersion.

D'abord, un certain nombre de ces codes sont bâtis sur des trajectoires particulaires ou Lagrangiennes. Cette approche semble relativement intuitive quand on a une source ponctuelle comme celle d'ETEX. Cependant, elle pose le problème que plus la dispersion devient importante, et plus le nombre de particules à injecter est grand si on veut pouvoir prédire les concentrations faibles observées loin des sources. Le traitement du transport turbulent est aussi une difficulté des codes particulaires. Certains font simplement l'impasse. D'autres utilisent des approches de type marche aléatoire. Par exemple, () commencent par estimer la hauteur de la couche limite puis déplacent aléatoirement les particules au sein de cette couche limite. On voit qu'il faut à nouveau ajouter un grand nombre de particule pour obtenir un traitement statistique correct de ce phénomène. Les codes eulériens bénéficient pour leur part des nombreux efforts de recherche développés dans les modèles météorologiques.

Le modèle LMDZT se distingue aussi des modèles de transport-chimie comme on l'a vu, en ce sens qu'il calcule sa propre météorologie. Pour la recherche climatique, les motivations pour cette approche sont claires. Le but ultime est de prendre en compte, de façon interactive, les rétroactions des espèces transportées sur la météorologie. Dans ce cadre, les versions guidées et débranchées du modèle ne sont qu'un mode particulier d'utilisation permettant la validation sur des campagnes d'observation. L'avantage potentiel pour simuler la dispersion dans un contexte de surveillance de l'environnement est plus subtile. Dans les modèles de transport classiques, des phénomènes physiques comme le mélange turbulent dans la couche limite ou le lessivage par les pluies doivent être calculés d'une façon ou d'une autre. Il faut alors essayer de diagnostiquer, à partir de champs météorologiques incomplets (disponibles toutes les 6 heures seulement par exemple), des coefficients de mélange turbulents ou des taux de précipitation dans l'atmosphère. Dans l'approche retenue ici, on effectue une simulation météorologique complète, dans laquelle ces différents aspects sont représentés, tout en guidant la simulation pour qu'elle colle au plus près aux champs de vents de grande échelle issus des analyses. A noter que pour répondre à la même préoccupation, les centres qui produisent les analyses et réanalyses météorologiques se sont mis petit à petit à archiver pour les modèles de transport-chimie des variables internes des paramétrisations comme les coefficients de diffusion turbulente ou les flux de masse convectifs.

Enfin, le modèle LMDZT se distingue des modèles régionaux par l'utilisation d'une grille globale à maille variable. Pour des applications où on utilise des observations réparties uniformément sur le globe, comme dans le cas du TICE exposé dans le Chapitre , le modèle global s'impose naturellement. Dans le cas d'évènements relativement localisés, les modèles à domaine limité présentent un meilleur rapport précision/coût. Cependant, même pour une source ponctuelle, le modèle global est intéressant en ce sens qu'on n'a pas à se poser à l'avance le problème du choix du domaine. Le zoom est effectué sur le point source. Tant que le panache est proche de la source, et donc relativement concentré, le calcul est très précis. Plus on s'éloigne de cette source et moins le calcul est précis. Mais cette perte de précision a aussi moins d'importance puisque le panache est de toutes façons beaucoup plus diffus.

Cette remarque s'applique de la même façon pour les rétro-simulations (cf. Chapitre ). Si on interprète une mesure à une station, par exemple la mesure d'une concentration élevée d'un polluant suggérant un accident industriel, une rétro-simulation à partir de la station avec un zoom aux abords de cette station permettra de bien décrire l'origine de l'air, à la fois finement près de la station mais également à l'autre bout du globe. Il ne sera donc pas nécessaire de faire des hypothèses a priori sur l'origine de la pollution.

Applications et perspectives

Applications "traceurs" avec le modèle LMDZ

Les développements présentés plus haut et mis en musique dans le modèle LMDZ sont impliqués dans un grand nombre d'études dont les plus importantes sont brièvement décrites ci-dessous.

Etude des couplages chimie-climat et aérosols-climat sur Terre

Au LOA, Olivier Boucher a été à l'origine d'une partie des développements concernant l'introduction de la composante traceurs dans LMDZ (notamment pour ce qui est du transport convectif). Il a depuis introduit dans LMDZ une représentation en ligne du cycle du soufre. Cette version contient 7 traceurs advectés (DMS, 2, H2S, DMSO, MSA, sulfates et 2O2). Les principaux oxydants (OH, 2, 3 et Ø3) sont pour leur part prescrits. Cette version du modèle a notamment été utilisée pour évaluer l'évolution entre l'époque pré-industrielle et l'époque actuelle du forçage radiatif des aérosol soufrés (, ). Une part importante de la problématique concerne l'évaluation de l'effet indirect de ces aérosols (, ) au travers des changements des propriétés microphysiques des nuages : diminution de la taille et augmentation du nombre des gouttes d'eau (premier effet indirect) à cause du nombre accru de noyaux de condensation en cas d'augmentation des aérosols, entraînant également une précipitation moindre et une durée de vie accrue des nuages (second effet indirect). Cette version a été récemment étendue aux autres composantes de l'aérosol et validée par rapport aux résultats de la campagne INDOEX (, ). Cette version avec aérosols a également été utilisée pour prédire la possible rétroaction d'un réchauffement climatique sur les aérosols naturels : notamment sur les sulfates à cause de la modification des DMS marins (, ,) mais aussi pour les sels marins et les aérosols désertiques, dont l'émission est sensible aux vents en surface.

Cette version du modèle avec cycle du soufre a récemment été adaptée aux hautes latitudes par l'équipe de Christophe Genthon au LGGE (, ) en vue notamment de l'interprétation des reconstitutions historiques effectuées à partir des calottes de glace (voir aussi , ).

Plus récemment, Didier Hauglustaine a développé pour LMDZ un module de chimie et aérosols interactifs, INCA. Les développements et études menés avec LMDZ-INCA se sont principalement portés sur une version 4-NOx-Ø3 troposphérique, visant principalement a étudier les gaz à effet de serre autres que le 2. Un très important travail de validation a été réalisé sur cette version dans laquelle une quarantaine d'espèces sont advectées (, ).

Cette version a été utilisée pour évaluer l'impact radiatif d'un changement des émissions de CO et des NOx via la modification de la distribution de l'ozone et des radicaux OH. Le modèle a également été utilisé pour étudier l'impact des émissions liées aux transport aérien sur la composition chimique de l'atmosphère et sur le climat.

Le modèle a également été couplé au module REPROBUS développé par Franck Lefèvre pour la chimie stratosphérique.

L'étude des grands cycles climatiques martiens

La version martienne du modèle LMDZ (, ,) s'est enrichie au fil des années. Les poussières, dont l'impact radiatif est très important même en dehors des grandes tempêtes de poussières globales, ont tout d'abord été incluse. Puis le cycle de l'eau (, ), très important pour essayer de contraindre les réservoirs d'eau actuels et attaquer les questions relatives aux climats passés de Mars. Une composante chimique a également été inclue dans le modèle, tout d'abord pour étendre le modèle jusqu'à la thermosphère dans le cadre du développement d'une base de données atmosphérique martienne pour l'ASE et le CNES (, ), puis, plus récemment avec le développement d'un module de chimie pour la basse atmosphère (0 à 120 km) (, ).

Autres applications

La version traceurs de LMDZ est également utilisée dans deux grands types d'applications qui font l'objet de deux chapitres spécifiques du présent document et ne sont donc que mentionnés ici.

Le modèle est d'abord utilisé en mode rétro-transport (Chapitre ) dans deux cadres principalement : l'inversion des puits et sources de 2 et la détection des essais nucléaires.

Le modèle a également été beaucoup utilisé pour l'étude des couplages entre dynamique atmosphérique, chimie et microphysique des brumes dans l'atmosphère de Titan (Chapitre ). A noter qu'une version vénusienne du modèle est actuellement en cours de développement.

Inclusion de la composante traceurs dans le modèle de convection d'Emanuel

Figure 2.19: Concentration moyenne de $^{222}$Rn (10$^{-21}$ mol/mol) en Janvier obtenue avec les paramétrisations de Tiedtke et Emanuel pour la convection profonde ainsi que la différence relative - (Emanuel - Tiedtke)/Tiedtke - en pourcentage.

Un travail a été effectué récemment par Marie-Angèle Filiberti et Jean-Yves Grandpeix pour introduire la composante traceurs dans le schéma convectif d'Emanuel. D'un point de vue du transport, la différence principale entre le schéma d'Emanuel et celui de Tiedtke est la possibilité pour une parcelle d'air sortant de l'ascendance convective à un certain niveau d'être détraînée dans l'environnement à n'importe quel niveau du modèle. La fermeture (permettant de calculer le flux à la base de la colonne convective en fonction d'autres paramètres du modèle de circulation) est également très différente et la convection pénètre généralement plus haut avec le schéma d'Emanuel. On n'entre pas ici dans les détails de la formulation mais on illustre simplement sur la Fig. 2.19 l'importance de la représentation du transport convectif pour le transport des traceurs.

La figure montre, pour un mois de janvier et en moyenne zonale, la structure méridienne de la concentration de $^{222}$Rn obtenue avec le schéma de Tiedtke (en haut), avec celui d'Emanuel (au milieu) et la différence relative. Dans la région de convergence intertropicale, localisée principalement sur l'Afrique, l'Amérique du sud et l'Indonésie à cette saison, entre 30S et l'équateur, le Radon est détraîné beaucoup plus haut en altitude avec le schéma d'Emanuel et beaucoup moins dans la troposphère moyenne. Aux autres latitudes, l'effet du schéma d'advection est sans doute moins directe. Les concentrations de Radon peuvent notamment être affectées par des variations de la circulation à grande échelle. On n'entre pas ici dans ces considérations.

Cette grande sensibilité au transport atmosphérique pourrait laisser penser que les traceurs peuvent fournir une contrainte forte sur la représentation du transport, et notamment la convection. Le $^{222}$Rn a d'ailleurs été beaucoup utilisé pour des inter-comparaisons de modèles globaux (cf. par exemple, , ,). Malheureusement, même avec les différences très importantes obtenues ici, les observations de la concentration restent souvent insuffisantes pour valider réellement les modèles et choisir entre deux paramétrisations. Des campagnes dédiées, comme la campagne AMMA d'Analyse Multi-échelles de la Mousson Africaine, avec déploiement d'instruments aéroportés autour des systèmes convectifs combinant détection passive et active et mesures météorologiques et chimiques devraient permettre d'avancer sur ce point.

Transport des variables actives

On peut bien sûr envisager d'utiliser, pour les variables météorologiques du modèle de circulation générale, les schémas en volumes finis introduits dans le modèle pour transporter les espèces traces. C'est en fait déjà le cas pour la vapeur d'eau dont l'advection est maintenant calculée dans la version standard du modèle avec le schéma I de Van Leer. L'advection de l'eau et de la température potentielle étaient déjà de fait écrites dans la version originale du modèle comme des schémas en volumes finis, avec un flux calculé simplement comme le produit du flux de masse par une interpolation linéaire de la quantité transportée à l'interface :

$$\frac{\partial \left(m\theta \right)}{\partial t} +{\cal F}\... ... +\delta_z \left({\overline{\theta }}^{ Z } W\right)=S_\theta$$ (2.76)

où $S_\theta$ est le terme source dû au chauffage diabatique et où le filtre longitudinal dans les hautes latitudes ${\cal F}\left[ \right]$^2.16s'applique aussi à l'équation de continuité

$$\frac{\partial m}{\partial t} + {\cal F}\left[ \delta_x U+ \delta_y V \right] +\delta_z W=0$$ (2.77)

Pour la température potentielle, on pourrait donc également facilement utiliser un autre schéma en volumes finis.

Ce qui est moins direct, c'est de remarquer que la propriété de conservation de l'enstrophie, importante pour la stabilité du modèle et pour la bonne représentation des transferts entre échelles pour les écoulements stratifiés (, ,), découle en fait d'une formulation volumes finis cachée de l'advection de la vorticité. Cette remarque permet d'envisager une réécriture complète du code dynamique, proche de l'esprit actuel, mais dans laquelle on remplacerait les schémas d'advection par des schémas en volumes finis présentant de meilleurs propriétés physiques comme le schéma I de Van Leer ou PPM.

En introduisant le facteur de Coriolis multiplié par l'aire de la maille, $f=2\Omega \sin{\phi} c_uc_v$ où $\Omega$ est la vitesse de rotation de la planète, ainsi que la vorticité potentielle absolue

$$Z=\frac{{\cal F}\left[ \delta_x \tilde{v}- \delta_y \tilde{u} \right]+ f}{m}$$ (2.78)

et l'énergie cinétique

$$K=\frac{1}{2} \left( {\overline{\tilde{u}\tilde{\tilde{u}}}}^{ Y } + {\overline{\tilde{v}\tilde{\tilde{v}}}}^{ X } \right)$$ (2.79)

les équations du mouvement discrétisées prennent la forme suivante:

$$\frac{\partial \tilde{u}}{\partial t} - {\overline{Z}}^{ Y }... ...^{ X } \delta_x {\cal F}\left[ \Pi \right] +WDU =S_{\tilde{u}}$$ (2.80)

$$\frac{\partial \tilde{v}}{\partial t} + {\overline{Z}}^{ X }... ...^{ Y } \delta_y {\cal F}\left[ \Pi \right] +WDV =S_{\tilde{v}}$$ (2.81)

où $\Pi =C_pp^\kappa$ est la fonction d'Exner au milieu de la maille.

Advection verticale

Les termes d'advection verticale de $u$ et $v$, respectivement $WDU$ et $WDV$, ont été modifiés dans le passé pour garantir partiellement la conservation du moment cinétique par le modèle (, ). Les formulations originales non conservatives et conservatives de ces termes sont données dans la Table 2.1.

En fait, Robert Sadourny (communication personnelle) avait remarqué fort justement que la conservation du moment cinétique était également garantie si la double moyenne en $Y$ sur le terme en $u$ était appliquée plutôt à $W$, ce qui conduisait à lisser le champ de vitesse verticale pour l'advection verticale. En faisant cette transformation, on peut revenir à la forme originelle avec une triple moyenne sur les champs de vitesse verticale (formulation 2 dans la table) ce qui est plus sympathique. En pratique, on peut commencer par calculer ${\overline{W}}^{ XY }$ avant la moyenne. Cette formulation 2 peut facilement être testée dans le modèle actuel.

Tableau 2.1: Différentes formulations possibles pour le terme d'advection verticale dans l'équation du mouvement (voir le texte pour les détails).

	$WDU$	$WDV$
non conservative	$\frac{{\overline{{\overline{W}}^{ X } \delta_z \tilde{u}}}^{ Z }}{{\overline{m}}^{ X }}$	$\frac{{\overline{{\overline{W}}^{ Y } \delta_z \tilde{v}}}^{ Z }}{{\overline{m}}^{ Y }}$
conservative 1	$- \frac{{\overline{\tilde{u}_{abs}}}^{ Y,Y } \delta_z {\overline{W}}^{ X } } {... ...W}}^{ X } {\overline{\tilde{u}_{abs}}}^{ Z } \right) } {{\overline{m}}^{ X } }$	$- \frac{{\overline{\tilde{v}}}^{ X,X } \delta_z {\overline{W}}^{ Y } } {{\over... ...erline{W}}^{ Y } {\overline{\tilde{v}}}^{ Z } \right) } {{\overline{m}}^{ X }}$
conservative 1b	$- \frac{\tilde{u}_{abs}\delta_z {\overline{W}}^{ XYY } } {{\overline{m}}^{ X }... ...{W}}^{ X } {\overline{\tilde{u}_{abs}}}^{ Z } \right) } {{\overline{m}}^{ X }}$	$- \frac{\tilde{v}\delta_z {\overline{W}}^{ XXY } } {{\overline{m}}^{ Y }} + \f... ...erline{W}}^{ Y } {\overline{\tilde{v}}}^{ Z } \right) } {{\overline{m}}^{ X }}$
conservative 2	$\frac{{\overline{{\overline{W}}^{ XYY } \delta_z \tilde{u}}}^{ Z }}{{\overline{m}}^{ X }}$	$\frac{{\overline{{\overline{W}}^{ XXY } \delta_z \tilde{v}}}^{ Z }}{{\overline{m}}^{ Y }}$

Vers les volumes finis

On réalise en fait facilement que l'astuce principale du schéma original de Sadourny, conservatif en enstrophie, a consisté à récrire la partie vorticité de l'équation du mouvement sous la forme d'un flux de vorticité. En effet, avec les notations ci-dessus, la version numérique de l'équation de bilan de vorticité s'écrit :

$$\frac{\partial \left({\overline{m}}^{ X,Y } Z\right)}{\partial t} = \frac{\partial {\cal F}\left[ \delta_x \tilde{v}- \delta_y \tilde{u} \right]}{\partial t}$$ (2.82)

$$= {\cal F}\left[ \delta_x \frac{\partial \tilde{v}}{\partial t}-\delta_y \frac{\partial \tilde{v}}{\partial t} \right]$$ (2.83)

$$= -{\cal F}\left[ \delta_x \left({\overline{Z}}^{ X } {\overline{U}... ...ight)+\delta_y \left({\overline{Z}}^{ Y } {\overline{V}}^{ X,Y }\right) \right]$$ (2.84)

$$+{\cal F}\left[ \delta_y \left({\overline{\theta }}^{ X } \delta_... ...line{\theta }}^{ Y } \delta_y {\cal F}\left[ \Pi \right]-WDV-S_v\right) \right]$$

Figure 2.20: Grilles volumes-finis équivalentes correspondant à l'advection des scalaires (à gauche) et de la vorticité (à droite) dans la formulation actuelle du modèle. Les tiretés correspondent aux délimitations des volumes de contrôle.

La disparition du terme $\Phi+K$ provient de l'égalité $\delta_x \delta_y =\delta_y \delta_x$.

Pour des écoulements barotropes ($W=0$ et $\delta_x \Pi =\delta_y \Pi =0$), on obtient une équation de bilan

$$\frac{\partial \left({\overline{m}}^{ X,Y } Z\right)}{\partial t}... ...ht)+\delta_y \left({\overline{Z}}^{ Y } {\overline{V}}^{ X,Y }\right) \right]=0$$ (2.85)

analogue à celle utilisée pour la température potentielle ou les traceurs à ceci près qu'on utilise une double moyenne des flux de masse pour travailler sur la grille en vorticité.

La figure 2.20 illustre le positionnement des variables dans la formulation actuelle si on les interprète en termes de volumes finis.

Il est donc tentant, de la même façon que pour les traceurs, ou la température potentielle, de remplacer les moyennes arithmétiques par des estimations amont à la Van Leer.

On n'obtiendra plus alors, comme dans l'ancien modèle, une conservation exacte de l'enstrophie. En revanche, en appliquant un schéma positif, monotone, etc ... on empêchera la création d'extrema numérique du champ de vorticité. On obtiendra également automatiquement de la dissipation à l'échelle de la maille en présence d'un extremum de vorticité important.

Changement de grille

On peut essayer d'appliquer directement une formulation volumes finis aux équations ci-dessus. Il y a un petit problème pratique : les schémas volumes finis utilisés pour la vorticité et la température potentielle travailleront sur des grilles décalées.

Cependant, on se rend compte qu'en décalant les variables scalaires sur la grille des points de vorticité, on obtient une écriture plus systématique.

Dans la nouvelle configuration, la position des points $\tilde{u}$ et $\tilde{v}$, des points de vorticité et des flux de masse utilisés pour l'advection est inchangée. Mais, les scalaires doivent être advectés, comme la vorticité, avec les flux de masse moyens en $X$ et en $Y$.

Il est en fait plus clair de redéfinir les variables $U$ et $V$ comme des doubles moyennes des flux de masse obtenus à partir de $\tilde{\tilde{u}}$ et $\tilde{\tilde{v}}$ :

$${\cal U}={\overline{\tilde{\tilde{u}}{\overline{m}}^{ Y }}}^... ...erline{\tilde{\tilde{u}}}}^{ XY }{\overline{m}}^{ Y }\mbox{ )}$$ (2.86)

et

$${\cal V}={\overline{\tilde{\tilde{v}}{\overline{m}}^{ X }}}^{ XY }$$ (2.87)

(les masses elle mêmes sont décalées comme les scalaires).

L'équation de continuité s'écrit alors formellement comme précédemment

$$\frac{\partial m}{\partial t} + {\cal F}\left[ \delta_x {\cal U}+ \delta_y {\cal V} \right] +\delta_z {\cal W}=0$$ (2.88)

à ceci près qu'on travaille maintenant sur la grille de vorticité.

La pression de surface est encore la somme verticale de la masse

$$\frac{{\cal A}}{g}{p_s} = \sum_{l=1}^N m$$ (2.89)

La vorticité potentielle absolue devient

$$Z=\frac{{\cal F}\left[ \delta_x \tilde{v}- \delta_y \tilde{u} \right]+ f}{m}$$ (2.90)

Pour l'énergie cinétique, on a le choix entre la calculer comme avant sur l'ancienne grille scalaire,

$$K=\frac{1}{2} \left( {\overline{\tilde{u}\tilde{\tilde{u}}}}^{ Y } + {\overline{\tilde{v}\tilde{\tilde{v}}}}^{ X } \right)$$ (2.91)

ou la calculer sur la grille de vorticité

$${\cal K}=\frac{1}{2} \left( {\overline{\tilde{u}\tilde{\tild... ...^{ Y } + {\overline{\tilde{v}\tilde{\tilde{v}}}}^{ X } \right)$$ (2.92)

puis en prendre la moyenne en $X$ et $Y$ pour la ramener aux points scalaires anciens

Les équations du mouvement prennent finalement la forme suivante:

$$\frac{\partial \tilde{u}}{\partial t} - {\overline{Z}}^{ Y }... ...elta_z \tilde{u}}}^{ Z }}{{\overline{m}}^{ Y }} =S_{\tilde{u}}$$ (2.93)

$$\frac{\partial \tilde{v}}{\partial t} + {\overline{Z}}^{ X }... ...elta_z \tilde{v}}}^{ Z }}{{\overline{m}}^{ X }} =S_{\tilde{v}}$$ (2.94)

avec

$$K^*=K\mbox{ ou }{\overline{\cal K}}^{ XY }$$ (2.95)

Pour le terme d'advection verticale, on utilise ici l'équivalent de la version "conservation 2" puisque les flux de masses horizontaux (et donc le flux de masse vertical) sont déjà moyennés dans les deux directions horizontales. De ce fait, ${\overline{W}}^{ Y }$ dans la nouvelle formulation correspond à ${\overline{W}}^{ XYY }$ dans l'ancienne. Il est intéressant de remarquer que le choix original de la forme non conservative était le choix le plus ``naturel" et que des modifications avaient dû être introduites après coup pour garantir la conservation du moment cinétique zonal. Ici, le choix ``naturel" correspond au choix conservatif.

L'équation thermodynamique

s'écrit formellement comme avant

$$\frac{\partial \left(m\theta \right)}{\partial t} +{\cal F}\... ... +\delta_z \left({\overline{\theta }}^{ Z } W\right)=S_\theta$$ (2.96)

Une maquette de cette nouvelle formulation a été développée par Phu LeVan sur une grille longitude-latitude. Elle est actuellement en cours de test. Si cette approche s'avère efficace et robuste, elle pourrait donner une excellente base pour développer une dynamique icosaédrique, ^2.17vieux rêve rentré de Robert Sadourny, qui hante encore les étagères de certains collègues au laboratoire.

À propos de ce document...

Représentation du transport direct et inverse dans les modèles globaux de climat et étude des couplages entre composition et dynamique atmosphérique sur Titan.

This document was generated using the LaTeX2HTML translator Version 2002-2-1 (1.71)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split 0 ttt

The translation was initiated by HOURDIN Frédéric on 2006-06-14

Notes

... finis^2.1

C'est à la suite d'un séminaire donné par Sandrine Edouard en thèse avec Bernard Legras sur l'advection d'ozone à la frontière du vortex polaire (, ,) que j'ai eu le bonheur de découvrir l'article de () sur les volumes finis et que je me suis retrouvé embarqué dans l'introduction des traceurs dans le modèle LMDZ et tout ce qui s'en suit...

... récemment.^2.2

Cette section, relativement technique, ne permettra pas à un non initié de se familiariser avec le monde des modèles de circulation générale.

... flexible^2.3

La flexibilité est ici une notion positive !

....^2.4

Pour résoudre cette équation, on commence par calculer la convergence de masse $\Omega_k$ cumulée depuis le sommet de l'atmosphère jusqu'au niveau $k$ :

$$\Omega_k= \sum_{l=k}^N \left(\delta_x U+ \delta_y V\right)$$ (2.7)

La convergence au sol $\Omega_1$ donne accès à l'évolution de la pression de surface

$$\frac{{\cal A}}{g}\frac{\partial p_s}{\partial t}=\Omega_1$$ (2.8)

Finalement, en remarquant que

$$\frac{\partial m}{\partial t}=-\frac{{\cal A}}{g}\delta_z B \frac{\partial p_s}{\partial t}=-\delta_z B \Omega_1$$ (2.9)

on obtient le flux de masse vertical $W$ par intégration depuis le haut de l'atmosphère

$$- \delta_z B \Omega_1 -\delta_z \Omega +\delta_z W=0$$ (2.10)

....^2.5

$\kappa=R/C_p$ où $R$ et $C_p$ sont les constantes thermodynamiques de l'air sec.

... savoir^2.6

A partir de l'équation hydrostatique 2.11, on peut écrire

$$\Phi_l=\Phi_s-\sum_{k=1}^l \left[{\overline{\theta }}^{ z } \delta_z \Pi \right]_k$$ (2.14)

avec la convention (utilisée dans le modèle) $\left[{\overline{\theta }}^{ z }\right]_1=\theta _1$. L'énergie potentielle de gravité totale de la colonne s'écrit donc

$$\sum_{l=1}^N\Phi_lm-\Phi_s = -\sum_{l=1}^Nm_l \sum_{1\le k \le l} \left[{\overline{\theta }}^{ z } \delta_z \Pi \right]_k$$ (2.15)

$$= -\sum_{k=1}^N\sum_{k\le l \le N} \left[{\overline{\theta }}^{ z } \delta_z \Pi \right]_k m_l$$ (2.16)

$$= -\sum_{k=1}^N\left[{\overline{\theta }}^{ z }\delta_z \Pi \right]_k \times \frac{{\cal A}}{g} p_k$$ (2.17)

$$= -\frac{{\cal A}}{g}\sum_{k=1}^N\theta _k\left[{\overline{p \delta_z \Pi }}^{ z }\right]_k$$ (2.18)

La dernière transformation, qui consiste à faire glisser la moyenne d'une demi maille sur la verticale, est valable si on choisit pour les conditions aux limites inférieure et supérieure :

$$\left[{\overline{p\delta_z \Pi }}^{ z }\right]_1=\frac{p_2 \left(\delta_z \Pi \right)_2}{2}+p_s\left(\delta_z \Pi \right)_1$$ (2.19)

et

$$\left[{\overline{p\delta_z \Pi }}^{ z }\right]_N=\frac{p_{N-1} \left(\delta_z \Pi \right)_{N-1}}{2}$$ (2.20)

On voit alors que l'équation 2.21 peut être satisfaite simplement si on choisit les niveaux $\Pi$ suivant la relation 2.12, cqfd.

... précédentes.^2.7

Un lourd exercice de ``convergence'' LMDZ/LMD-5ter sera mené à bien en 1996. Michèle Forichon y laissera pas mal d'énergie. Plus généralement, les acteurs de l'époque, qui tentent de mener à bien les décisions de Fontevraud, s'useront entre décisions, attaques et contre-décisions.

... dynamique^2.8

Cette séparation entre eau vapeur et liquide est en partie fictive car, juste après l'advection, en entrée des paramétrisations physiques, l'eau liquide est réévaporée pour travailler en eau totale. La séparation entre eau vapeur et eau liquide pour l'advection est donc uniquement numérique.

... conservé^2.9

A noter que cette séparation, très pratique pour le développement des modèles, peut cependant conduire a des problèmes numériques importants. On peut illustrer ce point sur un cas simple : la sédimentation des aérosols sous l'effet de leur poids est souvent traitée, comme les processus chimiques ou microphysiques, séparément de l'advection. Prenons le cas particulier ou le traceur est pris dans une ascendance avec une vitesse égale à la vitesse de sédimentation par rapport à l'air, de sorte que la vitesse réelle des aérosols est nulle. Dans ce cas, si on utilise pour l'advection un schéma en volumes finis comme ceux présentés plus loin, le traceur sera diffusé sur la verticale dans la succsession des mouvements vers le haut et vers le bas alors que l'application du même schéma avec une vitesse verticale nulle aurait évité ce problème.

... variable.^2.10

On fait souvent la confusion entre moyenne d'ensemble et moyenne spatiale ou temporelle. C'est par exemple le glissement qui s'opère ici quand on appuie le développement d'une paramétrisation sous-maille sur un raisonnement en moyenne d'ensemble. En toute rigueur, seule cette dernière peut permuter avec les dérivations spatiales et temporelles, condition indispensable pour les développements présentés ici.

... ^2.11

Ceux qui rechercheraient une présentation plus systématique et mathématique peuvent se reporter à ().

.... ^2.12

Le dernier terme de l'équation est une approximation numérique de $u\delta x/2 \partial^2c/\partial x^2$.

...Russ:81^2.13

La mise à jour de la pente, Eq. (11) de (), est la restriction exacte au cas unidimensionnel et non divergent de l'Eq (23) de ().

...eq:m1d)^2.14

Une certaine confusion sur ce point est entretenue dans la littérature (cf. e. g. , ) dans laquelle ont présente souvent les flux alternés en ne découpant que l'équation de transport des traceurs puis en introduisant comme des astuces numériques les corrections nécessaires qui, si elle sont bien choisies, consistent simplement à découper simultanément l'équation de conservation pour l'air.

... ordre.^2.15

Pour la formulation discrète, on introduit, par exemple pour l'ascendance, les quantités ${\hat{E}}_i\simeq \hat{e}\delta z \delta t$ et ${\hat{D}}_i\simeq \hat{d}\delta z \delta t$ (entraînement et détraînement vers et depuis l'ascendance pour la couche $i$ durant le pas de temps $\delta t$) et ${\hat{F}}_{i+1/2}\simeq \hat{f}\delta t$ (transfert de masse entre les couches $i$ et $i+1$). Les équations du modèle, discrétisées avec des schémas amont et en supposant que la concentration dans l'ascendance et dans la subsidence est en régime stationnaire, s'écrivent

$${\hat{E}}_i+{\hat{F}}_{i-1/2} = {\hat{D}}_i +{\hat{F}}_{i+1/2}$$ (2.66)

$${\check{E}}_i+{\check{F}}_{i+1/2} = {\check{D}}_i +{\check{F}}_{i-1/2}$$ (2.67)

$${\hat{E}}_i c_i+F_{i-1/2}\hat{c}_{i-1} = \hat{c}_i\left({\hat{D}}_i+{\hat{F}}_{i+1/2}\right)$$ (2.68)

$${\check{E}}_i c_i+F_{i+1/2}\check{c}_{i+1} = \check{c}_i\left({\check{D}}_i+{\check{F}}_{i-1/2}\right)$$ (2.69)

où $c_i$, ${\hat{c}}_i$ et ${\check{c}}_i$ sont les concentrations de traceur respectivement dans l'environnement (assimilée à la concentration moyenne dans la maille), l'ascendance et la subsidence. Si on note $m_i$ la masse de la maille $i$ et ${{c^*}_i}$ la concentration de traceur dans la maille $i$ au pas de temps $t+\delta t$, on a

$$m_i {{c^*}_i}- m_i c_i = {\hat{F}}_{i-1/2} \hat{c}_{i-1}-{\hat{F}}_{i+1/2}\hat{c}_{i} +{\hat{F}}_{i+1/2}c_{i+1}-{\hat{F}}_{i-1/2}c_{i}$$ (2.70)

$$+ {\check{F}}_{i+1/2} \check{c}_{i+1}-{\check{F}}_{i-1/2}\check{c}_{i} +{\check{F}}_{i-1/2}c_{i-1}-{\hat{F}}_{i+1/2}c_{i}$$ (2.71)

$$= {\hat{D}}_i \hat{c}_i -{\hat{E}}_i c_i +{\hat{F}}_{i+1/2}c_{i+1}-{\hat{F}}_{i-1/2}c_{i}$$ (2.72)

$$+ {\check{D}}_i \check{c}_i -{\check{E}}_i c_i +{\check{F}}_{i-1/2}c_{i-1}-{\hat{F}}_{i+1/2}c_{i}$$ (2.73)

...#tex2html_wrap_inline4167#^2.16

On a déjà mentionné plus haut, dans la présentation des schémas en volumes finis, le problème posé par le resserrement en longitude du maillage longitude-latitude à l'approche des pôles. Ce problème est contourné dans la partie dynamique du modèle en appliquant à partir de 60 degrés de longitude (en général) un filtre longitudinal qui ne retient, dans les régions polaires, que les échelles plus grandes que les échelles explicitement représentées à 60 degrés.

... icosaédrique,^2.17

L'icosaèdre est le polyèdre régulier le plus proche de la sphère. Il est constitué de 20 triangles équilatéraux, chacun des 12 sommets reliant entre eux 5 de ces triangles. Les triangles de base peuvent ensuite être redécoupés en triangles équilatéraux. Les singularités très fortes des deux pôles sont alors remplacées par 12 singularités beaucoup plus douces. L'icosaèdre produit le maillage le plus régulier possible de la sphère. Il est utilisé en particulier pour construire les géodes. Les icosaèdres sont revenus à la mode ces derniers temps dans le domaine météorologique, notamment parce qu'il permettent d'éviter le filtrage dans les haute latitudes, très pénalisant sur les ordinateurs vectoriels et encore davantage sur les ordinateurs parallèles.