Formulation générale et conditions aux limites

Problématique

L’étude des ondes internes et du mélange turbulent qu’elles induisent, nécessite une modélisation numérique à haute résolution. D’autant plus que dans le cas de simulations numériques directes à l’échelle du laboratoire, on souhaite résoudre la majeure partie de la cascade turbulente jusqu’aux échelles de la dissipation. De plus si l’on veut représenter correctement le déferlement de ces ondes ainsi que la génération et la propagation des ondes internes non-linéaires (solitons), il est nécessaire de prendre en compte les effets non-hydrostatiques. Si l’on souhaite représenter explicitement la génération de marée interne (forcées par des ondes de marée de surface), il est également nécessaire de représenter explicitement les ondes de surface et donc d’être dans des conditions dites de « surface libre ». La représentation numérique de ces processus nécessite donc des modèles de surface libre, très haute résolution et non-hydrostatiques.
On se dirige aujourd’hui vers de nouvelles infrastructures de calcul. Ainsi, l’avènement des systèmes massivement parallèles a permis d’améliorer significativement les performances de ce type de modèles numériques mais soulève de nouveaux problèmes nécessitant d’adapter les algorithmes.
Le développement récent des architectures GPUs ou Xeon Phi comme accélérateurs de calcul avec les CPU classiques pour des applications de calculs intensifs semble très prometteur. Ils présentent de nombreux intérêts que ce soit en termes d’accélération de calculs (nombre important de coeurs de calcul), de coût ou de consommation énergétique (greenHPC). Cependant, les performances de ce type de systèmes dépendent fortement des algorithmes à traiter et nécessitent donc le développement de nouveaux algorithmes plus adaptés.
Le calcul des différentes composantes de la pression (hydrostatique, non-hydrostatique et pression induite par les anomalies de surface libre) est un souci majeur des algorithmes nonhydrostatiques. Pour calculer ces différentes composantes, la correction de pression (CP) est une des méthodes les plus populaires dans les modèles de circulation océanique (Kanarska et al., 2007; Marshall et al., 1997a, 1997b). Cependant, ce schéma de correction de pression entraine la résolution coûteuse d’une équation de Poisson 3D. Pour minimiser les temps de calcul, une méthode de pas de temps séparés est souvent utilisée. Elle consiste à intégrer simultanément deux modes dynamiques : un mode externe dédié à la dynamique « barotope » à haute fréquence et un mode interne dédié à la dynamique basse fréquence barocline. La résolution explicite des ondes de surface « courtes » (surface libre) induit néanmoins des difficultés supplémentaires. Leur modélisation explicite requière généralement des pas de temps externe (mode « barotrope ») très petit. De plus, l’élévation de surface est couplée avec la composante non-hydrostatique de la pression, le pas de temps interne (mode « barocline ») est alors contraint par le pas de temps externe. Les gains de calcul que peut apporter une méthode de pas de temps séparés sont alors fortement limités. Ce type de modélisation nécessite donc des ressources de calculs très importantes. En ce sens, une des problématiques clés pour ce type d’approche est de trouver des schémas numériques plus performants.

Bathymétrie Mouvante

L’utilisation de topographies oscillantes est une des méthodes couramment employées dans les laboratoires pour étudier les interactions entre la marée et la topographie. Gerkema et Zimmerman (2008) ont montré qu’une telle représentation physique pouvait être considérée comme équivalente au mouvement de la marée au dessus d’un mont. Afin de faciliter la comparaison avec les expériences de laboratoire, dans lesquelles le forçage est réalisé par l’oscillation de la topographie, une version du modèle SNH à bathymétrie mouvante a été développée (Auclair et al., 2014). Cet algorithme ouvre de nouvelles opportunités de configurations océaniques dans lesquelles la bathymétrie évolue dans le temps à des échelles temporelles au moins comparables à celle de la dynamique océanique. On peut ainsi envisager des imbrications de modèles avec suivi Lagrangien, l’étude de catastrophes sismiques avec génération de tsunami ou encore l’étude de l’évolution sédimentaire près des côtes.

Condition au fond

Des conditions de non-glissement (no-slip) sont utilisées au fond, imposant des vitesses verticales ( ) et horizontales ( , α=1 (pour x) et 2 (pour y)) nulles en z = -H pour une bathymétrie fixe:

Modélisation non-hydrostatique

Mode non-hydrostatique « non-Boussinesq » (NBQ)

Le noyau non-hydrostatique basé sur une correction de pression n’étant pas directement utilisé dans le cadre de la présente thèse2, le lecteur est renvoyé pour plus d’informations sur sa description  proposée par Auclair et al. (2011). Le second noyau non-hydrostatique dit noyau « non-Boussinesq » est maintenant décrit plus en détails. Ce noyau résout le système général d’équations (SE1) et a donc nécessité une réécriture des modes internes et externes du modèle pour la quantité de mouvement (et non plus pour la vitesse). Seconde évolution importante : l’équation de continuité traite de la conservation de la masse et non plus de la conservation du volume comme cela était le cas sous l’hypothèse de Boussinesq. Enfin, ce n’est plus simplement la pression qui est écrite comme la somme d’une composante hydrostatique et une composante non-hydrostatique mais plus généralement la masse volumique qui est décomposée sous la forme :

Discrétisation numérique

Discrétisation temporelle

SNH est un modèle à surface libre, les ondes externes sont donc également explicitement représentées. Afin de représenter correctement la propagation de toutes ces ondes, pseudoacoustiques, externes et internes, qui se propagent à des vitesses très différentes, tout en minimisant les temps de calculs, une méthode de pas de temps séparés, ou « time-splitting » est utilisée. Pour résoudre correctement les ondes les plus rapides, la résolution temporelle du mode NBQ (ΔtNBQ) est plus fine que celle du mode externe (Δte), elle-même plus fine que celle du mode interne (Δti), comme présenté sur la Figure 2.20. Un schéma leap-frog, ou saute-mouton, est utilisé pour la discrétisation temporelle. La vitesse de phase des ondes pseudo-acoustiques, cs, peut être réduite artificiellement, rendant les processus hautes fréquences associés à la compressibilité sans signification physique, mais préservant une dynamique lente cohérente et permettant une relaxation des contraintes CFL sur les pas de temps. La représentation des ondes acoustiques nécessite donc des pas de temps de mode NBQ très petits, entrainant des coûts de calculs important mais peut être compensée par une physique simplifiée.

Discrétisation spatiale

Un schéma en différences finies centré sur une grille Arakawa-C est utilisé pour la discrétisation du modèle. Une discrétisation régulière est adoptée sur l’horizontale. Afin de s’adapter à des topographies complexes et aux variations de la surface libre, une coordonnée sigma s(x,z,t) est utilisée, particulièrement bien adaptée à l’étude des effets topographiques sur la dynamique océanique. 20 à 100 niveaux « s » sont utilisés selon la configuration. L’utilisation de cette coordonnée peut cependant poser des problèmes dans le cas d’une représentation simultanée de gradients topographiques et de gradients de masse volumique forts entrainant des erreurs de troncatures. La grille Arakawa-C peut être décomposée comme présenté sur la Figure 2.21. Les variables (ρ, T, S, P) sont calculées au centre de la maille (i,k), appelé « point de masse ». Les variables de vitesse sont calculées à des demi-niveaux horizontaux (u) et verticaux (w).

Configurations

Différentes configurations numériques réalisées avec le modèle SNH sont utilisées au cours de cette thèse. Le chapitre 3 consiste en une étude de régimes d’ondes internes et se base sur des simulations numériques directes (quasi-DNS) de résolution millimétrique à l’échelle du laboratoire qui permettent de résoudre explicitement toutes les structures tourbillonnaires potentiellement présentes. Ces configurations se basent sur les expériences physiques et numériques réalisées par Dossmann et al. (2013a) et Floor et al. (2009). A contrario, les configurations du chapitre 4 s’appuient sur des simulations de grande échelle (LES) à haute résolution permettant de résoudre les grandes échelles de l’écoulement mais n’allant pas jusqu’aux échelles de la turbulence.

Comparaison entre simulations numériques et expériences de laboratoire

La réalisation de simulations numériques directes et d’expériences de laboratoire effectuées par Floor (2009) sur des configurations physiques identiques a permis d’évaluer la pertinence des hypothèses sous-jacentes au modèle numérique. Cette évaluation a été effectuée, dans le cas d’une stratification linéaire, pour l’expérience présentée sur la Figure 2.23, sur les variables de masse volumique et de vitesse. La comparaison s’est révélée concluante, les approches numériques et expérimentales décrivent de manière comparable la dynamique des ondes internes, permettant l’utilisation du modèle dans cette configuration.

Configurations LES

Le chapitre 4 de cette thèse porte sur la dynamique des ondes internes à l’échelle océanique. Pour cela des simulations de grande échelle (LES) à haute résolution sont réalisées permettant de résoudre explicitement la dynamique des ondes internes. On passe donc à des domaines d’une largeur de plusieurs dizaines de kilomètres et à des résolutions métriques.

Force de Coriolis

Une version bidimensionnelle dans un plan vertical ( ) du modèle est utilisée. Des conditions cycliques dans la direction transverse sont appliquées aux configurations LES des cas d’études réalistes (détroit de Gibraltar et golfe du Maine). La vitesse transverse v est donc non nulle. On se place alors dans un référentiel tournant ( ).

Paramétrisation de la dissipation

Les simulations sont qualifiées de configurations LES dans la mesure où la dynamique des ondes internes est résolue explicitement bien que la résolution de l’ordre du mètre ne soit plus suffisante pour résoudre les petites échelles de la turbulence. Ces processus de petites échelles sont donc pris en compte en augmentant la viscosité. Pour cela les valeurs de la viscosité cinématique (ν = 2.10-5 m²/s) et du coefficient de diffusion de la masse volumique (Kρ = 2.10-7 m²/s) sont augmentées. Le coefficient multiplicatif du schéma d’advection “upstream” est également augmenté à une valeur de 0.015. L’objectif reste de représenter explicitement les processus physiques (ondes internes, instabilités …) induisant du mélange turbulent, et donc de calibrer ces paramètres au niveau le plus faible. Pour cela, une étude de sensibilité a été faite, permettant de fixer ces paramètres de dissipation à des valeurs minimales tout en assurant la stabilité du modèle. La seconde viscosité λ(2) est fixée à 102 m²/s dans cette configuration. Sa valeur a de nouveau peu d’impact sur la dynamique lente des ondes internes.

Diagnostics

Au cours de cette thèse, plusieurs diagnostics et post-traitement ont été développés pour l’analyse de ces configurations numériques.

Coordonnées isopycnales

Les ondes internes se propagent le long de surfaces isopycnales qui ne sont pas représentées explicitement par le modèle. Il est donc nécessaire de mettre au point un diagnostic permettent d’extraire la position des coordonnées isopycnales à partir des champs de masse volumique. Pour représenter l’évolution de ces surfaces isopycnales, un changement de variable est alors effectué via un post-traitement par une méthode de traçage (« tracking »): ρ(z)→z(ρ). Cette méthode de traçage permettra notamment d’effectuer des diagrammes Hovmöller représentant l’évolution spatiale et temporelle des ondes internes.

Décomposition en modes verticaux

Une des problématiques principales a été d’identifier les modes verticaux d’ondes internes et d’isoler chacun de leurs signaux à partir des champs de vitesses et de masse volumique résolus par le modèle.
La première difficulté pour identifier ces modes verticaux est d’avoir une référence analytique précise en termes de longueur d’onde (λn) , de vitesse de propagation (cn) et de structure verticale Wn(z). Pour cela l’approche modes normaux pour une stratification non-linéaire décrite dans la section (1.3.3) est utilisée. L’équation de Sturm-Liouville est résolue numériquement pour un profil de stratification donné : N(t = 0,x1,z). Ce profil de stratification est calculé à partir du champ de masse volumique extrait du modèle loin de la zone de génération et au dessus de celle-ci à l’instant t = 0.
Dans le cas de simulation avec une stratification complexe spatialement hétérogène comme celle du détroit de Gibraltar, il est nécessaire de prendre en compte la variabilité spatiale et temporelle du champ de masse volumique. L’équation de Sturm-Liouville est alors résolue en post-traitement pour chaque pas de temps externe et à chaque point de l’espace. Les vitesses de propagation théoriques déduites de cette méthode analytique permettront en autres d’analyser les diagrammes Hovmöller.
Les structures verticales associées servent de référence pour les deux méthodes statistiques présentées ci-dessous.

Outils statistiques

Pour pouvoir étudier la propagation de ces modes normaux et mesurer leurs amplitudes respectives à partir de nos sorties de modèles, divers outils statistiques ont été utilisés. Nous avons fait le choix d’utiliser le champ de vitesses verticales (w) pour identifier les modes normaux.

Principe de similitude : du laboratoire à l’échelle océanique

Les configurations physiques et numériques à l’échelle du laboratoire, présentées précédemment, sont inspirées de configurations océaniques. En ce sens, elles sont utilisées pour comprendre et décrire des processus physiques que l’on considère semblables à ceux observés dans l’océan. Cependant, il n’est pas possible de reproduire simultanément la totalité des caractéristiques physiques et géométriques des écoulements océaniques à l’échelle du laboratoire. La question qui se pose alors est de savoir si l’on peut réellement reproduire, à l’échelle du laboratoire, des régimes d’ondes internes similaires à ceux des échelles océaniques. Pour cela, il est nécessaire de mettre en place un principe de similitude permettant de conserver la dynamique des ondes internes tout en modifiant l’échelle de l’écoulement. De nombreux principes de similitudes existent dans la littérature, chacun d’eux étant façonnés pour répondre à des contraintes et des objectifs spécifiques.
Nous présenterons l’un d’entre eux, appliqué à la génération d’ondes internes solitaires dans le détroit de Luçon. Puis, nous présenterons, le principe de similitude que nous avons choisi pour l’étude de la propagation d’ondes internes dans un environnement supercritique.