Approche numérique et expérimentale pour une
meilleure description physique des processus de
subméso-échelle

Courant Nord : de l’échelle du bassin aux côtes varoises 

La circulation du bassin occidental méditerranéen a été étudiée à l’aide de nombreuses campagnes de mesures et d’observations satellites. Les travaux de Millot (1999) ont montré que le CN, qui est la part de la circulation cyclonique méditerranéenne au nord du bassin, est composé des MAW1 , LIW2 , WMDW3 et WIW4 quand elles existent (la description de la formation et écoulement de ces masses d’eau (figure 1.1) pourra être trouvée dans la bibliographie de C. Millot et ne sera pas reprise ici ). Le CN s’écoule le long du plateau continental, longeant ainsi les côtes françaises et espagnoles jusqu’au îles Baléares (Millot, 1991; Lopez-Garcia et al., 1994) où il se sépare en deux branches: une première longeant la côte jusqu’au détroit de Gibraltar, et une seconde remontant au nord des îles Baléares. La variabilité saisonnière des courants ECC 5 et WCC 6 qui composent le CN à leur jonction dans le bassin ligure, joue un rôle sur le positionnement et l’intensité du CN. Bien que la variabilité du WCC puisse être considérée négligeable (Astraldi and Gasparini, 1992), la circulation des LIW dans le canal Corse est trois fois plus forte en hiver (La Violette, 1994). De plus, l’ECC étant abrité des vents froids hivernaux, les LIW se refroidissent moins que les eaux qui composent le WCC. Cette variabilité saisonnière implique un CN plus chaud que les eaux du bassin de l’hiver au printemps (figure 1.2). Sur la figure 1.2, représentant la SST satellite moyennée sur un hiver, le CN longe la côte d’Azur, s’écoulant le long des côtes varoises avant de suivre le talus du GdL. Le rôle de barrière que le courant peut avoir entre les eaux du large et les eaux côtières est bien visible ici : une nette différence de température de surface est présente entre les eaux du plateau, qui sont plus froides, et celles de la veine de courant. Cette augmentation du transport de l’ECC en hiver est également à l’origine d’un resserrement (environ 25 km) et d’un approfondissement (450 m) de la veine, accompagnés d’une accélération (0.5 m.s-1) (Taupier-Letage and Millot, 1986; Millot, 1999) et d’une forte activité méso-échelle. De l’été à l’automne, le CN est plus large (40 km) avec une profondeur de 250 m (Albérola et al., 1995; Petrenko, 2003; Conan and Millot, 1995). Durant cette période où l’activité à méso-échelle est plus faible, le CN ralentit (0.2 m.s-1) et ses propriétés hydrologiques correspondent à celles des MAW. L’activité à méso-échelle du CN, comme les tourbillons et les méandres, est maximum en hiver (Sammari et al., 1995; Conan and Millot, 1995). En effet, s’écoulant très près des côtes en hiver, l’intéraction entre le fond et le courant est plus importante à cette saison, entraînant l’apparition 1. Mediterranean Atlantic Water 2. Levantine Intermediate Water 3. West Mediterranean Deep Water 4. Winter Intermediate Water 5. East Corsica Current 6. West Corsica Current 13 14 Chapitre 1. Courant Nord : de l’échelle du bassin aux côtes varoises Figure 1.1 – Circulation des masses d’eau dans le bassin occidental méditerranéen, de haut en bas : des eaux Atlantiques modifiées (MAW), des eaux intermédiaires Levantines (LIW) et des eaux profondes méditerranéennes (WMDW). (figures issues de (Millot, 1999)) 15 Figure 1.2 – Température (◦C) de surface moyenne observée par satellite durant l’hiver 2013- 2014. L’isobathe 1000 m est superposée en noir. d’instabilités (Sammari et al., 1995; Petrenko, 2003; Flexas et al., 2002) pouvant donner lieu à des méandres atteignant 10 à 100 km (Conan and Millot, 1995). Le CN est un courant géostrophique (Conan and Millot, 1995), dû à l’ajustement isopycnal entre les différentes masses d’eaux qui le composent, fortement influencé par la bathymétrie. Du bassin Ligure aux côtes espagnoles, le courant rencontre une zone où le plateau continental est quasiinexistant avec des fonds très escarpés couverts de profonds canyons jusqu’à l’entrée du GdL où la présence d’un vaste plateau continental peu profond (environ 150 m) va l’éloigner de la côte. Situées entre ces deux zones aux caractéristiques bathymétriques différentes, les côtes varoises sont une zone de transition. L’observation du CN au large de Toulon à l’aide de radar HF (par exemple sur la figure 1.3 où le courant de surface observé est présenté) montre un CN en moyenne bien établi au sud de l’isobathe 1000 m, avec une direction Ouest (figure 1.3a, haut) mais également une activité à petite échelle non négligeable observée sur des moyennes quotidiennes (figure 1.3b, panneaux inférieurs). En effet, sur 4 journées différentes, les observations radar montrent un courant de surface déstructuré par les effets d’un vent de Nord pour le 10/05/2013, par la présence de méandre (13/03/2013 et 10/01/2013) ou encore une structure tourbillonnaire au large de Toulon provoquant un décrochage du CN vers le Sud le 12/03/2013. Ces événements, modifiant la position du CN au large de Toulon, vont par la suite avoir un impact sur son écoulement en aval. De nombreuses études ont été réalisées sur le devenir du CN en amont (Taupier-Letage and Millot, 1986; Albérola et al., 1995; Sammari et al., 1995; Bethoux et al., 1982) et en aval (Albérola and Millot, 2003; Conan and Millot, 1995; Flexas et al., 2002; Petrenko, 2003; Rubio et al., 2009; Lapouyade and Durrieu de Madron, 2001) des côtes varoises. Bien que cette zone soit considérée comme importante dans l’accélération et l’orientation du CN (Ourmières et al., 2011) et un lieu de génération d’instabilités du courant dû au trait de côte (Guihou et al., 2013), peu d’études ont été menées au cœur de cette zone. Concernant la circulation dans la baie de Hyères, très peu d’observations sont disponibles. Cependant, cette baie peu profonde et semi-fermée par la présence de trois îles est sujette aux 16 Chapitre 1. Courant Nord : de l’échelle du bassin aux côtes varoises (a) Intensité et direction moyennées sur l’année 2013 (b) Direction du courant sur une moyenne horaire. L’échelle des flêches représentant l’intensité. Haut Gauche : 10/05/2013 – Haut Droite : 13/03/2013 – Bas Gauche : 10/01/2013 – Bas Droite : 12/03/2013 Figure 1.3 – Courants (m/s) de surface observés par radar HF au large de Toulon. 17 Figure 1.4 – Trajectoires de dérive de bouées lagrangiennes. Plusieurs lâchées ont été effectués entre le 16 Février et le 21 Juillet 2015 dans le cadre du projet CADOR. Les positions initiales des flotteurs sont symbolisées par des cercles. vents dominants dans la zone ainsi qu’au CN passant au Sud. L’utilisation de bouées dérivantes, dans le cadre d’une campagnes d’observations menée au laboratoire MIO, a permis de mettre en évidence l’influence du vent sur le transport dans cette région, comme par exemple les lâchés réalisés les 11/03 et 21/07/2015 par un vent à dominance Est (figure 1.4). L’impact des fleuves locaux, comme le Gapeau ou encore la Maravène, est négligeable dans la zone d’étude sauf en période de fortes pluies où la crue peut engendrer un apport considérable d’eau dans la baie. Ces crues lorsqu’elles sont suivies d’un vent à dominante Nord entraînent un export d’eau vers le Sud avec des vitesses pouvant atteindre en moyenne 40 cm/s (figure 1.4, lâchés de bouées les 16/02 et 18/03/2015). Le 18/02/2015 avec un vent de Nord (figure 1.4), la bouée a pris une direction Sud-Est avant de virer de cap et de partir plein Ouest vers la côte, mettant en évidence des cisailles horizontales de courant très proches des côtes, à moins de 2 km. Il ressort de ces observations une circulation dans la baie très influencée par les conditions météorologiques, interagissant aussi avec le courant de bord (CN) et sa variabilité.. En plus des observations, la modélisation numérique est un outil utilisé pour l’étude du CN le long des côtes françaises et espagnoles mais également pour l’étude de processus, comme la convection profonde en dans le bassin nord-ouest méditerranéen (Estournel et al., 2016; Madec et al., 1991), l’interaction du CN avec le talus du GdL (Echevin et al., 2003) ou encore l’activité à méso voire sub-méso-échelle (Bouffard et al., 2012; Garreau et al., 2011). Les travaux de Guihou et al. (2013) (figure 1.5) montrent la concordance des courants de surface modélisés à l’aide de la configuration GLAZUR64 (configuration décrite en section 2.2) avec la température et la concentration en chlorophylle-a issues d’observations satellites lors d’un événement tourbillonnaire au large des îles d’Hyères le 03/04/2011. Le positionnement et le phasage de ce tourbillon modélisé avec les observations montre que la modélisation numérique du CN dans cette zone est correctement établie à cette échelle. Comme le CN est un courant de densité, l’ordre de grandeur de sa méso-échelle est celui du rayon de déformation interne de Rossby Rd (1er mode barocline), déterminant l’échelle à partir de laquelle la force de Coriolis devient aussi importante que la force de flottabilité. À partir de travaux de modélisation, il a été montré que Rd est de l’ordre de 5-6 km en hiver pour le 18 Chapitre 1. Courant Nord : de l’échelle du bassin aux côtes varoises Figure 1.5 – Fond de carte : Données satellite pour le 03 Avril 2011. Gauche : Concentration en chlorophylle a (mg/m3 ) (produit L4, projet MyOcean). Droite : SST (◦C) (produit L4, projet MEDSPIRATION). Les courants de surface simulés à l’aide de la configuration NEMOGLAZUR64 sont superposés avec les vecteurs (en noirs) sur les deux cartes. (figures issues de Guihou et al. (2013)) CN (Guihou, 2013), définissant ainsi la résolution minimum des configurations, de l’ordre du kilomètre, nécessaire pour simuler cette méso-échelle. Dans la baie semi-fermée, il a été montré dans les mêmes travaux que le Rd est inférieur à 3 km en hiver, impliquant la nécessité d’une résolution numérique plus fine pour résoudre la méso-échelle, voire la subméso-échelle de la zone. Dans le but de poursuivre l’étude du CN ainsi que le couplage de celui-ci avec la circulation de la baie, il a été choisi dans ces travaux d’utiliser la modélisation numérique avec un couplage online de deux configurations modèles : une première couvrant la Côte d’Azur et le GdL à l’échelle kilométrique et une seconde, avec une résolution plus fine, couvrant la baie et le large des îles d’Hyères. Les simulations obtenues seront ensuite comparées aux observations disponibles dans les zones d’étude

Modélisation océanique 

Les modèles numériques sont des outils désormais incontournables et complémentaires aux observations océaniques pour la description des différents processus océaniques. Les modèles de circulation 3D sont largement utilisés pour étudier la dynamique méditerranéenne, différentes configurations ont été développées utilisant différents modèles : parmi les plus récentes on peut citer celles utilisant les codes NEMO (Brossier et al., 2014), MARS3D (Pairaud et al., 2011), SYMPHONIE (Marsaleix et al., 2008) et ROMS (Escudier et al., 2016). Les modèles numériques d’océan sont basés sur la résolution des équations primitives discrétisées sur des grilles de calculs, dont la résolution est liée aux dimensions des phénomènes étudiés. Les récentes innovations informatiques permettent de modéliser des processus physiques à méso- (voire sub-méso) échelles, une nécessité pour la description de la dynamique en zone côtière. Lors de ces travaux, il a été choisi d’utiliser le code de calcul NEMO, dont la description est faite ci-après. 2.1 Le code de calcul NEMO NEMO (Nucleus for European Modelling of the Ocean) permet la modélisation de la physique de l’océan, à l’aide du module OPA 1 , développé par Madec et al. (1998), mais comporte également un modèle de glace, LIM 2 , ainsi qu’un module de biogéochimie, TOP 3 . Comme son nom l’indique, NEMO est développé principalement en Europe, initialement au LOCEAN 4 , mais il est également utilisé par la communauté scientifique internationale, comme le NOC 5 et le Met Office qui sont les centres nationaux océanographique et météorologique de Grande Bretagne ou le groupe canadien MEOPAR modélisant la mer des Salish pour la gestion des risques dans les zones côtières lors de tempêtes. NEMO est également employé en configurations globales opérationnelles, par Mercator-Océan 6 par exemple. De plus, les récentes optimisations du code et l’ajout d’un module permettant l’imbrication de zoom online rendent ce code de calcul propice à la modélisation côtière réaliste à haute résolution. Le système d’équations primitives résolu par NEMO, ainsi que les paramétrisations propres aux configurations côtières utilisées dans ces travaux sont décrits ci-dessous. La documentation complète du code pourra éventuellement fournir au lecteur des détails sur le code (Madec, 2008). 1. Océan PArallélisé 2. Louvain-la-neuve Ice Model 3. Tracer in the Ocean Paradigm 4. Laboratoire d’Océanographie et de Climatologie, Expérimentation et Analyse Numérique 5. National Oceanography Centre 6. http ://www.mercator-ocean.fr : Centre français d’océanographie opérationnelle. 19 20 Chapitre 2. Modélisation océanique 

Équations primitives

L’océan est un fluide pouvant être décrit à l’aide des équations primitives, qui sont des équations d’état non linéaires reliant température, salinité et vitesse du fluide. Le système d’équations est composé de 7 inconnues : les 3 composantes de la vitesse, U~ , la température, T, la salinité, S, la pression, P, et la masse volumique, ρ. Ce système étant composé d’équations non-linéaires, il est nécessaire pour le résoudre de poser des hypothèses et des approximations simplificatrices qui sont liées aux échelles spatio-temporelles auxquelles on s’intéresse. 

 Équations générales

Les équations du mouvement, ou encore équations de Navier-Stokes, décrivent le mouvement d’un fluide. Elles sont basées sur les principes de la dynamique reliant l’accélération d’une particule fluide à la somme des forces qui lui sont appliquées. Dans un repère orthonormé (~i,~j,~k) relié à la Terre, ce système devient : dU~ dt = − 1 ρ ∇~ p + ~g − 2Ω~ × U~ + DU + F U (2.1) Où : — t est le temps — U~ représente le vecteur tridimensionnel de vitesse — ρ est la densité du fluide — p est la pression — ~g est l’accélération gravitationnelle — Ω~ représente la rotation de la Terre — F représente les forçages de surface — D représente la paramétrisation des processus à petite échelle — ∇~ est l’opérateur de dérivée vectorielle — L’opérateur dX dt = ∂X ∂t + U. ~ ∇~ X représente la dérivée particulaire. L’équation de conservation de la masse s’écrit : 1 ρ dρ dt + ∇.U~ = 0 (2.2) Les équations de conservation du sel et de la température sont respectivement : dS dt = −∇.(SU~ ) + DS + F S (2.3) dT dt = −∇.(TU~ ) + DT + F T (2.4)

Hypothèses simplificatrices — Approximation d’une Terre sphérique : les surfaces géopotentielles sont supposées être des sphères concentriques, donc la gravité est considérée parallèle au rayon de la Terre. — Approximation de couche fine : la profondeur de l’océan est considérée négligeable par rapport au rayon de la Terre. 2.1 Le code de calcul NEMO 21 — Hypothèse de fermeture turbulente : les flux turbulents, qui représentent l’effet des petites échelles sur les grandes, sont exprimés en terme de processus à grande échelle. — Hypothèse de Boussinesq : les variations horizontales de densité sont considérées négligeables devant la valeur moyenne. Seule la variation verticale de cette grandeur est prise en compte pour sa contribution aux forces de flottabilité. On estime la masse volumique comme étant la somme d’une masse volumique de référence, ρ0 et d’une perturbation, ρ 0 : ρ(p, T, S) = ρ0 + ρ 0 (x, y, z, t) (2.5) En utilisant cette approximation dans l’équation du mouvement (Eq.2.1) et une analyse en ordre de grandeur, on obtient : dU~ dt = − 1 ρ0 ∇~ p + ~g − 2Ω~ × U~ + DU + F U (2.6) — Hypothèse hydrostatique : l’équation du moment sur la verticale est réduite à l’équilibre entre le gradient vertical de pression et la force de flottabilité. Cette hypothèse induit que les processus convectifs ne sont pas pris en compte, ils doivent donc être paramétrés. En effet, la pression appliquée à une particule fluide est considérée comme la somme de la pression atmosphérique p0 et la pression de la colonne d’eau au-dessus, p 0 (x, y, z, t) : p = p0(x, y, z, t) + p 0 (x, y, z, t) (2.7) La pression atmosphérique est en général négligée, ce qui ramène l’expression de la pression totale p à la simple composante hydrostatique p 0 . Considérant également l’hypothèse d’incompressibilité, la variation de pression sur la verticale est associée à une variation de masse volumique, équilibrant ainsi le gradient de pression vertical et la force de flottabilité : ∇~ zp = ∂p ∂z = −ρg soit pz = Z z 0 −ρgdz (2.8) — Hypothèse d’incompressibilité : la divergence (3D) de la vitesse est nulle. ∇.U~ = 0 ou ∂u ∂x + ∂v ∂y + ∂w ∂z = 0 (2.9) L’équation du mouvement étant réduite à l’équation hydrostatique (Eq.2.8) pour la composante verticale, la non divergence de la vitesse du fluide permet de déduire la variation de vitesse verticale ∂w ∂z des variations de vitesses horizontales (Eq.2.9). Suivant ces hypothèses le système d’équations primitives précédent devient : ∂U~ h ∂t = −  (∇ × ~ U~ ) × U~ + 1 2 ∇(U~ 2 )  h − f ~k × U~ h − 1 ρ0 ∇~ hp + DU + F U (2.10) ∂p ∂z = −ρg (2.11) ∇.U~ = 0 (2.12) 22 Chapitre 2. Modélisation océanique ∂S ∂t = −∇.(SU~ ) + DS + F S (2.13) ∂T ∂t = −∇.(TU~ ) + DT + F T (2.14) ρ = ρ0 + ρ 0 (x, y, z, t) (2.15) où l’indice h représente les composantes horizontales des vecteurs, ρ0 est la densité de référence de l’eau de mer, f est le paramètre de Coriolis tel que f = 2Ω~ . ~k. Les termes D sont les paramètres des processus sous-mailles développés dans la section 2.1.4. 

Expression de la surface libre 

Suivant l’hypothèse hydrostatique (Eq.2.8), la pression est estimée comme étant la masse de la colonne d’eau présente au-dessus. La variation de la hauteur de mer doit être prise en compte dans le calcul du gradient de pression. Cette hauteur d’eau peut être déterminée de deux façons : (i) en introduisant l’élévation de la surface libre, notée η, ce qui permettra de calculer le gradient de pression de manière pronostique, (ii) en considérant l’interface entre l’air et l’océan comme étant un toit rigide, entraînant une détermination diagnostique du gradient de pression. Dans l’approximation de « toit rigide », on considère la surface de l’océan comme étant lisse, donc la surface est indexée en z=0, et qu’il n’y a pas de vitesses verticales en surface. Cette approximation est valable lorsque l’élévation de surface est négligeable devant la hauteur d’eau totale de l’océan. De moins en moins utilisée par la communauté, cette dernière formulation n’est plus présente dans le code NEMO à partir de la version 3.4. Dans la formulation de surface libre, l’élévation de surface est solution de l’équation pronostique : ∂η ∂t = −∇~ . h (H + η)U~ h i + P + R − E (2.16) avec H la profondeur de l’océan, P les précipitations, E l’évaporation et R l’apport d’eau terrestre. Les mouvements verticaux de la surface permettent la propagation d’ondes de gravité. Ces ondes sont barotropes avec une vitesse de phase élevée, ce qui implique un pas de temps de calcul très petit afin de préserver la condition CFL 7 (Courant et al., 1967) qui est ∆t < ∆x c . Afin de résoudre cette difficulté, différentes méthodes sont utilisées. La première, nommée timesplitting (Killworth et al., 1991; Zhang and Endoh, 1992), consiste à séparer la résolution des équations barotrope et barocline, avec un pas de temps barotrope plus petit que le barocline afin de résoudre la propagation des ondes de gravité. La part barotrope est ensuite moyennée sur le temps d’intégration barocline, avant d’être ajoutée à celle barocline. Une deuxième est l’utilisation d’un schéma implicite (Dukowicz and Smith, 1994) . Dans NEMO, une troisième méthode est utilisée, elle consiste à filtrer temporellement les ondes barotropes qui ne peuvent être résolues suivant la CFL par l’ajout d’une force supplémentaire −gTc∇  ρ˜ ∂η ∂t  qui agit comme un filtre passe-bas d’ordre 1 avec une période de coupure Tc (Roullet and Madec, 2000) : ∂U~ h ∂t = M~ − gT ∇~ (˜ρη) − gTc∇~  ρ˜ ∂η ∂t  (2.17) avec ρ˜ = ρ ρ0 la densité adimensionnée, et M~ qui représente les autres termes de l’équation du mouvement selon l’horizontale (Eq.2.10). La période de coupure, Tc est fixée à deux fois la valeur du pas de temps barocline afin de garantir une stabilité optimale.