Formation de glace et température
Sommaire
Calcul de l'effet de la formation de glace sur la température de la particule soulevée adiabatiquement
JYG & AJ Sep 27, 2018
Cette note a pour but de détailler les calculs permettant de déterminer la température d'une particule soulévée adiabatiquement lorsqu'il y a formation de glace. Elle se place dans la suite du travail d'Arnaud Jam et en particulier de son mémoire de 1998. Plus précisément, il s'agit de (re-)construire la procédure de calcul à mettre dans ündilute2" (ou dans "tlift" pour les anciennes versions) après l'appel de la procédure ïcefrac". L'idée générale dans tous ces calculs est que les particules d'air ascendant sont loin de l'équilibre et que leur phase condensée peut, à un même niveau, présenter une fraction englacée allant de 0% à 100%. La première phase du calcul ayant consisté à déterminer les caractéristiques des particules sans glace il s'agit ici de calculer la température des particules comportant la fraction englacée déterminée par la procédure ïcefrac".
On considère donc une particule de contenu en eau total $$q_t$$, dont on connait la température $$T_{p,ni}$$ et l'humidité massique $$q_{p,ni}$$ en l'absence de formation de glace (suffixe "ni" pour "no ice") et le contenu massique en glace, $$q_i$$. On se donne aussi la température $$T$$ de l'environnement au même niveau (cette température est supposée proche de la température de la particule avec glace).
On utilise les chaleurs latentes à la température $$T$$ :
| vaporisation | $$L_v = L_{v0} + (C_{pv}−C_l)(T−T_0)$$ |
| sublimation | $$ Ls = L_{s0} + (C_{pv}−C_i)(T−T_0)$$ |
| fusion | $$ L_f = L_s − L_v = L_{s0} − L_{v0} + (C_l−C_i)(T−T0) $$ |
où $$L_{v0}$$ et $$L_{s0}$$ sont les chaleurs latentes de vaporisation et de sublimation au point triple $$T_0$$ de l'eau ($$T_0$$ = 273.15 K). Les pressions de vapeur saturante, $$e_s$$ sur l'eau liquide et $$e_{si}$$ sur la glace sont données par les formules :
$$ e_s(T) = 6.112 exp[17.67 \frac{T−T_0} {T−29.65} ]$$
$$e_{si}(T)= exp[23.33086 − \frac{6111.72784} {T} + 0.15215 log(T)] $$ (2)
Les graphes de ces fonctions sont portés sur la figure. Ils se croisent évidemment au point triple (T = 273.15 K, $$e_s$$ = $$e_{si}$$ = 6.112 hPa).
Figure 1: Graphes des fonctions $$e_s$$ et $$e_{si}$$.
Avec $$e_s$$ et $$e_{si}$$ on peut écrire les humidités massiques à saturation :
$$q_{sat} = \frac{ϵe_s} {p − (1−ϵ)e_s} ≅ e_s \frac{ϵ} {p}$$ et $$ q_{sat,i} = \frac{ϵe_{si}} {p − (1−ϵ)e_{si}}≅ e_{si} \frac {ϵ}{p}$$
où $$ ϵ = R_d/R_v ≅ 0.62 $$ (noter que $$ q_{p,ni} = q_{sat}(T_{p,ni)})$$.
L'application de la formule de Clausius-Clapeyron donne en outre :
$$ \frac{d q_{sat,i}}{d T}(T) = \frac {L_s q_{sat,i}}{R_v T^2}$$
Les énergies statiques s'écrivent :
$$h_{p,ni} = [(1−q_t)C_{pd}+q_tC_l] T_{p,ni} + L_v q_{sat}(T_{p,ni}) + gz $$
pour la particule sans glace et
$$ h_p= [(1−q_t)C_{pd}+q_tC_l] T_p + L_v q_{sat,i}(T_p) − L_f q_i + gz $$
pour la particule avec glace.
En écrivant que l'énergie statique est conservée lors de la formation de glace, on obtient :
$$ [(1−q_t)C_{pd}+q_tC_l] T_{p,ni} + L_v q_{sat}(T_{p,ni}) = [(1−q_t)C_{pd}+q_tC_l] T_p + L_v q_{sat,i}(Tp) − L_f q_i $$
ce qui donne l'équation à résoudre en $$T_p$$ :
$$ [(1−q_t)C_{pd}+q_tC_l] (T_p − T_{p,ni}) = L_v (q_{sat}(T_{p,ni})−q_{sat,i}(T_p)) + L_f q_i$$ (3)
Résolution
On résoud l´équation par la méthode de Newton : on construit une suite $$T_j$$ partant de $$T_1 = ~T$$ et convergeant vers $$T_p$$ ; on passe de $$T_j$$ a $$T_{j+1}$$ par l'équation obtenue en linéarisant l'équation (3) :
$$ [(1−q_t)C_{pd}+q_tC_l] (T_{j+1} − T_{p,ni}) = L_v [q_{sat}(T_{p,ni}) − q_{sat,i}(T_j)− \frac{d q_{sat,i}} {d T} (T_j)(T_{j+1}−T_j) ] + L_f q_i$$ (4)
Ce qui donne :
| $$ [(1−q_t)C_{pd}+q_tC_l + L_v \frac {d q_{sat,i}} {d T}(T_j)] (T_{j+1}−T_j) $$ | $$ =[(1−q_t)C_{pd}+q-tC-l](T-{p,ni} − T_j) $$ |
| $$ + L_v [q_{sat}(T_{p,ni}) − q_{sat,i}(T_j)] $$ | |
| $$ + L_f q_i $$ |
(5)
soit :
$$ T_{j+1} = T_j + \frac { [(1−q_t)C_{pd}+q_tC_l](T_{p,ni} − T_j) + L_v [q_{sat}(T_{p,ni})−q_{sat,i}(T_j)] + L_f q_i} {(1−q_t)C_{pd}+q_tC_l + L_v \frac {d q_{sat,i}}{ d T} (T_j)} $$ (6)
et, finalement, en utilisant Clausius-Clapeyron :
$$ T_{j+1} = T_j + \frac { [(1−q_t)C_{pd}+q_tC_l](T_{p,ni} − T_j) + L_v [q_{p,ni}−q_{sat,i}(T_j)] + L_f qi} {(1−q_t)C_{pd}+q_tC_l +\frac {L_v L_s q_{sat,i}(T_j)}{ R_v T^2}} $$ (7)
Dans le code
La température de la particule avec glace $$T_p$$ ainsi que les approximations successives $$T_j$$ sont notées $$tp(i,k)$$. Les caractéristiques de la particule sans glace sont notées $$tg$$ pour $$T_{p,ni}$$ et $$qg$$ pour $$q_{p,ni}$$. L'eau totale est notées $$qnk(i)$$. Enfin la température $$T$$ de l'environnement est notée $$t(i,k)$$.
Tout cela donne la boucle :
alv=lv0-clmcpv*(t(i,k)-273.15)
alf=lf0+clmci*(t(i,k)-273.15)
als=alf+alv
tg=tp(i,k) ! On garde dans tg la temperature sans glace
tp(i,k) = t(i,k) ! L'iteration part de la temperature de l'environnement
do j=1,4
esi=exp(23.33086-(6111.72784/tp(i,k))
: +0.15215*log(tp(i,k)))
qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*als*qsat_new/
: (rrv*tp(i,k)*tp(i,k))
snew=1./snew
tp(i,k)=tp(i,k)+
: ( (cpd*(1.-qnk(i))+cl*qnk(i))*(tg-tp(i,k)) +
: alv*(qg-qsat_new) + alf*qi(i,k) )*snew
enddo
L'expérience montre que l'itération converge en moins de quatre coups à mieux que 0.01 K.
Remarques
- Les formules donnant les pressions de vapeur saturante devraient s'écrire de façon plus physique :
$$ e_{s0} $$ = 6.112 hPa Pression de vapeur saturante au point triple
$$e_s(T) = e_{s0} exp[17.67 \frac{ T−T_0}{T−29.65}] $$
$$ e_{si}(T) = e_{s0} exp[−6111.72784(\frac {1}{T} − \frac{1}{T0}) + 0.15215 log(\frac{T}{T0} )]$$ (8)
Le fait que l'itération converge prouve que l'on résoud bien l'équation (3). La question est alors de savoir si l'équation (3) est bien correcte. Les doutes portent sur l'écriture de l'énergie statique (dans le rapport d'Arnaud en 1998 il y a un un signe + devant le terme $$L_f q_i$$) et sur la dépendance en température de $$L_f$$ (signe contraire pour Arnaud). Il y a aussi un doute sur l'application de la relation de Clausius-Clapeyron à la sublimation mais, comme cela ne sert qu'à la résolution de l'équation (3), c'est sans importance.
La procédure Icefrac détermine un contenu total en glace. Comme ce contenu est déterminé avant le calcul de la température, on risque d'avoir un contenu en glace supérieur au contenu condensé. Il serait plus correct que Icefrac détermine la proportion de condensat qui est sous forme glace, puis qu'on résolve l'équation en température en tenant compte de la variation du contenu total en glace. En désignant par fi la fraction de condansat qui est sous forme glace, l'équation (3) deviendrait :
$$ [(1−q_t)C_{pd}+q_tC_l] (T_p − T_{p,ni}) = L_v (q_{sat}(T_{p,ni})−q_{sat,i}(T_p)) + L_f f_i (q_t − q_{sat,i}(T_p))$$ (9)
