Bifurcations dans un modèle climatique simple
Les questions

Marie-Cécile Selosse
(last modiﬁcation date: October 10, 2017)

Version pdf de ce document
Version sans bandeaux

1 Un modèle climatique simple (Ghil et Le Treut, 1981)
2 Écriture du modèle sous Scilab
3 Simulations et interprétations

1 Un modèle climatique simple (Ghil et Le Treut, 1981)
1.1 Variations de la température
1.2 Variations de l’extension méridionale de la calotte
2 Écriture du modèle sous Scilab
2.1 Données numériques
2.2 Fonctions physiques auxiliaires
2.3 Dynamique du système
3 Simulations et interprétations
3.1 Comportement global du système (pour une valeur ﬁxée du paramètre de bifurcation μ)
3.2 Sensibilité du système au paramètre de bifurcation μ

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 :

dT- CT dt = Ri − Ro, = Q {1 − [γα (L) + (1 − γ)α (T )]} − κ (T − T ). (1) land ocean κ

Dans cette équation, les paramètres sont :

capacité caloriﬁque 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éﬁnitions 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 :

αland(L ) = α0 + α1L, 0 ≤ L ≤ Lmax,

(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 aﬃne par morceaux de la température :

( || αmax, T ≤ T α,lower, { --αmin-−-αmax---- αocean = T − T (T − Tα,lower) + αmax, Tα,lower < T ≤ T α,upper, ||( α,upper α,lower αmin, Tα,upper < T.

(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 diﬀérence entre l’accumulation de neige et sa fonte.

On fait l’hypothèse que le proﬁl 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 :

1 1 HL (x ) = λ2(L − |x|)2.

(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 ﬁgure 1.

Figure 1: Modélisation de la calotte glaciaire

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 :

dV- ′ dt = aLa − a La′.

(5)

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

dV dV dL --- = --- ---, dt dL dt

(6)

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

∘ ------ dL Lmax --- = μ -----[(1 + 𝜀)La − L], dt L

(7)

où μ = a′ 3 1 _____ -------
√ λLmax 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éﬁnition de la zone d’accumulation

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

hiso(x,T, L) = h0(T ) + s(x + L ),

(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

h0(T) = β(T − T00),

(9)

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

-2s(αmax-−-α0-)- β = T − T α . 00 α,lower 1

(10)

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

[( ) 12 ( ) ] La = 1- 2s2L + sh0 + 1- − s2L + sh0 + 1- . s2 4 2

(11)

Déﬁnition du rapport accumulation/ablation

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

( ||{ 𝜀min, T ≤ T 𝜀,lower, 𝜀(T ) = ---𝜀max −-𝜀min--(T − T𝜀,lower) + 𝜀min, T𝜀,lower < T ≤ T𝜀,upper, || T𝜀,upper − T𝜀,lower ( 𝜀max, T𝜀,upper < T.

(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 eﬀet 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 suﬃsamment vigoureux pour amener les précipitations au-dessus de la calotte.

2 Écriture du modèle sous Scilab

En déﬁnitive, le modèle s’écrit

( || dT- -1- |{ dt = CT (Q {1 − [γ αland(L ) + (1 − γ)αocean(T)]} − κ(T − T κ)) ∘ ------ ||| dL Lmax ( --- = μ -----[(1 + 𝜀(T))La − L ] dt L

(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 dT-
dt et dL
dt . Pour cela, nous avons besoin de coder les diﬀérentes fonctions auxiliaires pour obtenir L_a, a_ocean, 𝜖.

2.1 Données numériques

  s=0.3E-3;

  T00=283;
  Talower=217;
  Taupper=283;

  a0=0.25;
  a1=4.1E-7;
  amax=0.85;
  amin=0.25;

  Tepslower=273;
  Tepsupper=283;

  epsmin=0.1;
  epsmax=0.5;

  Lmax=1.44E6;
  Q=362.2;
  CT=1;
  gama=0.3;//la variable gamma est prot\'eg\'ee
  kappa=1.74;
  Tkappa=154;

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

exec("nom_du_fichier_des_donnees.sce","c")

2.2 Fonctions physiques auxiliaires

On peut placer toutes les fonctions qui suivent dans le même ﬁchier nom_du_fichier.sci. Il suﬃra alors, pour les charger sous Scilab, de taper les commandes suivantes

exec("nom_du_fichier_des_donnees.sce","c")
getf("nom_du_fichier.sci","c")

Abscisse limite L_a de la zone neigeuse

  function x=La(T,L)
    beta=2*s*(amax-a0)/(T00-Talower)/a1;
    h0=beta*(T-T00);
    x=(sqrt(2*s^2*L+s*h0+1/4)-(s^2*L+s*h0+1/2))/s^2
  endfunction

Albédo de l’océan a_ocean

  function x=aocean(T)
    if T <= Talower then
      x=amax
    elseif T <= Taupper then
      x=amax+(T-Talower)*(amin-amax)/(Taupper-Talower)
    else
      x=amin
    end
  endfunction

Rapport accumulation / ablation 𝜖

  function eps=epsilon(T)
    if T <= Tepslower then
      eps=epsmin
    elseif T <= Tepsupper then
      eps=epsmin+(T-Tepslower)*(epsmax-epsmin)/.(Tepsupper-Tepslower)
    else
      eps=epsmax
    end
  endfunction

2.3 Dynamique du système

  function ydot=Ghil(t,y)
    //y(1)=T,temp\'erature
    //y(2)=L,extension
    ydot(1)=Q/CT*(1-(gama*(a0+a1*y(2))+(1-gama)*aocean(y(1))))-kappa*(y(1)-Tkappa);
    ydot(2)=mu*sqrt(Lmax/y(2))*((1+epsilon(y(1)))*La(y(1),y(2))-y(2));
  endfunction

3 Simulations et interprétations

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

Visualisation des points d’équilibre

Aﬁn 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 $dL- dt$ et $dT- dt$ s’annulent. Pour cela, nous déﬁnissons 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.

mu=0.5+(1.8-0.5)*rand

Écriture des composantes du champ de vecteurs

function z=Ghil1(T,L) z=Q/CT*(1-(gama*(a0+a1*L)+(1-gama)*aocean(T)))-kappa*(T-Tkappa);endfunction
function z=Ghil2(T,L) z=3*mu/2*sqrt(Lmax/L)*((1+epsilon(T))*La(T,L)-L);endfunction

Tracé des isoclines

  T=250:300;
  L=5E5:1E5:15E5;
  fcontour2d(T,L,Ghil1,[0,0],[9,9])
  fcontour2d(T,L,Ghil2,[0,0],[12,12],"000")

Question 1 Commenter vos observations sur les points d’équilibre.

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 :

portrait(Ghil)

Apparaissent ensuite plusieurs boîtes de dialogue. Aﬁn 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 :

portrait(Ghil,"default",[250,5E5,300,15E5],[10,0.1])

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.

Question 2 Commenter la forme du champ de vecteurs.

Retrouver les points d’équilibre. Identiﬁer leur nature (stabilité) ?

Examen de quelques trajectoires

Pour visionner une trajectoire correspondant à une condition initiale donnée, il suﬃt 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 diﬀérentielles ou récurrentes) dans la boîte de dialogue autant de fois que l’on veut.

Question 3 Tracer des trajectoires en cliquant un point initial dans la fenêtre ScilabGraphic.

Il arrive que l’intégration numérique avec la fonction Scilab ode échoue. Pourquoi ? (Ceci peut être relié à des propriétés mathématiques de stabilité du système (??), ou à des sorties hors du domaine de validité physique des formules comme (11)).

Observer quelques trajectoires. En déduire la nature des points d’équilibre.

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

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

Aﬁn 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 diﬀérentielles avec la fonction Scilab ode.

y=ode(y0,t0,t,Ghil)

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 diﬀé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ériﬁer que les deux vecteurs sont de même taille. En eﬀet, 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. Aﬁn 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.

  t=0:0.1:100;
  y=ode([278,9E5]',0,t,Ghil);

  //pour ne prendre que les valeurs de t pour lesquelles la temp\'erature a
  //effectivement \'et\'e calcul\'ee~:
  xbasc()
  si=size(y)
  plot(t(1):0.1:t(si(2)),y(1,:))
  //attention: le pas de temps dans plot doit \^etre le m\^eme que dans t

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 ≤ μ.

Question 4 Pour chaque valeur de μ choisie, observer :

le portrait de phase ;
quelques trajectoires, dont vous décrirez le comportement ;
la variation de la température au cours du temps pour une condition initiale que vous choisirez.

En étudiant cette dernière courbe, déterminer l’amplitude et la période des oscillations, ainsi que la convergence ou la divergence de la trajectoire choisie.

Y a-t-il des attracteurs ? Si oui, de quelles formes ?

Quelles sont les valeurs critiques de μ ? Quels types de bifurcation observe-t-on ?

Pour mieux observer les diﬀé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.

Question 5 Dans chaque cas, commenter vos observations.

Comment se comporte le point d’équilibre central ?

Décrire le comportement des trajectoires.

Bifurcations dans un modèle climatique simple Les questions

Contents