La troposphère : profil de température, intégration numérique, gradient adiabatique

La stratégie. L'équation locale et le gaz parfait donnent toujours $\mathrm{d}P/\mathrm{d}z=-MgP/(RT)$ ; seule change la fonction $T(z)$, désormais imposée par les mesures. Les variables se séparent : solution analytique dans la troposphère, solution numérique pour le profil complet, qui sert aussi à contrôler la première. La dernière question applique la loi de Laplace à une bulle et la poussée d'Archimède à sa stabilité.

1) L'équation et sa solution. Avec $\rho =\frac{PM}{RT(z)}$, l'équation locale devient $$\frac{\mathrm{d}P}{\mathrm{d}z}=-\frac{Mg}{RT(z)}P.$$ Dans la troposphère, $\frac{\mathrm{d}P}{P}=-\frac{Mg}{R}\frac{\mathrm{d}z}{T_0-\Gamma z}$, et l'on intègre de $0$ à $z$ : $$\ln \frac{P}{P_0}=\frac{Mg}{R\Gamma }\ln (1-\frac{\Gamma z}{T_0})⟹P(z)=P_0(1-\frac{\Gamma z}{T_0}{)}^{\alpha },\alpha =\frac{Mg}{R\Gamma }.$$ $$\alpha =\frac{29,0\times 10^{-3}\times 9,81}{8,314\times 6,5\times 10^{-3}}\approx 5,26.$$ $\alpha$ est sans dimension : $Mg/R$ est un gradient de température ($34,2\mathrm{K}\mathrm{k}{\mathrm{m}}^{-1}$), divisé par $\Gamma$.

2) Cas limite. Avec $u=\Gamma z/T_0\to 0$ : $$\ln \frac{P}{P_0}=\frac{Mg}{R\Gamma }\ln (1-u)\approx -\frac{Mg}{R\Gamma }\frac{\Gamma z}{T_0}=-\frac{Mgz}{R{T}_0},$$ c'est l'atmosphère isotherme à $T_0$. Le résultat général contient bien le cas particulier.

3) Everest et tropopause.

Par exemple $(1-6,5\times 8,849/288,15{)}^{5,26}=(0,8004{)}^{5,26}\approx 0,310$, et $e^{-8849/8417}\approx 0,349$. Le modèle isotherme surestime la pression en altitude : l'air réel, plus froid, est plus dense à pression égale, il pèse davantage par mètre, et la pression décroît plus vite. L'erreur croît avec l'altitude, parce que l'écart de température s'accumule. Au sommet de l'Everest, un tiers de la pression du sol : c'est la pression mesurée, à quelques hectopascals près selon la saison.

4) Capacité numérique. La méthode d'Euler remplace la dérivée par un taux d'accroissement sur un pas $\mathrm{d}z$ : $P(z+\mathrm{d}z)\approx P(z)-\frac{MgP(z)}{RT(z)}\mathrm{d}z$. L'échelle de variation de $P$ est $H\approx 8\mathrm{km}$ : un pas de $10\mathrm{m}$ est mille fois plus petit.

import numpy as np
import matplotlib.pyplot as plt

M, g, R, P0 = 29.0e-3, 9.81, 8.314, 1.013e5   # SI
T0, Gam = 288.15, 6.5e-3                     # K, K/m

def T(z):                                    # atmosphere type, 0-20 km
    return T0 - Gam * min(z, 11000.0)

def euler(dz):
    z = np.arange(0.0, 20000.0 + dz / 2, dz)
    P = np.empty_like(z); P[0] = P0
    for i in range(len(z) - 1):              # dP/dz = -M g P / (R T)
        P[i+1] = P[i] - M * g * P[i] / (R * T(z[i])) * dz
    return z, P

z, P = euler(10.0)
Piso = P0 * np.exp(-M * g * z / (R * 288.0))
def alt(Q, c):                               # 1re altitude ou Q <= c
    return z[np.argmax(Q <= c)] / 1000
print("P(11 km) =", round(P[1100]), "Pa ; analytique :",
      round(P0 * (1 - Gam * 11000 / T0) ** (M * g / (R * Gam))))
print("P(20 km) =", round(P[-1]), "Pa ; isotherme :", round(Piso[-1]))
for k in (2, 10):
    print("P0/" + str(k), ":", alt(P, P0 / k), "km ; isotherme :",
          alt(Piso, P0 / k), "km")
print("pas /10 : P(20 km) =", round(euler(1.0)[1][-1]), "Pa")

plt.plot(P / 100, z / 1000, label="profil type (Euler)")
plt.plot(Piso / 100, z / 1000, "--", label="isotherme 288 K")
plt.xlabel("P (hPa)"); plt.ylabel("z (km)"); plt.legend(); plt.show()

Le programme définit le profil à deux couches, intègre l'équation par Euler avec un pas de $10\mathrm{m}$, compare au modèle isotherme, recherche les altitudes demandées (première altitude où la pression passe sous le seuil), puis recommence avec un pas de $1\mathrm{m}$. Il affiche :

P(11 km) = 22553 Pa ; analytique : 22572
P(20 km) = 5437 Pa ; isotherme : 9410
P0/2 : 5.47 km ; isotherme : 5.84 km
P0/10 : 16.07 km ; isotherme : 19.38 km
pas /10 : P(20 km) = 5447 Pa

Contrôle. À $11\mathrm{km}$, Euler donne $22553\mathrm{Pa}$ contre $22572$ pour la formule de la question 1 : $0,08\%$ d'écart. Avec un pas dix fois plus petit, $P(20\mathrm{km})$ passe de $5437$ à $5447\mathrm{Pa}$, et la solution analytique (la troposphère, puis une exponentielle à $216,65\mathrm{K}$ de $11$ à $20\mathrm{km}$) donne $5448\mathrm{Pa}$ : l'erreur est divisée par dix avec le pas, comme il se doit pour la méthode d'Euler, dont l'erreur est proportionnelle au pas. (L'atmosphère type officielle, avec $P_0=101325\mathrm{Pa}$ et ses propres constantes, donne $22632\mathrm{Pa}$ à $11\mathrm{km}$ : notre $M$ arrondi et notre $g$ expliquent les trois pour mille.)

Comparaison. La moitié de la pression est atteinte à $5,47\mathrm{km}$ au lieu de $5,84$, le dixième à $16,1\mathrm{km}$ au lieu de $19,4$ ; à $20\mathrm{km}$, l'isotherme surestime la pression de $73\%$. Sur la courbe, les deux profils partent tangents (même pente au sol, $-\rho (0)g$), puis le profil réel passe sous l'isotherme ; au-dessus de $11\mathrm{km}$, il redevient une exponentielle, de hauteur caractéristique $RT/(Mg)\approx 6,3\mathrm{km}$ seulement.

5) Le gradient adiabatique. Le système est une bulle d'air sec, fermée. Sa pression égale celle de l'environnement à chaque altitude : $\frac{\mathrm{d}P}{\mathrm{d}z}=-\frac{MgP}{R{T}_{\text{env}}}$. La loi de Laplace, $\ln T+\frac{1-\gamma }{\gamma }\ln P=\text{cte}$, se différencie en $$\frac{\mathrm{d}T}{T}=\frac{\gamma -1}{\gamma }\frac{\mathrm{d}P}{P}=-\frac{\gamma -1}{\gamma }\frac{Mg}{R{T}_{\text{env}}}\mathrm{d}z.$$ Au départ, la bulle est à la température de l'air qui l'entoure, $T=T_{\text{env}}$, d'où $$\frac{\mathrm{d}T}{\mathrm{d}z}=-{\Gamma }_{\mathrm{a}},{\Gamma }_{\mathrm{a}}=\frac{\gamma -1}{\gamma }\frac{Mg}{R}=\frac{0,4}{1,4}\times \frac{29,0\times 10^{-3}\times 9,81}{8,314}\approx 9,8\times 10^{-3}\mathrm{K}{\mathrm{m}}^{-1},$$ soit $9,8\mathrm{K}\mathrm{k}{\mathrm{m}}^{-1}$. On peut l'écrire ${\Gamma }_{\mathrm{a}}=Mg/C_{P,\mathrm{m}}$, avec $C_{P,\mathrm{m}}=\gamma R/(\gamma -1)\approx 29,1\mathrm{J}{\mathrm{K}}^{-1}\mathrm{mo}{\mathrm{l}}^{-1}$ : l'énergie potentielle gagnée est prise sur l'enthalpie de la bulle.

Le point délicat. La pression de la bulle est celle de l'environnement : c'est $T_{\text{env}}$ qui figure dans $\mathrm{d}P/\mathrm{d}z$, et c'est seulement parce que la bulle part à la température ambiante que $T/T_{\text{env}}=1$ et que le gradient prend la valeur simple ${\Gamma }_{\mathrm{a}}$.

La stabilité. Soulevons une bulle de $\delta z$. Elle se refroidit de ${\Gamma }_{\mathrm{a}}\delta z$, l'air environnant n'est plus froid que de $\Gamma \delta z$ ; comme $\Gamma =6,5<{\Gamma }_{\mathrm{a}}=9,8\mathrm{K}\mathrm{k}{\mathrm{m}}^{-1}$, la bulle est plus froide que l'air qui l'entoure. À même pression, elle est plus dense : son poids l'emporte sur la poussée d'Archimède, qui est le poids de l'air déplacé, et elle redescend. Une bulle abaissée, plus chaude que son entourage, remonte. L'atmosphère type est stable vis-à-vis des mouvements verticaux d'air sec. Si le sol, chauffé par le soleil, rend le gradient réel supérieur à ${\Gamma }_{\mathrm{a}}$, une bulle soulevée continue de monter : c'est la convection, les ascendances des planeurs, les cumulus. (Dans l'air humide, la condensation libère de l'énergie et le gradient de la bulle tombe vers $5$ à $6\mathrm{K}\mathrm{k}{\mathrm{m}}^{-1}$ : la troposphère réelle, à $6,5$, est souvent instable pour l'air saturé.)

Ce que l'exercice installe. Résoudre $\mathrm{d}P/\mathrm{d}z=-MgP/(RT(z))$ pour un profil $T(z)$ donné : séparation des variables, cas limite $\Gamma \to 0$, puis méthode d'Euler contrôlée par la solution exacte et par la division du pas. Le modèle isotherme surestime la pression d'autant plus qu'on monte. Et la comparaison du gradient réel au gradient adiabatique, par la poussée d'Archimède, décide de la stabilité de l'atmosphère : on retrouvera cette étude de stabilité d'un équilibre avec les oscillateurs.