\section{Le modèle du thermique\label{sec:concept}}

\def\fm{f}

\subsection{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{cell}) instable
(typiquement de une à quelques centaines de mètres d'épaisseur),
une couche mélangée (CM sur la \fig{cell}) neutre, surmontée d'une zone stable (couche
d'entraînement plus troposphère libre). 
\begin{figure*}
\centerline{\includegraphics[width=14cm]{\local/FIGURES/cell.eps}}
\caption{Vue transverse idéalisée d'un rouleau convectif à la base du modèle
du thermique.
\label{fg:cell}}
\end{figure*}

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, $\tsbl$.\footnote{
Comme dans les sections précédentes, et afin d'alléger les notations,
on note $\theta$ la température potentielle virtuelle.}

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{equation}
\frac{d w}{d t}=w\frac{\partial w}{\partial z}= g \frac{\tsbl-\tmbl}{\tmbl}
\label{eq:dw}
\end{equation}
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 $\hmbl$.
A ce niveau, le carré de la vitesse verticale $\wmax$ est égal à deux fois
l'énergie potentielle disponible pour la convection ou CAPE
(pour Convective Available Potential Energy)
\begin{equation}
\frac{1}{2}w^2_{\mbox{max}}=\cape=\int_0^{\hmbl}  g \frac{\tsbl-\tmbl}{\tmbl} dz
\end{equation}
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 $\hth$ où
\begin{equation}
\int_0^{\hth} g\frac{\tsbl-\tmean}{\tmean} dz=0
\end{equation}
Cette intégrale est représentée par la zone grisée
sur la gauche de la \fig{cell}.
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
$\fm=\fraca \rho w$ où $\fraca$ 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 $\fm$ 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 $x$. 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{equation}
v\frac{\partial v}{\partial y}= -\frac{1}{\rho_s}
\frac{\partial \ps}{\partial y}
\end{equation}
où $\ps$ 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{equation}
\vmax^2\simeq \frac{2}{\rho_s}\Delta \ps
\end{equation}
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{equation}
\ps=p(z_i)+\int_0^{\hmbl} g\rho dz
\end{equation}
et donc, au premier ordre,
\begin{eqnarray}
\vmax^2&\simeq&\frac{2}{\rho_s}\Delta\int_0^{\hmbl} g \rho dz\\
&\simeq&2\int_0^{\hmbl} g\frac{\rho}{\rho_s}\frac{\tsbl-\tmbl}{\tmbl} dz
\end{eqnarray}
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
$\wmax$.

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{equation}
\wmax l(\hmbl) \rho(\hmbl)=\vmax \hsbl \rho(\hsbl/2)
\end{equation}
où $l(z)$ est la largeur du rouleau à la hauteur $z$.
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 $L$ 
entre la subsidence et l'ascendance est de l'ordre de 1.
En effet, si la largeur de la cellule $L$ était beaucoup plus grande que
sa hauteur, une particule partant de la distance $L$ 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
$\aspr=L/\hth$ est utilisé comme paramètre libre. La couverture fractionnaire
est minimum au niveau $\hmbl$ où elle vaut
$\fraca(\hmbl) =l(\hmbl)/(r \hth)$. On en tire une relation de fermeture
pour le calcul du flux de masse dans l'ascendance~:
\begin{equation}
\fm=\frac{l \wmax}{r \hth}=\frac{\rho(\hsbl/2) \hsbl \sqrt{2\cape}}{r \hth}
\label{eq:fm}
\end{equation}
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 $2-3$ \cite[]{LeMo:73,Moen:94,Atki:96} 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 à $2 r$.
La valeur $r=$1 qui sort de l'analyse ci-dessus est donc 
compatible avec les résultats des \LES.
Pour avoir méconnu ce facteur 2 dans un premier temps,
nous avons retenu $r=2$ 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).

\subsection{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ù $w$ 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 à $\hth$ (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}$
\cite[voir par exemple ][Ch.IV]{Pran:34}.
Dans notre cas,
$\lmix=\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$^2$~s$^{-1}$, $\lmix=10$~m.
On retient cette formulation pour réduire la section
du thermique en dessous de $\hmbl$, en considérant $\lmix$ comme un paramètre
ajustable.

\subsection{Les équations du modèle}

\def\e{\hat{e}}
\def\d{\hat{d}}
\def\ld{l_d}
\def\f{\hat{f}}
\def\phia{\hat{\phi}}
\def\vent{{\bf v}}
\def\va{\hat{\vent}}
\def\ve{{\vent}_e}
\def\fraca{\hat{\alpha}}
\def\wa{\hat{w}}
\def\ta{\hat{\theta}}


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 ($\wa$), sa température potentielle
($\ta$) et la couverture fractionnaire des thermiques ($\fraca$).
Le flux de masse dans les ascendances est $\f=\rho\fraca\wa$
($\rho$ est supposé identique à l'intérieur et à l'extérieur du thermique).
Notons
\begin{equation}
\mbox{I}(z,z')=\int_z^{z'} g \frac{\theta(z)-\theta(z'')}{\theta(z'')} dz''
\end{equation}
Pour de l'air instable provenant de l'altitude $z$
(si en $z$, $\partial \theta/\partial z < 0$), la
\CAPE\ est $I(z,hi(z))$ 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
$\hth$ comme la valeur maximum de $ht(z)$
avec
$ht(z)=\min\{z'\ \mbox{t.~q. }\ z'>z\ \mbox{et}\ \mbox{I}(z,z') = 0\}$.


Le flux de masse $\f$ est relié au flux d'air entrant $\e$ et au flux d'air
détraîné $\d$ par~:
\begin{equation}
\label{eq:f}
\frac{\partial \f }{\partial z}=\e  -\d
\end{equation}
En suivant l'analyse précédente (et l'\eq{fm}), on prescrit $\e$ suivant
\begin{equation}
\e(z)= {\rho(z)}\sqrt{2\mbox{CAPE}(z,h(z))}/\dep{r \hth}
\end{equation}
si $\partial \theta/\partial z<0$ et 0 sinon.
En dessous de $z_i$, le détraînement est calculé suivant
\begin{equation}
\d(z)=\frac{\partial}{\partial z} \dep{\rho \wa \sqrt{\lmix z}/(r \hth)}
\end{equation}
Au dessus de $z_i$, on réduit la largeur du thermique suivant
$[(\hth-z)/(\hth-\hmbl)]^\mu$. Comme $\e=0$ à cet endroit,  le détraînement
entre $z_i$ et $\hth$ s'écrit
\begin{equation}
\d(z)=-
\frac{\partial}{\partial z} \depb{\rho \wa \fraca(z_i)
\dep{\frac{\hth-z}{\hth-\hmbl}}^\mu}
\end{equation}
On utilise
$\mu=1$ et $\mu=2$ dans la version standard.

Le niveau de l'inversion $\hmbl$ est défini comme l'altitude où la température
à l'intérieur du thermique $\ta$ devient plus petite que celle de l'environnement
(début de l'overshoot).
\def\phim{\overline{\phi}}
On calcule alors les propriétés de l'air à l'intérieur du thermique
suivant l'équation de conservation
\begin{equation}
\frac{\partial \f \phia}{\partial z}=\e \phi_e -\d \phia\label{eq:fq} \label{eq:phi}
\end{equation}
où $\phi$ est un scalaire (y compris la température potentielle) et l'indice
$_e$ est relatif à l'environnement.
$\phi_e$ est relié à $\phia$ (dans le thermique) et à la valeur moyenne dans
la maille $\phim$ à un niveau donné
\begin{equation}
\phim=\fraca\phia+(1-\fraca)\phi_e
\end{equation}
On recalcule la vitesse verticale moyenne dans le thermique, en tenant compte
de sa flottabilité, à partir de l'équation
\begin{equation}
\frac{\partial \f \wa}{\partial z}=
-\d \wa+g\rho\frac{\ta-\theta_e}{\theta_e}\label{eq:fw}
\end{equation}
L'entraînement n'a pas de contribution à l'\eq{fw} dans la mesure où on
suppose que l'air entre dans le thermique avec une vitesse verticale
nulle.

\def\v{{\bf v}}
Pour le transport de la quantité de mouvement, l'idée la plus simple consiste
à traiter les deux composantes horizontales $u$ et $v$ 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 $L^2$ pour se ramener à une force par unité de surface
horizontale de sorte qu'on obtient finalement
\begin{equation}
\frac{\partial \f \va}{\partial z}=\e \ve -\d \va +\gamma\rho \frac{\sqrt{\fraca}}{L} ||\ve-\va|| (\ve-\va)
\end{equation}
$\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 $-\f$
\begin{equation}
\rho\overline{w'\phi'}=\f\dep{\phia-\phi_e}
\end{equation}

Une approximation supplémentaire est finalement introduite pour rendre
le schéma plus robuste. En effet,
si on intègre telle quelle l'\eq{phi}, 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{phi} prédirait
$\phia>0$  dans la couche mélangée alors qu'on a $\phim=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=\phim$). Cette approximation qui est valide
pour $\fraca\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 $r=2$ 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 \cite{Hour:02}.

\subsection{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.

\def\hmin{h_{\mbox{min}}}
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 à $\hmin=0,07/f$ où $f$ 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^*\vk z (1-z/\hmin)^2$.

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% TeX-master: "concept"
%%% TeX-master: "concept"
%%% TeX-master: "concept"
%%% TeX-master: "concept"
%%% TeX-master: "concept"
%%% End: