Back | Next


Chapitre 13 Un modèle global de haute précision : FES2000

13.1 Présentation du chapitre

Le chapitre précédent nous a permis de dégager les principales améliorations à apporter à notre modèle global afin d’obtenir des solutions de marée en s’affranchissant de toute utilisation de données de terrain. Nous allons maintenant appliquer toutes ces améliorations à une modélisation de la marée océaniques à l’échelle globale. En créant une nouvelle bathymétrie, en générant de nouveaux maillages éléments finis, en utilisant de nouveaux effets de charge et au moyen de la résolution par blocs introduite précédemment nous sommes à même de calculer de nouvelles solutions de marées libres à haute résolution. Nous expliquons dans ce chapitre comment nous avons obtenu ces nouvelles solutions dénommées FES2000.

13.2 Une nouvelle bathymétrie

La bathymétrie est un des paramètres d’entrée essentiels de notre modèle hydrodynamique pour améliorer le calcul d’un nouveau spectre de marée. C’est pourquoi, nous avons mis en œuvre le calcul d’une nouvelle bathymétrie mondiale adaptée à l’usage spécifique des modèles de marée FES.

13.2.1 Travaux mis en œuvre

Le calcul d’une nouvelle bathymétrie pour nos besoins spécifiques s’est décomposé en plusieurs étapes :

13.2.2 Choix des bases de données

Les besoins spécifiques de la formulation de notre modèle hydrodynamique de calcul de la marée océanique entraînent plusieurs nécessités sur la bathymétrie :

§ bonne bathymétrie en zones pélagiques afin de mieux modéliser les échanges océaniques de grandes longueurs d’ondes ;

§ bathymétrie très précise sur les plateaux continentaux et en petits fonds car de nombreuses ondes de marée y prennent naissance (ondes non linéaires) ;

§ bons positionnement des forts gradients de bathymétrie afin de mieux tenir compte de la génération des ondes internes ;

§ bonne bathymétrie en zones polaires car peu de bathymétries de ces zones sont disponibles et validées ;

§ bons positionnement des côtes car étant très précis, ils permettent d’ajuster les petits fonds au mieux.

C’est pourquoi la compilation de diverses bases de données est nécessaire.

13.2.3 Bathymétries mondiales

A ce jour, plusieurs bathymétries mondiales du fond des océans sont disponibles. Nous en avons sélectionné trois, basées sur des méthodes de calcul différentes : Elles ont chacune leurs caractéristiques propres, leurs qualités et leurs défauts.

13.2.4 ETOPO5

ETOPO5 (Earth Topography - 5 Minute) est une base de données digitales des élévations des fonds océaniques et des terres émergées de la Terre. Elle a été composée à partir de plusieurs bases de données grillées et digitalisées puis projetée sur une grille espacée de 5 minutes en latitude et 5 minutes en longitude. La bathymétrie océanique a été compilée par l’U.S. Naval Oceanographic Office, et révisée en 1987. Les topographies terrestre et océanique relevant des Etats Unis, du Japon et de l’Europe de l’Ouest ont été assemblées à partir de bases de données grillées fournies par la Defense Mapping Agency. Celles de la Nouvelle Zélande proviennent du Department of Scientific and Industrial Research of New Zealand; celles de l’Australie ont été fournies par le Bureau of Mineral Resources of Australia. Ces différentes topographies des zones continentales et océaniques ont été interpolées sur une grille de 10-minute par l’U.S. Navy Fleet Numerical. ETOPO5 a une couverture mondiale aussi bien pour les fonds océaniques que pour les terres émergées. L’édition révisée de 1988 contient de légères modifications afin de mieux prendre en considération les forts gradients de bathymétrie et de topographie. Quelques points de sondes ont été retirés afin d’enlever des erreurs d’interpolation.

13.2.5 DBDB5

DBDB5 (Digital Bathymetric Data Base - 5 minute grid) est composée de données bathymétriques digitalisées des fonds océaniques interpolée sur une grille de 5 minutes en latitudes par 5 minutes en longitude. Les données ont été assemblées par le Naval Ocean Research and Development Agency (NORDA) et le Naval Oceanographic Office (USNOO). Une version précédente était disponible sous le nom 'SYNBAPS II'. La base de données révisée contient une version retravaillée de la bathymétrie de l’hémisphère nord et inclut maintenant l’océan Arctique. Quelques légères modifications ont été apportées dans l’hémisphère sud. Les données ont été calculées par interpolation entre contours digitalisés (isobathes) à des intervalles réguliers de 200m en profondeur. Ces isobathes ont été digitalisées à partir de cartes marines (échelle : 1 pouce par degré en longitude soit des cartes aux 1 / 4.000.000) de l’US Navy datant d’avant 1983. La grille DBDB5 à 5 minutes contient environ 8 millions de points sondes. Les profondeurs ont été calculées en ayant considéré une vitesse de propagation du son dans l’eau de 1500 m/s.

13.2.6 Smith&Sandwell

Smith et Sandwell ont combiné des données d’échos sondeurs de navires avec des données de gravité dérivant de missions satellitaires afin de produire une estimation grillée des fonds océaniques pour toutes les surfaces océaniques dépourvues de glaces flottantes (entre +/-72 degrés de latitude). Les profondeurs proviennent de l’extraction de 6905 campagnes océanographiques compilées par le NGDC (Trkdas CD version 3.2), de la base de données de Scripps et Lamont, et de diverses autres données. La qualité des données à été contrôlé conformément à la méthode développée par Smith [1993]. Les champs de gravité satellitaires combinent toutes les données des satellites ERS1 et Geosat dont les données déclassifiées en 1995. Ces données ont une RMS de 3-5 mGal et une résolution de 20 km de longueur d’onde. La solution digitalisée de bathymétrie est projetée sur une grille Mercator espacée de 2 minutes en latitude et 2 minutes en longitude à l’équateur (grille non régulière). La méthode est similaire à celle développée par les travaux de Smith [1994]. Les données des campagnes océanographiques ont permis de calibrer la fonction de transfert de la gravité vers la topographie et des fournir les ondes longues périodes qui font partie intégrante de la solution. (les compensations isostatiques de la topographie limitent l’estimation de la gravité pour des longueurs d’ondes inférieures à 160 km). La corrélation des profondeurs estimées et mesurées est forte pour des longueurs d’ondes supérieures à 20 km, et les différences entre elles sont généralement inférieures à 100 m. Cependant des différences beaucoup plus grandes peuvent apparaître dans des zones à forts gradients de topographie et les profondeurs minimales aux sommets marins ne sont pas très précises. La dernière étape de calcul de la topographie a permis de lisser les données afin qu’elles soient en meilleur accord avec les échos sondeurs des navires et les fichiers de côtes digitalisées haute résolution de GMT. La carte de topographie calculée est légèrement différente de celles calculées au moyen de champs de gravité satellitaires. En effets les couches sédimentaires sont mal modélisées et certaines fractures ou bords de plateaux peuvent ne pas apparaître dans certaines zones.

13.2.7 Base de données de côtes digitalisées

En 1989, la US Defense Mapping Agency (DMA) a réalisé la World Vector Shoreline [WVS, 1989]. La WVS est une base donnée digitalisée à l’échelle 1:250,000 référencée par rapport au World Geodetic System (WGS-84). La couverture globale a été terminée en juillet 1989, avec comme spécification que 90% des côtes devaient être localisées à 500 mètres (i.e. 2mm au 1:250,000) de leur véritable position géographique. Des imprécisions demeurent toutefois pour 10% des côtes du monde, majoritairement en Antarctique. Cependant, une précision de 500 mètres est largement supérieure à la taille des mailles des maillages éléments finis de FES, ce qui en fait une base de données très précise et permis de positionner au mieux les zones en très petits fonds près des continents et de faire apparaître ou de mieux positionner de petites îles en plein océan.

13.2.8 Base de données d’isobathes

Le positionnement des marches de plateaux continentaux est essentiel pour tenir compte des propagations des ondes de marées. L’utilisation d’isobathes est donc une information supplémentaire et utile pour le calcul d’une meilleure bathymétrie. La base de données digitalisées GEBCO [GEBCO97, 1997] fournit de nombreuses isobathes échantillonnées à diverses profondeurs pour tous les océans du globe. La première version de la General Bathymetric Chart of the Oceans (GEBCO), qui est parue en 1905, établissait une tradition pour publier une série de cartes à couverture mondiale à l’échelle 1:10 million. Cette approche s’est poursuivie jusqu’à la publication de la Cinquième Edition de GEBCO par Canadian Hydrographic Service entre 1978 et 1982. En 1983, le Joint IOC/IHO Guiding Committee pour GEBCO a décidé que la Cinquième Edition de GEBCO devait être digitalisée Durant plus de dix ans, le Bureau Gravimétrique International à Toulouse et le British Oceanographic Data Centre (BODC) ont produit une base de données couvrant tous les contours bathymétriques (isobathes), les côtes et les sondes de campagnes océanographiques sous forme du GEBCO Digital Atlas (GDA). La Cinquième Edition GEBCO est composée de 16 cartes en projection Mercator couvrant le monde de 72°N à 72°S au 1:10 million à l’équateur et de deux cartes polaires en projection stéréographique polaire de 64°N à 90°N et de 64°S 90°S au 1:6 million à 75° de latitude. La base de données est régulièrement remise à jour sur des zones dépendant de la densité des échos sondeurs des navires. Ces zones sont à des échelles comprises entre 1:10 million et 1 250,000. Les contours fournis sont au minimum à 200m et à 500m de profondeur puis par intervalle de 500m. Bien souvent d’autres isobathes sont disponibles, accroissant ainsi la précision de la bathymétrie en petits fonds en particulier. Les profondeurs sont données en mètres corrigés.

13.2.9 Autres bases de données

Des bases de données bathymétriques locales ont aussi été utilisées. De nombreuses études scientifiques ou d’ingénierie sur certaines zones du globe ont nécessité une étude locale précise de détermination de la bathymétrie. Certaines bathymétries nous ont été fournies par des scientifiques, d’autres par des organismes publics, d’autres enfin sont issues de digitalisations de cartes réalisées par nos soins. Ces bases de données sont souvent issues de besoins locaux, ce qui leur confèrent une précision accrue.

Diverses zones ont ainsi été obtenues :

Les points de sondes que ces bathymétries locales ont ajoutés aux autres données déjà existantes ont fourni une importante information supplémentaire. En particulier dans les zones de petits fonds, elles ont permis d’améliorer considérablement les bases de données citées précédemment.

13.3 Algorithme de calcul

13.3.1 Objectifs de l’algorithme

Le but de l’algorithme est de calculer une nouvelle bathymétrie mondiale grillée de pas de grille 2 minutes en longitude et 2 minutes en latitude (afin de garder la précision de la bathymétrie Smith&Sandwell en particulier en plein océan).

13.3.2 Interpolation par éléments finis

Les différentes bases de données citées précédemment fournissent des points de sondes, c’est-à-dire un point de la surface du globe en latitude et en longitude avec soit la profondeur océanique, soit l’élévation terrestre. Si ce ne sont les données issues de grilles, les points de sondes ont des coordonnées exactes. Ce ne sont pas des interpolations. Afin de tenir compte au mieux de la qualité de ces données, un algorithme d’interpolation linéaire basé sur les éléments finis a été mis en place.

Une zone de calcul est choisie. Elle correspond en général aux points de sondes disponibles et provenant d’une même origine. Le calcul de la bathymétrie a été divisée en 35 zones (Figure 79) qui prennent en compte à la fois les zones océaniques et les possibilités de calcul des ordinateurs.


Figure 79 : Découpage en 35 zones pour le calcul de la nouvelle bathymétrie

Chaque zone a été subdivisée en petits domaines de 5°x5° (voire plus petits encore si le nombre de points de sondes à interpoler était très important) afin d’être en accord avec la puissance de calcul informatique disponible.

13.3.3 Exemple : calcul d’une bathymétrie sur un domaine de 5°x5

°Pour illustrer la méthode d’interpolation, un domaine de 5°x5° au nord de la Corse est utilisée.

Sur cette zone sont extraits tous les points de sondes de la base de données WVS qui correspondent à des profondeurs nulles. Ces points, très nombreux sont interpolés sur une grille de 0,5 x0,5° (supérieure à la précision de la base de données mais inférieure au pas de grille souhaité pour la bathymétrie mondiale).


Figure 80 : Points de sondes côtiers

Puis sont extraits tous les points de sondes disponibles dans toutes les autres bases de données disponibles sauf ETOPO5, DBDB5 et Smith&Sandwell. Ces points sont eux aussi interpolés sur une grille 0,5’x0,5’.


Figure 81 : Points de sondes côtiers, des isobathes et des bases de données personnelles

Ensuite le contour de la zone est déterminé pour compléter les points de sondes par des bathymétries mondiales.


Figure 82 : Détermination du contour de la zone pour l’interpolation ultérieure

DBDB5 est de bonne qualité le long des côtes, Smith&Sandwell est de bonne qualité en plein océan, donc pour compléter l’information des points de sondes définis précédemment, DBDB5 est utilisée de 0 à 200m, une interpolation linéaire entre DBDB5 et Smith&Sandwell est calculée entre 200 et 500m et Smith&Sandwell complète les fonds en dessous de 500m.


Figure 83 : Ensemble des points de sondes utilisés pour l’interpolation

Pour toutes les zones ou Smith&Sandwell n’est pas définie, ETOPO5 est utilisée. Pour les fonds océaniques, les points ajoutés sont interpolés sur une grille de 2’x2’. Pour les terres, une interpolation de 4’x4’ est suffisante. En effet après l’interpolation, toutes les élévations des terres émergées seront recopiées de la base de données Swith&Sandwell interpolée à 2’x2’ et qui est très précise. Ensuite, afin de réaliser l’interpolation linéaire de tous ces points de sondes, un maillage éléments finis est calculé. L’utilisation des éléments finis est très avantageuse dans notre cas car elle nous permet de tenir compte de la position exacte des points de sondes. En effet, si nous avions travaillé en différences finies nous aurions deux possibilités. Il aurait fallu prendre un pas de grille tellement petit pour ne pas perdre la précision de positionnement de ces points de sondes que nous aurions largement dépassé la puissance de calcul des ordinateurs à notre disposition (en l’occurrence deux stations Sun UltraSparc 1 cadencée à 177MHz avec 256 Mo de mémoire vive). Ou bien nous aurions du interpoler la position des points de sondes sur des points espacés d’une grille différences finies et nous aurions perdu l’information exacte apportée par chaque point de sonde.


Figure 84 : Maillage des points de sondes pour l’interpolation

Enfin, de ce maillage, une nouvelle bathymétrie a été interpolée linéairement Elle tient compte de la position réelle de tous les points de sondes, sur une grille 2’x2’.


Figure 85 : Nouvelle bathymétrie interpolée (profondeurs et altitudes en mètres)

Ainsi toutes les zones de 5°x5° sont calculées indépendamment puis elles sont assemblées en une bathymétrie mondiale. La continuité aux frontières des différentes zones est assurée car elles sont absolument identiques pour toutes les zones jointives.

Dans les zones polaires, Smith&Sandwell n’étant pas définie, ETOPO5 a été utilisée. Afin d’assurer une continuité entre les deux bathymétrie mondiales, une interpolation linéaire a été calculée entre –68° (Nord et Sud) et 72° (Nord et Sud).

Les maillages aux éléments finis des solutions FES antérieures étaient contraints par un critère basé sur la longueur d’onde théorique des ondes de marée de fréquence semi-diurne. Ils ne tenaient malheureusement pas compte de l’existence sur les pentes de la topographie d’ondes à plus courte échelle. Ils étaient également mal adaptés, sur ces même pentes, à la détermination précise des vitesses barotropes. Ils étaient sous-résolus pour les ondes quart-diurne. Du fait de la connaissance a priori que nous avons maintenant de la marée dans les bandes longues périodes, diurne, semi-diurne et bientôt quart-diurne, nous avons les moyens de corriger ce sous échantillonnage spatial en remaillant les domaines de modélisation. Par ailleurs, nous profiterons de cette remise à plat pour compléter notre modèle avec les mers adjacentes non incluses jusqu’à présent (Mer rouge, Golf Persique, Mer Noire...).

Nous présentons dans la suite la construction et les caractéristiques adaptées à l’usage spécifique des modèles de marée FES.

13.4 Création des maillages éléments finis

13.4.1 Choix d’une stratégie de construction

13.4.1.1 Définition de la stratégie des sous domaines

Des études énergétiques de dissipation de la marée sur une couverture océanique globale, (par exemple [Le Provost and Lyard, 1997]), ont permis de distinguer des zones de très forte dissipation :

§ Le plateau Européen ;

§ Le plateau de Patagonie ;

§ La baie d’Hudson ;

§ La Mer du Labrador

§ La zone de la Malaisie ;

§ La Mer Jaune et les Mers de Chine ;

Ces régions à forte dissipation sont étroitement corrélées avec les petits fonds océaniques. Elles nécessitent donc une résolution de maillage plus importante afin de mieux déterminer dénivellation et vitesses.

En outre des zones de petits fonds avec des bathymétries approximatives se trouvent dans les zones polaires :

§ L’Océan arctique proche des côtes ;

§ La Mer de Ross ;

§ La Mer de Weddell.

Le rôle fondamental joué par les phénomènes non-linéaires dans la dynamique de la marée (ainsi que le regain d’intérêt pour l’océanographie côtière) et la nécessité de bien représenter ces phénomènes nous imposent d’étudier avec une résolution plus importante que de coutume les zones citées précédemment. L’approche multi bassins de notre modèle s’adapte particulièrement bien à une définition hétérogène de la résolution des maillages en zones peu profondes. Les régions citées précédemment ont été détachées du modèle de bassin dans lequel elles étaient incluses auparavant, pour constituer un modèle indépendant qui sera au final traité de façon analogue aux modèles de bassins océaniques dans l’algorithme de résolution globale.

13.4.1.2 Choix des sous domaines

En considérant les remarques précédentes, le maillage global a été divisé en 12 domaines : Grâce à cette découpe de l’océan global, il est possible de reconstruire indépendamment un domaine éléments finis. Les maillages des zones à forte dissipation par exemple peuvent être retaillés sans que les maillages de plein océan ne soient modifiés.

13.4.2 Critères sur les éléments triangulaires des maillages

13.4.2.1 Choix des éléments triangulaires

Pour les calculs des solutions FES à l’échelle globale, les équations sont résolues sur un maillage global éléments finis. Chaque élément triangulaire du maillage est Lagrange P2, c’est-à-dire que dénivellation sont calculées pour les 3 sommets des triangles ainsi que pour les 3 milieux des côtés d’un triangle, soit 6 nœuds en tout. Les vitesses sont calculées aux 7 points de Gauss (points rouges) ce qui permet de calculer une interpolation des vitesses et de les dériver pour en déduire la dénivellation aux nœuds. L’ancien maillage global comprenait plus de 300000 nœuds Lagrange P2.

Afin de mieux modéliser la marée, le nouveau maillage est fortement raffiné suivant les trois critères présentés ci-dessous.

13.4.2.2 Précision le long des côtes

Les côtes constituent les frontières fermées du domaine de modélisation. Les maillages sont donc ancrés à ces limites. Le choix des fichiers de côtes et de leur discrétisation est donc primordial. Afin de ne pas outrepasser la puissance de calcul des super calculateurs disponibles, une distance de 15 km minimum de côté pour un triangle a été choisie. Les fichiers de côtes utilisés proviennent d’une base de donnée digitalisée, de haute précision. Ainsi entre deux nœuds Lagrange P2 la distance est de 7,5 km environ le long des côtes.

13.4.2.3 Critère sur la longueur d’onde

Le raffinement des maillages respecte le critère local définit par Le Provost et Vincent [1986]. Celui-ci fixe le maximum de distance  pour le côté d’un triangle afin d’obtenir une simulation correcte de la propagation d’une onde de marée à une valeur inférieure à :

(13.138)

13.4.2.4 Critère sur le gradient topographique

Cependant ce critère peut ne pas être suffisant dans des zones de forts gradients topographiques [Luettich and Westerink, 1995]. C’est pourquoi un autre critère de génération des maillages à été choisis. Il tient compte des gradients topographiques de la topographie. Ainsi la distance maximale pour le côté d’un élément triangulaire ne doit pas dépasser la valeur :

(13.139)

La bathymétrie utilisée dans le calcul de ce critère est issue de ETOPO5 et a été améliorée en petits fonds.

13.4.2.5 Choix du critère

Pour la génération des maillages, le plus contraignant de deux derniers critères a été sélectionné pour la création d’un élément triangulaire.

13.4.3 Algorithme de calcul des maillages éléments finis

13.4.3.1 Algorithme de calcul

Afin de construire un maillage élément finis, il faut définir une densité de taille de maille pour un élément triangulaire sur une zone bien définie. Cette zone est délimitée par les frontières fermées et la densité respecte les critères établis ci-dessus. De cette densité de probabilité est construit le maillage éléments finis par l’algorithme de Ruppert [1995] basé sur la méthode de Delaunay [1934].

La triangulation de Delaunay d’un ensemble de points est la triangulation de cet ensemble qui respecte la propriété stipulant qu’aucun point dans cet ensemble ne se trouve dans le cercle extérieur de tous les triangles de la triangulation. Ce cercle extérieur est le cercle qui passe par les trois sommets d’un élément triangulaire.

L’algorithme de Ruppert est peut être le premier algorithme de génération de maillages éléments finis qui fonctionne en pratique. Il génère des maillages sans angles trop petits en utilisant relativement peu de triangles (puisque la densité de triangle peut-être ajusté par l’utilisateur) et permet à la densité de triangle de varier rapidement sur de courtes distances.

13.4.3.2 Logiciel de génération des maillages

Utilisant l’algorithme de Ruppert, le logiciel Triangle [Shewchuk, 1996], génère un maillage éléments finis à partir des densités de tailles de mailles calculées d’après les trois critères définis précédemment.

13.4.4 Caractéristiques des maillages

Chaque maillage à ses caractéristiques propres. Ils ont été dimensionnés afin de répondre à la puissance de calculs des derniers super calculateurs accessibles. Les caractéristiques des maillages sont données Figure 86 et dans le Tableau 37.


Figure 86 : Caractéristiques des maillages

La Figure 86 indique le nombre d’éléments triangulaires, le nombre de nœuds Lagrange P2 et la largeur de bande qui composent chaque maillage éléments finis.
 

Nom
Triangles
Nœuds P2
Largeur bande
Frontières ouvertes P2
Taille tmpdir (bytes)
Arctic
37 086
79 302
1 061
120
4 037 423 424
Atlantic
78 385
161 615
1 407
821
10 973 761 728
Europe
10 113
21 992
254
145
267 774 592
Indian
73 554
151 181
1 635
379
11 999 666 800
Indonesia
66 982
140 428
1 252
461
8 415 613 760
Labrador
16 530
35 787
520
177
892 670 928
Mediter
22 797
49 505
676
25
1 605 546 160
PacificN
64 720
133 768
1 682
692
11 060 677 632
PacificS
67 374
137 672
1 355
833
8 951 984 128
Patagonia
11 045
23 424
488
267
548 308 992
Ross
3 190
6 725
218
61
70 262 800
Weddel
7 216
15 165
323
81
234 875 520
TOTAL
458 992
956 564
-
4 062

Tableau 37 : Caractéristiques des 12 maillages couvrant l’océan global

La taille du tmpdir est la taille du fichier qui doit être écrit sur disque dur au cours de la factorisation d’une matrice intermédiaire dans le code CEFMO calculant les solutions FES. Sa taille étant conséquente, elle conditionne les ressources nécessaires des machines de calculs. Elle est égale à :

(13.140)

avec :

Les différents éléments que nous avons développé précédemment (nouvelle bathymétrie, nouveaux critères de taille de maille et nouvel algorithme de calcul des maillages) nous ont permis de construire de nouveaux maillages éléments finis. En s’appuyant sur trois critères de tailles de mailles qui prennent en compte la bathymétrie des fonds océaniques, une densité de tailles de mailles des éléments triangulaire à permis de les créer grâce à un algorithme de calcul basé sur la méthode de Delaunay. 12 maillages ont été crées. Assemblés, ils forment un maillage triangulaire composé de presque un million de nœuds P2. La distance entre deux nœuds P2 est de l’ordre de 7,5 km en petits fonds et de l’ordre de 75 km en plein océan sur de faibles gradients de bathymétrie. Dans les zones de fort gradient bathymétrique la distance est de l’ordre de 10 km.


Figure 87 : Nouveau maillage éléments finis global

13.5 Calcul des nouvelles solutions libres

13.5.1 Ressources informatiques

Comme nous venons de le voir le nouveau maillage éléments finis est beaucoup plus important que celui utilisé pour calculer les solutions FES94.1, FES95.2, FES98 et FES99. Les ressources informatiques nécessaires sont donc aussi plus importantes. Nous avons alors utilisé le nouveau supercalculateur mis à disposition par l’IDRIS en remplacement des anciennes machines vectorielles Cray : le supercalculateur NEC SX-5.

Nous donnons pour information, les caractéristiques de ce supercalculateur utilisé dans le cadre de notre modélisation hydrodynamique.

13.5.2 Itération sur les deux ondes principales M2 et K1

Comme nous l’avons vu précédemment, notre code de calcul est itératif (Figure 88).


Figure 88 : Processus itératif simplifié du code CEFMO

Avec l’ancien maillage 6 itérations assuraient la convergence. Avec le nouveau maillage, nous avons fait des tests de convergence pour les ondes M2 et K1 afin de déterminer le nombre d’itérations nécessaires à la convergence. Ces solutions ont été calculées en mode forcé, c’est-à-dire que nous avons calculé les 12 bassins océaniques introduits dans le paragraphe 13.4.1.2, indépendamment les uns des autres. Les conditions aux limites sur les frontières fermées sont un flux nul Sur les frontières ouvertes, nous avons imposé des conditions aux limites en dénivellation extraites de la solution FES99. Imposer ces conditions aux limites nous permet de converger beaucoup plus vite vers un état stationnaire au cours du processus itératif. De plus, comme nous l’a montré l’étude locale du chapitre précédent, nous avons calculé nos solutions en mode mixte qui est plus proche de la réalité des phénomènes dissipatifs de la marée. Le mode mixte nous permet de tenir compte tout à la fois de la dominance des ondes diurnes et des ondes semi-diurnes pour la linéarisation du coefficient de frottement (M2 est dominante dans certaines zones de l’océan mondial et K1 dans les autres). En effet, nous avons vu par exemple que les marées à Dumont d’Urville sont diurnes (cf. mesures de ROSAME). Les marées semi-diurnes ne sont pas dominantes partout à la surface de la Terre. Afin de pouvoir quantifier ces itérations, nous avons comparé nos sorties de calculs avec la base de données ST95 composée de 95 données marégraphiques pélagiques (Figure 89 et Figure 90) et la base Topex composée de données altimétriques dans des fonds de plus de 200 m(Figure 91 et Figure 92).


Figure 89 : RMS (en cm) des solutions calculées pour M2 par rapport à ST95 en fonction du nombre d’itérations


Figure 90 : RMS (en cm) des solutions calculées pour K1 par rapport à ST95 en fonction du nombre d’itérations


Figure 91 : RMS (en cm) des solutions calculées pour M2 par rapport à la base Topex en fonction du nombre d’itérations


Figure 92 : RMS (en cm) des solutions calculées pour K1 par rapport à la base Topex en fonction du nombre d’itérations

Ces 4 figures nous montrent que nous avons la convergence assurée des solutions FES2000 en mode mixte forcé pour le nouveau maillage globalement atteinte en 10 itérations. Tous nos calculs en mode forcé avec le nouveau maillage nécessitent donc 10 itérations. Le saut de RMS que nous pouvons constater entre l’itération 0 et l’itération 1 s’explique par le fait que pour initialiser notre processus itératif, nous imposons initialement des vitesses arbitraires de 2 m.s-1 qui sont valables en milieu côtier, mais bien trop fortes en plein océan (plutôt de l’ordre du cm.s-1), il faut donc une à deux itérations supplémentaires pour obtenir un champ de vitesses de nature plus physique aussi bien en petits fonds qu’en plein océan.

13.5.3 Calcul des solutions libres

La résolution par blocs introduite dans la Partie 3, nous permet par linéarité de notre système, de retirer la contribution des conditions aux limites ouvertes que nous avons imposé pour calculer nos solutions. Nous obtenons alors ce que nous avons déjà nommé, les solutions libres en dénivellation. Grâce à la puissance de calcul du supercalculateur mis à notre disposition, nous avons pu calculer ces solutions libres. Le Tableau 38 donne les résultats des comparaisons de ces solutions libres pour M2 et K1 avec ST95 et la base de données Topex.
 
Solution libre
Comparaison avec ST95 (cm)
Comparaison avec la base Topex (cm)
M2
14,4
12,3
K1
1,8
1,8

Tableau 38 : Comparaison des solutions libres avec ST95 et la base de données Topex

Comme nous pouvons le voir, nous sommes loin des précisions des solutions actuelles telles FES99. Sachant que la RMS sur les données de ST95 est de 33,5 cm (11,1 cm) la comparaison de notre nouvelle solution libre M2 (K1) avec ST95 représente un écart aux observations d’environ 33% (16%). Avec la solution assimilée FES99, il est de 4% (10%)

Si nous considérons les données de la base Topex, la RMS sur les données est de 26,1 cm (10,5 cm) et donc la comparaison de la solution libre M2 (K1) avec ces données est d’environ 47% (17%). Avec la solution assimilée FES99, elle est de 3% (9%). La Figure 93 présente la différence complexe entre la solution libre FES2000 et la solution FES99 pour l’onde M2. Globalement en plein océan, nous constatons des différences de l’ordre de plusieurs centimètres. Mais, dans la Mer du Labrador, la Mer de Weddell, la Mer de Tasmanie, le détroit de Béring, sur les plateaux Européens et Patagonien, ces différences atteignent 50 cm. Après une étude sur les phases, nous avons constaté que les différences en phases étaient minimes : elles sont donc principalement dues aux amplitudes de marées. Nous avons tracé Figure 94, la différence en amplitude pour l’onde M2 de la solution FES2000 moins la solution FES99, ce qui nous permet de déterminer les différences relatives de FES2000 par rapport à un modèle de marée en accord avec les mesures in situ. Sur la majeure partie de l’océan mondial, la nouvelle solution libre en dénivellation est plus élevée que la solution validée FES99 (à part le détroit de Madagascar la Mer de Timor et le plateau de Patagonie). Nous pouvons donc soupçonner un manque de dissipation dans notre modèle. Nous verrons plus loin que cette hypothèse est confirmée par un bilan énergétique de la marée.


Figure 93 : Comparaison complexe pour M2 de la solution FES2000 et de la solution FES99 (différences en cm)


Figure 94 : Différence en amplitude pour M2 de la solution FES2000 et de la solution FES99 (différences en cm)

De même que pour l’onde M2, la Figure 95 présente la différence en amplitude de la solution libre FES2000 moins la solution FES99 pour l’onde K1. Les différences sont moins marquées que pour M2. Si, majoritairement, l’amplitude de FES2000 est plus élevée dans le Pacifique et une grande partie de l’Indien, ce n’est pas le cas dans l’Atlantique Nord. Pour K1, nous n’avons pas comme pour M2 une amplitude uniformément trop forte de FES2000 par rapport à FES99.


Figure 95 : Différence en amplitude pour K1 de la solution FES2000 et de la solution FES99 (différences en cm)

Nous constatons le même phénomène que pour les solutions libres de FES98 et FES99 : la comparaison des nouvelles solutions libres n’est pas en bon accord avec les données de terrain et les modèles validés de marées. Nous pouvons soupçonner notre modèle de faire des approximations trop restrictives quant à la représentation physique du phénomène des marées. Par exemple, CEFMO ne modélise que les ondes de marées barotropes et ne tient pas compte des ondes baroclines. Il manque donc au moins une paramétrisation dans notre résolution. Le processus itératif présenté Figure 88 montre bien que le calcul des dénivellations de marée est issu du calcul des vitesses des marées et donc des coefficients de frottement que nous avons sur le fond des régions océaniques. Il est donc légitime de penser que si nous avons un problème de paramétrisation dans la résolution de nos équations, il faut regarder de près les vitesses et plus particulièrement le frottement qui est la base même de la physique de la marée océanique. Nous allons tenter d’identifier ce manque par une étude du bilan énergétique des marées issues de nos nouvelles solutions libres FES dans le paragraphe suivant. En effet, nous nous permettons de le rappeler, le but ultime des études en cours au LEGOS sur ces nouvelles solution FES est d’obtenir une solution de marée sans introduire d’information de mesures de terrain, donc sans condition aux limites ouvertes et sans assimilation de données.

13.6 Bilan énergétique

13.6.1 Intérêt du bilan énergétique des marées

Plusieurs études sur la dissipation des marées océaniques ont été entreprises [Munk, 1997]. La plupart de ces études ont été réalisées pour l’onde M2 qui est, du point de vue des courants de marées, largement prédominante. Elle est donc bien représentative à elle seule des bilans énergétiques pouvant être fait pour la marée globale. Nous nous intéresserons dans la suite à la seule onde M2. Pour calculer les bilans énergétiques de la marée M2, il faut disposer des champs de dénivellations et de vitesses de cette marée à l’échelle globale. En général ces champs utilisés pour calculer des bilans énergétiques sont issus de modèles hydrodynamiques [Le Provost and Lyard, 1997] de modèles empiriques utilisant l’altimétrie [Church et al., 2000; Egbert, 1997] ou encore de modèles hydrodynamiques assimilant des données altimétriques [Kantha et al., 1995]. En effet, grâce à ces champs et grâce aux équations de la marée ces bilans peuvent être calculés. Si ces champs sont bons physiquement et si les équations sont exactes, nous devrions avoir un équilibre entre les flux d’énergie que les marées génèrent en déplaçant les masses d’eaux océaniques et l’énergie qui est apportée par les astres générant les marées plus l’énergie que ces mêmes marées dissipent sur le fond océanique. A l’heure actuelle, les estimations pour l’onde M2 sont que la Lune apporte environ 2,5 TW aux masses océaniques [Munk, 1997].

13.6.2 Equations de l’énergie dans le modèle hydrodynamique éléments finis

13.6.2.1 Energie de marée

En reprenant les notations de Le Provost et Lyard [1997], l’énergie  par unité de surface d’une colonne d’eau élémentaire est donnée par :

(13.141)

avec :

Le bilan énergétique est obtenu des équations de la marée [Lyard, 1997] :

(13.142)

avec :

Ainsi le seul phénomène dissipatif qui est considéré ici est dû au frottement de fond. La dissipation causée par la diffusion latérale et par les ondes internes est ignorée. Nous verrons plus loin les conséquences possibles de cette approximation.

En intégrant sur un espace  délimité par ces frontières  et en moyennant sur une période T, nous obtenons :

F = P + D (13.143)

avec :

(13.144)

Due à la périodicité des marées, seuls les termes intégraux pairs sont non-nuls :

(13.145)

Ainsi, le flux d’énergie au travers des frontières ouvertes  du domaine océanique  doit correspondre au travail fourni par le forçage astronomique et l’élévation de la surface des océans due à la marée moins la dissipation.

13.6.2.2 Dissipation de la marée

Dans le modèle CEFMO, le frottement  est paramétrisé au moyen d’un coefficient adimensionnel de type Chézy par la relation :

(13.146)

La forme non-linéaire et non analytique nous pose un problème car nous ne pouvons pas adopter une approche linéaire pour le calcul des bilans énergétiques. En effet le spectre de marée est composé d’un très grand nombre d’ondes périodiques qui interagissent entre elles. Comme nous l’avons vu précédemment, dans le cadre de CEFMO, le coefficient de frottement est quasi linéarisé grâce à la supposition d’une onde dominante (cf. [Le Provost and Lyard, 1997] pour plus de détails). Pour chaque onde d’indice k, le frottement s’écrit sous la forme :

(13.147)

avec :

Nous donnons l’expression de ces coefficients dans la Partie 1.

Cependant Le Provost et Lyard [1997] ont montré que l’onde M2 n’est pas dominante partout (du point de vu des vitesses). L’onde K1 est plus importante dans certaines zones (Mer d’Okhotsk, plateau de Sunda en Malaise et quelques zones autour de l’Antarctique). C’est pourquoi, le coefficient de frottement utilisé pour nos calculs est plus complexe car il prend en compte les effets de M2 et de K1 simultanément. Le frottement est ainsi évalué par :

(13.148)

sachant que :

Le premier (second) terme est prépondérant quand M2 (K1) est dominant.

Ainsi, pour l’onde M2, l’énergie dissipée par le frottement de fond est la somme de deux termes :

D = D1 + D2 (13.149)

13.6.3 Bilan énergétique des solutions libres

Au moyen des équations présentées dans le paragraphe précédent, nous pouvons calculer les différentes énergie mises en jeu dans le processus dynamique des marées océaniques pour l’onde M2 (Tableau 39). Dans le Tableau 39, P1 représente l’énergie relative au forçage astronomique et aux marées terrestres, P2 l’énergie relative aux effets de charges et d’auto-attraction (P=P1+P2), D1 l’énergie relative au frottement de M2 sur les semi-diurnes, D2 l’énergie relative au frottement de M2 sur les diurnes (D=D1+D2).
 
Bassin
Surface 106km2
P1 
P2
P
D1
D2
D
Balance
Arctic
11,863
-3,5
-6,2
-9,7
-132,2
-37,0
-169,2
-178,9
Atlantic
78,551
1439,9
-37,4
1402,5
-366,7
-51,2
-417,9
984,6
Europe
2,22
0,3
-16,2
-15,9
-459,5
-15,6
-475,0
-491,0
Indian
69,239
379,4
-164,8
214,6
-86,6
-28,2
-114,8
99,7
Indonesia
25,923
154,6
-47,4
107,2
-422,9
-256,2
-679,1
-571,9
Labrador
5,46
24,3
-23,6
0,7
-669,7
-58,6
-728,3
-727,6
Mediter
3,074
-0,1
-0,2
-0,2
-2,3
-0,6
-2,9
-3,1
PacificN
82,929
930,2
-69,7
860,5
-227,1
-184,1
-411,3
449,3
PacificS
75,434
310,0
-118,8
191,2
-257,1
-40,4
-297,5
-106,3
Patagonia
4,144
5,7
2,3
8,1
-29,2
-8,8
-38,0
-29,9
Ross
0,904
0,2
0,0
0,2
-66,6
-35,0
-101,7
-101,5
Weddel
2,407
10,4
-0,3
10,0
-88,7
-55,1
-143,9
-133,8
Total
362,148
3251,4
-482,4
2769,0
-2808,7
-770,8
-3579,5
-810,5

Tableau 39 : Bilan d’énergie de la solution libre M2 (énergie en Giga watts)

Comparons ces résultats avec ceux obtenus par Le Provost et Lyard [1997] grâce aux champs de vitesses calculés avec FES94.1 pour les grands bassins océaniques :
 

Bassin
Surface 106km2
P
D
Balance
Arctic
11,6
-7
-98
-105
AtlanticN
42,5
194
-850
-656
AtlanticS
48,3
680
-261
419
Indian
73,6
541
-239
302
PacificN
73,4
518
-441
77
PacificS
107
281
-120
161

Tableau 40 : Bilan énergétique obtenus par Le Provost et Lyard pour M2 avec le modèle FES94.1

Il faut noter que nos bassins et ceux de Le Provost et Lyard ne sont pas rigoureusement identiques (cf. surfaces) car les maillages utilisés pour calculer FES94.1 et FES2000 sont différents. C’est pourquoi, nous avons essayé dans la mesure du possible de reconstruire les quatre principaux bassins entre les deux solutions en sommant les bilans énergétiques. Nous avons obtenu les quatre bassins suivants : Arctique, Atlantique, Indien et Pacifique. Afin de quantifier les différences entre les énergies, nous avons calculé les rapports de l’apport énergétique et de la dissipation déduites des vitesses de FES2000 avec ceux calculés à partir de FES94.1 sur les quatre bassins (Tableau 41).
 

Bassin
P (FES94)
P (FES2000)
Rapport P
D (FES94)
D (FES2000)
Rapport D
Arctique
-7
-10
1,43
-98
-169
1,72
Atlantique
874
1405
1,61
-1111
-1806
1,63
Indien
541
215
0,40
-239
-115
0,48
Pacifique
799
1159
1,45
-561
-1490
2,66

Tableau 41 : Rapport des bilans énergétiques des solutions libres de FES2000 par rapport aux solutions de FES94.1

Les résultats du Tableau 41 sont assez difficiles à interpréter car les zones de comparaisons entre FES94.1 et FES2000 ne sont pas rigoureusement identiques. Cependant, il ressort que pour trois des quatre bassins (Arctique, Atlantique et Pacifique) les apports astronomiques et les dissipations calculées avec les solutions FES2000 sont plus élevées que pour celles calculées avec FES94.1. La dissipation obtenue avec FES94.1 a été validée par d’autres travaux indépendants [Church et al., 2000; Egbert, 1997; Kantha et al., 1995; Munk, 1997]. C’est donc la dissipation obtenue avec les solutions libres FES2000 qui est trop forte (3579,5 GW ce qui est environ 78% trop fort comparé aux 2010 GW obtenus avec FES94.1). Cette constatation est confirmée par notre étude sur les amplitudes de FES2000 (paragraphe 13.5.3) qui nous a montré que les amplitudes de marées sont trop fortes et donc que le phénomène dissipatif est trop faible. Nous savons que la dissipation est proportionnelle au cube de la vitesse. Les vitesses obtenues avec FES2000 sont donc trop fortes, ce qui conduit à penser que les coefficients de frottement qui sont introduit dans le calcul itératif de CEFMO sont trop faibles. Sous forme tensorielle, les coefficients de frottement pour l’onde M2 s’écrivent sous la forme de (13.147). Les coefficients r, r’, r’’ et r’’’ sont des fonctions de R et R’ (deux termes introduits dans la Partie 1 de cette thèse) dans les équations quasi-linéarisées [Le Provost and Vincent, 1986]. Dans le cas de l’onde dominante M2 nous avons :

(13.150)

Pour avoir un ordre d’idée de la distribution des coefficients de frottements issus de notre code de calcul, nous avons tracé les valeurs de R et R’ après 10 itérations du code CEFMO pour l’onde M2 (Figure 96 et Figure 97).


Figure 96 : Coefficient de frottement R


Figure 97 : Coefficient de frottement R’

Nous voyons bien sur les deux figures que R est plus important que R’ d’un facteur 10 environ. Dans une première approximation, nous pouvons considérer le frottement dans notre code comme linéaire et le comparer à des frottements linéaires obtenus par d’autres méthodes. Ainsi Church et al. [2000] ont déduit une carte du coefficient de frottement sur le fond océanique (Figure 98) à partir de la dissipation globale déduite du modèle FES99 (Church, communication personnelle, 2000).


Figure 98 : Coefficient de frottement linéaire calculé par Church et al.

La Figure 96 comparée à la Figure 98 rend bien compte de la trop faible valeur de notre coefficient de frottement R en particulier en plein océan où les différences peuvent atteindre un facteur 100. Par quoi peuvent donc s’expliquer ces différences ? Dans notre modélisation des marées, nous ne prenons en compte que les vitesses de la marée barotrope. Or, nous ne tenons pas compte de l’effet barocline de la marée alors que la génération d’ondes internes due à l’effet barocline n’est pas à négliger [Munk, 1997; Ray and Mitchum, 1996]. Les vitesses de marée doivent donc être moins importantes que celles que nous calculons dans CEFMO. Les coefficients de frottement sont donc sous estimés. Il est cependant difficile de dire de combien est cette sous estimation. Le facteur 100 n’est-il pas trop élevé ?

Nous avons voulu estimer ce facteur en utilisant dans notre code les coefficients de frottement déduits par Church et al. Le champ de Church et al. est limité aux zones de grands fonds et aux zones non polaires. C’est pourquoi, nous avons remplacé notre coefficient R calculés précédemment pour obtenir les dénivellations libres de FES2000, par ce nouveau champ de frottement sur tous les fonds de profondeurs plus importantes que 500 m et dans les zones polaires. De ce nouveau champ de frottement, nous avons déduit des vitesses à partir de nos solutions libres en dénivellation. Mais la simple technique de recollement des deux champs de frottement que nous avons mis en place, pose des problèmes de raccordement (pas de continuité) entre ces deux champs. C’est pourquoi les calculs que nous avons effectués sont entachés d’erreurs numériques trop importantes pour pouvoir faire des comparaisons plus avancées. Par manque de temps nous n’avons pas pu pousser plus loin nos investigations dans cette voie.

13.7 Bilan sur FES2000

Les apports de l’étude locale dans la Mer Jaune et les Mers de Chines pour le modèle CEFMO nous ont permis d’améliorer notre code à l’échelle globale. Ainsi, nous avons régénéré une nouvelle bathymétrie mondiale adapté à nos besoins et nous avons recréé un nouveau maillage global éléments finis trois fois plus précis que celui utilisé pour calculer FES94.1, FES95.2, FES98 et FES99. L’adaptation de CEFMO à ce nouveau maillage couplé aux nouvelles puissances de calculs fournies par l’IDRIS a permis de calculer de nouvelles solutions libres de marées FES2000. Le but de notre étude est d’obtenir une solution de marée globale indépendante de toutes données de terrain. Cependant, ces nouvelles solutions libres sont trop élevées en dénivellation. Le processus itératif de CEFMO qui lie les calculs des vitesses, de coefficients de frottement et des dénivellations, nous permet d’avancer que les vitesses sont trop élevées et donc que les coefficients de frottement sont trop faibles. Cette constatation est confirmée par l’étude de Church et al. [2000]. L’hypothèse de vitesse barotrope dans le calcul de CEFMO est mise en cause. En effet, les calculs ne tenant pas compte des ondes internes, ils sont biaisés par rapport au phénomène physique de la dissipation des marées océaniques. Le manque de temps ne nous a permis d’explorer plus avant cet aspect énergétique afin d’identifier le manque à apporter dans la paramétrisation de nos équations.
Back | Next

Title: Thèse de Fabien Lefèvre
Issue: Version 1.0
Date: 29/09/2000