Modèle mathématique réchauffement climatique

Bernard Parisse Institut Fourier UMR 5582 du CNRS Universit� de Grenoble I

R�sum�: Ce texte reprend certains th�mes des sciences du climat o� interviennent des math�matiques au niveau Licence. On abordera la d�finition des climats, les causes naturelles de variation du climat (influence des param�tres orbitaux de la Terre), et les causes anthropiques (effet de serre).

Table des mati�res

  • 1��Introduction
  • 2��Le rayonnement solaire.
    • 2.1��L’insolation au cours de l’ann�e.
    • 2.2��Les saisons
    • 2.3��L’orbite de la Terre.
      • 2.3.1��Calcul en utilisant le vecteur excentricit�.
      • 2.3.2��Calcul par l’�quation diff�rentielle.
      • 2.3.3��Lois de K�pler.
    • 2.4��Quelques propri�t�s de l’ellipse
    • 2.5��Influence de l’ellipse sur les saisons
    • 2.6��L’�quation du temps, la dur�e des saisons.
    • 2.7��Les variations des param�tres orbitaux
  • 3��Les climats du pass�
    • 3.1��La transform�e de Fourier discr�te.
    • 3.2��La transform�e de Fourier rapide
    • 3.3��Application � la temp�rature
  • 4��L’effet de serre
    • 4.1��La d�finition du climat.
      • 4.1.1��Lois
      • 4.1.2��Ind�pendance
      • 4.1.3��Moyenne et �cart-type
      • 4.1.4��La loi normale
      • 4.1.5��Th�or�me central-limite
    • 4.2��Corr�lation et r�gression
      • 4.2.1��Corr�lation entre deux variables al�atoires.
      • 4.2.2��R�gression lin�aire avec un mod�le lin�aire
      • 4.2.3��Autres r�gressions lin�aires
      • 4.2.4��Corr�lations entre gaz � effet de serre et temp�rature
  • 5��Mod�lisation num�rique (climat, �nergie)
    • 5.1��Exemple�: un mod�le tr�s simplifi� du climat
    • 5.2��Exemple�: sc�narios d’�mission de CO2
      • 5.2.1��CO2 et consommation d’�nergie
      • 5.2.2��Sc�narios de consommation d’�nergie et d’�missions de CO2
    • 5.3��Mod�les discrets.
    • 5.4��Mod�les continus�: �quations et syst�mes diff�retiels
      • 5.4.1��Le mod�le logistique
      • 5.4.2��Syst�mes diff�rentiels lin�aires
      • 5.4.3��Les mod�les non lin�aires
      • 5.4.4��Exemples de mod�les de ce type
    • 5.5��Les mod�les de type EDP
  • 6��Energies renouvelables.
    • 6.1��Le solaire thermique
    • 6.2��Le solaire � concentration.
    • 6.3��Les aspects �conomiques.
  • 7��Conclusions
  • 8��R�f�rences sur le web
  • A��Pr�cession.

1��Introduction

La temp�rature du globe d�pend d’une part de l’intensit� et de la r�partition du rayonnement solaire, d’autre part de la chaleur r��mise vers l’espace.

Lorsque le Soleil est � la verticale on re�oit en effet beaucoup plus d’�nergie que s’il est bas sur l’horizon. La position du Soleil et la dur�e de l’ensolleillement journalier d�pend d’une part de l’endroit de la Terre o� l’on se situe (la latitude: pr�s des poles, des zones temp�r�es, des tropiques ou de l’�quateur), et d’autre part de l’inclinaison de l’axe de rotation de la Terre et de la position de la Terre sur son orbite, responsables des saisons, cf. figure 1.

Modèle mathématique réchauffement climatique

Figure 1: Saisons et obliquit� (source http://www.mnhn.fr/mnhn/lop/DIOGENE/morphologie/orbite.html)

Le calcul de l’inclinaison des rayons du Soleil en fonction de la date et de la latitude fait intervenir principalement de la trigonom�trie et un peu de g�om�trie dans l’espace.

La puissance du rayonnement solaire est inversement proportionnelle � la distance au carr� entre la Terre et le Soleil (en admettant que la puissance du Soleil soit constante � la source), la forme de l’orbite de la Terre influe donc sur le climat. Si on n�glige les autres plan�tes et on assimile que la Terre est un point mat�riel, l’orbite de la Terre est une ellipse (c’est-�-dire un cercle aplati dans une direction) dont le Soleil est un foyer, cf. figure 2.4 Au cours de son orbite la Terre se rapproche ou s’�loigne du Soleil ce qui fait l�g�rement varier le contraste et la dur�e des saisons (actuellement les saisons sont moins contrast�es dans lh�misph�re Nord que Sud et l’hiver bor�al est plus court que l’�t� bor�al ce qui explique que f�vrier n’ait que 28 jours, plus pr�cis�ment l’hiver bor�al dure 88.99 jours, le printemps 92.75 jours, l’�t� 93.65 jours, l’automne 89.85 jours). Le calcul pr�cis de la puissance de l’ensolleillement n�cessite donc la connaissance d’un peu de g�om�trie de l’ellipse, ainsi que de la r�solution num�rique d’�quations.

En r�alit� l’ellipse d�crite par la Terre autour du Soleil se d�forme au cours du temps, son excentricit� varie ainsi que les dates de passage au plus pr�s et au plus loin du Soleil, cf. figure 1 et voir pourquoi en annexe

Modèle mathématique réchauffement climatique

Figure 2: Pr�cession des �quinoxes

De plus l’angle d’inclinaison de l’axe de rotation de la Terre varie lentement ce qui r�partit l’�nergie solaire diff�remment, plus l’inclinaison est grande plus les poles en re�oivent, cf. figure 1.

Modèle mathématique réchauffement climatique

Figure 3: Variation de l’obliquit� (source Cyril Langlois)

Ces variations sont � la base de la th�orie astronomique des climats de Milankovitch. L’�tude num�rique des enregistrements du climat du pass� sur plusieurs centaines de milliers d’ann�es (en particulier les forages en Antarctique, par exemple Vostok, Dome C, projet Epica) est en bon accord avec cette th�orie, on peut ainsi mettre en �vidence (par transform�e de Fourier) les p�riodicit�s principales de ces enregistrements et les comparer avec celles calcul�es par les astronomes pour l’excentricit�, la pr�cession des �quinoxes et l’inclinaison de l’axe de rotation de la Terre, cf. figure 1.

Modèle mathématique réchauffement climatique

Figure 4: Temp�rature et gaz � effet de serre (Vostok)

L’arriv�e au sol et la r��mission de chaleur par la Terre est elle influenc�e par les nuages et la pr�sence dans l’air de mol�cules ayant des fr�quences de vibration propres proches de celles du rayonnement �mis par le sol (qui peuvent rendre l’air opaque � ce type de rayonnement) provoquant l’effet de serre. Les mol�cules principalement responsables de l’effet de serre sont le dioxyde de carbone, le m�thane, les oxydes d’azote, l’ozone ainsi que les gaz de type halofluorocarbone. Ainsi un doublement de concentration de dioxyde de carbone �quivaut � une augmentation d’environ 4 watt par m�tre carr� au niveau du sol (� comparer aux 1360/4 watt par m�tre carr� re�u au sommet de l’atmosph�re). Leur concentration a vari� dans le pass� de mani�re naturelle en phase avec les enregistrement de temp�ratures relev�s dans les glaces, cf. ci-dessus, on peut le mettre en �vidence de mani�re statistique, ici par une r�gression lin�aire de la temp�rature en fonction du CO2, figure 1

Modèle mathématique réchauffement climatique

Figure 5: Temp�rature en Antarctique et taux de CO2

Depuis le d�but de l’�re industrielle, cette concentration a beaucoup augment�, largement au-del� des limites naturelles atteintes au cours des 750 mille derni�res ann�es (par exemple 380 parties par million pour le CO2 a Mauna Loa, largement en-dehors de l’intervalle 180-290, en hausse de 2 ppm par an, figure 1)

Modèle mathématique réchauffement climatique

Modèle mathématique réchauffement climatique

Figure 6: �volution r�cente du taux de CO2 (absolue et moyenne mobile de la hausse)

faisant craindre un r�chauffement de la plan�te au 21�me si�cle et au cours des si�cles suivants, r�chauffement trs probablement d�j� commenc�, figure 1.

Modèle mathématique réchauffement climatique

Figure 7: Temp�rature actuelle (source IPCC)

La cause principale de l’augmentation de ces concentrations est l’utilisation d’�nergie fossiles pour le CO2 (p�trole, gaz naturel et charbon) et d’engrais azot�s (agriculture intensive). Une simple r�gle de trois permet de faire la correspondance

Modèle mathématique réchauffement climatique
Modèle mathématique réchauffement climatique
Modèle mathématique réchauffement climatique

Figure 8: Carbone �mis par tonne �quivalent p�trole (source J.M. Jancovici) et exemple de sc�nario �nerg�tique et de population pour le 21i�me si�cle (IPCC et ONU)

Ainsi, en un an (d�cembre 2004-novembre 2005), la France a consomm� 12.9tonnes �quivalent p�trole (tep) de charbon, 95.1 tep de p�trole et 40.5 tep de gaz, dont la combustion a �mis dans l’atmosph�re 122 millions de tonnes de carbone (sous forme de CO2), soit 1.8 tonne de carbone par habitant, si le reste du monde consommait au m�me rythme cela correspondrait � 6 parties par millions de CO2 �mises par an.

La quantification de l’ampleur du r�chauffement (et des autres param�tres climatiques en particulier la r�partition et l’intensit� des pr�cipitations et des vents) se fait en deux temps, d’une part par des sc�nari d’emission de gaz � effet de serre, d’autre part par mod�lisation num�rique du climat en fonction des sc�nari d’�mission. La mod�lisation num�rique donne des r�sultats relativement pr�cis pour la temp�rature (3 � 4 degr� pour un doublement du CO2), plus incertains pour les pr�cipitations.

Par contre la mod�lisation des �missions de gaz � effet de serre est beaucoup plus difficile, d’une part parce que les r�serves de combustibles fossiles utilisables sont controvers�es (par exemple sur la date du pic de production p�trolier), d’autre part parce que le mix �nerg�tique du 21�me si�cle est largement inconnu (fossile, nucl�aire, renouvelables?), mais aussi parce que la r�ponse de la biosph�re au r�chauffement est encore mal quantifi�e (par exemple on estime que l’oc�an absorbe actuellement environ la moiti� des �missions de CO2, mais la solubilit� du gaz carbonique dans l’eau diminue lorsqu’on la r�chauffe). Les math�matiques interviennent ici d’abord pour d�finir le climat (par exemple pourquoi on utilise une s�rie de 30 ann�es pour d�finir une temp�rature moyenne en un lieu donn�), dans la r�alisation des mod�les num�riques (nous ne parlerons pas de ce sujet qui d�passe le niveau licence scientifique), dans les mod�les statistiques de r�chauffement dus aux gaz � effet de serre (corr�lation), ainsi que dans les mod�les d’�mission de gaz � effet de serre, par l’interm�diaire des mod�les de production et consommation de combustibles fossiles.

Modèle mathématique réchauffement climatique
Modèle mathématique réchauffement climatique
Modèle mathématique réchauffement climatique

Figure 9: Production de p�trole aux �tats-Unis (source EIA et Oildrum) mettant en �vidence le ph�nom�ne de pic p�trolier et recherche du meilleur mod�le de pr�diction 1976-2004 connaissant l’avant 1976. En millions de barils par jour (1 baril=159 litres). Mod�lisation des ressources ultimes en p�trole en milliards de barils et de production d’�nergie en combustibles fossiles en milliards de tonne �quivalent p�trole par an (source Lah�rr�re)

2��Le rayonnement solaire.

2.1��L’insolation au cours de l’ann�e.

Pour connaitre la quantit� d’�nergie recue � un moment donn�, il faut calculer l’angle entre la verticale du lieu et la direction du Soleil. Plus g�n�ralement, on va calculer les composantes du vecteur Terre-Soleil et les composantes des vecteurs de la base locale (verticale locale, direction du Sud et direction du parall�le). On choisit d’abord comme r�f�rence le plan Txy de l’�cliptique (plan de l’orbite de la Terre autour du Soleil), avec Ty orthogonal � l’axe de rotation de la Terre (donc Tx la projection de l’axe de rotation de la Terre sur ce plan). Soit θ l’angle que fait la Terre avec la direction du passage au p�rih�lie, et θ0 l’angle de la position de la Terre au solstice d’hiver avec la direction du p�rih�lie, l’angle entre la direction Terre-Soleil et Tx est donc θ−θ0

Modèle mathématique réchauffement climatique

Dans Txyz le vecteur unitaire s de la direction Terre-Soleil a pour coordonn�es�:

s=−(cos(θ−θ0),sin(θ−θ0),0)

On effectue ensuite une rotation autour de Ty d’angle i l’inclinaison de l’axe de rotation de la Terre. On obtient ainsi un rep�re TXyZ (TX et TZ se d�duisent de Tx et Tz par rotation d’angle i). Dans ce rep�re le vecteur untaire s a pour coordonn�es�:

s=−(cos(θ−θ0)cos(i),sin(θ−θ0), cos(θ−θ0)sin(i))�

Calculons maintenant dans ce rep�re TXyZ les coordonn�es des vecteurs de la base locale. On se place en un point de latitude l et de longitude φ, on note J la dur�e d’une p�riode de r�volution de la Terre sur elle-m�me (23 heures 56 minutes, c’est un peu moins d’un jour car il faut encore en moyenne 4 minutes pour compenser le d�placement de la Terre sur son orbite autour du Soleil). La verticale locale a pour coordonn�es�:

v=(cos(l)cos(φ+2π�t/J),cos(l)sin(φ+2π�t/J),sin(l))

L’�nergie solaire recue au lieu donn� (sur une surface horizontale ; pour un panneau solaire, il faudrait calculer les coordonn�es d’un vecteur perpendiculaire au panneau) est proportionnelle �

o� ρ(θ)=a(1−e2)/(1+ecos(θ)) d�signe la distance Terre-Soleil. Le calcul de s.v donne, en notant ϕ=φ+2π t/J�:

−s.v��=��cos(l)�cos(ϕ)�cos(θ−θ0)cos(i) +�cos(l)sin(ϕ)�sin(θ−θ0)� +�sin(l)�cos(θ−θ0)sin(i)�

On rassemble les deux premiers termes qui d�pendent rapidement du temps par l’interm�diaire de ϕ (le 3�me terme n’en d�pend que par θ qui ne varie que d’environ 1 degr� pendant une journ�e) et on applique la formule de trigonom�trie:

A�cosα�+�B�sinα�=� cos(α−α0),����











Ici, apr�s avoir factoris� cos(l), on a�:

=
cos(θ−θ0)2cos(i)2+sin(θ−θ0)2�
=�

On peut aussi calculer

qui donne α0 modulo π et compl�ter en regardant le quadrant o� se trouve (A,B), ici α0 et θ−θ0 sont tous deux dans [0,π] ou tous deux dans [−π,0]. Finalement, on obtient le

Theorem�1�� L’�nergie solaire recue au sol est proportionnelle �

o� s.v est donn� par�:

s.v�=�−�cos(l)� cos(ϕ−ϕ0) −�sin(l)�cos(θ−θ0)sin(i)

et�

  • e est l’excentricit� de l’orbite elliptique (environ 0.0167 actuellement)
  • i est l’obliquit� (inclinaison de l’axe de rotation de la Terre, environ 23 degr� 27 minutes actuellement)
  • θ est l’angle fait par la direction Terre-Soleil avec la direction du demi grand axe (Soleil-p�rih�lie), θ0 le m�me angle au solstice d’hiver de l’h�misph�re Nord (environ -13 degr�s). En premi�re approximation, on peut faire varier θ proportionnellement au temps, voir la fin de la section 2.6 pour un calcul plus pr�cis.
  • l est la latitude, ϕ la longitude tenant compte de la rotation de la Terre (somme de la longitude g�ographique φ et du terme d�pendant du temps 2π t/J)
  • ϕ0 ∈ [−π,π] est de m�me signe que θ−θ0 et est d�fini par�:

Variations de s.v au cours d’une journ�e dans l’approximation o� θ ne varie pas�:
On obtient une sinusoide entre les deux valeurs extr�mes�:

��cos(l)� �− sin(l)�cos(θ−θ0)sin(i)�

Le maximum est atteint pour

ϕ=ϕ0+�π��⇒�� 2�π�t/J�=�ϕ0�−�φ�+�π,����J=23h56m�

le moment correspondant est appel� culmination (c’est le midi solaire si le maximum est positif) et ne d�pend pas de la latitude (bien entendu la valeur du maximum en d�pend). Si le maximum est n�gatif ou nul, la nuit dure 24h. Si le minimum est positif ou nul, le jour dure 24h. Par exemple au solstice d’hiver, θ=θ0, selon la latitude on obtient un maximum n�gatif pour l=π/2 (pole Nord), positif pour l=−π/2 (pole Sud), le minimum et le maximum croissent entre ces 2 valeurs. Si le maximum est positif et le minimum est n�gatif, il y a 2 instants ou s.v=0 (lever et coucher du soleil).

L’�nergie solaire recue pendant une journ�e par une surface horizontale est proportionnelle � l’int�grale entre le lever et le coucher de s.v/ρ2. S’il n’y a pas de lever/coucher, soit on ne recoit rien (nuit polaire), soit on recoit l’int�grale entre 0 et 24h de s.v (jour polaire).

L’intervalle entre 2 culminations n’est pas constant au cours de l’ann�e, car ϕ0 n’est pas une fonction lin�aire de θ (qui lui m�me n’est pas lin�aire en fonction du temps sauf en premi�re approximation avec une orbite terrestre circulaire). On peut le calculer en d�rivant (1). Par exemple dans l’approximation d’une excentricit� nulle, au solstice d’hiver (θ=θ0), on obtient

ϕ0=�0,����(1+02)�dϕ0�=�(1+02)

avec dθ qui correspond � 4 minutes, on trouve dϕ0 correspondant � 4.36 minutes. L’�cart entre 2 culminations est donc d’environ 24h 20secondes. Au moment du solstice, le Soleil se l�ve et se couche donc environ 20 secondes plus tard entre un jour et son lendemain, dans l’hypoth�se d’un mouvement circulaire de la Terre autour du Soleil. En r�alit�, l’orbite terrestre �tant faiblement elliptique, l’�cart est un peu moins de 30 secondes en hiver et de 15 secondes en �t�, le mouvement de la Terre autour du Soleil �tant plus rapide d’environ 3% au solstice d’hiver et moins rapide d’environ 3% au solstice d’�t�. Comme 3% de l’�cart moyen entre 2 culminations (4 minutes=240 secondes) correspond � 7 secondes cela explique la diff�rence.

2.2��Les saisons

Les deux figures (2.2, 2.2) qui suivent repr�sentent la Terre en 4 points de son orbite, les 2 solstices et les 2 �quinoxes. Les solstices sont d�finis par les 2 points de l’orbite o� la projection de l’axe de rotation terrestre est parall�le � l’axe Terre-Soleil. Les �quinoxes sont d�finis par les 2 points o� il y a perpendicularit�.

Modèle mathématique réchauffement climatique

Figure 10: Saisons (h�misph�re Nord), source wikipedia

Modèle mathématique réchauffement climatique

Figure 11: Saisons (h�misph�re Sud), source wikipedia

Au solstice d’hiver, on voit que les parall�les situ�s aux hautes latitudes Nord ne sortent jamais de l’obscurit�. Aux latitudes int�rm�diaires, le morceau de parall�le situ� au jour est nettement plus petit que celui situ� dans l’obscurit�. � l’�quinoxe de printemps, chaque parall�le est � moiti� au jour et � moiti� dans l’obscurit� (derri�re la grille). Au printemps, la situtation est analogue. Au solstice d’�t�, on est dans la situation inverse de l’hiver.

2.3��L’orbite de la Terre.

En premi�re approximation, l’orbite de la Terre est uniquement influenc�e par la force de gravitation entre la Terre et le Soleil, ce dernier pouvant �tre consid�r� comme fixe en raison de sa masse (on peut �viter cette approximation en rempla�ant le Soleil par le centre de gravit� du syst�me Terre-Soleil). La force de gravitation qui d�rive d’un potentiel inversement proportionnel � la distance Terre-Soleil est de la forme

o� r d�signe le vecteur Terre-Soleil, mT est la masse de la Terre, �=GmS est le produit de la constante de gravitation universelle par la masse du Soleil. Le moment cin�tique de la rotation de la Terre autour du Soleil est d�fini par�:

On v�rifie que sa d�riv�e est nulle, donc L garde une direction fixe k, orthogonale � r, l’orbite de la Terre reste donc dans le plan d�fini � un instant donn� par l’axe Terre-Soleil et le vecteur vitesse de la Terre. De plus la conservation de L entraine la loi des aires, l’aire balay�e par le rayon Soleil-Terre est proportionnelle au temps.

On utilise un rep�re en coordonn�es polaires centr� au Soleil, ρ d�signant la distance Terre-Soleil et θ l’angle fait par rapport � une direction fixe, on a alors

car si on calcule en coordonn�es polaires d r/dt, la composante sur le vecteur radial er est dρ/dt, et la composante sur le vecteur perpendiculaire eθ est ρ dθ/dt.

2.3.1��Calcul en utilisant le vecteur excentricit�.

Montrons que le vecteur

est aussi conserv� (o� on rappelle que � provient de la force de gravitation F′=F/mT=−� r/r3). Le deuxi�me terme est proportionnel au vecteur radial −er, dont la d�riv�e est le vecteur orthogonal −dθ/dt eθ. Comme L est constant, la d�riv�e du premier terme est

F′�∧�L�=� �∧�L�k =� eθ�=�−� �ρ2

Notons que E est dans le plan de l’orbite, prenons comme origine des angles pour rep�rer la Terre par rapport au Soleil la direction de E. En faisant le produit scalaire de E avec r, on obtient en notant e la norme de E

e�ρ�cos(θ) = E.r
  =
  =
  =
  =
  =

d’o��:

2.3.2��Calcul par l’�quation diff�rentielle.

On a les �quations de conservation de l’�nergie et du moment cin�tique�:

+











+�


ρ�








=C1,����ρ2� �=�L,����K<0,�m>0

On change de variable d�pendante pour ρ, en prenant θ au lieu de t, comme dρ/dt=ρ′ dθ/dt (o� ρ′=dρ/dθ), on a�:

+





(ρ′2+ρ2)�











=�C1

On effectue ensuite le changement de variable ρ=1/u, ρ′=−u′/ u2, d’o��:

K�u�+�





+


L2�u4�


=�C1

soit�:

donc�:

K′�u�+�u2�+�u′2�=�C3,���� K′= �<0,�C3=

On pose maintenant v=u+K′/2, d’o��:

On montre (en exprimant v′ en fonction de v puis en s�parant les variables), que cette �quation diff�rentielle a pour solution g�n�rale

D’o��:

Comme K′<0 et comme la trajectoire de la Terre autour du Soleil passe par tous les angles (donc ρ est d�fini pour tout θ, le d�nominateur ne peut pas s’annuler), on a�:

On d�finit ensuite a par 2/−K′=a(1−e2), et on obtient finalement l’�quation d’une ellipse dont l’origine (le Soleil) est un des foyers�:

On suppose maintenant quitte � faire pivoter l’axe des x que θ0=0.

2.3.3��Lois de K�pler.

L’orbite de la Terre est donc une ellipse dont le Soleil occupe un des foyers (1�re loi de K�pler). On a aussi vu que L=ρ2 dθ/dt est constant, ceci entraine la loi des aires, infinit�simalement on a�:

ce qui se traduit par l’aire balay�e par le rayon vecteur Soleil-Terre est proportionnelle au temps (2�me loi de K�pler). Au cours d’une p�riode T, l’aire parcourue est celle de l’ellipse, donc

En prenant le carr�, et en appliquant

on en d�duit la troisi�me loi de K�pler�:

o� on rappelle que � est le produit de la constante de gravitation universelle par la masse du Soleil. (On peut �videmment faire le m�me calcul pour la Lune autour de la Terre).

2.4��Quelques propri�t�s de l’ellipse

D�finition
L’ellipse E de foyers F1 et F2 de demi-grand axe a est l’ensemble des points M du plan tels que

Cf. figure 2.4.

Modèle mathématique réchauffement climatique

Figure 12: ellipse (source Florence Kalfoun)

On note 2c=F1F2 la distance entre les deux foyers, qui doit �tre plus petite que 2a pour que l’ellipse soit non vide. L’excentricit� de l’ellipse est d�finie par e=c/a < 1. Si e=0, on obtient un cercle de centre F1=F2 et de rayon a. Si e≠ 0, on va voir qu’il s’agit d’un cercle contract� selon l’axe perpendiculaire � F1F2 dans un rapport de √1−e2. On va �galement calculer l’�quation en coordonn�es polaires de E pour montrer que l’�quation obtenue ci-dessus est bien celle d’une ellipse dont le Soleil occupe un foyer.

Soit O le milieu de F1 et F2, on se place dans le rep�re orthonorm� dont le premier axe Ox contient F1 et F2 donc les coordonn�es de F1 sont (c,0) et celles de F2 sont (−c,0). Soit M(x,y) un point de l’ellipse, on a d’une part�:

MF12�−�MF22�=�(x−c)2−(x+c)2�=�−4cx�

et d’autre part�:

MF12�−�MF22�=�(MF1�+�MF2)(MF1�−�MF2�)�=�2a�(MF1�−�MF2�)

donc�:

en additionnant avec MF1+MF2=2a et en appliquant c=ea, on en d�duit�:

� MF1�=�a�−� �=�a−ex� ����(2)

En prenant le carr�, on a�:

d’o��:

y2�+�x2�(1−e2)�=�a2(1−e2)�

finalement�:

qui est bien la contraction selon Oy de rapport √1−e2 du cercle de centre O et de rayon a (appel� grand cercle de l’ellipse).

En coordonn�es polaires, on note ρ la distance de F1 � M, et θ l’angle entre l’axe Ox et F1M. L’abscisse de M est donc�:

que l’on combine avec (2) pour obtenir�:

ρ�=�a−ex�=a(1−e2)�−�e�ρ�cos(θ)�

donc�:

ce qui nous permet d’affirmer que l’orbite de la Terre dans l’approximation du point mat�riel soumis uniquement au Soleil suppos� fixe est une ellipse dont le Soleil occupe un foyer.

2.5��Influence de l’ellipse sur les saisons

Il faut prendre garde � ne pas confondre les solstices et �quinoxes avec le moment o� la Terre coupe le grand axe de son ellipse autour du Soleil. Il n’y a aucune raison que la projection de l’axe de rotation de la Terre sur le plan de l’ellipse soit parall�le ou perpendiculaire au grand axe de l’ellipse, et actuellement ce n’est pas le cas, le solstice d’hiver a lieu le 21 d�cembre alors que le passage au plus proche du Soleil a lieu vers le 3 janvier (donc pendant l’hiver de l’h�misph�re Nord) et le passage au plus loin du Soleil a lieu d�but juillet (pendant l’�t�). C’est pour cette raison que les saisons sont moins marqu�es dans l’h�misph�re Nord que dans l’h�misph�re Sud. De plus la loi des aires oblige la Terre a se d�placer plus vite lorsqu’elle est proche du Soleil que lorsqu’elle en est �loign�e ce qui diminue la dur�e de l’hiver bor�al et augmente la dur�e de l’�t� bor�al (c’est pour cette raison que f�vrier n’a que 28 jours alors que juillet et aout ont 31 jours).

2.6��L’�quation du temps, la dur�e des saisons.

Modèle mathématique réchauffement climatique

Figure 13: Ellipse et �quation du temps

La trajectoire elliptique E de la Terre autour du Soleil est repr�sent�e sur la figure 2.6 en bleu, l’excentricit� de l’orbite a �t� �norm�ment exag�r�e, il s’agit d’une ellipse de foyers S (le Soleil) et S′. Le point A d�signe le p�rih�lie de l’orbite (passage de la Terre au plus proche du Soleil), qui a lieu vers le 4 janvier. En noir, on a dessin� le grand cercle de l’ellipse (l’ellipse s’obtient par contraction du grand cercle de rapport √1−e2 o� e est l’excentricit� de l’orbite). L’aire d�crite par le rayon Soleil-Terre (ST) est proportionnelle au temps (loi des aires qui d�coule de la conservation du moment cin�tique), il en est donc de m�me de l’aire (en vert) du d�crite par le rayon SM. Si on ajoute � cette aire verte l’aire en rouge du triangle OSM, on obtient l’aire de l’arc de cercle OAM. Donc

est proportionnel au temps �coul� depuis le passage au p�rih�lie. Comme HM=OM sin(V) et OS= e � OA, on en d�duit que

� V�−�e�sin(V)�=�C�t�=�2π� ����(3)

o� la constante C s’obtient en faisant varier V de 0 � 2π ce qui correspond � la dur�e T d’une r�volution de la Terre autour du Soleil (1 an).

La relation entre θ (not� t sur la figure) et V s’obtient par exemple en calculant l’abscisse de M

x = a�cos(V)�
  = ea�+�ρ�cos(θ)�
  =

Les angles V et θ sont de m�me signe et

et r�ciproquement�:

Dur�e des saisons�:
Il suffit de connaitre l’angle θ lors du solstice d’hiver et de lui ajouter kπ/2 pour k=1,2,3 pour connaitre l’angle θ au printemps, en �t� et � l’automne, on en d�duit V par (4) puis le temps �coul� depuis le p�rih�lie avec (3).

Calcul de θ en fonction du temps �coul� depuis le passage au p�rih�lie�:
Il faut calculer V par des m�thodes num�riques (point fixe ou m�thode de Newton) en appliquant (3), on en d�duit θ avec (5). En r�sum�, on a le�:

Theorem�2�� Soit θ l’angle entre le demi grand axe de l’ellipse et la direction Soleil-Terre, t∈ [−T/2,T/2] le temps �coul� depuis le passage au p�rih�lie (t=0 lorsque θ=0, T=1 an). Soit V∈ [−π,π] la solution de

o� e est l’excentricit� de l’ellipse. Alors θ est donn� par

2.7��Les variations des param�tres orbitaux

La Terre n’est pas une sph�re id�ale, elle a un renflement au niveau de l’�quateur, due � rotation de la Terre sur elle-m�me (la force centrifuge y est plus importante). Ce renflement est dans un plan qui fait un angle avec le plan de l’�cliptique, le Soleil exerce donc un couple sur ce renflement. Ce ph�nom�ne est � l’origine de la pr�cession des �quinoxes, le passage au p�rih�lie de la Terre se d�cale dans le temps. De plus, la Terre n’est pas seulement soumise � l’influence du Soleil, mais aussi des autres plan�tes, en particulier Jupiter. Cela modifie sur de tr�s longues p�riodes tous les param�tres de l’orbite terrestre, en particulier l’excentricit�, la pr�cession des �quinoxes, mais aussi l’obliquit� (inclinaison de l’axe de rotation terrestre par rapport � la perpendiculaire au plan de l’�cliptique). Le calcul de ces variations est bien au-del� des pr�tentions de ce texte, le lecteur int�ress� pourra se r�f�rer par exemple aux publications de Laskar (chercher ce mot-clef ou des mots comme orbite, perturbation, symplectique, hamiltonien, ...). On se bornera ici � indiquer que le demi-grand axe ne varie pas, ce qui donne une relation entre les variations de la constante des aires et de l’excentricit�

L� est proportionnel � �

Les variations des param�tre orbitaux modifient � long terme l’ensolleillement de la Terre (la valeur de l’�nergie re�ue en un lieu sur une surface horizontale s.v/ρ2 d�pend de la latitude, de la position de la Terre sur son orbite mais aussi de l’excentricit� de l’orbite, de l’obliquit� et de la date du p�rih�lie par rapport aux saisons) et sa r�partition sur le globe par latitude, il est naturel de supposer qu’elles influent sur le climat de la Terre. Par exemple, l’�nergie moyenne recue par la Terre au cours d’une p�riode T de une ann�e est donn�e par

est proportionnelle � 1/(T L) donc � (1−e2)−1/2 (car T est aussi constant d’apr�s la 3�me loi de K�pler). Au premier ordre, la variation de e entraine donc une variation de l’ensolleillement global de

Pour la Terre, cela repr�sente au plus 2.5 pour mille (la p�riode la plus favorable aux glaciations �tant celle o� l’orbite est circulaire), soit, sans r�troactions, une variation globale de 0.2 degr�s Kelvin.

3��Les climats du pass�

Les carottages de la glace accumul�e � Vostok (Antarctique) et ailleurs permettent de reconstituer la temp�rature en Antarctique (ou ailleurs) jusqu’� 700 mille ans. On peut effectuer une �tude num�rique des valeurs obtenues pour faire apparaitre les principales p�riodicit�s de ce signal num�rique, et comparer avec les p�riodicit�s de l’orbite terrestre. On utilisera aussi les donn�es de ces carottes pour mettre en �vidence les corr�lations entre gaz � effet de serre et temp�rature dans le pass�.

3.1��La transform�e de Fourier discr�te.

Soit N un entier fix�. Une suite x p�riodique de p�riode N est d�termin�e par le vecteur x=[x0,x1,...xN−1]. La transform�e de Fourier discr�te (DFT) not�e FN fait correspondre � une suite x p�riodique de p�riode N une autre suite y p�riodique de p�riode N, d�finie pour k=0..N−1 par :

o� ωN est une racine N-i�me primitive de l’unit�, on prend ω=e2iπ/N si x est � coefficients r�els ou complexes. Cette transformation est lin�aire, la transform�e de la somme de 2 suites est la somme des transform�es, et la transform�e du produit par une constante d’une suite est le produit par cette constante de la transform�e de la suite.

La transform�e de Fourier discr�te FN est une transformation bijective dont la r�ciproque est donn�e par�:

FN−1= ,��� (FN−1(y))k= yjωNk��j

On le prouve en rempla�ant y par sa valeur�:

(FN−1(y))k =
  =
  =
  =
�xl�





si ωN(k−l)≠�1�
N si ωN(k−l)�=�1�

Or si ωNk−l=e2iπ(k−l)/N=1 si et seulement si k=l d’o� le r�sultat.

Propri�t�
La transform�e de Fourier discr�te d’une suite r�elle v�rifie yN−k=yk.
La preuve est imm�diate en appliquant la d�finition.

Un des int�r�ts de la DFT est de mettre en �vidence rapidement d’�ventuelles p�riodicit�s de x divisant N. Plus pr�cis�ment soit j est un entier divisant N. Consid�rons une suite r�elle x dont la DFT y est nulle sauf yl et yN−l. Par lin�arit�, on peut se ramener � 2 cas yl=yN−l=1 et yl=i, yN−l=−i. Dans le premier cas, on obtient xk=ωNlk+ωN−lk=2cos(2π kl/N), dans le deuxi�me cas, on obtient xk=−2sin(2π kl/N), qui sont p�riodiques de p�riode N/l.

R�ciproquement, si x a comme p�riode T=N/l, alors en posant j=T m + r avec m∈[0,l[ et r∈[0,T−1], on a xj=xr donc�:

si ωN−kT ≠ 1. Comme (ωN−kT)l=ωN−klT=1, yk=0 si kT=kN/l n’est pas un multiple de N. Finalement si k n’est pas un multiple de l, alors yk=0.

Voyons maintenant le cas de “pseudo-p�riodes”, supposons donc que x est p�riodique de p�riode N mais que de plus pour un T>0 quelconque (ne divisant pas forc�ment N), on ait

xj+T=xj,����∀�j�∈[0,N−T]�

On peut refaire le raisonnement ci-dessus, modulo des erreurs. plus pr�cis�ment�:

yk�−� xj�ωN−k��j� =� �xr�ωN−k�(T�m+r)�

On calcule donc yk � une erreur de ceil(N/T)T−N termes major�s par |xj| pr�s. Et le membre de droite vaudra�:

�xr�ωN−kr�
1−ωN−kceil(N/T)�T
1−ωN−kT

Le module de la fraction est �gal �

|
sin(π�k�ceil(N/T)�T/N)
sin(π�k�T/N)
�| =�|
sin(π�k�(ceil(N/T)�T/N−1))
sin(π�k�T/N)
�|

il est petit si k n’est pas proche d’un multiple de ceil(N/T). Par exemple, prenons N=216=65536 et T ≈ N/10 =6554. Dans ce cas ceil(N/T)T=10 � 6554=65540, il y a donc une erreur de 4 termes sur le calcul de yk. Si k n’est pas proche d’un multiple de 10, on doit trouver yk proche de 0 relativement � la valeur des |xj|.

Les p�riodes et pseudo-p�riodes de x correspondent donc aux valeurs de yk grandes par la r�gle k * p�riode =N.

3.2��La transform�e de Fourier rapide

Le calcul de la DFT est relativement lent, il n�cessite de l’ordre de N2 op�rations, car il revient � calculer la valeur du polyn�me de degr� N−1:

aux N points 1, ωN, ..., ωNN−1 (on a yk=P(ωNk)). Mais si N est une puissance de 2, on peut calculer de mani�re plus astucieuse et r�duire le nombre d’op�rations � un ordre N ln(N). En effet N=2M, on d�coupe P en 2 parties de m�me longueur�:

on a alors

P(ω2k)=�(Q+R)�((ω2)k),��� P(ω2k+1)=�(−Q�+R)ω((ω2)k)��� Sω(x)=sk�wk�xk

On est donc ramen� � deux additions de 2 polyn�mes de degr� M, une multiplication coefficient par puissances (4M op�rations), et au calcul des deux DFT de Q+R et R−Q. Si N=2n on v�rifie que cela n�cessite O(n2n) op�rations, donc Nln(N) op�rations. On appelle alors FFT cette m�thode de calcul (DFT=FFT si N2n).

3.3��Application � la temp�rature

Les enregistrements ne sont pas faits � des dates r�guli�rement espac�es. Pour se ramener � une suite, il faut interpoler la valeur de la temp�rature � une date donn�e en fonction de la temp�rature aux dates les plus proches de l’�chantillon. Le calcul de la FFT de l’�cart de temp�rature avec la temp�rature actuelle sur les 400 000 derni�res ann�es donne le graphe figure 3.3

Modèle mathématique réchauffement climatique

Figure 14: Norme de la FFT de la temp�rature de Vostok

On observe que la FFT yk est de module relativement grand pour k=4, k=10 et k proche de 20. La valeur de k=4 correspond � une p�riodicit� de 100 000 ans, celle pour k=10 � 40 000 ans. Il est plus difficile de tirer des conclusions du maximum relatif pour k proche de 20 car 20 est un multiple des p�riodes pr�c�dentes 4 et de 10.

On peut comparer ces 2 p�riodes tir�es de l’enregistrement num�rique des p�riodes de variation des param�tres orbitaux de la Terre�:

  • l’excentricit� de la Terre varie avec des p�riodes principales de l’ordre de 100 000 et 400 000 ann�es
  • l’obliquit� de la Terre varie avec une p�riode de 41 000 ans.
  • la pr�cession des �quinoxes a des p�riodicit�s principales de 19 000 et 23 000 ans.

L’excentricit� et l’obliquit� ont donc des p�riodicit�s que l’on retrouve dans les enregistrements de la temp�rature du pass� ce qui confirme l’influence suppos�ee de ces param�tres sur le climat de la Terre.

La figure 3.3 montre la reconstruction par FFT inverse en tronquant la FFT aux 30 premiers entiers (et sym�triquement)

Modèle mathématique réchauffement climatique

Figure 15: Reconstruction de la temp�rature en tronquant la FFT aux 30 premiers entiers

4��L’effet de serre

Nous nous sommes int�ress�s jusqu’� pr�sent � la quantit� d’�nergie qui arrive du Soleil, sans regarder ce qui peut se passer entre le sommet de l’atmosph�re et le sol. Or la composition chimique de l’atmosph�re a une influence sur le pi�geage des radiations solaires, donc sur le climat. Les mol�cules des gaz � effet de serre ont des fr�quences propres de vibrations qui sont proche des fr�quences de radiation �mises par la Terre et les rend plus ou moins opaques (alors qu’ils sont quasiment transparents dans le domaine de fr�quence des radiations provenant du Soleil). Il s’agit sur Terre principalement du CO2, CH4, des chlorofluorocarbone, de l’ozone (O3) et d’oxydes d’azote. Leur concentration a vari� dans le pass�, de mani�re naturelle corr�l�e avec la temp�rature (cf. la section 4.2.2). Elle varie aujourd’hui du fait de l’homme, et les scientifiques �laborent des mod�les pour pr�voir la r�action du climat � ces variations.

4.1��La d�finition du climat.

Le domaine des math�matiques concern� sont les probabilit�s et les statistiques. Il s’agit de savoir quel sens attribuer � une phrase comme “la hausse de temp�rature de la Terre au 20�me si�cle a �t� de 0.6 degr�”. Il s’agit aussi de savoir pourquoi on consid�re que 30 ans est un nombre raisonnable d’ann�es pour d�finir le climat en un lieu donn�. On peut faire l’hypoth�se que le temps un jour donn� n’est pas ind�pendant de celui du jour pr�c�dent (m�moire du temps), mais qu’il est ind�pendant de celui du m�me jour l’ann�e pr�c�dente�: le temps suit une loi qui d�pend de la saison. Cette loi subit une lente d�rive au cours du temps (variations du climat). Un des r�sultats importants qui justifie la p�riode de 30 ann�es est le th�or�me central-limite. Nous allons rappeler quelques d�finitions avant de l’�noncer.

4.1.1��Lois

Lorsqu’on n’est pas capable de pr�dire de mani�re d�terministe plusieurs �v�nements du m�me type, on les mod�lise souvent en disant qu’ils d�pendent du hasard mais pas de mani�re anarchique, ils suivent une loi. Ainsi, le r�sultat du lancer d’un d� � 6 faces non pip� ne peut pas �tre pr�dit, mais en revanche il suit une loi (uniforme), chaque face a autant de chances d’apparaitre qu’une autre. La valeur d’un lancer est une variable al�atoire � valeurs dans {1,2,3,4,5,6} suivant la loi uniforme (la probabilit� d’observer chacune des faces de 1 � 6 est de 1/6). Il s’agit ici d’une loi discr�te, on ne peut observer que 6 valeurs. L’analogue existe pour une variable al�atoire dont la valeur est r�elle (par exemple la temp�rature), on utilise alors une densit� de probabilit�, c’est une fonction f de ℝ � valeurs dans ℝ+ telle que la probabilit� d’observer la variable al�atoire dans un intervalle [a,b] soit ∫ab f(x) dx. L’int�grale sur ℝ de f est donc �gale � 1. Par exemple pour la loi uniforme sur [1,6], f est constante et vaut 1/5 sur [1,6] et 0 ailleurs. Un autre exemple important est la loi normale d�velopp�e plus bas.

4.1.2��Ind�pendance

Lorsqu’on effectue plusieurs observations successives, une observation peut d�pendre de l’observation pr�c�dente ou non. Par exemple, la temp�rature observ�e un jour donn� est en g�n�ral proche de celle qui sera observ�e le lendemain. Mais ell est ind�pendante de celle observ�e l’ann�e pr�c�dente � la m�me date. De m�me qu’un jet de d� donne un r�sultat ind�pendant du jet de d� pr�c�dent. Math�matiquement parlant, l’ind�pendance se d�finit par le fait que la probabilit� d’observer les �v�nements A et B est le produit des probabilit�s d’observer A et d’observer B.

4.1.3��Moyenne et �cart-type

Il y a en fait deux notions de moyenne, la premi�re consiste � faire la moyenne d’un �chantillon (n observations ind�pendantes) de la variable al�atoire, la deuxi�me, appel�e esp�rance de la variable al�atoire, est la moyenne de la valeur de la variable al�atoire pond�r�e par la probabilit� d’observer chacune des valeurs (intuitivement on pense s’en rapprocher lorsqu’on fait la moyenne d’un �chantillon assez grand).

=
�=E(X) = =�xi�P(X=xi)�(loi discr�te)�
�=E(X) =
x�f(x)�dx�(loi continue)�

Le lien entre les deux notions de moyenne est�:

L’esp�rance de la moyenne est �gale � la moyenne, lorsqu’on ne connait pas la loi d’une variable al�atoire, on peut estimer son esp�rance de mani�re non biais�e par la moyenne de l’�chantillon. Intuitivement, si les observations sont ind�pendantes et plus il y en a, plus la moyenne devrait s’approcher de l’esp�rance. (on verra avec le th�or�me central-limite qu’on en sait plus sur la distribution de la moyenne lorsque le nombre d’observations ind�pendantes pour calculer la moyenne est grand). Pour quantifier cela, il faut pouvoir �valuer l’�cart � la moyenne. La moyenne ne donne en effet pas d’informations sur les fluctuations de part et d’autre de la moyenne. Pour cela, on utilise l’�cart au carr� � la moyenne, dont on fait la moyenne puis on en prend la racine carr�e�: c’est l’�cart-type. On peut calculer l’�cart-type pour un �chantillon (par rapport � la moyenne de cet �chantillon ou par rapport � l’esp�rance de la loi si elle est connue) et l’�cart-type d’une loi (N.B. on appelle variance le carr� de l’�cart type).

�σ2(X)��(�chantillon) =
σ2��(loi discr�te) = E(�(X−E(X))2�) =�(xi−E(X))2�P(X=xi)��dx�
  = E(X2)−E(X)2�= =�xi2��P(X=xi)��−�E(X)2�
σ2��(loi continue) =
E(�(X−E(X))2�)� =� (x−E(X))2�f(x)��dx�
  =
E(X2)−E(X)2=� x2��f(x)��dx�−�E(X)2�

On observe que si m est une constante�:

C’est faux si m n’est pas constant, on a tout de m�me la relation

lorsque X et Y sont ind�pendants.

Contrairement au cas de la moyenne, il y a un ajustement � faire lorsqu’on estime l’�cart type d’une loi � partir d’un �chantillon si on ne connait pas la moyenne de la loi. En effet, on estime alors la moyenne de la loi � partir de la moyenne de l’�chantillon, ce qui conduit � sous-estimer la somme des carr�s des �carts � la moyenne. Lorsque les n valeurs de Xi de l�chantillon sont ind�pendants, on a�:

E(σ2� echantillon )=� �E(σ2(X))�

En effet�:

E(σ2(X)) =
  =
  =
�E(X2) −� �E(� �Xk2�+� �Xk�Xj�)�
  =
E(X2)�−� �(n�E(X2)�+�n(n−1)�E(X)2)�
  =
  =

On utilise donc

au lieu de

pour estimer l’�cart-type d’une loi � partir d’un �chantillon.

La d�finition de ces deux param�tres va nous permettre d’�noncer deux r�sultats sur la dispersion�:

  • σ�≠�0�⇒�P(|X− |�≥���σ)� ≤�
    En effet, on minore la variance de X en restreignant aux xi v�rifiant la condition ci-dessus�:
    σ2�≥� �(xi− )2 P(X=xi)�≥��2�σ2�P(|X− |�≥���σ)
  • Si on a n variables ind�pendantes, alors la variance de la somme de ces n variables est la somme des variances de chaque variable. En particulier si elles ont le m�me �cart-type σ, l’�cart-type de la moyenne est σ/√n.

La d�finition de la loi normale (cf. infra) permet de donner un sens quantitatif au fait qu’une moyenne d’observations ind�pendantes de m�me loi tend vers l’esp�rance de la loi (c’est le th�or�me central-limite qui sera admis).

4.1.4��La loi normale

Il s’agit d’une loi continue, d�finie par la fonction

o� la constante devant l’exponentielle est d�termin�e par la condition ∫−∞+∞ N�,σ2(x) dx=1.

On montre que la loi normale de param�tres � et σ2 a pour esp�rance � et variance σ2. Pour l’esp�rance, on fait le changement de variable x=� +y, on obtient

E(X)�=� �x e � =� (y+�)e =���+�� y e =��

la derni�re �galit� provenant de l’imparit� de la fonction � int�grer. Pour la variance, il faut faire une int�gration par parties pour enlever le terme en x2.

Plage de normalit� au niveau de confiance α ∈ ]0,1[: il s’agit d’un intervalle centr� autour de la moyenne I=]� − t σ, � +t σ[, tel que la probabilit� d’�tre dans I soit de α. Le calcul pratique se fait avec la fonction sp�ciale erf (error function) d�finie par

On montre par un changement de variable que t v�rifie l’�quation

Ainsi, comme erf(√(2))=0.954..., on utilise souvent la plage de normalit� au niveau de confiance de 95% qui est approximativement [�−2σ,�+2σ] (la probabilit� d’�tre dans cet intervalle est l�g�rement sup�rieure � 95%)

4.1.5��Th�or�me central-limite

Ce th�or�me dit que la variable al�atoire obtenue en faisant la moyenne de N variables al�atoires ind�pendantes et identiquement distribu�es (de m�me loi inconnue) suit une loi gaussienne de moyenne la moyenne de la loi inconnue et d’�cart type, l’�cart type de la loi inconnue divis� par √N. On estime alors la moyenne de la gaussienne par la moyenne M des observations, et l’�cart type S de la gaussienne par l’�cart type des observations divis� par √N. Pour une variable al�atoire suivant une loi gaussienne de moyenne m et d’�cart type σ, la probabilit� de tirer une valeur en-dehors de l’intervalle [m−2σ,m+2σ] est de 5%, on peut estimer que la moyenne “r�elle” pour le climat est dans l’intervalle [M−2S/√N,M+2S/√N]. Si S/√N est suffisamment petit, on a une d�finition raisonnable du climat. Plus N est grand, plus S/√N sera petit, mais d’une part on ne dispose pas de s�ries de mesure du temps sur des p�riodes tr�s grandes, d’autre part la loi d�finissant le climat d’un lieu �volue au cours du temps, on ne peut donc pas prendre N trop grand.

On peut faire le calcul de S/√N pour certaines stations m�t�o (les chiffres sont disponibles gratuitement sur Internet pour de nombreuses stations m�t�o am�ricaines), on trouve pour N=30 des valeurs de l’ordre de 0.1 degr�. Lorsqu’on moyenne sur le globe, les �carts-type diminuent au moins d’un ordre de grandeur, ce qui permet de parler de mani�re raisonnable d’augmentation de temp�rature (qui est largement sup�rieure aux incertitudes dues aux S/√N).

4.2��Corr�lation et r�gression

4.2.1��Corr�lation entre deux variables al�atoires.

Si deux variables X et Y sont ind�pendantes, alors

Dans le cas g�n�ral, cette relation n’est pas vraie, on d�finit la covariance entre X et Y par

C(X,Y)�=�E(X�Y�)�−�E(X)�E(Y)

Cette valeur n’est pas “normalis�e” (si X et Y ont des unit�s physiques par exemple ce nombre a une dimension), on d�finit donc la corr�lation par

ce nombre est sans unit�. On montre qu’il appartient � l’intervalle [−1,1]. Si Y=X alors c=1, si Y=−X alors c=−1. Si les variables sont ind�pendantes, alors c=0. Mais la r�ciproque est fausse. On utilise toutefois ce calcul pour estimer si deux variables sont corr�l�es ou ind�pendantes. Plus pr�cis�ment, si deux variables paraissent corr�l�es, on peut essayer d’exprimer l’une en fonction de l’autre, c’est ce qu’on va faire dans la section suivante.

4.2.2��R�gression lin�aire avec un mod�le lin�aire

Soit (Xi)1 ≤ i ≤ n et Yi deux s�ries statistiques. On cherche � savoir s’il est raisonnable de pr�voir Yi en fonction de Xi, avec par exemple une relation lin�aire Yi = a + b Xi. Pour cela on minimise f(a,b)=∑(a + b Xi −Yi)2 par rapport � a et b, ce qui donne les �quations :

�quivalentes � :









����(6)

Si on note U la matrice de premi�re colonne remplie de 1, et de deuxi�me colonne remplie par Xi (en ligne i),

et si on note U* la transpos�e de U,

on observe que

U*�U�=�















le syst�me (6) devient

U*�U�



=�















= U*�Y,����o� Y=







Donc la solution est donn�e par

A=



=�inv(U*�U)�(U*�Y)

En pratique avec Xcas, on place dans un tableur deux colonnes, une avec les Xi, une avec les Yi, on rajoute une colonne de 1 avant la colonne des Xi, on s�lectionne les colonnes de 1 et de Xi, on retourne dans l’historique (hist), on �crit U:= puis on copie (bouton du milieu de la souris ou touche copie) puis Entr�e. On revient dans mtrw, on s�lectionne la colonne des Yi, on repasse dans hist, on �crit Y:= et on recopie puis Entr�e. On effectue alors A:=inv(tran(U)*U)*(tran(U)*Y).

On peut ensuite juger de la r�gression en comparant f(a,b)/n avec la variance de Yi ou en comparant √f(a,b)/n avec l’�cart-type de Yi. On voit que f(a,b) est la somme des carr�s des composantes du vecteur U A − Y, donc le calcul pratique de f(a.b)/n se fait par la formule l2norm(U*A-Y)^2/size(Y), � comparer avec stddev(Y)^2. Plus ce rapport est proche de 0, plus on explique les variations de Yi par la r�gression lin�aire.

On peut montrer que les coefficients de la droite sont donn�s par�:

et

f(a,b)�=�n�(1−r2)�σ(Y)2,� ���r=

On montre que cette m�thode se g�n�ralise (dans sa formulation matricielle), si on veut pr�voir par exemple Yi � partir de deux s�ries statistiques xi et Xi avec une relation lin�aire de la forme Yi=a+bxi+cXi, alors (a,b,c) est obtenu en calculant inv(U* U) U* Y, o� cette fois la matrice U poss�de 3 colonnes, la colonne de 1, la colonne des xi et la colonne des Xi. Les calculs pratiques se font de mani�re identique.

4.2.3��Autres r�gressions lin�aires

On peut aussi supposer que Y d�pend de mani�re non lin�aire de X, par exemple logarithmique, exponentielle, polynomiale, etc. mais toujours de mani�re lin�aire des param�tres � d�terminer et calculer ces param�tres en utilisant la m�me m�thode (minimiser l’erreur quadratique). Ceci peut servir par exemple � mod�liser la consommation d’une denr�e non renouvelable (par exemple les combustibles fossiles) par une fonction logistique ou par une gaussienne.

4.2.4��Corr�lations entre gaz � effet de serre et temp�rature

Si on fait une corr�lation temp�rature-CO2 sur les 400 000 derni�res ann�es, on trouve un rapport f(a,b)/n sur variance de 0.35, soit une qualit� de corr�lation de 0.65 (65% de la variance de la temp�rature est expliqu�e par le taux de CO2) et un coefficient d’augmentation de temp�rature de 9 degr�s pour 100 ppm de CO2. On peut aussi faire une corr�lation temp�rature en fonction du CO2 et du CH4, on trouve alors 6 degr�s pour 100 ppmv de CO2 et 1.5 degr� pour 100 ppbv de CH4 avec une qualit� de corr�lation de 0.68.

Bien entendu, une corr�lation statistique ne suffit pas � conclure qu’une des variables d�pend (au moins partiellement) de l’autre, mais ici la physique nous montre que l’augmentation du taux de CO2 augmente les radiations solaires captur�es au niveau du sol, la corr�lation permet de quantifier cet effet (calcul que l’on pourra comparer � celui obtenu par d’autres m�thodes).

5��Mod�lisation num�rique (climat, �nergie)

Dans cette section, on pr�sente diff�rents types de mod�les math�matiques qui id�alisent la situation r�elle. Les principaux types s’appliquant ici sont les mod�les discrets, les mod�les continus se ramenant � des �quations ou syst�mes diff�rentiels et les mod�les de type �quations aux d�riv�es partielles (qui ne seront pas pr�sent�s ici).

5.1��Exemple�: un mod�le tr�s simplifi� du climat

On suppose ici que la Terre est assimilable � un oc�an (pas de glace modifiant l’alb�do) recevant un rayonnement solaire constant et �mettant un rayonnement dont la puissance est compos�e d’une partie proportionnelle � T4 (comme un corps noir) dont on retranche le pi�gage du au gaz carbonique (proportionnel au log de la concentration de CO2). La diff�rence entre les deux rayonnements incident et renvoy�s sert � r�chauffer l’eau. On obtient l’�quation (o� le coefficient 6 et le ln pour le taux de CO2 proviennent de calculs de physique)

�=�k�(�6�ln(CO2/280)�−�σ�(T4−T04)�)

Le calcul de k donne (86400 � 365)/(3e6 � 4.18e3) =0.0025 degr�s par an (on chauffe l’eau sur 3000 m�tres en moyenne avec une puissance de 1W/m2). On peut prendre T0=288 (15 degr� celsius), et pour σ la constante de Stefan-Boltzman (5.67e−8, ceci aura tendance � sous-estimer le ph�nom�ne car l’eau de mer n’est quand m�me pas un corps noir). On peut prendre par exemple une variation du taux de CO2 constante et �gale � 2 ppm par an, soit si t est exprim�e en ann�es

On ne sait pas r�soudre exactement ce type d’�quation diff�rentielle, on utilise alors des m�thodes de r�solution approch�e. La plus simple consisterait dans notre cas � supposer la temp�rature du globe constante pendant une ann�e et � calculer la diff�rence entre l’ann�e en cours et l’ann�e suivante

� T(t+1)=T(t)�+�k�(�6�ln(CO2/280)�−�σ�(T4−T04)�) ����(7)

Il s’agit de la m�thode d’Euler (avec un pas de 1 an), graphiquement cela revient � tracer le graphe de la solution d’une �quation diff�rentielle en suivant les tangentes.

L’avantage de cette discr�tisation est qu’on peut mod�liser de mani�re plus r�aliste la variation du taux de CO2.

5.2��Exemple�: sc�narios d’�mission de CO2

5.2.1��CO2 et consommation d’�nergie

L’homme �met dans l’atmosph�re des gaz � effet de serre principalement par combustion des p�trole, gaz, charbon, �galement par la d�forestation et les �missions de m�thane de l’agriculture. Un petit calcul �l�mentaire permet d’�tablir la correspondance entre 1 ppm de CO2 (une partie de CO2 par million en volume) et 2 Gt (gigatonnes) de carbone.

Le m�tre a �t� pendant longtemps d�fini comme 1/40 millioni�me de la circonf�rence de la Terre, on a donc 2π R=4e7 m o� R d�signe le rayon de la Terre. On en d�duit la surface de la Terre

4�π�R2=�(4e7)2/π=5e14�m2�

La pression de l’air au sol est �quivalente � une colonne d’eau de 10m (ou encore 76cm de mercure dont la densit� est 13). La masse de l’atmosph�re terrestre est donc �gale � celle de 5e15m3 d’eau soit 5e18kg. L’air est pour trois quarts compos� d’azote (N2 masse molaire 214=28), et pour le reste d’oxyg�ne (O2 masse molaire 216=32), alors que le carbone C a une masse molaire de 12. Comme il y a 1 atome de carbone par mol�cule de CO2, on en d�duit que la masse de carbone correspondant � une concentration de 1ppm de CO2 est de�:

�*�5e18�*�1e−6=�2.1�e12�kg�

soit 2.1Gt de carbone.

La consommation mondiale de combustibles fossiles s’�tablit pour 2003 � 3.5 Gtep (giga-tonne �quivalent p�trole) de p�trole, 2.3 Gtep de gaz naturel et 2.4 Gtep de charbon. Avec un taux moyen en tonne de carbone par tep de 0.856 pour le p�trole, 0.651 pour le gaz et 1.123 pour le charbon, on a �mis

3.5�נ0.856�+�2.3�נ0.651�+�1.123�נ2.4 =�7.1�Gt�C

ce qui correspond � 3.5 ppm de carbone. Auquel il faut ajouter le carbone �mis par d�forestation, retrancher le CO2 absorb� par l’oc�an et les �cosyst�mes, on a mesur� une augmentation de 2ppm.

5.2.2��Sc�narios de consommation d’�nergie et d’�missions de CO2

Les incertitudes sur la sensibilit� du climat � une augmentation fix�e des taux de gaz � effet de serre sont assez faibles. Par contre les pr�visions de taux de gaz � effet de serre sont tr�s incertaines puisque par exemple la concentration de CO2 � la fin du si�cle d�pend de la quantit� de combustibles fossiles que nous aurons brul�e. Celle-ci est tr�s incertaine, car�:

  • nous ne connaissons pas pr�cis�ment la quantit� de p�trole, gaz, charbon � notre disposition. Si comme certains g�ologues le pensent, nous avons consomm� la moiti� des p�troles conventionnels, la consommation de p�trole va diminuer dans les ann�es � venir, au profit d’autres �nergies fossiles plus ou moins polluantes en CO2 par tep, d’�nergies renouvelables ou d’�conomies d’�nergie.
  • nous ne connaissons pas l’avenir des fili�res nucl�aires
  • nous ne pouvons pas pr�dire les conditions �conomiques et d�mographiques qui pr�vaudront dans 20 ou 50 ans

Plusieurs mod�les existent donc, par exemple les sc�narios du GIEC ou d’experts (par exemple ceux d�riv�s des projections de l’ASPO). Quelques principes math�matiques sous-tendent ces mod�les�:

  • Des mod�les de type suite g�om�trique lorsque la croissance n’est pas contrainte ou pour certains mod�les simples de d�pletion d’une quantit� (on suppose que la quantit� donn�e va �voluer comme une suite g�om�trique de raison >1 dans le premier cas, de raison < 1 dans le deuxi�me cas).
  • Des mod�les malthusiens lorsqu’on s’int�resse � une ressource dont la quantit� va limiter l’accroissement (de population, de consommation, ...), par exemple logistique�: on suppose que la variation relative de la quantit� est lin�aire, avec un coefficient directeur n�gatif pour limiter la croissance. Ce type de mod�le est utilis� par exemple pour estimer la quantit� ultime r�cup�rable de p�trole dans le monde.
  • Des mod�les lin�aires lorsque plusieurs quantit�s sont en interaction, on suppose que la variation de consommation d’une quantit� donn�e d�pend lin�airement de toutes les quantit�s. On r�sout alors un syst�me diff�rentiel lin�aire.
  • etc.

5.3��Mod�les discrets.

Les mod�les les plus utilis�s sont les suites arithm�tiques, les suites g�om�triques, et plus g�n�rallement les suites r�currentes. Historiquement, c’est Malthus qui soutint la th�se que la production croit comme une suite arithm�tique, alors que la population croit comme une suite g�om�trique, pr�disant de gros probl�mes lorsque la suite g�om�trique d�passait la suite arithm�tique. C’est Verhulst qui proposa alors le mod�le logistique pour am�liorer la mod�lisation (cf. infra dans le cas continu).

Plus g�n�rallement, les suites r�currentes un+1=f(un) permettent de mod�liser des ph�nom�nes o� la relation ne d�pend pas explicitement du temps, il est alors int�ressant d’�tudier les positions d’�quilibre aussi appel�s points fixes, et leur stabilit� (si u0 est proche d’un point fixe l, la suite un converge-t-elle vers l�?). On peut aussi �tudier la sensibilit� de l’�quilibre � une variation des param�tres d�finissant la fonction f (d�placement de l’�quilibre, passage d’un �quilibre stable � un �quilibre instable). Le r�sultat principal sur la stabilit� est que si |f′(l)|<1 alors l’�quilibre est stable (on peut le comprendre intuitivement en faisant une repr�sentation graphique de la suite en dimension 1).

Si f d�pend explicitement du temps (i.e. de n), par exemple pour (7), il n’y a pas de r�sultats simples. On peut r�soudre explicitement quelques cas particuliers, par exemple un+1=Aun+P(n), o� A est une matrice constante et P un polyn�me.

5.4��Mod�les continus�: �quations et syst�mes diff�retiels

On a vu dans les sections pr�c�dentes que de nombreux mod�les se ramenent � la r�solution d’�quations diff�rentielles ou de syst�mes diff�rentiels�:

o� X est la quantit� recherch�e et t le temps. La plupart du temps ces �quations et syst�mes ne peuvent �tre r�solues que de mani�re num�rique, mais certains cas particuliers importants peuvent �tre r�solus exactement. Nous verrons le cas du mod�le logistique comme exemple d’�quation r�soluble (qui est utilis� comme mod�le de production d’�nergie non renouvelable) et les syst�mes diff�rentiels affines (qui aident � comprendre une am�lioration du mod�le tr�s simple, 5.1, tenant compte des interactions entre temp�rature, CO2 et alb�do). Pour les cas non r�solubles, on essaie de d�terminer le comportement qualitatif de l’�quation � partir de cas connus, en g�n�ral en lin�arisant. On mettra en �vidence sur ces syst�mes la notion d’�quilibre stable et instable, et la notion de r�troaction.

5.4.1��Le mod�le logistique

Le mod�le exponentiel, historiquement propos� par Malthus et encore tr�s (trop?) utilis� par les �conomistes, correspond � l’�quation diff�rentielle�:

o� y(t) d�signe la quantit� totale extraite au temps t. La solution de cette �quation est une exponentielle qui tend vers +∞ lorsque t tend vers l’infini, ce qui n’est pas possible pour une ressource non renouvelable. Pour avoir un mod�le plus r�aliste qui tienne compte d’une ressource finie, Verhulst (1840) passe d’un taux d’accroissement constant � un taux lin�aire et qui diminue lorsque y augmente

�=�−ay+b�=−a(y−c),����a>0,�c=

La solution de cette �quation est appel�e fonction logistique. On s�pare les variables�:

on d�compose en �l�ments simples�:

a�dt�=�


�−�


soit�:

b�dt�=��


�−�


d’o� en int�grant�:

puis, en enlevant les valeurs absolues et en posant K=� eC�:

Donc K>0, on peut poser K=e−bt0 et finalement�:

On peut �tudier le comportement de la fonction y en fonction du temps, en t=−∞, la limite de y est nulle, la fonction croit pour t n�gatif comme une exponentielle

(avec le m�me taux de croissance que lorsque a=0) et est convexe. Pour t=t0, y=c/2, la fonction a un point d’inflexion et devient concave pour rejoindre sa limite c en t=+∞. Si on calcule y′ qui mod�lise la production instantann�e, on obtient�:

On observe que y′ a un graphe en forme de cloche avec deux limites nulles aux extr�mit�s et un maximum en t=t0�:

Modèle mathématique réchauffement climatique

En pratique, pour mod�liser y on effectue une r�gression lin�aire du quotient de la production annuelle par la production cumul�e (par exemple pour le p�trole depuis 1900) en ordonn�e (ce qui repr�sente dy/ydt) en fonction du temps. Un grand int�r�t de ce mod�le est qu’il permet de pr�dire la quantit� totale que l’on pourra extraire. Sa principale faiblesse est la sensibilit� de cette estimation en fonction de la qualit� des donn�es initiales et des irr�gularit�s (par exemple chocs p�troliers). Pour la production de p�trole, on peut un peu l’am�liorer en mod�lisant plusieurs cycles de p�trole donc en additionnant plusieurs courbes logisitiques (p�trole on-shore facile, p�trole offshore, deepwater, ultra-deepwater, polaire, p�troles lourds).

5.4.2��Syst�mes diff�rentiels lin�aires

On s’int�resse plus g�n�rallement � des syst�mes affines:

o� X(t)=(x1,...,xn) est un vecteur de dimension n, A une matrice n,n et C=(c1,...,cn) est un vecteur constant Plus g�n�rallement, on peut r�soudre exactement ce type de syst�me lorsque C ne d�pend que du temps.

Pour les syst�mes qui nous int�ressent, on suppose en g�neral que les coefficients sur la diagonale sont n�gatifs, de sorte que si les interactions entre param�tres �taient nulles, la solution du syst�me convergerait. Pour une mod�lisation r�aliste, A n’est pas diagonale, les coefficients non diagonaux traduisent des r�troactions entre les param�tres, ces r�troactions peuvent �tre positives si ces coefficients sont positifs (l’interaction amplifie l’action de la variation originale du param�tre) ou n�gatives dans le cas contraire.

Si la matrice A est inversible, il existe une unique position d’�quilibre, en effet si dX/dt=0, on a X=Xe=−A−1 C et r�ciproquement cette valeur de X est solution de (8). Pour une condition initiale X0 quelconque, on peut alors se demander si X(t) va se rapprocher de la position d’�quilibre Xe (�quilibre stable) ou non (instabilit� si on s’�loigne de l’�quilibre). Si le syst�me est stable, on peut aussi se demander comment la position d’�quilibre varie lorsque C varie (par exemple pour les param�tres du climat si le for�age radiatif est modifi� par variation des param�tres orbitaux de la Terre).

Il ne faut pas confondre la question de la stabilit� du syst�me et l’existence de r�troactions positives ou n�gatives. Un syst�me avec des r�troactions positives peut tr�s bien �tre stable.

On se limitera ici au cas de matrices A r�elles et diagonalisables (sur ℂ). Soient λ1,...,λn les valeurs propres de A (r�elles ou par paires de complexes conjugu�s) et P la matrice de passage vers une base propre correspondante

P−1�A�P=D�=�diag(λ1,...,λn)

On pose X=PY (Y=P−1X), on a alors

=P−1 �=�P−1(AX+C)��=�P−1APY�+�P−1C =�D�Y�+�P−1C�

Soit B=P−1C, on est ainsi ramen� � n �quations diff�rentielles lin�aires d’ordre 1 � coefficients constants�:

Lorsque bk=0, l’�quation admet pour solution g�n�rale

En faisant varier la constante c=c(t), et en rempla�ant dans l’�quation v�rifi�e par yk, on obtient�:

c′�eλkt�=�bk��⇒��c′�=�e−λkt��bk�

donc c=∫e−λkt bk, et si bk est constant,

o� la valeur de K se calcule en posant t=0, finalement la solution est donn�e par�:

yk(t)=−bk/λk�+�(yk(0)+ )�e−λkt�

c’est la somme d’un terme constant (k-i�me coordonn�e de la position d’�quilibre dans la base propre) et d’un multiple d’une exponentielle dont le comportement lorsque t tend vers +∞ permet de d�terminer si le syst�me est stable.

Si ℜ(λk)>0 pour au moins une des valeurs de k, alors la composante correspondante yk(t) tend vers l’infini, et le syst�me est instable. Si ℜ(λk)=0, la composante reste born�e, il n’y a pas convergence vers l’�quilibre mais le syst�me fait des tours autour de la position d’�quilibre. Si ℜ(λk)<0 pour toutes les valeurs de k, il y a convergence (exponentiellement rapide) vers l’�tat d’�quilibre, le syst�me est stable, la vitesse de convergence est caract�ris�e par l’inverse de la valeur propre dont la partie r�elle est la plus proche de 0.

On prend souvent comme exemple de cette th�orie le mod�le proie-pr�dateur, on peut ici dans le contexte �nerg�tique proposer un mod�le baril (prix)-demande-production (sans contrainte g�ologique). Plus le prix du baril est �lev�, plus la production augmente et la demande baisse. Plus la production est �lev�e, plus le baril baisse. Plus la demande est �lev�e, plus le baril monte. On est donc en dimension 3.







=











+�C

On remarque que 0 est valeur propre, la trace de la matrice est nulle, il y a 2 valeurs propres imaginaires pures � i √bd+ac conjugu�es et le syst�me va effectuer des cycles (dont la p�riode est 2π/√bd+ac).

5.4.3��Les mod�les non lin�aires

Un syst�me diff�rentiel g�n�ral d’ordre 1 peut s’�crire

Si F ne d�pend pas explicitement du temps (syst�r�me autonome), on peut commencer par chercher les positions d’�quilibre possibles, i.e. les solutions de F(X)=0, puis si Xe est une solution, lin�ariser, c’est-�-dire r�soudre

pour voir si en partant de conditions initiales proches de l’�quilibre, le syst�me s’en approche (stable) ou s’en �loigne (instable).

Un autre cas est int�ressant, celui d’une perturbation d�pendant du temps

par exemple pour �tudier la perturbation introduite par l’ajout de CO2 � un syst�me lin�aire mod�lisant le syst�me climatique non perturb�. Dans le cas du CO2, on peut par exemple supposer un taux d’�mission constant (p(t)=p) et voir comment l’�quilibre du syst�me climatique est affect�. Plus g�n�rallement, pour trouver l’�quilibre on souhaitera trouver Xe(t) tel que

c’est le th�or�me des fonctions implicites qui permet de d�terminer Xe(t) en fonction de t, sous r�serve que ∂X G soit inversible. Si le syst�me diff�rentiel lin�aris� correspondant � G(Xe(t)) est stable et que la perturbation p(t) varie lentement par rapport aux constante de temps caract�risant la stabilit� ou est d’amplitude faible, alors la position d’�quilibre va suivre. Pour le syst�me climatique soumis aux variations externes des param�tres orbitaux, on peut consid�rer que c’est le cas. Par contre ce n’est pas vrai pour la perturbation anthropique du CO2, on peut toutefois �tudier le syst�me lin�aris�, en supposant que G′(Xe) reste � peu pr�s constant.

5.4.4��Exemples de mod�les de ce type

On consid�re 3 param�tres pour mod�liser le climat, la temp�rature moyenne T, la masse de glace G (qu’on suppose reli� � l’alb�do par la surface de glace) et le taux de CO2 dans l’air C. La temp�rature moyenne influe sur la puissance rayonn�e (proportionnelle � T4), sur le CO2 dissous dans l’oc�an (si T augmente, C augmente) et sur la glace (si T augmente, G diminue). Le CO2 influe sur la puissance rayonn�e en fonction de −ln(C) et sur le taux de CO2 dissous dans l’eau de mer. La glace influe sur la puissance rayonn�e, en fonction de la surface qu’on peut mod�liser par G(2/3). Enfin la puissance rayonn�e est reli�e � la temp�rature par la capacit� calorifique de l’oc�an.

Il y a �videmment encore beaucoup de simplifications dans un tel mod�le � 3 param�tres, la temp�rature de l’oc�an n’est pas homog�ne et la temp�rature des eaux profondes ne varie que tr�s lentement, il y a d’autres gaz � effet de serre, ...

On a alors







= F�=�


cc�(σ(T04−T4)+α�ln(C)�−�β�G2/3�)�
f(T)�
g(C,T)



o� cc est proportionnel la capacit� calorifique de l’oc�an, f est une fonction d�croissante, ∂T g est positif.

On peut aussi utiliser ce type de mod�le pour obtenir le pourcentage de CO2 �mis par l’homme restant dans l’atmosph�re au bout de t ann�es

CO2(t)=18+26�e−t/3.4�+�24�e−t/21�+�18�e−t/70�+�14�e−t/420

5.5��Les mod�les de type EDP

Les fluides de l’atmosph�re ob�issent aux �quations physique de la dynamique des fluides. Contrairement au probl�me � 2 corps, on ne peut pas r�soudre ces �quations de mani�re exacte. On utilise alors des m�thodes num�riques, en discr�tisant l’�tat de l’atmosph�re. On obtient alors une solution approch�e des �quations, c’est le principe des pr�visions num�riques m�t�orologiques. Ce principe s’applique aussi au climat. On peut s’�tonner que le m�me principe qui ne permet pas des pr�vision fiables � plus de quelques jours puisse permettre de faire des pr�visions � l’�ch�ance du si�cle. La contradiction n’est qu’apparente, parce que les mod�les ne pr�voient pas le temps pr�cis un jour donn�, mais le temps moyen, en fonction des grandes contraintes �nerg�tiques (ensoleillement, effet de serre) de m�me qu’on peut pr�voir qu’il fera en moyenne plus chaud en �t� qu’en hiver parce que l’ensoleillement est plus important en �t� qu’en hiver. L’int�r�t des mod�les est de pouvoir quantifier le “il fera plus chaud si l’effet de serre est renforc�”.

Les probl�mes math�matiques li�s � la mod�lisation du climat sont li�s � la stabilit� des m�thodes de r�solution num�rique. On peut en avoir une petite id�e en �tudiant les m�thodes num�riques d’int�gration ou la m�thode d’Euler pour r�soudre des �quations diff�rentielles de type y′=f(x,y) (mais ici on a 3 dimensions d’espace et 1 de temps et beaucoup plus qu’une variable, il y a la pression, la temp�rature, l’humidit�, ...). On peut aussi s’inspirer d’exemple de suites it�ratives pour essayer de comprendre certains ph�nom�nes.

6��Energies renouvelables.

6.1��Le solaire thermique

Le calcul du rendement d’un panneau solaire thermique d�coule du calcul de l’ensolleillement fait au th�or�me 1, mais il faut tenir compte du fait que la normale n au panneau solaire n’est plus la verticale locale. On exprime les coordonn�es de n dans le rep�re local orthonorm� form� de la verticale locale, du Sud local, et de la direction obtenue par produit vectoriel des deux vecteurs pr�c�dents. La verticale locale v est connue (cf. section 2.1), le vecteur norm� port� par la direction du Sud local est �gal � l’oppos� de la d�riv�e de v par rapport � la latitude, apr�s normalisation.

Lorsque l’angle avec la verticale est fix�, on a bien sur int�r�t � aligner le panneau avec le Sud local. On peut calculer l’angle optimal pour recevoir le maximum d’�nergie solaire tout au long de l’ann�e, mais il est plus int�ressant de calculer l’angle optimal pour maximiser le rayonnement re�u pendant la p�riode de chauffe pour minimiser l’utilisation d’un appoint (par exemple sur la p�riode automne-hiver, on trouve un angle voisin de 60 degr�s � nos latitudes sans tenir compte de la m�t�o). On peut aussi chercher � optimiser l’angle au solstice d’hiver si on souhaite �viter l’utilisation d’un appoint (on utilisera par contre plus de surface de panneaux, ceux-ci devant alors �tre toujours exc�dentaires).

6.2��Le solaire � concentration.

Il s’agit de g�n�rer de la chaleur (four solaire) ou du courant �lectrique comme avec une centrale thermique classique, la chaleur �tant fournie par les rayons solaires au lieu de la combustion de charbon, gaz ou fuel. Les rayons solaires ne sont bien sur pas suffisamment chauds pour obtenir une chaleur ou un rendement suffisant sans concentration. La concentration se fait � partir d’un ou de miroirs qui r�fl�chissent les rayons en une zone la plus petite possible.

Etude id�alis�e�: on suppose que le miroir peut prendre la forme d’une courbe r�guli�re (ou que les miroirs sont suffisament petits pour �tre assimilable � des tangentes � une courbe). Quelle forme doit prendre la courbe pour focaliser les rayons solaires incidents en un seul point�? Une r�ponse est connue depuis longtemps, il s’agit d’une parabole (ou plus pr�cis�ment d’un paraboloide de r�volution) qui focalise les rayons perpendiculaires � la directrice de la parabole vers le foyer de la parabole. Montrons qu’effectivement cette contrainte impose une parabole, en effectuant une coupe par un plan passant par le foyer.

Rappel�: une parabole est d�finie g�om�triquement par une droite D appel�e directrice et un point O appel� foyer, c’est le lieu des points �quidistants � D et O.

Modèle mathématique réchauffement climatique

Si on calcule l’�quation de ce lieu, en prenant O comme origine, M(x,y) et l’�quation de D sous la forme y=−m, on a

x2+y2=OM2=d(M,D)2=(y+m)2=y2+2ym+m2

donc

c’est bien l’�quation usuelle d’une parabole. La normale � la parabole est obtenue en prenant le gradient de l’�quation qui la d�finit

c’est donc

qui est le vecteur directeur de la bissectrice de la verticale et de OM (plus pr�cis�ment la bissectrice int�rieure du vecteur OM et de −j ou de MO et de j). La tangente est donc la bissectrice ext�rieure de cet angle, ou encore la bissectrice int�rieure de (MO,MK). On en d�duit qu’un rayon vertical venant du haut se r�fl�chit sur la parabole sym�triquement � la normale donc en passant par le foyer O.

Consid�rons maintenant une courbe ayant cette propri�t�: les rayons provenant du haut verticalement se r�fl�chissent dessus en passant par O. Soit M(x,y) un point de cette courbe. Prolongeons le rayon vertical passant par M vers le bas de la longueur OM, on obtient un point K d’ordonn�e

Montrons que lorsque x varie, yk a une d�riv�e nulle. On a�:

D’autre part, la tangente en M � la courbe est orthogonale � OK puisque c’est la bissectrice de l’angle (MO,MK) (hypoth�se que le rayon r�flechi passe par O) et que MO=MK (par construction de K), et le vecteur directeur de la tangente est (1,y′), on a donc

(1,y′).(x,yk)=0�⇒�x+y′y−y′ =0

donc yk′ est bien nul et K est donc sur une droite horizontale fixe D, d’o� l’on d�duit que M est sur la parabole de foyer O et directrice D, la courbe est donc forc�ment une parabole.

Rendement d’un tel syst�me�: on peut le d�finir comme le rapport de l’�nergie focalis�e sur la surface de miroir. Pour un miroir cylindrique avec une coupe parabolique, la surface est proportionnelle � la longueur de la courbe en coupe. On peut faire le calcul sur la parabole y=x2/2 pour d�terminer par exemple quel portion de parabole on peut raisonnablement utiliser. Le cout de l’�nergie est proportionnel �

�dx� =�


ln( −l)


Il ne faut pas prendre l trop grand (sinon le cout augmente trop vite), mais pas trop petit non plus (sinon on ne concentrera pas assez d’�nergie pour avoir par exemple un bon rendement de Carnot). On peut aussi jouer sur le param�tre de la parabole (pour l’aplatir pr�s de son point le plus bas, ce qui permet d’augmenter l � rendement constant), mais il faut alors augmenter l’altitude du foyer.

Ce syst�me est bien adapt� s’il reste petit, par exemple pour un four solaire pour faire la cuisson d’un repas, mais n’est pas adapt� � de grandes installations, la taille du miroir parabolique � r�aliser �tant trop grande et trop difficile � manoeuvrer (toute l’installation doit bouger pour que la directrice de la parabole reste perpendiculaire aux rayons solaires). Pour une centrale solaire � concentration, le foyer doit �tre � une position fixe, et les miroirs doivent pouvoir �tre inclin�s pour faire suivre les rayons du Soleil vers le foyer au cours de la journ�e. La direction miroir-foyer est donc fixe, et la direction miroir-Soleil impos�e par la position du Soleil, le plan du miroir est d�termin� par le point du miroir, la bissectrice ext�rieure aux 2 directions (dans le plan les contenant) et la perpendiculaire � ce plan. Ceci impose bien d’avoir des miroirs amovibles. L’�nergie capt�e par ce syst�me est proportionnel au cosinus de l’angle entre la normale au miroir et la direction miroir-Soleil. Comme la normale au miroir est la bissectrice int�rieure de l’angle miroir-Soleil, miroir-foyer, l’�nergie capt�e (par miroir) est proportionnelle au cosinus de la moiti� de l’angle miroir-Soleil, miroir-foyer. On est donc ramen� � un calcul analogue � celui du rendement de panneaux solaires thermiques, la normale au panneau �tant remplac�e par la direction miroir-foyer, mais on divise ensuite l’angle par 2.

6.3��Les aspects �conomiques.

Il y a actuellement deux grands freins au d�veloppement des �nergies renouvelables�: l’intermittence (conjointement avec la possibilit� de stocker l’�nergie) et la facon dont le co�t est calcul�. Le premier frein n�cessite des avanc�es technologiques pour pouvoir monter en puissance, le deuxi�me ne me semble pas intrins�que mais uniquement fond� sur la facon de financer un projet �nergie renouvelables. En effet, un projet d’installation de panneaux solaires ou d’�oliennes n�cessite un investissement important au d�part (qui peut repr�senter jusqu’� plus de 90% du co�t total par exemple pour des panneaux solaires thermiques) contrairement � un projet de centrale thermique au charbon ou au gaz. Comment doit-on r�partir ce cout de mani�re juste sur la dur�e de fonctionnement de l’installation�? Les �conomistes utilisent un taux d’actualisation, ils comparent avec ce qu’une somme d’argent �quivalente aurait pu rapporter en euros constants. Plus pr�cis�ment le cout annuel estim� correspond � la moyenne des couts de fonctionnement plus le cout annuel qu’il faudrait payer pour rembourser en euros constants l’investissement initial. Le choix du taux d’actualisation est alors un param�tre tr�s important qui peut faire varier du simple au double (voire plus) le co�t estim� des �nergies renouvelables. Les taux actuels de 5% � 10%, calcul�s sur la base d’une croissance soutenue de l’�conomie mondiale rendue possible par la hausse de la consommation mondiale d’�nergie fossiles et l’am�lioration de l’efficacit� �nerg�tique, paraissent peu compatibles avec une p�riode de “plateau fossiles” voire de d�pl�tion.

7��Conclusions

Les variations de concentration de gaz � effet de serre au cours du pass� agissent sur le bilan �nerg�tique de la Terre de mani�re logarithmique de l’ordre de quelques W/m2 (Watt par m�tre carr�), une valeur qui est du m�me ordre de grandeur que celles dues aux variations orbitales de la Terre ou que celles dues aux variations des volumes de glace (modification de l’alb�do terrestre). La corr�lation entre gaz � effet de serre et temp�rature calcul�e sur les 400 000 derni�res ann�es en Antarctique donne une sensibilit� de 9 degr�s pour 100 ppm de CO2 au niveau de l’Antarctique. En divisant ce chiffre par 2 pour tenir compte de la plus grande variabilit� de la temp�rature en Antarctique puis par 3 pour tenir compte des 3 causes principales de variation de m�me importance (gaz � effet de serre, variations orbitales, volume des glaces), on obtient une sensibilit� empirique de 1.5 degr� pour 100ppm de CO2. Ce qui veut dire que si nous conservions la concentration actuelle en CO2 (100 ppm de plus qu’avant le d�but de l’�re industrielle en 2005), ce mod�le pr�voit une hausse de 1.5 degr�s de la temp�rature, ce qui est parfaitement plausible avec l’�l�vation de 0.6 degr�s du 20�me si�cle (l’oc�an jouant le role de r�gulateur thermique de par son inertie sur des �chelles de temps entre le si�cle et le mill�naire, par exemple le for�age radiatif actuel est estim� � plus de 2W/m2, or pour r�chauffer l’oc�an de 1K avec une telle puissance il faut attendre plus de 200 ans, pour faire ce calcul, on consid�re qu’il y a une �paisseur moyenne de 3000 m�tres d’eau, donc chaque heure, 2Wh servent � r�chauffer 3000 m�tre cubes d’eau, or il faut un peu plus d’un kWh pour r�chauffer un m�tre cube d’eau d’un degr�, il faut donc environ 1,5 millions d’heures pour r�chauffer l’oc�an � cette puissance).

Bien �videmment, une corr�lation seule ne suffit pas � justifier l’existence d’une relation de cause � effet entre gaz � effet de serre et temp�rature, ce sont les lois de la physique qui justifient cette relation, la corr�lation nous donne toutefois des informations quantitatives sur l’effet que l’on peut rapprocher des calculs issus de la physique. Ceux-ci donne un for�age radiatif de 4W par doublement du CO2 soit 0.3% de l’�nergie re�ue au sommet de l’atmosph�re perpendiculairement au Soleil, ou encore 1.2% de l’�nergie re�ue au sommet de l’atmosph�re en moyenne (pour tenir compte de la surface de la Terre 4π R2 qui est 4 fois plus grande que sa section perpendiculaire au Soleil π R2). En tenant compte de la partie du rayonnement arrivant effectivement dans les basses couches de l’atmosph�re, c’est un for�age de 2 � 3%. Si la Terre �tait un corps noir � une temp�rature de 288K, ce for�age se traduirait par une hausse de plus de 0.5% de la temp�rature pour compenser (car la puissance �mise par un corps noir est proportionnelle � T4, il faut donc diviser par 4 le for�age), donc de 2K environ. Ce chiffre est raisonnablement en accord avec les mod�les du climat qui pr�voient une hausse de 3 � 4K pour un doublement du CO2, en tenant compte des diverses r�troactions (modification de l’alb�do par exemple).

Si nous voulons �viter une variation trop brusque de la temp�rature et plus encore des pr�cipitations, il faut donc agir d�s maintenant, ce qui ne sera pas facile car 1ppm de CO2 correspond � 2Gt de C soit � peine plus de 300kg de C par Terrien (ou un peu plus de 400 litres de p�trole)�: rappelons qu’un Fran�ais consomme (en gros pour moiti� directement par les transports et le chauffage, et l’autre moiti� indirectement par la consommation de biens) en moyenne 1.5 tep de p�trole, 0.67 tep de gaz et 0.2 tep de charbon par an, il �met donc 2 tonnes de C par an (ce qui g�n�ralis� � l’�chelle de la plan�te correspondrait � une �mission de 6 ppmv de CO2 par an).

8��R�f�rences sur le web

  • Astronomie
    • calcul de l’insolation en fonction de la latitude au cours des mill�naires �coul�s (Laskar et al)
      http://www.imcce.fr/fr/presentation/equipes/ASD/insola/earth/earth.html
    • La th�orie de Milankovitch: rechercher Milankovitch sur google par exemple sur le site de l’ENS Lyon
      http://www.ens-lyon.fr/Planet-Terre/Infosciences/Histoire/Paleoclimats/
  • Climat
    • IPCC http://www.ipcc.ch/
    • Realclimate www.realclimate.org
    • Reconstructions climatiques (nombreux liens)
      http://www.ncdc.noaa.gov/paleo/recons.html,
    • Donn�es num�riques de Vostok
      ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/vostok/, le fichier ch4nat.txt contient les donn�es pour le m�thane, co2nat.txt pour le CO2, deutnat.txt pour la temp�rature. (remonter d’un cran dans l’arborescence du site pour d’autres lieux).
    • donn�es du projet EPICA les plus r�centes
      ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/epica_domec/
    • Donn�es climatiques au 20 i�me si�cle
      http://data.giss.nasa.gov/
      http://hydrolab.arsusda.gov/nicks/nicks.htm
  • �nergie, �mission de gaz � effet de serre:
    • Concentration de CO2 � Mauna Loa (Hawai)
      http://cdiac.esd.ornl.gov/ftp/trends/co2/maunaloa.co2
    • Statistiques du gouvernement fran�ais
      http://www.industrie.gouv.fr/energie/sommaire.htm
    • The Oil Drum http://www.theoildrum.com/
    • ASPO-France http://aspofrance.org/,
    • Institut Francais du p�trole (IFP) http://www.ifp.fr
    • site de J.M. Jancovici (�nergie et effet de serre), http://www.manicore.com/
  • En compl�ment, vous pouvez aussi faire des recherches sur wikipedia, google et d’autres moteurs de recherche.

A��Pr�cession.

(D’apr�s Fr�d�ric Faure). La rotation de la Terre a pour cons�quence qu’elle n’est pas sph�rique. En effet, la force centrifuge vaut mω2 r, o� r est le vecteur issu de la projection sur l’axe de rotation et ω=2π/1 jour est la vitesse angulaire. Cette force d�rive du potentiel −m ω2 r2/2 qu’il faut ajouter au potentiel d� � la gravit�, en premi�re approximation mgh o� h+R est la distance au centre de la Terre, (R �tant la distance moyenne de la surface de la Terre � son centre). Si l’on ne tient pas compte de la force centrifuge, h=0 pour avoir un potentiel identiquement nul sur toute la surface, on a une sph�re. Si l’on en tient compte, ce ne sera pas le cas sauf aux poles (o� r=0). Ainsi � l’�quateur, pour avoir un potentiel nul, il faut prendre h= ω2 R2/(2g), soit h=11km environ, il y a donc un bourrelet un peu plus �pais au niveau de l’�quateur. En r�alit�, on observe un aplatissement de l’ordre du double, l’approximation du potentiel de gravitation n’est en effet pas tr�s bonne, parce qu’elle suppose une Terre sph�rique, ce qui n’est pas le cas comme on vient de le voir.

L’existence d’un bourrelet �quatorial et l’inclinaison de l’axe de rotation de la Terre par rapport au plan de l’orbite terrestre entraine l’existence d’un couple exerc� par le Soleil et la Lune sur la Terre. Ce couple s’exerce dans des plans contenant l’axe de rotation de la Terre (qui est axe de sym�trie) et a donc un moment perpendiculaire � l’axe de rotation. Les �quations du mouvement (d�riv�e du moment cin�tique donc de la vitesse angulaire = moment du couple) entrainent alors que la d�riv�e du vecteur rotation est perpendiculaire � l’axe de rotation (et cette d�riv�e est petite), le vecteur rotation conserve donc sa norme au cours du temps mais sa direction bouge lentement perpendiculairement � lui-m�me. Cela entraine en un an une petite rotation du vecteur rotation de la Terre sur un cone autour de la perpendiculaire au plan de l’orbite. On peut faire le calcul en mod�lisant la Terre par un ellipsoide de r�volution (cf. par exemple www-fourier.ujf-grenoble.fr/~faure/enseignement/ puis meca_analytique_L3/index.html#Travaux_dirig%E9s), et on montre que l’action du Soleil et de la Lune s’additionne pour faire faire un tour complet � l’axe de rotation de la Terre autour de la perpendiculaire au plan de l’orbite en un peu moins de 26 000 ans. C’est la pr�cession astronomique.

Le probl�me se complique si on consid�re les autres plan�tes du syst�me solaire, en particulier Jupiter. En effet l’orbite elliptique de la Terre se d�forme au cours des mill�naires, en particulier la position de l’axe des foyers, donc la direction de l’axe aph�lie-p�rih�lie, or c’est l’angle entre cet axe et la projection de l’axe de rotation de la Terre qui importe pour les variations des saisons. C’est pour cette raison que la p�riodicite de la pr�cession “climatique” est un peu plus courte que la pr�cession astronomique.


Ce document a �t� traduit de LATEX par HEVEA

Pourquoi un modèle climatique est une modélisation mathématique du climat ?

Un modèle climatique vise ainsi à représenter le climat et son évolution. Comme le climat est complexe, les modèles climatiques peuvent prendre en compte un nombre fixé de variables et donc se rapprocher plus ou moins de la réalité. Mais la prise en compte d'un grand nombre de phénomènes rallonge le temps de calcul.

Quelle méthode permet de mesurer le réchauffement climatique ?

Découper la terre pour évaluer le réchauffement global Ces observations servent ensuite à faire des projections pour l'avenir grâce à des modélisations. Pour réaliser ces projections, les scientifiques découpent la planète en différentes “mailles” de 100 à 150 km, explique Météo France.

Comment un modèle de simulation climatique est construit ?

Les modèles climatiques sont des systèmes d'équations différentielles basées sur les lois fondamentales de la physique, du mouvement des fluides et de la chimie. Pour «exécuter» un modèle, les scientifiques divisent la planète en une grille tridimensionnelle, appliquent les équations de base et évaluent les résultats.

Comment modéliser l'évolution de la température mondiale ?

L'analyse du système climatique, réalisée à l'aide de modèles numériques, repose sur des mesures et des calculs faisant appel à des lois physiques, chimiques, biologiques connues.