OSCILLATIONS DÉCENNALLES THERMIQUES DANS UN MODÈLE COUPLÉ OCÉAN-ATMOSPHÈRE.

Guillaume MAZE

Thierry HUCK

Laboratoire de Physique des Océans

SUJET DU STAGE:

La variabilité du climat de l'Atlantique Nord est associée à des variations de l'étendue de la glace de mer dans les hautes latitudes. Néanmoins les relations de cause à effet ne sont pas très claires bien que les changements de circulation atmosphériques semblent précéder de quelques mois les déplacements de glace (Deser et Blackmon 1993). La glace de mer joue un rôle important dans les interactions air-mer en réduisant à néant les flux de chaleur entre l'océan et l'atmosphère qui dépassent facilement les 100 W$m^{-2}$ dans ces régions en l'absence de glace. C'est un élément interactif majeur du climat qui joue un rôle essentiel à toutes les échelles de temps, en gardant la mémoire du système d'une année à l'autre.

La dynamique de la glace (influence du vent et des courants marins, rhéologie) est néanmoins suffisament complexe pour qu'elle reste un problème majeur dans la plupart des modèles les plus réalistes.

À l'aide de modèles thermodynamiques simplifiés de glace, plusieurs auteurs ont mis en évidence des oscillations dues aux interactions océan - glace de mer (Zhang et al. 1995, Yang et Neelin 1997) et leur mécanisme parait bien compris (Welander 1977). Les interactions avec l'atmosphère dans ces modèles ne sont toutefois pas très satisfaisantes d'un point de vue de conservation de l'énergie. Dans la plupart des problèmes climatiques, il est nécessaire de satisfaire les bilans radiatifs au sommet de l'atmosphère et un des modèles les plus simples respectant cet équilibre est utilisé depuis longtemps pour étudier la réponse du climat à des changements d'insolation par exemple (``Energie Balance Model''). L'objectif est donc ici de prendre en compte la formation de la glace de mer dans un modèle couplé océan - atmosphère simplifié en équilibre énergétique et d'étudier son influence sur la stabilité et la variabilité du système.


Table des matières

Introduction

La nécessité de coupler un modèle thermodynamique de glace à un modèle de circulation océanique dont au moins une zone recouvre les latitudes polaires est évidente. La présence ou non de glace à la surface de l'océan joue un rôle majeur dans les quantités de chaleur échangées à l'interface océan-atmosphère (la glace est quasiment un isolant thermique) et donc dans la détermination de la température océanique de surface (dont le rôle de conditions aux limites dans énormément de modèles est primordiale). Cette influence se retrouve également dans la salinité car la glace de mer joue un rôle de source (formation de glace) ou de puit (fonte de glace) de sel pour l'océan de surface. Cette influence de la glace de mer sur la température et la salinité et donc surtout sur la densité des eaux superficielles constitue l'un des moteurs de la circulation océanique. Tous les phénomènes s'y rapportant ont donc un rôle à jouer dans la compréhension globale du système.

Or les modèles (2D, 3D ou en boîtes) de circulation océanique, qu'ils soient couplés ou non à l'atmosphère, présentent une complexité de plus en plus croissante. Ce degré de perfectionnement, qui tend à reproduire de plus en plus fidèlement les phénomènes naturels, engendre parallèlement une certaine difficulté d'identification précise des processus physiques mis en cause. C'est pourquoi reconnaître des phénomènes simples établis après diverses hypothèses simplificatrices (mais souvent reproduits en laboratoire), au milieu de modèles toujours plus complets (les phénomènes physico-chimiques commencent même à y faire leur entrée) n'est pas très aisé. L'intérêt d'une telle démarche peut être de déterminer si des phénomènes qui apparaissent dans des hypothèses physiques relativement simples, se maintiennent dans un milieu plus réaliste et dans quelles mesures les divers paramètres naturels l'influencent.

Divers régimes d'oscillation de la circulation thermohaline (THC) et de l'étendue de la glace de mer, recouvrant une gamme de période très large (de quelques dizaines d'années à plusieurs millénaires) ont été découverts et très bien étudiés. Toutefois ce n'est pas le cas pour les oscillations décennales thermiques d'un fluide réchauffé par le bas et refroidi par le haut découvertes par Welander et exposées en 1977 dans la revue Dynamics of Atmospheres and Oceans.

Welander montre qu'un modèle simple thermodynamique glace-océan-atmosphère posséde trois états d'équilibres:

Dans ce modèle binaire (présence ou non de glace sur la colonne de fluide) le regime oscillatoire peut s'expliquer comme suit: l'atmosphère refroidi l'océan jusqu'à sa température de congélation, il y a alors formation de glace. La faible conductivité thermique de celle-ci isole bientôt complétement l'océan de l'atmosphère et arrive un stade où l'océan réchauffe plus la glace que l'atmosphère ne la refroidit et permet son extension (dans ce cas uniquement en épaisseur et pas en surface). La glace commence alors à fondre et quand elle a complètement disparue, l'atmosphère reprend à refroidir l'océan jusqu'à sa température de congélation et la reprise de la formation de glace.

Les états ice-covered et ice-free ne nous intéressent pas ici. Nous allons par contre tenter d'analyser dans quelles mesures le régime d'oscillations auto-entretenues (d'une période de l'ordre de 5 à 10 ans) peut être retrouvé et maintenu dans un modéle thermodynamique plus complet (Bjornsson 2000). Nous l'appliquerons ensuite à un modèle en boîte simple océan-atmosphère et équatorial-pôlaire pour étudier les éventuelles oscillations.

Le modèle thermodynamique

Étape 1: Welander 1977 (EBM0)

Équations

On considère une colonne verticale de surface unité composée de haut en bas: d'eau de mer sur une hauteur $H_e=400\,m$, éventuellement de glace d'épaisseur h ($h< <H_e$ ) et d'air sur une hauteur $H_a=8000\,m$.

Les hypothèses formulées par Welander sont:

Il est clair que de considérer la température atmosphérique comme fixe constitue un cas particulier dont il est indispensable de sortir pour émettre une quelconque hypothèse sur la réalité. Et même lorsque $T_a$ ne sera pas constante, les modèles trop simples utilisés dans cette étude ne permettront toujours pas de donner des informations plausibles sur le monde réel. L'objectif ici est de discuter de la stabilité et de la variabilité de la THC à travers des modèles simples en équilibre thermodynamique où sont facilement identifiables l'influence de termes dont la signification physique est claire.

Les flux de chaleur océan-atmosphère $Q_{oa}$ (W$m^{-2}$), glace-atmosphère $Q_{ia}$ (W$m^{-2}$) et océan-glace $Q_{oi}$ (W$m^{-2}$) ont donc pour expression (EBM0 considère que l'océan en contact direct avec la glace est à la température de fusion $T_{fr}$):


(2.1)

avec $T_{fr}=-1.9\,^oC$ et $E_{oa}$, $E_{oi}$ et $E_{ia}$ les coefficients d'échange thermodynamique respectivement aux interfaces océan-atmosphère, océan-glace et glace-atmosphère. Nous discuterrons plus loin de leur valeur.

D'autre part l'EBM0 suppose que si la chaleur perdue à la surface de la glace est supérieure à celle gagnée à sa base, le déficit est employé à la formation de glace; la variation de h dépend donc directement de la différence entre les flux entrant et sortant $Q_{ia}$ et $Q_{oi}$ comme (conservation de la masse):


$\displaystyle {\rho_{ice}}L_f\frac{{\partial}h}{{\partial}t} = Q_{ia}-Q_{oi}$     (2.2)

Avec $L_f=3.34{\times}10^{5}\,JK^{-1}$ la chaleur latente de fusion de la glace et $\rho_{ice}=910\,kg\,m^{-3}$ sa masse volumique.

L'EBM0 suppose également que la glace est en équilibre thermodynamique et que la diffusion de la chaleur à travers elle est linéaire. Le flux de chaleur $Q_{ia}$ à sa surface s'ecrit donc aussi:


$\displaystyle Q_{ia}$ $\textstyle =$ $\displaystyle -K_{ice}\frac{{\partial}T}{{\partial}z}$  
$\displaystyle Q_{ia}$ $\textstyle =$ $\displaystyle -K_{ice}\frac{T_{z=h}-T_{z=0}}{h}$  
$\displaystyle Q_{ia}$ $\textstyle =$ $\displaystyle K_{ice}\frac{T_{fr}-T_i}{h}$  

D'où:


\begin{displaymath} T_i=T_{fr}-\frac{hQ_{ia}}{K_{ice}} \end{displaymath} (2.3)

avec $K_{ice}=2\,Wm^{-1}K^{-1}$ la conductivité thermique de la glace. On peut également exprimer $T_i$ en fonction de $T_a$ en introduisant $Q_{ia}$ de (3) :


\begin{displaymath} T_i=\frac{K_{ice}T_{fr}+hE_{ia}T_a}{K_{ice}+hE_{ia}} \end{displaymath} (2.4)

Le bilan océanique de flux de chaleur s'écrit finalement:


\begin{displaymath}Q_{oce}= {\rho_e}Cp_eH_e\frac{{\partial}T_o}{{\partial}t} = ... ... glace} \\ Q_{deep}-Q_{oa} & \mbox{sinon} \end{array}\right.\end{displaymath}

avec ${\rho_e}=1025\,kg\,m^{-3}$ la masse volumique et $Cp_e=4180\,J\,kg^{-1}\,K^{-1}$ la chaleur spécifique de l'eau de mer.

On note ${\mu}={\rho_e}Cp_eH_e$ ($JK^{-1}m^{-2}$) la capacité calorifique de la colonne d'eau de surface unité et on la considère comme constante. En effet les variations de la masse volumique de l'eau de mer sont faibles pour les plages de températures envisagées ici.

Les expressions de tous les flux et du bilan océanique forment donc le système d'équations suivants:

Analyse

Nous n'allons pas refaire ici l'analyse complète d'un tel modèle, nous invitons pour cela le lecteur à consulter directement l'article de Welander. Nous rappelerons juste que le système possède 0, 1 ou 2 états d'équilibres: un équilibre ice-cover stable et un équilibre ice-free instable aux faibles perturbations.

Une bonne représentation du système peut se faire dans le plan ($T_o$,h) où l'on peut placer les deux points représentants les équilibres ice-cover ($S_0$) et ice-free ($S_1$). Lorsque les paramètres imposés au système empêchent celui-ci d'atteindre l'un ou l'autre des états d'équilibres, le système oscille dans des états intermédiaires. Il s'agit du cas où les deux points $S_0$ et $S_1$ sont en dehors du cadrant physique ($h>0$ et $T_o>T_{fr}$) du plan ($T_o$,h) et c'est le cas qui nous intéresse ici. En effet, on peut se demander si de telles oscillations sont uniquement dues à la configuration du système où si on peut les retrouver dans un système thermodynamique plus réaliste. Pour l'instant, cherchons à caractériser le regime où l'on obtient des oscillations.

Le système se met donc à osciller lorsqu'il ne peut atteindre les états d'équilibre, c'est à dire lorsque les cas impossibles physiquement: $T_o<T_{fr}$ dans le cas sans glace (SG) et $T_i>T_{fr}$ dans le cas avec glace (AG) sont vérifiés.

A l'équilibre le bilan océanique $Q_{oce}$ est nul (dérivées temporelles nulles), donc de l'équation (10) on sort:

\begin{displaymath} T_o=\frac{1}{E_{oa}}Q_{deep}+T_a \end{displaymath} (2.10)

et des équations (7) et (8):
$\displaystyle Q_{ia}$ $\textstyle =$ $\displaystyle Q_{oi}$  
$\displaystyle E_{ia}\,(T_i\,-\,T_a)$ $\textstyle =$ $\displaystyle Q_{deep}$  
$\displaystyle T_i$ $\textstyle =$ $\displaystyle \frac{1}{E_{ia}}Q_{deep}+T_a$ (2.11)

D'où l'expression des deux critères à verifier pour obtenir des oscillations:

\fbox{ \begin{minipage}{0.7\textwidth} \begin{eqnarray} \frac{1}{E_{oa}}Q_{dee... ...}}Q_{deep}+T_a > T_{fr} &:& {Crit\\lq ere\, 2\, (AG)} \end{eqnarray}\end{minipage}}

Dans le plan ($T_a$,$Q_{deep}$) il est donc possible d'établir une zone dans laquelle le système donnera des oscillations. Nous avons mené une exploration numérique du modèle afin de le vérifier, les résultats sont présentés sur la figure 1. Le maillage de tous les ``runs'' effectués n'est pas très fin mais suffit à mettre en évidence que les critéres sont bons. Notons que les oscillations ne sont possibles que dans le plan $Q_{deep}>0$. En effet, dans le cas (AG) l'équilibre n'est pas atteint si $T_o>T_{fr}$; or d'après l'équation (7) à l'équilibre cela reviens à dire que $Q_{deep}$ doit être positif.

Notons également l'importance des coefficients d'échanges thermodynamique $E_{oa}$ et $E_{ia}$. En effet, si à termes $T_a$ et $Q_{deep}$ sont amenés à disparaître de la liste des constantes ($T_a$ dès la prochaine étape et $Q_{deep}$ dans le modèle en boîte), les coefficients d'échanges sont des paramètres physiques constants. Leur détermination est difficile c'est pourquoi nous utiliserons ici les valeurs les plus répandues dans la littérature.

Dans un régime d'oscillations, le bilan de flux de chaleurs $Q_{oce}$ n'arrive jamais à s'équilibrer. Pour indication, au point ($-5^oC$,$50Wm^{-2}$) de la figure 1, les valeurs moyennes des flux sont: $\overline{Q_{ia}}=21Wm^{-2}$, $\overline{Q_{oi}}=17Wm^{-2}$ et $\overline{Q_{oa}}=170Wm^{-2}$.

Figure: Vérification numérique des critéres d'oscillations d'EBM0. Le nombre en indice des points d'oscillations indique la période en années ( $E_{oa}=35WK^{-1}m^{-2}$, $E_{ia}=10WK^{-1}m^{-2}$ et $T_{fr} = -1.9 ^oC$).
\includegraphics[scale=0.7]{Images/fig_20b.eps}

Étape 2: $T_a$ variable (EBM1)

L'hyphothèse d'un forçage constant par l'atmosphère sur l'océan est très forte parce qu'elle implique une température atmosphèrique constante. Déplacer le forçage de l'interface océan - atmosphère à l'interface atmosphère - espace par l'intermédaire du flux solaire fixe $Q_s$ permet de rendre l'atmosphère et l'océan interactifs.

Équations

Afin de se rapprocher des conditions réelles, on utilise pour déterminer $Q_s$ la formulation (Toggweiler et al. 1998): $Q_s = \frac{S_0}{4}(1+0.7(cos^2({\lambda})-\frac{2}{3})) $ où $S_0=1360Wm^{-2}$ est la constante solaire et ${\lambda}$ la latitude. En fait, notre $Q_s$ représente le flux entrant dans l'atmosphère, il faut donc prendre en compte l'albédo au sommet de l'atmosphère ${\alpha}$. Thierry Huck dans la formulation d'un EBM pour coupler au modèle MOM3 (Huck 2001) propose la formulation: ${\alpha} = 0.2 + 0.36sin^2({\lambda})$. L'expression définitive de notre $Q_s$ est donc:

\begin{displaymath}Q_s = (1-{\alpha})\frac{S_0}{4}\left(1+0.7\left(cos^2({\lambda})-\frac{2}{3}\right)\right) \end{displaymath}

On suppose que notre colonne se situe aux hautes latitudes et en moyennant $Q_s$ entre $70^oN$ et $80^oN$ on obtient $Q_s=92Wm^{-2}$.

D'autre part, on doit prendre en compte l'émission dans les ondes longues de l'atmosphère vers l'espace. Ce nouveau flux de chaleur est de la forme $Q^{lw}={\sigma}{\epsilon}_p{T_a}^4$ avec ${\sigma}=5.67{\times}10^{-8}\,Wm^{-2}K^{-4}$ la constante de Stefan-Boltzman et ${\epsilon}_p=0.51$ l'émissivité planétaire. Afin de faciliter l'étude du système, ce terme est linéarisé en: $Q^{lw}={\sigma}{\epsilon}_p(A^4+4A^3{T_a})$ avec $A=273.15K$. Pour des facilités de notation avec les considérations qui surviendront dans les prochains paragraphes, on pose $a_T={\sigma}A^4=315.64Wm^{-2}$ et $b_T={\sigma}4A^3=4.62Wm^{-2}K^{-1}$. D'où: $Q^{lw}={\epsilon}_p(a_T+b_T{T_a})$.

Finalement, le bilan de flux de chaleur atmosphèrique $Q_{atm}$ s'écrit:

\begin{displaymath}Q_{atm}= {\rho_a}Cp_aH_a\frac{{\partial}T_a}{{\partial}t} = ... ...} \\ Q_s+Q_T+Q_{oa}-Q^{lw} & \mbox{sinon} \end{array}\right.\end{displaymath}

avec ${\rho_a}=1.225 \,kg\,m^{-3}$ la masse volumique et $Cp_a=1004\,J\,kg^{-1}\,K^{-1}$ la chaleur spécifique de l'atmosphère, $Q_T$ ($Wm^{-2}$) un transport de chaleur horizontal (représentant par exemple un apport de chaleur des latitudes équatoriales aux latitudes p ô laires) et $Q_s$ ($Wm^{-2}$) le flux de chaleur net entrant au sommet de l'atmosphère. On note ${\gamma}={\rho_a}Cp_aH_a$ en $JK^{-1}m^{-2}$ la capacité calorifique de la colonne d'air de surface unité et on la considère comme constante.

Les expressions de tous les flux et des bilans océaniques et atmosphèrique forment donc le nouveau système d'équations (à une inconnue suplémentaire $T_a$):

Analyse

L'équation (17) est indépendante du reste du système, on peut donc déterminer la température océanique d'équilibre dans le cas (AG):

\begin{displaymath}T_o = T_{fr} + \frac{Qdeep}{E_{oi}}\end{displaymath}

Le systéme (AG) possède donc 2 équations, 2 inconnues ($T_a$,h); tout comme le système (SG) où les inconnues sont ($T_a$,$T_o$).

Pour déterminer la nouvelle formulation des critères pour obtenir des oscillations entre les deux états d'équilibres, il faut dans chaque cas calculer $T_a$ puis $T_o$ (SG) ou $T_i$ (AG):

La température atmosphérique d'équilibre ne dépend donc pas de la présence ou non de glace à la surface océanique. Elle dépend uniquement des flux entre la totalité de la colonne et l'extérieur.

L'expression des deux critères (SG:$T_o<T_{fr}$ et AG:$T_i>T_{fr}$) à verifier est donc:


$\displaystyle \left(\frac{1}{E_{oa}}+\frac{1}{{\epsilon}_pb_T}\right)Q_{deep}+ ... ...}{{\epsilon}_pb_T}Q_T + \frac{1}{{\epsilon}_pb_T}Q_s - \frac{a_T}{b_T} < T_{fr}$ $\textstyle :$ $\displaystyle {Crit\\lq ere\, 1\, (SG)}$ (2.22)
$\displaystyle \left(\frac{1}{E_{ia}}+\frac{1}{{\epsilon}_pb_T}\right)Q_{deep}+ ... ...1}{{\epsilon}_pb_T}Q_T +\frac{1}{{\epsilon}_pb_T}Q_s - \frac{a_T}{b_T} > T_{fr}$ $\textstyle :$ $\displaystyle {Crit\\lq ere\, 2\, (AG)}$ (2.23)

De la même manière qu'avec l'EBM0, nous avons exploré numériquement ce modèle pour vérifier les critères d'oscillations (figure 2). Comme nous supposons que la colonne de fluide est située dans les hautes latitudes, les flux de chaleurs $Q_{deep}$ et $Q_T$ doivent être positifs pour signifier un apport de chaleur et donc un réchauffement de la colonne par les eaux et l'atmosphère adjacentes. Nous avons donc ``zoomé''dans la figure 2 pour explorer le plan $Q_T>0$,$Q_{deep}>0$ (figure 3). L'exploration numérique confirme bien la véracité des critères (26) et (27).

Il peut être intéressant de regarder comment évolue la période dans la zone des oscillations en fixant soit $Q_{deep}$ soit $Q_T$. On se place donc dans une petite zone et on relève la période des oscillations (figure 4). On trace ensuite la période en fonction de $Q_T$ pour $Q_{deep}=30Wm^{-2}$ (figure 5). Le caractère non linéaire des oscillations apparaît clairement: de manière non symétrique, plus on s'éloigne des limites au delà desquelles on a équilibre, plus la période est courte. Le profil de la période se fait de la même manière suivant que l'on soit à $Q_{deep}$ fixé ou $Q_T$ fixé.

Comparer les oscillations de l'EBM1 avec celles de l'EBM0 n'est pas évident puisqu'un changement fondamental s'est opéré de l'un à l'autre ($T_a$ varie) et que l'on n'exprime plus les critères d'oscillations dans le même plan. Toutefois si on regarde pour le même forçage $Q_{deep}=30Wm^{-2}$, EBM0 donne des oscillations de période environ 5.5 ans pour $T_a$ compris entre $-4.9^oC$ et $-2.8^oC$ et EBM1 donne des oscillations de période environ 50 ans pour $Q_T$ compris entre $27Wm^{-2}$ et $32Wm^{-2}$ où la température atmosphèrique moyenne est $\overline{T_a}=-3.8^oC$. On constate que $\overline{T_a}$ est situé au milieu de la plage d'oscillations de EBM0, on peut donc supposer que l'on a approximativement les mêmes forçages dans les deux modèles ce qui permet une comparaison.

La différence la plus frappante est dans l'augmentation de la période qui est presque multipliée par un facteur 10. Rendre la température atmosphèrique réactive à celle de l'océan ou de la glace engendre une forte inertie dans les oscillations. Dans les différentes phases (augmentation puis diminution de $T_o$) la réaction de l'atmosphère à la température de surface du sol (glace ou océan) ralentie le cheminement vers l'équilibre impossible. Le forçage constant de EBM0 ne présente pas cette inertie donc la période est plus courte.

En supposant donc que l'on se situe dans des forçages plus ou moins identiques dans les deux modèles, on peut tenter de comparer les oscillations de manière adimensionnelle. Sur la figure 6 on a représenté l'évolution de la température océanique durant une période pour les deux modèles. La phase d'augmentation de $T_o$ (du point A aux points B et C: qui correspond à la présence de glace sur la colonne) est plus rapide dans l'EBM1 que l'EBM0 et inversement, la phase de diminution de $T_o$ (des points B et C à D: absence de glace) est plus rapide dans l'EBM0 que l'EBM1, ces vitesse étant relative à une période. Cela signifie qu'avec les mêmes paramètres physiques, l'interactivité de la température atmosphèrique a le même effet virtuel qu'une diminution du coefficient d'échange thermodynamique ocean-atmosphère $E_{oa}$. En effet en diminuant celui-ci on rend plus difficile les échanges et donc on augmente les temps de transition entre chaque états du système. Pour confirmer cela nous avons tracé la courbe obtenue avec l'EBM0 en augmentant $E_{oa}$ (courbe en tiret-point). Augmenter les échanges océan-atmosphère, diminue le temps relatif sur une oscillation où $T_o$ diminue parce que l'océan est directement refroidie par l'atmosphère. Pour retrouver des périodes d'oscillations du même ordre de grandeur (5 à 15 ans) entre EBM1 et EBM0, il faudrait donc augmenter $E_{oa}$.

Figure: Vérification numérique des critéres d'oscillations d'EBM1. ( $Q_s=92Wm^{-2}, E_{oi}=10WK^{-1}m^{-2}, E_{ia}=10WK^{-1}m^{-2}, E_{oa}=35WK^{-1}m^{-2}, et T_{fr} = -1.9 ^oC$).
\includegraphics[scale=0.6]{Images/fig_21ba.eps}

Pour illustrer le fait que les oscillations sont obtenues lorsque le système ne peut atteindre aucun état d'équilibre, nous avons tracé sur la figure 7 l'évolution du bilan radiatif océanique $Q_{oce}$ en fonction du temps. Les artefacts situés à environ 15 et 65 ans sont dues à la transition discontinue entre chaques modes (ice-free et ice-cover). S'il éxistait un équilibre quelconque, une fois celui-ci atteint le bilan deviendrait nul mais ce n'est pas le cas lors des oscillations.

Du point A au point B, le bilan est négatif et il n'y a pas de glace ce qui signifie que l'océan se refroidit et qu'il tente d'atteindre un équilibre ice-free ($Q_{oce}$ augmente vers zero). Mais il atteint la température de congélation avant la température d'équilibre supposée, il se forme donc de la glace et le bilan devient aussitôt positif (on rappelle que le modèle est binaire). Alors l'océan se réchauffe (du point C au point D) vers une nouvelle température d'équilibre ($Q_{oce}$ diminue vers zero) mais la glace disparaît entièrement avant et le bilan devient négatif. En excluant les équilibres ice-free et ice-cover de la réalité physique on fait osciller le système d'un état à l'autre indéfiniment.
Enfin, Welander a montré que l'équilibre ice-cover est stable et l'équilibre ice-free instable (l'ajout momentané d'une petite quantité de glace perturbe le système ice-free) ce qui suppose une dépendance de l'état final avec les conditions initiales. Qu'en est-il dans un régime oscillatoire ? Nous avons tracé dans le plan ($T_o$,h) ce que donnait l'EBM1 pour des forçages identiques mais des conditions initiales différentes (figure 8). Chaque cycle obtenue se superpose parfaitement aux autres, le régime oscillatoire est donc complètement indépendant des conditions initiales.

Figure: Vérification numérique des critéres d'oscillations d'EBM1 dans le plan $Q_T>0$ et $Q_{deep}>0$ ( $Q_s=92Wm^{-2}, E_{oi}=10WK^{-1}m^{-2}, E_{ia}=10WK^{-1}m^{-2}, E_{oa}=35WK^{-1}m^{-2} \,et\, T_{fr} = -1.9 ^oC$).
\includegraphics[scale=0.6]{Images/fig_21bb.eps}

Figure: Période des oscillations en années avec EBM1( $Q_s=92Wm^{-2}, E_{oi}=10WK^{-1}m^{-2}, E_{ia}=10WK^{-1}m^{-2}, E_{oa}=35WK^{-1}m^{-2} \,et\, T_{fr} = -1.9 ^oC$).
\includegraphics[scale=0.6]{Images/fig_22ba.eps}

Figure: Evolution de la période sur la coupe $Q_{deep}=30Wm^{-2}$ avec EBM1( $Q_s=92m^{-2}, E_{oi}=10WK^{-1}m^{-2}, E_{ia}=10WK^{-1}m^{-2}, E_{oa}=35WK^{-1}m^{-2} \,et\, T_{fr} = -1.9 ^oC$).
\includegraphics[scale=0.7]{Images/fig_22bb.eps}

Figure: Evolution adimensionnelle de $T_o$ pour EBM0 et EBM1 pendant une période ( $Q_{deep}=30Wm^{-2}, Q_s=92.38Wm^{-2}, E_{oi}=10WK^{-1}m^{-2}, E_{ia}=10WK^{-1}m^{-2} \,et\, T_{fr} = -1.9 ^oC$).
\includegraphics[scale=0.7]{Images/fig_24b.eps}

Figure: Evolution temporelle de $Q_{oce}$ de EBM1 dans un régime oscillatoire. ( $Q_{deep}=30Wm^{-2}$,$Q_T=30Wm^{-2}$, $Q_s=92.38Wm^{-2}$, $E_{oi}=10WK^{-1}m^{-2}$, $E_{ia}=10WK^{-1}m^{-2}$, $E_{oa}=35WK^{-1}m^{-2}$et $T_{fr} = -1.9 ^oC$).
\includegraphics[scale=0.7]{Images/fig_25.eps}

Figure: Diagramme ($T_o$,h) obtenus avec EBM1 pour différentes conditions initiales ( $Q_{deep}=28Wm^{-2}, Q_T=31Wm^{-2}, Q_s=92Wm^{-2}, E_{oi}=10WK^{-1}m^{-2}, E_{ia}=10WK^{-1}m^{-2}, E_{oa}=35WK^{-1}m^{-2} \,et\, T_{fr} = -1.9 ^oC$).
\includegraphics[scale=0.7]{Images/fig_23.eps}

Étape 3: les ondes courtes (EBM2)

Maintenant que l'atmosphère réagit au reste de la colonne de fluide, il convient de prendre en compte les réflexions et absorptions des ondes courtes afin d'établir un bilan énergétique plus satisfaisant.

Équations

Si $Q_s$ est le rayonnement solaire incident entrant au sommet de la colonne, alors ${\tau}Q_s$ est la partie qui atteind effectivement le sol après absorption par l'atmosphère (${\tau}=0.65$: coefficient d'absorption de l'atmosphère).

On ne prend pas en compte les multiples réflexions entre le sol et l'atmosphère, on suppose simplement que le flux $Q_s$ entre dans l'atmosphère, se réfléchie contre le sol et repart vers l'espace. Dans ces conditions les flux d'ondes courtes absorbés par la glace et l'océan sont:

\begin{displaymath}Q_{ia}^{sw}=(1-{\alpha}_i){\tau}Q_s\end{displaymath}


\begin{displaymath}Q_{oa}^{sw}=(1-{\alpha}_o){\tau}Q_s\end{displaymath}

avec ${\alpha}_o=0.06$ l'albédo de l'océan et ${\alpha}_i=0.7$ celui de la glace. La partie réfléchie par le sol ${\alpha}{\tau}Q_s$ vers le haut est de nouveau absorbée par l'atmosphère et les flux sortant vers l'espace sont:

\begin{displaymath}Q^{esc}_o={\alpha}_o{\tau}^2Q_s\end{displaymath}


\begin{displaymath}Q^{esc}_i={\alpha}_i{\tau}^2Q_s\end{displaymath}

Ainsi, les bilans aux interfaces océan-atmosphère et glace-atmosphère changent mais pas le bilan à l'interface océan-glace. Avec l'expresion du nouveau bilan atmosphèrique, on obtient le nouveau système d'équations:

Analyse

La résolution des deux systèmes donne pour les températures d'équilibre:

D'où les deux nouveaux critères d'oscillations:


$\displaystyle \left(\frac{1}{E_{oa}}+\frac{1}{{\epsilon}_pb_T}\right)Q_{deep}+ ... ...o)}{E_{oa}}\right)Q_s + \frac{1}{{\epsilon}_pb_T}Q_T - \frac{a_T}{b_T} < T_{fr}$ $\textstyle :$ $\displaystyle (SG)$ (2.34)
$\displaystyle \left(\frac{1}{E_{ia}}+\frac{1}{{\epsilon}_pb_T}\right)Q_{deep}+ ... ..._i)}{E_{ia}}\right)Q_s +\frac{1}{{\epsilon}_pb_T}Q_T - \frac{a_T}{b_T} > T_{fr}$ $\textstyle :$ $\displaystyle (AG)$ (2.35)

Le changement par rapport à l'EBM1 se situe dans les coefficients de $Q_s$ où l'on voit apparaître l'absorption de l'atmosphère et les albédos de l'océan et de la glace.

Nous avons exploré numériquement le modèle pour vérifier les critéres dans le plan $Q_{deep}$ et $Q_T$ positifs. Les critères sont bien vérifiés mais une zone supplémentaire d'oscillations apparaît (figure 9). Il semble que le deuxième critère correspondant à l'équilibre avec glace soit plus faible que prévu. La délimitation de cette nouvelle zone semble indiquer que le critère 2 est bon dans son coefficient directeur mais nécessite une correction dans son ordonnée à l'origine.

Contrairement à l'EBM1, la température atmosphérique d'équilibre (dans le cas d'absence d'oscillations) dépend cette fois de la présence où non de glace à l'interface.

Dans les régimes d'équilibre, chaque variable vérifie une équation différentielle du deuxième ordre à coefficients constants et sans second membre. Il serait donc possible d'obtenir une formulation analytique de la période des oscillations en calculant le temps mis pour chaque régime transitoire pour atteindre la valeur critique qui fait changer le système d'équations (avec ou sans glace). Par exemple dans le cas sans glace, la température océanique vérifie l'équation:


\begin{displaymath}\frac{{\partial}^2{T_o}}{{\partial}t} + \left( \frac{E_{oa}}{... ...\partial}t} + \frac{E_{oa}{{\epsilon}_pb_T}}{{\gamma}{\mu}}T_o \end{displaymath}

On calcule ainsi le temps mis par $T_o$ pour atteindre une température d'équilibre $T_o^{eq}$:


\begin{displaymath}t &=& \left\vert\frac{-1}{r_1}ln\left(\frac{T_o-T_o^{eq}}{T_o(0)-T_o^{eq}}\right)\right\vert \end{displaymath}

Avec $r_1$ la racine positive du polynôme caractéristique de l'équation différentielle et en prenant $T_o^{eq} = T_{fr}$ et $T_o(0)$ la température océanique au moment où la glace avait disparu au précédant changement d'équations. Nous avons testé cette formule et elle semble fonctionner correctement; elle donne en fait la semi-période correspondant à la zone de diminution de $T_o$ de la figure 6.

Toutefois, même s'il est possible (quoique bien plus difficile) d'obtenir une formulation de l'autre semi-période (correspondant au cas avec glace), l'expression finale ne serait aisément analysable du point de vue de l'influence des différents paramètres.

Nous avons ensuite voulu comparer l'EBM1 et l'EBM2 du point de vue des périodes d'oscillations pour voir quel était l'effet des ondes courtes. Pour $Q_{deep}=30Wm^{-2}$ nous avons donc relevé les périodes obtenues pour chaques modèles, les résultats sont présentés sur la figure 10. Le caractère non linéaire de la période semble conservé à l'identique. La différence vient d'une baisse sensible des valeurs de la période. L'introduction de la température atmosphérique dans l'EBM1 avait introduit une certaine inertie dans les oscillations (et ainsi avait augmenté les périodes) mais la prise en compte explicite des réflexions-absorptions des ondes courtes paraît rediminuer ces périodes. Les interactions d'ondes courtes favorisent les échanges aux interfaces et accélèrent donc les oscillations.

Figure: Vérification numérique des critéres d'oscillations d'EBM2 dans le plan $Q_T>0$ et $Q_{deep}>0$ ( $Q_s=92.38Wm^{-2}, E_{ia}=4WK^{-1}m^{-2}$).

\includegraphics[scale=0.6]{Images/fig_28b.eps}
Figure: Comparaison des périodes obtenues avec l'EBM1 et l'EBM2 pour $Q_{deep}=30Wm^{-2}$ et $E_{ia}=4WK^{-1}m^{-2}$.

\includegraphics[scale=0.6]{Images/fig_27a.eps}

Étape finale 4: les ondes longues (EBM3)

Nous l'avons vu dans le paragraphe 2.2.1 nous modélisons les émissions d'ondes longues comme celles d'un corps noir suivant la loi de Stefan.

Équations

Aux interfaces, la glace et l'océan rayonnent (verticalement vers le haut) vers l'atmosphère respectivement: $Q_i^{lw}={\sigma}{\epsilon}_i{T_i}^4$ et $Q_o^{lw}={\sigma}{\epsilon}_o{T_o}^4$ et l'atmosphère rayonne (verticalement vers le bas): $Q_a^{lw}={\sigma}{\epsilon}_a{T_a}^4$. Les bilans radiatifs aux interfaces donnent donc: $Q_{ia}^{lw}={\sigma}({\epsilon}_i{T_i}^4-{\epsilon}_a{T_a}^4)$ et $Q_{oa}^{lw}={\sigma}({\epsilon}_o{T_o}^4-{\epsilon}_a{T_a}^4)$. Ce qui, avec les notations déjà employées, s'écrit:


\begin{displaymath}Q_{ia}^{lw} &=& \left\{\begin{array}{ll} {\epsilon}_i(a_T+b_... ...}_a) + b_T({\epsilon}_iT_i-{\epsilon}_aT_a) \end{array}\right.\end{displaymath}


\begin{displaymath}Q_{oa}^{lw} &=& \left\{\begin{array}{ll} {\epsilon}_o(a_T+b_... ...}_a) + b_T({\epsilon}_oT_o-{\epsilon}_aT_a) \end{array}\right.\end{displaymath}

Avec ${\epsilon}_o=0.96$, ${\epsilon}_i=0.97$ et ${\epsilon}_a=0.84$ les émissivités de l'océan, de la glace et de l'atmosphère. L'émission de l'amosphère vers l'espace a déjà été prise en compte comme: $Q^{lw}={\epsilon}_p(a_T+b_T{T_a})$.

Les nouveaux bilans océanique et atmosphérique sont donc:

Analyse

On calcule les critères d'oscillations comme précédemment et on obtient:


$\displaystyle \left(\frac{1}{E_{oa}+b_T{\epsilon}_o} + \frac{1}{{\epsilon}_pb_T}\frac{E_{oa}+b_T{\epsilon}_a}{E_{oa}+b_T{\epsilon}_o}\right)Q_{deep}+$      
$\displaystyle \left(\frac{1-{\tau}^2{\alpha}_o}{{\epsilon}_pb_T}\frac{E_{oa}+b_... ...T{\epsilon}_o}+ \frac{{\tau}(1-{\alpha}_o)}{E_{oa}+b_T{\epsilon}_o}\right)Q_s +$      
$\displaystyle \left(\frac{1}{{\epsilon}_pb_T}\frac{E_{oa}+b_T{\epsilon}_a}{E_{oa}+ b_T{\epsilon}_o}\right)Q_T-$      
$\displaystyle \left(\frac{a_T}{b_T}\frac{E_{oa}+b_T{\epsilon}_a}{E_{oa}+b_T{\ep... ...+ \frac{a_T({\epsilon}_o-{\epsilon}_a)}{E_{oa}+b_T{\epsilon}_o}\right) < T_{fr}$ $\textstyle :$ $\displaystyle (SG)$ (2.43)


$\displaystyle \left(\frac{1}{E_{ia}+b_T{\epsilon}_i} + \frac{1}{{\epsilon}_pb_T}\frac{E_{ia}+b_T{\epsilon}_a}{E_{ia}+b_T{\epsilon}_i}\right)Q_{deep}+$      
$\displaystyle \left(\frac{1-{\tau}^2{\alpha}_i}{{\epsilon}_pb_T}\frac{E_{ia}+b_... ...{\epsilon}_i} + \frac{{\tau}(1-{\alpha}_i)}{E_{ia}+b_T{\epsilon}_i}\right)Q_s +$      
$\displaystyle \left(\frac{1}{{\epsilon}_pb_T}\frac{E_{ia}+b_T{\epsilon}_a}{E_{ia}+ b_T{\epsilon}_i}\right)Q_T-$      
$\displaystyle \left(\frac{a_T}{b_T}\frac{E_{ia}+b_T{\epsilon}_a}{E_{ia}+b_T{\ep... ...+ \frac{a_T({\epsilon}_i-{\epsilon}_a)}{E_{ia}+b_T{\epsilon}_i}\right) > T_{fr}$ $\textstyle :$ $\displaystyle (AG)$ (2.44)

Nous avons biensûr vérifié ces critères et ils semblent être corrects. Toutefois la zone où l'on obtient des oscillations (pour les forçages classiques: $Q_s=92.38Wm^{-2}$, $E_{ia}=4WK^{-1}m^{-2}$) n'est pas située dans le plan ( $Q_{deep}>0,Q_T>0$) habituel mais dans le plan ( $Q_{deep}>0,Q_T<0$) (figure 11). Un transport atmosphérique de chaleur négatif signifie que la colonne d'air exporte de la chaleur ce qui vient contredire l'hypothèse suivant laquelle la colonne est située dans les hautes latitudes. La manière la plus simple pour ramener les oscillations dans le plan ( $Q_{deep}>0,Q_T>0$) est d'abaisser la quantité $Q_s$ jusqu'à $35Wm^{-2}$. Cette valeur est discutable mais l'objectif ici étant de faire une première approche de ces oscillations sans atteindre un degré de réalisme extrème, nous l'accepterons.

Pour ce nouveau forçage les critères d'oscillations sont biensûr toujours vérifiés (figure 12). L'exploration numérique du modèle a montré que les oscillations avaient une période beaucoup plus longue que celle des modèles précédents, de l'ordre de la centaine d'années (97.63 ans pour $Q_T=31Wm^{-2},Q_{deep}=86Wm^{-2}$ ce qui constitue presque un minimum pour l'ensemble des oscillations). Ainsi, l'ajout des rayonnements basses fréquences introduit dans les oscillations une très grande inertie.

Figure: Vérification numérique des critéres d'oscillations d'EBM3 pour le forçage habituel: $Q_s=92.38Wm^{-2}, E_{ia}=4WK^{-1}m^{-2}$.
\includegraphics[scale=0.7]{Images/fig_30.eps}

Figure: Vérification numérique des critéres d'oscillations d'EBM3 pour le forçage: $Q_s=35Wm^{-2}, E_{ia}=4WK^{-1}m^{-2}$.
\includegraphics[scale=0.5]{Images/fig_31.eps}

Extension du modèle à la fraction de glace

Jusqu'ici nous avons considéré un modèle binaire, c'est à dire que la surface de l'océan était entièrement ou pas du tout recouverte de glace. Pour plus de réalisme nous avons introduit une fraction de glace qui permet à la surface de l'océan de n'être recouverte que partiellement. Nous nous sommes appuyé sur le modèle de Bjornsson (H.Bjornsson, Mars 2000) ce pourquoi nous ne donnerons ici que les résultats sans entrer dans le détails des calculs.

L'introduction de cette fraction permet de s'affranchir des distinctions entre les cas avec ou sans glace, nous n'avons plus qu'un seul jeu d'équations.

Les flux de chaleur aux interfaces et les émissions/absorption d'ondes longues/courtes sont inchangés:


\begin{displaymath}Q_{ia} = E_{ia}(T_i-T_a)-{\tau}(1-{\alpha}_i)Q_s+a_T({\epsilon}_i-{\epsilon}_a)+b_T({\epsilon}_iT_i-{\epsilon}_aT_a)\end{displaymath}


\begin{displaymath}Q_{oa} = E_{oa}(T_o-T_a)-{\tau}(1-{\alpha}_o)Q_s+ a_T({\epsilon}_o-{\epsilon}_a)+b_T({\epsilon}_oT_o-{\epsilon}_aT_a)\end{displaymath}


\begin{displaymath}Q_{oi} = E_{io}(T_o-T_{fr}) \end{displaymath}


\begin{displaymath}Q^{lw} = {\epsilon}_p(a_T+b_T{T_a}) \end{displaymath}


\begin{displaymath}Q_i^{esc} = {\alpha_i}{\tau}^2Q_s \end{displaymath}


\begin{displaymath}Q_o^{esc} = {\alpha_o}{\tau}^2Q_s \end{displaymath}

Le changement intervient dans les bilans de flux:


$\displaystyle Q_{atm}={\gamma}\frac{{\partial}T_a}{{\partial}t}$ $\textstyle =$ $\displaystyle Q_s+Q_T+ A\,Q_{ia} + (1-A)\,Q_{oa} - \left(A\,Q_i^{esc}+(1-A)Q_o^{esc}\right)-Q_{lw}$ (2.45)
$\displaystyle Q_{oce}={\mu}\frac{{\partial}T_o}{{\partial}t}$ $\textstyle =$ $\displaystyle Q_{deep}-A\,Q_{oi}-(1-A)\,Q_{oa}\,{\aleph}(T_o-T_{fr})$ (2.46)
$\displaystyle {\rho_{ice}}L_f\frac{{\partial}h}{{\partial}t}$ $\textstyle =$ $\displaystyle \left\{\begin{array}{ll} A\,(Q_{ia}-Q_{oi})+(1-A)Q_{oa} & {si (T_o=T_{fr}\,et\,Q_{oa}>0)} \\ A\,(Q_{ia}-Q_{oi}) & sinon \end{array}\right.$ (2.47)

Il faut également prendre en compte les variations de A:
$\displaystyle {\rho_{ice}}L_f\frac{{\partial}A}{{\partial}t}$ $\textstyle =$ $\displaystyle \left\{\begin{array}{ll} (1-A)\frac{Q_{oa}}{H_0} & {si (T_o=T_{fr... ...}t} & {si \frac{{\partial}h}{{\partial}t}<0} \\ 0 & {sinon} \end{array}\right.$ (2.48)

Avec $H_0=0.1m$ la limite entre la glace fine et épaisse, ${\aleph}(T_o-T_{fr})$ la fonction de Heavyside qui est nulle si $T_o<T_{fr}$ (dans ce cas, toute la chaleur perdue par l'océan va à former de la glace) et A la fraction de glace comprise entre 0 et 1.

La manipulation de ce système d'équations étant beaucoup moins facile que celle des précédents nous n'avons ici que sommairement exploré le modèle. Nous avons retrouvé des régimes d'oscillations avec tous les modèles et les résultats comparatifs entre modèles sont quasiment les mêmes. Toutefois nous avons observé des oscillations particulières avec l'EBM3. Il semble exister une zone où les oscillations interviennent dans l'état ice-cover sans le quitter. Le système oscille mais n'atteind jamais l'état sans glace (voir figure 13). Le régime semble dépendre des conditions initiales, du moins pour la partie transitoire avant établissement du régime stable d'oscillations.

Nous avons cherché numériquement à résoudre le système tout en introduisant des marqueurs dans les équations pour étudier l'influence des paramètres (voir annexe). Il apparaît bien une zone où les valeurs propres sont positives mais les résultats ne sont pas satisfaisants et ce types d'oscillations mériteraient une étude plus approfondie que nous n'avons pas eu le temps de faire ici. Notre première approche analytique nous a toutefois permis de superposer les résultats théoriques à ceux expérimentaux du modèle. Dans le plan ($h,T_o$) nous avons tracé les oscillations du modèle et celles obtenues après calcul par la méthode des perturbations (figure 14) et la concordance semble prometteuse pour une étude plus pousser en introduisant la non-linéarité des équations.

Figure: Oscillations dans l'état ice-cover avec l'EBM3 et la fraction de glace pour les paramètres: $E_{ia}=10WK^{-1}m^{-2},E_{io}=10WK^{-1}m^{-2},E_{oa}=35WK^{-1}m^{-2}$. La courbe noir représente h, la bleue $T_o$, la verte $T_i$ et la rouge $T_a$.
\includegraphics[scale=0.7]{Images/fig_12.eps}

Figure: $Q_T=60Wm^{-2},Q_{deep}=40Wm^{-2},Q_s=35Wm^{-2},E_{ia}=10WK^{-1}m^{-2},E_{io}=10WK^{-1}m^{-2},E_{oa}=35WK^{-1}m^{-2}$.
\includegraphics[scale=0.7]{Images/fig_29.eps}

Première approche du modèle en boîte

Nous avons obtenu des oscillations avec tous les modèles, que ce soit en mode binaire ou à fraction de glace. Mais jusqu'ici nous ne considérions qu'une colonne de fluide sans réelles interactions avec son environnement immédiat. Nous avons donc voulu tester nos équations dans le cas d'un modèle en boîte simple.

Equations

Nous considérons 4 boîtes: 2 représentants le pôle et 2 autres l'équateur. Dans chaque zone nous avons une boîte océanique et une atmosphérique. Les interactions entre les boîtes se font comme:

Tous les flux et les dimensions des boîtes sont représentés sur la figure 15.

On modélise la diffusion et l'advection comme:

\begin{displaymath} Q^{dif}=\frac{K_H}{L^2}{\Delta}T \\ Q^{adv}=\frac{q}{M}{\Delta}T \end{displaymath} (3.1)

Avec $ K_H(m^2s^{-1}) $ le coefficient de diffusion turbulente (qui dépend du fluide), ${\Delta}T$ la différence de température, L(m) la largeur en lattitude de la boîte, q( $kg{\cdot}s^{-1}$) le débit massique et M(kg) la masse de la boîte.

Les dérivées temporelles des différentes variables s'écrivent sous la forme: