\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}.

