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. RELU JUSQUE PAGE 110 INCLU

Frédéric Hourdin

Mémoire présenté pour obtenir une Habilitation à Diriger des Recherches

soutenu le 6 avril 2005

devant un jury composé de:

M. Michel CABANE Président

M. Pierre DROSSART Rapporteur

M. Ray PIERREHUMBERT Rapporteur

M. Patrick MASCART Rapporteur

M. Alain HAUCHECORNE Examinateur

M. Daniel GAUTIER Examinateur

M. Olivier TALAGRAND Examinateur

Travaux réalisés :

au Laboratoire de Météorologie Dynamique du CNRS

Institut Pierre Simon Laplace

NOLOGO

Table des matières

Merci ...

Un grand merci d'abord aux étudiants que j'ai eu le plaisir d'encadrer ces dernières années et qui ont contribué à des titres divers à l'avancée des travaux présentés ici. Plutôt que de les nommer dans ce préambule, je les ai mentionnés au fil des chapitres. Beaucoup de ces étudiants ont depuis trouvé leur place, dans la recherche ou ailleurs, mais quelques uns sont venus malheureusement grossir durablement les bataillons de la précarité que nous avons laissé prospérer autour de nos instituts de recherche.

Merci ensuite aux collègues qui gravitent autour du noyau dur du "LMD-Jussieu" et avec qui nous avons construit au fil des ans un environnement de travail sérieux et néanmoins fort sympathique sans lequel ce travail n'aurait pas été possible. Quand je pense à cet environnement, me vient souvent l'image d'une niche écologique où compétences (depuis la physique fondamentale jusqu'au secrétariat), intelligence et originalité s'additionnent dans la bonne humeur, à l'abri du carriérisme et de la surenchère programmatique et technocratique qui semblent le lot quotidien de la recherche "moderne" ; où l'on sait encore qu'avant d'aller jouer les moulins à vent de mitinges en projets européens, il faut avoir du grain à moudre et donc commencer par labourer, semer, moissonner. Alors, à vous tous qui vous reconnaîtrez, un immense merci.

Je veux enfin remercier ceux des responsables qui m'ont épaulé à un moment ou un autre au cours de ces travaux. Comme ils auront peut-être plus de difficulté à se reconnaître, je citerai en particulier Olivier Talagrand, Jean-Paul Huot, Gérard Mégie, Sylvie Joussaume, Hervé Le Treut et Claude Basdevant. Merci enfin aux voyageurs, proches ou lointains, qui ont accepté de participer au jury.

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, 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, 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 de l'atmosphère (Dufresne, 2002). 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 (Bopp et al., 2004).

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 Pierrehumbert, 1998), 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, 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' $\oe$ il d'un observateur extérieur. La circulation est alors profondément modifiée (Hourdin et al., 1993). 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 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 (Forget et al., 1998; Montmessin, 2004; Hourdin et al., 1995a,1993). Pour Titan, les couplages avec la photochimie et la brume (Lebonnois et al., 2003; Rannou et al., 2002; Hourdin, 2004) ou le méthane (Tokano, 2001) 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 $\oe$ uvre dans le modèle LMDZ des schémas en volumes finis développés à l'origine par Van Leer (1977). 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 .

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 2. 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 et 2 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 (Hourdin et al., 1995b) 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 (Toublanc et al., 1995) et la microphysique des brumes (Cabane et al., 1992; Rannou et al., 1995). 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 (Hourdin et al., 1998) : ``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 (Rannou, 2004). Le modèle a également permis d'expliquer les contrastes latitudinaux observés dans la composition chimique (Hourdin, 2004; Lebonnois et al., 2001). 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 (Rannou et al., 2005). 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.

Le ``modèle du thermique"

Le transport turbulent dans la couche limite est essentiel à la fois bien sûr pour la météorologie et le climat, mais également pour toutes les espèces traces ayant des sources et puits en surface (vapeur d'eau, CO, espèces chimiques et aérosols de pollution, aérosols désertiques). Les évolutions récentes des paramétrisations des modèles de climat se sont portées davantage sur la représentation de la convection nuageuse profonde que sur la couche limite. Pour modéliser la convection nuageuse, les modèles récents utilisent en particulier souvent des schémas dits "en flux de masse" dans lesquels on explicite des flux ascendants, souvent intenses et concentrés, des descentes également concentrées associées aux précipitations, et des flux compensatoires plus lents dans l'environnement.

Comparativement, les paramétrisations utilisées dans les modèles grande échelle pour la couche limite restent souvent rudimentaires. Elles sont la plupart du temps basées sur des adaptations de formulations locales. On entend par là que le flux vertical turbulent d'une quantité ne dépend que des caractéristiques locales de l'écoulement. Cette approche locale s'appuie sur une analogie avec la diffusion moléculaire, les mouvements turbulents aléatoires jouant le rôle des mouvements browniens des molécules pour la diffusion moléculaire. Ce flux s'exprime alors comme le produit du gradient local de par un coefficient de mélange turbulent ne dépendant lui-même que des conditions météorologiques locales. On parle de diffusion turbulente ou super-viscosité. On sait depuis longtemps que cette vision diffuse de la turbulence ne permet pas de rendre compte d'un certain nombre de phénomènes atmosphériques et notamment du transport de chaleur en remontant le gradient (du froid vers le chaud en termes de température potentielle) très souvent observé dans la couche limite convective. Ce transport à contre-gradient est effectué en fait par des structures organisées méso-échelles qui emmènent directement l'air chaud de la couche de surface vers le haut de la couche limite. Ce sont ces structures - thermiques, pompes, rouleaux convectifs - qui font la joie des pilotes de planeurs et autres engins volants.

Depuis longtemps, cette difficulté est contournée dans les modèles de circulation générale en utilisant un "contre-gradient" (Deardorff, 1972b): au lieu d'utiliser directement le gradient vertical de la température potentielle $\theta$ dans le calcul du flux de chaleur, on soustrait à ce gradient une valeur qui permet d'avoir un flux de chaleur vers le haut, même dans une atmosphère légèrement stable. Des développements récents (Troen et Mahrt, 1986; Abdella et McFarlane, 1997; Holtslag et Boville, 1993) on permis de dériver des expressions moins arbitraires pour le contre-gradient, en prenant en compte indirectement l'existence de ces structures méso-échelles. Ces formulations permettent également d'introduire cette composante non-locale pour le transport des espèces traces ce qui n'étaient pas le cas dans la version originale de Deardorff (1972b).

Stull (1984) avait souligné l'importance des aspects non locaux du transport vertical dans la couche limite et proposé un formalisme général basé sur des matrices d'échange (baptisées matrices de ``transilience'') pour traiter ce problème dans les modèles numériques, afin de rompre de façon radicale d'avec la diffusion turbulente. Sur la base des matrices de transiliences, Pleim et Chang (1992) ont proposé un ``modèle de convection asymétrique" basé sur l'image d'un transport par une ascendance concentrée et une subsidence compensatoire lente. Le caractère non local verticalement du mélange, et la dissymétrie entre ascendances convectives concentrées (et descentes précipitantes pour les cumulo-nimbus et congénères) et subsidences compensatoires sont à la base des schémas dit ``en flux de masse" (Arakawa et Schubert, 1974; Emanuel, 1991; Tiedtke, 1989) qui se sont largement répandus dans les modèles de circulation générale. Ces idées de modèles ``en flux de masse" ont également été appliquées pour la couche limite convective (notamment pour les cumulus d'alizés) mais en utilisant généralement des modèles dits de "couche mélangée" (``mixed layer formulation'' ou aussi ``bulk models'' en anglais).

Le "modèle du thermique" proposé ici s'inspire plus directement des modèles de la convection nuageuse, en reprenant notamment, comme dans le schéma d'Emanuel, l'idée d'un panache ascendant non mélangé, épluché progressivement au cours de son ascendance. Sans entrer dans le détail, on détermine pour chaque couche instable (surmontée par de l'air plus froid en termes de température potentielle virtuelle), un profil vertical de vitesse ascendante à partir de la flottabilité d'une parcelle d'air entraînée depuis cette couche en conservant sa température potentielle. Ces calculs sont ensuite utilisés pour décrire une ascendance (le thermique) alimentée en air par les couches instables près de la surface, et compensée par une subsidence plus lente dans l'environnement. Ce schéma tient compte de la structure géométrique de ces cellules convectives, notamment pour relier vitesses verticales et flux de masse.

Le modèle du thermique n'est pas un modèle de couche limite complet. Il ne représente que la partie meso-échelle de la dynamique de la couche limite. On conserve dans le modèle de circulation générale, en parallèle du modèle du thermique, une formulation en diffusion turbulente, active notamment dans la couche limite de surface. Il faut donc considérer qu'on a rajouté dans le modèle climatique une échelle entre l'échelle turbulente et l'échelle de la convection nuageuse.

Dans ce chapitre, on revient assez largement sur les approches classiques de la paramétrisation de la couche limite (Section 2.1). On présente aussi les spécificités de la couche limite convective (Section 2.2) et les approches qui ont été proposées pour la représenter (Section 2.3). On présente ensuite en détail le modèle du thermique (Section 2.4). Ce modèle est enfin testé, aussi bien pour les aspects météorologiques - en utilisant à la fois les résultats de simulations des grands tourbillons (Section 2.5) et des mesures in-situ (Section 2.6) - que pour le transport des espèces traces (Section 2.7).

Nous avions commencé avec Richard Fournier à nous intéresser à ces questions de couche limite sur Mars (en adaptant notamment dans le modèle de circulation martien le modèle de couche limite de Mellor et Yamada). Fleur Couvreux, Camille Risi et Catherine Rio ont contribué lors de différents stages à développer, affiner ou valider la paramétrisation. Abderrahmane Idelkadi a effectué des comparaisons systématiques sur le transport du Radon. Enfin les discussions avec Alain Lahellec, Anne Mathieu et Jean-Yves Grandpeix ont beaucoup contribué au mûrissement des idées.

Bases physiques des paramétrisations en diffusion

Comme on l'a dit plus haut, les paramétrisations en diffusion turbulente sont construites par analogie avec la diffusion moléculaire. Cette approche s'est révélée particulièrement fructueuse pour expliquer certaines caractéristiques de la couche limite atmosphérique et dériver des paramétrisations pour les modèles de circulation atmosphérique. Avant de parler de mise en défaut de cette théorie dans le cas des couches limites convectives, on retrace les grandes lignes des différentes approches en diffusion turbulente. Cette section permet également de décrire des paramétrisations qui sont par la suite comparées, sur des cas académiques ou réalistes de couches limites convectives, au modèle du thermique. On introduit les différentes formulations sous des hypothèses classiques de couche limite qui supposent que les quantités moyennes varient moins vite horizontalement que verticalement ( $\vert\partial \tilde{X}/\partial x\vert \ll \vert\partial \tilde{X}/\partial z\vert$ et $\vert\partial^2\tilde{X}/\partial x^2\vert \ll \vert\partial^2\tilde{X}/\partial z^2\vert$ où $\tilde{X}$ est la moyenne d'ensemble de pondérée par l'air définie dans la Section ).

La longueur de mélange

Dans le chapitre précédent, on a introduit rapidement l'utilisation de la diffusion turbulente pour paramétriser le transport turbulent d'une quantité dans la couche limite

$\begin{displaymath} \overline{\rho w' c'}=-\rho K_z \frac{\partial c}{\partial z} \end{displaymath}$

(2.1)

Une façon physique d'introduire cette formulation et d'estimer le coefficient est l'approche de la longueur de mélange introduite à l'origine par Prandtl en 1925. L'image physique sous-tendue est l'existence de petites structures turbulentes avec une taille caractéristique, ou longueur de mélange, .

En se plaçant pour fixer les idées dans une configuration où la concentration moyenne du traceur croît verticalement ( $\partial\overline{c}/\partial z>0$ ), un mouvement descendant sera associé en moyenne à une fluctuation positive de . La fluctuation associée à ce mouvement turbulent sera d'autant plus grande que les contrastes verticaux de sont grands. En supposant que la particule qui descend a conservé les propriétés qu'avait l'air à une distance au-dessus, et en supposant petit devant la hauteur caractéristique des variations de $\overline{c}$ (c'est cette dernière hypothèse qui n'est pas valide dans les cas de couches limites convectives) on peut écrire

$\begin{displaymath} c'=l \frac{\partial c}{\partial z} \end{displaymath}$

(2.2)

De même, les mouvement turbulents ascendants seront associés en moyenne à une fluctuation négative de

de sorte qu'il y a dans les deux cas une corrélation négative entre

. On peut finalement écrire

$\begin{displaymath} \overline{w'c'}=-l\overline{\vert w'\vert} \frac{\partial c}{\partial z} \end{displaymath}$

(2.3)

Cette équation peut être prise comme définition de la longueur de mélange

Dans ce cadre, le coefficient est simplement $l\vert w'\vert$ .

Si on suppose que la turbulence est isotrope (par exemple en atmosphère neutre et loin du sol), on peut aller un cran plus loin en remarquant qu'on a

$\begin{displaymath} \vert w'\vert\simeq\vert u'\vert\simeq l\frac{\partial u}{\partial z} \end{displaymath}$

(2.4)

ce qui conduit à choisir

$\begin{displaymath} K_z=l^2\vert\vert\frac{\partial {\bf v}}{\partial z}\vert\vert \end{displaymath}$

(2.5)

La théorie de la diffusion turbulente a remporté un de ses plus francs succès dans l'explication de la structure de la couche limite de surface.

Près de la surface, on suppose que la longueur caractéristique des échanges turbulents est proportionnelle à la distance à la surface, $l=\kappa z$ . En choisissant un repère local tel que $\overline{u}$ soit dans la direction du vent moyen près de la surface, et en remarquant que $\overline{u}$ (noté ci-après) et son gradient en sont du même signe près de la surface - puisque $\overline{u}$ doit s'annuler en -, il vient

$\begin{displaymath} \overline{w'u'}=-K_z \frac{\partial u}{\partial z}=-{\left(\kappa z\right)}^2{\left(\frac{\partial u}{\partial z}\right)}^2 \end{displaymath}$

(2.6)

Si on définit la couche de surface comme la partie de la couche limite dans laquelle les flux turbulents ne diffèrent pas trop (par moins de 10% par exemple) du flux en surface - c'est à dire qu'on peut écrire $\overline{w'u'}\simeq\overline{u'v'}_0=-\tau/\rho$ où $\tau$ est le module de la tension de vent en surface - on obtient

$\begin{displaymath} u_*=\kappa z \frac{\partial u}{\partial z} \end{displaymath}$

(2.7)

où $u_*^2=\tau/\rho$ est la vitesse de friction. A noter que les dérivations ci-dessus aboutissent également à $\vert w'\vert\simeq u_*$ dans la couche de surface.

Les mesures expérimentales de la constante de Von Karman, $\kappa$ , supposée universelle, donnent des valeurs comprises entre 0,35 et 0,43.

Sous ces hypothèse, près de la surface, le vent varie de façon logarithmique avec la verticale. La singularité en surface est résolue en supposant que le vent s'annule non pas en mais à une altitude , appelée longueur de rugosité, telle que

$\begin{displaymath} \overline{u}(z)=\frac{u_*}{\kappa }\ln{\frac{z}{z_0}} \end{displaymath}$

(2.8)

Physiquement, les dérivations précédentes ne sont plus valables très près de la surface, région dans laquelle l'atmosphère échange de la quantité de mouvement au travers de couples de pression sur les obstacles. La longueur de rugosité est typiquement de l'ordre de la fraction de mm sur mer, de quelques centimètres sur les prairies ou les déserts caillouteux, jusqu'à quelques mètres dans les régions boisées ou urbaines.

De la même façon, à partir du flux de chaleur sensible en surface $C_p\rho\overline{w'\theta'}_0$ , on peut introduire une échelle des fluctuations turbulentes de température potentielle, $\theta^*=\overline{w'\theta'}_0/u^*$ , reliée au gradient de température potentielle par

$\begin{displaymath} \theta^*=\kappa ' z \frac{\partial \theta}{\partial z} \end{displaymath}$

(2.9)

ce qui aboutit également à une forme logarithmique

$\begin{displaymath} \theta-\theta_S=\frac{\theta^*}{\kappa '}\ln{\frac{z}{z_0}} \end{displaymath}$

(2.10)

où $\theta_S$ est la température de surface. Le rapport $R=\kappa /\kappa '$ est le nombre de Prandtl turbulent, rapport entre les diffusivités de la quantité de mouvement et de la température. Ce nombre est de l'ordre de 1 pour les gaz. Une valeur de 0,7 est couramment utilisée pour l'air. Les mesures dans la couche limite de surface dans les grandes plaines américaines (Businger et al., 1971) suggèrent

en conditions neutres ou stables (Deardorff, 1972a). La hauteur

peut également être différente pour la quantité de mouvement et la température. Dans ce dernier cas,

correspond plutôt à la hauteur à partir de laquelle on passe d'un transport conductif à un transport turbulent.

Les fermetures basées sur une équation pronostique de l'énergie cinétique turbulente

La décomposition de Reynolds, qui a permis de faire apparaître des termes croisés du type $\overline{w'u'}$ dans les équations des variables moyennes (Section ), peut bien sûr être poussée plus loin. On écrit alors des équations d'évolution pour les quantités turbulentes , , $\theta '$ en soustrayant la moyenne d'ensemble à l'équation complète. On peut ainsi obtenir des équations d'évolution temporelle des termes croisés $\overline{w'u'}$ , $\overline{{w' }^2}$ , $\overline{{\theta'}^2}$ , ... On peut donc imaginer des fermetures où, au lieu de spécifier directement $\overline{w'u'}$ en fonction des grandeurs moyennes (de type diffusion turbulente), on considère le flux turbulent lui-même comme une variable indépendante du modèle suivant sa propre évolution. Mais ces nouvelles équations font elles-mêmes apparaître des termes du troisième ordre.

Une littérature très savante a été consacrée à ces fermetures à des ordres plus élevés. Cette histoire, qui semble avoir commencé dans les années 50, a produit les premiers modèles utilisables comme fermetures turbulentes au début des années 70. C'est le cas notamment du travail de Mellor et Yamada (1974). Dans cette approche, la fermeture (dite d'ordre 2) s'effectue au niveau des termes du troisième ordre, en introduisant notamment une mesure de l'anisotropie de la turbulence. De façon générale, on peut calculer les termes croisés à partir d'équations pronostiques qui font apparaître trois types de termes : des termes de transport, des corrélations pression-vitesse et des termes de dissipation. On aboutit typiquement alors à une dizaine d'équations pour prédire les différents termes croisés.

Les fermetures utilisées en pratique dans les modèles atmosphériques, comme celle de Mellor et Yamada (1974) ou les fermetures dites $K-\epsilon$ , sont des versions simplifiées des fermetures d'ordre 2 où on se focalise sur l'équation d'évolution de l'énergie cinétique turbulente

$\begin{displaymath} e=\frac{1}{2}\left[\overline{u'^2}+\overline{v'^2}+\overline{w'^2}\right] \end{displaymath}$

(2.11)

cette énergie étant ensuite utilisée pour le calcul du coefficient de mélange

En effet, si on revient à la présentation de la longueur de mélange faite précédemment, il est naturel de prendre $\sqrt{e}$ comme amplitude des mouvements turbulents verticaux $\overline{\vert w'\vert}$ dans le cas d'une turbulence isotrope en atmosphère neutre. Dans le cas d'une atmosphère stable, où l'on s'attend à des mouvements anisotropes plutôt horizontaux, $\sqrt{e}$ fournira plutôt une surestimation de $\overline{\vert w'\vert}$ (respectivement une sous-estimation pour une atmosphère instable).

Sous les hypothèses de couche limite mentionnées plus haut, on peut montrer (voir e. g. Stull, 1988) que l'évolution de l'énergie cinétique turbulente s'écrit sous la forme

$\begin{displaymath} \frac{\partial e}{\partial t}=P-D -\frac{1}{\rho}\frac{\part... ...artial z} -\frac{\partial \overline{w'e}}{\partial z}-\epsilon \end{displaymath}$

(2.12)

Le terme

$\begin{displaymath} P=-\overline{w'u'}\frac{\partial u}{\partial z}-\overline{w'v'}\frac{\partial v}{\partial z} \end{displaymath}$

(2.13)

est généralement positif (en particulier si $\overline{w'u'}=-K_z\partial u/\partial z$ ) et correspond à la production mécanique de turbulence. Le terme

$\begin{displaymath} D=-\frac{g}{\theta}\overline{w'\theta'} \end{displaymath}$

(2.14)

est en général positif dans une atmosphère stable. Il correspond alors à l'inhibition (ou destruction) de la turbulence par les effets de stratification. Il peut devenir un terme de production dans les atmosphères instables. Les termes suivants sont un terme de pression, le transport vertical turbulent d'énergie cinétique turbulente (le transport par la grande échelle est négligé car

est supposé petit), et $\epsilon$ est le terme de dissipation mécanique de la turbulence. Physiquement, l'énergie contenue principalement dans les plus grosses structures turbulentes ``cascade" vers les échelles moléculaires où elle est dissipée.

Mellor et Yamada, à partir d'une paramétrisation des termes du troisième ordre (à la base il s'agit donc d'une fermeture à l'ordre 2) proposent une série de simplifications, avec une hiérarchie de schémas. Nous présentons ici le schéma de niveau 2.5 (qui malgré cette nomenclature originale n'est qu'une approximation relativement grossière d'un schéma d'ordre 2), le plus largement testé dans des applications météorologiques, et dont certains résultats sont présentés plus loin. Dans le modèle 2.5, les termes croisés s'écrivent formellement, comme pour la diffusion turbulente,

$\begin{displaymath} \overline{w'\phi'}=-K_\phi \frac{\partial \phi}{\partial z} \end{displaymath}$

(2.15)

On écrit la diffusivité turbulente $K_\phi=lqS_\phi$ en fonction d'une longueur de mélange,

, d'une vitesse typique $q=\sqrt{2e}$ et d'une fonction de stabilité, qui peut prendre des valeurs différentes pour les traceurs ( $S_\phi$ ), la température (

) ou la quantité de mouvement (

Un des tours de force de ce travail est de faire sortir directement de simplifications successives d'une fermeture du second ordre une formulation analytique pour les fonctions et . Ces grandeurs sont fonctions du seul nombre de Richardson de flux, , mesure de l'importance relative du forçage mécanique de la turbulence et de l'inhibition par stratification. On retient ici les valeurs numériques données par Yamada (1983). Cette version a été utilisée dans une étude d'intercomparaison de fermetures turbulentes par Ayotte et al. (1996). Elle est également utilisée plus loin dans des tests numériques. Pour la quantité de mouvement, on obtient :

$\begin{displaymath} S_m=1,96\frac{\left(0,1912-Ri_f\right)\left(0,2341-Ri_f\right)}{\left(1-Ri_f\right)\left(0,2231-Ri_f\right)}, Ri_f< 0,16 \end{displaymath}$

(2.16)

$\begin{displaymath} S_m=0,085, Ri_f\ge 0,16 \end{displaymath}$

(2.17)

Pour la température, on obtient pour l'inverse du nombre de Prandtl turbulent $\omega=K_h/K_m=S_h/S_m$ , l'expression

$\begin{displaymath} \omega=1,1318\frac{0,2231-Ri_f}{0,2341-Ri_f}, Ri_f< 0,16 \end{displaymath}$

(2.18)

$\begin{displaymath} \omega=1,12, Ri_f\ge 0,16 \end{displaymath}$

(2.19)

En remarquant que $Ri_f=\omega Ri$ , où

$\begin{displaymath} Ri=\frac{g}{\theta}\frac{\frac{\partial \theta}{\partial z}}... ...t\frac{\partial {\bf v}}{\partial z}\right\vert\right\vert}^2} \end{displaymath}$

(2.20)

est le nombre de Richardson gradient, on obtient finalement

$\begin{displaymath} Ri_f= 0,6588[\ \ Ri+0,1776- \sqrt{Ri^2-0,3221Ri+0,03156}\ \ ]\mbox{, } Ri<Ri_c \end{displaymath}$

(2.21)

$\begin{displaymath} Ri_f=Ri_{fc}\mbox{, } Ri\ge Ri_c \end{displaymath}$

(2.22)

où $Ri_{fc}=0,191$ et

sont des nombres de Richardson critiques. On obtient ainsi des expressions explicites de

et $\omega$ en fonction des seules variables de grande échelle.

Dans le modèle 2.5, le terme de pression est négligé et le transport vertical est traité comme une diffusion turbulente de sorte que l'évolution de l'énergie cinétique turbulente ou de s'écrit finalement

$\begin{displaymath} \frac{1}{2}\frac{\partial q^2}{\partial t}= q l S_m{\left\ve... ...al z} \left[l q S_q \frac{\partial q^2 / 2}{\partial z}\right] \end{displaymath}$

(2.23)

avec

. Yamada (1983) utilise

à la fois pour l'eau et pour les traceurs passifs.

La paramétrisation de Mellor et Yamada, même dans ses versions plus sophistiquées, laisse en fait une grosse zone d'ombre sur la spécification de la longueur de mélange . Dans certains articles, les auteurs ont suggéré une équation pour analogue à l'équation pour mais cette équation est beaucoup moins fondée que sa grande s $\oe$ ur (selon les auteurs eux-mêmes). Le plus souvent, la longueur de mélange est spécifiée, par exemple en utilisant la formule de Blackadar (1962)

$\begin{displaymath} l=l_0\frac{\kappa z}{\kappa z+l_0} \end{displaymath}$

(2.24)

asymptotique à $\kappa z$ près de la surface et à une longueur

au-dessus de la couche de surface. La longueur

peut elle-même être fixée à une constante (de l'ordre de 100-200 m la plupart du temps) ou être calculée en fonction d'autres quantités. Yamada (1983) utilise par exemple pour

une altitude moyenne pondérée par l'intensitée de la turbulence

$\begin{displaymath} l_0=0,2\ \frac{\int_0^\infty z q dz}{\int_0^\infty q dz} \end{displaymath}$

(2.25)

C'est cette formulation qui est retenue pour la version introduite dans le modèle LMDZ et dont on montre certains résultats plus loin. ^2.1

Figure 2.1: Cycle diurne de la couche limite sur Mars. Comparaison d'observations par les sondes Viking (qui ont fonctionné plusieurs années à la surface de Mars dans les années 70) et les résultats de simulations numériques effectuées avec une version unidimensionnelle du modèle de circulation générale martien du LMD. a : Cycle diurne de la température de l'air mesurée au bout du bras Viking (à 1,6 m au-dessus du sol) et simulée dans la première couche du modèle (altitude de 4 m). b : Module du vent horizontal (m s $^{-1}$ ). c : Fluctuations turbulentes du vent (m s $^{-1}$ ). Pour b) et c), on montre, à gauche, les résultats du modèle en fonction de l'heure locale et de l'altitude et, à droite, les évolutions comparées des mêmes grandeurs près de la surface pour le modèle et les observations Viking. Les axes horizontaux correspondent aux heures locales. Pour les fluctuations turbulentes (c), les données sont calculées à partir de mesures haute fréquence du vent et les valeurs simulées sont estimées à partir de l'énergie cinétique turbulente prédite par la paramétrisation de Mellor et Yamada.

a) cycle diurne de la température
$\includegraphics[width=7cm]{therm/FIGURES/MARS/pbl_t.eps}$
b) cycle diurne du vent
$\includegraphics[width=7cm,clip]{therm/FIGURES/MARS/pbl_u_ref.eps}$	$\includegraphics[width=7cm]{therm/FIGURES/MARS/pbl_v.eps}$
c) cycle diurne de l'énergie cinétique turbulente
$\includegraphics[width=7cm,clip]{therm/FIGURES/MARS/pbl_q_ref.eps}$	$\includegraphics[width=7cm]{therm/FIGURES/MARS/pbl_q.eps}$

Avant de la tester couplée au modèle du thermique dans la version terrestre de LMDZ (comme on l'explique plus loin), nous avons avec Richard Fournier introduit la fermeture 2.5 de de Mellor et Yamada dans la version martienne du modèle (Forget et al., 1999). Avec sa fine atmosphère de CO, l'immense désert martien connaît des cycles diurnes très marqués. Dans les tropiques, ou l'été dans les moyennes latitudes, la température de surface peut varier de plusieurs dizaines de degrés entre la nuit et le jour. On montre sur la Fig. 2.1 un exemple de comparaison de résultats de simulations numériques avec des observations par les sondes Viking. On voit sur cette figure un cas typique de jet nocturne comme il en existe dans les déserts terrestres. Le mélange vertical intense pendant la journée s'éteint subitement en fin d'après-midi (voir la brusque chute de sur les figures c). Le vent qui était dans l'après-midi en équilibre entre un gradient de pression et ce terme de mélange se trouve subitement en déséquilibre et entre en oscillation inertielle. Dans les premières couches du modèle, le vent diminue rapidement sous l'effet du frottement turbulent sur la surface tandis que les couches supérieures, découplées de la surface, accélèrent. On obtient alors en sommet de couche limite une couche fortement cisaillée qui génère à nouveau de la turbulence. Les comparaisons aux données sont bonnes en générale pour les quantités moyennes ainsi que pour l'énergie cinétique turbulente la nuit. On note cependant une très forte sous-estimation des fluctuations turbulentes du vent le jour, dues sans doute à la non prise en compte des mouvement convectifs (comme on l'explique plus loin).

Les fermetures basées sur un équilibre de l'énergie cinétique turbulente

On peut aller un cran plus loin dans les simplifications en supposant que les termes de production ou destruction d'énergie turbulente et la dissipation sont constamment à l'équilibre $P-D=\epsilon$ , d'où l'on tire

$\begin{displaymath} q^2=l^2{\left\vert\left\vert\frac{\partial {\bf v}}{\partial z}\right\vert\right\vert}^2S_m B_1\left(1-\omega Ri\right) \end{displaymath}$

(2.31)

et donc

$\begin{displaymath} K_m=l^2\left\vert\left\vert\frac{\partial {\bf v}}{\partial ... ...ight\vert\right\vert\sqrt{{S_m}^3 B_1\left(1-\omega Ri\right)} \end{displaymath}$

(2.32)

On voit qu'on retombe exactement sur l'approche de Prandtl avec un terme correctif ( $\sqrt{{S_m}^3 B_1\left(1-\omega Ri\right)}$ ) rendant compte des effets de la stratification.

La formulation utilisée dans le modèle original du LMD (Laval et al., 1981) est également basée sur un modèle stationnaire de l'énergie cinétique turbulente. Seuls les coefficients diffèrent. Le coefficient de mélange s'écrit simplement sous la forme

$\begin{displaymath} K_m=l\ Max\left[ l \left\vert\left\vert\frac{\partial {\bf v... ...\right\vert\right\vert\sqrt{1- Ri/Ri_c},\sqrt{e_{min}} \right] \end{displaymath}$

(2.33)

où

est un nombre de Richardson critique. De façon un brin arbitraire, la longueur de mélange utilisée dans la version standard du modèle décroît quand on s'éloigne de la surface comme

avec

35 m.

Dans le monde des sciences de l'ingénieur, les fermetures dites $K-\epsilon$ sont davantage utilisées que le modèle de Mellor et Yamada. Ces fermetures font intervenir deux équations pronostiques, l'une pour l'énergie turbulente et l'autre pour la dissipation $\epsilon$ .

Les paramétrisations basées sur des relations de similitude

Les méthodes de similitude ont remporté un grand succès dans l'explication des observations des grandeurs turbulentes dans la couche limite de surface. Dans cette approche, on s'intéresse à une couche limite en régime stationnaire, on adimensionalise les équations et on dérive des relations ou modèles à partir des seuls paramètres dont dépendent les équations. La couche logarithmique, présentée plus haut à partir de la longueur de mélange, peut déjà être présentée à partir des relations de similitude si on remarque que le gradient vertical du vent près de la surface dans une atmosphère neutre ne peut dépendre que de et .

Monin et Obukov ont introduit les effets de la stratification dans cette description de la turbulence mécanique. L'hypothèse de base de leur théorie est de supposer que le cisaillement adimensionnel, $\kappa u_*\partial u/\partial z$ , égal à 1 pour une atmosphère neutre, ne dépend que d'une mesure de l'importance relative des flux de moment et de chaleur, le nombre de Richardson de flux introduit plus haut. Dans la couche de surface, ce nombre s'écrit

$\begin{displaymath} {Ri}_f=\frac{g}{\theta}\frac{\overline{w'\theta'}}{\overline... ... z {\overline{w'\theta'}}_0}{\kappa {u_*}^3\theta}=\frac{z}{L} \end{displaymath}$

(2.34)

où

$\begin{displaymath} L=\frac{\kappa {u_*}^3 \theta}{g{\overline{w'\theta'}}_0} \end{displaymath}$

(2.35)

est la longueur d'Obukov.

On suppose donc, dans la couche de surface, que le gradient de vent peut s'écrire sous la forme

$\begin{displaymath} \frac{\partial u}{\partial z}=\frac{u_*}{\kappa z}\phi_m\left(\frac{z}{L}\right) \end{displaymath}$

(2.36)

où $\phi_m$ est une fonction universelle, dite fonction de stabilité, qui vaut 1 dans les conditions neutres ou quand $z \rightarrow 0$ , c'est à dire quand la turbulence est dominée par les effets mécaniques. Cette équation peut s'intégrer verticalement pour fournir directement le vent à un niveau donné dans la couche de surface :

$\begin{displaymath} u=\frac{u_*}{\kappa }\left[\ln\left(\frac{z}{z_0}\right)-\psi_m\left(\frac{z}{L}\right)\right] \end{displaymath}$

(2.37)

où

$\begin{displaymath} \psi_m\left(\xi\right)=\int_{z_0/L}^{\xi}\left[1-\phi_m\left... ...0 }^{\xi}\left[1-\phi_m\left(\xi\right)\right]\frac{d\xi}{\xi} \end{displaymath}$

(2.38)

La seconde forme de $\psi_m$ permet en général des calculs analytiques relativement simples et suffisamment précis.

Comme pour la quantité de mouvement, on suppose que le gradient de température adimensionnel est relié au gradient neutre par une fonction de

$\begin{displaymath} \frac{\partial \theta}{\partial z}=\frac{R \theta_*}{\kappa z}\phi_h\left(\frac{z}{L}\right) \end{displaymath}$

(2.39)

avec $\theta_*= {\overline{w'\theta'}}_0/u_*$ qui conduit comme pour le vent à

$\begin{displaymath} \theta-\theta_s=\frac{R \theta_*}{\kappa }\left[\ln\left(\frac{z}{z_0}\right)-\psi_h\left(\frac{z}{L}\right)\right] \end{displaymath}$

(2.40)

où $\theta_s$ est la température potentielle associée à la température du sol. L'introduction du nombre de Prandtl $R\simeq 0.7$ dans la formule suit la présentation par Deardorff (1972b) et permet de choisir la fonction $\phi_h$ égale à l'unité pour les conditions neutres ou dominées par la turbulence mécanique.

Des campagnes de mesures ont été dédiées à la mesure de ces fonctions. Les formules proposées par Businger et al. (1971) et ajustées sur les résultats d'une campagne dans le Kansas sont encore largement utilisées. Ces formules sont données dans la Table 2.1. De nombreuses versions modifiées ont été proposées depuis, utilisant parfois des données plus récentes (cf. par exemple Högström, 1988).

Tableau 2.1: Fonctions de Businger Dyer telles qu'elles sont utilisées dans les simulations présentées plus loin.

	instable	stable $L \ge 0$


$\phi_m\left(\xi\right)$	$\left(1-\gamma_1\xi\right)^{-\frac{1}{4}}$	$1 + \gamma_3\xi$

$\psi_m\left(\xi\right)$	$\ln\left[\left(\frac{1+x^2}{2}\right){\left(\frac{1+x}{2}\right)}^2\right]-2 \tan^{-1}x +\frac{\pi}{2}$	$-\gamma_3\xi$
	avec $x=1/\phi_m$


$\phi_h\left(\xi\right)$	$\left(1-\gamma_2\xi\right)^\frac{1}{2}$	$1+\gamma_3\xi/R$

$\psi_h\left(\xi\right)$	$2 \ln\left(\frac{1+x}{2}\right)$	$-\gamma_3\xi/R$
	avec $x=1/\phi_h$

Notons que Mellor (1973) a appliqué son modèle de fermeture à la couche limite de surface, en supposant des flux constants, et réussi à interpréter les mesures de Businger et al. (1971) avec des jeux de coefficients compatibles avec des mesures de souffleries en conditions neutres.

Certains auteurs ont essayé d'extrapoler les approches en similitude à la couche mélangée. Il faut alors introduire au moins un paramètre supplémentaire : la hauteur de la couche limite . Nous testons par exemple plus loin une formulation analytique du coefficient de diffusion turbulente en fonction de l'altitude proposée à l'origine par Brost et Wyngaard (1978)

$\begin{displaymath} K_m=u_*\kappa z \phi_m^{-1}\left[1-\frac{z}{h}\right]^p \end{displaymath}$

(2.41)

Cette formulation est choisie de façon a être asymptotique à la prédiction de Monin-Obukov dans la couche limite de surface et à s'annuler en

. Brost et Wyngaard (1978) retiennent un exposant

Pour appliquer cette formule, il faut commencer par estimer la hauteur de la couche limite . Une approche maintenant classique et relativement robuste consiste à utiliser un nombre de Richardson non local :

$\begin{displaymath} R_{ib}(z)=\frac{g\,(z-z_{1})}{\theta(z)}\frac{(\theta(z)-\theta(z_{1}))} {u(z)^{2}+v(z)^{2}} \end{displaymath}$

(2.42)

où

est une altitude de référence proche de la surface. On définit la hauteur de la couche limite comme l'altitude à laquelle ce nombre dépasse une valeur seuil, typiquement de l'ordre de 0,2-0,25. Cette approche est par exemple retenue dans les travaux de Troen et Mahrt (1986).

Spécificités de la couche limite convective

A la base, les formulations - plus ou moins sophistiquées - en diffusion turbulente font l'hypothèse que la longueur caractéristique des mouvements turbulents est petite devant les échelles spatiales typiques, et notamment devant la hauteur de la couche limite.

Les limites de cette approche sont reconnues depuis longtemps, en particulier dans le cas des couches limites convectives, où les ascendances thermiques, résultant de l'accumulation de chaleur près de la surface, s'organisent sous forme de panaches ou de rouleaux à des échelles comparables aux échelles de la couche limite. Dans la suite, on appellera meso-échelle cette échelle des structures convectives de couche limite. Dans la couche limite convective, le flux de chaleur, dirigé vers le haut pour évacuer l'énergie accumulée à la surface, est souvent associé à un profil neutre ou même marginalement stable de température potentielle, c'est à dire que le flux d'énergie remonte le gradient, du froid vers le chaud, ce qui est incompatible avec une approche en diffusion.

Les couches limites convectives se caractérisent plus précisément en trois régions :
- une couche de surface instable chauffée directement par le sol,
- une couche mélangée épaisse typiquement de 1 à 2 km dans les régions tempérées mais qui peut atteindre 3 km aux jours les plus chauds de l'été même en région parisienne et plus de 5 km sur les déserts ou sur la planète Mars.
- une couche d'inversion très stable, épaisse de quelques dizaines à quelques centaines de mètres. La hauteur de cette inversion est souvent utilisée comme hauteur de couche limite. ^2.2

Ce sont les particules d'air de la couche de surface, particules plus chaudes donc plus légères que celles de la couche mélangée, qui s'élèvent dans la couche mélangée pour s'organiser en ascendances thermiques sous forme de rouleaux, de cellules ou de panaches isolés. L'accélération d'une particule de la couche de surface dans l'environnement est donnée par

$\begin{displaymath} \gamma= g\frac{{\theta_v}_P - {\theta_v}_e}{{\theta_v}_e} \end{displaymath}$

(2.43)

où

$\begin{displaymath} \theta_v=\theta (1+0.061 q) \end{displaymath}$

(2.44)

est la température potentielle virtuelle et

est l'humidité spécifique. Cette température potentielle tient compte, pour le calcul de la flottabilité, des changements de masse molaire de l'air dus aux changements de contenu en vapeur d'eau.

Avant de présenter quelques approches pour paramétriser la couche limite dans ces conditions particulières, et de décrire en détail le ``modèle du thermique", on présente dans cette section une analyse d'échelle de la couche limite convective ainsi que les grandes lignes des connaissances sur le sujet, que ce soit au travers d'observations ou de simulations dites des grands tourbillons (ou Large Eddy Simulation en anglais).

Parce que les développements proposés ici concernent essentiellement la couche limite convective en ciel clair, on ne parlera pratiquement pas de nuages, même s'il est clair que la capacité de la nouvelle paramétrisation à prédire les caractéristiques statistiques des nuages (couverture nuageuse, contenu en eau des nuages) sera un élément essentiel de sa possible adoption comme paramétrisation de base d'un modèle de climat.

Organisation à meso-échelle

**Figure 2.2:** Rues de nuages observées dans la mer de Bering.
$\includegraphics[width=16cm]{therm/FIGURES/Bering.eps}$

**Figure 2.3:** Vue schématique de l'organisation de la convection de couche limite en rouleaux, le long de l'axe du vent. D'après Brown (1980).
$\includegraphics[width=14cm,clip]{therm/FIGURES/SCAN/w1.eps}$

L'existence de structures organisées dans la couche limite convective est bien connue des amateurs de vol libre qui utilisent les ``pompes" thermiques pour gagner de l'altitude. Les vitesses verticales rapportées par ces amateurs sont typiquement de 1 à 4 m s $^{-1}$ en plaine et plutôt de 5 à 10 m s $^{-1}$ en montagne.

Ces structures peuvent prendre la forme de panaches isolés ou s'organiser en forme de cellules ou de rouleaux. Un travail d'investigation systématique de ces structures, notamment à partir de vols avions, a été entrepris depuis une trentaine d'années, à partir des travaux pionniers de LeMone (1973). Les mesures in-situ à bord d'avions, les photos satellites à haute résolution, les instruments de détections active (radar et lidar) ainsi que les simulations numériques dites ``des grands tourbillons" (Large Eddy Simulations en anglais) ont permis de mieux comprendre et caractériser les structures organisées de la couche limite. On se contente ici de montrer quelques illustrations issues de ces études.

Les rues de nuages constituent une des réalisations les plus spectaculaire de l'organisation de la convection de couche limite. Les structures de rues s'observent à toutes les latitudes et en toutes saisons, mais les arrivées sur la mer d'air très froid ayant séjourné un moment sur des glaciers ou des banquises offrent souvent des photos spectaculaires comme celle montrée sur la Fig. 2.2. L'air froid et sec, en arrivant sur la mer plus chaude (ici le vent souffle du nord au milieu de la mer de Bering) donne naissance à une couche limite convective. Au début, l'air est encore clair. Il se charge petit à petit en humidité et des cumulus se mettent à bourgeonner en sommet de la partie ascendante de grands rouleaux convectifs, créant ces grandes rues de nuages alignées le long du vent dominant. L'image faisant un millier de kilomètres de large environ, on voit que les rues de nuages sont typiquement espacées de 5 km dans la partie nord et jusqu'à une vingtaine dans le sud.

Au sud de la zone, la structure en rouleaux disparaît au profit d'une organisation en cellules. Mais, dans un cas comme dans l'autre, on distingue nettement une organisation à une échelle de quelques kilomètres à quelques dizaines de kilomètres.

A noter également, dans le sud-ouest de la photo, au sud (en aval) des Aléoutiennes, des structures transversales associées très vraisemblablement à des ondes de gravité piégées dans le sillage des reliefs que constituent les îles.

Figure 2.4: Simulations des grands tourbillons de la couche limite convective d'après Moeng et Sullivan (1994). Coupes horizontales instantanées à

pour deux simulations (B à gauche et SB1 à droite). On montre, en haut, le vent vertical

(m s $^{-1}$ ) et, en bas, les perturbations de température potentielle virtuelle $\theta _v$ (K).

Simulation B	Simulation SB1
$\includegraphics[height=7cm]{therm/FIGURES/SCAN/moengBw02.eps}$	$\includegraphics[height=7cm]{therm/FIGURES/SCAN/moengSBw02.eps}$
$\includegraphics[height=7cm]{therm/FIGURES/SCAN/moengBT02.eps}$	$\includegraphics[height=7cm]{therm/FIGURES/SCAN/moengSBT02.eps}$

Caractéristiques et valeurs des iso-contours pour les deux simulations.
B : domaine de 3 $\times$ 3 km

, $\overline{w'\theta_0'}$ =0.24 m s $^{-1}$ K,

=10 m s $^{-1}$

: (-2;-1,5;-1;-0,5;0,5;1;1,5;2;2,5;3), gris [sombre/clair] pour [

$\theta_v'$ : (-0,3;-0,2;-0,1;0,1;0,2;0,3;0,4;0,5), [ $\theta_v'>0,1/\theta_v'<-0,1$ ]
SB1 : domaine de 5 $\times$ 5 km

, $\overline{w'\theta_0'}$ =0.05 m s $^{-1}$ K,

=15 m s $^{-1}$

: (-1,8;-1,5;-1,2;-0,9;-0,6;-0,3;-0,1;0,1;0,3;0,6;0,9;1,2;1,8), [

]
$\theta_v'$ : (-0,2;-0,1;-0,05;0,05;0,1;0,2), [ $\theta_v'>0,1/\theta_v'<-0,1$ ]

**Figure 2.5:** Perturbation du vent zonal à deux altitudes dans la simulation SB1 de Moeng et Sullivan (1994).
$\includegraphics[width=14cm]{therm/FIGURES/SCAN/moengu.eps}$ Valeur des iso-contours : (-3;-2,5;-2;-1,5;-1;-0,5;-0,1;0,1;0,5;1;1,5;3), gris [sombre/clair] pour [ ]

De nombreux travaux théoriques et numériques ont été consacrés à l'étude de ces structures convectives de la couche limite.

Une étude théorique des instabilité de la couche d'Eckman a permis de prédire l'orientation des rouleaux par rapport aux vents dominants (Brown, 1972). Les rouleaux sont alignés à 30 $^{o}$ à gauche du vent pour des couches stables, 18 $^{o}$ pour des couches limites neutres et essentiellement alignées avec le vent pour des couches instables. Le modèle de la couche d'Eckman n'est pas vraiment applicable à la couche limite instable mais la prédiction est cependant relativement proche de l'observation (LeMone, 1973) même si les rouleaux sont plutôt orientés également à 10-20 $^{o}$ du vent dans la couche convective. On montre sur la Fig. 2.3 une vision schématique de cette organisation en rouleaux.

Après le travail pionnier de Sommeria et LeMone (1978), de nombreuses simulations des grands tourbillons ont été consacrées à l'organisation de la couche limite convective. Dans les simulations des grands tourbillons, on résout explicitement des équations dynamiques non-hydrostatiques (différentes approximations sont cependant utilisées pour filtrer les modes acoustiques les plus rapides) jusqu'à une échelle typique de 20 à 100 m suivant les cas. On suppose à cette échelle que la turbulence est bien représentée par des idées de cascades vers les petites échelles et par des fermetures locales de type $K-\epsilon$ ou Mellor et Yamada. A partir de telles simulations, Moeng et Sullivan (1994) et d'autres ont par exemple montré que la sélection entre les différents modes d'organisation était en grande partie contrôlée par l'importance relative des forçages thermiques (par le chauffage en surface) et mécanique (par le cisaillement de vent) de la turbulence.

Les simulations de Moeng et Sullivan (1994) sont relativement académiques, avec une turbulence de couche limite forcée par un flux de chaleur imposé en surface et un forçage géostrophique engendrant des cisaillements de vent et donc de la turbulence mécanique près de la surface. Il s'agit de couches limites non nuageuses. Le calcul est effectué sur un domaine carré de 3 ou 5 km avec une maille d'une cinquantaine de mètres dans les trois directions d'espace et des conditions aux limites périodiques horizontalement. En faisant varier indépendamment l'intensité du forçage thermique $\overline{w'\theta'}_0$ en surface et l'intensité du vent géostrophique , on trouve dans les simulations deux modes d'organisation : une organisation en rouleaux quand le cisaillement est important et des panaches isolés sans organisation apparente quand le flux de chaleur domine. Sur la Fig. 2.4, on montre des coupes instantanées à l'altitude , où est la hauteur de l'inversion, des perturbations de la vitesse verticale (en haut) et de la température potentielle virtuelle (en bas) pour deux simulations. Dans la première - appelée B pour ``buoyant" par les auteurs - le flux en surface vaut $\overline{w'\theta'}_0=0,24$ m K s $^{-1}$ avec un vent géostrophique de 10 m s $^{-1}$ . Dans la seconde simulation, le forçage mécanique est plus important avec $\overline{w'\theta'}_0=0,05$ m K s $^{-1}$ et =15 m s $^{-1}$ . Cette simulation avec cisaillement et flottabilité (shear and buoyancy) est appelée SB1 par les auteurs.

Pour les deux simulations, on voit clairement les structures thermiques, avec de l'air chaud associé à des vitesses ascendantes. Ces structures thermiques couvrent dans les deux cas une fraction relativement faible de la surface. La simulation B ne présente pas de structure bien marquée et on a plutôt l'impression de voir des panaches isolés. La simulation SB1, avec un forçage mécanique important, présente une organisation en rouleaux. La relative faible étendue du domaine fait qu'on ne simule que deux rouleaux. Des simulations plus récentes, utilisant des domaines plus grands par rapport aux structures représentées, confirment ces résultats (cf. par exemple Weckwerth et al., 1997).

Les subsidences amènent vers la surface de l'air provenant du haut du domaine et associé de ce fait à un excès de quantité de mouvement (le vent géostrophique est positif dans la direction ) comme on le voit sur la Fig. 2.5 pour la simulation SB1. En , sur la gauche de la figure, les structures organisées sont encore perturbées par les organisations à plus petite échelle dans la couche limite de surface. Un peu plus haut dans l'atmosphère et à droite sur la même figure, en , la structure en rouleaux domine encore davantage l'écoulement.

Dans ces simulations, le rapport d'aspect - rapport entre la séparation des rouleaux et la hauteur de la couche limite - est compris entre 2 et 3.

**Figure 2.6:** En haut, champs de réflectivité radar pour 2 situations particulières observées en Floride (a) le 6 août 1991 à 1700 UTC et (b) le 12 août 1991 à 2000 UTC. La température (nombre du haut) et le point de rosée (nombre du bas) ainsi que le vent sont superposés pour certaines stations d'observation. La figure de gauche montre une couche limite régulièrement organisée en rouleaux. La figure de droite montre un front de brise de mer sur la droite et des cellules en haut. Les figures du dessous correspondent à un lissage des échos radar pour le sous-domaine repéré par un carré dans les figures du haut, avec des contours tous les 4 dBZ à partir de 0. Les valeurs plus grandes que 4 et 8 dBZ sont grisées respectivement en gris clair et gris foncé. D'après Weckwerth et al. (1997).
$\includegraphics[width=14cm]{therm/FIGURES/SCAN/w2gb.eps}$ $\includegraphics[width=14cm]{therm/FIGURES/SCAN/w3b.eps}$

L'importance de l'organisation en cellules ou en rouleaux de la couche limite convective, même en l'absence de nuages, a été confirmée avec l'utilisation de plus en plus systématique de la télédétection active, lidar ou radar, pour observer l'atmosphère. Les échos lidar ou radar sont en effet souvent capables de distinguer, dans la couche limite, l'air montant depuis la couche de surface de son environnement. Pour les lidars, c'est la présence d'aérosols dans les panaches qui permet en général de les visualiser alors que, pour les radars, on pense qu'on voit souvent des insectes.

On montre sur la Fig. 2.6 deux exemples d'observations radar issus d'une étude de Weckwerth et al. (1997) montrant à gauche une organisation en rouleaux et, à droite, une organisation en cellules. Cette étude assez systématique d'observation de la couche limite convective en Floride a permis de confirmer certains résultats obtenus avec les simulations des grands tourbillons, comme l'apparition systématique de rouleaux dans certaines gammes de cisaillement et de flux de chaleur ou l'estimation du rapport d'aspect des rouleaux.

Le cycle diurne de la couche limite continentale

**Figure 2.7:** Cycle diurne de la couche limite continentale.
$\includegraphics[width=16cm,clip]{therm/FIGURES/sirtagb.eps}$ Profils verticaux de température potentielle et d'humidité relative enregistrés par les radiosondages de Trappes à midi pour trois jours successifs, les 26, 27 et 28 mai 2003, et les échos enregistrés au cours de la journée par le lidar aérosol LNA du SIRTA.

**Figure 2.8:** Observations de la couche limite en Oklaoma le 14 juin 2002 pendant la campagne IHOP. Transects avions à différentes heures de la journées avec des mesures radar en dessous et au-dessus de l'avion. On voit bien la croissance des thermiques qui atteignent, en milieu de journée, environ 1,5 km.
$\includegraphics[width=15cm]{therm/FIGURES/FLEUR/Bart_fig.eps}$

Des couches limites convectives particulièrement développées sont observées l'après-midi sur les continents, notamment par beau temps. Sur les déserts, elles peuvent atteindre plus de 5 km d'altitude. On montre sur la Fig. 2.7 un exemple d'observations de couches limites convectives sur trois jours consécutifs au SIRTA, le site instrumenté atmosphérique de l'IPSL. Les figures de droite montrent l'écho lidar (le lidar LNA^2.3) observé au cours du temps à la verticale de l'école Polytechnique, à Palaiseau. Sur ces figures, la couche limite est matérialisée grossièrement comme la zone gris clair, correspondant à une réflexion sur des aérosols. La couche limite nocturne très fine (quelques centaines de mètres) se développe dans la matinée, entre 9 heures et midi. Dans cette phase de croissance, on voit très clairement des panaches ascendants, plus clairs que l'air environnant. En début d'après-midi, la couche limite convective est bien développée et on voit se former, au sommet des thermiques, des cumulus qui réfléchissent totalement le signal.

Les profils de température potentielle associés, observés à Trappes par radiosondage à midi, montrent tous les trois une légère instabilité dans la couche de surface et un profil très bien mélangé sur 1 à 1,5 km suivant les jours. Dans cette région, l'humidité spécifique est relativement bien mélangée également. Ceci correspond à une humidité relative qui croît avec l'altitude, pour approcher les 100 $\%$ en sommet de couche limite, là où les nuages sont observés.

Cette pulsation de la couche limite entre couche limite nocturne stable et couche limite convective développée dans l'après-midi conditionne au premier chef les concentrations observées pour les espèces émises en surface.

On montre également sur la Fig. 2.8 des mesures radar aéroportées obtenues pendant la campagne IHOP qui s'est déroulée pendant l'été 2002 dans l'Oklaoma. Là encore, on voit se développer les thermiques en cours de matinée. On voit également que les petits thermiques du matin s'organisent peu à peu en ascendances plus importantes et plus espacées.

Caractérisation des grandeurs turbulentes dans la couche limite convective

**Figure 2.9:** Principe d'une analyse en composites d'évènements chauds de mesures avions de la turbulence de couche limite. La partie du haut montre une séquence de mesures de la température potentielle et du vent vertical. Une courbe lissée des températures potentielles est utilisée pour identifier les évènements chauds. Chaque évènement est associé à un segment qui est ensuite étiré dans l'espace pour ramener tous les segments à une longueur identique. Pour chaque variable, on peut alors construire des moyennes ou des écart-types pour un thermique moyen. La figure et l'approche sont issues d'une très jolie étude de de Williams et Hacker (1992).
$\includegraphics[width=15cm]{therm/FIGURES/SCAN/williams1crop.eps}$

**Figure 2.10:** Analyse en composites d'évènements chauds d'une série de vols avion effectués dans la couche limite convective en Australie. Les rangées correspondent à différentes gammes de valeurs de . Les deux rangées du haut correspondent au milieu de la couche mélangée. De gauche à droite sont représentés, la température potentielle $\theta$ , le flux de chaleur $H=\rho C_p \overline {w'\theta '}$ , l'humidité spécifique et la vitesse verticale , normalisées par les échelles $\theta ^*$ , $H_0=\rho C_p\overline {w'\theta '}_0$ , et . Les courbes pleines correspondent aux moyennes des grandeurs et les courbes pointillées aux écart-types associés. Le nombre d'évènements associés à chaque mesure est donné à gauche (22ev par exemple veut dire qu'on a fait des statistiques avec 22 segments). D'après Williams et Hacker (1992).
$\includegraphics[width=15cm]{therm/FIGURES/SCAN/williams2b.eps}$

Depuis les travaux de LeMone (1973), de nombreux travaux ont porté sur la caractérisation des ascendances thermiques à partir des mesures avions.

Pour quantifier les fluctuations turbulentes, on effectue des vols en avion aussi stables que possible en altitude avec un échantillonnage rapide, et on analyse les fluctuations de vent, température et humidité. Un exemple de séquence de mesure de la température potentielle et de la vitesse verticale est donné en haut de la Fig. 2.9. Pour cet exemple, issu du travail de Williams et Hacker (1992), la fréquence d'acquisition était de 13 Hz ce qui correspond, pour un avion qui volait en moyenne à 40 m s $^{-1}$ , à un pas d'échantillonnage d'environ 3 m.

On voit clairement apparaître sur ce cas particulier des évènements chauds, d'une longueur d'une centaine de mètres, associés à une vitesse verticale plutôt positive mais très bruitée.

La façon la plus classique d'analyser de telles observations consiste à calculer les flux par corrélation entre fluctuations de vent et de température ( $\overline {w'\theta '}$ ). Cette approche permet effectivement d'estimer les flux mais en perdant toute l'information sur les structures organisées. Il est en plus délicat de restituer la géométrie des structures méso-échelles traversées à partir des vols avions. En effet, les variations verticales ne peuvent être reconstituées qu'au travers de vols horizontaux successifs et qui n'explorent donc pas les mêmes panaches thermiques.

Williams et Hacker (1992) ont proposé une approche très éclairante sur la nature du transport dans la couche limite convective à partir de la construction d'un thermique moyen, défini comme un composite des évènements chauds. Les thermiques sont en général beaucoup plus facilement identifiables sur les mesures de $\theta$ , qui montrent une grande asymétrie entre un fond un peu froid et des évènements chauds intenses et relativement bien isolés, que sur celles du vent. Pour caractériser ces thermiques, Williams et Hacker (1992) commencent donc par identifier les segments chauds sur les mesures de $\theta$ - après un lissage - comme les portions où $\theta '$ dépasse 1 $\sigma$ où $\sigma$ est l'écart-type des fluctuations. Tous les segments sont ensuite ramenés par homothétie sur un segment de longueur unité. On peut alors, à partir de tout ces segments, calculer des moyennes ou écart-types de toutes les grandeurs mesurées pour construire une image d'un thermique moyen. Cette méthode est illustrée sur la Fig. 2.9.

Sur la Fig. 2.10, on montre les résultats obtenus par Williams et Hacker (1992) à différents niveaux dans la couche limite convective. Les différentes variables sont normalisées par des échelles caractéristiques. Pour le vent vertical, l'échelle utilisée est l'échelle convective (proposée à l'origine par Deardorff, 1970) construite à partir du flux de chaleur au sol et de la hauteur de la couche limite :

$\begin{displaymath} w_*=\left[\frac{g}{\theta}z_i\overline{w'\theta'}_0\right]^{\frac{1}{3}} \end{displaymath}$

(2.45)

On revient un peu plus loin sur le sens de cette échelle. L'échelle de température est

$\begin{displaymath} \theta_*=\frac{\overline{w'\theta'}_0}{w_*}= \left[\frac{g}{... ...\right]^{-\frac{1}{3}} {\overline{w'\theta'}_0 }^{\frac{2}{3}} \end{displaymath}$

(2.46)

On retrouve sur ces composites des caractéristiques déjà suggérées par l'observation directe des séquences avions. Les thermiques, sélectionnés à partir d'un excès de température, sont aussi associés à un vent vertical positif en moyenne. La valeur du vent ascendant au milieu du thermique est très proche de l'échelle de vitesse alors que l'écart type (en pointillé) est de l'ordre de grandeur du vent moyen dans et à proximité du thermique pour la vitesse ascendante. Ceci suggère que le thermique et ses abords immédiats sont le siège d'une turbulence de petite échelle importante, le thermique n'expliquant à lui seul que la moitié environ de $\overline{w'^2}$ . En regard, l'écart-type est deux à trois fois plus faible que la moyenne pour la température.

Le flux de chaleur dans l'ascendance est 4 fois supérieur au chauffage par la surface en bas et encore 2 à 3 fois supérieur au milieu de la couche mélangée alors qu'il est beaucoup plus faible autour. Les thermiques qui occupent typiquement 20 $\%$ de la surface dans les observations, semblent donc, à eux seuls, capables d'expliquer l'essentiel du transport de chaleur.

Dans le cas considéré, la vapeur d'eau est légèrement plus abondante dans le panache, mais avec une grande dispersion.

A noter qu'en regard de cette analyse en composites particulièrement éclairante, il existe une littérature relativement abondante dans laquelle les thermiques sont caractérisés sur la base de seuils sur les vitesses verticales. Ces études suggèrent qu'une fraction seulement (typiquement une bonne moitié) du flux de chaleur est contenu dans les plus grandes structures. On citera en particulier pour ces questions le travail de Schumann et Moeng (1991) qui compare des tris en vitesse verticale à la fois dans des simulations des grands tourbillons et dans des observations et le travail de Wang (1996) qui étudient l'influence du choix du seuil sur la caractérisation des structures organisées. Le tri sur les vitesses verticales est en fait particulièrement peu sélectif des thermiques car les variations turbulentes du vent sont dues aussi bien aux structures méso-échelles qu'à la turbulence de petite échelle, active dans toute la couche limite. Les fluctuations de petite échelle, en mélangeant un air déjà bien mélangé, n'affectent que peu les fluctuations turbulentes de la température. Les excès de température sont donc davantage caractéristiques de l'origine de l'air.

**Figure 2.11:** Distributions et distributions croisées des vitesses verticales et des fluctuations de la température potentielle virtuelle dans des simulations des grands tourbillons effectuées avec le modèle mésoNH pour un cas de couche limite convective observé dans les grandes plaines américaines pendant la campagne IHOP le 14 juin 2002. Des observations radar de ce cas particulier sont présentées sur la Fig. 2.8.
$\includegraphics[width=12cm]{therm/FIGURES/FLEUR/fleur.eps}$

On voit donc que les structures thermiques sont associées à des distributions de fluctuations turbulentes de ou $\theta$ fortement asymétrique. Ces distributions peuvent être calculées dans les simulations des grands tourbillons. On montre pour illustration sur la Fig. 2.11, les distributions de et $\theta_v'$ obtenues avec le modèle Meso-NH pour une simulation d'un jour particulier de la campagne IHOP (Couvreux, 2005), correspondant aux mesures radar aéroportées de la Fig. 2.8. La distribution croisée de $\theta '$ et montre un maximum important pour des températures basses et des vitesses verticales légèrement négatives et des températures plus chaudes associées à des vitesses verticales positives et plus importantes, les thermiques. Les distributions individuelles des deux variables montrent une forte asymétrie avec une queue de distribution du côté des vitesses positives et des températures chaudes. Pour , la forme de la distribution est en très bon accord avec celle déduite de vols avions (non montrées).

Analyse d'échelle de la couche limite convective

**Figure 2.12:** Représentation schématique de la couche limite convective
$\includegraphics[width=15cm]{therm/cbl.eps}$

Dans la couche mélangée, le transport de chaleur est donc effectué principalement par les structures de grande échelle avec des ascendances, associées à de l'air plus chaud provenant de la couche limite de surface, compensées par des subsidences plus froides. Si on idéalise l'ascendance thermique en supposant qu'elle est associée à un excès constant, ${\theta'}_a$ , de température potentielle virtuelle^2.4et une vitesse verticale (=) et si on suppose que les particules d'air dans le thermique montent sous l'effet de leur flottabilité, on a, avec l'approximation classique $p'/p\ll \theta'/\theta$ ,

$\begin{displaymath} \frac{dw_a}{dt}=g\frac{{\theta'}_a}{\overline{\theta}} \end{displaymath}$

(2.47)

On obtient, pour une ascendance stationnaire,

$\begin{displaymath} w_a\frac{\partial w_a}{\partial z}=g\frac{{\theta'}_a}{\overline{\theta}} \end{displaymath}$

(2.48)

qui s'intègre, en supposant que ${\theta'}_a$ ne dépend pas de

, en

$\begin{displaymath} w_a^2=2g\frac{{\theta'}_a}{\overline{\theta}}z \end{displaymath}$

(2.49)

Si on suppose que les ascendances couvrent une fraction $\alpha$ de la surface et qu'elles sont compensées par une subsidence de vitesse moyenne , associée à un déficit moyen de température potentielle ${\theta'}_d<0$ avec :

$\begin{displaymath} \alpha w_a+ (1-\alpha) w_d=0 \end{displaymath}$

(2.50)

$\begin{displaymath} \alpha {\theta'}_a+ (1-\alpha) {\theta'}_d=0 \end{displaymath}$

(2.51)

on voit que le flux de chaleur total s'écrit

$\begin{displaymath} \overline{w'\theta'}=\alpha w_a{\theta'}_a+ (1-\alpha) w_d{\theta'}_d =\frac{\alpha}{1-\alpha} w_a{\theta'}_a \end{displaymath}$

(2.52)

On obtient donc (en utilisant l'Eq. 2.49) une relation entre le flux de chaleur et la vitesse dans les ascendances

$\begin{displaymath} \overline{w'\theta'}=\frac{\alpha}{1-\alpha}\frac{\overline{\theta}{w_a}^3}{2gz} \end{displaymath}$

(2.53)

La couche limite convective est par nature très brassée et elle se caractérise donc par un profil de température potentielle très homogène. De ce fait un excès de température en surface se répartit très rapidement dans l'ensemble de la couche mélangée. Le flux de chaleur fourni par la surface produit donc un réchauffement quasi uniforme de la couche limite, jusqu'à la couche d'inversion. En conséquence aussi, la divergence du flux de chaleur doit être nulle dans la couche mélangée, c'est à dire que le flux doit décroître linéairement depuis la surface.

Cette idée est illustrée sur la Fig. 2.12. La courbe noire à gauche montre un profil typique de température potentielle dans la couche limite convective, avec une couche de surface instable, une couche mélangée neutre et une inversion en sommet de couche limite. Le flux de chaleur en surface brasse la couche mélangée provoquant un chauffage homogène de la couche mélangée et un léger refroidissement au niveau de l'inversion. Le schéma du milieu montre le flux moyen de chaleur associé à cette évolution de la température potentielle. La dérivée du flux de chaleur s'annule à l'altitude où la température n'évolue pas.

Pour l'analyse d'échelle, on voit donc que le flux de chaleur décroît linéairement pour s'annuler à une altitude légèrement inférieure à . Si on se place au milieu de la couche limite en , le flux calculé par l'Eq. 2.53 doit être proche de $\overline{w'\theta'_0}/2$ d'où

$\begin{displaymath} \frac{\overline{w'\theta'_0}}{2}\simeq \frac{\alpha}{1-\alpha}\frac{\overline{\theta}{w_a}^3}{gz_i} \end{displaymath}$

(2.54)

d'où l'on déduit

$\begin{displaymath} w_a\simeq w_*{\left[\frac{1-\alpha}{2\alpha}\right]}^{\frac{1}{3}} \end{displaymath}$

(2.55)

où

est l'échelle de vitesse convective définie plus haut (Eq. 2.45). Cette formule est relativement peu sensible à la fraction $\alpha$ avec $w_a\simeq w_*$ pour $\alpha=30\%$ ou $w_a\simeq 2w_*$ pour $\alpha=5\%$ .

Les fermetures non locales et la couche limite convective

Avant de présenter le modèle du thermique, on donne ci-dessous un aperçu des différentes approches qui ont été proposées ou utilisées pour pallier les déficiences de la diffusion turbulente dans les cas de couches limites convectives.

Contre-gradient et modèles non locaux

Une première approche simple pour pallier le problème du transport de chaleur en remontant le gradient dans la couche limite, consiste à calculer la diffusion, non pas par rapport à un profil neutre de température potentielle mais par rapport à un profil légèrement stable.

Cette approche a été proposée par Deardorff (1966) et consiste à prescrire le flux de chaleur sous la forme

$\begin{displaymath} \overline{w'\theta'}=-K_z\left(\frac{\partial \theta}{\partial z}-\gamma_c\right) \end{displaymath}$

(2.56)

avec un contre-gradient $\gamma_c$ positif.

Plusieurs études (notamment par Deardorff lui-même) ont tenté de donner une expression physique de ce contre-gradient mais les modèles de circulation utilisent souvent une valeur constante de l'ordre de 0.5 K/km. La prescription d'un tel contre-gradient est encore de mise dans le modèle du LMD.

Plus récemment, Holtslag et Boville (1993) ont introduit, dans le modèle du NCAR, un schéma de couche limite non local qui inclut un terme de contre-gradient relié directement aux caractéristiques de la couche limite, en suivant une approche développée à l'origine par Troen et Mahrt (1986). Dans cette approche, le coefficient de mélange est prescrit en suivant un profil à la Brost et Wyngaard (1978)

$\begin{displaymath} K_h=\kappa w_h z \left[1-\frac{z}{h}\right]^{1.5} \end{displaymath}$

(2.57)

en utilisant comme vitesse caractéristique pour la couche mélangée une combinaison de la vitesse caractéristique de la couche de surface et de la vitesse convective

$\begin{displaymath} w_m^3=u_*^3+c_1w_*^3 \end{displaymath}$

(2.58)

avec

. Le contre-gradient est calculé de façon à avoir le bon flux convectif au milieu de la couche limite, là où le gradient de température est quasiment nul :

$\begin{displaymath} \frac{\overline{w'\theta'}_0}{2}=\kappa w_h \frac{h}{2} 0.5^{1.5} \gamma_c \end{displaymath}$

(2.59)

D'où l'on tire

$\begin{displaymath} \gamma_c=C \frac{\overline{w'\theta'}_0}{h w_h} \end{displaymath}$

(2.60)

La hauteur de couche limite

elle même est estimée avec une formule en Richardson non local (Eq. 2.42) en rajoutant un excès de température à la température de l'air près de la surface ( $\theta(z_1)$ ). Cet excès de température dépend lui-même du flux de chaleur en surface.

La paramétrisation de Holtslag et Boville (1993) permet d'étendre le transport non local aux traceurs.

Bougeault et Lacarrere (1989) ont également développé une paramétrisation en partie non locale en utilisant une longueur de mélange reliée à la distance qu'une particule d'air, issue du niveau considéré, peut parcourir après qu'on lui a donné au départ une impulsion correspondant à l'énergie cinétique turbulente moyenne dans la couche en question. Notons que cette approche est non locale mais qu'elle ne permet pas le transport en remontant le gradient.

Abdella et McFarlane (1997) ont pour leur part proposé de modifier la fermeture du second ordre de Mellor et Yamada (1974) en paramétrisant les moments du second ordre en utilisant des images en flux de masses comme celles qui sont développées plus loin dans le modèle du thermique. Comme dans le travail de Troen et Mahrt (1986), on aboutit alors à une expression physique d'un terme de contre-gradient.

Toutes ces approches ont en commun qu'elles essaient d'introduire des aspects non locaux dans un formalisme hérité des formulations locales en diffusion turbulente ou dans des développements à des ordre successifs d'équations locales de la turbulence.

Matrices de transilience

Depuis les années 1980, Stull insiste, parfois lourdement, sur la nécessité de rompre radicalement avec la diffusion turbulente. Stull a introduit le concept de matrice de ``transilience", matrice contenant les taux d'échange d'air entre les différentes mailles d'une colonne atmosphérique. Le coefficient $C_{i,j}$ de la matrice est par exemple la concentration massique dans la maille au temps $t+\delta t$ d'un traceur passif injecté uniformément dans la maille au temps en quantité unitaire. Plus que d'une paramétrisation, il s'agit en fait d'un cadre formel dans lequel développer une paramétrisation. Dans ce formalisme, la diffusion turbulente sera représentée par une matrice d'échange tridiagonale (en un pas de temps, une maille n'échange qu'avec les mailles immédiatement au-dessus et au-dessous).

**Figure 2.13:** Calcul et interprétation d'une matrice de transilience à partir d'une simulation des grands tourbillons dans une couche limite convective. D'après Ebert et al. (1989).
$\includegraphics[width=16cm]{therm/FIGURES/SCAN/Erbert.eps}$ L'axe horizontal correspond à l'origine du traceur et l'axe vertical à la destination. Les iso-contours montrent la distribution des traceurs au temps (compté à partir de l'instant auquel les traceurs sont introduits dans la simulation). Le temps caractéristique est de l'ordre de 18 minutes dans cette simulation.

Stull et ses collaborateurs ont analysé les résultats de simulations des grands tourbillons en termes de matrices de transilience. Les résultats montrés sur la Fig. 2.13 ont été obtenus pour une simulation des grands tourbillons d'un cas académique de couche limite convective similaire aux simulations de Moeng et Sullivan (1994) décrites plus haut. Dans ce cas particulier, le temps caractéristique d'advection par les thermiques, , est de l'ordre de 18 minutes. Pour calculer les matrices, on commence par intégrer le modèle sans traceurs jusqu'à atteindre un état de régime de la turbulence. On introduit alors un traceur uniformément dans chaque couche du modèle. Après un certain temps, on peut tracer le profil vertical moyen de la concentration de ce traceur. Les matrices de transilience associées, avec en abcisse la couche d'origine et en ordonnée la couche d'arrivée du traceur, sont montrées sur la Fig. 2.13, ainsi que l'interprétation physique de la forme de la matrice, en haut à gauche. Au-dessus de 1,2 , l'air est très peu affecté par la turbulence. Seuls les termes diagonaux sont donc non nuls. Au bout d'un temps , l'air de la partie haute de la couche limite est descendu lentement ce qui donne de grands termes juste en dessous de la diagonale. A l'opposé, on voit sur la gauche de la matrice que l'air provenant de la couche de surface se répartit rapidement dans la couche mélangée. A , l'air de la couche de surface se retrouve principalement en haut de la couche mélangée. A la fois à et , la matrice de transilience est fortement asymétrique. A , l'air a eu le temps de bien se mélanger et la matrice prend alors une forme plus symétrique, avec des coefficients relativement uniformes dans la couche limite : de l'air originaire de la couche limite se retrouve à peu près équi-réparti dans la couche limite au bout de .

**Figure 2.14:** Forme prise par les matrices de transilience (ou d'échange) dans le cas d'une formule en diffusion turbulente à gauche et dans le cas du modèle de convection asymétrique de Pleim et Chang (1992) à droite. On montre par des flèches sur des colonnes verticales les échanges mis en jeu dans ces paramétrisations et en grisé les éléments non nuls de la matrice associée.
$\includegraphics[width=15cm]{therm/FIGURES/acm.eps}$

La seule fermeture à proprement parler qui semble avoir été proposée dans le cadre des matrices de transilience est le modèle asymétrique de Pleim et Chang (1992). Dans ce modèle, une forme particulière est imposée à la matrice pour rendre compte de l'opposition entre les ascendances thermiques plus concentrées et les subsidences plus lentes. La première couche au-dessus de la surface échange de l'air avec toutes les mailles situées au-dessus tandis que toutes les autres mailles transfèrent de l'air à la maille immédiatement en dessous. On illustre le principe de cette paramétrisation sur la Fig. 2.14. Cette paramétrisation permet d'obtenir un transport de chaleur remontant le gradient de température potentielle, sans autre traitement, dés lors que la première couche du modèle est plus chaude (en température potentielle) que les couches au-dessus.

Les modèles en flux de masse

La dernière catégorie d'approches, celle à laquelle se rattache le modèle du thermique présenté ci-après, est la catégorie des modèles en flux de masse.

Les approches en flux de masse remontent aux travaux de Lilly (1968) et Arakawa et Schubert (1974) et étaient principalement motivées à leur origine par la paramétrisation des mouvements dans les cumulus.

On entend de façon générale par "flux de masse" une approche qui tend à expliciter le transport vertical par des mouvement (ou flux de masse) dans des sous colonnes de la colonne atmosphérique.

Les schémas de Tiedtke (1989) et Emanuel (1991) appartiennent tous les deux à cette catégorie des schémas en flux de masse.

Dans le schéma de Tiedtke, on isole par exemple trois sous colonnes, une pour les ascendances rapides (au c $\oe$ ur des nuages convectifs), une pour les descentes précipitantes (entretenues par l'évaporation des précipitations) et une pour l'environnement. Dans les panaches ascendants et descendants, on calcule explicitement la température, l'eau etc... à chaque niveau, sous des hypothèses de stationnarité, pour ensuite calculer le bilan total dans la colonne.

Dans le cadre de l'étude de la couche limite et de la convection peu profonde, le concept des flux de masse a été beaucoup utilisé d'une façon relativement différente, en développant en général des modélisations simple à fin d'analyse plutôt que de simulation proprement dite, en utilisant une idéalisation des profils verticaux. Ces approches remontent aux travaux de Betts (1973).

Dans ces approches, la couche limite convective est décrite en une ou quelques couches dans les quelles ont suppose que les quantités physiques sont bien mélangées, et des couches infiniment fines de transition au travers des quelles les différentes variables subissent des discontinuités.

Pour une couche limite sèche comme celle décrite plus haut, on aura typiquement une discontinuité de température en surface, une température potentielle unique sur toute la hauteur de la couche limite, une nouvelle discontinuité au niveau de l'inversion, pour se raccorder au profil de la troposphère libre.

Le modèle du thermique

Idéalisation d'une cellule thermique

Le schéma développé ici s'appuie sur une vue idéalisée d'une cellule convective ou d'un thermique proche de celle présentée plus haut.

Considérons un profil de température potentielle typique de la couche limite convective avec une couche de surface (CS sur la Fig. 2.15) instable (typiquement de une à quelques centaines de mètres d'épaisseur), une couche mélangée (CM sur la Fig. 2.15) neutre, surmontée d'une zone stable (couche d'entraînement plus troposphère libre).

**Figure 2.15:** Vue transverse idéalisée d'un rouleau convectif à la base du modèle du thermique.
$\includegraphics[width=14cm]{therm/FIGURES/cell.eps}$

Dans cet environnement idéalisé, le thermique est vu comme un panache d'air chaud montant de la couche de surface sous l'effet de sa flottabilité. On suppose que l'air dans le thermique conserve sa température potentielle virtuelle qui est donc celle de la couche de surface, $\theta_{\mbox{\footnotesize CS}}$ .^2.5

On suppose de plus que le mouvement dans le panache est stationnaire et sans friction, que la vapeur d'eau ne condense pas et que l'air monte sous le seul effet de sa flottabilité. Dans ce cas, l'accélération verticale de l'air, dans le thermique, s'écrit

$\begin{displaymath} \frac{d w}{d t}=w\frac{\partial w}{\partial z}= g \frac{\the... ...a_{\mbox{\footnotesize CM}}}{\theta_{\mbox{\footnotesize CM}}} \end{displaymath}$

(2.61)

Dans ces conditions, l'air est accéléré de manière uniforme jusqu'au niveau où la température potentielle virtuelle dans l'environnement est supérieure à celle de la couche de surface. On retient ce niveau comme définition de $z_{\mbox{\footnotesize i}}$ . A ce niveau, le carré de la vitesse verticale $w_{\mbox{max}}$ est égal à deux fois l'énergie potentielle disponible pour la convection ou CAPE (pour Convective Available Potential Energy)

$\begin{displaymath} \frac{1}{2}w^2_{\mbox{max}}=\mbox{CAPE}=\int_0^{z_{\mbox{\fo... ...\mbox{\footnotesize CM}}}{\theta_{\mbox{\footnotesize CM}}} dz \end{displaymath}$

(2.62)

Au dessus de ce niveau, la vitesse verticale est encore positive (on parle alors d'``overshoot") mais décroît pour finalement s'annuler au niveau $z_{\mbox{max}}$ où

$\begin{displaymath} \int_0^{z_{\mbox{max}}} g\frac{\theta_{\mbox{\footnotesize CS}}-\bar{\theta}}{\bar{\theta}} dz=0 \end{displaymath}$

(2.63)

Cette intégrale est représentée par la zone grisée sur la gauche de la Fig. 2.15. Dans la logique de cette vision très théorique, les particules d'air qui atteignent le point où leur vitesse s'annule peuvent ensuite redescendre puisqu'elles sont alors plus lourdes que l'environnement. On négligera cette possibilité dans les développements qui suivent en supposant que cette zone d'overshoot est aussi une zone de mélange important. En réalité, il semble qu'on observe souvent des descentes d'air assez marquées sur les bords des panaches ascendants, sans doute alimentées en grande partie par l'air provenant de ces overshoots.

Pour calculer le transport convectif, il faut en fait calculer le flux de masse par unité de surface $f=\hat{\alpha}\rho w$ où $\hat{\alpha}$ est la fraction de la surface horizontale couverte par des panaches ascendants (on parlera de couverture fractionnaire). Dans un premier temps, on va supposer que le flux de masse ne dépend pas de l'altitude (panache conservatif) et on va essayer de déterminer ce flux. Pour cela, il est nécessaire de rentrer un peu plus en détail dans la description de la géométrie du thermique.

On va considérer ici la configuration en rouleau, plus simple à analyser, c'est à dire une situation où on a une direction horizontale invariante, par exemple selon l'axe . L'ascendance est alimentée par convergence horizontale dans la couche de surface. Si les effets de friction et de rotation sont négligés, et si on suppose que l'écoulement est stationnaire on a

$\begin{displaymath} v\frac{\partial v}{\partial y}= -\frac{1}{\rho_s} \frac{\partial p_s}{\partial y} \end{displaymath}$

(2.64)

où

est la pression de surface. Si on intègre cette équation sur la largeur de la cellule, on obtient une estimation de la vitesse maximale

$\begin{displaymath} v_{\mbox{max}}^2\simeq \frac{2}{\rho_s}\Delta p_s \end{displaymath}$

(2.65)

où la notation $\Delta$ indique la différence entre les conditions à la base de l'ascendance thermique et dans l'environnement. Si on suppose qu'au sommet de la couche limite, la pression est la même dans l'environnement et dans le panache, on obtient en surface

$\begin{displaymath} p_s=p(z_i)+\int_0^{z_{\mbox{\footnotesize i}}} g\rho dz \end{displaymath}$

(2.66)

et donc, au premier ordre,

$\displaystyle v_{\mbox{max}}^2$	$\textstyle \simeq$	$\displaystyle \frac{2}{\rho_s}\Delta\int_0^{z_{\mbox{\footnotesize i}}} g \rho dz$	(2.67)
	$\textstyle \simeq$	$\displaystyle 2\int_0^{z_{\mbox{\footnotesize i}}} g\frac{\rho}{\rho_s}\frac{\t... ...ize CS}}-\theta_{\mbox{\footnotesize CM}}}{\theta_{\mbox{\footnotesize CM}}} dz$	(2.68)

en supposant toujours que $\Delta \rho/\rho \simeq \Delta \theta/\theta$ . Au coefficient $\rho/\rho_s$ près, de l'ordre de 1 à 0.5 dans la couche limite, cette estimation est exactement celle de la vitesse verticale maximale $w_{\mbox{max}}$ .

D'autre part, en géométrie bidimensionnelle, le flux vertical dans le thermique doit être égal à la convergence horizontale dans la couche de surface :

$\begin{displaymath} w_{\mbox{max}}l(z_{\mbox{\footnotesize i}}) \rho(z_{\mbox{\f... ...}}z_{\mbox{\footnotesize s}}\rho(z_{\mbox{\footnotesize s}}/2) \end{displaymath}$

(2.69)

où

est la largeur du rouleau à la hauteur

. Avec cette approche simple, la largeur de l'ascendance à l'inversion est l'épaisseur de la couche de surface.

Le fait que les vitesses maximales horizontale et verticale soient du même ordre de grandeur implique en retour que la distance entre la subsidence et l'ascendance est de l'ordre de 1. En effet, si la largeur de la cellule était beaucoup plus grande que sa hauteur, une particule partant de la distance atteindrait le sommet de la couche limite avant d'atteindre le panache ascendant et créerait ainsi une ascendance secondaire. Cette isotropie est effectivement observée dans les expériences de Rayleigh-Benard. Dans la paramétrisation, le rapport $r=L/z_{\mbox{max}}$ est utilisé comme paramètre libre. La couverture fractionnaire est minimum au niveau $z_{\mbox{\footnotesize i}}$ où elle vaut $\hat{\alpha}(z_{\mbox{\footnotesize i}}) =l(z_{\mbox{\footnotesize i}})/(r z_{\mbox{max}})$ . On en tire une relation de fermeture pour le calcul du flux de masse dans l'ascendance :

$\begin{displaymath} f=\frac{l w_{\mbox{max}}}{r z_{\mbox{max}}}=\frac{\rho(z_{\m... ...{\mbox{\footnotesize s}}\sqrt{2\mbox{CAPE}}}{r z_{\mbox{max}}} \end{displaymath}$

(2.70)

Les observations et les simulations des grands tourbillons suggèrent des rapports d'aspect compris entre 1 et 10 et le plus souvent de l'ordre de

(Moeng et Sullivan, 1994; Atkinson et Zhang, 1996; LeMone, 1973) pour les rouleaux de la couche limite convective sèche. En fait, le rapport d'aspect en question est le rapport entre la distance séparant deux rouleaux et la hauteur de la couche limite. Ce rapport d'aspect est donc égal à

. La valeur

1 qui sort de l'analyse ci-dessus est donc compatible avec les résultats des simulations des grands tourbillons. Pour avoir méconnu ce facteur 2 dans un premier temps, nous avons retenu

comme valeur nominale pour la paramétrisation et pour la plupart des tests présentés plus loin. La paramétrisation est en fait relativement peu sensible à la valeur de ce paramètre (ce qui veut dire que l'instabilité dans la couche de surface s'ajuste pour obtenir le bon flux de masse dans l'ascendance).

Détraînement et environnement du thermique

On a supposé jusqu'ici que le flux de masse était constant dans le thermique. Cette hypothèse conduit à une couverture fractionnaire infinie au sommet de la couche limite, là où est nul. En pratique, ceci signifie que le thermique ne peut pas rester à flux constant et qu'il doit restituer son air dans l'environnement. Pour la paramétrisation, on suppose que, au-dessus de l'inversion, la largeur du thermique décroît pour s'annuler à $z_{\mbox{max}}$ (en pratique on teste ci-dessous une décroissance linéaire et une décroissance quadratique).

En plus, on permet que l'air soit détraîné du thermique sous l'effet du mélange turbulent en dessous de l'inversion. Ce détraînement est pris en compte sous forme d'un épluchage du thermique. Une analyse d'échelle montre que l'épaisseur de la couche limite d'un jet pénétrant dans un environnement au repos se développe en $\sqrt{\lambda z}$ (voir par exemple Prandtl et Tietjens, 1934, Ch.IV). Dans notre cas, $\lambda =\nu_g /w$ où $\nu_g$ est la viscosité turbulente du gaz. Pour une vitesse typique de 1 m/s et une viscosité turbulente de 10 m s $^{-1}$ , $\lambda =10$ m. On retient cette formulation pour réduire la section du thermique en dessous de $z_{\mbox{\footnotesize i}}$ , en considérant $\lambda$ comme un paramètre ajustable.

Les équations du modèle

La présentation précédente partait d'un profil de température idéalisé, avec notamment une température potentielle constante dans la couche mélangée. Dans la formulation complète, qui doit pouvoir être utilisée pour un profil de température quelconque, le thermique est encore caractérisé par sa vitesse verticale ( $\hat{w}$ ), sa température potentielle ( $\hat{\theta}$ ) et la couverture fractionnaire des thermiques ( $\hat{\alpha}$ ). Le flux de masse dans les ascendances est $\hat{f}=\rho\hat{\alpha}\hat{w}$ ( $\rho$ est supposé identique à l'intérieur et à l'extérieur du thermique). Notons

$\begin{displaymath} \mbox{I}(z,z')=\int_z^{z'} g \frac{\theta(z)-\theta(z'')}{\theta(z'')} dz'' \end{displaymath}$

(2.71)

Pour de l'air instable provenant de l'altitude

(si en

, $\partial \theta/\partial z < 0$ ), la CAPE est

avec $hi(z)=\min\{z'\ \mbox{t.q. }\ z'>z\ \mbox{et}\ \theta(z)=\theta(z')\}$ . On définit de même le sommet du thermique $z_{\mbox{max}}$ comme la valeur maximum de

avec $ht(z)=\min\{z'\ \mbox{t.~q. }\ z'>z\ \mbox{et}\ \mbox{I}(z,z') = 0\}$ .

Le flux de masse $\hat{f}$ est relié au flux d'air entrant $\hat{e}$ et au flux d'air détraîné $\hat{d}$ par :

$\begin{displaymath} \frac{\partial \hat{f}}{\partial z}=\hat{e}-\hat{d} \end{displaymath}$

(2.72)

En suivant l'analyse précédente (et l'Eq. 2.70), on prescrit $\hat{e}$ suivant

$\begin{displaymath} \hat{e}(z)= {\rho(z)}\sqrt{2\mbox{CAPE}(z,h(z))}/\left(r z_{\mbox{max}}\right) \end{displaymath}$

(2.73)

si $\partial \theta/\partial z < 0$ et 0 sinon. En dessous de

, le détraînement est calculé suivant

$\begin{displaymath} \hat{d}(z)=\frac{\partial}{\partial z} \left(\rho \hat{w}\sqrt{\lambda z}/(r z_{\mbox{max}})\right) \end{displaymath}$

(2.74)

Au dessus de

, on réduit la largeur du thermique suivant $[(z_{\mbox{max}}-z)/(z_{\mbox{max}}-z_{\mbox{\footnotesize i}})]^\mu$ . Comme $\hat{e}=0$ à cet endroit, le détraînement entre

et $z_{\mbox{max}}$ s'écrit

$\begin{displaymath} \hat{d}(z)=- \frac{\partial}{\partial z} \left[\rho \hat{w}\... ...}{z_{\mbox{max}}-z_{\mbox{\footnotesize i}}}\right)^\mu\right] \end{displaymath}$

(2.75)

On utilise $\mu=1$ et $\mu =2$ dans la version standard.

Le niveau de l'inversion $z_{\mbox{\footnotesize i}}$ est défini comme l'altitude où la température à l'intérieur du thermique $\hat{\theta}$ devient plus petite que celle de l'environnement (début de l'overshoot). On calcule alors les propriétés de l'air à l'intérieur du thermique suivant l'équation de conservation

$\begin{displaymath} \frac{\partial \hat{f}\hat{\phi}}{\partial z}=\hat{e}\phi_e -\hat{d}\hat{\phi}\end{displaymath}$

(2.76)

où $\phi$ est un scalaire (y compris la température potentielle) et l'indice

est relatif à l'environnement. $\phi_e$ est relié à $\hat{\phi}$ (dans le thermique) et à la valeur moyenne dans la maille $\overline{\phi}$ à un niveau donné

$\begin{displaymath} \overline{\phi}=\hat{\alpha}\hat{\phi}+(1-\hat{\alpha})\phi_e \end{displaymath}$

(2.77)

On recalcule la vitesse verticale moyenne dans le thermique, en tenant compte de sa flottabilité, à partir de l'équation

$\begin{displaymath} \frac{\partial \hat{f}\hat{w}}{\partial z}= -\hat{d}\hat{w}+g\rho\frac{\hat{\theta}-\theta_e}{\theta_e} \end{displaymath}$

(2.78)

L'entraînement n'a pas de contribution à l'Eq. 2.78 dans la mesure où on suppose que l'air entre dans le thermique avec une vitesse verticale nulle.

Pour le transport de la quantité de mouvement, l'idée la plus simple consiste à traiter les deux composantes horizontales et comme des scalaires. C'est cette approche qui est retenue dans la version standard. On teste aussi une version dans laquelle on inclut un échange de quantité de mouvement par les forces de pression exercées entre le thermique et l'environnement. Cet échange va tendre à incliner le thermique dont la vitesse horizontale va progressivement s'approcher de celle de l'environnement. Pour ce test, on utilise la géométrie tridimensionnelle et on applique la formule utilisée pour calculer par exemple la traînée d'un ballon dans l'atmosphère. Cette traînée s'exprime comme le produit du carré du vent apparent par la section apparente et un coefficient de traînée. Dans la géométrie tridimensionnelle, la section apparente du thermique est $L\sqrt{\alpha}$ . Il faut diviser le terme de traînée par pour se ramener à une force par unité de surface horizontale de sorte qu'on obtient finalement

$\begin{displaymath} \frac{\partial \hat{f}\hat{{\bf v}}}{\partial z}=\hat{e}{{\b... ...{\bf v}}_e-\hat{{\bf v}}\vert\vert ({{\bf v}}_e-\hat{{\bf v}}) \end{displaymath}$

(2.79)

$\gamma$ est utilisé comme paramètre ajustable. A noter qu'on peut s'attendre à un freinage moins efficace dans la configuration en rouleaux.

Le transport total d'une variable $\phi$ est finalement donné par la somme du transport vers le haut dans le thermique et du transport vers le bas dans l'environnement avec un flux compensatoire $-\hat{f}$

$\begin{displaymath} \rho\overline{w'\phi'}=\hat{f}\left(\hat{\phi}-\phi_e\right) \end{displaymath}$

(2.80)

Une approximation supplémentaire est finalement introduite pour rendre le schéma plus robuste. En effet, si on intègre telle quelle l'Eq. 2.76, on peut aboutir à des valeurs aberrantes de $\phi_e$ . Considérons l'exemple d'un traceur initialement confiné dans la couche de surface dans une situation convective. L'Eq. 2.76 prédirait $\hat{\phi}>0$ dans la couche mélangée alors qu'on a $\overline{\phi}=0$ de sorte qu'on aurait $\phi_e<0$ . Ce problème provient de l'hypothèse de stationnarité, hypothèse qu'on n'a pas envie de faire tomber, en particulier pour des raisons de simplicité. Une alternative consiste à assimiler l'environnement et les grandeurs moyennes ( $\phi_e=\overline{\phi}$ ). Cette approximation qui est valide pour $\hat{\alpha}\ll 1$ , est en général utilisée dans les schémas de convection profonde. Elle est moins évidente dans notre cas où la couverture fractionnaire des ascendances est typiquement comprise entre 5 et 30 $\%$ . Cette nouvelle approximation est validée a posteriori par le bon comportement de la paramétrisation quand on la compare à des simulations des grands tourbillons.

Le schéma dépend finalement de 4 paramètres ajustables. Pour choisir la valeur du rapport d'aspect on se base sur les observations dans les configurations en rouleaux. Le paramètre $\gamma$ est fixé à 0 par soucis de simplicité. Finalement, les paramètres relatifs au détraînement sont ajustés pour obtenir un bon accord avec les résultats de simulations des grands tourbillons présentés dans la section suivante. Bien que différentes combinaisons soient possibles, nous retenons $\lambda =20$ m et $\nu=2$ .

La discrétisation des équations est donnée par Hourdin et al. (2002).

Fermeture locale

La description ci-dessus ne traite que de la partie méso-échelle du transport de couche limite. Le mélange de petite échelle, important notamment dans la couche de surface, doit être traité avec une autre paramétrisation. Dans les tests montrés par la suite, on utilise pour ces échelles soit la paramétrisation d'origine du modèle du LMD soit celle de Mellor et Yamada. Comme on le verra par la suite, ces paramétrisations en diffusion turbulente, sensées à l'origine traiter l'ensemble de la turbulence de couche limite, voient leur rôle se cantonner naturellement à la couche de surface une fois la paramétrisation du thermique activée, cette dernière stabilisant le profil de température dans la couche mélangée.

Un traitement spécial doit également être appliqué pour les couches limites nocturnes. Les formulations classiques en énergie cinétique turbulente sont connues pour s'effondrer en conditions très stables. Les tests effectués plus loin pour des cas de fort cycle diurne continental confirment ce point. En l'absence de traitement spécifique, on observe la nuit un découplage de la surface et de la première couche d'atmosphère. Pour pallier ce problème, on adopte une formule citée dans un certain nombre de travaux et qui dit que la hauteur de la couche limite doit être au moins égale à $h_{\mbox{min}}=0,07/f$ où est le facteur de Coriolis en s $^{-1}$ . Cette formule ayant la fâcheuse propriété d'être singulière à l'équateur, on prend arbitrairement $f=10^{-4}$ s $^{-1}$ . On en déduit une estimation d'un coefficient de diffusion turbulente minimum ${K_z}_{\mbox{min}}=u^*\kappa z (1-z/h_{\mbox{min}})^2$ .

Comparaison avec les simulations des grands tourbillons

Pour valider la nouvelle paramétrisation, on se base tout d'abord sur des résultats de simulations des grands tourbillons. L'avantage d'utiliser des résultats de simulations des grands tourbillons plutôt que des observations réside dans le fait que les conditions expérimentales sont connues exactement et peuvent être modifiées, ainsi que dans la plus grande facilité d'accès à des diagnostics spécifiques.

Ici, on a retenu les résultats de simulations des grands tourbillons utilisées par Ayotte et al. (1996) dans une étude d'intercomparaison de paramétrisations de la couche limite en conditions neutres ou convectives.

Description des simulations des grands tourbillons

On ne donne ici qu'une description succincte des simulations qui sont décrites en détail par Ayotte et al. (1996). Le code utilisé a été développé à l'origine par Moeng (1984). Il est pseudo-spectral sur l'horizontale et en différences finies sur la verticale. Le domaine est supposé périodique horizontalement. La paramétrisation des tourbillons sous-mailles repose sur une équation pronostique de l'énergie cinétique turbulente.

Les simulations sont effectuées dans des conditions sèches et pour des profils de température neutres ou instables, avec un flux de chaleur prescrit en surface ${\overline{w'\theta'}}_0$ . Neuf simulations, 00WC, 05WC, 00SC, 03SC, 05SC, 24SC, 24F, 15B, 24B, sont utilisées. Elle diffèrent d'abord par la valeur du flux de chaleur en surface (le 05 de 05SC correspond par exemple à ${\overline{w'\theta'}}_0=$ 0.05 K m s $^{-1}$ ) et par le profil initial de température potentielle. Les simulations WC (weakly caped) présentent une inversion peu marquée en sommet de couche limite contrairement aux cas SC (strongly caped). Pour les simulations SC et WC, on impose un vent géostrophique constant de 15 m s $^{-1}$ suivant . La simulation 24F (Free) correspond à un cas de convection libre sans forçage géostrophique alors que les cas 15B et 24B (Baroclines) correspondent à un vent géostrophique de 10 m s $^{-1}$ suivant et variant de 0 en surface à 20 m s $^{-1}$ à 2000 m suivant . Les simulations 24SC et 05SC correspondent à peu près aux simulations B et SB de l'article de Moeng et Sullivan (1994) dont les résultats sont montrés sur la Fig. 2.4.

Toutes les simulations sont initialisées à partir de champs homogènes horizontalement auxquels on ajoute des petites perturbations aléatoires. Une première simulation est effectuée sur une période de 5 $\tau$ où $\tau$ est un temps caractéristique (l'échelle de temps convective de l'ordre de 500 à 1000 s pour les simulations avec un flux non nul en surface). L'état final de cette première simulation est alors utilisé pour initialiser les simulations des grands tourbillons et les simulations uni-colonnes au temps . Les analyses présentées ci-dessous portent sur la période $\left[t_1=t_0+4\tau,t_2=t_0+10\tau\right]$ .

En plus des variables météorologiques, deux traceurs sont pris en compte dans les simulations des grands tourbillons. Les deux sont initialisés au temps avec une discontinuité à l'inversion comme suit:

$\displaystyle B$	$\textstyle =$	$\displaystyle 13.5,\ 0 < z \le 1.01 z_i$
	$\textstyle =$	$\displaystyle 3.0,\ z > 1.01 z_i$	(2.81)
$\displaystyle C$	$\textstyle =$	$\displaystyle 0.0,\ 0 < z \le 1.01 z_i$
	$\textstyle =$	$\displaystyle 1.0,\ z > 1.01 z_i$	(2.82)

Le traceur C est uniquement transporté alors que le traceur B, un analogue à la vapeur d'eau dans l'idée des auteurs, est en plus forcé par un flux en surface. Ce flux est calculé en fonction du contraste entre la concentration dans l'atmosphère au-dessus du sol et une valeur de surface de 15 (unités arbitraires).

Pour ces différentes simulations, l'auteur de l'étude - retrouvé au fil de la toile derrière les rangées d'éoliennes qu'il commercialise en Australie - nous a gracieusement communiqué les profils de début et de fin des simulations des grands tourbillons. Nous disposons à la fois des moyennes des différentes variables du modèle et des moyennes des termes croisés (variances et covariances). En revanche, il n'a pas été possible de récupérer ou de recalculer d'autres grandeurs qui auraient pourtant été intéressantes pour la validation de la paramétrisation tels les moments d'ordre 3 ou les facteurs d'asymétrie.

Les simulations uni-colonnes

On compare les résultats des simulations de Ayotte et al. (1996) avec les résultats d'une version uni-colonne du modèle LMDZ. Dans le modèle uni-colonne, restreint aux processus de couche limite, l'équation donnant l'évolution temporelle des variables , , $\theta$ , et peut s'écrire de façon générique

$\begin{displaymath} \frac{\partial \phi}{\partial t} = {\cal S}_\phi -\frac{1}{\rho}\frac{\partial {\rho\overline{w'\phi'}}}{\partial z} \end{displaymath}$

(2.83)

Le terme source ${\cal S}_\phi$ est le forçage géostrophique pour le vent ( ${\cal S}_u=f \left(v-v_g\right)$ et ${\cal S}_v=-f \left(u-u_g\right)$ ). Il est nul pour les autres variables.

Le flux de surface est prescrit pour le traceur (il est nul) et pour la température potentielle. Pour les vents et le traceur , la condition en surface suit les relations de similitude de Businger-Dyer :

$\begin{displaymath} \frac{\vert\vert{{\bf v}}_1\vert\vert}{u_*}=\frac{1}{\kappa ... ...\frac{z+z_0}{z_0}\right)-\Psi_m\left(\frac{z}{L}\right)\right] \end{displaymath}$

(2.84)

$\begin{displaymath} \frac{B_1-B_s}{B_*}=\frac{1}{\kappa } \left[\ln\left(\frac{z+z_0}{z_0}\right)-\Psi_h\left(\frac{z}{L}\right)\right] \end{displaymath}$

(2.85)

où ${{\bf v}}_1$ et

sont le vent et la concentration du traceur

dans la première couche du modèle. $\kappa =0.4$ est la constante de Von Karman et

est la longueur de Monin-Obukov. Les fonctions $\psi_m$ et $\psi_h$ sont calculées avec les formules données dans la Table 2.1 pour $\gamma_1=\gamma_2=15$ , $\gamma_3=5$ et

Au sein de l'atmosphère, le flux turbulent dans les différents jeux de paramétrisations testés peut s'écrire de façon générique comme

$\begin{displaymath} \rho\overline{w'\phi'}=-\rho K_\phi\left(\frac{\partial\phi}... ...ial z}-\Gamma_\phi\right) +\hat{f}\left(\hat{\phi}-\phi\right) \end{displaymath}$

(2.86)

où $\Gamma_\phi$ est un terme de contre-gradient.

On compare différentes paramétrisations de la couche limite.

MY: $K_\phi$ est calculé suivant le schéma de MY et $\Gamma_\phi=0$ pour toutes les variables. $\hat{f}$ =0.
MY+TH : idem mais en activant le modèle du thermique pour le calcul de $\hat{f}$ et $\hat{\phi}$ . C'est cette version qui est considérée comme le modèle standard de couche limite avec thermiques.
LMD : $K_\phi$ est calculé avec la paramétrisation standard du LMD avec $\Gamma_\phi=0$ et $\hat{f}=0$ .
LMD+CG : comme LMD mais avec en plus un terme de contre-gradient pour la température potentielle.
LMD+AJS : comme LMD à ceci près qu'un ajustement convectif est appliqué à $\theta$ .
LMD+TH : comme LMD avec activation du modèle du thermique.
HB: la paramétrisation de Holtslag et Boville, telle qu'elle est décrite par Ayotte et al. (1996). $\hat{f}=0$ . $K_\phi$ est prescrit en fonction de l'altitude en incluant des aspects non locaux. Le contre-gradient $\Gamma$ est non nul pour $\theta$ , et .

Les simulations de HB et MY sont des duplications de celles présentées par Ayotte et al. (1996) et donnent des résultats très proches de ceux présentés dans l'article.^2.6

Les simulations sont effectuées avec la grille verticale des simulations des grands tourbillons, avec une taille de mailles de 10 ou 20 m suivant les cas. Suivant le cas également, le pas de temps varie entre 15 et 100 s.

Résultats numériques

Figure 2.16: Comparaison des résultats des simulations des grands tourbillons avec ceux obtenus avec le modèle du thermique (MY+TH) en mode uni-colonne dans sa configuration nominale ( $\lambda =20$ , $\mu =2$ ,

et $\gamma =0$ ) pour les cas (05WC, 05SC, 24SC et 24F). Pour chaque cas, on montre le profil initial (pointillés), la moyenne entre les temps

pour le simulations des grands tourbillons (courbe fine) et le modèle du thermique (courbe épaisse). Il n'y a pas de valeur initiale pour le flux de chaleur. Pour le cas 24F, les vents sont nuls.

05WC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05WC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05WC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05WC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05WC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05WC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05WC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )
05SC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05SC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05SC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05SC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05SC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05SC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/05SC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )
24SC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24SC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24SC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24SC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24SC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24SC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24SC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )
24F	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24F/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24F/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24F/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY+TH/24F/C.epsi}$
	$\theta K$	$w'\theta'$	B	C

Sur la Fig. 2.16, on commence par montrer les profils verticaux de température potentielle, de flux de chaleur, de traceurs et de vents pour 4 cas, à la fois pour les simulations des grands tourbillons et pour le modèle du thermique dans sa configuration nominale (MY+TH avec $\lambda =20$ m, $\nu=2$ , et $\gamma =0$ ). On montre les profils initiaux ainsi que les moyennes entre $t_1=t_0+4\tau$ et $t_2=t_0+10\tau$ .

L'accord apparemment très bon entre la paramétrisation et les simulations des grands tourbillons est dû pour une grande part au forçage très contraignant de ces simulations académiques. Comme on l'a expliqué plus haut, avec un un flux de chaleur prescrit en surface et une température homogène dans la couche mélangée, le modèle est essentiellement contraint à chauffer uniformément la couche mélangée et le flux ne peut donc que décroître à peu près uniformément depuis la surface jusqu'à l'inversion.

Figure 2.17: Comparaison des résultats de la paramétrisation HB avec les simulations des grands tourbillons pour les cas 05WC et 24SC.

05WC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/05WC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/05WC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/05WC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/05WC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/05WC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/05WC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )
24SC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/24SC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/24SC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/24SC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/24SC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/24SC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/HB/24SC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )

Figure 2.18: Comparaison des résultats de la paramétrisation MY avec les simulations des grands tourbillons pour les cas 05WC et 24SC.

05WC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/05WC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/05WC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/05WC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/05WC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/05WC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/05WC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )
24SC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/24SC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/24SC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/24SC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/24SC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/24SC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/MY/24SC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )

Figure 2.19: Comparaison des résultats de la paramétrisation LMD avec les simulations des grands tourbillons pour les cas 05WC et 24SC.

05WC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/05WC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/05WC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/05WC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/05WC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/05WC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/05WC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )
24SC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/24SC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/24SC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/24SC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/24SC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/24SC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD/24SC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )

Figure 2.20: Comparaison des résultats de la paramétrisation LMD+CG avec les simulations des grands tourbillons pour les cas 05WC et 24SC.

05WC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/05WC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/05WC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/05WC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/05WC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/05WC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/05WC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )
24SC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/24SC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/24SC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/24SC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/24SC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/24SC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+CG/24SC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )

Figure 2.21: Comparaison des résultats de la paramétrisation LMD+AJS avec les simulations des grands tourbillons pour les cas 05WC et 24SC.

05WC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/05WC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/05WC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/05WC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/05WC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/05WC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/05WC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )
24SC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/24SC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/24SC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/24SC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/24SC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/24SC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+AJS/24SC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )

Figure 2.22: Comparaison des résultats de la paramétrisation LMD+TH avec les simulations des grands tourbillons pour les cas 05WC et 24SC.

05WC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/05WC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/05WC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/05WC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/05WC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/05WC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/05WC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )
24SC	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/24SC/temp.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/24SC/flux.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/24SC/B.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/24SC/C.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/24SC/u.epsi}$	$\includegraphics[width=3.cm,angle=-90]{therm/AYOTTE/SIMUS/LMD+TH/24SC/v.epsi}$
	$\theta K$	$w'\theta'$	B	C	u (m s $^{-1}$ )	v (m s $^{-1}$ )

Les Fig. 2.17 et 2.18 montrent les résultats pour les cas 05WC et 24SC pour les paramétrisations de HB et MY. La paramétrisation de HB se comporte globalement bien dans ce type de conditions (Ayotte et al., 1996). Cette paramétrisation a cependant tendance à exagérer l'inversion. Elle surestime aussi l'entraînement en sommet de couche limite dans le cas 05WC alors qu'elle le sous-estime légèrement pour le cas 24SC. Le schéma de MY sous-estime systématiquement l'entraînement. A noter que le bon accord relatif pour le flux n'est obtenu dans ce cas que grâce à une déstabilisation du profil de température potentielle sans laquelle on ne peut avoir un flux vers le haut. Les figures 2.19, 2.20, 2.21 et 2.22, correspondant respectivement aux simulations LMD, LMD+CG, LMD+AJS et LMD+TH, illustrent l'importance du traitement des aspects non locaux. La paramétrisation du LMD sans aucun traitement aboutit à des atmosphères très instables avec peu d'entraînement au sommet dans le cas 24SC. L'introduction d'un contre-gradient (LMD+CG) résout seulement en partie le problème. Le contre-gradient, imposé ici indépendamment de l'intensité de l'instabilité, est trop fort dans le cas 05WC mais trop faible pour 24SC. L'ajustement sec (LMD+AJS) améliore d'une certaine façon les profils de température potentielle mais sans affecter beaucoup l'entraînement des traceurs en sommet de couche limite. En revanche, l'activation du modèle du thermique (LMD+TH) conduit à des résultats très proches de MY+TH.

Afin de quantifier cette intercomparaison, nous retenons les diagnostics proposés par Ayotte et al. (1996) et qui se concentrent sur l'entraînement en sommet de couche limite. En particulier, ils proposent de calculer pour un scalaire $\phi$ ( $\theta$ , B ou C), le coefficient

$\displaystyle A1$	$\textstyle =$	$\displaystyle - \frac{1}{t_f-t_1}\int_{z_i(t_0)}^H \left[\phi(z,t_f)-\phi(z,t_0)\right] dz$
	$\textstyle =$	$\displaystyle \frac{1}{t_f-t_1}\int_{t_0}^{t_f} \overline{w'\phi'}(z_i(t_0),t) dt$	(2.87)

où

est une altitude située au-dessus du sommet de couche limite. Ce coefficient est facile à calculer en utilisant directement les profils aux instants

. Pour les estimations montrées sur les figures on a retenu comme profil final la moyenne entre

avec

. Pour les traceurs,

est l'aire qui sépare les courbes pleine et pointillée au-dessus de la discontinuité.

Figure 2.23: Coefficients

pour la température potentielle et pour les traceurs

ainsi que l'évolution de la température en surface (valeur moyenne entre

moins celle à

) pour les 9 simulations. On compare les résultats des simulations des grands tourbillons avec la version nominale du modèle du thermique ( $\lambda =20$ m, $\mu =2$ ,

et $\gamma =0$ ) et avec HB et MY.

A1 $_\theta$	A1 $_{\mbox{B}}$
$\includegraphics[width=5cm,angle=-90]{therm/AYOTTE/a1t.eps}$	$\includegraphics[width=5cm,angle=-90]{therm/AYOTTE/a1b.eps}$
A1 $_{\mbox{C}}$	$T_1(\frac{t_1+t_2}{2}) - T_1 (t_0)$
$\includegraphics[width=5cm,angle=-90]{therm/AYOTTE/a1c.eps}$	$\includegraphics[width=5cm,angle=-90]{therm/AYOTTE/ts.eps}$

Figure 2.24: Coefficients

pour la température potentielle et pour le traceur

pour les simulations utilisant la couche limite diffuse du LMD avec différents traitements pour les aspects non locaux (LMD, LMD+CG, LMD+AJS et LMD+TH).

A1 $_\theta$	A1 $_{\mbox{B}}$
$\includegraphics[width=5cm,angle=-90]{therm/AYOTTE/a1tlmd.eps}$	$\includegraphics[width=5cm,angle=-90]{therm/AYOTTE/a1blmd.eps}$

Sur la Fig. 2.23 on montre, pour les neuf cas étudiés par Ayotte et al. (1996), le coefficient pour les simulations des grands tourbillons, pour le modèle du thermique et pour les schémas de HB et MY. La paramétrisation de MY tend à sous-estimer systématiquement l'entraînement. Le schéma de HB surestime l'entraînement pour le cas 05WC et et le sous-estime pour les autres cas. Les deux schémas ne réussissent pas à prédire un entraînement significatif dans le cas de convection libre. Dans l'intercomparaison de Ayotte et al. (1996), on retrouve le même comportement pour tous les schémas testés sauf pour un modèle de ``couche mélangée" qui tend, lui, à surestimer systématiquement l'entraînement. L'accord avec les simulations des grands tourbillons est généralement meilleur pour le modèle du thermique pour les coefficients ainsi que pour la température de l'air près de la surface (dernier panneau de la Fig. 2.23).

Sur la Fig. 2.24, on voit que la couche limite du LMD tend à surestimer l'entraînement pour les cas intermédiaires (00SC, 03SC et 05SC). Cette surestimation est peu sensible au traitement des aspects non locaux. Pour les cas très convectifs, seul le modèle du thermique arrive à effectuer un transport raisonnable. A nouveau, l'entraînement calculé avec LMD+TH est très proche de celui obtenu avec MY+TH (excepté dans les cas intermédiaires pour la température potentielle ou on conserve la surestimation par le modèle du LMD).

Sensibilité aux paramètres

Figure 2.25: Coefficient

normalisé pour la température potentielle et le traceur B. On montre les résultats des simulations des grands tourbillons et le modèle du thermique à la fois dans sa configuration nominale ( $\lambda =20$ m, $\mu =2$ ,

et $\gamma =0$ ) et pour les expériences de sensibilité dans lesquelles on a changé la valeur d'un ou deux paramètres. Le cas ``incliné" correspond à $\gamma =0.5$ (se reporter au texte pour plus de détails). Les lignes épaisses (simulations des grands tourbillons et modèle du thermique en configuration nominale) sont déjà montrées sur la Fig. 2.23.

A1 $_\theta$	A1 $_{\mbox{B}}$
$\includegraphics[width=6cm,angle=-90]{therm/FIGURES/a1t2.eps}$	$\includegraphics[width=6cm,angle=-90]{therm/FIGURES/a1b2.eps}$

Comme on l'a dit plus haut, les résultats des simulations des grands tourbillons ont été utilisés pour sélectionner les valeurs nominales des paramètres pour le détraînement. La figure Fig. 2.25 montre le paramètre pour la température potentielle et le traceur B et pour différentes valeurs des paramètres ( $\lambda =80$ m, $\lambda =0$ m, 5 et $\mu=1$ ). Quand on augmente $\lambda$ , on diminue dans les cas convectifs. A l'opposé, sans épluchage en dessous de l'inversion ( $\lambda =0$ ), ou quand on utilise une décroissance moindre au-dessus de l'inversion ( $\mu=1$ ), l'entraînement est surestimé. Cependant, une simulation avec $\lambda =80$ m et $\mu=1$ produit un accord avec les simulations des grands tourbillons aussi bon que le cas nominal. L'entraînement est relativement peu sensible au rapport d'aspect . Enfin, quand on prend en compte l'échange de quantité de mouvement entre le thermique et l'environnement ( $\gamma =0.5$ ), on ne modifie pas non plus beaucoup l'entraînement.

Pour finir, remarquons que, même pour des valeurs des paramètres assez éloignées des valeurs nominales, le modèle du thermique conserve une bonne sensibilité aux flux de surface et aux profils initiaux. Par exemple, pour $\theta$ , le coefficient pour est toujours un peu plus grand pour le cas 24F que pour le cas 24SC; est aussi systématiquement plus grand pour le cas 05SC que pour le cas 05WC. Enfin, pour le traceur , les résultats sont encore moins sensibles aux paramètres du modèle (les résultats sont similaires pour ).

Dans le thermique

Figure 2.26: Structure du thermique pour le cas 24SC. Pour le premier graphe, la courbe fine (resp. épaisse) correspond à la surface fractionnaire avant (resp. après) épluchage. Pour $\overline {w'\theta '}$ (courbe épaisse sur le dernier graphe) on montre la décomposition entre le thermique (courbe fine) et la partie diffuse (tireté).

$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/24SCfraca.epsi}$	$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/24SCwam.epsi}$	$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/24SCfm.epsi}$	$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/24SCwt.epsi}$
$\alpha$	$\hat{w}$	$\hat{f}$	$\overline {w'\theta '}$
$\%$	m s $^{-1}$	kg m $^{-2}$ s $^{-1}$	K m s $^{-1}$

La Fig. 2.26 montre la structure du thermique telle qu'elle est simulée avec la paramétrisation nominale ainsi que la décomposition du flux de chaleur entre la partie du thermique à proprement parler (courbe fine) et la partie diffuse paramétrisée avec le schéma de MY (courbe pointillée). La chaleur est d'abord transférée dans la couche de surface par le schéma de MY puis répartie dans la couche limite par le thermique. Dans la couche mélangée, le thermique transporte la chaleur vers le haut malgré une atmosphère légèrement stable.

A remarquer que, quand elle est utilisée seule, la paramétrisation de MY est active également dans la couche mélangée du fait de la déstabilisation du profil de température.

Les thermiques couvrent environ 10 $\%$ de la surface. On obtient des valeurs similaires pour les autres cas. A noter que le bon accord en termes d'entraînement qu'on obtient avec correspond à des couvertures beaucoup plus faibles (3-5 $\%$ ). Ces valeurs sont en accord avec les analyses des simulations des grands tourbillons (Moeng et Sullivan, 1994) (qui montrent des fractions de 10-20 $\%$ avec un rapport d'aspect de 2-3). Les mesures avions dans des conditions similaires conduisent à des résultats semblables (Williams et Hacker, 1993). Les valeurs sont cependant dans la tranche basse des observations ce qui peut être relié au fait que nous ne considérons comme ascendance que l'air qui provient directement de la couche de surface sous l'effet de sa flottabilité.

Moments du second ordre

Figure 2.27: Moments d'ordre deux de la température potentielle et du vent vertical pour les cas 05WC et 24SC pour les simulations des grands tourbillons (courbes fines) et la paramétrisation nominale (courbes épaisses). Les courbes pointillées correspondent à la contribution du thermique seule.

$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/05WCtt.epsi}$	$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/05WCww.epsi}$
$\overline{{\theta'}^2}/{\theta_*}^2$ , 05WC	$\overline{{w'}^2}/{w_*}^2$ , 05WC
$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/24SCtt.epsi}$	$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/24SCww.epsi}$
$\overline{{\theta'}^2}/{\theta_*}^2$ , 24SC	$\overline{{w'}^2}/{w_*}^2$ , 24SC

La comparaison des moments d'ordre 2 aux résultats des simulations des grands tourbillons permet de mieux comprendre le comportement de la paramétrisation et la physique de la couche limite convective.

Cette comparaison pour les cas 05WC et 24SC est montrée sur la Fig. 2.27. Dans les résultats de simulations des grands tourbillons fournis par Ayotte (courbes épaisses), les moments d'ordre deux contiennent la partie explicite plus la partie paramétrisée pour le vent et seulement la partie explicite pour la température. Pour la paramétrisation, on montre à la fois la contribution des thermiques (tiretés)

$\displaystyle \overline{{\phi'}^2}$	$\textstyle =$	$\displaystyle \alpha {(\hat{\phi}-\phi)}^2+(1-\alpha ){(\phi_e-\phi)}^2$	(2.88)
	$\textstyle =$	$\displaystyle \frac{\alpha }{1-\alpha }{(\hat{\phi}-\phi)}^2$	(2.89)

et la somme de cette contribution avec la contribution de la fermeture locale (courbes fines).

Dans la couche de surface, les fluctuations de vent et de température sont bien reproduites. La légère surestimation de $\overline{{\theta'}^2}$ est sans doute due au fait que les résultats des simulations des grands tourbillons n'incluent pas la partie paramétrisée.

Dans la couche mélangée, pour $0.1 \le z/z_i \le 0.6$ , la prédiction de $\overline{{\theta'}^2}$ par le modèle du thermique est aussi très proche des résultats des simulations des grands tourbillons. En revanche, $\overline{{w' }^2}$ est sous-évalué d'un facteur 2 environ. Ceci peut se comprendre de la façon suivante. Dans la couche mélangée, nous tenons compte seulement de la partie liée aux structures méso-échelles, en nous basant sur une vue idéalisée d'un thermique homogène en vent et en température comme l'illustrent les graphes du haut de la Fig. 2.28. Dans la réalité, l'air dans le thermique (et dans l'environnement mais dans une moindre mesure) est également turbulent à plus petite échelle. Cependant, comme on l'a dit plus haut, les fluctuations de température associées avec ces fluctuations de petite échelle de sont petites parce que $\theta$ est relativement homogène verticalement, à la fois dans le thermique et dans l'environnement. Avec une vue lagrangienne des choses, on peut dire qu'une fluctuation positive de la température est attribuée à une particule de fluide qui provient de la couche de surface. Dans une limite non visqueuse, que la trajectoire empruntée par la particule soit une ligne droite ou qu'elle s'apparente davantage à une marche aléatoire importe peu. La vue un peu moins simpliste du thermique dans laquelle on superpose de la turbulence de petite échelle au structures organisées est illustrée sur les graphes du milieu de la Fig. 2.28. Puisque les fluctuations à haute fréquence de sont seulement partiellement associées à des fluctuations de $\theta$ , les deux visions conduisent à des flux de chaleur comparables. C'est bien ce que suggèrent les analyses en composite de Williams (1991) présentées plus haut, dont deux exemples sont reproduits en bas de la Fig. 2.28.

Figure 2.28: Des vues schématiques de sondages horizontaux de

et $\theta '$ pour le thermique idéalisé sous-tendant la présente paramétrisation (figures du haut) et pour une vision un peu plus réaliste incluant des tourbillons de petite échelle (figures du milieu) sont comparées à des ``composites" obtenues à partir de vols avions (figures du bas). Les ``composites" sont ceux de Williams et Hacker (1992), Fig. 5, et correspondent à $0.2\le z/z_i \le 0.3$ . La courbe pleine représente la structure moyenne des thermiques. La courbe tiretée montre la variance associée aux différents thermiques échantillonnés pour construire l'image moyenne. Dans l'analyse en composites, les thermiques sont superposés sur la même échelle de longueur.

$\includegraphics[width=9.cm]{therm/FIGURES/scnd.eps}$

En revanche, dans le haut de la couche mélangée et dans la zone d'entraînement, le modèle ne reproduit correctement ni les moments d'ordre deux de la température ni ceux du vent. En fait, les processus en jeu dans cette région sont beaucoup plus compliqués. Ils sont résumés ici dans l'épluchage du thermique. La réduction de l'aire fractionnaire au-dessus du niveau d'inversion cache en fait des particules qui, après overshoot, redescendent pour être mélangées dans l'environnement plus bas que leur niveau d'excursion maximale. Si on prenait en compte ces aller-retours de façon plus explicite, on pourrait aboutir à des flux de chaleur identiques mais pour des valeurs de $\overline{\theta'}^2$ très supérieures. On peut penser de même que les ondes de gravité, excitées en particulier par la percussion des panaches thermiques dans la couche d'inversion, contribuent peu au flux de chaleur malgré des fluctuations importante de température (dues au fort gradient vertical de température potentielle à cette endroit). Le fait que, dans les simulations des grands tourbillons, on ait encore des fluctuations importantes de la température pour , à une altitude où $\overline{w'\theta'}=0$ , souligne le fait que seulement une part des fluctuations turbulentes est impliquée dans des processus d'échange.

Transport de quantité de mouvement

A titre d'illustration, on compare sur la Fig. 2.29 le vent horizontal simulé avec la paramétrisation nominale et avec celle incluant un terme de traînée (on choisit $\gamma =0.5$ ). Dans la version standard, le vent est constant dans le thermique au-dessus de la couche de surface. Quand on introduit la traînée, le thermique s'incline et le vent converge vers la valeur dans l'environnement. L'impact n'est pas vraiment positif. C'est le genre de raffinements auxquels on pourrait penser si on avait accès à des diagnostics plus détaillés des simulations des grands tourbillons. Il faut cependant se souvenir que l'impact sur les flux de chaleur et de traceurs est faible.

Figure 2.29: Composante horizontale du vent

pour le modèle du thermique dans sa configuration nominale et quand on inclue une traînée (respectivement $\gamma =0$ et 0.5). Les figures montrent, pour les cas 05WC et 24SC, le profil initial, le profil moyen obtenu avec les simulations des grands tourbillons, le résultat des paramétrisations ainsi que la valeur simulée dans le thermique ( $\hat{u}$ ).

$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/u05WCKref.epsi}$	$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/u05WCKrefV05.epsi}$
(m s $^{-1}$ ), 05WC, $\gamma =0$	(m s $^{-1}$ ), 05WC, $\gamma =0.5$
$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/u24SCKref.epsi}$	$\includegraphics[width=4.5cm,angle=-90]{therm/FIGURES/u24SCKrefV05.epsi}$
(m s $^{-1}$ ), 24SC, $\gamma =0$	(m s $^{-1}$ ), 24SC, $\gamma =0.5$

Simulations uni-colonnes de la POI 2 d'ESQUIF

On illustre l'impact de la nouvelle paramétrisation dans une configuration moins académique en utilisant le modèle de circulation générale en version unicolonne. Dans cette configuration, la couche limite est couplée aux autres processus comme le rayonnement et la résolution verticale est plus proche des configurations classiques des modèles de circulation. Pour ce faire, on s'appuie sur un cas de convection sèche bien développée, observé dans la région parisienne pendant la campagne ESQUIF, plus précisément lors de la seconde Période d'Observation Intensive de cette campagne (POI 2), en août 1998.

La POI 2 d'ESQUIF

ESQUIF (``Étude et Simulation de la QUalité de l'air en Ile de France") est un projet dédié à l'étude et à la modélisation de la qualité de l'air dans la région parisienne (Menut, 2000). La POI 2 d'ESQUIF correspondait à une période très chaude, avec des vents faibles et sans nuages, c'est-à-dire un cas idéal pour tester le modèle du thermique. La POI 2 s'est déroulée du 7 au 11 août 1998, et précédait de peu un très important pic de pollution à l'ozone dans la région parisienne.

Pour la validation du modèle, on a utilisé les observations de la température, de l'humidité relative et de la pression à 2 m au-dessus du sol enregistrées toutes les heures par le réseau Mesonet de Météo-France déployé autour de Paris. On a aussi utilisé les sondages de Trappes. En plus des sondages effectués en routine à 00:00 et 12:00 (heures TU comme toutes les heures utilisées dans cette section), nous disposions de sondages additionnels effectués toutes les 3 heures. La Fig. 2.30 illustre l'évolution de la structure de la couche limite pendant la POI. La hauteur de la couche limite est marquée par des points sur des profils de température potentielle. Cette hauteur est déterminée à partir d'un seuil sur un nombre de Richardson non local, $R_{ib}> 0.21$ , calculé suivant l'Eq. 2.42. Pour les sondages, on considère comme altitude de référence le premier point disponible sur les profils (situé entre 0 et 50 m suivant les cas).

Figure 2.30: Température potentielle virtuelle enregistrée à Trappes les 8 et 9 août 1998 et hauteur de la couche limite (points; voir texte). Les profils sont décalés suivant l'heure de lâché du ballon : le temps correspond au point de surface de chaque profil et la relation entre temps et température est de 2.75 K/h.

$\includegraphics[width=8cm]{therm/FIGURES/rstzidecal.epsi}$

Les 8 et 9 août correspondent à des conditions typiques de couche limite convective développée avec, la nuit, une couche limite écrasée dans les 300 premiers mètres et, l'après-midi, une température bien mélangée jusque 2300 m le 8 et 2800 m le 9 août.

Les simulations unidimensionnelles

Les simulations sont effectuées sur la période du 7 au 9 août en partant des sondages de 23:30, le 6. Le 7, qui ne présentait pas une couche convective bien développée, est utilisé comme période de mise en régime pour les paramétrisations. Les champs météorologiques sont à nouveau initialisés avec les sondages du 8 août à 5:30. Les résultats sont présentés pour les 8 et 9.

Dans les simulations, l'albédo de surface est fixé à 0.19 et la longueur de rugosité à 0.4 m. L'inertie thermique de la surface est ajustée à 3000 J m $^{-2}$ s $^{-1/2}$ K $^{-1}$ , valeur typique des surfaces continentales, de manière à reproduire correctement l'amplitude du cycle diurne de la température de surface. La température initiale du sol est fixée à 292 K dans les 11 niveaux du modèle de sub-surface. Enfin, l'aridité du sol $\beta$ , qui relie l'évaporation à l'humidité spécifique selon $E=1.35 \beta \rho {u_*}^2 [q_s(T_s)-q]$ (Laval et al., 1981), ne varie pas au cours du temps et est utilisée comme paramètre ajustable des simulations. On prend finalement $\beta=0.08$ pour obtenir une évolution de l'humidité en surface comparable à celle observée sur les 3 jours.

Figure 2.31: Comparaison des mesures mesonet (grisés, 2 m au-dessus de la surface), des sondages (point) et du modèle (courbes) concernant l'évolution sur les 8 et 9 août de la température d'une part (K, graphique du haut) et de l'humidité spécifique (en g/kg, graphique du milieu). Pour les sondages, on montre la valeur moyenne sur une hauteur correspondant à la première couche du modèle. Le graphique du bas montre, avec les mêmes conventions, la hauteur de la couche limite calculée comme pour la Fig. 2.30.

$\includegraphics[width=10cm]{therm/FIGURES/ts_modrst.epsi}$

Malgré des vents faibles, il faut prendre en compte le forçage grande échelle dans ces simulations. Pour la vapeur d'eau, et parce que la POI était dépourvue de nuages, l'humidité spécifique du modèle n'est affectée que par l'évaporation en surface et la paramétrisation de la couche limite. Dans l'atmosphère libre, on observe un assèchement systématique pendant le 8 août qui ne peut s'expliquer que comme un effet de la grande échelle. Il peut s'agir par exemple d'une advection d'air sec et chaud depuis le sud. Afin de garder un forçage aussi simple que possible, on prescrit ce forçage comme une subsidence

$\begin{displaymath}\frac{\partial q}{\partial t}= -w \frac{\partial q}{\partial p} \end{displaymath}$

(2.90)

avec $w=w_0 \times \sin\left(2 \pi p/p_s\right)$ et

=0.6 Pa s $^{-1}$ - correspondant à une vitesse d'environ 1cm/s dans la troposphère moyenne -. Avec ce forçage, on obtient un assèchement correct au court de la journée du 8 août.

En principe, la dérivation est moins directe pour la température du fait de l'importance du forçage radiatif. Cependant, comme le forçage radiatif est sans doute assez bien simulé pendant cette période sans nuages, on retient la même approche. Les résultats montrent qu'il faut rajouter un terme de chauffage par la grande échelle pour obtenir un accord satisfaisant dans la moyenne atmosphère. On retient la même forme que précédemment mais avec Pa s $^{-1}$ pour le 8 août et pour le 9.

Dans les simulations que nous présentons ici, le vent horizontal est directement forcé par des valeurs interpolées linéairement entre deux radiosondages. On utilise une discrétisation verticale avec 40 couches, typique des résolutions qui seront utilisées avec le modèle de circulation générale dans les années qui viennent. Le premier niveau est centré à environ 40 m au-dessus de la surface et le 15 $^{\mbox{ème}}$ niveau est situé vers 4.4 km avec une résolution verticale de 500 m entre 1.5 et 3.5 km. Le modèle est intégré avec un pas de temps de 3 minutes.

Résultats

nous présentons tout d'abord des comparaisons effectuées avec le modèle du thermique dans sa configuration nominale et les schéma de MY et HB. Sur la Fig. 2.31 on compare la température et l'humidité spécifique simulées dans le premier niveau du modèle avec les sondages de Trappes au même niveau et les mesures Mesonet à 2 m. A noter que le forçage du modèle a été réglé de façon à reproduire correctement la température de surface (dérive et amplitude du cycle diurne) et la dérive de l'humidité de surface observée à Trappes. Les mesures Mesonet, à 2 m au-dessus du sol, donnent une idée de la stabilité de la couche de surface et des fluctuations des champs météorologiques près de la surface.

La Fig. 2.32 montre une comparaison des sondages de Trappes avec les résultats du modèle pour la température potentielle virtuelle $\theta_{v}$ (K) et l'humidité spécifique (g/kg) pour le 9 août à 5:30 et 17:30.

**Figure 2.32:** Température potentielle virtuelle et humidité spécifique pour le 9 août.
$\includegraphics[width=8cm]{therm/FIGURES/s1prof_9A98.epsi}$

Pour les sondages de 5:30, les profils montrent deux couches d'inversion. La plus basse correspond au sommet de la couche limite nocturne avec un très fort gradient à 500 m au-dessus de la surface. La couche limite nocturne est bien représentée avec les 3 paramétrisations. Ceci n'est pas surprenant dans la mesure où on a modifié la paramétrisation de MY (aussi bien quand elle est utilisée seule qu'avec le modèle du thermique) pour les conditions très stables en suivant Holtslag et al. (1990) comme on l'explique plus haut. La couche résiduelle (vers 2200 m) est mieux simulée avec le modèle du thermique du fait d'une meilleure simulation de la couche convective le 8. Au-dessus de cette inversion résiduelle, les profils sont essentiellement déterminés par le forçage de grande échelle.

A 17:30 les profils sont le reflet des processus convectifs du début d'après-midi. Le modèle du thermique montre un accord bien meilleur avec les sondages, avec une hauteur de couche limite de 2500 m alors que les schémas de HB et MY la trouvent vers 2000 m. La surestimation de au-delà de 2700 m n'est due qu'au forçage de grande échelle.

Figure 2.33: Sensibilité aux paramètres : température potentielle et humidité spécifique à 17:30 pour le 9 août obtenues pour différentes valeurs des paramètres du modèle du thermique.

$\includegraphics[width=8cm]{therm/FIGURES/s1sty_prof_therm_9A98_1730.epsi}$

On effectue les mêmes tests de sensibilité que ceux présentés pour la comparaison aux simulations des grands tourbillons (Fig. 2.33). La sensibilité est généralement moindre. A noter en particulier que le détail de la façon dont on décrit le détraînement au-dessus de l'inversion importe peu, ce qui ce comprend aisément vue la résolution verticale très grossière à ce niveau.

Sensibilité du transport des traceurs à la paramétrisation de la couche limite

**Figure 2.34:** Grille utilisée pour les simulations guidées et zoomées sur la France avec l'emplacement du SIRTA et des stations radon de HD et JFJ.
$\includegraphics[width=14cm]{therm/FIGURES/carte.eps}$

Figure 2.35: Evolution, pour la période du 1 au 20 août 1998, de la concentration de radon en surface à Heidelberg (en haut, Bq m $^{-3}$ ), et, dans la première couche du modèle, de la température (K), de l'humidité relative (%) et de l'humidité spécifique (g/kg). Pour les trois courbes du bas, les données correspondent aux réanalyses ERA40 interpolées dans la couche correspondante.

$\includegraphics[angle=-90,width=15cm]{therm/FIGURES/evtrac01.eps}$

$\includegraphics[angle=-90,width=15cm]{therm/FIGURES/evtemp.eps}$

$\includegraphics[angle=-90,width=15cm]{therm/FIGURES/evrhum.eps}$

$\includegraphics[angle=-90,width=15cm]{therm/FIGURES/evovap.eps}$

Jour (août 1998)

On présente dans cette section des tests relatifs au transport des traceurs avec les nouveaux schémas de couche limite. On se concentre sur le mois d'août 1998, pour lequel on dispose à la fois des sondages de la POI2 d'ESQUIF et de mesures de Radon en continu en Europe, montrant un fort cycle diurne pour la même période. On utilise pour ce faire des simulations tridimensionnelles, zoomées sur le nord de la France et guidées par les réanalyses ERA40. Pour le vent, on utilise des constantes de temps de guidage différentes à l'extérieur (2h30) et à l'intérieur (1 jour) du domaine zoomé - en pratique cette constante est spécifiée comme une fonction de la taille de la maille considérée -. Pour la température et l'humidité relative, on applique une constante de temps uniforme de 1 jour. Le guidage est donc très peu contraignant pour le modèle à l'intérieur du maillage. Le guidage plus important du vent à l'extérieur du domaine garantit une bonne représentation de l'advection à grande échelle, en phase avec la situation synoptique observée. Sur la Fig. 2.34 on montre le maillage, la position du SIRTA et des deux stations de mesure du radon utilisées par la suite, l'une située à Heidelberg (HD), à basse altitude, et l'autre située au sommet du Jungfraujoch (JFJ), à 3400 m d'altitude.^2.7 ^2.8

Comme pour les simulations unidimensionnelles présentées ci-dessus, il est tout d'abord nécessaire d'ajuster l'inertie thermique et l'humidité du sol. Des tests discutés dans la conclusion de ce chapitre montrent que les modèles de sol à notre disposition (y compris le modèle ORCHIDEE) ne sont pas capables de maintenir des concentrations d'eau correctes en l'absence de pluie. On choisit donc d'imposer l'humidité du sol, ce qui revient à imposer le facteur $\beta=E/E_{\mbox{pot}}$ . On choisit pour le mois d'août $\beta=0,133$ , c'est-à-dire 10/75 où 10 mm est la hauteur d'eau dans le sol et 75 mm la hauteur à partir de laquelle $\beta=1$ (se reporter à la description des paramétrisations physiques de LMDZ, Section ).

On montre sur la Fig. 2.35 les résultats d'une simulation de référence utilisant la version nominale du modèle du thermique (MY+TH). On voit que le modèle reproduit correctement à la fois l'évolution sur le mois et le cycle diurne de la concentration de radon mesurée à Heidelberg (HD) et des variables météorologiques en région parisienne. Pour les variables météorologiques, on compare sur cette figure les valeurs obtenues dans la première couche du modèle aux sorties des réanalyses ERA40 du ECMWF interpolées à la même altitude. Le réglage de l'inertie thermique (à 1700 USI) et de l'humidité du sol (à 10 mm d'eau) a été fait en privilégiant la première dizaine de jours. On voit qu'avec ce réglage, on sous-estime le cycle diurne après le 13 août, et ce à la fois pour la température et le radon.

**Figure 2.36:** Cycle diurne moyen du radon en surface à Heidelberg et, dans la première couche du modèle, de la température, de l'humidité relative et spécifique dans la première couche du modèle (période du 7 au 10 août).
$\includegraphics[angle=-90,width=8cm]{therm/FIGURES/trac01.eps}$ $\includegraphics[angle=-90,width=8cm]{therm/FIGURES/temp.eps}$ $\includegraphics[angle=-90,width=8cm]{therm/FIGURES/rhum.eps}$ $\includegraphics[angle=-90,width=8cm]{therm/FIGURES/ovap.eps}$

**Figure 2.37:** Evolution de la température et de l'humidité spécifique à Trappes et pour trois des simulations pour les 7, 8 et 9 août, journées pour lesquelles on dispose des radiosondages à haute fréquence. Les valeurs simulées correspondent à la première couche du modèle et sont donc à comparer aux données de Trappes à 35 m.
$\includegraphics[width=14cm,clip]{therm/FIGURES/trappes.eps}$

**Figure 2.38:** Influence de l'humidité du sol. Cycles diurnes du radon, de l'humidité et de la température (période du 7 au 10 août 1998) pour la simulation nominale (MYTH, contenu en eau du sol de 10 mm) et des simulations avec un contenu divisé (MYTHQS5) ou multiplié (MYQS20) par 2.
$\includegraphics[angle=-90,width=8cm]{therm/FIGURES/trac01QS.eps} \includegraphics[angle=-90,width=8cm]{therm/FIGURES/tempQS.eps}$ $\includegraphics[angle=-90,width=8cm]{therm/FIGURES/rhumQS.eps} \includegraphics[angle=-90,width=8cm]{therm/FIGURES/ovapQS.eps}$

On privilégie dans ce qui suit la période du 7 au 10 janvier qui correspondait à un fort cycle diurne estival, période au milieu de laquelle on bénéficie des radiosondages de la POI2 d'ESQUIF. La Fig. 2.36 montre pour cette période le cycle diurne moyen des différentes quantités montrées sur la Fig. 2.35 pour les simulations LMD(+CG), LMD+TH, MY et MY+TH. La température est bien réglée pour l'ensemble des simulations. Le cycle diurne du radon est relativement bien reproduit avec les différents modèles. Le radon, accumulé près de la surface la nuit dans la couche limite nocturne, se dilue au cours de la matinée en même temps que se développe la couche limite convective. On voit cependant que le cycle diurne des concentrations de surface est davantage sensible à la paramétrisation en diffusion qu'à l'utilisation ou non du modèle du thermique (qui tend à creuser un peu plus ce cycle diurne). A noter que des différences plus importantes étaient rapportées par Idelkadi (2002) mais qu'on s'est rendu compte par la suite que ces différences étaient dues en fait à des différences cachées dans la formulation de la diffusion turbulente.

Le cycle diurne observé se situe quelque part entre celui obtenu avec les schémas de MY et du LMD. Comme pour le radon, l'humidité est minimum près de la surface dans l'après-midi du fait du mélange vertical dans la couche limite convective. Cet effet domine largement celui du renforcement de l'évaporation dans la journée. L'apparent désaccord avec les données ERA40 semble en fait plutôt imputable aux réanalyses. En effet, la comparaison directe aux observations, pour les trois jours où nous disposons de radiosondages toutes les 3 heures (Fig. 2.37) est nettement plus favorable. Le cycle diurne est un peu trop faible pour le modèle LMD mais très bien représenté pour la simulation nominale MY+TH.

La Fig. 2.38 montre le cycle diurne moyen du radon et des variables météorologiques pour cette période et pour trois simulations avec différentes valeurs de l'humidité du sol. On voit que l'air est à la fois trop sec et trop chaud si on divise par deux la valeur du contenu en eau du sol. On obtient de façon symétrique un air trop humide et trop froid pour une humidité du sol deux fois plus grande. La version nominale est bien réglée à la fois en humidité et en température.

**Figure 2.39:** Cycle diurne moyen (du 7 au 10 août 1998) du profil vertical (axe des ordonnées en Pa) de radon (Bq m $^{-3}$ ) simulé à Heidelberg avec pour la diffusion turbulente, soit la paramétrisation LMD(+CG) (à gauche) soit la paramétrisation MY (à droite). On montre les simulations sans (en haut) et avec (en bas) thermiques ainsi que la différence relative entre les deux (avec moins sans).
$\includegraphics[width=16cm]{therm/FIGURES/RNz.eps}$

Si le cycle diurne des concentrations de surface est finalement relativement peu sensible à l'introduction du modèle du thermique, il n'en va pas de même en altitude. On montre sur la Fig. 2.39 les cycles diurnes moyens sur la verticale. On voit bien la montée plus rapide du radon sur les différences relatives, entre 9:00 et 15:00 UTC. En moyenne aussi, les simulations avec thermiques montrent des profils verticaux plus marqués avec moins de radon près de la surface et plus en sommet de couche limite (signature d'un épaississement de cette couche limite). Il est intéressant de noter que l'impact est relativement similaire pour les deux modèles de diffusion turbulente (LMD et MY).

**Figure 2.40:** Evolution sur le mois d'août 1998 de la concentration de radon observée à la station Jungfraujoch, située à 3400 m d'altitude, et simulée au même point et dans la couche du modèle située à cette même altitude.
$\includegraphics[width=17cm]{therm/FIGURES/JFJ.eps}$

Figure 2.41: Profils d'humidité spécifique simulés et observés à Trappes.

8 août,	9 août,	9 août,	9 août
$\includegraphics[width=3.7cm]{therm/FIGURES/8A17h.eps}$	$\includegraphics[width=3.7cm]{therm/FIGURES/9A05h.eps}$	$\includegraphics[width=3.7cm]{therm/FIGURES/9A14h.eps}$	$\includegraphics[width=3.7cm]{therm/FIGURES/9A17h.eps}$

Il existe malheureusement apparemment peu de données permettant de départager les modèles, même pour des différences aussi importantes. Les campagnes d'observation du radon en altitude sont relativement limitées et ont souvent été effectuées dans des conditions particulières, notamment dans des régimes de brises de mer sur les côtes. Les résultats de modélisation (des tests, non montrés, ont été effectués avec LMDZ pour certains de ces cas) montrent que les autres sources d'incertitudes sont souvent trop importantes pour apporter des réponses.

On montre cependant deux indications du meilleur comportement du modèle avec thermiques. C'est d'abord (Fig. 2.40) l'évolution de la concentration de radon pour le mois d'août à la station de Jungfraujoch. On peut penser que cette station, située sur un sommet relativement isolé, à 3400 m d'altitude, est souvent sensible au niveau moyen de radon à cette altitude plutôt qu'à des effets locaux. On compare sur la figure observations et résultats de simulations dans la couche du modèle correspondant à l'altitude de JFJ. On voit d'abord que le modèle reproduit raisonnablement l'évolution au cours du mois de l'ordre de grandeur des concentrations de radon. Les simulations sont très proches les unes des autres. On observe cependant des différences pour les 8, 9 et 14 août (et dans une moindre mesure le 7). Pour ces trois jours, les modèles avec thermiques (LMD+TH et MY+TH) prédisent des niveaux de radon plus élevées, en meilleur accord avec les observations. Il faut cependant prendre ces résultats plus comme une indication que comme une validation (faible nombre de jours concernés, non prise en compte des montagnes dans la prévision du transport vertical).

La Fig. 2.41 montre également les profils de vapeur d'eau à Trappes pour les 8 et 9 août. On retrouve des résultats similaires à ceux des simulations unidimensionnelles. L'introduction du modèle du thermique permet de prédire des couches limites plus étendues, en meilleur accord avec les sondages observés. On retrouve aussi le fait que l'effet des thermiques est relativement similaire pour les paramétrisations MY et LMD. La paramétrisation MY+TH se comporte globalement très bien.

Il est à noter que les résultats sont obtenus ici avec le modèle tridimensionnel guidé. Dans les simulations unidimensionnelles montrées précédemment, il avait fallu une phase d'ajustement pour prédire les forçages à grande échelle du modèle, essentiels pour obtenir un bon accord sur les profils verticaux simulés. Ici, les simulations sont évidemment beaucoup plus lourdes d'un point de vue informatique, mais le forçage grande échelle est calculé automatiquement par le modèle zoomé et guidé. Cette approche permet donc facilement d'effectuer des simulations contraintes sur de longues périodes de temps.

**Figure 2.42:** Influence de la paramétrisation de la couche de surface. Cycles diurnes du radon, de l'humidité et de la température (période du 7 au 10 août 1998) pour la simulation nominale (MYTH), la simulation LMD et une simulation avec le même modèle mais avec un seuil minimum sur la diffusivité verticale introduit pour éviter les inversions trop fortes dans les régions polaires l'hiver.
$\includegraphics[angle=-90,width=8cm]{therm/FIGURES/trac01KS.eps} \includegraphics[angle=-90,width=8cm]{therm/FIGURES/tempKS.eps}$ $\includegraphics[angle=-90,width=8cm]{therm/FIGURES/rhumKS.eps} \includegraphics[angle=-90,width=8cm]{therm/FIGURES/ovapKS.eps}$

Pour finir, notons que si le cycle diurne des concentrations en surface semble relativement peu sensible à la paramétrisation des thermique, ces mesures peuvent cependant permettre de relever quelques grosses déficiences des paramétrisations. On montre pour illustration sur la Fig. 2.42 l'impact sur le cycle diurne de l'introduction d'un seuil minimum sur la diffusivité verticale quand on utilise le modèle du LMD. L'introduction de ce seuil dégrade considérablement la simulation du cycle diurne du radon. Ce seuil est effectivement activé dans la configuration nominale du modèle de climat afin d'éviter que ne se crée un découplage irréaliste entre la première couche d'atmosphère et la surface, notamment l'hiver dans les hautes latitudes. L'introduction de ce seuil, grandement bénéfique pour le climat des hautes latitudes, dégrade considérablement la simulation du cycle diurne du radon.

Conclusions

Résumé des principaux résultats

On a donc introduit dans le modèle LMDZ une nouvelle paramétrisation qui concerne une échelle intermédiaire entre les échelles ``diffuses" de la couche limite et les échelles de la convection nuageuse. Cette paramétrisation est inspirée des schémas en flux de masse de la convection nuageuse mais diffère fortement de ceux-ci par la fermeture employée, appuyée ici sur une analyse géométrique des structures méso-échelles organisées de la couche limite convective. Cette paramétrisation présente également de fortes analogies avec le schéma en convection asymétriques de Pleim et Chang (1992). Le flux dans l'ascendance est cependant spécifié très différemment. Dans le modèle de Pleim et Chang (1992), ce flux est inspirée des approches en similitude avec un profil de vitesse calculé à partir de la hauteur de la couche limite et des flux en surface. Ce modèle suppose également que l'ascendance est alimentée uniquement par la première couche du modèle. Comme le modèle en convection asymétrique, le modèle du thermique permet de rendre compte de façon relativement physique du transport de chaleur en remontant le gradient dans les couches limites convectives.

Avec la nouvelle paramétrisation, on arrive à reproduire assez fidèlement les résultats des simulations des grands tourbillons pour des cas de couche limite convective claire. Chose importante, le modèle permet de bien rendre compte de la sensibilité du transport à l'intensité relative des forçages thermique et mécanique de la turbulence.

Le modèle du thermique permet de prédire un certain nombre de variables comme la largeur des thermiques, les fluctuations turbulentes du vent ou de la température dans la couche mélangée. Les thermiques négligent complètement la composante de turbulence de petite échelle dans la couche limite. Cette hypothèse pourrait être levée en utilisant le modèle du thermique pour advecter l'énergie cinétique turbulente depuis la couche de surface, en prenant éventuellement en compte la génération de turbulence par cisaillement sur le bord du panache ascendant. Un tel raffinement ne sera envisagé que si il répond à des problèmes particuliers.

Certaines améliorations ont été entreprises ou sont envisagées sur le modèle. Le modèle d'origine commençait par calculer les caractéristiques d'un thermique pour chaque couche instable pour ensuite les regrouper en un thermique unique. En fait, on peut en ne changeant que peu les résultats ne calculer qu'un thermique depuis le début, en se donnant a priori le profil vertical de l'entrainement vers le thermique (ce profil peut par exemple être pris comme fonction de et de $\partial \theta/\partial z$ ). On allège ainsi les calculs de façon significative. Les résultats de la Section 2.7 ont d'ailleurs été obtenus avec cette version modifiée du schéma (le travail de modification a été effectué par Catherine Rio).

On pourrait également penser à rajouter un excès de température dans cette couche d'alimentation, tenant compte de la dispersion des températures prédite par exemple par la fermeture en diffusion. Le calcul de cet excès de température pourrait aussi inclure des aspects liés aux hétérogénéités de surface (albédo, îlots de chaleur). On pourrait également entraîner de l'air par mélange turbulent dans la couche mélangée, ce qui aurait tendance à faire décroître la flottabilité du thermique. Les comparaisons favorables avec les simulations des grands tourbillons suggèrent soit qu'une tel calcul ne serait pas opportun soit que les résultats seraient peu sensibles à une telle sophistication.

Les résultats sont finalement assez frustrants en ce qui concerne le transport des traceurs. Les simulations donnent des résultats relativement différents, mais le nombre de degrés de liberté et le faible nombre de données (en particulier en altitude) sont tels qu'ils ne suffisent pas à départager les résultats. La vapeur d'eau, malgré sa source plus complexe, est peut-être finalement un meilleur traceur que le radon du fait du grand nombre de données disponibles. A noter également l'importance de la validation simultanée de la météorologie en surface quand on s'intéresse au transport vertical dans la couche limite. Dans l'avenir, la combinaison des profils aérosols du SIRTA et des mesures de vapeur d'eau et éventuellement de CO obtenus au décollage et à l'atterrissage par un certain nombre d'avions de lignes dans le cadre du programme MOZAIC, pourrait permettre d'avancer sur ce point.

On peut également penser dans l'avenir à utiliser de façon beaucoup plus systématique les traceurs pour valider le comportement des paramétrisations par rapport aux simulations des grands tourbillons. Un tel travail est actuellement entrepris en collaboration entre le LMD et le CNRM. Les collaborateurs du CNRM devraient en particulier refaire tourner les cas de Ayotte et al. (1996) avec des traceurs émis dans chaque couche (quelque chose d'analogue au travail de Ebert et al., 1989, sur les matrices de transiliences). C'est aussi une approche qu'on compte promouvoir dans le cadre du projet AMMA d'étude de la mousson africaine.

Pour les simulations des grands tourbillons, il serait également particulièrement intéressant d'effectuer des analyses en composites d'évènements chauds similaires à celles réalisées par Williams et Hacker (1992) pour les observations avions.

Les nuages

L'étape suivante du travail sur la paramétrisation concerne les nuages. Les rétroactions nuageuses sont à l'heure actuelle une des plus grandes sources d'incertitude des modèles de climat. On sait aussi que ces modèles de climat soufrent encore souvent de gros défauts concernant la représentation des nuages.

Il est couramment admis par exemple que les modèles de climat on tendance à sous-estimer la couverture nuageuse basse et moyenne (même si les jeux de données sont souvent insuffisant pour quantifier une telle allégation). Une intercomparaison récente de modèles basée sur l'utilisation du simulateur de radiances ISCCP^2.9, à la quelle le modèle LMDZ4 a participé (Zhang et al., 2004), confirme en tous cas la grande dispersion des modèles.

Le modèle du thermique est a priori particulièrement adapté à la simulation des cumulus de couche limite. Un travail dans cette direction a déjà commencé. On a notamment essayé de simuler les cas de cumulus montrés sur la Fig. 2.7. Plusieurs sophistications du modèle peuvent être envisagées pour ce faire. D'abord, on peut, comme pour la convection profonde, coupler le schéma de nuages au modèle du thermique. Comme on l'a rapidement expliqué dans la présentation de la partie physique de LMDZ, la fraction nuageuse, , et le contenu en eau condensé, ${q_c}$ , du nuage étaient précédemment prédits dans le modèle à partir de l'eau totale dans la maille, $\overline{q}$ , et de l'eau à saturation, ${q_{\mbox{sat}}}$ (calculée à partir de la température à grande échelle), en utilisant une distribution de probabilité de forme imposée pour la distribution sous-maille de l'eau totale. Avec ces notations, la fraction nuageuse (fraction de la maille où l'eau totale dépasse la saturation) est donnée par

$\begin{displaymath} f=\int_{{q_{\mbox{sat}}}}^{\infty} P(q)dq \end{displaymath}$

(2.91)

Le contenu en eau condensée de la maille vaut pour sa part

$\begin{displaymath} {q_c}=\int_{{q_{\mbox{sat}}}}^{\infty}(q- {q_{\mbox{sat}}}) P(q)dq \end{displaymath}$

(2.92)

Dans la formulation originale de Le Treut et Li (1991), la distribution sous maille d'eau totale

est prescrite à parti d'une distribution carrée de largeur $\sigma=r\overline{q}$ autour de $\overline{q}$ où

est un paramètre imposé.

Au cours du réglage du modèle de climat, et après l'introduction du modèle d'Emanuel pour la convection profonde, on a modifié ce schéma en suivant le travail de Bony et Emanuel (2001). D'une part, on remplace les distributions carrées par des fonctions log-normales généralisées bornées en zéro. Ainsi bornées, ces fonctions présentent une asymétrie vers les valeurs positives, et cette asymétrie est d'autant plus grande que le rapport $r=\sigma/\overline{q}$ (où $\sigma$ est comme précédemment la largeur de la distribution) est grand. Des grandes valeurs de correspondent donc à la fois à une grande dispersion des humidités et à une forte asymétrie positive. C'est typiquement ce qui est observé dans les cumulus qui présentent des humidités très fortes mais sur une étendue relativement faible, au milieu d'un environnement très sec. Le deuxième grand changement consiste à ne plus imposer la largeur de la distribution mais à la calculer en partant de l'eau nuageuse prédite par la paramétrisation de la convection. Ici, typiquement, on prendra pour ${q_c}$ la valeur de $\hat{q}-{q_{\mbox{sat}}}$ ou $\hat{q}$ est l'eau totale dans le thermique et on en déduira le paramètre de largeur de la distribution nuageuse en inversant l'Eq. 2.93. L'introduction de ce couplage pour le schéma de convection d'Emanuel a joué un rôle déterminant dans l'amélioration de la représentation des forçages radiatifs dans le modèle LMDZ. Les premiers tests de couplage avec le modèle du thermique sont également encourageants.

La deuxième modification possible en présence de nuages concerne la fermeture. L'idée la plus simple serait de remplacer la CAPE sèche par une CAPE humide (on tient compte du dégagement de chaleur latente dans le calcul de la flottabilité). Mais se pose alors immédiatement la question de la transition vers la convection profonde et de l'articulation avec les paramétrisations utilisées pour cette convection. En effet, la fermeture du modèle du thermique est basée sur une image de la couche limite convective très éloignée des images de la convection profonde où les processus de condensation, les descentes précipitantes ou les fronts de rafales créés devant les proches froides jouent un rôle déterminant. Des tests on commencé dans le modèle de circulation pour regarder l'articulation entre le modèle du thermique et la paramétrisation de la convection profonde. On revient sur ces aspects dans la conclusion générale du document.

Utilisation du modèle zoomé et guidé

**Figure 2.43:** Evolution de la température près de la surface à Trappes pour des simulations guidées et zoomées avec différentes versions du modèle. Les simulations sont comparées à deux jeux de données : les températures extraites des réanalyses ERA40 (et des analyses opérationnelles pour 2003 et 2004) et températures à 2 m à Trappes.
$\includegraphics[angle=-90,width=13cm]{therm/FIGURES/OLIVIA/C/T0b.eps}$

Pour clore ce chapitre, soulignons le potentiel qu'offre pour le travail sur les paramétrisations la version guidée et zoomée du modèle. Pour illustrer davantage ce point, on montre des simulations réalisées par Coindreau et al. (2006) avec la configuration de grille présentée plus haut et avec différentes versions du modèle. Ces simulations étaient réalisées pour valider et ajuster une version régionale de LMDZ adaptée à l'étude des paramétrisations au SIRTA et à la surveillance de l'environnement.

Les simulations sont effectuées sur une période de 6 ans entre 1998 et 2004. Pour le vent, l'humidité relative et la température, les constantes de temps de rappel sont fixées à 30 minutes à l'extérieur du domaine zoomé et 10 jours à l'intérieur afin de laisser le maximum de degrés de liberté au modèle dans la région d'intérêt et éviter les dérives à l'extérieur.^2.10On montre tout d'abord sur la Fig. 2.43 l'évolution simulée de la température à 2 m à Trappes pour les années 2000 à 2004. On voit que le modèle surestime de façon significative le cycle saisonnier, avec des températures trop chaudes de 5 degrés environ l'été. Chose intéressante à noter, ce biais est tout à fait similaire à celui observé dans le modèle climatique standard qui surestime largement les températures estivales sur les continents de l'hémisphère nord. On retrouve ce biais chaud estival pour toutes les simulations montrées ici, en utilisant la couche limite originale du LMD (LMD) ou le modèle du thermique couplé à MY (MY+TH), en utilisant le schéma de Tiedtke (TI) ou d'Emanuel (KE) pour la convection profonde, ou encore en utilisant le modèle de sol en saut d'eau (BUCKET) ou le modèle ORCHIDEE. On voit cependant que l'utilisation du modèle BUCKET aboutit à un cycle saisonnier moins en phase.

**Figure 2.44:** Cycle saisonnier moyen de la température moyenne à 2 m, de l'amplitude du cycle diurne de la même température (en fait l'écart-type des températures sur la journée) et de l'humidité relative. Le dernier graphique montre l'évolution du contenu en eau du sol sur les 5 dernières années de simulation.
$\includegraphics[angle=-90,width=8cm]{therm/FIGURES/OLIVIA/C/T3b.eps} \includegraphics[angle=-90,width=8cm]{therm/FIGURES/OLIVIA/C/T4b.eps}$ $\includegraphics[angle=-90,width=8cm]{therm/FIGURES/OLIVIA/C/q3b.eps} \includegraphics[angle=-90,width=8cm]{therm/FIGURES/OLIVIA/C/qsol0b.eps}$

Ce comportement est bien visible si on regarde le cycle saisonnier moyen de la température (en haut à gauche de la Fig. 2.44). Ce comportement est fortement corrélé à celui de l'humidité relative. On voit en effet que la surestimation des températures estivales est associée à une trop grande sécheresse. Le comportement spécifique du modèle en saut d'eau s'explique aussi à partir de l'humidité. En hiver et au début du printemps, quand le bilan précipitation moins évaporation est positif, le sol est gorgé d'eau et le modèle surestime l'évaporation (égale à l'évaporation potentielle dans ces conditions). Les températures sont donc trop froides. Quand le bilan devient négatif, le sol se vide rapidement et on retrouve finalement en août un air aussi sec et aussi chaud que dans les simulations avec ORCHIDEE. On voit en fait que le contenu en eau du sol est au plus de l'ordre de quelques millimètres l'été (en bas à droite de la figure) alors qu'on a vu dans les tests plus haut qu'il fallait conserver des humidités de l'ordre de 10 mm en août avec la formule utilisée dans le modèle BUCKET.

Le phasage saisonnier est meilleur quand on utilise ORCHIDEE. Cependant, le modèle n'arrive pas à suffisamment évaporer l'été, malgré un contenu en eau qui est loin de s'annuler. Le modèle ORCHIDEE contient en fait deux réservoirs d'eau. Un réservoir superficiel et un réservoir profond. C'est le contenu de ce dernier qui est montré sur la figure. Il semble que la capacité à évaporer l'eau stockée dans ce réservoir profond soit très insuffisante dans le modèle.

On voit que le modèle a suffisamment de liberté pour véritablement évaluer les paramétrisations du modèle de climat. A noter que des rétroactions complexes du modèle complet de climat sont ici mises en jeu. Par exemple, la surestimation de l'humidité en hiver avec le BUCKET aboutit à des pluies plus importantes qui renforcent en retour l'humidité du sol. On voit aussi que le modèle, même biaisé, permet de retrouver des éléments de la variabilité inter-annuelle comme la canicule de 2003. On peut donc envisager de décortiquer certains des mécanismes mis en jeu avec cet outil.

En attendant de disposer d'un modèle de surface plus performant (le modèle à dix couches développé par Patricia De Rosnay devrait être testé prochainement), on compte imposer un cycle saisonnier de l'humidité du sol pour éviter les gros biais saisonniers montrés ici. Cette étape est essentielle si l'on veut travailler plus finement sur la dynamique de la couche limite et des nuages associés, en tirant avantage de l'instrumentation sophistiquée installée au SIRTA (mesures météo classiques, mat, radiomètres, lidars, radars).

Bibliographie

Abdella, K., et N. McFarlane, 1997: A new second-order turbulence closure scheme for the planetary boundary layers, J. Atmos. Sci., 54, 1850-1867.
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.
Angelats i Coll, M., F. Forget, M. A. López-Valverde, P. L. Read, et S. R. Lewis, 2004: Upper atmosphere of Mars up to 120 km: Mars Global Surveyor accelerometer data analysis with the LMD general circulation model, J. Geophys. Res., 1011.
Arakawa, R. A., et W. H. Schubert, 1974: Interaction of a cumulus cloud ensemble with the large scale environment. part I, J. Atmos. Sci., 31, 674-701.
Atkinson, B. W., et J. W. Zhang, 1996: Mesoscale shallow convection in the atmosphere, Reviews of Geophysics, 34, 403-431.
Ayotte, K. W., P. P. Sullivan, A. Andrén, S. C. Doney, A. A. Holtslag, W. G. Large, J. C. McWilliams, C.-N. Moeng, M. J. Otte, J. J. Tribbia, et J. C. Wyngaard, 1996: An evaluation of neutral and convective planetary boundary-layer parameterizations relative to large eddy simulations, Boundary-layer Meteorol., 79, 131-175.
Betts, A. K., 1973: Non-precipitating convection and its parameterization, Q. J. R. Meteorol. Soc., 99, 178-196.
Blackadar, A. K., 1962: The vertical distribution of wind and turbulent exchange in neutral atmosphere, J. Geophys. Res., 67, 3095-3102.
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.
Bopp, L., O. Boucher, O. Aumont, S. Belviso, J.-L. Dufresne, M. Pham, et P. Monfray, 2004: Will marine dimethylsulphide emissions alleviate global warming? a model study, Canadian Journal of Fisheries and Aquatic Sciences, 61(5), 826-835.
Boucher, O., C. Moulin, S. Belviso, O. Aumont, L. Bopp, E. Cosme, R. von Kuhlmann, M. Lawrence, M. Pham, M. S. Reddy, J. Sciare, et C. Venkataraman, 2003: Sensitivity study of the dms atmospheric concentrations and the sulfate aerosol indirect radiative forcing to the DMS source representation, Atmos. Chem. Phys., 3, 49-65.
Boucher, O., et M. Pham, 2002: History of sulfate aerosol radiative forcings, Geophys. Res. Lett., 29, 22-1.
Bougeault, P., et P. Lacarrere, 1989: Parameterization of orography induced turbulence in a mesobeta-scale model, Mon. Wea. Rev., 117, 1872-1890.
Brost, R. A., et J. C. Wyngaard, 1978: A model study of the stably stratified planetary boundary layer, J. Atmos. Sci., 1427-1440.
Brown, R. A., 1972: On the Inflection Point Instability of a Stratified Ekman Boundary Layer., J. Atmos. Sci., 29, 850-859.
Brown, R. A., 1980: Longitudinal instabilities and secondary flows in the planetary boundary layer: A review., Rev. Geophys. Space Phys., 18, 638-697.
Businger, J. A., J. C. Wyngaard, Y. Izumi, et E. F. Bradley, 1971: Flux-pofiles relationships in the atmospheric surface layer, J. Atmos. Sci., 181-189.
Cabane, M., E. Chassefière, et G. Israel, 1992: Formation and growth of photochemical aerosols in Titan's atmosphere, Icarus, 96, 176-189.
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.
Coindreau, O., F. Hourdin, M. Haeffelin, A. Mathieu, et C. Rio, A global climate model with strechable grid and nudging: a tool for assesment of physical parametrizations, Submitted to Monthly Weather Rev., 2006.
Colella, P., et P. R. Woodward, 1984: The piecewise parabolic method (PPM) for gaz-dynamical simulations, J. Computational Phys., 54, 174-201.
Cosme, E., C. Genthon, P. Martinerie, O. Boucher, et M. Pham, 2002: The sulfur cycle at high-southern latitudes in the LMD-ZT General Circulation Model, J. Geophys. Res., 107(D23), 4690.
Couvreux, F., F. Guichard, J. L. Redelsperger, C. Kiemle, V. Masson, J. P. Lafore, et C. Flamant, 2005: Water-vapour variability within a convective boundary-layer assessed by large-eddy simulations and IHOP_2002 observations, Q. J. R. Meteorol. Soc., 131, 2665-2693.
Cressman, G. P., 1959: An operational objective analysis system, Mon. Wea. Rev., 87, 367-374.
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.
Deardorff, J. W., 1966: The counter-gradient heat-flux in the lower atmosphere and in the laboratory, J. Atmos. Sci., 23, 503-506.
Deardorff, J. W., 1970: Convective velocity and temperature scales for the unstable planetary boundary layer and for Rayleigh convection, J. Atmos. Sci., 27, 1211-1213.
Deardorff, J. W., 1972a: Parameterization of the boundary layer for use in general circulation models, Mon. Wea. Rev., 77, 93-106.
Deardorff, J. W., 1972b: Theoretical expression for the countergradient vertical heat flux, J. Geophys. Res., 77, 5900-5904.
Dufresne, J.-L., L. Fairhead, H. Le Treut, M. Berthelot, L. Bopp, P. Ciais, P. Friedlingstein, et P. Monfray, 2002: On the magnitude of positive feedback between future climate change and the carbon cycle, Geophys. Res. Lett., 29, 43-1.
Ebert, E. E., U. Schumann, et R. B. Stull, 1989: Non local turbulent mixing in the convective boundary layer evaluated from Large-Eddy Simulation, J. Atmos. Sci., 46, 2178-2207.
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., 1991: A scheme for representing cumulus convection in large-scale models, J. Atmos. Sci., 48, 2313-2335.
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.
Forget, F., F. Hourdin, et O. Talagrand, 1998: CO snow fall on Mars: Simulation with a general circulation model, Icarus, 131, 302-316.
Forget, F., F. Hourdin, R. Fournier, C. Hourdin, O. Talagrand, M. Collins, S. R. Lewis, P. L. Read, et J.-P. Huot., 1999: Improved general circulation models of the Martian atmosphere from the surface to above 80 km, J. Geophys. Res., 104, 24,155-24,176.
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 cloud-resolving models and single-column models, Q. J. R. Meteorol. Soc., 130, 3139-3172.
Hauglustaine, D. A., F. Hourdin, L. Jourdain, M.-A. Filiberti, S. Walters, J.-F. Lamarque, et E. A. Holland, 2004: Interactive chemistry in the Laboratoire de Météorologie Dynamique general circulation model: Description and background tropospheric chemistry evaluation, J. Geophys. Res., 109, 4314-4357.
Holtslag, A. A. M., E. DeBruijn, et H. Pan, 1990: A high resolution air mass transformation model for short-range weather forecasting, Mon. Wea. Rev., 118, 1561-1575.
Holtslag, A. A. M., et B. A. Boville, 1993: Local versus non-local boundary-layer diffusion in a global climate model, J. Climate, 6, 1825-1842.
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., R. Corutin, D. Gautier, C. P. McKay, P. Rannou, et O. Talagrand, Dynamique de l'atmosphère de titan, in Actes du 2 $^{\mbox{ème}}$ Colloque National de Planétologie de l'INSU, vol 1, edited by W. Kofman, J. Lilensten, B. Pibaret, et I. Raynaud, 1998.
Hourdin, F., F. Couvreux, et L. Menut, 2002: Parameterisation of the dry convective boundary layer based on a mass flux representation of thermals, J. Atmos. Sci., 59, 1105-1123.
Hourdin, F., S. Lebonnois, D. Luz, et P. Rannou, 2004: Titan's stratospheric composition driven by condensation and dynamics, J. Geophys. Res., 109.
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, et P. LeVan, The lmdz4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection, Clim. Dyn., in press, 2006.
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.
Högström, U., 1988: Non-dimensional wind and temperature profiles in the atmospheric surface layer : a re-evaluation, Boundary-layer Meteorol., 55-78.
Idelkadi, A., Validation du transport atmosphérique direct et inverse des espèces traces aux échelles continentales dans un modèle de circulation générale atmosphérique à grille variable, thèse, Université Paris-6, 2002.
Jacob, D. J., M. J. Prather, P. J. Rasch, R. Shia, Y. J. Balkanski, S. R. Beagley, D. J. Bergmann, W. T. Blackshear, M. Brown, M. Chiba, M. P. Chipperfield, J. de Grandpré, J. E. Dignon, J. Feichter, C. Genthon, W. L. Grose, P. S. Kasibhatla, I. Köhler, M. A. Kritz, K. Law, J. E. Penner, M. Ramonet, C. E. Reeves, D. A. Rotman, D. Z. Stockwell, P. F. J. Van Velthoven, G. Verver, O. Wild, H. Yang, et P. Zimmermann, 1997: Evaluation and intercomparison of global atmospheric transport models using $^{222}$ Rn and other short-lived tracers, J. Geophys. Res., 102, 5953-5970.
Jacob, D., et M. Prather, 1990: Radon 222 as a test of convective transport in a general circulation model, Tellus, 42, 118-134.
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.
Krinner, G., et C. Genthon, 2003: Tropospheric transport of continental tracers towards Antarctica under varying climatic conditions, Tellus, 53, 54-70.
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.
Lebonnois, S., D. Toublanc, F. Hourdin, et P. Rannou, 2001: Seasonal variations of titan's atmospheric composition, Icarus, 152, 384-406.
Lebonnois, S., F. Hourdin, P. Rannou, D. Luz, et D. Toublanc, 2003: Impact of the seasonal variations of composition on the temperature field of Titan's stratosphere, Icarus, 163, 164-174.
Lefèvre, F., S. Lebonnois, F. Montmessin, et F. Forget, 2004: Three-dimensional modeling of ozone on Mars, J. Geophys. Res. (planets), 109, 7004-+.
LeMone, M. A., 1973: The structure and dynamics of horizontal roll vorticies in the planetary boundary layer, J. Atmos. Sci., 30, 1077-1091.
Lilly, D. K., 1968: Models of cloud-topped mixed layers under a strong inversion, Q. J. R. Meteorol. Soc., 94, 292-309.
Lin, S.-J., et R. B. Rood, 1996: Multi-dimensional flux form semi-lagrangian transport schemes, Mon. Wea. Rev., 124, 2046-2070.
Mahowald, N. M., P. J. Rasch, B. E. Eaton, S. Whittlestone, et R. G. Prinn, 1997: Transport of $^{222}$ radon to the remote troposphere using the Model of Atmospheric Transport and Chemistry and assimilated winds from ECMWF and the National Center for Environmental Prediction/NCAR, J. Geophys. Res., 102, 28,139-28,151.
Mellor, G. L., et T. Yamada, 1974: A hierarchy of turbulence closure models for planetary boundary layers, J. Atmos. Sci., 31, 1791-1806.
Mellor, G. L., 1973: Analytical prediction of the properties of stratified planetary surface layers, J. Atmos. Sci., 30, 1061-1069.
Menut, L., R. Vautard, M. Beekmann, et C. Honoré, 2000: Sensitivity of photochemical pollution using the adjoint of a simplified chemistry-transport model, J. Geophys. Res., 105, 15379-15402.
Moeng, C., et P. P. Sullivan, 1994: A comparison of shear- and buoyancy-driven planetary boundary layer flows, J. Atmos. Sci., 51, 999-1022.
Moeng, C., 1984: A large-eddy-simulation model for the study of planetary boundary-layer turbulence, J. Atmos. Sci., 41, 2052-2062.
Montmessin, F., F. Forget, P. Rannou, M. Cabane, et R. M. Haberle, 2004: Origin and role of water ice clouds in the Martian water cycle as inferred from a general circulation model, J. Geophys. Res., 10004-+.
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.
Pierrehumbert, R. T., et R. Roca, 1998: Evidence for control of Atlantic subtropical humidity by large-scale advection, Geophys. Res. Lett., 25, 4537-4540.
Pleim, J. E., et J. S. Chang, 1992: A non-local closure model for vertical mixing in the convective boundary layer, Atmosph. Environ., 26A, 965-981.
Prandtl, L., et O. G. Tietjens, Applied hydro- and aeromechanics, Engineering Societies Monographs, Dover Publications, Inc., New-York, 1934.
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.
Quaas, J., O. Boucher, et F.-M. Bréon, 2004: Aerosol indirect effects in polder satellite data and the laboratoire de météorologie dynamique-zoom (lmdz) general circulation model, J. Geophys. Res., 109, D08205, doi:10.1029/2003JD004317.
Rannou, P., M. Cabane, E. Chassefière, R. Botet, C. P. McKay, et R. Courtin, 1995: Titan's geometric albedo : Role of the fractal structure of the aerosols, Icarus, 96, 355-372.
Rannou, P., F. Hourdin, et C. P. McKay, 2002: A wind origin for Titan's haze structure, Nature, 418, 853-856.
Rannou, P., F. Hourdin, C. P. McKay, et D. Luz, 2004: A coupled dynamics-microphysics model of Titan's atmosphere, Icarus, 170, 443-462.
Rannou, P., S. Lebonnois, F. Hourdin, et D. Luz, 2005: Titan atmosphere database, Adv. Space Res., 36, 2194-2198.
Reddy, M. S., O. Boucher, C. Venkataraman, S. Verma, J.-F. Léon, et a. M. P. N. Bellouin, 2004: Gcm estimates of aerosol transport and radiative forcing during indoex, J. Geophys. Res., 109(D16), D16205, doi:10.1029/2004JD004557.
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.
Sadourny, R., 1975a: Compressible model flows on the sphere, J. Atmos. Sci., 32, 2103-2110.
Sadourny, R., 1975b: The dynamics of finite-difference models of the shallow-water equations, J. Atmos. Sci., 32, 680-689.
Schumann, U., et C. Moeng, 1991: Plume fluxes in clear and cloudy boundary layers, J. Atmos. Sci., 48, 1746-1770.
Sommeria, G., et M. A. LeMone, 1978: Direct testing of a three dimensional model of the planetary boundary layer against experimental data, J. Atmos. Sci., 35, 25-39.
Stull, R. B., 1984: Transilient turbulence theory. Part I: The concept of eddy-mixing across finite distances, J. Atmos. Sci., 41, 3351-3367.
Stull, R. B., An introduction to boundary layer meteorology, Kluwer Academic Publishers, 1988.
Tiedtke, M., 1989: A comprehensive mass flux scheme for cumulus parameterization in large-scale models, Mon. Wea. Rev., 117, 1179-1800.
Tokano, T., F. M. Neubauer, M. Laube, et C. P. McKay, 2001: Three-Dimensional Modeling of the Tropospheric Methane Cycle on Titan, Icarus, 153, 130-147.
Toublanc, D., J. P. Parisot, J. Brillet, D. Gautier, F. Raulin, et C. P. McKay, 1995: Photochemical modelling of Titan's atmosphere, Icarus, 113, 2-26.
Troen, I., et L. Mahrt, 1986: A simple model of the atmospheric boundary layer: Sensitivity to surface evaporation, Boundary-layer Meteorol., 37, 129-148.
H. Van Dop et K. Nodop (Eds.), ETEX, A european tracer experiment, vol. 32 of Atmos. Environn. special issue, pp. 4089-4378, 1998.
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.
Wang, S., et B. Stevens, 1996: Top-Hat Representation of Turbulence Statistics in Cloud-Topped Boundary Layers: A Large Eddy Simulation Study., Nature, 382, 528-531.
Weckwerth, T. M., J. W. Wilson, R. M. Wakimoto, et N. A. Crook, 1997: Horizontal convective rolls: determining the environmental conditions supporting their existence and characteristics, Mon. Wea. Rev., 125, 505-526.
Williams, A. G., et J. M. Hacker, 1992: The composite shape and structure of coherent eddies in the convective boundary layer, Boundary-layer Meteorol., 61, 213-245.
Williams, A. G., et J. M. Hacker, 1993: Interactions between coherent eddies in the lower convective boundary layer, Boundary-layer Meteorol., 64, 55-74.
Williams, A. G., Internal structure and interactions of coherent eddies in the lower Convective Boundary Layer, Ph.D. thesis, Flinders University of South Australia, 1991.
Williamson, D., et P. J. Rasch, 1989: Two dimensional semi-lagrangian transport with shape-preserving 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.
Yamada, T., 1983: Simulations of nocturnal drainage flows by a turbulence closure model, J. Atmos. Sci., 40, 91-106.
Zhang, M. H., W. Y. Lin, S. A. Klein, J. T. Bacmeister, S. Bony, R. T. Cderwall, A. D. Del Genio, J. J. Hack, N. G. Loeb, U. Lohmann, P. Minnis, I. Musat, R. Pincus, P. Stier, M. J. Suarez, M. J. Webb, J. B. Wu, M.-S. Xie, et J. H. Zhang, Comparing clouds and their seasonal variations in 10 atmospheric general circulation models with satellite measurements, submitted to JGR, 2004.

À 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.70)

The command line arguments were:
latex2html -split 1 hdrt

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

HOURDIN Frédéric 2006-06-15

M. Michel CABANE	Président
M. Pierre DROSSART	Rapporteur
M. Ray PIERREHUMBERT	Rapporteur
M. Patrick MASCART	Rapporteur
M. Alain HAUCHECORNE	Examinateur
M. Daniel GAUTIER	Examinateur
M. Olivier TALAGRAND	Examinateur