Le modèle numérique de circulation océanique régionale SYMPHONIE

Les équations du modèle

Les équations du mouvement et de continuité permettent de calculer les trois composantes du courant :où ρ0 est la masse volumique de l’eau de mer de référence, f est le paramètre de Coriolis, p est la pression, Kh est le coefficient de diffusion horizontale, et Kv est le coefficient de viscosité turbulente verticale, calculé à partir du schéma de fermeture de la turbulence proposé par Gaspar et al. [1990], détaillé plus loin. La pression p est donnée par l’équation de l’équilibre hydrostatique :

Le schéma de fermeture de la turbulence

Le schéma de fermeture de la turbulence proposé par Gaspar et al. [1990] permet de calculer le coefficient de viscosité turbulente verticale Kv à partir de l’énergie cinétique turbulente E = 1 2 (u′2+v′2+w′2), où l’apostrophe désigne les composantes turbulentes de la vitesse : où Ck est une constante empirique et lk est la longueur de mélange, caractéristique de la dimension des tourbillons associés à l’écoulement turbulent. L’équation de l’énergie cinétique turbulente est déduite des équations de quantité de mouvement.
Les valeurs des constantes empiriques Ck et Cε, respectivement 0.15 et 0.5, ont été légèrement modifiées par rapport à celles déterminées par Gaspar et al. [1990] (0.1 et 0.7) dans le but de favoriser le mélange dans le cas où la colonne d’eau devient instable. Une valeur seuil Emin est appliquée à l’énergie cinétique turbulente. Les échelles de longueur turbulentes de mélange et de dissipation lk et lε sont calculées en fonction des échelles de longueur lu et ld qui représentent la conversion d’énergie turbulente en énergie potentielle respectivement vers le haut et vers le bas Bougeault and Lacarrère [1989].
Par ailleurs, des instabilités statiques sont susceptibles de se développer à certains endroits du domaine. Les processus convectifs qui rétablissent la stabilité dans l’océan ne sont pas représentés par le modèle numérique en raison de l’hypothèse hydrostatique. On a donc recours à un algorithme d’ajustement convectif non pénétratif pour paramétriser ces processus [Andrich, 1998], suivant la méthode adoptée et validée par Madec [1990] : dans le cas où la colonne d’eau devient instable, on part de la surface (z = 0), et on calcule la moyenne entre la densité du point de grille et celle du point de grille inférieur. On continue en descendant jusqu’à ce que la colonne soit stable.

Les conditions aux limites

Pour résoudre les équations présentées plus haut, il est nécessaire de préciser les conditions aux limites, à savoir la surface, le fond, les embouchures de fleuves et les frontières latérales.

A la surface libre

Pour les composantes horizontales de la vitesse et l’énergie cinétique turbulente, les conditions à la surface libre s’expriment en fonction de la tension du vent −→τs = (τsx, τsy) :

A l’embouchure des fleuves

A l’embouchure du fleuve, la salinité est nulle, la température varie tout au long de l’année, et une vitesse horizontale u est appliquée dans la direction horizontale la plus proche de l’axe du fleuve, vérifiant :

La discrétisation des équations

Les valeurs utilisées dans les études présentées ici pour les différents paramètres et constantes intervenant dans les équations précédentes sont détaillées dans le Tab. La formulation mathématique complète des schémas numériques utilisés dans SYMPHONIE pour discrétiser les équations du modèle peut être trouvée dans Marsaleix et al. [2008].

La discrétisation spatiale

Les équations du modèle sont résolues par la méthode des différences finies sur une grille C [Arakawa and Suarez , 1983] présentée sur la Fig. 4.1, où sont représentées en haut la maille du mode barotrope et en bas la maille du mode barocline. Les variables baroclines sont définies un point sur deux sur l’horizontale et la verticale.
La température et la salinité sont définies au centre de la maille du mode barocline à chaque demi niveau. Les composantes horizontales de la vitesse sont obtenues au milieu des côtés, de façon alternée, à chaque demi niveau. L’énergie cinétique, les échelles de longueur turbulentes et la composante verticale de la vitesse sont calculées au centre de la maille du mode barocline à chaque niveau vertical entier. Les composantes du transport sont calculées au milieu des cotés de la maille barotrope, et l’élévation de la surface est calculée au centre de la maille du mode barotrope.
Par ailleurs, un système de coordonnées hybrides sigma-z généralisées est utilisé. La conversion de la coordonnée sigma en coordonnée z s’écrit : tel que σ vaut 1 en surface et 0 au fond. Ce système de coordonnées verticales permet une meilleure représentation des effets bathymétriques. Le système hybride permet de réduire le nombre de niveaux près des côtes, diminuant ainsi les coûts de calculs, sans sacrifier pour autant la résolution verticale au large. Il réduit également les erreurs de troncature associées au fort cisaillement vertical de la densité dans les zones de forte pente bathymétrique [Auclair et al., 2000b], par exemple le long du talus qui sépare le plateau continental et l’océan profond.

Le modèle numérique de circulation océanique régionale : SYMPHONIE

Le schéma d’advection est un schéma centré classique pour la vitesse [Arakawa and Suarez , 1983] et un schéma hybride centré/upstream pour la température et la salinité , adapté de Beckers [1995].

La discrétisation temporelle

Un schéma leapfrog ou saute-mouton explicite à l’ordre 2 est utilisé pour la discrétisation temporelle des équations. La valeur des variables au temps t + 1 est calculée en fonction des variables au temps t − 1 et t :

Aux frontières latérales

La composante normale de la vitesse à la frontière est nulle aux frontières latérales fermées. Aux frontières latérales ouvertes, les conditions aux limites ont un double objectif : la radiation des ondes sortantes, et le forçage de la solution intérieure au domaine par les champs externes, fournis dans nos simulations par les résultats d’un modèle de plus grande échelle. Comme suggéré par Blayo and Debreu [2005] et Marsaleix et al. [2006], ces objectifs sont réalisés en appliquant les schémas de conditions aux limites ouvertes à la différence entre les variables modélisées et les variables forçantes plutôt qu’à la variable elle-même. Considérons les frontières ouvertes en x = 0 et x = xm. On distingue le courant barotrope (u, v) et le courant barocline (u, v). Pour le courant barotrope, on applique une condition de type Flather [Flather , 1976] à l’élévation de surface :

La séparation des pas de temps

Le modèle calcule explicitement les ondes de gravité de surface. Les ondes de gravité externes se propagent plus rapidement que les ondes de gravité internes.
Une résolution temporelle plus fine est donc requise pour les résoudre et assurer la stabilité numérique du modèle. Une technique de séparation des pas de temps [Blumberg and Mellor , 1987] est utilisée afin de calculer séparément le cisaillement vertical de courant et le courant moyenné sur la verticale avec des pas de temps appropriés.
Le mode externe se calcule à partir de l’équation 4.2 et des équations du mouvement intégrées sur la verticale :

Implémentation du modèle SYMPHONIE enMéditerranée Nord-Occidentale

Dans toutes les simulations océaniques présentées dans la suite de ce document, nous utilisons le modèle SYMPHONIE avec une grille horizontale de 260 par 155 points, 40 niveaux sur la verticale, et une résolution spatiale de 3 km. Le rayon de déformation de Rossby étant égal à ∼ 10 km dans la région étudiée, cette résolution nous permet de représenter correctement les processus de méso-échelle, comme nous le verrons dans le chapitre 7. Le pas de temps du mode interne est de 225 s et celui du mode externe de 5.625 s. Le Rhône est la seule rivière prise en compte, avec un débit mensuel obtenu à partir de la base de donnée RivDis de l’UNESCO [Vörösmarty et al., 1996]. Comme nous le verrons dans la suite, les champs utilisés pour les conditions aux limites latérales ainsi que pour les conditions initiales sont obtenus à partir des résultats de simulations effectuées sur toute la Méditerranée au moyen du modèle OPAMED8 [Somot, 2005; Somot et al., 2006]. Ces champs sont interpolés sur la grille de SYMPHONIE au moyen de l’outil de méthode variationnelle inverse VIFOP présenté par [Auclair et al., 2000a, 2006]. Les flux atmosphériques sont obtenus à partir de différents modèles, que nous préciserons dans chaque chapitre.
La bathymétrie de la zone modélisée est représentée sur la Fig. 4.2. Cette bathymétrie est obtenue à partir de la base de données de Smith and Sandwell [1997], et a été raffinée aux endroits où la topographie est fortement accidentée (le plateau et le talus du Golfe du Lion et le Canal Corse) en utilisant les points de sonde indiqués sur les cartes topographiques du SHOM. Les frontières de la zone représentée sur la Fig. 4.2 ont été choisies de manière à ce que les processus physiques qui régissent le fonctionnement de la circulation et des écosystèmes en Méditerranée nord-occidentale, décrits dans le chapitre 2, soient contenus dans cette région : la convection profonde au large, le cascading d’eau dense le long du talus du Golfe du Lion, le courant Nord et plus largement la circulation cyclonique caractéristique de la zone.