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].

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.

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.

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é.

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.