Chapitre 1

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 Joussaume (1990) 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 (Genthon et Armengaud, 1995; Preiss et Genthon, 1997).

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 (Hourdin et al., 1993, 1995a) et Titan (Hourdin, 1992; Hourdin et al., 1995b). 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^{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 Van Leer (1977). 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. Van Leer (1977) 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 Russell et Lerner (1981) (schéma des pentes du NASA/GISS) et Prather (1986). 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 (Hourdin et Armengaud, 1999). 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.

1C’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 (Edouard et al., 1996; Edouard, 1997) que j’ai eu le bonheur de découvrir l’article de Van Leer (1977) 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...

1.1 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}

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 Sadourny et Laval (1984). Ce modèle est plus modulaire et flexible^{3} que 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.

1.1.1 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 ϕ – où a est le rayon de la planète et ϕ la latitude – plutôt que (a+z) cos ϕ, 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.

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

3La flexibilité est ici une notion positive !

Fig. 1.1 – Disposition des variables pour la grille du modèle LMDZ.

1.1.2 Discrétisation des équations primitives

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. Kasahara, 1977). On note X et Y les coordonnées horizontales : X (resp. Y) est une fonction biunivoque de la longitude λ (resp. de la latitude ϕ). Les variables scalaires (la température potentielle θ, le géopotentiel Φ 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. 1.1. De façon à pouvoir modifier la distribution des longitudes et latitudes de la grille, on utilise en fait les composantes covariantes (ũ et ṽ) et contravariantes (˜ũ et ˜ ṽ) du vent définies par

ũ = c_{u}u et ˜ ũ = u/c_{u} avec c_{u} = a cos ϕ (dλ/dX)
ṽ = c_{v}v et ˜ ṽ = v/c_{v} avec c_{v} = a (dϕ/dY) (1.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 σ. La coordonnée σ = 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 σ − 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_{l}p_{s}  (1.2)

On choisit près de la surface A ∼ 0 et B ∼ 1. Vers le haut du modèle, on a A ∼ 0 et B ∼ 0 avec A ≫ Bp_{s}. On retrouve ainsi la coordonnée σ 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 1.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} = A\frac{p_{k} − p_{k+1}}{g}  (1.3)

(A = c_{u}c_{v} est l’aire de la maille) ou encore

m_{k} = \frac{A}{g} [A_{k} − A_{k+1} + (B_{k} − B_{k+1}) p_{s}]  (1.4)

On introduit alors les trois composantes du flux de masse :

U = \overline{m}^{X} ˜ũ, V = \overline{m}^{Y} ˜ ṽ et W  (1.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{∂m}{∂t} + δ_{x}U + δ_{y}V + δ_{z}W = 0  (1.6)

la notation δ_{X} correspondant à la différence entre deux points consécutifs suivant la direction X.^{4}

On ne présente pas ici la discrétisation des équations du mouvements qui n’a pas changé par rapport à la version σ 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 Π = p^{κ}.^{5} L’équilibre hydrostatique est intégré verticalement suivant le schéma :

δ_{z}Φ = −\overline{θ}^{z}δ_{z}Π  (1.11)
Ω_{k} = ∑_{l=k}^{N} (δ_{x}U + δ_{y}V)  (1.7)
\frac{A}{g} \frac{∂p_{s}}{∂t} = Ω_{1}  (1.8)
\frac{∂m}{∂t} = −\frac{A}{g}δ_{z}B \frac{∂p_{s}}{∂t} = −δ_{z}BΩ_{1}  (1.9)

4Pour résoudre cette équation, on commence par calculer la convergence de masse Ω_{k} cumulée depuis le sommet de l’atmosphère jusqu’au niveau k :

La convergence au sol Ω_{1} donne accès à l’évolution de la pression de surface Finalement, en remarquant que on obtient le flux de masse vertical W par intégration depuis le haut de l’atmosphère

−δ_{z}BΩ_{1} − δ_{z}Ω + δ_{z}W = 0  (1.10)

5κ = R/C_{p} où R et C_{p} sont les constantes thermodynamiques de l’air sec.

6  CHAPITRE 1. REPRÉSENTATION DU TRANSPORT DES ESPÈCES TRACES

où Φ est le géopotentiel, avec \overline{θ}_{1}^{z} = θ_{1} et δ_{z}Π_{1} = p_{s}^{κ} − Π_{1}. Au lieu de calculer Π_{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 Π :

\frac{A}{g}\overline{pδ_{z}Π}^{z} = −κΠm  (1.12)

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

∫_{0}^{∞} Φρdz = ∫_{0}^{∞} RT ρdz  (1.13)

entre énergie potentielle et énergie interne, à savoir^{6}

∑_{l=1}^{N} Φ_{l}m_{l} − Φ_{s} = ∑_{l=1}^{N} RT_{l}m_{l} = ∑_{l=1}^{N} κθ_{l}Π_{l}m_{l}  (1.21)

6A partir de l’équation hydrostatique 1.11, on peut écrire

Φ_{l} = Φ_{s} − ∑_{k=1}^{l} [\overline{θ}^{z}δ_{z}Π]_{k}  (1.14)

avec la convention (utilisée dans le modèle) [\overline{θ}^{z}]_{1} = θ_{1}. L’énergie potentielle de gravité totale de la colonne s’écrit donc

∑_{l=1}^{N} Φ_{l}m − Φ_{s} = − ∑_{l=1}^{N} m_{l} ∑_{1≤k≤l} [\overline{θ}^{z}δ_{z}Π]_{k}  (1.15)

= − ∑_{k=1}^{N} ∑_{k≤l≤N} [\overline{θ}^{z}δ_{z}Π]_{k} m_{l}  (1.16)

= − ∑_{k=1}^{N} [\overline{θ}^{z}δ_{z}Π]_{k} × \frac{A}{g}p_{k}  (1.17)
= −\frac{A}{g} ∑_{k=1}^{N} θ_{k} [\overline{pδ_{z}Π}^{z}]_{k}  (1.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 :

[\overline{pδ_{z}Π}^{z}]_{1} = \frac{p_{2} (δ_{z}Π)_{2}}{2} + p_{s} (δ_{z}Π)_{1}  (1.19)

et

[\overline{pδ_{z}Π}^{z}]_{N} = \frac{p_{N−1} (δ_{z}Π)_{N−1}}{2}  (1.20)

On voit alors que l’équation 1.21 peut être satisfaite simplement si on choisit les niveaux Π suivant la relation 1.12, cqfd.

1.1.3 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 (Hourdin et al., 1993). 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.^{7} C’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 (Hourdin et al., 2006).

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 Laval et al. (1981). 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 (Morcrette, 1991). La partie concernant l’absorption et la diffusion du rayonnement solaire est une version raffinée du modèle de Fouquart et Bonnel (1980). La partie infrarouge thermique a été développée par Morcrette et al. (1986). 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 Tiedtke (1989).

7Un 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.

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 (Le Treut et Li, 1991). 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^{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 (Hourdin et al., 1993). Si on note ξ 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∂T/∂z′ et la conduction dans le sol ∂T/∂t = ∂^{2}T/∂z′^{2}. I = \sqrt{ξC} est appelé inertie thermique et z′ est une pseudo profondeur définie par z′ = z\sqrt{C/ξ}.

Suivant Laval et al. (1981), l’évaporation à la surface est calculée comme E = 1, 35βρ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é β 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, dq_{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é β 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.

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’Emanuel (1993). 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 Fig. 1.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.

8Cette 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.

Fig. 1.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. 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. Guichard et al., 2004; Hourdin et al., 2006). Un travail important est réalisé au LMD sur cette paramétrisation d’Emanuel (Grandpeix et al., 2004).

2. le schéma statistique de nuages qui a été couplé à la paramétrisation de la convection suivant l’approche de Bony et Emanuel (2001). 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 (de Rosnay et al., 2002; Krinner et al., 2005). 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 souscolonnes 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. 1.2 et Fig. 1.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.

Fig. 1.4 – Séparation entre les parties dynamique 3D et physique 1D du code LMDZ.

1.1.4 Organisation informatique - version unidimensionnelle - mode guidé

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. 1.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{∂X}{∂t} = M(X) + \frac{X_{a} − X}{τ}  (1.22)

où X est une variable météorologique (u, v, θ, 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 τ 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 (Jeuken et al., 1996).

Une approche intermédiaire a été mise en œuvre dans LMDZ (Bonazzola, 2001). 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.

1.2 L’équation de transport : séparation d’échelles

1.2.1 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é^{9}.

9A 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.

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{dc}{dt} = 0  (1.23)

ou

\frac{∂c}{∂t} + vgrad c = 0  (1.24)

Dans la présentation qui suit, on suppose le champ de vent 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.

1.2.2 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.^{10}

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

Si on note 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{∂c}{∂t} + vgrad c = 0  (1.25)

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

\frac{∂ρ}{∂t} + div(ρv) = 0,  (1.26)

sous sa forme conservative ou flux

\frac{∂ρc}{∂t} + div(ρvc) = 0  (1.27)

10On 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.

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

\overline{ρvc} = \overline{ρ}ṽc̃ + \overline{ρv′c′},  (1.28)

on obtient l’équation

\frac{∂\overline{ρ}c̃}{∂t} + div(\overline{ρ}ṽc̃) + div(\overline{ρv′c′}) = 0  (1.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 ρ, u, v et c pour les variables grande échelle, on obtient finalement

\frac{∂c}{∂t} + v.grad c + \frac{1}{ρ}div(\overline{ρv′c′}) = 0  (1.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 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.

1.3 Le transport grande échelle

1.3.1 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 1.30 ou 1.29 formellement équivalentes aux Eqs 1.25 et 1.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 + δt en remontant à t “le long du vent” : q(x, t + δt) = q(x − δtv, t). Le problème est que, si le point x se trouve sur le maillage, le point x − δtv 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 “isoniveaux” 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 Rood, 1987). 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é.

1.3.2 Les schémas en volumes finis

11 L’intégration de l’équation de transport (1.27) et de l’équation de conservation de la masse (1.26) sur un polyèdre (volume de contrôle) à N faces conduit aux formulations :

\frac{∂m}{∂t} = ∑_{n=1}^{N} U_{n}  (1.31)

et

\frac{∂cm}{∂t} = ∑_{n=1}^{N} F_{n}  (1.32)

où m est la masse totale d’air dans le polyèdre, c est la concentration massique moyenne du traceur dans le volume, U_{n} et 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 Godunov (1959), F est simplement estimé comme le produit de 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

11Ceux qui rechercheraient une présentation plus systématique et mathématique peuvent se reporter à Roux

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.

Van Leer (1977) 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 Russell et Lerner (1981) – le schéma III dans la classification de Van Leer – et Prather (1986) – 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.

1.3.3 Description des schémas en dimension un

En dimension 1 et après intégration sur un pas de temps, les Eq. 1.31 et 1.32 s’écrivent simplement

δ_{t}(m)_{i} = U_{i−1/2} − U_{i+1/2}  (1.33)

et

δ_{t}(cm)_{i} = F_{i−1/2} − F_{i+1/2}  (1.34)

où les indices correspondent à la position sur l’axe des x, δ_{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} = ∫_{Δx} ρ(x, t)dx  (1.35)

c_{i} = ∫_{Δx} cρ(x, t)dx/ ∫_{Δx} ρ(x, t)dx  (1.36)
U_{i+1/2} = ∫_{t}^{t+Δt} ρ(x_{i+1/2}, t)u(x_{i+1/2}, t)dt  (1.37)

et

F_{i+1/2} = ∫_{t}^{t+Δt} ρ(x_{i+1/2}, t)c(x_{i+1/2}, t)u(x_{i+1/2}, t)dt  (1.38)

Fig. 1.5 – Principe du schéma de Van Leer et notations. On montre le cas où la distribution sousmaille 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.

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 c̆_{i+1/2} dans l’air qui traverse l’interface au cours du 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 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. 1.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}^{∗} = (c_{i}m_{i} + U_{i−1/2}c̆_{i−1/2} − U_{i+1/2}c̆_{i+1/2}) /m_{i}^{∗}  (1.39)

et

m_{i}^{∗} = m_{i} + U_{i−1/2} − U_{i+1/2}  (1.40)

Pour des valeurs positives de U_{i−1/2} et U_{i+1/2}, l’Eq. 1.39 peut se récrire

c_{i}^{∗} = [U_{i−1/2}c̆_{i−1/2} + (m_{i} − U_{i+1/2}) c̆_{i}] /m_{i}^{∗}  (1.41)

où c̆_{i} est la concentration massique moyenne dans l’air qui reste dans le volume i pendant le pas de temps (cf. Fig. 1.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 c̆ (les croix sur la Fig. 1.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 Van Leer (1977) 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 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} ≤ c_{i} ≤ c_{i+1} et U ≥ 0. La condition ci-dessus s’écrit simplement

c_{j} ≤ c̆_{j+1/2} ≤ c_{j+1} pour j = i − 1 et j = i  (1.42)

et

c_{i−1} ≤ c̆_{i} ≤ c_{i}.  (1.43)

En introduisant la condition 1.42 dans l’Eq. 1.39, on assure que c_{i}^{∗} ≤ c_{i} tandis que l’Eq. 1.41 et les conditions 1.42 et 1.43 garantissent que c_{i−1} ≤ 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 Van Leer (1977), cette proposition implique la monotonie.

Si on rajoute en plus la condition que 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. Rood, 1987).

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 c̆_{i−1/2} et c̆_{i+1/2} par une valeur constante dans l’Eq. 1.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.

1.3.4 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). c̆ est alors simplement la valeur de c dans la maille amont (c̆_{i+1/2} = c_{i} si U_{i+1/2} > 0 et c_{i+1} sinon). Ce schéma proposé à l’origine par Godunov (1959) 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δx/(mδt) =cste à une dimension), le schéma 1.39 s’écrit

c_{i}^{∗} = c_{i} + \frac{U}{m}c_{i−1} − \frac{U}{m}c_{i}  (1.44)

et peut se récrire

\frac{c_{i}^{∗} − c_{i}}{δt} = u\frac{c_{i−1} − c_{i}}{δx}  (1.45)
= u\frac{c_{i−1} − c_{i+1}}{2δx} + \frac{uδx}{2} \frac{c_{i−1} − 2c_{i} + c_{i+1}}{δx^{2}}  (1.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δx/2.^{12}

1.3.5 Schémas du second ordre

Pour passer à un ordre supérieur, on suppose que la distribution sous-maille est linéaire avec une pente (δc)_{i} donnant les valeurs aux bords gauche et droit de la maille i comme c_{i}^{±} = c_{i} ± (δc)_{i}/2 (voir l’illustration de la Fig. 1.5).

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

c̆_{i+1/2} = c_{i} + \frac{1}{2} (1 − \frac{U_{i+1/2}}{m_{i}} ) (δc)_{i} si U_{i+1/2} > 0  (1.47)
= c_{i+1} − \frac{1}{2} (1 + \frac{U_{i+1/2}}{m_{i+1}} ) (δc)_{i+1} , sinon.  (1.48)

Différents schémas ont été proposés (Van Leer, 1977, 1979) correspondant à différentes estimations de (δc), dont deux schémas particulièrement intéressants, l’un pour son faible coût et l’autre pour sa précision.

Dans le schéma I de Van Leer, la pente est simplement estimée à chaque pas de temps par différences finies (δc_{i} = (c_{i+1} − c_{i−1}) /2) comme l’illustre le graphique a de la Fig. 1.6. Ce schéma peut être considéré comme une version volumes-finis du schéma de From (1968).

Dans le second schéma, la pente au pas de temps t + δ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. 1.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 Russell et Lerner (1981)^{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. 1.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.

12Le dernier terme de l’équation est une approximation numérique de uδx/2 ∂^{2}c/∂x^{2}.

13La mise à jour de la pente, Eq. (11) de Van Leer (1977), est la restriction exacte au cas unidimensionnel et non divergent de l’Eq (23) de Russell et Lerner (1981).

Fig. 1.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 Russell et Lerner (1981) 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.

1.3.6 Limiteurs de pentes

L’un des intérêts principaux du travail de Van Leer (1977) 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 1.3.3.

Pour les schémas du second ordre, Van Leer (1977) 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

(δc)_{i} = sign(c_{i+1} − c_{i}) ×  (1.49)
min(\frac{|c_{i+1} − c_{i−1}|}{2}, 2|c_{i+1} − c_{i}|, 2|c_{i} − c_{i−1}|)

si c_{i} est compris entre c_{i−1} et c_{i+1} et 0 sinon.

Fig. 1.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. 1.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. 1.42 et Eq. 1.43), qui correspondent respectivement aux Eq. (66) et (74) données par Van Leer (1977) est illustré sur la Fig. 1.7a. Un test numérique de l’impact de ces limiteurs est montré sur les Fig. 1.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. 1.49 qui consiste à utiliser la moyenne géométrique des deux pentes voisines

(δc)_{i} = \frac{2 (c_{i+1} − c_{i}) (c_{i} − c_{i−1})}{c_{i+1} − c_{i−1}}  (1.50)

quand (c_{i+1} − c_{i}) (c_{i} − c_{i−1}) > 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.

1.3.7 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. Woodward et Colella (1981) 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 à Woodward et Colella, 1984; Colella et Woodward, 1984). 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 Carpenter et al. (1990), Lin et Rood (1996) ou Vautard et al. (2001). 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 Prather (1986). 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.

1.3.8 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. 1.33 ou l’Eq. 1.6 donnée plus loin.

Flux alternés

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. 1.8 pour une configuration bidimensionnelle. Si on considère un pas de temps d’advection, le calcul direct ne tient pas Fig. 1.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.

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. 1.39) mais aussi l’advection de transport de l’air (Eq. 1.40)^{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 (Lin et Rood, 1996, montrent que ce n’est pas garanti par toutes les formulations en flux alternés).

14Une certaine confusion sur ce point est entretenue dans la littérature (cf. e. g. Carpenter et al., 1990) 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.

Cette approche en flux alternés est utilisée assez systématiquement dans les calculs d’advection en volumes finis (Van Leer, 1979; Allen et al., 1991; Russell et Lerner, 1981; Prather, 1986).

On utilise ici la séquence proposée par Russell et Lerner (1981) :

Direction pas de temps
X  δt/2
Y  δt/2
Z  δt
Y  δt/2
X  δt/2

Sauter des mailles

Si on veut ne travailler qu’avec des nombres de Courant (δt u/δ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

∑_{j=i−n+1}^{i} m_{j} < U_{i+1/2} ≤ ∑_{j=i−n}^{i} m_{j}  (1.51)

on prend

F_{i+1/2} = ∑_{j=i−n+1}^{i} m_{j}c_{j} + (U_{i+1/2} − ∑_{j=i−n+1}^{i} m_{j}) c̆_{i−n+1/2}  (1.52)

avec

c̆_{i−n+1/2} = c_{i−n} + \frac{1}{2} (1 − \frac{U_{i+1/2} − ∑_{j=i−n+1}^{i} m_{j}}{m_{i−n}} ) (δc)_{i}  (1.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 Lin et Rood (1996) à 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. Yabe et al., 2001).

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. 1.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 isolignes 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).

Fig. 1.9 – Test numérique d’advection transpolaire. Suivant Williamson et Rasch (1989) et Allen et al. (1991), la distribution initiale de traceur est donnée par c(λ, ϕ) = (1 + cos[min(r[λ, ϕ]/R, 1)])/2 avec r = arccos(cos λ cos ϕ) et R = 7 × (2π)/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.

1.3.9 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 Prather, 1986, 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 (λ, ϕ), les composantes zonale (u) et méridienne (v) du vent satisfont la relation suivante :

\frac{∂u}{∂λ} + \frac{∂v cos ϕ}{∂ϕ} = 0  (1.54)

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

Ψ = aU_{0} cos^{2} \frac{λ}{2} cos^{2} ϕ  (1.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{∂Ψ}{∂ϕ} = −2U_{0} cos ϕ sin ϕ cos^{2} \frac{λ}{2}  (1.56)
v = \frac{1}{a cos ϕ} \frac{∂Ψ}{∂λ} = −U_{0} cos ϕ cos \frac{λ}{2} sin \frac{λ}{2}  (1.57)

En suivant une trajectoire (valeur constante de Ψ), la vitesse méridienne v = adϕ/dt peut être récrite en combinant les équations 1.55 et 1.57 :

adϕ/dt = −\sqrt{ΨU_{0}/a} signλ\sqrt{1 − cos^{2}(λ/2)}  (1.58)

Si on introduit α = \sqrt{ΨU_{0}/a^{3}} et qu’on remplace cos^{2}(λ/2) selon l’Eq. 1.55, cette équation devient

−α ∫ signλdt = ∫ \frac{cos ϕdϕ}{\sqrt{cos^{2} ϕ − Ψ/(aU_{0})}} = [arcsin (\frac{sin ϕ}{\sqrt{1 − Ψ/(aU_{0})}})]  (1.59)

Les trajectoires, iso-valeurs de Ψ, 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π/α = 2πa/(U_{0} cos ϕ cos λ/2). Ce temps est plus court au centre du domaine et infini sur les bords.

Résultats numériques

Sur le graphique en haut à gauche de la Fig. 1.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 × 60 (résolution pleine), 60 × 30 (résolution 1/2) et 40×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.

Fig. 1.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 × 60, 60 × 30 et 40 × 20 points.

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.

1.3.10 Remarques pour conclure

Les schémas en volumes finis proposés par Van Leer (1977) 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. Roux, 2002).

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é.

32  CHAPITRE 1. REPRÉSENTATION DU TRANSPORT DES ESPÈCES TRACES

Bibliographie

Allen, D. J., A. R. Douglass, R. B. Rood, et P. D. Guthrie, 1991 : Application of a monotonic upstream-biased transport scheme to three dimensional constituent transport calculations, Mon. Wea. Rev., 119, 2456–2464.

Bonazzola, M., Analyse de la convection et de la circulation atmosphérique de grande échelle pendant la campagne INDOEX, thèse, Université Paris-6, 2001.

Bony, S., et K. A. Emanuel, 2001 : A parameterization of the cloudiness associated with cumulus convection ; evaluation using TOGA COARE data, J. Atmos. Sci., 58, 3158–3183.

Carpenter, R. L., K. K. Droegemeier, P. R. Woodward, et C. E. Hane, 1990 : Application of the piecewise parabolic method (PPM) to meteorological modeling, Mon. Wea. Rev., 118, 586– 612.

Colella, P., et P. R. Woodward, 1984 : The piecewise parabolic method (PPM) for gaz-dynamical simulations, J. Computational Phys., 54, 174–201.

de Rosnay, P., J. Polcher, M. Bruen, et K. Laval, 2002 : Impact of a physically based soil water flow and soil-plant interaction representation for modeling large scale land surface processes, J. Geophys. Res., 107, 10.1029/2001JD000634.

Edouard, S., B. Legras, et V. Zeitlin, 1996 : The effect of dynamical mixing in a simple model of the ozone hole, J. Geophys. Res., 101, 16771–16779.

Edouard, S., Dynamique et chimie du vortex polaire : le problème de l’ozone, Thèse, Université Paris-6, 1997.

Emanuel, K. A., 1993 : A cumulus representation based on the episodic mixing model : the importance of mixing and microphysics in predicting humidity, A.M.S. Meteorol. Monographs, 24, 185–192.

Fouquart, Y., et B. Bonnel, 1980 : Computations of solar heating of the Earth’s atmosphere : A new parametrization, Contrib. Atmos. Phys., 53, 35–62.

From, J. E., 1968 : A method for reducing dispersion in convective difference schemes, J. Computational Phys., 3, 176–189.

Genthon, C., et A. Armengaud, 1995 : Radon 222 as a comparative tracer of transport and mixing in two general circulation models of the atmosphere, J. Geophys. Res., 100, 2849–2866.

Godunov, S. K., 1959 : Finite-difference methods for the numerical computations of equations of gas dynamics, Math. Sb., 7, 271–290.

Grandpeix, J. Y., V. Phillips, et R. Tailleux, 2004 : Improved mixing representation in Emanuel’s convection scheme, Q. J. R. Meteorol. Soc., 130, 3207–3222.

Guichard, F., J. C. Petch, J. L. Redelsperger, P. Bechtold, J. P. Chaboureau, S. Cheinet, W. Grabowski, H. Grenier, C. G. Jones, M. Köhler, J. M. Piriou, R. Tailleux, et M. Tomasini, 2004 : Modelling the diurnal cycle of deep precipitating convection over land with cloudresolving models and single-column models, Q. J. R. Meteorol. Soc., 130, 3139–3172.

Hourdin, F., P. Le Van, F. Forget, et O. Talagrand, 1993 : Meteorological variability and the annual surface pressure cycle on Mars, J. Atmos. Sci., 50, 3625–3640.

Hourdin, F., F. Forget, et O. Talagrand, 1995a : The sensitivity of the Martian surface pressure to various parameters : A comparison between numerical simulations and Viking observations, J. Geophys. Res., 100, 5501–5523.

Hourdin, F., O. Talagrand, R. Sadourny, C. Régis, D. Gautier, et C. P. McKay, 1995b : General circulation of the atmosphere of Titan, Icarus, 117, 358–374.

Hourdin, F., I. Musat, S. Bony, P. Braconnot, F. Codron, J.-L. Dufresne, L. Fairhead, M.-A. Filiberti, P. Friedlingstein, J.-Y. Grandpeix, G. Krinner, P. Levan, Z.-X. Li, et F. Lott, 2006 : The LMDZ4 general circulation model : climate performance and sensitivity to parametrized physics with emphasis on tropical convection, Climate Dynamics, 27, 787–813.

Hourdin, F., et A. Armengaud, 1999 : The use of finite-volume methods for atmospheric advection of trace species. part i : Test of vairious formulations in a general circulation model, Mon. Wea. Rev., 127, 822–837.

Hourdin, F., Conservation du moment angulaire dans le modèle de circulation générale du LMD, Note Interne 175, Lab. de Météorol Dyn., Cent. Nat. de la Rech. Sci., Paris, 1992.

Jeuken, A. B. M., P. C. Siegmund, H. L. C., J. Feichter, et L. Bengtsson, 1996 : On the potential of assimilating meteorological analyses in a global climate model for the purpose of model validation, J. Geophys. Res., 101, 939.

Joussaume, S., 1990 : Three-dimensional simulations od the atmospheric cycle of desert dust particles using a general circulation model, J. Geophys. Res., 95, 1909–1941.

Kasahara, A., Computational aspects of numerical models for weather prediction and climate simulation, in Methods in computational physics, edited by J. Chang, vol. 17, 1–66, Academic press, inc., 1977.

Krinner, G., N. Viovy, N. de Noblet-Ducoudré, J. Ogée, J. Polcher, P. Friedlingstein, P. Ciais, S. Sitch, et C. Prentice, 2005 : A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system, Glob. Biogeochem. Cyc., 19, GB1015, doi :10.1029/2003GB002199.

Laval, K., R. Sadourny, et Y. Serafini, 1981 : Land surface processes in a simplified general circulation model, Geophys. Astrophys. Fluid Dyn., 17, 129–150.

Le Treut, H., et Z. X. Li, 1991 : Sensitivity of an atmospheric general circulation model to prescribed SST changes : Feedback effects associated with the simulation of cloud optical properties., Clim. Dyn., 5, 175–187.

Lin, S.-J., et R. B. Rood, 1996 : Multi-dimensional flux form semi-lagrangian transport schemes, Mon. Wea. Rev., 124, 2046–2070.

Morcrette, J. J., L. Smith, et Y. Fouquart, 1986 : Pressure and temperature dependence of the absorption in longwave radiation parametrizations, Contrib. Atmos. Phys., 59, 455–469.

Morcrette, J., 1991 : Radiation and cloud radiative properties in the European Centre for Medium Range Weather Forecasts forecasting system, J. Geophys. Res., 96, 9121–9132.

Prather, M. J., 1986 : Numerical advection by conservation of second order moments, J. Geophys. Res., 91, 6671–6681.

Preiss, N., et C. Genthon, 1997 : Use of a new database of lead 210 for global aerosol model validation, J. Geophys. Res., 102, 25347–25358.

Rood, R. B., 1987 : Numerical advection algorithms and their role in atmospheric transport and chemistry models, Rev. Geophys., 25, 71–100.

Roux, J., http ://gershwin.ens.fr/houches2002/Cours/Roux/leshouches.pdf, 2002.

Russell, G. L., et J. A. Lerner, 1981 : A new finite-differencing scheme for the tracer transport equation, J. Appl. Met., 20, 1483–1498.

Sadourny, R., et K. Laval, January and July performance of the LMD general circulation model, in New perspectives in Climate Modeling, edited by A. Berger et C. Nicolis, Elsevier, 173–197, Amsterdam, 1984.

Tiedtke, M., 1989 : A comprehensive mass flux scheme for cumulus parameterization in largescale models, Mon. Wea. Rev., 117, 1179–1800.

Van Leer, B., 1977 : Towards the ultimate conservative difference scheme : IV. a new approach to numerical convection, J. Computational Phys., 23, 276–299.

Van Leer, B., 1979 : Towards the ultimate conservative difference scheme. V. a second-order sequel to Godunov’s method, J. Computational Phys., 32, 101–136.

Vautard, R., M. Beekmann, J. Roux, et D. Gombert, 2001 : Validation of a hybrid forecasting system for the ozone concentrations over Paris area, Atmosph. Environ., 35, 2449–2461.

Williamson, D., et P. J. Rasch, 1989 : Two dimensional semi-lagrangian transport with shapepreserving interpolation, Mon. Wea. Rev., 117, 102–127.

Woodward, P. R., et P. Colella, 1981 : High resolution difference schemes for compressible gaz dynamics, Lect. Notes Phys., 141, 434–441.

Woodward, P. R., et P. Colella, 1984 : The numerical simulation of two-dimensional fluid flow with strong shocks, J. Computational Phys., 54, 115–173.

Yabe, T., R. Tanaka, T. Nakamura, et F. Xiao, 2001 : An exactly conservative semi lagrangian scheme (CIP-CSL) in one dimension, Mon. Wea. Rev., 129, 332–344.