Modèles simples du climat 2

Modèles dynamiques

Pista negra El 6 mayo 2020  - Escrito por  Nils Berglund, avec la participation de Marc Monticelli pour les simulations Ver los comentarios

Dans cette deuxième partie d’une série de trois articles sur des modèles simples du climat, nous allons introduire une version dynamique du modèle à deux boîtes de la circulation thermohaline atlantique. Cela nécessitera l’introduction de mathématiques un peu plus sophistiquées, à savoir d’équations différentielles au lieu d’équations algébriques. En retour, nous obtiendrons des informations sur la stabilité des équilibres déterminés dans la première partie. De plus, nous pourrons étudier la réponse du système à des changements périodiques dans l’apport d’énergie solaire, et à d’autres perturbations.

Dans la première partie de cette série d’articles, nous avons introduit une version d’un modèle à deux boîtes de la circulation thermohaline atlantique nord (c’est-à-dire du Gulf Stream et de son extension, la dérive nord-atlantique). La version initiale est due à Henry Stommel [St], et la version modifiée à Paola Cessi [Ce].

Le modèle est décrit dans la figure du logo de l’article. L’Atlantique nord est représenté par deux boîtes de même volume, l’une plus au sud, à température $T_1$ et salinité $S_1$, l’autre plus au nord, à température $T_2$ et salinité $S_2$. En utilisant la conservation de la quantité totale de sel et de la chaleur, nous avons déterminé des relations satisfaites par les différences $\Delta T = T_1-T_2$ et $\Delta S = S_1 - S_2$ lorsque le système est à l’équilibre. Ceci nous a permis de montrer que pour certaines valeurs du flux $F$ d’eau douce transmise via l’atmosphère par les dépressions, il peut exister plusieurs états d’équilibre, qui se distinguent par l’intensité de la circulation thermohaline.

Cette étude ne permet toutefois pas de décrire ce qui se passe si l’on n’est pas dès le départ dans un état d’équilibre. Afin de s’affranchir de cette limitation, et d’obtenir une version dynamique de ce modèle, nous allons écrire deux équations donnant la variation temporelle, c’est-à-dire la dérivée par rapport au temps, des différences de température et de salinité. La variation de la température de chaque boîte est proportionnelle à la différence de la chaleur reçue et cédée par cette boîte. La variation de la salinité est déduite de manière similaire d’un bilan entre l’eau reçue et le sel cédé via le transfert de masse.

Rappelons que le terme d’échange de masse d’eau $Q(\Delta\rho)$ est donné par l’expression
\[ Q(\Delta\rho) = \frac1{t_{\text{d}}} + C \Delta\rho^2 \]
où $\Delta\rho$ est la différence de densité entre les deux boîtes, qui s’écrit
\[ \Delta\rho = \alpha_S \Delta S - \alpha_T \Delta T \]
Ici $\alpha_S$ et $\alpha_T$ sont des coefficients positifs liés aux propriétés thermodynamiques de l’eau, et $t_{\text{d}}$ est un temps de diffusion thermique, de l’ordre de 200 ans. Celui-ci est beaucoup plus grand que le temps de relaxation $t_{\text r}$ de la différence de température vers sa valeur d’équilibre $\theta$, qui est de l’ordre de 25 jours.

Prenant la différence des équations de bilan de chaleur et de quantité de sel pour chaque boîte, on obtient le système d’équations différentielles
\[ \begin{cases} \displaystyle \frac{\text{d}}{\text{d}t}\Delta T &=& - \displaystyle \frac1{t_{\text r}} (\Delta T - \theta) - Q(\Delta\rho)\Delta T \\ \displaystyle \frac{\text{d}}{\text{d}t}\Delta S &=& F - Q(\Delta\rho) \Delta S \end{cases} \]
Remarquons que les membres de droite des deux égalités s’annulent dans les états d’équilibre déterminés dans la première partie. Dans ce cas, les dérivées de $\Delta T$ et de $\Delta S$ sont nulles, donc ces différences restent bien constantes.

Variables sans dimension

Comme nous l’avons fait dans la première partie, nous allons commencer par introduire des quantités sans dimension afin de simplifier ce système d’équations différentielles. Soit donc $x=\Delta T/\theta$ une température sans dimension, et $y=\alpha_S\Delta S/(\alpha_T\theta)$ une salinité sans dimension. De plus, il s’avère utile de prendre le temps de diffusion $t_{\text{d}}$ comme unité de temps.
Notant $\dot x$ et $\dot y$ les dérivées par rapport à ce nouveau temps, on aboutit à
\[ \begin{cases} \dot x &=& -\displaystyle \frac{1}{\varepsilon}(x-1) - [1+a(y-x)^2] x \\ \dot y &=& p - [1+a(y-x)^2] y \end{cases} \]
Comme dans la première partie, nous avons introduit trois paramètres sans dimension :

  • $\varepsilon = t_{\text{r}}/t_{\text{d}}$ est le rapport entre le temps de relaxation $t_{\text r}$ de la différence de température $\Delta T$ et le temps de diffusion $t_{\text d}$, qui est de l’ordre de $3\times10^{-4}$.
  • $p = \alpha_s t_{\text{d}}F/(\alpha_T\theta)$ mesure le flux d’eau douce via l’atmosphère.
  • $a = C\alpha_T^2\theta^2 t_{\text{d}}$ mesure la sensibilité du taux de transfert de masse à la différence de densité. Sa valeur utilisée dans [Ce] est $a = 7,\!5$.

L’applet java ci-dessous permet de simuler les solutions de ce système d’équations différentielles. Vous pouvez choisir les valeurs initiales de $x$ et $y$ en cliquant en un point du graphique. Vous pouvez également choisir tout un ensemble de valeurs initiales, en bougeant la souris tout en maintenant pressé son bouton. La simulation démarre alors une fois que vous aurez relâché le bouton. Les paramètres peuvent être modifiés à l’aide des régleurs. Par exemple, augmenter la valeur de $\varepsilon$ accélère la dynamique.

En observant le comportement des solutions au bout d’un temps assez long, on pourra notamment vérifier que les points d’équilibre correspondent bien à ceux trouvés dans la première partie, dans la limité $\varepsilon \to 0$. En particulier, nous y avons vu qu’il existe trois points d’équilibre si et seulement si $a>3$ et $p$ appartient à un intervalle dépendant de $a$, que l’on peut lire sur la figure suivante :

Par exemple, pour $a=7,\!5$, cet intervalle est donné approximativement par $]0,\!96;1,\!48[$. Que se passe-t-il lorsque $p$ s’approche d’un bord de cet intervalle ? [1]

Réduction à une dimension

Dans la suite, nous allons étudier une version unidimensionnelle du système ci-dessus, obtenue en supposant le paramètre sans dimension $x$ reflétant l’écart de température est constant et égal à $1$. Comme on le constate sur les simulations, cette approximation est justifiée, lorsque $\varepsilon$ est petit, par le fait que la variable $x$ tend très rapidement vers $1$, puis en reste très proche [2].

Justification mathématique de l’approximation

Un premier argument est basé sur la fonction
\[ U(x) = \frac12 (x-1)^2 \]
qu’on appelle une fonction de Lyapunov. Cette fonction est nulle si $x=1$, et strictement positive autrement. Calculons alors la dérivée de $U(x(t))$ lorsque $x(t)$ est une solution du système d’équations différentielles. Elle vaut
\[ \begin{array}{rcl} \displaystyle \frac{\text{d}}{\text{d}t} U(x(t)) &=& (x(t)-1)\dot x(t) \\ &=& -\displaystyle\frac{1}{\varepsilon}(x-1)^2 - [1+a(y-x)^2] x (x-1) \end{array} \]
Si $x\not\in]0,1[$, alors $x(x-1)$ est positif, et l’on a
\[ \frac{\text{d}}{\text{d}t} U \le -\frac{2}{\varepsilon}U \]
Ceci montre que $U$ décroît exponentiellement vite. Pour $x\in]0,1[$, l’inégalité triangulaire fournit les majorations
\[ \begin{array}{rcl} \displaystyle \frac{\text{d}}{\text{d}t} U(x(t)) &\le& \displaystyle -\frac{2}{\varepsilon}U + \bigl| [1+a(y-x)^2] x (x-1) \bigr| \\ &\le&-\displaystyle \frac{2}{\varepsilon}U + [1+a(|y|+|x|)^2] |x| |x-1| \end{array} \]
Comme $|x-1| = \sqrt{2U}$ et donc $|x| \le 1 + \sqrt{2U}$, on obtient
\[ \frac{\text{d}}{\text{d}t} U(x(t)) \le -\frac{2}{\varepsilon}U + [1+a(|y|+1+\sqrt{2U})^2] \sqrt{2U} (1+\sqrt{2U}) \]
Tant que $U$ est sensiblement plus grand que $\varepsilon^2$, sa dérivée est négative, et donc $U$ décroît également. Par conséquent, quelle que soit la valeur initiale de $x$, $U(x(t))$ va finir par atteindre des valeurs d’ordre $\varepsilon^2$, ce qui signifie que $x(t)-1$ est d’ordre $\varepsilon$.

Un second argument, donnant des informations plus précises, est basé sur la notion de courbe invariante. Supposons qu’il existe une courbe d’équation $x = 1 + \varepsilon h(y)$ telle que toute solution du système démarrant sur cette courbe y reste indéfiniment. Cela revient à exiger que
\[ \varepsilon h'(y) = \frac{\text{d}x}{\text{d}y} = \frac{\dot{x}}{\dot{y}} \]
En effet, cette condition exprime le fait que le champ de vecteurs est partout tangent à la courbe. En remplaçant $\dot{x}$ et $\dot{y}$ par leurs expressions en fonction de $x$ et $y$, avec $x$ remplacé par $1 + \varepsilon h(y)$, on obtient une équation différentielle pour $h(y)$. On peut alors montrer que cette équation admet une solution passant par tous les points d’équilibre. À l’aide d’une fonction de Lyapunov, on montre ensuite que la courbe invariante attire également toute solution de l’équation.

L’équation différentielle de dimension $1$, obtenue dans la limite $\varepsilon\to0$, a la forme
\[ \dot y = p - [1+a(y-1)^2] y \]
Notons qu’elle équivaut à $\dot y = p - f(y)$, où la fonction $f(y) = ay^3 - 2ay^2 + (a+1)y$ a été étudiée dans la première partie. Il est particulièrement utile d’écrire cette équation sous la forme [3]
\[ \dot y = - V'(y) \]
où $V$ est une fonction potentiel donnée par
\[ V(y) = \frac14 a y^4 - \frac23 a y^3 + \frac12(a+1) y^2 - py \]
L’équation $\dot y = - V'(y)$ implique que $V(y)$ tend à diminuer au cours du temps (des précisions sont données dans le bloc dépliant un peu plus bas).

Ci-dessous, on a représenté le potentiel pour $a=7,\!5$, et pour $p=0,\!9$ à gauche, $p=1,\!2$ au centre, et $p=1,\!5$ à droite. Remarquez que seule la valeur $1,\!2$, qui se trouve dans l’intervalle où plusieurs équilibres coexistent, donne lieu à un potentiel à deux puits. Les autres valeurs de $p$ sont en-dehors de cet intervalle, et conduisent à un potentiel avec un seul minimum.

PNG - 462.9 KB
PNG - 462.9 KB
PNG - 462.9 KB

C’est naturel, puisque les points stationnaires de $V$ (c’est-à-dire ses points critiques, et en particulier les extrema locaux) sont exactement les points où $V'$ s’annule. Or, la condition $V'(y)=0$ équivaut à l’équation $f(y) = p$, dont les solutions sont les points d’équilibre.

L’intérêt de cette approche est qu’elle permet de visualiser la dynamique de l’équation différentielle $\dot y = -V'(y)$ comme le mouvement d’une bille évoluant dans ce potentiel, et soumise à une force de frottement visqueux. La bille reste immobile si elle se trouve en un minimum ou en un maximum local de $V$, mais sinon elle tend à rouler le long de la pente du potentiel, pour s’approcher d’un minimum local. Il suit de cette observation que, comme nous l’avons annoncé dans la première partie, s’il y a trois états d’équilibre, celui du milieu est instable, puisqu’il correspond à un maximum local du potentiel. Dès que la bille est déplacée légèrement de ce maximum, elle tombe dans l’un des puits de potentiel.

Détails sur l’analogie de la bille

Observons d’abord que la dérivée du potentiel le long d’une trajectoire vaut
\[ \frac{\text{d}}{\text{d}t} V(y(t)) = V'(y(t)) \dot y(t) = -V'(y(t))^2 \le 0 \]
Ainsi, $V(y(t))$ décroît strictement en tout point où la dérivée de $V$ est non nulle, et est constante autrement. En d’autres termes, $V$ est à nouveau une fonction de Lyapunov. Ceci confirme que $V(y(t))$ ne peut jamais augmenter, comme pour une bille amortie roulant dans le puits de potentiel.

Pour mieux comprendre le lien avec le mouvement d’une bille, appliquons la deuxième loi de Newton à son mouvement. Celle-ci affirme que
\[ m \ddot{y} = - \gamma \dot{y} - V'(y) \]
où $m$ est la masse de la bille, et $\gamma$ est un coefficient de frottement visqueux. Le terme $-\gamma y'$ décrit une force de frottement visqueux, proportionnelle à la vitesse, alors que $-V'(y)$ est la force dérivant du potentiel. Il s’agit d’une équation différentielle du second ordre, dont les solutions peuvent présenter des oscillations. Toutefois, si $\gamma$ est suffisamment grand par rapport à la masse $m$, on peut vérifier que les solutions sont proches de celles de l’équation du premier ordre
\[ \gamma \dot{y} = - V'(y) \]
qui est équivalente à $\dot{y} = -V'(y)$ via un changement d’unité de temps.
Plus précisément, on écrit d’abord l’équation de second ordre de Newton comme un système de deux équations du premier ordre, puis on montre, un peu comme dans le premier bloc dépliant, qu’il existe une courbe invariante sur laquelle le mouvement est décrit par l’équation du premier ordre.

Flux d’eau douce dépendant du temps

Maintenant que nous avons bien compris la dynamique de l’équation du premier ordre, nous pouvons modifier le modèle, afin de le rendre plus réaliste. Jusqu’ici, nous avons supposé le flux $p$ d’eau douce via l’atmosphère constant. Or, comme nous l’avons signalé dans la première partie, le climat terrestre est soumis à différentes influences qui peuvent varier au cours du temps. Parmi elles, on trouve des influences d’origine astronomique, liées à des changements dans l’orbite terrestre, en particulier de son excentricité. Celle-ci influence la quantité de chaleur reçue de la part du soleil, qui modifie à son tour la quantité d’eau évaporée et l’intensité des précipitations.

Il est donc pertinent de considérer l’équation différentielle dépendante du temps
\[ \dot y = p(t) - [1+a(y-1)^2] y \]
Pour rester simples, supposons que le flux d’eau douce est périodique en temps, de la forme
\[ p(t) = p_0 + A \cos(\omega t) \]
Rappelons que nous utilisons le temps de diffusion $\text{t}_d$, qui est de l’ordre de 200 ans, comme unité de temps. Nous avons vu qu’un des mécanismes proposés pour expliquer l’alternance entre périodes glaciaires et interglaciaires est une variation de l’insolation sur une échelle d’environ 100 000 ans. Cela implique que la pulsation (ou «fréquence angulaire») $\omega$ est de l’ordre de $2\times 10^{-3}$.

L’applet java suivant permet de simuler les solutions de cette équation différentielle, pour différentes valeurs des paramètres $p_0$, $A$ et $\omega$.

La figure ci-dessous montre deux exemples de solutions périodiques (courbes bleues), dans le plan $(p,y)$, obtenues pour un forçage symétrique ($p_0$ est la moyenne des deux valeurs limites $p_\pm$ de $p$ pour lesquelles il existe plusieurs équilibres), des amplitudes $A$ différentes, et une valeur de $\omega$ de l’ordre de $0,\!1$. La courbe noire en forme de S indique les points d’équilibre statiques, les points stables en trait plein, et les points instables en pointillé.

PNG - 583.5 KB
Solutions de l’équation dépendante du temps $\dot{y} = p(t)) - [1+7,\!5(y-1)^2]$ avec $p(t) = p_0 + A \cos(\omega t)$, représentées dans le plan $(p,y)$. Les valeurs des paramètres sont $p_0 = 1,\!22$, $\omega=0,\!1$, et $A = 0,\!2$ à gauche, $A = 0,\!75$ à droite.

Comme on le voit sur la figure de gauche, si $A$ est trop petit, alors $p(t)$ reste toujours dans le domaine où il existe trois équilibres, et la solution périodique reste voisine du même équilibre statique. En revanche, sur la figure de droite on voit que si l’amplitude $A$ est assez grande, la solution périodique a un comportement d’hystérésis (dynamique) : deux fois par période, la différence de salinité fait une transition assez abrupte entre un régime de fort taux d’échange de masse et un régime de faible taux d’échange. Cela signifie que l’intensité de la circulation thermohaline change périodiquement entre un régime très actif et un régime peu actif. Dans l’image de la particule dans un potentiel, le cycle se comprend ainsi :

On notera que les transitions sont tout de même moins abruptes que dans le cas statique. En effet, les solutions de l’équation différentielle dépendante du temps présentent une certaine inertie, et ne réagissent donc pas immédiatement à la disparition de l’état d’équilibre qu’elles suivent. Pontriaguine a montré dans [Po] que ce retard est d’ordre $\omega^{2/3}$ pour $\omega$ petit [4].

Flux d’eau douce dépendant de l’état

Outre les variations de flux causées par des influences externes, il peut également exister des variations causées par l’état même du système. En effet, la différence de salinité détermine l’intensité de la circulation thermohaline, qui peut à son tour influencer le flux d’eau douce. On aboutit alors à un système de la forme
\[ \begin{cases} \dot y &=& p - [1+a(y-1)^2] y \\ \dot p &=& F(y,p) \end{cases} \]
où $F(y,p)$ est une fonction qui modélise l’effet de la différence de salinité et du flux sur ses variations temporelles.

Pour que ce système ait des solutions périodiques, il faut que $F$ soit négative sur la branche d’équilibre supérieure, correspondant à une forte circulation thermohaline, et positive sur la branche inférieure, correspondant à une circulation faible. On peut par exemple imaginer que comme une circulation faible augmente la différence de température entre les deux boîtes, elle augmente également le flux d’eau douce - c’est ce qu’on appelle un effet de rétroaction négative. Mais un bon modèle de la fonction $F$ nécessite une meilleure compréhension des différents effets impliqués. Si $F$ a effectivement les bonnes propriétés, on aboutit à nouveau à des solutions sous forme de cycles d’hystérésis, semblables à celles décrites ci-dessus.

Remarquons pour terminer que les phénomènes d’hystérésis sont fréquents dans le système climatique. Le schéma ci-dessus en donne un autre exemple : l’épaisseur de la calotte glaciaire du Groenland résulte de l’équilibre entre l’accumulation de neige dans les régions élevées et la fonte des glaces dans les régions plus basses. Si la limite pluie-neige s’élève en raison d’un réchauffement climatique, la fonte l’emporte sur l’accumulation de neige, et la calotte glaciaire diminue d’épaisseur. Même si la limite pluie-neige revient ensuite à son niveau initial, la surface de la glace étant plus basse qu’avant, la calotte glaciaire continue de fondre. Pour que la tendance s’inverse, il faut que la limite pluie-neige descende suffisamment pour qu’il y ait à nouveau prédominance de l’accumulation de glace.

Conclusion

Nous avons vu plusieurs manières d’enrichir le modèle à deux boîtes statiques introduit dans la première partie, d’abord en le rendant dynamique, puis en permettant des variations du flux d’eau douce. Cela nous a permis d’observer des comportements plus riches, comme l’évolution vers un équilibre stable à partir d’un état initial qui n’est pas un équilibre, ou encore des solutions variant périodiquement au cours du temps, lorsque le flux d’eau douce varie.

Un inconvénient majeur du modèle reste que ses solutions sont asymptotiquement périodiques. Or, nous avons vu dans la première partie que les événements de Daansgard-Oeschger, ou réchauffements subits durant la dernière ère glaciaire, ont un comportement beaucoup plus irrégulier. Dans la dernière partie de cette série d’articles, nous verrons comment on peut obtenir de tels comportements en ajoutant un terme aléatoire à notre équation différentielle.

Bibliographie

[Ce] Cessi, P. (1994), A Simple Box Model of Stochastically Forced Thermohaline Flow. J. Phys. Oceanogr. 24: 1911-1920.

[Po] Pontryagin, L. S. (1957), Asymptotic behavior of solutions of systems of differential equations with a small parameter in the derivatives of highest order. J. Izv. Akad. Nauk SSSR Ser. Mat. 21 (5): 605-626.

[St] Stommel, H. (1961). Thermohaline Convection with Two Stable Regimes of Flow. Tellus. 13 (2): 224-239.

Post-scriptum :

Merci aux relecteurs Lalanne, Jean-Romain et Himynameisarno pour leurs suggestions, qui ont permis d’améliorer la lisibilité de cet article.

Article édité par Jérôme Buzzi

Notas

[1Les lignes violettes délimitant les régions à un et trois équilibres dans le plan $(p,a)$ s’appellent des lignes de bifurcation. En franchissant l’une de ces lignes, un point d’équilibre stable et un point d’équilibre instable «entrent en collision» et «s’annihilent».

[2Cette question n’est pas anodine, car même si deux équations différentielles sont proches, cela ne signifie pas automatiquement que leurs solutions resteront proches pour tous les temps.

[3Rappelons qu’ici $\dot y$ désigne la dérivée de $y$ par rapport au temps, alors que $V'(y)$ désigne la dérivée de $V$ par rapport à la variable $y$.

[4Le résultat de Pontriaguine nécessite une étude de l’équation différentielle dépendante du temps au voisinage des «points de rebroussement» de la courbe en forme de S décrivant les états d’équilibre, c’est-à-dire des points où les états stables (en traits pleins) et instables (en pointillé) se rejoignent. Cela fait intervenir des mathématiques assez sophistiquées. Mais l’idée de base est qu’il faut zoomer sur le point de rebroussement, en changeant les échelles d’espace et de temps. Cela permet de se ramener à l’étude d’équations différentielles bien connues, appelées équations de Riccati, qui dans ce cas particulier peuvent être transformées en une équation d’Airy.

Comparte este artículo

Para citar este artículo:

Nils Berglund, avec la participation de Marc Monticelli pour les simulations — «Modèles simples du climat 2» — Images des Mathématiques, CNRS, 2020

Comentario sobre el artículo

Dejar un comentario

Foro sólo para inscritos

Para participar en este foro, debe registrarte previamente. Gracias por indicar a continuación el identificador personal que se le ha suministrado. Si no está inscrito/a, debe inscribirse.

Conexióninscribirse¿contraseña olvidada?

Registros

Este artículo es parte del registro «Modèles simples de climat» consulte el registro
La traducción del sitio del francés al castellano se realiza gracias al apoyo de diversas instituciones de matemáticas de América Latina.