\section{Bases physiques des paramétrisations en diffusion\label{sec:kdiff}} 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 ($|\partial \mave{X}/\partial x| \ll |\partial \mave{X}/\partial z|$ et $|\partial^2\mave{X}/\partial x^2| \ll |\partial^2\mave{X}/\partial z^2|$ où $\mave{X}$ est la moyenne d'ensemble de $X$ pondérée par l'air définie dans la \sec{traceurs}). \subsection{La longueur de mélange} \def\q{c} 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é $\q$ dans la couche limite \begin{equation} \overline{\rho w' \q'}=-\rho K_z \der{\q}{z} \end{equation} Une façon physique d'introduire cette formulation et d'estimer le coefficient $K_z$ 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, $l$. En se plaçant pour fixer les idées dans une configuration où la concentration moyenne du traceur croît verticalement ($\partial\overline{\q}/\partial z>0$), un mouvement descendant sera associé en moyenne à une fluctuation positive de $\q$. La fluctuation $\q'$ associée à ce mouvement turbulent sera d'autant plus grande que les contrastes verticaux de $\q$ sont grands. En supposant que la particule qui descend a conservé les propriétés qu'avait l'air à une distance $l$ au-dessus, et en supposant $l$ petit devant la hauteur caractéristique des variations de $\ave{\q}$ (c'est cette dernière hypothèse qui n'est pas valide dans les cas de couches limites convectives) on peut écrire \begin{equation} \q'=l \der{\q}{z} \end{equation} De même, les mouvement turbulents ascendants seront associés en moyenne à une fluctuation négative de $\q$ de sorte qu'il y a dans les deux cas une corrélation négative entre $w'$ et $\q'$. On peut finalement écrire \begin{equation} \ave{w'\q'}=-l\ave{|w'|} \der{\q}{z} \end{equation} Cette équation peut être prise comme définition de la longueur de mélange $l$. Dans ce cadre, le coefficient $K_z$ est simplement $l|w'|$. 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{equation} |w'|\simeq|u'|\simeq l\der{u}{z} \end{equation} ce qui conduit à choisir \begin{equation} K_z=l^2||\der{\V}{z}|| \end{equation} 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=\vk 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 $\ave{u}$ (noté $u$ ci-après) et son gradient en $z$ sont du même signe près de la surface -- puisque $\overline{u}$ doit s'annuler en $z=0$ --, il vient \begin{equation} \overline{w'u'}=-K_z \der{u}{z}=-{\dep{\vk z}}^2{\left(\der{u}{z}\right)}^2 \end{equation} 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{equation} \ustar=\vk z \der{u}{z} \end{equation} où $\ustar^2=\tau/\rho$ est la vitesse de friction. A noter que les dérivations ci-dessus aboutissent également à $|w'|\simeq \ustar$ dans la couche de surface. Les mesures expérimentales de la constante de Von Karman, $\vk$, 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 $z=0$ mais à une altitude $z=z_0$, appelée longueur de rugosité, telle que \begin{equation} \overline{u}(z)=\frac{\ustar}{\vk}\ln{\frac{z}{z_0}} \end{equation} 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\ave{w'\theta'}_0$, on peut introduire une échelle des fluctuations turbulentes de température potentielle, $\theta^*=\ave{w'\theta'}_0/u^*$, reliée au gradient de température potentielle par \begin{equation} \theta^*=\vk' z \der{\theta}{z} \end{equation} ce qui aboutit également à une forme logarithmique \begin{equation} \theta-\theta_S=\frac{\theta^*}{\vk'}\ln{\frac{z}{z_0}} \end{equation} où $\theta_S$ est la température de surface. Le rapport $R=\vk/\vk'$ 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 \cite[]{Busi:71} suggèrent $R=0,74$ en conditions neutres ou stables \cite[]{Dear:72a}. La hauteur $z_0$ peut également être différente pour la quantité de mouvement et la température. Dans ce dernier cas, $z_0$ correspond plutôt à la hauteur à partir de laquelle on passe d'un transport conductif à un transport turbulent. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{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 (\sec{traceurs}), peut bien sûr être poussée plus loin. On écrit alors des équations d'évolution pour les quantités turbulentes $u'$, $w'$, $\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 \cite{Mell:74}. 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 \cite{Mell:74} 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{equation} e=\frac{1}{2}\depb{\overline{u'^2}+\overline{v'^2}+\overline{w'^2}} \end{equation} cette énergie étant ensuite utilisée pour le calcul du coefficient de mélange $K_z$. 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 $\ave{|w'|}$ 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 $\ave{|w'|}$ (respectivement une sous-estimation pour une atmosphère instable). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Sous les hypothèses de couche limite mentionnées plus haut, on peut montrer \cite[voir \eg][]{Stul:88} que l'évolution de l'énergie cinétique turbulente s'écrit sous la forme \begin{equation} \der{e}{t}=P-D -\frac{1}{\rho}\der{\ave{w'p'}}{z} -\der{\ave{w'e}}{z}-\epsilon \end{equation} Le terme \begin{equation} P=-\overline{w'u'}\der{u}{z}-\overline{w'v'}\der{v}{z} \end{equation} 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{equation} D=-\frac{g}{\theta}\overline{w'\theta'} \end{equation} 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 $w$ 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{equation} \overline{w'\phi'}=-K_\phi \der{\phi}{z} \end{equation} On écrit la diffusivité turbulente $K_\phi=lqS_\phi$ en fonction d'une longueur de mélange, $l$, 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 ($S_h$) ou la quantité de mouvement ($S_m$). \def\ca{1,96} \def\cb{0,1912} \def\cc{0,2341} \def\cd{1} \def\ce{0,2231} \def\cf{1,318} \def\ba{0,6588} \def\bb{0,1776} \def\bc{0,3221} \def\be{0,03156} \def\ce{0,2231} \def\cf{1,1318} 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 $S_m$ et $S_h$. Ces grandeurs sont fonctions du seul {\bf nombre de Richardson de flux}, $\rif=D/P$, 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 \cite{Yama:83}. Cette version a été utilisée dans une étude d'intercomparaison de fermetures turbulentes par \cite{Ayot:96}. Elle est également utilisée plus loin dans des tests numériques. Pour la quantité de mouvement, on obtient~: \begin{equation} S_m=\ca\frac{\dep{\cb-\rif}\dep{\cc-\rif}}{\dep{\cd-\rif}\dep{\ce-\rif}}, \rif< 0,16 \end{equation} et \begin{equation} S_m=0,085, \rif\ge 0,16 \end{equation} 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{equation} \omega=\cf\frac{\ce-\rif}{\cc-\rif}, \rif< 0,16 \end{equation} \begin{equation} \omega=1,12, \rif\ge 0,16 \end{equation} En remarquant que $\rif=\omega \ri$, où \def\shear{\left|\left|\der{\V}{z}\right|\right|} \begin{equation} \ri=\frac{g}{\theta}\frac{\der{\theta}{z}}{{\shear}^2} \end{equation} {\bf est le nombre de Richardson gradient}, on obtient finalement \begin{equation} \rif= \ba [\ \ \ri+\bb - \sqrt{\ri^2-\bc\ri+\be}\ \ ]\mbox{, } \ri<\ric \end{equation} \begin{equation} \rif=\rifc\mbox{, } \ri\ge\ric \end{equation} où $\rifc=0,191$ et $\ric=0,195$ sont des nombres de Richardson critiques. On obtient ainsi des expressions explicites de $S_m$ 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 $q$ s'écrit finalement \def\v{{\bf v}} \begin{equation} \label{eq:dqorig} \frac{1}{2}\frac{\partial q^2}{\partial t}= q l S_m{\shear}^2 \depb{1 - \omega \ri} - \frac{q^3}{lB_1} + \frac{\partial}{\partial z} \depb{l q S_q \frac{\partial q^2 / 2}{\partial z}} \end{equation} avec $B_1=16,6$. \cite{Yama:83} utilise $S_q=0,2$ à 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 $l$. Dans certains articles, les auteurs ont suggéré une équation pour $q^2l$ analogue à l'équation pour $q^2$ 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 \cite{Blac:62} \begin{equation} l=l_0\frac{\vk z}{\vk z+l_0} \end{equation} asymptotique à $\vk z$ près de la surface et à une longueur $l_0$ au-dessus de la couche de surface. La longueur $l_0$ 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. \cite{Yama:83} utilise par exemple pour $l_0$ une altitude moyenne pondérée par l'intensitée de la turbulence \begin{equation} l_0=0,2\ \frac{\int_0^\infty z q dz}{\int_0^\infty q dz} \end{equation} 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. \footnote{A noter que l'intégration numérique du modèle de Mellor et Yamada s'avère souvent délicate. Une intégration naïve de l'équation d'évolution de l'énergie cinétique turbulente avec un schéma temporel explicite (on calcule les termes sources et puits du membre de droite au temps $t$ qu'on ajoute à l'énergie cinétique au temps $t$ pour obtenir la nouvelle valeur à $t+\delta t$) contraint à prendre des pas de temps de quelques secondes, même avec les discrétisations grossières utilisées dans le modèle de circulation générale. Dans la version développée pour LMDZ, on contourne en partie cette difficulté en récrivant formellement l'équation d'évolution de l'énergie cinétique turbulente (sans diffusion) sous la forme \begin{equation} \frac{1}{2}\frac{\partial q^2}{\partial t}=q^3 \chi \end{equation} ou encore \begin{equation} \label{eq:dq} \frac{\partial}{\partial t}\dep{\frac{1}{q}}=- \chi \end{equation} avec \begin{equation} \chi = \frac{l S_m}{q^2} M^2\dep{1-\rif}-\frac{1}{l B_1} \end{equation} Si on suppose que $\chi$ ne varie pas au court d'un pas de temps, la solution de l'\eq{dq} est \begin{equation} q^{(t+\delta t)} = \frac{q^{(t)}}{1-\chi^{(t)} q^{(t)} \delta t} \end{equation} On retient directement cette solution quand $\chi\le 0$. En revanche, quand $\chi>0$, on utilise une forme approchée \begin{equation} q^{(t+\delta t)} = q^{(t)} \dep{1+\chi^{(t)}q^{(t)}\delta t} \end{equation} Cette formulation numérique produit des résultats numériques presque indiscernables de l'intégration temporelle explicite de l'équation d'origine mais avec des pas de temps de typiquement quelques minutes à dizaines de minutes pour les configurations classiques du modèle de circulation. La diffusion verticale de l'énergie cinétique turbulente est calculée a posteriori. }% Fin de la footnote \begin{figure} %\begin{tabular}{cc} %\multicolumn{2}{c}{{\bf a)} cycle diurne de la température}\\ %\multicolumn{2}{c}{$\includegraphics[width=7cm]{\local/FIGURES/MARS/pbl_t.eps}$} \\ %\multicolumn{2}{c}{{\bf b)} cycle diurne du vent}\\ %$\includegraphics[width=7cm,clip]{\local/FIGURES/MARS/pbl_u_ref.eps}$ & %$\includegraphics[width=7cm]{\local/FIGURES/MARS/pbl_v.eps}$\\ %\multicolumn{2}{c}{{\bf c)} cycle diurne de l'énergie cinétique turbulente}\\ %$\includegraphics[width=7cm,clip]{\local/FIGURES/MARS/pbl_q_ref.eps}$ & %$\includegraphics[width=7cm]{\local/FIGURES/MARS/pbl_q.eps}$ %\end{tabular} \begin{tabular}{c} {{\bf a)} cycle diurne de la température}\\ \includegraphics[width=7cm]{\local/FIGURES/MARS/pbl_t.eps}\\ {{\bf b)} cycle diurne du vent}\\ $\includegraphics[width=7cm,clip]{\local/FIGURES/MARS/pbl_u_ref.eps}\includegraphics[width=7cm]{\local/FIGURES/MARS/pbl_v.eps}$\\ {{\bf c)} cycle diurne de l'énergie cinétique turbulente}\\ $\includegraphics[width=7cm,clip]{\local/FIGURES/MARS/pbl_q_ref.eps}\includegraphics[width=7cm]{\local/FIGURES/MARS/pbl_q.eps}$ \end{tabular} \caption{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. {\bf 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). {\bf b :} Module du vent horizontal (m~s$^{-1}$). {\bf c~:} Fluctuations turbulentes du vent (m~s$^{-1}$). Pour b) et c), on montre, {\bf à gauche,} les résultats du modèle en fonction de l'heure locale et de l'altitude et, {\bf à 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 ({\bf 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. \label{fg:mymars}} \end{figure} 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 \cite[]{Forg:99}. Avec sa fine atmosphère de CO$_2$, 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{mymars} 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 $q$ sur les figures {\bf 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). \subsection{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{equation} q^2=l^2{\shear}^2S_m B_1\dep{1-\omega\ri} \end{equation} et donc \begin{equation} K_m=l^2\shear \sqrt{{S_m}^3 B_1\dep{1-\omega\ri}} \end{equation} On voit qu'on retombe exactement sur l'approche de Prandtl avec un terme correctif ($\sqrt{{S_m}^3 B_1\dep{1-\omega\ri}}$) rendant compte des effets de la stratification. La formulation utilisée dans le modèle original du LMD \cite[]{Lava:81} 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{equation} K_m=l\ Max\depb{ l \shear \sqrt{1- Ri/Ri_c},\sqrt{e_{min}} } \end{equation} où $Ri_c$ 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 $l=l_0 (p/p_s)^2$ avec $l=l_0=$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 $K$ et l'autre pour la dissipation $\epsilon$. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{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 $u^*$ et $z$. 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, $\vk\ustar\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 $\rif$ introduit plus haut. Dans la couche de surface, ce nombre s'écrit \begin{equation} {\ri}_f=\frac{g}{\theta}\frac{\overline{w'\theta'}}{\overline{w'u'}\der{u}{z}}\simeq\frac{g z {\overline{w'\theta'}}_0}{\vk {\ustar}^3\theta}=\frac{z}{L} \end{equation} où \begin{equation} L=\frac{\vk {\ustar}^3 \theta}{g{\overline{w'\theta'}}_0} \end{equation} 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{equation} \der{u}{z}=\frac{\ustar}{\vk z}\phi_m\dep{\frac{z}{L}} \end{equation} 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{equation} u=\frac{\ustar}{\vk}\depb{\ln\dep{\frac{z}{z_0}}-\psi_m\dep{\frac{z}{L}}} \end{equation} où \begin{equation} \psi_m\dep{\xi}=\int_{z_0/L}^{\xi}\depb{1-\phi_m\dep{\xi}}\frac{d\xi}{\xi} \simeq \int_{0 }^{\xi}\depb{1-\phi_m\dep{\xi}}\frac{d\xi}{\xi} \end{equation} La seconde forme de $\psi_m$ permet en général des calculs analytiques relativement simples et suffisamment précis. \def\thetastar{\theta_*} Comme pour la quantité de mouvement, on suppose que le gradient de température adimensionnel est relié au gradient neutre par une fonction de $z/L$ \begin{equation} \der{\theta}{z}=\frac{R \thetastar}{\vk z}\phi_h\dep{\frac{z}{L}} \end{equation} avec $\thetastar = {\overline{w'\theta'}}_0/\ustar$ qui conduit comme pour le vent à \begin{equation} \theta-\theta_s=\frac{R \thetastar}{\vk}\depb{\ln\dep{\frac{z}{z_0}}-\psi_h\dep{\frac{z}{L}}} \end{equation} 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 \cite{Dear:72} 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 \cite{Busi:71} 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 \tb{Businger}. De nombreuses versions modifiées ont été proposées depuis, utilisant parfois des données plus récentes \cite[cf. par exemple][]{Hogs:88}. \def\phimi{\dep{1-\gamma_1\xi}^{-\frac{1}{4}}} \def\phims{1 + \gamma_3\xi} \def\psimi{\ln\depb{\dep{\frac{1+x^2}{2}}{\dep{\frac{1+x}{2}}}^2}-2 \tan^{-1}x +\frac{\pi}{2}} \def\psims{-\gamma_3\xi} \def\phihi{\dep{1-\gamma_2\xi}^\frac{1}{2}} \def\phihs{1+\gamma_3\xi/R} \def\psihi{2 \ln\dep{\frac{1+x}{2}}} \def\psihs{-\gamma_3\xi/R} \def\saut{&&\\} \begin{table} \centerline{ \begin{tabular}{ccc} & instable $L<0$ & stable $L \ge 0$ \\\saut\hline\saut $\phi_m\dep{\xi}$ & $\phimi$ & $\phims$ \\\saut $\psi_m\dep{\xi}$ & $\psimi$ & $\psims$ \\ & avec $x=1/\phi_m$ & \\\saut\hline\saut $\phi_h\dep{\xi}$ & $\phihi$ & $\phihs$ \\ \saut $\psi_h\dep{\xi}$ & $\psihi$ & $\psihs$ \\ & avec $x=1/\phi_h$ & \\\saut\hline \end{tabular} } \caption{Fonctions de Businger Dyer telles qu'elles sont utilisées dans les simulations présentées plus loin.\label{tb:Businger}} \end{table} Notons que \cite{Mell:73} 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 \cite{Busi:71} 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 $h$. Nous testons par exemple plus loin une formulation analytique du coefficient de diffusion turbulente en fonction de l'altitude proposée à l'origine par \cite{Bros:78} \begin{equation} K_m=\ustar\vk z \phi_m^{-1}\depb{1-\frac{z}{h}}^p \end{equation} 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 $z=h$. \cite{Bros:78} retiennent un exposant $p=1,5$. Pour appliquer cette formule, il faut commencer par estimer la hauteur de la couche limite $h$. Une approche maintenant classique et relativement robuste consiste à utiliser un nombre de Richardson non local~: \begin{equation}\label{eq:bulkri} R_{ib}(z)=\frac{g\,(z-z_{1})}{\theta(z)}\frac{(\theta(z)-\theta(z_{1}))} {u(z)^{2}+v(z)^{2}} \end{equation} où $z_1$ 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 \cite{Troe:86}.