\section{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 ($|\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|$ pour toute grandeur $X$). \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 de tenter 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-tendu 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 fluctuations de $\ave{\q}$ (c'est cette hypothèse qui n'est pas valide dans les cas de couche limite convective) 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} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{L'énergie cinétique turbulente et les fermetures d'ordre 2} 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, 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 1970. 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 dizaines d'équations pour prédire les différents termes croisés. Les fermetures utilisées en pratique dans les modèles atmosphériques, comme celui 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} 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). Le modèle de \cite{Mell:74} donne une base théorique à ces considérations ainsi que des solutions "abordables" pour la paramétrisation du transport turbulent. Sous les hypothèses de couche limite mentionnées plus haut, on peut montrer \cite[vor \eg][]{Stul:84} 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 négatif dans une atmosphère stable et correspond alors à l'inhibition de la turbulence par les effets de stratification. Ce terme 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ù elles 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 $\omega$. Ces grandeurs sont fonctions du seul nombre de Richardson de flux, $\rif=D/P$, mesure de l'importance relative des forçages mécanique et thermique de la turbulence. On retient ici les valeurs numériques données par \cite{Yama:83} et utilisées dans une étude d'intercomparaison de fermetures turbulentes par \cite{Ayot:96}. 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 spécifie le nombre de Prandtl turbulent, $\omega=S_h/S_m$ sous la forme \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ù $\ri$ est le nombre de Richardson gradient \def\shear{\left|\left|\der{\V}{z}\right|\right|} \begin{equation} \ri=\frac{g}{\theta}\frac{\der{\theta}{z}}{{\shear}^2} \end{equation} 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é, le terme de transport vertical est paramétrisé en K-diffusion 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 on 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}~: \afaire{Verifier que c'est la bonne ref a Blakadar} \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. \cite{Yama:83} utilise par exemple pour $l$ 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 rajoute à 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\ge 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 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 et où la longueur $l$ est spécifiée en fonction de la pression comme \afaire{Donner la formule de $l$ pour le modèle du LMD}. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{La couche limite de surface} 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. Cette couche limite de surface est généralement définie comme la région au-dessus de la surface dans laquelle les flux turbulents de moment et de chaleur restent proches de leurs valeurs en surface (on prend souvent un critère de 10$\%$). Cette couche de surface peut être aussi définie comme la couche dans laquelle les flux turbulents dominent largement tous les autres termes. Dans cette couche, la direction du vent est presque constante et on a l'habitude, pour simplifier, de prendre l'axe $x$ selon la direction du vent. Dans cette couche aussi, on suppose que la longueur caractéristique des échanges turbulents est proportionnelle à la distance à la surface, $l=\vk z$. Dans ces conditions, le flux turbulents dans la couche limite de surface s'écrit \begin{equation} \overline{w'u'}={\dep{\vk z}}^2\left|\der{\overline{u}}{z}\right|\der{\overline{u}}{z} \end{equation} En prenant $\overline{u}$ dans le sens du vent, en remarquant que le gradient de $\overline{u}$ est du signe de $\overline{u}$ près de la surface puisque $\overline{u}$ doit s'annuler en $z=0$, et en remarquant que, dans la couche de surface, $\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=\tau/\rho$ est la vitesse de friction. L'équation précédente montre aussi qu'on a implicitement $|w'|\simeq \ustar$ dans la couche de surface. Les mesures expérimentales de la constante de Von Karman, $\vk$, supposée universelle, donne 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 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. Cette hauteur dépend donc du caractère plus ou moins accidenté de la surface et est appelée de ce fait longueur de rugosité. Elle 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. 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. Dans la couche de surface, se nombre s'écrit \begin{equation} {\ri}_f=\frac{D}{P}=\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ù $\rho C_p {\overline{w'\theta'}}_0$ est le flux de chaleur en surface et \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 \longrightarrow 0$, c'est à dire quand la turbulence est dominée par les effets mécaniques. Cette équation peut s'intégrer verticalement pour donner 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 permettant en général des calculs analytiques relativement simples et suffisamment précis. \def\thetastar{\theta_*} Pour la température, l'approche de la longueur de mélange dans la couche limite de surface neutre donne aussi une expression du flux de la forme \begin{equation} \overline{w'\theta'}\simeq-l|w'|\der{\theta}{z}\simeq-k'z\ustar\der{\theta}{z} \end{equation} Dans les conditions neutres, par analogie avec \afaire{demander à Jean-Louis}, on a $k'\simeq 0.8 k$ \afaire{à vérifier}. Comme pour la quantité de mouvement, on suppose dans le cas général que le gradient de température potentielle s'écrit \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$ et $\phi_h(0)/\phi_m(0)\simeq 0.74 ou 1.35$. qui conduit comme pour le vent à \begin{equation} \theta=\frac{R \thetastar}{\vk}\depb{\ln\dep{\frac{z}{z_0}}-\psi_h\dep{\frac{z}{L}}} \end{equation} L'introduction du nombre de Prandtl $R$ 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. L'altitude $z_0$ jusqu'à laquelle il faut intégrer l'équation n'est pas forcément physiquement la même que celle impliquée dans les équations pour la quantité de mouvement. En effet, pour la quantité de mouvement il s'agit de la hauteur sous laquelle le transfert de moment cinétique s'effectue par couple de pression sur les obstacles alors que, pour la température, il s'agit plutôt de la hauteur sous laquelle les échanges sont dominés par la conduction. Des campagnes de mesures ont été dédiées à la mesure de ces fonctions. Les fits proposés par \cite{Busi:71} pour une campagne dans le Kansas sont encore largement utilisés. Une exemple de fit issu des travaux de \cite{Busi:71} basés sur une campagne menée dans une région très plate du Kansas est donné dans la \tb{Businger}. De nombreuses versions modifiées ont été proposée 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 Bussinger Dyer telles qu'elles sont utilisées dans les simulations présentées plus loin.\label{tb:Businger}} \end{table} Pour conclure avec cette partie sur la couche limite de surface, notons que \cite{Mell:73} a appliqué sont modèle de fermeture à la couche limite de surface, en supposant des flux constants, et réussi à interpréter les mesure 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 résultats de la couche limite de surface à la couche mélangée. Il faut alors introduire au moins un paramètre supplémentaire~: la hauteur de la couche limite $h$. Mais, si les hypothèses de flux constant (variant peu devant le cisaillement de vent), permet de dériver des expressions théoriques pour la couche limite de surface, la situation est beaucoup moins claire pour le reste de la couche limite. Nous testons cependant 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$. Une telle formulation requiert évidemment d'estimer la hauteur de la couche limite $h$. Une approche classique 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}.