4 Une modélisation simple de l’anomalie de température globale en fonction du rayonnement solaire

Symboles utilisés dans les formules

\(t\)	Le temps
\(\Delta T(t)\)	L’anomalie de température au temps \(t\) , c-à-d la différence entre la température et une température moyenne, en général calculée sur une période de 30 ans. La série de données hadsst3 (global monthly) a été retenue.
\(a\)	Un coefficient
\(\tau\)	Une durée correspondant à l’inertie thermique des océans
\(\bar S_{t-\tau}^{\ \ t}\)	Le rayonnement solaire moyen entre le temps \(t\) et \(t-\tau\)
\(S_{ref}\)	Un rayonnement solaire de référence.
\(S(t)\)	Le rayonnement solaire au temps \(t\)

4.1 Un modèle statique

L’approche classique et statique peut s’exprimer comme suit:

\[ \Delta T(t)\ =\ a(S(t)\ -\ S_{ref}) \tag{4.1}\]

Elle suppose que l’équilibre thermique s’établit instantanément. Cette approche produit de très mauvaises corrélations et ne peut expliquer que de très faibles variations de température. Voir Figure 4.1.

Figure 4.1: Régression directe de la température océanique (hadsst3gm) par l’ITS (nrl2Tsi) Cette approche statique classique ne permet d’expliquer que des variations de température de quelques dixièmes de °C, avec une piètre corrélation. C’est sur cette base que le GIEC s’appuie pour conclure que l’ impact de l’activité solaire sur le climat est négligeable.

4.2 Un modèle dynamique

En remplaçant le rayonnement solaire instantané par le rayonnement solaire moyen des \(\tau\) dernières années (Voir Figure 2.4 et Figure 2.5 tracées avec \(\tau\) = 110 ans) on obtient un modèle dynamique.

\[ \Delta T(t)\ =\ a(\bar S^{\ \ t}_{t-\tau }\ -\ S_{ref})\ avec\ \bar S^{\ \ t}_{t-\tau } =\ \frac 1{\tau }\ \int _{t-\tau }^tS(t)dt \tag{4.2}\]

Utiliser le rayonnement solaire moyen des \(\tau\) dernières années équivaut à considérer que la variation de température consécutive à une variation de rayonnement solaire apparaît progressivement et n’est pleinement effective qu’après \(\tau\) années. Il faut du temps pour que la température océanique s’adapte à la variation de l’activité solaire. Cette approche simple est en fait la solution approchée d’un modèle thermique du climat. Voir Section 4.3.

Dans ce modèle, la planète est perpétuellement en état de déséquilibre thermique. (Abdussamatov 2013) a bien décrit les mécanismes physiques responsables de ce déséquilibre.

L’anomalie de température est nulle lorsque le rayonnement solaire moyen est égal au rayonnement solaire de référence. Lorsque le rayonnement solaire moyen est supérieur au rayonnement solaire de référence, l’anomalie de température est positive. Elle est négative dans le cas contraire.

La température ne se stabilise que si le rayonnement solaire reste constant pendant \(\tau\) années. Le modèle ne fait que lisser le rayonnement solaire, mais il peut produire des phases chaudes ou froides de grande amplitude et de longue durée, parce que de faibles variations du rayonnement solaire peuvent se cumuler via l’énorme inertie thermique des océans. De plus, l’intensité solaire peut diminuer, comme constaté lors des dernières années, sans que la température ne diminue immédiatement. Ce qui répond à une objection du GIEC réfutant sur cette base, le rôle moteur du soleil dans le mécanisme climatique.

Une parenthèse: l’océan comme volant d’inertie

Si dans l’Équation 4.2 on remplace \(\Delta T(t)\) par \(\dot \omega\), \(a\) par \(\frac{1}{J}\), \(\bar S_{t-\tau}^{\ \ \ t}\) par \(C_m(t,\tau)\) et \(S_{ref}\) par \(C_{r0}\), on obtient l’équation du mouvement d’un volant d’inertie rotationnelle \(J\) soumis à un couple moteur \(C_m(t,\tau)\) et à un couple résistant \(C_{r0}\) :

\[J \dot \omega = C_m(t,\tau) - C_{r0}\]

L’équation qui lie l’anomalie de température et le rayonnement solaire est donc analogue à celle qui décrit l’accélération angulaire d’un volant d’inertie.

Lorsqu’on utilise l’Équation 4.2 dans une analyse de régression, l’examen des résidus révèle la présence d’un terme harmonique, indépendant du rayonnement solaire, et qui pourrait résulter des boucles convectives verticales et horizontales (les « courants », le « conveyor belt », le gulf stream, les phénomènes NAO, PDO, etc.) existant au sein des océans et circulant respectivement entre couches superficielles et profondes, ou entre zones chaudes et froides; et dont on détecte la trace de la composante principale dans les résidus. En l’ajoutant à l’Équation 4.2 , on obtient :

\[ \Delta T(t)\ =\ a(\bar S^{\ \ t}_{t-\tau }\ -\ S_{ref})\ +\ A\cos (2\pi \frac{(t-t_0)}{per}) \tag{4.3}\]

Dans cette expression,

\(A\)	représente l’amplitude,
\(per\)	la période et
\(2 \pi \frac{t_0}{per}\)	la phase de l’oscillation

Les analyses de régression multiple effectuées optimisent tous les paramètres de l’Équation 4.2 ou de l’Équation 4.3 par une technique classique de minimisation des moindres carrés, les équations non linéaires étant résolues en optimisant des itérations successives.

Pour quantifier le rayonnement solaire on peut utiliser le nombre de taches solaires ou l’irradiance totale solaire (ITS).

4.3 Un modèle thermique équivalent au modèle dynamique

Il est intéressant de noter que la modélisation de l’anomalie de température par une moyenne du rayonnement solaire n’est qu’une autre formulation du modèle thermique utilisé par Stockwell et exposé dans un de ses autres articles (Stockwell 2011) :

\[ \frac{dT}{dt}\ =\ \frac 1 C\ (S\ -\ F) \tag{4.4}\]

Dans ce modèle, \(\frac{dT}{dt}\) représente la dérivée de la température par rapport au temps, S le rayonnement solaire, F une perte radiative et C une constante. La perte radiative est supposée proportionnelle à la température:

\[ F\ =\ \frac T{\lambda } \tag{4.5}\]

Cette hypothèse semble bizarre à première vue, la contribution radiative comprenant normalement un terme en T⁴ . Cette formule correspond cependant à une démarche relativement courante, lorsqu’on travaille aux environs de la température ambiante. Voir (McCabe, Smith, et Harriott 1993).

Elle revient à considérer que le terme radiatif est un complément au coefficient de transfert convectif hc, la formule de transfert de chaleur combinée radiative-convective devenant :

\[ \begin{gathered} \rho \ c_p\ \frac{dT}{dt}\ =\ (h_c\ +\ h_r)\ S_e\ (T_0\ -\ T)\\\mathit{o\text{ù}}\ h_r\ =\ 4\sigma \epsilon T^3 \end{gathered} \tag{4.6}\]

Dans ces formules :

\(\rho\)	masse spécifique [kg/m³]
\(c_p\)	chaleur spécifique [Joules /(kg°K)]
\(T\)	température [°K]
\(t\)	temps (h)
\(h_c\)	coefficient de transfert convectif (empirique) [W/(m² °K)]
\(h_r\)	coefficient de transfert radiatif “linéarisé” [W/(m²°K)]
\(S_e\)	surface d’échange [m²]
\(T_a\)	température de l’autre milieu avec lequel la chaleur est échangée [°K]
\(\sigma\)	la constante de Stefan-Boltzmann [\(\sigma\) = 5,670 374 × 10^-8 W m^-2 K^-4 ]
\(e\)	l’émissivité [0]

Moyennant cette linéarisation, obtient alors

\[ \frac{\mathit{dT}}{\mathit{dt}}\ +\ \frac T{\lambda C}\ =\ \frac S C \tag{4.7}\]

Posons \(\lambda \ C\ =\ k\)

Si l’on considère \(t=0\) lorsque \(\frac{dT}{dt}=0\) et qu’à cet instant,\(S=S_0\) et \(T=T_0=\lambda S_0\), une solution de l’équation différentielle est donnée par

\[ T(t)\ -\ T_0\ =\ \frac{e^{-\frac t k}} C\ \int _0^te^{\frac t k}\ \Delta S(t)\ \mathit{dt} \tag{4.8}\]

Dans cette formule, \(\Delta S(t)=S(t)-S_0\)

Si \(\Delta S(t)\) est une fonction escalier (\(\Delta S=\Delta S_0\) pour t >= 0), l’intégrale se calcule facilement. On obtient

\[ T(t)-T_0\ =\ \frac{e^{-\frac t k}} C\ k\ \Delta S_0\ (e^{\frac t k}\ -\ 1)\ =\ \lambda \ \Delta S_0\ (1\ -\ e^{-\frac t k}) \tag{4.9}\]

Après une longue période, la variation de température se stabilise :

\[ T(t)-T_0\ =\ \lambda \ \Delta S_0 \tag{4.10}\]

Si l’on discrétise la fonction S(t) en plusieurs escaliers, en la supposant constante dans les intervalles i, i+1 on obtient

\[ T(t)\ -\ T_0\ =\ \lambda \ \sum _i\ \Delta S_i\ (1\ -\ e^{-\frac{(t-i)} k})\mathit{\ avec\ }\Delta S_i=S_{i+1}-S_i \tag{4.11}\]

Si l’on approche \((1\ -\ e^{-\frac t k})\) par une fonction linéaire-palier égale à \(\frac t \tau\) pour t compris entre 0 et \(\tau\), et égale à 1 pour t supérieur à \(\tau\) (voir Figure 4.2),

Figure 4.2: Exemple d’approximation d’une fonction exponentielle par une fonction linéaire-escalier. Dans cet exemple. La valeur de k a été déduite de \(\tau\) par une analyse de régression non linéaire, pour une durée d’observation de 270 ans. Voir Section 4.4

la formule peut alors s’écrire sous forme de différences finies

\[ T(t)\ -\ T_0\ =\ \lambda \ (T_1\ +\ T_2) \tag{4.12}\]

Avec

\(T_1\ =\ \sum _i\Delta S_i\) pour i=0, 1, … , \(t-\tau\) et

\(T_2\ =\ \sum _i\Delta S_{t-\tau +i}\ \frac{\tau \ -\ i}{\tau }\) pour i=0, 1, … , \(\tau\)

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

\(T_1\ =\ S_{t-\tau }\ -\ S_0\)

\(T_2\ =\ \bar S^{\ \ t}_{t-\tau }\ -\ S_{t-\tau }\)

avec \(\bar S^{\ \ t}_{t-\tau }\) = la moyenne de S au cours des \(\tau\) dernières années

Il vient finalement

\[ T(t)\ -\ T_0\ =\ \lambda \ (\bar S^{\ \ t}_{t-\tau }\ -\ S_0) \tag{4.13}\]

Ce qui est tout à fait équivalent à l’Équation 4.2 introduite à la Section 4.2.

Remarquons au passage que \(T(t)\ -\ T_0\) représente une anomalie de température.

Dans le cas où \(\Delta S(t)\ =\ \Delta S_a\ \sin (2\pi \frac t{\mathit{per}})\), la solution de l’Équation 4.8 est la suivante :

\[ T(t)\ -\ T_0\ =\ \frac{\mathit{per}}{\sqrt{(4\pi ^2k^2\ +\ \mathit{per}^2)}}\ \lambda \ \Delta S_a\ \sin (2\pi \ \frac t{\mathit{per}}\ -\ \phi ) \tag{4.14}\]

En remplaçant \(per\) par \(2\pi/\omega\), l’Équation 4.14 s’écrit :

\[ T(t)\ -\ T_0\ =\ \frac 1{\sqrt{(1\ +\ (k\omega )^2)}}\ \lambda \ \Delta S_a\ \sin (\omega \ t\ -\ \phi ) \tag{4.15}\]

La variation de température est un signal sinusoïdal dont l’amplitude vaut

\[ \Delta T_a\ =\ \frac 1{\sqrt{(1\ +\ (k\omega )^2)}}\ \lambda \ \Delta S_a \tag{4.16}\]

et qui est déphasée d’un angle \(\phi\) par rapport au rayonnement solaire :

\[ \phi \ =\ \mathit{atan}(2\pi \frac k{\mathit{per}})\ =\ \mathit{atan}(k\omega ) \tag{4.17}\]

Pour des oscillations lentes, \(k \omega\) tend vers zéro, l’amplitude de la variation de température tend vers \(\lambda \Delta S_a\) et \(\phi\) 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 \(\phi\) tend vers \(\pi/2\).

4.4 Estimation des valeurs numériques des paramètres du modèle thermique

L’analyse TS.5 (voir Section 4.7.6) fournit directement la valeur de \(\lambda\) du modèle. C’est la valeur du coefficient nrl2TsiAvg qui vaut 1.66. Elle fournit également la valeur de \(\tau\) qui est égale au coefficient nrl2TsiAvg_tau, soit 106 ans. Une analyse de régression non linéaire fournit une valeur de k égale à 56.6 ans. C’est avec ces valeurs que la Figure 4.2 a été tracée.

Le coefficient k peut également être trouvé en résolvant l’Équation 4.15 par rapport à k :

\[ k\ =\ \frac{\sqrt{(\frac{\lambda \ \Delta S_a}{\Delta T_a})^2\ -\ 1}}{\omega } \tag{4.18}\]

L’examen des ondulations de la courbe rouge de la Figure 4.27 autour des années 1960 permet d’estimer pour un cycle solaire de 11 ans, \(\Delta S_a{\simeq}\ 0.4W/M^2\) et \(\Delta T_a{\simeq}\ 0.02{}^{\circ}C\).

Avec \(\lambda\)= 1.66 et \(\omega=2\pi/11\), il vient \(k\simeq 58\) ans. Les estimations basées sur la Figure 4.27 étant assez grossières, le hasard vient vraisemblablement à la rescousse pour donner une concordance aussi bonne.

Avec k=56.6 ans, le déphasage du cycle solaire par rapport au rayonnement est donné par l’Équation 4.17, et vaut 88.2 degrés.

Ceci confirme le décalage d’un quart de période entre la température et l’activité solaire (voir Figure 2.16) et valide donc également l’analogie de l’océan avec un filtre passe-bas « intégrateur ».

Avec \(\lambda=1.66\ °K/(W/m^2)\) et \(k=56.6 \ ans\) et un saut constant de rayonnement \(\Delta S=1\ W/m^2\), l’Équation 4.10 donne \(\Delta T=1.66\ °K\).

Pour un rayonnement sinusoïdal \(\Delta S_a=1\ W/m^2\) de période \(per\), avec les mêmes valeurs de \(\lambda\) et de \(k\), l’amplitude de la variation de température a été calculée pour différentes valeurs de \(per\). Voir Table 4.1.

Table 4.1: Amplitudes et déphasages de la variation de température en réponse à un rayonnement solaire périodique ayant une amplitude de 1 W/m².
\(per\)	1	11	100	200	500	1000	10000	year
\(\omega\)	6.2832	0.5712	0.0628	0.0314	0.0126	0.0063	0.000063	1/year
\(\Delta T_a\)	0.0047	0.0513	0.4494	0.8137	1.3527	1.5640	1.6600	°K
\(\phi_{deg}\)	89.8	88.2	74.3	60.6	35.4	19.6	0.2	angular °
\(\phi_{years}\)	0.25	2.70	20.64	33.69	49.20	54.38	56.60	year

L’amplitude de la variation de température est faible lorsque la période est courte, et tend vers la valeur correspondant à une variation constante du rayonnement solaire lorsque la période est très longue. Le déphasage passe progressivement de 90 degrés vers 0 degré. Pour une même variation de rayonnement solaire, la variation de température est environ 16 fois plus forte pour une période de 200 ans que pour une période de 11 ans. Les différences de température entre les valeurs max et min sont égales au double des valeurs mentionnées. Par exemple, pour un cycle de 200 ans, la variation totale de température est de 1.6 °K.

4.5 Analogie du modèle thermique avec un circuit électrique de type RC

L’Équation 4.7 est équivalente à celle qui décrit le fonctionnement d’un circuit électrique de type RC. Voir ici.

Figure 4.3: Circuit RC

La différence de potentiel aux bornes du condensateur est donnée par l’équation suivante :

\[ C\frac{\mathit{dV}_C}{\mathit{dt}}+\frac{V_C} R=V_{in} \tag{4.19}\]

L’équivalence des symboles des 2 équations est donnée par la table suivante :

Modèle thermique	Circuit RC
\(T\)	\(V_C\)
\(\lambda\)	\(R\)
\(C\)	\(C\)
\(S\)	\(V_{in}\)

4.6 Rayonnement solaire, températures et anomalies de température

On considère souvent que la température étant une grandeur intensive qui n’a qu’une signification thermodynamique locale, moyenner spatialement des températures locales n’a pas de sens, et qu’il faudrait par conséquent se limiter à ne faire que des analyses climatiques locales.

Nous allons voir que ce n’est pas le cas pour les températures de surface océanique : le modèle thermique décrit précédemment explique à la fois les températures et leurs anomalies.

Il se fait qu’en moyenne, à chaque latitude, on peut définir une anomalie de température locale qui est proportionnelle au rayonnement solaire incident local qui ne dépend également que de la latitude.

Ces anomalies de températures locales peuvent dès lors être moyennées spatialement pour en dériver une anomalie de température globale qui est proportionnelle au rayonnement solaire global au sommet de l’atmosphère qui n’est autre que l’irradiance totale solaire.

Les analyses suivantes sont basées sur la série hadisst du Met Office – Hadley Centre. Les données peuvent être téléchargées ici (« Hadisst data »).

Cette série reprend pour chaque mois depuis 1870 les températures de surface océanique pour un maillage en latitude et en longitude de 1° x 1° , soit 64800 points par mois, et plus de 100 millions de données au total.

4.6.1 Principes des calculs

Les températures de surface océaniques sont relativement homogènes à une latitude donnée, quelle que soit la saison. Il y a quand même plusieurs exceptions :

l’Atlantique Nord est plus chaud que le Pacifique Nord
au Sud de l’équateur, les côtes Est de l’Afrique et de l’Amérique du Sud sont plus chaudes que les côtes Ouest
les phénomènes El Niño – La Niña induisent des perturbations.

Il y a une variation saisonnière, mais aux niveaux de l’équateur et des cercles polaires, les températures sont assez stables. Voir Figure 4.4.

Figure 4.4: Températures de surface de la série hadisst pour les mois de mars, juin, septembre et décembre de l’année 2020. Les températures varient de -1.8 °C aux pôles à un peu plus de 30 °C à l’équateur.

Cette distribution des températures n’est pas sans rappeler celle du rayonnement solaire incident à la surface de la planète qui présente également un maximum au niveau de l’équateur, et des minima dans les zones polaires.

Le rayonnement solaire incident a été calculé pour chacun des 64800 points de la surface de la planète, toutes les 10 minutes pour une année complète, avec des moyennes par jour, par mois et pour l’année entière. On considérera que l’orbite de la planète est suffisamment stable entre 1870 et nos jours pour que la comparaison avec n’importe quelle période ou année de la série hadisst soit valide.

Il faut appliquer un flag ( 1=océan ou 0=terres) pour ne garder que les points correspondants aux océans. Ce flag est déduit de la série hadisst pour chaque mois. Dans la série hadisst, la banquise est en effet affectée d’une température conventionnelle égale à -1000 °C qui doit être ignorée pour ne pas fausser les calculs. La banquise évoluant au fil des saisons et des années, cette information n’est évidemment pas constante. Les terres quant à elles, ont une température égale à NA (not available) qui doit également être ignorée. Cela peut éventuellement être le cas aussi pour des données manquantes.

Pour chaque mois analysé, les températures et le rayonnement solaire incident ont été moyennés par latitude et comparés par analyse de régression.

La température moyenne globale a également été calculée pour chaque mois. Elle a fait l’objet de quelques analyses statistiques et été comparée avec la série mensuelle hadsst3.

4.6.2 Températures globales hadisst et comparaison avec les anomalies de température de la série hadsst3

La température moyenne globale de la série hadisst a été calculée en moyennant les 64800 mesures mensuelles pondérées par le cosinus de leur latitude. Son graphe est repris à la Figure 4.5 qui présente également celui de la série hadsst3 à titre de comparaison.

Figure 4.5: En haut, la température globale dérivée de la série hadisst. En bas, l’anomalie de température de la série hadsst3. La série hadisst présente une variabilité plus importante.

Les deux graphiques sont superposés dans la Figure 4.6, après ajout à la série hadsst3 de la température moyenne de la série hadisst au cours de la période de référence (1961 – 1990).

Figure 4.6: En rouge, la série hadisst. En vert, la série hadsst3 + la température moyenne de hadisst pendant la période de référence de hadsst3. En bleu, la température moyenne de hadisst pendant la période de référence de hadsst3. La concordance entre les 2 séries est très bonne jusqu’en 1980, mais un peu moins bonne par la suite.

Afin d’affiner la comparaison, les deux séries ont fait l’objet d’une décomposition additive en termes de tendance (trend), de variation saisonnière (seasonal) et le reste (remainder). Voir Figure 4.7 et Figure 4.8.

Figure 4.7: Décomposition de la série hadisst en composants de tendance (trend), de variation saisonnière (seasonal) et le reste (remainder) La variation saisonnière crête à crête est importante, environ 0.5 °C, soit à peu près 10 fois plus que celle de la série hadsst3. La variabilité augmente fortement après 1980, avec des sauts importants aux alentours de 2010. C’est difficile à expliquer par rapport à la variabilité de la série hadsst3 qui diminue au fil du temps.

Figure 4.8: Décomposition de la série hadsst3 en composants de tendance (trend), de variation saisonnière (seasonal) et le reste (remainder). Il subsiste une faible variation saisonnière crête à crête d’environ 0.05 °C, à peu près 10 fois plus faible que celle de la série hadisst. C’est tout à fait logique, parce que pour les anomalies de température il n’y a pas une seule température de référence, mais 12 (une par mois), ce qui ramène la variabilité saisonnière quasiment à 0. La variabilité diminue au fil du temps à mesure que le nombre de données augmente.

4.6.3 Le rayonnement solaire incident

La Terre tourne autour du Soleil en parcourant en 365.25 jours une orbite elliptique dont le soleil occupe un des foyers. La distance moyenne Terre – Soleil définit l’unité astronomique (1 UA = 1.496 108 km). Elle passe par un maximum au cours de l’été de l’hémisphère Nord, et par un minimum au cours de l’hiver de l’hémisphère Nord. Voir Figure 4.9.

Figure 4.9: Déclinaison du Soleil au cours de l’année. Voir chapitre 3 dans (Stine et Geyer 2001)

L’altitude aussi appelée élévation du soleil en un point de la surface terrestre est l’angle entre la direction du soleil et le plan horizontal local. Le rayonnement solaire en ce point est proportionnel au sinus de l’altitude, en tenant compte du fait que le soleil est sous l’horizon lorsque l’altitude est négative, et que dans ce cas, le rayonnement incident est bien évidemment nul. Voir Figure 4.10.

Figure 4.10: Système de coordonnées local à la surface de la Terre montrant l’azimuth du Soleil (angle A), l’altitude (ou élévation) du Soleil (angle \(\alpha\)), et l’angle du Soleil par rapport au zénith (\(\theta_z\)). La direction du Soleil est représentée par le vecteur S. Le point Q représente un observateur à la surface de la Terre. Voir chapitre 3 dans (Stine et Geyer 2001)

Quelques exemples de trajectoires du Soleil sur la sphère céleste sont affichés à la Figure 4.11.

Figure 4.11: La trajectoire du Soleil sur la sphère céleste au cours de la journée pour un observateur à 56° de latitude nord. La trajectoire du Soleil change avec sa déclinaison au cours de l’année. Les intersections des lignes avec l’axe horizontal montrent des azimuts en degrés par rapport au nord où le soleil se lève et se couche aux solstices d’été et d’hiver. Par Deditos — Travail personnel, CC BY-SA 4-0, https://commons.wikimedia.org/w/index.php?curid=37904452

L’altitude du Soleil dépend de sa direction définie par la droite joignant les centres de la Terre et du soleil, ainsi que de la position de l’observateur à la surface terrestre et du temps. Voir Figure 4.12.

Figure 4.12: Système de coordonnées au centre de la Terre montrant la déclinaison du Soleil (angle \(\delta\)) et l’angle horaire (\(\omega\)) qui est égal à la différence entre la longitude de l’observateur Q et la longitude du méridien contenant le vecteur de la direction du Solaire indiquée par le vecteur S’. Voir chapitre 3 dans (Stine et Geyer 2001)

Figure 4.13: Vue composite des deux schémas précédents. L’angle ϕ représente la latitude de l’observateur Q. Les directions du Soleil S et S’ peuvent être considérées parallèles parce que le rayon de la Terre est négligeable par rapport à la distance Terre – Soleil. Voir chapitre 3 dans (Stine et Geyer 2001)

Après quelques manipulations mathématiques, on trouve que le sinus de l’altitude du soleil peut être calculé par l’Équation 4.20, et le rayonnement incident par l’Équation 4.21

\[ \sin (\alpha )\ =\ \sin (\delta )\ \sin (\phi )\ +\ \cos (\delta )\ \cos (\omega )\ \cos (\phi ) \tag{4.20}\]

\[ R_u^{\mathit{inc}}\ =\ \sin (\alpha )\ /\ (d_{\mathit{TS}}/\bar d_{\mathit{TS}})^2\ =\ \sin (\alpha )\ /\ (d_{\mathit{TS}}^{\mathit{UA}})^2 \tag{4.21}\]

Les symboles repris dans ces formules, sont décrits ci-après :

\(\alpha\)	L’altitude du Soleil en un point de la surface terrestre à un instant donné.
\(\delta\)	La déclinaison du Soleil à un instant donné.
\(\phi\)	La latitude du point d’observation sur la surface terrestre.
\(\omega\)	L’angle horaire, qui est la différence entre la longitude du point d’observation et la longitude du méridien contenant le vecteur de la direction du Soleil à un instant donné.
\(R_u^{inc}\)	Le rayonnement incident par unité de surface en un point de la surface terrestre à un instant donné, pour un rayonnement unitaire au sommet de l’atmosphère.
\(d_{TS}\)	La distance Terre - Soleil à un instant donné.
\(\bar d_{TS}\)	La distance Terre – Soleil moyenne = 1 unité astronomique
\(d_{TS}^{UA}\)	La distance Terre – Soleil à un instant donné, exprimée en unité astronomique

La déclinaison du soleil, la longitude du vecteur de sa direction et la distance Terre-Soleil dépendent du temps. Elles sont calculées en utilisant des tables d’éphémérides et des programmes publiés par le Jet Propulsion Laboratory. La matière est très touffue, et il n’est pas facile de trouver son chemin. Il existe heureusement des packages en langage Python (SpiceyPy et Skyfield qui facilitent quelque peu les calculs.

Quelques résultats de calcul sont repris à la Figure 4.14.

Figure 4.14: Le rayonnement incident moyen pour quelques dates et l’année entière 2020. Il est exprimé en fraction du rayonnement au sommet de l’atmosphère. Les images de la colonne de gauche représentent la distribution spatiale du rayonnement en latitude et longitude. Les tracés de la colonne de droite reprennent la distribution du rayonnement à 0 degré de longitude. De haut en bas, les journées de l’équinoxe de printemps, du solstice d’été, de l’équinoxe d’automne, du solstice d’hiver, et pour terminer, l’année complète 2020. Sur une journée entière, le rayonnement incident moyen est quasiment indépendant de la longitude.

4.6.4 Comparaison des températures hadisst et du rayonnement solaire incident

Pour tous les mois de l’année 2020, la température de surface océanique hadisst a été moyennée par latitude. Le rayonnement incident a également été moyenné par latitude, en se limitant aux points considérés comme « océan » dans la série hadisst. Une analyse de régression a chaque fois été menée entre la température hadisst et le rayonnement incident. Deux séries de calculs ont été effectuées. Dans la première série, le rayonnement incident moyen du mois a été utilisé (Figure 4.15). Dans la seconde série, le rayonnement incident moyen de l’année a été utilisé (Figure 4.16).

Finalement, la régression de la température océanique moyenne annuelle par latitude de hadisst par le rayonnement incident moyen annuel a également été calculée (Figure 4.17). La corrélation est quasiment parfaite, et permet de conclure que les observations montrent que la répartition de la température océanique découle de celle du rayonnement incident annuel. Ceci a d’importantes conséquences quant à l’usage de moyennes spatiales et temporelles des températures océaniques dans les études climatiques. Ceci est discuté dans la section suivante.

Figure 4.17: Régression de la température moyenne annuelle par latitude de hadisst par le rayonnement incident moyen annuel. La corrélation est excellente.

4.6.5 Réhabilitation des anomalies de température

L’examen de la série de température hadisst révèle qu’à une latitude donnée, les températures de surface océaniques varient peu à latitude constante, et qu’en moyenne annuelle, en formalisant l’analyse de régression ( Figure 4.17), on a :

\[ T_{\mathit{ANN}}^{\mathit{LOC}}(\mathit{lat})\ -\ T_{\mathit{REF0}}\ =\ \alpha \ U_{\mathit{ANN}}^{\mathit{LOC}}(\mathit{lat}) \tag{4.22}\]

Dans cette formule,

\(lat\)	= la latitude.
\(T_{ANN}^{LOC}(lat)\)	= la température annuelle locale moyenne à la latitude considérée.
\(T_{REF0}\)	= une température de référence.
\(\alpha\)	= une constante.
\(U_{ANN}^{LOC}(lat)\)	= le rayonnement incident annuel moyen à la latitude considérée, pour une irradiance totale solaire unitaire au sommet de l’atmosphère.

En intégrant cette expression sur la surface des océans, et en divisant par la surface totale océanique, on obtient une anomalie globale moyenne :

\[ \frac{\iint _{\mathit{ocean}}[T_{\mathit{ANN}}^{\mathit{LOC}}(\mathit{lat})\ -\ T_{\mathit{REF0}}]\ \mathit{dS}}{\iint _{\mathit{ocean}}\mathit{dS}}\ =\ \alpha \ \frac{\iint _{\mathit{ocean}}U_{\mathit{ANN}}^{\mathit{LOC}}(\mathit{lat})\ \mathit{dS}}{\iint _{\mathit{ocean}}\mathit{dS}} \tag{4.23}\]

En considérant une Terre sphérique de rayon R, l’élément de surface dS en coordonnées sphériques est égal à :

\[ \mathit{dS}\ =R^2\cos (\mathit{lat})\ \mathit{dlat}\ \mathit{dlon} \tag{4.24}\]

Pour calculer la moyenne d’une grandeur sur une sphère, il faut la pondérer par le cosinus de sa latitude.

L’Équation 4.23 peut s’écrire

\[ \Delta T_{\mathit{ANN}}^{\mathit{GLO}}(T_{\mathit{REF0}})\ =\ T_{\mathit{ANN}}^{\mathit{GLO}}\ -\ T_{\mathit{REF0}}\ =\ \alpha \ U_{\mathit{ANN}}^{\mathit{GLO}} \tag{4.25}\]

Dans cette formule,

\(\Delta T_{\mathit{ANN}}^{\mathit{GLO}}(T_{\mathit{REF0}})\)	= l’anomalie globale de température annuelle moyenne pour une température de référence \(T_{REF0}\)
\(T_{\mathit{ANN}}^{\mathit{GLO}}\ =\ \frac{\iint _{\mathit{ocean}}T_{\mathit{ANN}}^{\mathit{LOC}}(\mathit{lat})\ \mathit{dS}}{\iint _{\mathit{ocean}}\mathit{dS}}\)	= la température globale annuelle moyenne
\(U_{\mathit{ANN}}^{\mathit{GLO}}\ =\ \frac{\iint _{\mathit{ocean}}U_{\mathit{ANN}}^{\mathit{LOC}}(\mathit{lat})\ \mathit{dS}}{\iint _{\mathit{ocean}}\mathit{dS}}\)	= le rayonnement global incident annuel moyen. Ce terme ne dépend que de la géométrie de la planète et de son orbite autour du Soleil.

Moyennant un choix adéquat de sa température de référence, l’anomalie globale de température annuelle moyenne est proportionnelle au rayonnement incident annuel moyen.

Bien que la température de surface océanique soit une grandeur intensive, il est tout à fait légitime d’en faire des moyennes spatiales pour calculer une anomalie globale, parce qu’en procédant ainsi, on ne fait que moyenner des rayonnements, ce qui est tout à fait correct. C’est l’Équation 4.22 qui justifie cette pratique.

L’Équation 4.25 s’adapte facilement si une anomalie de température utilise une autre température de référence \(T_{REF1}\) :

\[ \Delta T_{\mathit{ANN}}^{\mathit{GLO}}(T_{\mathit{REF1}})\ =\ T_{\mathit{ANN}}^{\mathit{GLO}}\ -\ T_{\mathit{REF1}}\ =\ \alpha \ (U_{\mathit{ANN}}^{\mathit{GLO}}\ -\ U_{\mathit{REF}}) \tag{4.26}\]

Dans cette relation,

\(U_{\mathit{REF}}\ =\ \frac{(T_{\mathit{REF1}}\ -\ T_{\mathit{REF0}})}{\alpha }\) = un rayonnement unitaire de référence.

Pour passer du rayonnement unitaire au rayonnement au sommet de l’atmosphère, il suffira de modifier la constante de proportionnalité . On aura alors

\[ \Delta T_{\mathit{ANN}}^{\mathit{GLO}}(T_{\mathit{REF1}})\ =\ \ a\ (S_{\mathit{ANN}}\ -\ S_{\mathit{REF}}) \tag{4.27}\]

Dans cette formule,

\(a\)	= une constante
\(S_{ANN}\)	= le rayonnement annuel au sommet de l’atmosphère
\(S_{REF}\)	= un rayonnement au sommet de l’atmosphère de référence

l’Équation 4.27 est valable lorsque le rayonnement solaire au sommet de l’atmosphère est constant. Remarquons que dans ce cas, \(S_{ANN}\) est égal à sa valeur moyenne au cours des \(\tau\) dernières années \(\bar S^{\ \ t}_{t-\tau }\), quelle que soit la valeur de \(\tau\).

L’Équation 4.27 se transforme alors en

\[ \Delta T^{\mathit{GLO}}(t,T_{\mathit{REF}})\ =\ T^{\mathit{GLO}}(t)\ -\ T_{\mathit{REF}}\ =\ a\ (\bar S^{\ \ t}_{t-\tau }\ -\ S_{\mathit{REF}}) \tag{4.28}\]

On retrouve ainsi la solution approchée de l’équation différentielle du modèle thermique simple (Équation 4.13), telle que décrite à la Section 4.3. Cette solution est valable lorsque le rayonnement solaire au sommet de l’atmosphère n’est pas constant, à condition de faire un choix approprié de \(\tau\).

4.6.6 Synthèse

A chaque latitude, et à chaque instant, les moyennes annuelles des températures de surface océaniques sont très bien modélisées par le rayonnement solaire incident.

Il est tout à fait légitime d’effectuer des moyennes spatiales et temporelles des températures de surface océaniques.

Travailler avec des températures ou avec leurs anomalies est équivalent, à une translation près.

Lorsqu’on détermine une tendance temporelle, on risque d’introduire un biais si la durée de la période utilisée pour déterminer la tendance n’est pas sensiblement plus longue que la durée du cycle, mais le biais sera le même avec la température ou avec son anomalie.

Indépendamment de ce qui précède, il serait préférable de travailler avec la température plutôt qu’avec son anomalie, parce que dans ce dernier cas, on perd de l’information :

La variabilité saisonnière annuelle disparaît.
La variabilité spatiale disparaît également. Sur base des anomalies de température, il est impossible de voir que le rayonnement incident est responsable de la distribution de la température par latitude.

4.7 Analyses effectuées au moyen du modèle dynamique

4.7.1 Sélection des variables indépendantes dans les analyses de régression TS (température - soleil)

Les matrices de corrélation de la Section 2.6.1.1 montrent que l’activité solaire moyennée nrl2TsiAvg est meilleure que la série histTsiAvg et que la série sunspots. Les tests de causalité de la Section 2.6.2.2 montrent une causalité bidirectionnelle entre la température de surface océanique (hadsst3) et le logCO2, ainsi qu’une causalité unidirectionnelle de l’activité solaire moyennée vers hadsst3. Un terme harmonique (harm), comme dans l’Équation 4.3 a également été pris en compte.

Toutes les combinaisons possibles des variables indépendantes harm, nrl2TsiAvg, et logCO2 ont été utilisées pour modéliser la variable dépendante hadsst3 par analyse de régression.

Les résultats de ces analyses sont résumées dans la Table 4.2.

Table 4.2: Capacités prédictives de toutes les combinaisons possibles d’un terme harmonique (harm), de l’ITS moyenne (nrl2TsiAvg) et du log(CO2).
harm	nrl2TsiAvg	logCo2	CV	AIC	AICc	BIC	AdjR2	%SE harm	%SE nrl2TsiAvg	%SE logCo2
x	x		0.0157	-8452.2	-8452.1	-8452.2	0.793	4.05	1.24
x	x	x	0.0157	-8450.8	-8450.8	-8450.8	0.793	4.04	5.16	80.15
x		x	0.0165	-8346.5	-8346.4	-8346.5	0.782	3.62		1.17
	x		0.0183	-8131.8	-8131.8	-8131.8	0.758		1.25
	x	x	0.0184	-8129.8	-8129.8	-8129.8	0.758		4.56	665.40
		x	0.0227	-7700.0	-7700.0	-7700.0	0.701			1.45
x			0.0627	-5631.2	-5631.2	-5631.2	0.173	4.85

La table est triée par capacité prédictive décroissante. Les trois dernières colonnes reprennent le pourcentage de l’erreur standard sur le coefficient de chaque terme de régression utilisé. La meilleure solution est donnée par la combinaison d’un terme harmonique et de l’ITS moyenne.

L’ajout d’un terme log(CO2) détériore la situation. Seul, ou en combinaison avec un terme harmonique, le log(CO2) donne un résultat correct, mais moins bon. Lorsqu’il est combiné avec l’ITS moyenne, l’erreur sur son coefficient devient très importante, alors que les erreurs sur les coefficients du terme harmonique et de l’ITS moyenne sont peu affectées.

De plus nous avons montré à la Section 2.6.3.3 que l’évolution du CO₂ suit celle de la température avec un décalage de 8 mois.

On peut en conclure que rien ne justifie l’introduction d’un terme log(CO₂) dans les analyses de régression TS.

Nous verrons plus tard au Chapitre 6, qu’en dehors des émissions anthropiques, le CO₂ naturel dépend de la température suivant l’équilibre de Henry – van ’t Hoff. Suivant la même relation d’équilibre, le CO₂ anthropique a un petit effet sur la température, mais il est suffisamment faible pour qu’on puisse le négliger dans les analyses.

Plusieurs analyses sont détaillées ci-après.

Les analyses TS.1, TS.2 et TS.3, basées uniquement sur l’activité solaire moyenne, n’ont été reprises que pour montrer de manière détaillée que des 3 variables sunspotsAvg, histTsiAvg et nrl2TsiAvg, c’est effectivement cette dernière qui est la meilleure, ce qui confirme ce qui a été trouvé à la Section 2.6.1.2.

4.7.2 Analyse TS.1: taches solaires et température

Cette analyse utilise les taches solaires (sunspots) comme mesure du rayonnement solaire. Voir Figure 4.18 et Figure 4.19

Figure 4.18: En bleu, l’anomalie de température (variable dépendante). En vert, la seule variable indépendante basée sur le nombre moyen de taches solaires, masquée par le modèle (en rouge). La corrélation est bonne. Le problème des taches solaires, c’est que les minima sont tous au même niveau (zéro). Cela limite leur potentiel explicatif. Le nombre moyen de taches solaires est passé par un maximum et a commencé à diminuer.

Table 4.3: TS.1 Regression Summary
Parameter	Value	95% interval
sunspotsAvg	0.0387408	3.32%
sunspotsAvg_a	82.8024000	0.25%
sunspotsAvg_tau	110.2870000	0.31%

Figure 4.19: Les résidus de l’analyse TS.1 Il reste une tendance et un biais assez important au début et à la fin. On décèle également une périodicité. Une analyse détaillée de ces résidus a été faite à la Section 4.7.9.2.1

4.7.3 Analyse TS.2: irradiance totale solaire et température

Cette analyse utilise la reconstruction historical_tsi comme mesure du rayonnement solaire. Voir Figure 4.20 et Figure 4.21

Figure 4.20: En bleu, l’anomalie de température (variable dépendante). En vert, la seule variable indépendante basée sur l’irradiance totale solaire, masquée par le modèle (en rouge). La corrélation est similaire à celle de l’analyse TS.1

Table 4.4: TS.2 Regression Summary
Parameter	Value	95% interval
histTsiAvg	4.99454	3.39%
histTsiAvg_ref	1361.17000	0.00%
histTsiAvg_tau	109.77500	0.39%

Figure 4.21: Les résidus de l’analyse TS.2 Il y a une tendance et un biais assez important au début et à la fin. Une analyse détaillée de ces résidus a été faite à la Section 4.7.9.2.2

4.7.4 Analyse TS.3: irradiance totale solaire et température

Cette analyse utilise la reconstruction nrl2_tsi comme mesure du rayonnement solaire. Voir Figure 4.22 et Figure 4.23.

Figure 4.22: En bleu, l’anomalie de température (variable dépendante). En vert, la seule variable indépendante basée sur le nombre moyen de taches solaires, masquée par le modèle (en rouge). La corrélation est sensiblement meilleure que dans les analyses TS.1 et TS.2

Table 4.5: TS.3 Regression Summary
Parameter	Value	95% interval
nrl2TsiAvg	1.83245	2.77%
nrl2TsiAvg_ref	1360.74000	0.00%
nrl2TsiAvg_tau	111.97500	1.39%

Figure 4.23: Les résidus de l’analyse TS.3 Il n’y a pas de tendance dans les résidus, mais on peut déceler un terme harmonique. Une analyse détaillée de ces résidus a été faite à la Section 4.7.9.2.3

L’analyse spectrale de la série hadsst3gm après detrending cubique confirme la présence d’un pic important et de quelques pics secondaires de moindre amplitude. Voir Figure 4.24

Figure 4.24: Périodogramme de la série hadsst3gm

Les périodes des 4 premiers pics principaux sont égales à 66,7 ans, 21,4 ans, 9,0 ans et 3,6 ans. La figure 4 dans l’article de (Scafetta 2021) contient des pics très voisins.

Sur cette base, un terme harmonique a été ajouté dans l’analyse suivante conformément à l’Équation 4.3.

4.7.5 Analyse TS.4: TS.3 + un terme harmonique

Cette analyse utilise la reconstruction nrl2_tsi et un terme harmonique dont la période calculée par analyse de régression multiple non linéaire vaut 66.28 ans. Cette valeur est très voisine de la période de 66.7 ans trouvée dans l’analyse spectrale précédente. Voir Figure 4.25 et Figure 4.26

Figure 4.25: En bleu, l’anomalie de température (variable dépendante). En vert, le terme harmonique avec une amplitude d’environ un dixième de degré et une période d’environ 66 ans. Cette période correspond très bien au premier pic de l’analyse spectrale montrée à la Fig. 27. En orange, le terme correspondant au rayonnement solaire. En rouge, le modèle (somme des 2 termes précédents). Le terme basé sur le rayonnement solaire moyen et le terme harmonique viennent de passer par un maximum et commencent à diminuer maintenant. La corrélation est meilleure que dans l’analyse TS.3

Table 4.6: TS.4 Regression Summary
Parameter	Value	95% interval
harm	9.70287e-02	9.27%
harm_phase	1.88140e+03	0.07%
harm_period	6.62807e+01	1.99%
nrl2TsiAvg	1.56256e+00	3.07%
nrl2TsiAvg_ref	1.36076e+03	0.00%
nrl2TsiAvg_tau	1.01040e+02	1.77%

Figure 4.26: Les résidus de l’analyse TS.4 Il n’y a plus aucune tendance dans les résidus. Les résidus sont plus importants au début parce que le maillage des mesures était beaucoup plus faible dans les premières décennies des relevés de température. Une analyse détaillée de ces résidus a été faite à la Section 4.7.9.2.4

4.7.6 Analyse TS.5: TS.4 + lissage de la température

C’est la même analyse que l’analyse TS.4, sauf que le relevé des températures a fait l’objet d’un lissage par spline de degré 3 pour éliminer les épisodes El Niño/La Niña et l’impact des éruptions volcaniques. Ces éléments qui « perturbent » l’analyse n’ont pas d’influence sur le climat parce qu’ils sont plus ou moins cycliques ou ont une brève durée. Le lissage affecte peu les paramètres du modèle, mais révèle la véritable corrélation.

Figure 4.27: Les intervalles de confiance des paramètres sont beaucoup plus faibles. Le coefficient de corrélation est impressionnant. Le rayonnement solaire et un terme harmonique modélisent parfaitement le relevé des températures de la série hadsst3 lissée.

Table 4.7: TS.5 Regression Summary
Parameter	Value	95% interval
harm	8.85452e-02	2.72%
harm_phase	1.88116e+03	0.02%
harm_period	6.67980e+01	0.58%
nrl2TsiAvg	1.66317e+00	0.86%
nrl2TsiAvg_ref	1.36075e+03	0.00%
nrl2TsiAvg_tau	1.05546e+02	0.45%

Figure 4.28: Les résidus de l’analyse TS.5. Notez l’échelle des résidus qui sont compris entre plus ou moins 0.07 °C. Il y a apparemment encore un terme harmonique, mais si on l’ajoute dans l’analyse, on trouve un terme d’amplitude inférieure à 0.01 °C avec une période d’environ 20 ans qui correspond assez bien au deuxième pic de l’analyse spectrale montrée à la Figure 4.24. Ce terme ne fournirait qu’une amélioration dérisoire du coefficient de corrélation qui passerait à 0.982.

4.7.7 Analyse TS.6: une projection rétrospective

C’est la même analyse que l’analyse TS.4, mais en se plaçant rétrospectivement en l’an 2000. Cette analyse a été menée afin d’évaluer les capacités prédictives du modèle et de voir comment l’anomalie de température évolue dans le futur. Les mesures antérieures à 2000 servent à calibrer le modèle, tandis que les mesures plus récentes servent à le valider.

Il est difficile d’utiliser le modèle pour faire des prévisions parce qu’il est basé sur l’irradiance totale solaire qui est encore plus difficile à prédire que les taches solaires.

L’option choisie considère que le rayonnement solaire futur est constant et égal à la valeur du rayonnement solaire moyen en 2000 (modèle prédictif simple dit “naïf”), soit 1360.9 W/m2 (Voir Figure 4.29). Cette hypothèse est en dessous de la réalité puisque le rayonnement solaire moyen a continué à augmenter après l’an 2000, avant de commencer à décroître légèrement (Cfr Figure 2.5 ).

Figure 4.29: Le rayonnement solaire utilisé dans l’analyse TS.6 En trait plein, la valeur instantanée, et en pointillé la valeur mobile des 110 dernières années. Bien que la projection du rayonnement solaire instantané soit constante (1360.9 W/m2), le rayonnement solaire moyen continue d’augmenter environ jusqu’en 2040 parce que pendant environ 40 ans les années qui sortent du tampon ont en moyenne un rayonnement instantané plus faible. Par la suite, c’est le contraire, et le rayonnement solaire moyen diminue progressivement pour s’établir à la valeur constante en 2110.

Figure 4.30: A comparer avec la Figure 4.25. En bleu (trait épais), la série hadsst3 entre 1850 et 2000. C’est la variable dépendante de l’analyse de régression. En bleu (trait fin) les températures réellement relevées après 2000. En vert, le terme harmonique (trait plein) et sa projection (trait pointillé). En orange, le terme correspondant au rayonnement solaire (trait plein) et sa projection (trait pointillé). En rouge, le modèle (trait plein) et sa projection (trait pointillé). Le coefficient de corrélation est moins bon (0.657) que dans l’analyse TS.4 (0.794) mais reste tout à fait correct. Les valeurs des paramètres et leurs intervalles de confiance sont voisins de ceux de l’analyse TS.4. Une projection qui aurait été effectuée en 2000 correspondrait très bien aux observations pendant 20 ans. L’explication est assez simple et parfaitement logique. Dans cette analyse, l’inertie thermique calculée vaut environ 100 ans. Au fil du temps, 20% d’anciennes données sont progressivement remplacés par 20% de nouvelles données. Après 20 ans, 80% des données sont inchangées. La longue inertie thermique permet de faire des analyses de tendance au-delà de quelques années. C’est également ce qui est mentionné dans (Abdussamatov 2013).

Table 4.8: TS.6 Regression Summary
Parameter	Value	95% interval
harm	0.107502	10.31%
harm_phase	1881.760000	0.06%
harm_period	65.497600	2.04%
nrl2TsiAvgProj	1.625270	3.78%
nrl2TsiAvgProj_ref	1360.750000	0.00%
nrl2TsiAvgProj_tau	100.235000	2.01%

Figure 4.31: Les résidus de l’analyse TS.6 sont semblables à ceux de l’analyse TS.4

Ce résultat valide le modèle.

4.7.8 Analyse TS.7: sortie du petit âge glaciaire

C’est exactement la même analyse que l’analyse TS.5. Les données de la série nrl2 remontent jusqu’en 1610. En tenant compte de l’inertie thermique de 105 ans, l’anomalie de température peut être calculée à partir de 1715.

Figure 4.32: Ceci donne une idée de la sortie du petit âge glaciaire. On retrouve le minimum de Maunder entre 1700 et 1750, ainsi que le minimum de Dalton entre 1800 et 1850.

Ce résultat constitue une autre validation du modèle.

4.7.9 Analyses complémentaires

4.7.9.1 Mesures de la capacité prédictive des modèles

Le coefficient de détermination R², historiquement le plus répandu, n’est pas la seule mesure utilisée pour comparer la capacité prédictive des modèles. Les statistiques CV, AIC, AICc, BIC sont également couramment utilisées. Un modèle est d’autant meilleur que son R² est grand et que ses statistiques CV, AIC, AICc et BIC sont petites.

Ces mesures permettent de sélectionner le modèle qui a la meilleure capacité prédictive. Elles sont expliquées au chapitre 5 dans (Hyndman et Athanasopoulos 2018).

Elles ont été résumées dans la Table Table 4.9.

Table 4.9: Models predictive capacity
Anal	CV	AIC	AICc	BIC	AdjR2	BG Test
TS.1	0.0276524	-6929.969	-6929.957	-6913.270	0.6441065	❌
TS.2	0.0277407	-7246.240	-7246.228	-7229.404	0.6239604	❌
TS.3	0.0183328	-8131.797	-8131.785	-8114.944	0.7580459	❌
TS.4	0.0156629	-8452.167	-8452.147	-8429.696	0.7934072	❌

TS.5	0.0011843	-13703.430	-13703.410	-13680.960	0.9810580	❌

TS.6	0.0162153	-7417.273	-7417.251	-7395.291	0.6561721	❌

Des modèles TS.1 à TS.4, il apparaît clairement que celui de l’analyse TS.4 est le meilleur, quelle que soit la statistique utilisée pour effectuer la comparaison.

Le modèle TS.5 ne peut être directement comparé avec les modèles précédents, parce qu’il a été calculé avec une anomalie de température lissée.

Le modèle TS6, qui est le modèle TS4 calculé avec des données arrêtées en l’an 2000, est repris à titre indicatif.

Les résidus des analyses échouent au test de Breusch-Godfrey (colonne BG Test). Les résidus sont analysés en détail à la Section 4.7.9.2.

4.7.9.2 Résidus des analyses de régression

Les résidus d’une analyse de régression sont égaux à la différence entre les valeurs observées et les valeurs modélisées. La moyenne des résidus doit être nulle, et ils ne devraient pas présenter d’auto-corrélation.

Les séries utilisées dans les analyses de régression présentent souvent une tendance et de l’auto-corrélation dans leurs résidus. Cela ne signifie pas que les prédictions de ces modèles soient biaisées, mais ces prédictions pourraient être améliorées en exploitant l’information contenue dans l’auto-corrélation.

Pour apprécier la qualité des résidus, on les compare à un bruit blanc qui ne contient aucune information et aucune auto-corrélation.

95 % des coefficients du graphe d’auto-corrélation des résidus doivent être situés dans l’intervalle \(\pm 2/\sqrt{(T)}\) où T représente le nombre d’observations. On peut également effectuer un test de type portemanteau (Breusch-Godfrey ou Ljung-Box) pour évaluer l’absence d’auto-corrélation dans les premiers résidus des analyses.

En cas d’auto-corrélation dans les résidus, on peut recourir à des modèles de régression dynamiques beaucoup plus sophistiqués, qui combinent des modèles de régression classiques et des modèles de type ARIMA, en espérant obtenir des résidus proches d’un bruit blanc qui satisfont aux tests de type portemanteau. Ces techniques de calcul sont décrites aux chapitres 8 et 9 dans (Hyndman et Athanasopoulos 2018).

Le test de Breusch-Godfrey est préféré pour les modèles de régression classiques, et celui de Ljung-Box pour les modèles de régression dynamiques. On considère qu’il y a absence d’auto-corrélation lorsque la probabilité de l’hypothèse nulle du test est supérieure à 0.05.

Tous les modèles de régression classique TS.1 à TS.4 présentent de l’auto-corrélation, et échouent au test de type portemanteau. Ils ont dès lors été complétés par un modèle de type ARIMA. Les projections de ces modèles ont été calculées en supposant que pour les 200 mois suivants, les valeurs des sunspots, irradiances totales solaires et du terme harmonique éventuel sont égales à leurs dernières observations (« naïve forecast »).

Les graphiques des analyses reprennent

en cyan, la modélisation par régression classique.
en rouge, la variable dépendante (l’anomalie de température hadsst3).
en noir, la modélisation « régression + ARIMA », qui se confond quasiment avec la courbe rouge, tellement l’ajustement est élevé.
en bleu, la projection « régression + ARIMA », avec ses intervalles de confiance à 80 % (en bleu foncé) et à 95 % en bleu clair.

Les capacités prédictives des modèles avec ARIMA sont résumées dans la Table Table 4.10.

Table 4.10: Models predictive capacity (ARIMA)
Anal	AIC	AICc	BIC	LJ Test
TS.1	-4529.11	-4529.06	-4495.71	❌
TS.2	-4588.48	-4588.45	-4560.42	❌
TS.3	-4661.13	-4660.95	-4588.10	✓
TS.4	-4679.82	-4679.64	-4606.79	✓

Il apparaît clairement que le modèle TS.4 est le meilleur, quelle que soit la statistique utilisée pour effectuer la comparaison.

Le détail de ces analyses est commenté ci-dessous.

4.7.9.2.1 Analyse résidus TS.1

4.7.9.2.1.1 Sans ARIMA

Figure 4.33: Les résidus de l’analyse TS1.

## 
##  Breusch-Godfrey test for serial correlation of order up to 10
## 
## data:  Residuals
## LM test = 1539.1, df = 10, p-value < 2.2e-16

Il y a une tendance dans les résidus qui sont dissymétriques. Le graphe d’auto-corrélation se stabilise à une valeur trop élevée pour satisfaire au test portemanteau.

4.7.9.2.1.2 Avec ARIMA

## regr_ar = auto.arima(df$hadsst3gm, xreg = df$sunspotsAvg)
##
## Series: df$hadsst3gm 
## Regression with ARIMA(2,1,1) errors 
## 
## Coefficients:
##          ar1     ar2      ma1  drift    xreg
##       0.5649  0.2677  -0.9891  3e-04  0.5928
## s.e.  0.0225  0.0224   0.0042  1e-04  0.1308
## 
## sigma^2 estimated as 0.005585:  log likelihood=2270.55
## AIC=-4529.11   AICc=-4529.06   BIC=-4495.71
## 
## Training set error measures:
##                        ME       RMSE       MAE MPE MAPE      MASE         ACF1
## Training set 5.346618e-06 0.07461492 0.0550197 Inf  Inf 0.9295593 -0.002741821

Figure 4.34: Les résidus de l’analyse TS1 avec ARIMA.

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(2,1,1) errors
## Q* = 87.45, df = 25, p-value = 7.426e-09
## 
## Model df: 5.   Total lags used: 30

Il n’y a pas de tendance dans les résidus qui sont quasiment symétriques. Les coefficients d’auto-corrélation sont légèrement trop élevés pour pour être validés par le test portemanteau. Les intervalles de confiance des projections augmentent avec le temps. Le coefficient sur le nombre moyen de taches solaires (sunspotsAvg) est quasiment égal à 0.5928, ce qui invalide l’analyse correspondante sans ARIMA.

Figure 4.35: Projections de l’analyse TS1 avec ARIMA.

4.7.9.2.2 Analyse résidus TS.2

4.7.9.2.2.1 Sans ARIMA

Figure 4.36: Les résidus de l’analyse TS2.

## 
##  Breusch-Godfrey test for serial correlation of order up to 10
## 
## data:  Residuals
## LM test = 1583.1, df = 10, p-value < 2.2e-16

Il y a une tendance dans les résidus qui sont dissymétriques. Le graphe d’auto-corrélation se stabilise à une valeur trop élevée pour satisfaire au test portemanteau.

4.7.9.2.2.2 Avec ARIMA

## regr_ar = auto.arima(df$hadsst3gm, xreg = df$histTsiAvg)
##
## Series: df$hadsst3gm 
## Regression with ARIMA(1,1,2) errors 
## 
## Coefficients:
##          ar1      ma1     ma2    xreg
##       0.8948  -1.3227  0.3300  0.6015
## s.e.  0.0144   0.0273  0.0265  0.1198
## 
## sigma^2 estimated as 0.006025:  log likelihood=2299.24
## AIC=-4588.48   AICc=-4588.45   BIC=-4560.42
## 
## Training set error measures:
##                       ME       RMSE        MAE MPE MAPE      MASE         ACF1
## Training set 0.002490818 0.07752253 0.05691908 Inf  Inf 0.9259642 -0.005702849

Figure 4.37: Les résidus de l’analyse TS2 avec ARIMA.

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(1,1,2) errors
## Q* = 81.36, df = 26, p-value = 1.292e-07
## 
## Model df: 4.   Total lags used: 30

Il n’y a pas de tendance dans les résidus qui sont quasiment symétriques. Les coefficients d’auto-corrélation sont légèrement trop élevés pour être validés par le test portemanteau. Les intervalles de confiance des projections augmentent avec le temps. Le coefficient sur le nombre moyen de taches solaires (histTsiAvg) est égal à 0.6015, ce qui invalide l’analyse correspondante sans ARIMA.

Figure 4.38: Projections de l’analyse TS2 avec ARIMA.

4.7.9.2.3 Analyse résidus TS.3

4.7.9.2.3.1 Sans ARIMA

Figure 4.39: Les résidus de l’analyse TS3.

## 
##  Breusch-Godfrey test for serial correlation of order up to 10
## 
## data:  Residuals
## LM test = 1377.9, df = 10, p-value < 2.2e-16

Il n’y a pas de tendance dans les résidus qui sont symétriques. Le graphe d’auto-corrélation se stabilise à une valeur trop élevée pour satisfaire au test portemanteau.

4.7.9.2.3.2 Avec ARIMA

## regr_ar = auto.arima(df$hadsst3gm, xreg = df$nrl2TsiAvg)
##
## Series: df$hadsst3gm 
## Regression with ARIMA(5,0,5) errors 
## 
## Coefficients:
##          ar1     ar2      ar3      ar4     ar5     ma1      ma2     ma3     ma4
##       0.0653  0.7917  -0.8496  -0.1234  0.8227  0.5105  -0.2809  0.7922  0.6033
## s.e.  0.3650  0.1537   0.0650   0.0743  0.3516  0.3736   0.3607  0.3726  0.1146
##           ma5  intercept    xreg
##       -0.2623     0.0011  1.0080
## s.e.   0.2686     0.0147  0.0562
## 
## sigma^2 estimated as 0.005875:  log likelihood=2343.56
## AIC=-4661.13   AICc=-4660.95   BIC=-4588.1
## 
## Training set error measures:
##                        ME      RMSE        MAE MPE MAPE      MASE         ACF1
## Training set 5.937753e-05 0.0764199 0.05646927 Inf  Inf 0.9215865 -0.001125632

Figure 4.40: Les résidus de l’analyse TS3 avec ARIMA.

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(5,0,5) errors
## Q* = 24.086, df = 18, p-value = 0.1522
## 
## Model df: 12.   Total lags used: 30

Il n’y a pas de tendance dans les résidus qui sont quasiment symétriques. Les coefficients d’auto-corrélation restent dans les limites admises et sont validés par le test portemanteau. Les intervalles de confiance des projections sont constants. La déviation standard de la projection est égale à celle des données historiques. Le coefficient sur l’ITS moyenne (nrl2TsiAvg) est quasiment égal à 1, ce qui valide l’analyse correspondante sans ARIMA.

Figure 4.41: Projections de l’analyse TS3 avec ARIMA.

4.7.9.2.4 Analyse résidus TS.4

4.7.9.2.4.1 Sans ARIMA

Figure 4.42: Les résidus de l’analyse TS4.

## 
##  Breusch-Godfrey test for serial correlation of order up to 10
## 
## data:  Residuals
## LM test = 1272.5, df = 10, p-value < 2.2e-16

Il n’y a pas de tendance dans les résidus qui sont symétriques. Le graphe d’auto-corrélation finit par rentrer dans les limites admises, mais le test portemanteau n’est pas satisfait.

4.7.9.2.4.2 Avec ARIMA

## regr_ar = auto.arima(df$hadsst3gm, xreg = cbind(df$harm,df$nrl2TsiAvg))
##
## Series: df$hadsst3gm 
## Regression with ARIMA(5,0,5) errors 
## 
## Coefficients:
##          ar1     ar2      ar3      ar4     ar5     ma1      ma2     ma3     ma4
##       0.0671  0.8154  -0.8652  -0.1736  0.8043  0.4997  -0.3293  0.7787  0.6393
## s.e.  0.0203  0.0334   0.0350   0.0329  0.0285  0.0296   0.0574  0.0264  0.0397
##           ma5   xreg1   xreg2
##       -0.2422  1.0041  1.0019
## s.e.   0.0326  0.1620  0.0451
## 
## sigma^2 estimated as 0.005821:  log likelihood=2352.91
## AIC=-4679.82   AICc=-4679.64   BIC=-4606.79
## 
## Training set error measures:
##                        ME       RMSE        MAE MPE MAPE      MASE         ACF1
## Training set 7.192395e-05 0.07607203 0.05621399 Inf  Inf 0.9174202 -0.00112353925632

Figure 4.43: Les résidus de l’analyse TS4 avec ARIMA.

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(5,0,5) errors
## Q* = 24.389, df = 18, p-value = 0.1427
## 
## Model df: 12.   Total lags used: 30

Il n’y a pas de tendance dans les résidus qui sont quasiment symétriques. Les coefficients d’auto-corrélation restent dans les limites admises et sont validés par le test portemanteau. Les intervalles de confiance des projections sont constants. La déviation standard de la projection est égale à celle des données historiques. Les coefficients sur le terme harmonique (harm) et sur l’ITS moyenne (nrl2TsiAvg) sont quasiment égaux à 1, ce qui valide l’analyse correspondante sans ARIMA.

Figure 4.44: Projections de l’analyse TS4 avec ARIMA.

4.7.10 Projections d’activité solaire et de température

Pour tester la capacité prédictive du modèle, il faut disposer de projections de l’intensité solaire, ou de son proxy, l’irradiance. Différents modèles théoriques d’irradiance ont été développés.

4.7.10.1 Modèles projectifs d’activité solaire

Le modèle de (Salvador 2013) est fondé sur les conjonctions des planètes qui influencent les orbites terrestres et l’activité solaire, qui dépendrait de la position du soleil par rapport au centre de gravité du système solaire. Les constantes du modèle sont disponibles à cette adresse. Sa reconstruction des taches solaires depuis l’an 1000 reconstitue fidèlement les minima historiques (Voir Figure 1.2). Salvador considère que son modèle est précis pour les 2 ou 3 prochains cycles solaires et devrait donner une bonne projection pour les 100 prochaines années qui devraient voir un grand minimum solaire, et donc aussi une mini-période glaciaire.

Voir Figure 4.45.

Figure 4.45: Un grand minimum solaire à l’horizon. Voir (Salvador 2013).

Le modèle de (Zharkova et al. 2015) est basé sur une analyse détaillée du champ magnétique solaire des cycles solaires 21 à 24. Le champ magnétique solaire découlerait de la somme de 2 composantes ayant des fréquences proches produisant une modulation dont la période est d’environ 350 ans qui se superpose à un cycle d’environ 22 ans (deux fois la durée du cycle de Schwabe). Cette période de 22 ans correspond à une inversion des pôles magnétiques du soleil. L’allure des variations de ce champ magnétique présente une ressemblance remarquable avec les taches solaires et reproduit également leurs minima historiques. Voir Figure 1.3. Zharkova considère que son approche ouvre une nouvelle ère dans la projection de l’activité solaire à l’échelle du millénaire. Elle pense également que nous nous dirigeons vers un grand minimum solaire.

Dans (Velasco Herrera, Mendoza, et He 2015) il y a deux reconstructions de l’irradiance totale solaire (ITS) depuis l’an 1000 avec des projections jusqu’en 2100. Ces reconstructions ont été calibrées sur base des composites PMOD et ACRIM qui utilisent des algorithmes différents pour calculer l’ITS. L’ITS PMOD est la plus récente. Le satellite à la base de l’ITS ACRIM n’est plus opérationnel. Il n’y a pas eu de période de recouvrement permettant de calibrer les mesures du dernier satellite avec le précédent. Les spécialistes ne sont pas d’accord entre eux quant à la manière de raccorder les séries de mesures. Les périodes et durées des grand minima historiques se retrouvent dans chaque reconstruction de Velasco, avec peu de différences entre elles.

4.7.10.2 Projections de température

Dans le présent travail, des calculs de régression analogues à ceux des analyses TSx ont été effectués pour les modèles de Salvador et Velasco. Les données publiées du modèle de Zharkova contiennent trop peu d’informations pour le faire. Un terme harmonique a été utilisé dans chacune des analyses.

Figure 4.46: En bleu, trait fin : l’anomalie de température brute (hadsst3gm) et en trait épais, son lissage par spline. En rouge, la modélisation de Velasco (PMOD) En mauve, la modélisation de Velasco (ACRIM) En orange, la modélisation de Salvador. La modélisation PMOD de Velasco reproduit le mieux les minima de Maunder et de Dalton.

Les modélisations présentent toutes un coefficient de corrélation fort élevé, ce qui ne les empêche pas d’être sensiblement différentes lorsqu’on sort de l’intervalle d’observation. Ce résultat confirme qu’il est beaucoup plus difficile d’extrapoler des données que de les interpoler. Les projections que l’on peut faire à un siècle doivent donc être considérées avec toutes les réserves d’usage, et cela quelle que soit la méthode utilisée pour les calculer.

Moyennant ces réserves, on constate que les trois simulations prédisent, à des degrés divers, un refroidissement pour les décennies à venir. La diminution de température jusqu’en 2070 devrait être (au moins) aussi rapide que l’augmentation de température encourue à la fin du siècle dernier. A cet horizon, la température serait redevenue semblable à celle du début du siècle passé, voire beaucoup plus basse, et cela de manière purement naturelle, sans aucune intervention humaine.

Ces résultats sont conformes aux avis d’une longue liste d’auteurs (Voir (Mörner 2015)) pour lesquels un refroidissement significatif, voire un nouveau petit âge glaciaire doit être envisagé, plutôt qu’un réchauffement catastrophique, tel qu’annoncé par le GIEC.