Si l’on prend en exemple l’atmosphère, son évolution est régie par des équations de la mécanique des fluides, en l’occurrence l’équation de Navier-Stokes (ou l’équation d’Euler, suivant les situations), comprenant les équations de conservation bien classiques en physique : conservation de la masse, du moment cinétique, ou de l’énergie. Ces équations sont des équations d’évolution dans le temps, et si on connaît la solution de l’équation à un instant donné, on peut en déduire la solution à tout instant dans le futur.

Les variables de ce type de modèle sont généralement la vitesse de l’air (ou de l’eau, dans le cas d’un océan), la pression (de l’air, atmosphérique, ou de l’eau pour un océan), la température, la densité, ...

En supposant toutes ces variables connues à un instant donné (typiquement l’instant présent), on peut espérer qu’en résolvant ces équations à partir de la donnée initiale connue, on obtienne l’évolution dans le futur de ces variables, et donc du système météorologique considéré : les fameuses prévisions météo !

Malheureusement, ces systèmes sont instables, voire chaotiques : c’est le fameux effet papillon de Lorenz [1]. Ce météorologue fut l’un des premiers à remarquer que le système météorologique était fortement instable et non-linéaire. Il a développé un modèle (connu sous le nom de modèle de Lorenz) d’apparence très simple, mais en réalité extrêmement complexe, permettant d’étudier le caractère chaotique d’un système. En prenant deux points de départ très proches, ce modèle conduit souvent à des situations radicalement différentes en un temps assez court. Selon ses propos en 1963, le battement d’aile d’une mouette à un certain endroit du globe pourrait modifier inéluctablement l’évolution du système météorologique dans le futur [2]. Un peu plus tard, une version plus connue de sa théorie posera la question suivante : est-ce que le battement d’aile d’un papillon au Brésil peut causer une tornade au Texas ? La réponse étant théoriquement oui, vu l’instabilité du système météorologique, mais fort heureusement pour nous (et surtout pour le Texas !!), cela n’arrive pas à chaque fois, c’est même extrêmement rare à cette échelle de perturbation atmosphérique.

Regardons l’image ci-dessous. Nous avons choisi deux points de départ, matérialisés par un petit rond bleu et un petit rond rouge, très très proches l’un de l’autre. Comme vous pouvez le constater, au départ ($t=0$), les deux points semblent confondus. Et nous avons résolu le système de Lorenz à partir de ces deux points, pour obtenir deux trajectoires - l’une en bleu, l’autre en rouge. Si vous regardez attentivement, à partir de $t=9$, les deux trajectoires commencent à être légèrement différentes. Et à partir de $t=15$, elles sont clairement très différentes, et chacune semble vivre sa vie.

Illustration du caractère chaotique (et imprévisible !) du système de Lorenz

C’est l’illustration parfaite du caractère chaotique, et donc imprévisible, du système de Lorenz : à partir de deux points de départ, aussi proches soient-ils, on obtiendra à partir d’un certain temps deux trajectoires totalement différentes !

Eh bien le système météorologique est lui aussi imprévisible, en un certain sens. Mais à une échelle de temps qui nous intéresse le plus souvent (1 à quelques jours), il ne l’est pas tant que ça. Fort heureusement, sinon il serait impossible de savoir quel temps il fera demain !! Par ailleurs, la dimension du système météorologique (ou océanographique) est bien supérieure à la dimension 3 du système de Lorenz, et en pratique, cela lui confère une plus grande stabilité — certes toute relative [3].

De ce fait, l’évolution dans le temps de ces systèmes est très sensible à la condition initiale. Par exemple, il arrive que le temps soit le même (ou semble être le même !) à plusieurs jours ou mois d’intervalle (ou à la même période, certaines années), mais malgré cela, on voit bien que l’évolution de la météorologie les jours suivants diffère et ne conduit pas aux mêmes situations. Ce sont les petites différences qui existent entre les deux situations de départ, même si elles ne semblent pas visibles pour le commun des mortels avec par exemple des différences de température inférieures à 0,1 degré, qui conduisent à des scenarios totalement différents quelques jours après [4].

Cela est particulièrement vrai dans certaines zones très sensibles, comme les massifs de montagnes, ou les zones côtières, pour lesquels les variations brutales de relief, de vent, ou d’humidité (voire les trois en même temps !), augmentent considérablement l’instabilité. La végétation a également un impact souvent plus important qu’on pourrait le croire, influant localement sur l’humidité et la température.

Comme il n’est pas possible de connaître très précisément l’état de l’atmosphère à l’instant présent (il y a finalement très peu de zones sur la Terre, et encore moins dans le ciel, où des mesures de température ou pression sont réalisées, mais en plus, dans les zones observées, les erreurs de mesure ne sont pas toujours négligeables), la donnée initiale ne permet pas d’obtenir une solution satisfaisante.

Pour permettre de faire une prévision fiable sur une longue période de temps, il convient donc de connaître la condition initiale du système (par exemple l’état précis de l’atmosphère au jour présent) avec une grande précision, bien supérieure à celle des appareils de mesure, et surtout il faut connaître l’état de l’atmosphère à beaucoup plus d’endroits, sur toute la surface du globe et en altitude, que simplement les zones observées.

Pour cela, on va utiliser les observations sur une certaine période de temps, passée, de façon à trouver la condition initiale (au temps présent) telle que la trajectoire calculée par le modèle passe le plus près possible de ces observations. On pourra alors espérer que la trajectoire continuera à être proche de la réalité dans le futur, et ainsi que les prévisions numériques de l’évolution du système météorologique ou océanique seront les plus fiables et précises possible.

L’histoire des prévisions météorologiques est évidemment très riche, et on retrouve trace de prévisions — ou au moins de mesures récurrentes — de phénomènes particuliers (mousson, sécheresse, ...) dans l’antiquité. Mais c’est principalement après la seconde guerre mondiale que les choses se sont accélérées avec l’apparition des ordinateurs. A partir de 1944, John von Neumann (entre autres) participe au programme ENIAC (Electronic Numerical Integrator and Computer), consistant à concevoir et fabriquer un supercalculateur. Une fois l’ordinateur opérationnel, il s’agissait de lui trouver une utilité ! Et von Neumann eut l’idée de le mettre à contribution pour les prévisions météorologiques, et c’est ainsi qu’en 1950, avec des météorologistes américains, il réalise la première simulation numérique pour prédire la météo.

Depuis, la situation a bien évolué, grâce à l’utilisation d’ordinateurs de plus en plus performants et la disponibilité croissante d’observations (notamment avec l’envoi massif de satellites d’observation de la Terre dans les années 90). Et nous allons faire un petit tour d’horizon des méthodes et algorithmes les plus souvent utilisés pour les prévisions météorologiques et océanographiques ces trente dernières années.

Algorithme 4D-VAR

L’algorithme 4D-VAR (VAR comme variationnel, et 4D car il prend en compte les 3 dimensions de notre espace ainsi que la dimension temporelle) est le plus connu parmi les méthodes variationnelles, et il est actuellement utilisé de façon opérationnelle dans plusieurs centres de prévision, comme Météo-France. Il est issu de la théorie mathématique du contrôle optimal, principalement développée par Jacques-Louis Lions dans les années 60, et introduite dans le contexte de l’assimilation de données en météorologie et océanographie dans les années 80 par François-Xavier Le Dimet et Olivier Talagrand.

La méthode variationnelle est basée sur la minimisation de la distance entre les observations et les résultats du modèle numérique, sur une période couvrant les heures ou jours précédents. L’état estimé du système est ajusté jusqu’à ce que la simulation numérique parvienne à approcher au mieux les observations passées (de quelques heures à jours dans le passé, jusqu’à l’instant présent). Cet ajustement se fait principalement par la condition initiale du modèle numérique : en changeant la condition initiale, on change la trajectoire correspondante, qu’on cherche à approcher au mieux des données.

Dans la figure ci-dessous, on peut comprendre le principe de l’algorithme 4D-VAR. Les données, ou observations, sont représentées par les petites croix rouges [5]. Elles sont étalées dans le temps (axe des abscisses). De l’autre côté, nous avons le modèle. Si on fournit au modèle une condition initiale (que l’on appelle $X_0$), c’est-à-dire un point de départ à l’instant $t=0$ (à gauche, sur l’axe des abscisses), alors le modèle nous permet d’en déduire la trajectoire correspondante, représentée par la courbe bleue.

Principe de l’assimilation variationnelle

Le but du jeu est d’identifier la courbe bleue : il s’agit de la meilleure solution du modèle. Pourquoi est-elle optimale ? Parce qu’elle commence par un point de départ $X_0$ proche de $X_b$, qu’on appelle l’ébauche - nous y reviendrons un peu plus tard, pas d’inquiétude ! Et parce qu’elle passe près de toutes les observations (croix rouges), L’ébauche, c’est une estimation de ce que la condition initiale devrait être. Elle a plusieurs rôles : d’abord, cela permet de plus facilement identifier la trajectoire qui passe le plus près de toutes les observations, car sans idée du point de départ, cela peut être difficile. Ensuite, elle est généralement le résultat d’une prévision passée, et elle fournit une bonne idée de ce que la solution devrait être. Imaginez les prévisions météorologiques : si aujourd’hui, elle vous disent qu’il devrait faire beau demain, alors demain, les météorologistes utiliseront cette information (issue d’une prévision de la veille) pour faire la nouvelle prévision. Car même si la météo peut changer très vite, elle n’en reste pas moins corrélée avec ce qui s’est passé la veille.

Revenons au 4D-VAR. L’identification de cette trajectoire optimale se fait grâce à une approche « moindres carrés ». Le but du jeu est de rendre la plus petite possible la quantité suivante : la distance au carré entre le point de départ et l’ébauche, plus toutes les distances au carré entre les observations et la solution du modèle (ce sont les petits traits rouges verticaux entre les données et la trajectoire bleue).

La difficulté réside dans l’ajustement entre la trajectoire du modèle et les observations, rendu particulièrement complexe en raison (notamment) des points suivants :

• caractère instable du système : une petite variation de l’état estimé du système peut conduire à une simulation très différente ;
• dimension du problème : comme vu précédemment, l’état à estimer comprend de plusieurs dizaines de millions à plusieurs milliards de valeurs ;
• incertitudes : ni le modèle, ni les observations, ne sont corrects, et la présence de nombreuses erreurs dans les données, et incertitudes dans les modèles utilisés, accroît la complexité du problème [6].

Sur la figure suivante, on peut voir un exemple de calcul, sur un « petit modèle » (modèle simplifié, qui est assez facile à utiliser, mais qui simule assez bien les écoulements géophysiques). Imaginez une très grande piscine, de 2000 km de côté (ce serait sans doute plus simple de penser à l’océan Atlantique !). La hauteur de l’eau n’est pas la même partout (à cause des courants marins, notamment), et les zones en rouge représentent une eau plus élevée que la moyenne, les zones en bleu correspondent à une eau moins élevée que la moyenne, la moyenne étant en vert.

La figure de gauche montre l’ébauche du système, il s’agit donc d’une estimation de ce que la solution devrait être. La figure du milieu représente l’état réel de la piscine, que l’on cherche à retrouver grâce aux observations. Et la figure de droite représente la solution optimale fournie par le 4D-VAR. Pour cette simulation assez simple, nous avons découpé la piscine suivant une grille avec environ 200000 points, et le modèle utilise 3 variables physiques qui sont la hauteur de l’eau (qui est représentée ici), et les deux composantes de la vitesse horizontale de l’eau. Nous avons utilisé uniquement des observations portant sur la hauteur de l’eau (aucune donnée sur la vitesse), et réparties de façon assez éparse en temps et en espace : 1 mesure tous les 25 points de grille en espace, et tous les 50 pas de temps. Cela correspond à environ 10000 données chaque jour, à comparer avec les 200000 valeurs à identifier - 1 à chaque point de la grille (ou plutôt 3, puisqu’il y a 3 variables physiques).

Comme on peut le voir, la solution identifiée est très proche de la solution exacte, alors que sans assimilation, nous aurions du nous contenter de l’ébauche (à gauche). Celle-ci paraît très éloignée de la solution réelle, mais elle ne l’est pas autant qu’on le croit, car même si les hauteurs d’eau ne sont pas correctes, les principaux courants (notamment celui qui part du milieu du bord gauche en direction du bord droit, entre le bleu et le rouge) est assez bien placé.

Filtre de Kalman

Une autre grande classe d’algorithmes d’assimilation de données repose sur le filtre de Kalman (qui, comme son nom l’indique, a été développé par Rudolf Kalman, dans les années 60). Cette méthode est généralement qualifiée de séquentielle, par opposition à l’approche variationnelle précédente. En effet, l’idée est de traiter les observations ’séquentiellement’ dans le temps, c’est-à-dire au fur et à mesure qu’elles sont disponibles. L’état estimé du système est alors corrigé ’en temps réel’, avec les observations. En l’absence de données disponibles, le modèle numérique d’évolution du système est utilisé pour progresser dans le temps, et lorsque des données arrivent, le système les ’filtre’ et corrige l’état estimé du système, dans une sorte de compromis entre le précédent état et les observations.

Principe du filtre de Kalman

La figure ci-dessus présente le principe du filtre de Kalman. Partons de la gauche : la courbe en noire qui arrive au point $x^f$ représente le résultat d’une prévision passée. A l’instant matérialisé par $k-2$, des observations $y^o$ sont disponibles (triangle vert). On « filtre » alors la prévision $x^f$ en utilisant les données $y^o$, pour obtenir une moyenne (ou un compromis) $x^a$. L’état a donc été corrigé grâce aux observations. Puis on utilise le modèle pour repartir de ce point $x^a$, et avancer dans le temps, jusqu’à l’instant $k-1$. De nouvelles données $y^o$ sont disponibles (triangle bleu), et on corrige alors la solution $x^f$ pour obtenir le meilleur compromis possible, $x^a$, entre la trajectoire que l’on vient de calculer, et les données disponibles. Et ainsi de suite !

Cette méthode présente l’avantage d’une mise en oeuvre a priori plus simple, mais une difficulté réside dans la gestion des erreurs. En effet, les données et le modèle étant ’faux’ (comme évoqué précédemment, avec de nombreuses incertitudes), le compromis optimal entre les données nouvellement disponibles et l’état estimé du système (afin d’obtenir le nouvel état estimé) n’est optimal que si l’on dispose d’informations suffisantes sur les incertitudes. La propagation de ces incertitudes, au même titre que celle de l’état estimé du système, est responsable d’un coût de calcul extrêmement élevé.

Nudging

Il existe des méthodes plus simples à mettre en oeuvre pour identifier la condition initiale d’un système.

Le terme anglais ’nudging’ signifie au sens propre ’donner un coup de coude’ et au sens figuré ’encourager’ ou ’pousser’. Pour pousser les résultats du modèle vers les résultats des observations (mesures de vent, de température, d’humidité, de pression, de salinité, ...), les physiciens (météorologues ou océanographes) ont imaginé d’ajouter aux équations de leur modèle un terme de rappel entre les observations et les quantités correspondantes venant des calculs à partir du modèle (équations des fluides géophysiques, conservation de la matière et de l’énergie, lois de comportement, ...). Ainsi ce terme dit de ’relaxation newtonienne’ va permettre de ’nudger’ le modèle vers les observations et donc de satisfaire un compromis entre modèle et observations. On peut démontrer mathématiquement que cet observateur (dit de Luenberger) est, sous certaines conditions, asymptotique et donc que, quand le temps va tendre vers l’infini, l’état du modèle converge vers la réalité expérimentale.

La figure ci-dessous montre le principe du nudging. Les observations sont représentées par les croix rouges. La courbe verte correspond à une solution du modèle, pour laquelle on a utilisé un point de départ différent des observations. La trajectoire (courbe verte) reste donc éloignée des observations. On voit que même après un temps assez long ($t=10$), la solution reste éloignée des données. En rajoutant un terme de rappel vers les observations (terme de nudging), on obtient alors la courbe bleue, qui se rapproche lentement - mais sûrement !! - des observations. On voit toutefois qu’il faut attendre un temps relativement long pour se rapprocher suffisamment des données.

Mais dans la réalité on ne peut pas se permettre d’attendre un temps long pour avoir cette concordance.

C’est pour cela qu’on a mis au point le ’back and forth nudging’, qui consiste à itérer des résolutions directes et rétrogrades du modèle avec un terme de nudging. On tombe alors sur la notion de réversibilité des équations, c’est-à-dire la possibilité de les résoudre aussi bien en avançant qu’en remontant le sens du temps. Ce dernier point est bien sûr indispensable pour remonter à la condition initiale. Les équations des ondes sont réversibles, par contre celles modélisant la diffusion de la chaleur ne le sont pas. Les équations des fluides géophysiques ne sont pas totalement réversibles, mais, ô joie, le terme de nudging permet de stabiliser la résolution rétrograde de ces équations (celle qui remonte le temps).

En pratique, on part d’une certaine condition initiale — par exemple l’ébauche, comme dans le 4D-VAR —, on résout les équations du modèle avec le terme de nudging, puis on effectue une résolution rétrograde (de la fin vers le début de la fenêtre temporelle), avec le nudging, en partant de la solution obtenue à la fin de la résolution directe. Et on obtient ainsi à la fin de la résolution rétrograde une nouvelle estimation de la condition initiale. On itère ce processus jusqu’à convergence, ce qui dans la pratique (numérique) est obtenu en quelques itérations [7]. Cette méthode présente l’avantage d’être très simple à mettre en oeuvre et d’être également très rapide.

La figure ci-dessus reprend la simulation précédente : les observations (croix rouges), la trajectoire sans nudging (courbe verte) et la trajectoire avec nudging (courbe bleue) sont les mêmes que dans l’exemple précédent. On a juste zoomé sur un temps un peu plus court, afin de mieux voir ce qui se passe au début.

Si on ne peut travailler que sur un temps assez court (nous avons choisi $t=2$ ici), on voit que la solution avec nudging (courbe bleue) est encore assez éloignée des observations. On applique alors le nudging direct et rétrograde, et on obtient la courbe noire. Au début, elle est similaire à la courbe bleue, c’est du nudging ’direct’ (classique). Mais une fois qu’on arrive à l’instant $t=2$, on voit bien que la solution est encore très éloignée des observations, alors on repart en sens inverse, autrement dit, on remonte le temps !! On retourne de $t=2$ à $t=0$, toujours avec un terme de rappel vers les observations. On voit alors que la courbe en noire se rapproche encore un peu des observations quand on revient vers l’instant initial. Une fois arrivé en $t=0$, on repart dans le sens « normal », et ainsi de suite. Après 3 allers-retours, on constate que la solution s’est fortement rapproché des observations, et cela permet d’obtenir une solution très proche des observations. Et même sans terme de nudging après $t=2$, la solution en noir est (et reste) plus proche des observations que la solution du nudging classique (en bleu). On voit bien qu’on a contourné la difficulté induite par un temps « court » : en faisant plusieurs allers-retours dans la petite fenêtre de temps $[0;2]$, on s’est redonné du temps, afin de mieux se rapprocher des observations !

Cet algorithme a fait ses preuves sur des modèles allant des équations de Lorenz à un vrai modèle complet d’océan.

Le film ci-dessous montre un exemple avec 5 allers-retours de nudging direct et rétrograde sur notre grande piscine (souvenez-vous, 2000 km de côté !!), pendant lesquels la trajectoire que l’on construit se rapproche des observations, permettant ainsi d’obtenir un état plus précis de cet océan, et ainsi de réaliser de bien meilleures prévisions !

Regardez bien dans quel sens les tourbillons tournent, le modèle fonctionne alternativement dans un régime « normal » (direct) et un mode « rétrograde » (les tourbillons vont dans l’autre sens). Le modèle change de régime toutes les 5 secondes dans le film (5 secondes en direct, puis 5 secondes en rétrograde, etc), et nous avons représenté 5 allers-retours (pendant les 50 premières secondes) avant d’utiliser le modèle en mode direct, sans nudging, afin de simuler des prévisions (pendant les 20 dernières secondes du film).

Conclusion

Les prévisions météorologiques et océanographiques représentent un défi encore difficile, et d’actualité, même si les progrès ces dernières années ont été phénoménaux, que ce soit en termes de moyens de calcul, d’observations disponibles (grâce aux satellites) ou de méthodologie. Mais il reste encore beaucoup de travail, et ce problème est au coeur des préoccupations de nombreux instituts de recherche, universitaires ou privés, et d’organismes de prévision météorologique dans le monde [8].