Modélisation des courants marins

Pour toutes les équations, nous travaillerons en coordonnées cartésiennes, avec $\vec{e_x}$ et $\vec{e_y}$ les axes horizontaux et $\vec{e_z} $ l'axe vertical.

Les hypothèses

Nous allons prendre en compte plusieurs hypothèses, afin de simplifier les équations.

Approximation de Boussinesq

On considère que les variations de la masse volumique en fonction du temps et des variables spatiales  sont très faibles, soit :

$ \rho = \rho_0 + \rho' (x,y,z,t) $ avec  $\rho' << \rho_0 $

Cette approximation nous permet d'obtenir une divergence de la vitesse nulle : $ \nabla \mathbf U = 0 $.

Approximation hydrostatique

On considère que les variations verticales sont négligeables par rapport aux variations horizontales (facteur 1000). Les termes de Coriolis, d'advection et d'accélération ainsi négligeables devant le gradient de pression, ce qui nous permet d'obtenir l'équilibre hydrostatique : \[ \frac {\delta P}{\delta z}= - \rho g \]

Incompressibilité 

La masse volumique de l'eau ne dépend pas de la pression.

Hypothèse sur la profondeur 

On suppose également que la profondeur de l'océan est négligeable par rapport au rayon de la Terre, on peut donc se placer en coordonnées cartésiennes et non en coordonnées sphériques, ce qui simplifie grandement les équations.

 

 

 

Les équations régissant les courants

Les logiciels que nous allons utiliser pour les modélisations se servent des équations mises en place pour décrire les courants marins. Ces équations, dites primitives, sont obtenues après simplifications grâce aux hypothèses et donnent une description des courants suffisantes pour permettre de les modéliser. Elles permettent d'établir la circulation océanique générale, mais ne prennent pas compte les comportements plus locaux, difficilement prévisibles et modélisables, car ayant un caractère aléatoire.

Différents paramètres agissent sur l'océan et entraînent la formation des courants marins. Il y a deux types de forces influant sur les courants. D'une part, les forces actives, qui sont les forces productrices de courant. Elles comprennent les forces internes à l'océan qui sont dues aux variations de densité provenant des échanges énergétiques à l'interface entre l'eau et l'atmosphère et les forces externes dues aux vents, à la variation de la pression atmosphériques, ...

D'autre part, nous avons les forces passives, qui sont modificatrices du mouvement. Ce sont la force de Coriolis, la force de frottement due aux vents, ...

Tous ces paramètres permettent d'établir les équations primitives de circulation océanique, qui sont composées des équations de Navier-Stokes, de continuité, de conservation de la température et de la salinité ainsi que l'équation d'état de l'eau de mer.

Équations de conservation de la quantité de mouvement 

\[\frac{\delta u }{\delta t } + \mathbf {U} \nabla u = - \frac{1}{\rho_0} \frac{\delta P}{\delta x}+ f v   +A_x \frac{\delta^2 u}{\delta^2 x } + A_y \frac{\delta^2 u}{\delta^2 y } + A_z \frac{\delta^2 u}{\delta^2 z } \]

\[\frac{\delta v }{\delta t } + \mathbf {U} \nabla v= - \frac{1}{\rho_0} \frac{\delta P}{\delta y} - f u   +A_x \frac{\delta^2 v}{\delta^2 x } + A_y \frac{\delta^2 v}{\delta^2 y }+ A_z \frac{\delta^2 v}{\delta^2 z } \]

Dans ces équations, on voit l'influence de la pression, de la force de Coriolis et de la turbulence.
Pour la pression, on peut remarquer que le fluide se déplace de pressions hautes vers les basses pression (gradient négatif).

La force de Coriolis est la force due à la rotation de la Terre. En effet, en tournant, cette dernière entraîne une déviation de la trajectoire initiale. Dans l'hémisphère Nord, cette variation se fait dans le sens anti-trigonométrique.

Équation de continuité 

\[ \nabla \mathbf{U} = \frac{\delta u }{\delta x }+\frac{\delta v }{\delta y }+\frac{\delta w }{\delta z } = 0 \]

Équation hydrostatique 

\[dP = - \rho g dz \]

Équation de conservation de la salinité 

\[ \frac {\delta S}{ \delta t } + \mathbf{U}  \nabla S = K_{S,h} \nabla_h^2 S + K_{S,z} \frac{ \delta^2 S}{\delta z^2} \]

Équation de conservation de la température 

\[\frac{\delta T }{\delta t } + \mathbf{U} \nabla T = K_{T,h} \nabla_h^2 T + K_{T,z} \frac{ \delta^2 T}{\delta z^2} \]

Équation internationale d'état de l'eau de mer (IES80) 

\[ \rho = \rho (T,S,P) \]

Cette équation a été établie en 1980. Elle comporte 11 polynômes auxquels sont associés 41 coefficients numériques.

 

Les variables utilisées 

- $\mathbf {U} $: vitesse non turbulente du fluide de composantes u, v, w suivant les axes x, y , z. (m.s-1)

- $\mathbf {U'} $ : vitesse turbulente de composantes u', v', w' (m.s-1)

on a les relations suivantes entre la vitesse non turbulente et la vitesse turbulente :

$ \overline{u' u'} = - A_x \frac {\delta u }{\delta x} $          $ \overline{u' v'} = - A_y \frac {\delta u }{\delta y} $          $ \overline{u' w'} = - A_z \frac {\delta u }{\delta z} $

$ \overline{v' u'} = - A_x \frac {\delta v }{\delta x} $           $ \overline{v' v'} = - A_y \frac {\delta v }{\delta y} $          $ \overline{v' w'} = - A_z \frac {\delta u }{\delta z} $

$ \overline{w' u'} = - A_x \frac {\delta w }{\delta x} $          $ \overline{w' v'}= - A_y \frac {\delta w }{\delta y} $          $ \overline{w' w'} = - A_z \frac {\delta w }{\delta z} $

- $ \mathbf {A} $ : coefficient de viscosité turbulente ou d'Austaucht (m².s-1) de composantes Ax, Ay, Az.

- t : temps (s)

- $\rho_0 $ : masse volumique ( kg.m-3)

- P : pression (Pa)

- $f = 2 \Omega sin \theta $ : paramètre de Coriolis avec $\Omega$ le vecteur rotation de la Terre et $\theta $ la latitude.

- S : salinité (g.kg-1)

- T : température (K)

- $ \mathbf {K_S} $ : diffusivité cinématique (m².s-1) du sel de composantes KS,x, KS,y, KS,z. On a KS,x = KS,y = KS,h.

- $ \mathbf {K_T} $ : diffusivité cinématique (m².s-1) du sel de composantes KT,x, KT,y, KT,z. On a KT,x = KT,y = KT,h.

Les conditions limites

Conditions limites en surface

Aux équations primaires,  s'ajoutent les apports des éléments extérieurs en surface : le vent, le climat...influent également sur les courants. En effet, à la surface de l'océan, il y a des échanges permanents entre l'eau et l'atmosphère. Ces échanges sont de plusieurs types :

  • Échanges de moment, par les frottements dus à la différence de vitesse entre le vent et les courants en surface.
  • Échanges de chaleur, dus à la convection du vent et à la propre émission thermique de l'océan.
  • Échanges de masse, par la pluie et l'évaporation.

Ces influences sont en partie prises en compte grâce aux conditions limites en surface.

A la surface, la hauteur d'eau n'est pas constante, on définie $\eta$ tel que à la surface $ z= \eta(x,y,t) $.

Jusqu'à une certaine profondeur, les vents vont avoir une influence sur les courants. on va alors avoir ce qu'on appelle les transports d'Ekman. Expliquons rapidement le principe de ces transports : à la surface, les vents entraînent le fluide dans leur direction, cependant cette influence diminue rapidement avec la profondeur et  une autre force intervient, la force de Coriolis. Elle entraîne une déviation de la trajectoire des courants. Dans la couche limite ou les vents influencent les courants, ces derniers se retrouvent donc déviés par rapport à la direction du vent d'en moyenne 45°.

source : oceanmotion.org

L'influence du vent sur la surface de l'océan est décrit par les équations suivantes :

\[ A_z \frac{\delta u}{\delta z}=\ \frac{\tau_{s,x}}{\rho_0} \]

\[ A_z \frac {\delta v}{\delta z }=\ \frac{\tau_{s,y}}{\rho_0} \]

et \[ \frac {\delta \eta}{ \delta z} =w \]

avec $ \mathbf {\tau_s} = ( \tau_{s,x},\tau_{s,y})  $ la tension de suface en Pa. on a $ \mathbf {\tau_s} = \rho_{air} C_v || \mathbf V_{vent} || \mathbf V_{vent} $ avec :

  • $\rho_{air}$ la masse volumique de l'air
  • Cv le coefficient de trainée
  • $ \mathbf V_{vent} $ le vitesse du vent à 10 mètres d'altitude.

L'équation de conservation de la quantité de mouvement s'écrira alors en régime stationnaire:

\[ 0 = - \frac{1}{\rho_0} \frac{\delta P}{\delta x}+ f v  + \frac{1}{\rho}_0 \frac{\delta \tau_{s,x}}{\delta z} \]

\[ 0 = - \frac{1}{\rho_0} \frac{\delta P}{\delta y}- f u  + \frac{1}{\rho_0} \frac{\delta  \tau_{s,y}}{\delta z} \]

En résolvant ces équations, on obtient une formulation de la vitesse en surface, qui nous permet de définir l'épaisseur d'Ekman telle que : \[ D_E= \sqrt{ \frac {2 A_z}{f}}\]

Au dela de $\pi$ fois cette épaisseur, on considère que les vents n'ont plus d'influence sur les courants.

Pour la température 

\[ K_{T,z} \frac {\delta T}{\delta z}|_{z=\eta} = \frac {Q}{\rho_0 C_p}\]  avec Q le flux de chaleur de l'océan vers l'atmosphère en W.m-2 et Cp la capacité calorifique en J.kg-1.

Pour la salinité 

\[ K_{S,z} \frac {\delta S}{\delta z}|_{z=\eta}= \frac { S(E-P)}{\rho_0}\]

avec E l'évaporation et P la pluie.

 

Conditions limites en profondeur

Considérons une profondeur h. Au fond, z=-h. Ici, les conditions limites vont être dues aux frottements du fluide contre la paroi du fond.

On a les équations suivantes :

\[ A_z \frac {\delta u}{\delta z }= \frac{1}{\rho_0}\tau_{b,x} \]

\[ A_z \frac {\delta v}{\delta z }= \frac{1}{\rho_0}\tau_{b,y} \]

et \[ w= -u \mathbf { \nabla H} \]

 avec la contrainte pariétale en Pa $\mathbf{\tau_{b}} = \rho_0 C_d || \mathbf {V_b}|| \mathbf {V_b} $ avec Cd le coefficient de trainée et Vb la vitesse horizontale proche de la frontière.

Nous pouvons à nouveau obtenir des formulations de la vitesse dans une couche limite.

Pour la température 

\[ K_{T,z} \frac {\delta T}{\delta z}|_{z=-h} = 0 \]

Pour la salinité

\[ K_{S,z} \frac {\delta S}{\delta z}|_{z=-h} = 0 \]

Présentation du logiciel ROMS_AGRIF

Le logiciel

Le logiciel ROMS (Regional Oceanic Modeling System) est un outil de modélisation de circulation océanique. Cet outil a été développé par l'universioté Rutgers et l'université de Californie à Los Angeles ainsi que d'autres contributaires à travers le monde. Au sein de cette étude, nous allons utiliser ROMS_AGRIF, qui n'est autre que la version française du code ROMS, mise au point par l'IRD (Institut de Recherche pour le Développement) et l'INRIA (Institut National de Recherche en Informatique et en Automatique).

ROMS_AGRIF se base sur les équations que nous avons établi précédemment, en ajoutant l'hypothèse de turbulence isotrope horizontale soit Ax = Ay =Ah avec Ax le coefficient de viscosité turbulente suivant l'axe x et Ay le coefficient de viscosité turbulente suivant l'axe y.

Les équations de conservation de la quantité de mouvement deviennent donc :

\[\frac{\delta u }{\delta t } + \mathbf {U} \nabla u = - \frac{1}{\rho_0} \frac{\delta P}{\delta x}+ f u   +A_h \nabla_h^2 u + A_z \frac{\delta^2 u}{\delta^2 z } \]

\[\frac{\delta v }{\delta t } + \mathbf {U} \nabla v= - \frac{1}{\rho_0} \frac{\delta P}{\delta y} + f v  +A_h \nabla_h^2 v + A_z \frac{\delta^2 v}{\delta^2 z } \]

Le logiciel utilise les autres équations telles que nous les avons décrites précédemment.

Il prend également en compte les forçages en surface suivant :

  • Bilan de sel qui est équivalent à la soustraction évaporation - précipitation à laquelle on ajoute l'apport d'eau douce par le biais des fleuve.
  • Radiation solaire.
  • Tension du vent en surface.​

​​

Le maillage​​

Le modèle utilise la méthode des volumes finis pour résoudre les équations. Le maillage horizontal est basé sur une grille décalée d'Arawaka C, ou les flux sont décalés par rapport aux centres.

Le maillage vertical est fait de telle sorte que l'on ai toujours le même nombre de maille quelque soit la profondeur de l'océan. C'est une grille dont les niveaux suivent la topographie,la taille des mailles varient donc suivant la bathymétrie.

Source : Grids in Numerical Weather and Climate Models (intechopen.com/books/climate-change-and-regional-local-responses/grids-in-numerical-weather-and-climate-models)

 

 

Modélisation

Notre étude se fera au niveau du gyre Atlantique Nord. Nous avons donc choisi la zone d'étude de coordonnées : E = -85° , W =+11°; N= +50°,  S= +10°, afin d'avoir des modélisations du gyre entier.

Zone d'étude et sa bathymétrie

Procédure

Pour obtenir une modélisation des courants avec ROMS_AGRIF, nous devons suivre une procédure que nous allons résumer ici.

La première partie de la procédure se fait sous Matlab. Il faut tout d'abord régler les paramètres correspondant à notre zone d'étude et à la modélisation que l'on désire : les coordonnées de la zone, le maillage, le temps de modélisation, les forçages, les frontières, ainsi que les conditions initiales et limites. Nous pouvons indiquer la valeur de tous ces paramètres dans le fichier romstools_param.

Nous créons ensuite le maillage (make_grid) puis faisons intervenir les équations primaires et les forçages (make_forcing, make_clim).

Il nous faut ensuite compiler puis exécuter le programme.

Enfin, on peut visualiser les résultats à l'aide du logiciel roms_gui (sous Matlab).

 

Critère de stabilité

Avant de compiler, on va vérifier si le critère de stabilité CFL (courant Friedrichs-Lewy) est bien respecté. Ce dernier est défini ainsi :

\[ CFL =  C \Delta t_{externe} \sqrt{\frac{1}{dx^2}+\frac{1}{dy^2}} \]

Il doit être inférieur à 1 pour respecter la stabilité, ce qui nous donne :

\[ \Delta t_{externe} \leqslant \frac{1}{C} [ \frac {1}{ dx ^2} + \frac{1}{dy^2} ]^{\frac{-1}{2} } \]

avec $ C = \sqrt{g h_{max} }$ la vitesse des ondes de gravité en m.s-1 (g la gravité et $ h_{max} $ la profondeur maximum), dx et dy les pas suivant l'axe x et y en m.

Pour le pas de temps interne :

\[ 30 \leqslant \frac{\Delta t_{interne}}{\Delta t_{externe}} \leqslant 80 \]

 

Maillage

Pour nos modélisations nous avons choisi une résolution de 2/3 afin d'avoir un maillage assez fin tout en gardant des temps de calcul raisonnables.

 Nous choisissons également un nombre de mailles verticales : N=20.

Avec ces paramètres, nous obtenons un maillage avec L= 143 mailles suivant l'axe x et M= 71 mailles suivant l'axe y. Ce qui nous fait un total d'environ 2*105 mailles.

En prenant CFL=1, on a $\Delta t_{externe} = 80 s $.

On prend $\frac{\Delta t_{interne}}{\Delta t_{externe}} = 60 $, donc $ \Delta t_{interne} = 4800s$.

 

Temps de modélisation

Nous choisissons d'effectuer d'abord d'effectuer une modélisation sur un an.

Cependant, au centre  du gyre, il faut entre 4 et 6 ans pour faire un tour complet. D'autre part, il faut quelques années au logiciel pour se stabiliser. Pour une bonne représentation du trajet des plastiques, nous allons donc réaliser une modélisation plus longue, sur 10 ans.

Résultats

A la fin de la modélisation, nous obtenons des résultats sur de nombreuses variables. Nous nous intéresserons d'abord aux paramètres de l'océan qui influencent les courants, c'est-à-dire :

  • La salinité
  • La température de l'eau
  • La masse volumique.

En deuxième partie, nous analyserons les résultats obtenues pour les vecteurs vitesse.

Première modélisation : 1 an

Nous allons d'abord étudier les différents paramètres qui influent sur les équations des courants. Pour cela nous avons décidé de prendre leurs modélisations à deux périodes bien différentes de l'année : en hiver (décembre) et en été (juillet). 
 

La salinité

  

Image 1 : Salinité des mois de julllet (à gauche) et de décembre (à droite) en g de sel par kg de liquide

Nous pouvons observer que la salinité est plus forte au centre du gyre. Le long des grands courants (Gulf Stream) elle est plus basse. En effet, comme le la gyre est subtropicale, la température de l'air est plus élevée, ce qui entraîne une forte évaporation et donc une augmentation de la salinité.

 

Température de l'eau

Image 2 : Température à la surface de l'eau pour les mois de juillet (à gauche) et de décembre (à droite) en degré

​Plus on se dirige vers le sud, plus la température est élevée, ce qui est logique vu qu'on se rapproche de l'équateur. De plus, la température est plus élevée en été qu'en hiver . On peut également apercevoir les courants : au niveau du Gulf Stream, courant chaud, par exemple, la température est plus élevée.

 

La masse volumique

​​

​Image 4 : masse volumique de l'eau en kg/m​3 en juillet (à gauche) et en décembre (à droite).

La masse volumique dépend de la température de l'eau, de sa pression et de sa salinité. Lorsque la température augmente, la masse volumique diminue. Au sud de la carte, donc proche de l'Equateur, la masse volumique est plus faible que plus au Nord. C'est également le cas le long des côtes américaines, révélant bien la présence du Gulf Stream, courant chaud.

 

A présent, nous allons nous intéresser aux vecteurs vitesses.

Image 5 : Vitesse des courants suivant l'axe horizontale x en  juillet (à gauche) et en décembre (à droite).

Cette modélisation nous permet de visualiser les courants qui constituent le gyre Atlantique Nord. Tout d'abord, à l'Ouest de la carte, nous observons un courant fort passant par la pointe de la Floride et remontant le long des côtes américaines avant de se diriger vers le Nord-Est et donc vers l'Europe. Il s'agit du Gulf Stream. Ensuite, au Sud, se dessine un courant qui traverse l'océan Atlantique d'Est en Ouest. C'est le courant Nord Équatorial.

Cependant, les deux autres courants du gyre, à savoir le courant des Canaries et la dérive Nord Atlantique, sont moins bien représentés. Le courant des Canaries longe les côtes marocaines du Nord vers le Sud. Sur la carte, on voit que les vecteurs vitesse vont dans cette direction mais ils ne sont pas aussi marqués que pour le Gulf Stream par exemple. Cela est en partie dû au fait que le courant des Canaries est moins fort. En ce qui concerne la dérive Atlantique Nord, prolongement du Gulf Stream, elle est également peu visible car d'une part la zone d'étude choisie ne permet de la visualiser entièrement et d'autre part pour les mêmes raisons que le courant des Canaries.

Nous obtenons dans l'ensemble une bonne visualisation qualitative des courants du gyre, qui va nous permettre par la suite d'expliquer les trajectoires des déchets plastiques.

Néanmoins, cette modélisation a été faite sur un an, temps insuffisant pour la stabiliser. Nous ne pouvons donc pas être surs de l'exactitude des informations qu'elle nous fournit. Nous allons donc réaliser une modélisation plus longue, sur 10 ans, afin de stabiliser le modèle.

Seconde modélisation : 10 ans

En ce qui concerne la salinité, la température et la masse volumique de l'eau, nous observons que les tendances sont identiques chaque année. Au bout de 10 ans, lorsque nous avons un modèle stable, les résultats concernant ces paramètres sont pratiquement identiques à ceux de la première année. Nos interprétations précédantes restent donc valables, ce qui était déjà prédictibles, lorsque l'on compare avec les données réelles observées.

Intéressons nous maintenant aux courants.

Modélisation sur 10 ans

Nous observons toujours les courants décrits précédemment, cependant en 10 ans il y a une nette évolution du modèle. Sur la vidéo, il semble se stabiliser vers la cinquième année. Comparons maintenant la première et la dernière année de la modélisation. 

Modélisation des mois de juillet pour l'année 1 (à gauche) et l'année 10 (à droite)

Tout d'abord, nous remarquons que les vitesses modélisées lors de la première année sont plus importantes que la dernière année mais l'on observe tout de même que les tendances sont identiques. Les grands courants sont les mêmes. La stabilisation du modèle joue donc un rôle essentiellement sur la valeur des champs de vitesse et le modèle nous fournit donc dès la première année une représentation qualitative tout à fait correcte. 

Toutefois, pour obtenir une meilleure modélisation, il faudrait augmenter la précision des calculs (maillage plus fin, raffinage du maillage aux niveaux des courants principaux, ...)

De plus, cette représentation n'est que qualitative. Elle nous donne uniquement une idée des tendances des courants océaniques.