Un modèle thermique simple à la surface de la Planète

1 Introduction

La Terre absorbe de l’énergie en provenance du Soleil. De l’énergie s’accumule sous sa surface. Elle évacue de l’énergie vers le Cosmos par différents mécanismes de transfert de la chaleur qui fonctionnent en parallèle: conduction, évaporation, convection et rayonnement. Les parts de chacun de ces mécanismes dans le total de l’énergie évacuée sont mal connues. Tout ce que l’on sait, c’est que plus la surface de la Terre est chaude, plus l’énergie évacuée est importante.

Il n’est pas inutile de rappeler que le Soleil est de très loin la principale source d’énergie de la Terre. On considère que le rayonnement solaire représente 99,97% du budget énergétique de la Terre. Voir ici et ici.

2 Equation différentielle du modèle

Pour faire le bilan énergétique de ce schéma, il suffit de constater que la variation d’énergie accumulée est égale à l’énergie absorbée moins l’énergie évacuée.

Plaçons nous en un point quelconque de la surface de la Terre. Supposons qu’à l’instant $t$, la température de surface y soit égale à $T$, et qu’elle ait la valeur $T+dT$ à l’instant $t+dt$.

En ignorant les échanges d’énergie avec les régions voisines, pour un petit élément de surface autour de ce point, on peut considérer que la variation d’énergie accumulée soit proportionnelle à la variation de température dT. Soient $E_{abs}$ et $E_{ev}$ les énergies absorbée et évacuée pendant l’intervalle de temps $dt$, on a:

$$\begin{array}{}\text{(1)}& CdT=E_{abs}-E_{ev}\end{array}$$

Si $P_{abs}$ est la puissance absorbée et $P_{ev}$ la puissance évacuée ,

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{E}_{abs}& =P_{abs}dt\\ E_{ev}& =P_{ev}dt\end{array}\end{array}$$

L’Eq. 1 devient

$$\begin{array}{}\text{(3)}& \frac{dT}{dt}=\frac{1}{C}(P_{abs}-P_{ev})\end{array}$$

En utilisant la notation de Stockwell

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{P}_{abs}& =S\\ P_{ev}& =F\end{array}\end{array}$$

On obtient

$$\begin{array}{}\text{(5)}& \frac{dT(t)}{dt}=\frac{1}{C}(S(t)-F(T))\end{array}$$

Si la Terre est soumise à un rayonnement constant, un équilibre thermique finit par s’établir. A cet instant, la puissance évacuée est égale à la puissance absorbée, et la température reste constante:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\frac{dT(t)}{dt}& =0\\ S& =F\end{array}\end{array}$$

On ne peut pas tirer grand chose de plus de cette équation. $S(t)$ dépend directement du rayonnement solaire au sommet de l’atmosphère, auquel il faut appliquer un coefficient mal connu, pour tenir compte du passage dans l’atmosphère et de la réflexion au niveau de la surface. La loi $F(T)$ est très mal connue. Elle contient un terme en $T^4$ correspondant à la puissance évacuée par rayonnement thermique, et un ou plusieurs termes de puissance inférieure pour les autres modes d’évacuation de la chaleur.

Par contre, si on s’intéresse aux perturbations par rapport à une situation d’équilibre, on peut en tirer beaucoup d’enseignements.

Supposons donc qu’au temps $t=0$, $T=T_0$, $S=S_0$ et $F=F_0=S_0$. Linéarisons $F(T)$ par un développement de Taylor au premier ordre. Si on pose $\frac{1}{\lambda }=F^{\prime }(T_0)$

$$\begin{array}{}\text{(7)}& F(T)=F_0+\frac{T-T_0}{\lambda }\end{array}$$

Posons encore $T_a=T-T_0$ (l’anomalie de température par rapport à $T_0$), $S_a=S-S_0$ (l’anomalie de rayonnement par rapport à $S_0$) et en notant que $d{T}_a=dT$ et que $F_0=S_0$ , l’Eq. 7 s’écrit

$$\begin{gathered} \begin{array}{}\text{(8)}& eq-05 \\ \frac{d{T}_a(t)}{dt}=\frac{1}{C}(S_a(t)-\frac{T_a(t)}{\lambda })\end{array} \end{gathered}$$

En posant $k=\lambda C$, la solution de cette équation est donnée par

$$\begin{gathered} \begin{array}{}\text{(9)}& eq-06 \\ {T}_a(t)=\frac{\lambda }{k}{e}^{-t/k}{\int }_0^t{e}^{u/k}{S}_a(u)du\end{array} \end{gathered}$$

Il est mathématiquement prouvé que cette équation est convergente vers une solution asymptotique lorsque le coefficient k est positif.

3 Solutions analytiques

3.1 Fonction en escalier

Si le système est à l’équilibre à l’instant $t=0$ à la température $T_0$ et que $S_a(t)=\mathrm{\Delta }{S}_0$ pour tout $t>0$ , alors

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathrm{\Delta }T(t)& =T(t)-T_0\\ & =\frac{\lambda }{k}{e}^{-t/k}k\mathrm{\Delta }{S}_0(e^{t/k}-1)\\ & =\lambda \mathrm{\Delta }{S}_0(1-e^{-t/k})\end{array}\end{array}$$

La variation de température évolue asymptotiquement vers la valeur $\lambda \mathrm{\Delta }{S}_0$. Voir Fig. 1.

Fig. 1: Eq. 10 avec $\lambda =\mathrm{\Delta }{S}_0=1$, pour différentes valeurs de k

3.2 Fonction générale

Une fonction générale $S_a(t)$ peut être discrétisée en plusieurs escaliers. Si on la suppose constante dans les intervalles $i,i+1$, en utilisant l’Eq. 10, on obtient

$$\begin{gathered} \begin{array}{}\text{(11)}& T(t)-T_0=\lambda \sum _i\mathrm{\Delta }{S}_i(1-e^{-(t-i)/k}) \\ \mathit{avec} \\ \mathrm{\Delta }{S}_i=S_{i+1}-S_i=S_{a_{i+1}}-S_{a_i}\end{array} \end{gathered}$$

Si l’on approche $(1-e^{-t/k})$ par une fonction linéaire-palier égale à $t/\tau$ pour $t$ compris entre 0 et $\tau$, et égale à 1 pour t supérieur à $\tau$ (voir Fig. 2), l’Eq. 11 peut être fortement simplifiée. En effet, la formule peut alors s’écrire sous forme de différences finies:

$$\begin{gathered} \begin{array}{}\text{(12)}& T(t)-T_0=\lambda (T_1+T_2) \\ {T}_1=\sum _{i=0,1,...,t-\tau }\mathrm{\Delta }{S}_i \\ {T}_2=\sum _{i=0,1,...,\tau }\mathrm{\Delta }{S}_{t-\tau +i}\frac{\tau -i}{\tau }\end{array} \end{gathered}$$

En notant

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\overline{S}}_{t-\tau }^t& =lamoyennedeSentrelesinstantstett-\tau \\ {\overline{S}}_{a_{t-\tau }}^t& =lamoyennede{S}_{a}entrelesinstantstett-\tau \end{array}\end{array}$$

En développant les sommes et en regroupant les termes , on obtient facilement

$$\begin{gathered} \begin{array}{}\text{(14)}& T_1=S_{t-\tau }-S_0 \\ {T}_2={\overline{S}}_{t-\tau }^t-S_{t-\tau }\end{array} \end{gathered}$$

Il vient finalement

$$\begin{array}{}\text{(15)}& T(t)-T_0=\lambda ({\overline{S}}_{t-\tau }^t-S_0)=\lambda {\overline{S}}_{a_{t-\tau }}^t\end{array}$$

Fig. 2: Exemple d’approximation d’une fonction exponentielle par une fonction linéaire au départ et constante par la suite.

3.3 Variation sinusoïdale

Supposons que $S_a(t)$ soit une sinusoïde d’amplitude $\mathrm{\Delta }{S}_a$ et de période $per$ :

$$\begin{array}{}\text{(16)}& S_a(t)=\mathrm{\Delta }{S}_{a}sin(2\pi \frac{t}{per})\end{array}$$

Dans ce cas la solution est donnée par

$$\begin{array}{}\text{(17)}& \mathrm{\Delta }T(t)=T(t)-T_0=e^{-t/k}Asin(\varphi )+Asin(\omega t-\varphi )\end{array}$$

Dans cette dernière formule,

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\omega & =\frac{2\pi }{per}\\ \varphi & =atan(k\omega )\\ A& =\frac{1}{\sqrt{1+(k\omega {)}^2}}\lambda \mathrm{\Delta }{S}_a\end{array}\end{array}$$

Le premier terme de l’Eq. 17 est un terme transitoire. Il tend vers 0 lorsque $t$ tend vers l’infini, et $\mathrm{\Delta }T(t)$ tend vers une sinusoïde déphasée d’un angle $\varphi$ par rapport au rayonnement solaire.

Pour des oscillations lentes, $k\omega$ tend vers zéro, l’amplitude de la variation de température tend vers $\lambda \mathrm{\Delta }{S}_a$ et $\varphi$ tend vers zéro.

Pour des oscillations rapides, $k\omega$ tend vers l’infini, l’amplitude de la variation de température tend vers zéro et $\varphi$ tend vers $\pi /2$.

Voir Fig. 3.

Fig. 3: Anomalie de la température pour une anomalie solaire sinusoïdale (Eq. 17, sans le terme transitoire), avec $\lambda =\mathrm{\Delta }{S}_a=1,etk=5$, pour différentes valeurs de la période $per$. En abscisse, le temps et en ordonnée l’anomalie de température, sans unités particulières.
On notera que l’amplitude de l’anomalie de température est très faible lorsque la période du signal d’entrée est courte. C’est l’inertie thermique qui est à l’oeuvre.

3.4 Variation sinusoïdale déphasée

Si $\mathrm{\Delta }S(t)$ est déphasée d’un angle $\alpha$ :

$$\begin{array}{}\text{(19)}& \mathrm{\Delta }S(t)=\mathrm{\Delta }{S}_{a}sin(2\pi \frac{t}{per}-\alpha )\end{array}$$

Dans ce cas la solution de l’Eq. 9 devient

$$\begin{array}{}\text{(20)}& \mathrm{\Delta }T(t)=T(t)-T0=e^{-t/k}Asin(\varphi +\alpha )+Asin(\omega t-(\varphi +\alpha ))\end{array}$$

Les paramètres de cette formule sont également calculés selon l’Eq. 18.

3.5 Fonction périodique quelconque

Si $S_a(t)$ est une fonction périodique, elle peut être décomposée en série de Fourier:

$$\begin{array}{}\text{(21)}& S_a(t)=\mathrm{\Delta }{S}_0+\sum _{i=1..n}\mathrm{\Delta }{S}_{i}sin({\omega }_{i}t+{\beta }_i)\end{array}$$

Le terme $\mathrm{\Delta }{S}_0$ correspond à l’harmonique 0 de la décomposition et est égal à la moyenne de $S_a(t)$ sur sa période. La solution s’obtiendra en lui appliquant l’Eq. 10 et l’Eq. 19 aux autres termes.

3.5.1 Exemple: rayonnement solaire incident à l’équinoxe du printemps ou d’automne

Aux latitudes tempérées, à un facteur d’échelle près, le rayonnement est très bien modélisé par une demi sinusoïde pendant les 12 heures qui suivent le lever du soleil, et une valeur nulle les 12 heures suivantes. La moyenne de cette fonction est égale à $1/\pi$. L’anomalie de ce rayonnement peut s’exprimer comme suit:

$$\begin{array}{}\text{(22)}& \begin{array}{rlr}{S}_a(t)& =\mathrm{\Delta }{S}_0(sin(\omega t)-1/\pi )& 0<=t<=\frac{\pi }{\omega }\\ & =-\mathrm{\Delta }{S}_0/\pi & \frac{\pi }{\omega }<=t<=\frac{2\pi }{\omega }\end{array}\end{array}$$

Pour $\mathrm{\Delta }{S}_0=1$ , la transformée de Fourier de cette fonction possède une solution analytique donnée par la formule suivante:

$$\begin{array}{}\text{(23)}& fft(S_a(t))=\frac{1}{2}sin(\omega t)+\frac{2}{\pi }\sum _{n=2,4,6,....}\frac{1}{n^2-1}sin(n\omega t-\frac{\pi }{2})\end{array}$$

Il suffira d’appliquer les formules précédentes à chaque terme. Voir Fig. 4.

Fig. 4: Elle reprend le rayonnement solaire théorique un jour d’équinoxe. Ce rayonnement est approché par une demi-sinusoïde pendant 12 heures après le lever du soleil, suivie par une valeur nulle pendant la nuit. Il est représenté par la courbe Sa.
L’anomalie de température a été calculée en appliquant les formules de l’Eq. 18 et de l’Eq. 20 à chaque terme de l’Eq. 23. Les calculs ont été menés en considérant $\lambda =1$, pour quelques valeurs de $k$. En abscisse, le temps exprimé en heures depuis le lever du soleil. En ordonnée, l’anomalie du rayonnement solaire et les anomalies de température (sans unités particulières).

4 Analogies au modèle thermique

4.1 Vidange d’un récipient alimenté en eau

Considérons un récipient de section constante alimenté en eau par un robinet, et muni à sa base d’un orifice de vidange. Le bilan de l’eau dans le récipient est facile à établir.

En notant:

$S$	La section du récipient
$s$	La section du trou de vidange
$h$	La hauteur d’eau au-dessus du trou de vidange
$D_{in}$	Le débit d’alimentation
$D_{out}$	Le débit évacué

On a:

$$\begin{array}{}\text{(24)}& \frac{dh}{dt}=\frac{1}{S}(D_{in}-D_{out})\end{array}$$

C’est analogue à l’Eq. 3. Il suffit de considérer qu’un volume d’eau est de l’énergie, un débit d’eau est une puissance, et la hauteur d’eau est la température.

En utilisant la formule de Torricelli, on trouve facilement que

$$\begin{array}{}\text{(25)}& D_{out}=s\sqrt{2gh}\end{array}$$

A l’équilibre, pour un débit d’alimentation constant, on obtient

$$\begin{array}{}\text{(26)}& h=\frac{1}{2g}(D_{in}/s{)}^2\end{array}$$

Finalement, l’Eq. 24 s’écrit

$$\begin{array}{}\text{(27)}& \frac{dh}{dt}=\frac{1}{S}(D_{in}-s\sqrt{2gh})\end{array}$$

Comme précédemment, cette équation peut être linéarisée par rapport à une situation d’équilibre. Supposons qu’au temps $t=0$, $h=h_0$, $D_{in}=D_{i{n}_0}$ et $D_{out}=D_{ou{t}_0}=s\sqrt{2g{h}_0},D_{ou{t}_0}^{\prime }=s\sqrt{g/2h_0}$ .

En notant $h^a=h-h_0,D_{in}^a=D_{in}-D_{i{n}_0}$, il vient finalement

$$\begin{array}{}\text{(28)}& \frac{d{h}^a}{dt}=\frac{1}{S}(D_{in}^a-h^{a}s\sqrt{g/2h_0})\end{array}$$

Ceci permettra de juger dans quelle mesure l’approximation par linéarisation s’écarte de la solution exacte.

4.1.1 Un exemple concret pour un débit d’alimentation sinusoïdal

4.1.1.1 Solution exacte

L’Eq. 24 a été intégrée numériquement pour un récipient d’une section $S=100c{m}^2$, muni d’un trou de section $s=0.5c{m}^2$ alimenté par un débit variable $D_{in}=D_0(1+ksin(2\pi t/per))$ avec $D_0=3l/min$. La hauteur d’eau initiale $h_0$ a été calculée par l’Eq. 26 avec $D_{in}=D_0$.

Fig. 5: Solution de l’Eq. 24 par intégration numérique pour un cas concret, pour $k=0.5$ et quelques valeurs de $per$. On peut observer la convergence vers une sinusoïde qui oscille autour de la valeur moyenne du débit en entrée. L’amplitude et le déphasage de la sinusoïde sont analogues à ce que l’on observe sur la Fig. 3

4.1.1.2 Solution approchée par linéarisation

Fig. 6: Solution de l’Eq. 28 par intégration numérique pour un cas concret, pour $k=0.5$ et quelques valeurs de $per$. Les différences avec la Fig. 5 sont quasiment imperceptibles.

Fig. 7: Les erreurs d’approximation par linéarisation sont négligeables.

4.1.2 Réponse à une impulsion

Supposons que le même récipient que précédemment ait atteint l’équilibre pour un débit d’alimentation $D_0=3l/min$ , et que pendant une brève période, ce débit soit doublé.

Fig. 8: Solution exacte de l’Eq. 28 par intégration numérique, pour un doublement du débit pendant 5 secondes, en partant d’une situation à l’équilibre. Le niveau augmente rapidement pendant le doublement avant de décroître beaucoup plus lentement vers l’équilibre initial.

Un phénomène analogue est vraisemblablement à l’oeuvre dans les événements de Dansgaard-Oeschger pendant lesquels on observe une augmentation de la température de 5 à 8 °C en quelques décennies suivie d’une période de refroidissement s’étendant sur plusieurs centaines d’années.

4.2 Equation du mouvement d’un volant d’inertie

Les vélos d’intérieur sont souvent équipés d’un volant d’inertie qui dissipe par courants de Foucault la puissance produite par le cycliste. L’équation différentielle d’un tel système soumis à un couple moteur et à un couple résistant est facile à poser en se basant sur le principe de conservation de l’énergie.

En notant:

éééééé

$$\begin{array}{}\text{(29)}& \begin{array}{rl}\theta & =Positionangulaireduvolant\\ \omega =\dot{\theta }& =Vitesseangulaireduvolant\\ \dot{\omega }=\ddot{\theta }& =Accélérationangulaireduvolant\\ J& =Inertierotationnelleduvolant\\ E_v=J{\omega }^2/2& =Energiecinétiqueduvolant{d}^{\prime }inertie\\ C_m& =Couplemoteur\\ C_r& =Couplerésistant\\ P_m=C_m\omega & =Puissancemotrice\\ P_r=C_r\omega & =Puissancerésistante\\ E_m=P_{m}dt& =Energieinjectéedanslevolantpendant{l}^{\prime }intervalledt\\ E_r=P_{r}dt& =Energieperdueparlevolantpendant{l}^{\prime }intervalledt\end{array}\end{array}$$

Le bilan énergétique du volant s’écrit:

$$\begin{array}{}\text{(30)}& \begin{array}{rl}d{E}_v& =E_m-E_r\\ J\omega d\omega & =P_{m}dt-P_{r}dt\\ J\omega \dot{\omega }& =C_m\omega -C_r\omega \\ \\ \dot{\omega }=\ddot{\theta }& =\frac{1}{J}(C_m-C_r)\end{array}\end{array}$$

Cette équation est équivalente à l’Eq. 5. Elle permet de mieux comprendre le sens qu’il faut donner à l’inertie thermique.

5 Détermination des paramètres du modèle sur base des observations

On dispose de plusieurs séries de température. Elles sont globales ou locales. Certaines séries sont publiées sous forme de grille qui couvre tout ou partie de la surface terrestre. Si les températures sont absolues, il sera facile d’en dériver une anomalie par différence avec une moyenne.

Le rayonnement absorbé est mal connu. Dans certaines stations, le rayonnement incident est mesuré depuis de longues années, mais pas le rayonnement réfléchi.

En l’absence de mesures locales du rayonnement incident, il faudra se baser sur l’irradiance totale solaire (ITS). Celle-ci n’est mesurée par satellite que depuis quelques décennies, mais il existe des reconstructions historiques de l’ITS qui remontent bien plus loin dans le passé. Lorsque le calcul est fait en un point précis, on pourra passer de l’ITS au rayonnement incident théorique local par des calculs d’astronomie. Il existe des librairies spécialisées qui facilitent grandement ces calculs.

On peut utiliser les Eq. 8, Eq. 9, Eq. 15 dans tous les cas, que l’activité solaire soit périodique ou non.

L’Eq. 8 requiert de passer par une intégration numérique qui peut être assez lourde.

L’Eq. 9 est déconseillée parce que la multiplication de très petits nombres par de très grands nombres provoquera généralement des instabilités numériques liées.

L’Eq. 15 est assez légère à manipuler, mais nécessite de connaître l’historique de l’activité solaire sur une période suffisamment longue avant le début de la série de température.

Pour les activités solaires périodiques (journalière ou saisonnière), le plus simple sera d’effectuer une décomposition du signal en série de Fourier.

Dans tous les cas, il faudra passer par une régression multiple non-linéaire pour obtenir les paramètres du modèle.

Si l’on dispose de données d’évaporation, il faudra ajouter un terme correspondant dans les régressions. Sinon l’évaporation va se retrouver de manière implicite dans le coefficient du rayonnement et/ou dans le terme indépendant. L’évaporation est loin d’être négligeable.