Contents

1 Un modèle climatique simple (Ghil et Le Treut, 1981)

Dans (Ghil M. & Tavantzis J. Global Hopf Bifurcation in a Simple Climate Model. SIAM J. Appl. Math., 43(5):1019–1041, (October 1983)), Ghil et Tavantzis s’inspirent d’un modèle de Ghil et Le Treut (1981), et représentent simplement l’évolution de la calotte glaciaire à l’aide de deux variables :

• la température moyenne de l’atmosphère T ;
• l’extension méridionale de la calotte polaire L.

1.1 Variations de la température

Les variations de la température T dépendent :

• de l’énergie solaire incidente R_i reçue par le système ;
• de l’énergie R_o réémise directement (à cause de l’albédo de la surface, et sous forme de rayonnement infrarouge).

Bilan d’énergie

Le bilan d’énergie correspondant s’écrit alors :

Dans cette équation, les paramètres sont :

• capacité calorifique du système climatique ; on normalise le temps t en mettant cette valeur à 1 ;
• constante solaire, que l’on peut faire varier avec les périodes de Milankovitch ;
• albédo terrestre ;
• albédo de l’océan ;
• paramètre de proportion entre terre et océan ;
• rayonnement infra-rouge de la planète.

Définitions des albédos

L’albédo terrestre α_land dépend essentiellement de celui de la calotte de glace ; on supposera donc qu’il suit une loi donnée par :

(2)

où L_max représente l’extension maximale de la glace (continent recouvert par la calotte).

L’albédo de l’océan α_ocean est donné par une fonction affine par morceaux de la température :

(3)

La dépendance de α_ocean en la température T tient au fait que, à ces échelles de temps, la glace de mer s’allonge ou se rétracte pratiquement instantanément avec les variations de température, ce qui n’est pas le cas des calottes de glace continentales.

1.2 Variations de l’extension méridionale de la calotte

L’extension de la calotte glaciaire est la différence entre l’accumulation de neige et sa fonte.

On fait l’hypothèse que le profil géométrique H_L(x) de la calotte possède une forme constante, constitué de deux demi-paraboles avec un terme correcteur pour la déformation du sol :

(4)

On supposera dans le modèle que la calotte est limitée au nord par un océan (Arctique) et ne peut s’étendre que vers le sud. Ce modèle de calotte est décrit par la figure 1.

---

---

La calotte possède une zone active, divisée en une partie d’accumulation a de longueur L_a et une partie d’ablation a′ de longueur L_a′, telles que L_a + L_a′ = L.Ces deux parties sont séparées par une droite d’isotherme 0 au nord de laquelle l’accumulation est prépondérante, et qui est dominée par l’ablation au sud.

Conservation de la masse

Le bilan de masse de la calotte s’écrit :

(5)

avec V = ∫ _−L^LH _L(x)dx. Or, en dérivant l’équation de volume, on trouve

(6)

ce qui permet de trouver une équation pour l’évolution de L

(7)

où μ = a′ 3 1 _____ et 𝜀 = a∕a′ est le rapport accumulation / ablation.¹

Le modèle ainsi décrit est donc une déduction ¡¡ naïve ¿¿ des premiers principes de bilan d’énergie et de conservation de la masse.

Définition de la zone d’accumulation

Dans ce modèle, L_a est définie par l’intersection du profil de la calotte avec l’isotherme 0 d’équation :

(8)

où s est la pente de l’isotherme, h₀ est la hauteur de l’isotherme sur la côte arctique (x = −L) paramétrisée par

(9)

et T₀₀ est la température pour laquelle l’isotherme intercepte la côte arctique ; β est défini par :

(10)

Ceci permet d’estimer la zone d’accumulation L_a :

(11)

Définition du rapport accumulation/ablation

Le rapport 𝜀 est paramétrisé par une fonction affine par morceaux de T, pour exprimer la rétroaction du cycle hydrologique sur la calotte de glace :

(12)

La croissance de 𝜀 avec T exprime le fait que ce sont les hivers les plus doux qui possèdent les plus forts taux de précipitations neigeuses. Ainsi, lorsque la température est “élevée”, l’évaporation dans les zones océaniques de basses latitudes est plus forte, l’humidité est advectée vers le nord où elle précipite, en accentuant le bilan de masse de la calotte polaire. On suppose alors que cet effet dépasse l’augmentation de la fonte de la calotte durant des étés chauds. A l’inverse, si la température est trop faible, l’océan est couvert de glace de mer et le cycle hydrologique n’est pas suffisamment vigoureux pour amener les précipitations au-dessus de la calotte.

2 Écriture du modèle sous Scilab

En définitive, le modèle s’écrit

(13)

avec les formules (2), (3), (11) et (12), et les constantes données ci-dessous.

Le modèle (13) peut être codé sous Scilab en écrivant la fonction Ghil(t,y), qui associe au vecteur y = [T,L] le vecteur formé des expressions de et . Pour cela, nous avons besoin de coder les différentes fonctions auxiliaires pour obtenir L_a, a_ocean, 𝜖.

2.1 Données numériques

Pour récupérer les données utiliser le lien : [onnées]données.code/donnees.sce Les différentes constantes du problème peuvent être placées à part dans un fichier
nom_du_fichier_des_donnees.sce. Par la suite, ces valeurs pourront être chargées sous Scilab par la commande

2.2 Fonctions physiques auxiliaires

On peut placer toutes les fonctions qui suivent dans le même fichier nom_du_fichier.sci. Il suffira alors, pour les charger sous Scilab, de taper les commandes suivantes

Abscisse limite L_a de la zone neigeuse

Albédo de l’océan a_ocean

Rapport accumulation / ablation 𝜖

2.3 Dynamique du système

3 Simulations et interprétations

3.1 Comportement global du système (pour une valeur fixée du paramètre de bifurcation μ)

Visualisation des points d’équilibre

Afin de repérer les points d’équilibre du système, nous allons traçer ses isoclines, c’est-à-dire les courbes sur lesquelles les dérivées et s’annulent. Pour cela, nous définissons les fonctions correspondant aux expressions de ces deux dérivées, que nous appelons Ghil_1 et Ghil_2.

Valeur du paramètre de bifurcation μ La fonction Ghil fait appel à un paramètre μ dont la valeur mu doit être précisée avant de lancer les programmes précédents. On prendra pour mu un réel compris entre 0.5 et 1.8.

Écriture des composantes du champ de vecteurs

Tracé des isoclines

Visualisation du portrait de phase global

Nous obtenons le portrait de phase des solutions du système par la fonction Scilab portrait. L’appel le plus simple est :

Apparaissent ensuite plusieurs boîtes de dialogue. Afin d’avoir une vue globale, vous pouvez choisir les options suivantes :

• la température peut être comprise entre 250 K et 300 K ;
• l’extension peut être choisie entre 5.10⁵ m et 15.10⁵ m ;
• choisir de dessiner le champ de vecteurs ;
• choisir 15 points est satisfaisant pour visualiser l’ensemble du comportement ;
• la boîte de dialogue qui suit concerne le pas d’intégration du système ; un pas de 0.1 et un nombre de points de 10 sont satisfaisants.

Vous pouvez aussi entrer directement ces données lors de l’appel de la fonction de la façon suivante :

Il vous sera alors demandé si vous voulez ou non dessiner le champ de vecteurs.

À ce stade, répondre Quit à la question Choose dans la boîte de dialogue Scilab Choose Panel.

Examen de quelques trajectoires

Pour visionner une trajectoire correspondant à une condition initiale donnée, il suffit de choisir l’option New initial point dans la boîte de dialogue Scilab Choose Panel et de pointer sur le portrait de phase à l’endroit de la condition initiale voulue. Le programme trace alors une partie de la trajectoire ; pour en voir une plus grande partie, on peut choisir l’option Continue ode (ode est un solveur de systèmes d’équations différentielles ou récurrentes) dans la boîte de dialogue autant de fois que l’on veut.

Tracé d’une courbe d’évolution de température

On peut tracer des courbes de température grâce aux fonctions ode et plot.

Afin de tracer l’évolution de la température avec le temps, pour une condition initiale donnée, on résoud le sytème d’équations différentielles avec la fonction Scilab ode.

où y₀ est le vecteur (colonne) des conditions initiales, t₀ le temps initial et t le vecteur des temps pour lesquels la solution est calculée. Ici, la première ligne du vecteur solution y contiendra les différentes températures, tandis que la deuxième contiendra les extensions méridionales.

Attention. Lorsque vous utilisez la fonction plot pour tracer la température en fonction du temps, pensez à bien vérifier que les deux vecteurs sont de même taille. En effet, lorsqu’il y a divergence, la calcul de la température s’arrête : le vecteur des températures est donc plus court que celui des pas de temps. Afin de toujours tracer des vecteurs de même taille, même en cas de divergence du solveur ode, il est recommandé de s’inspirer de l’exemple qui suit.

3.2 Sensibilité du système au paramètre de bifurcation μ

Nous allons étudier la sensibilité du système au paramètre de bifurcation μ : μ < 1.6, 1.6 ≤ μ < 1.7, 1.7 ≤ μ.

Pour mieux observer les différences de comportement du système, on peut se placer au niveau de la ¡¡ spirale ¿¿ observée précédemment, en choisissant une échelle de température comprise entre 268 et 290 dans l’appel de la fonction portrait.