<figure><img src="image_1.png" /></figure>

HAL
open science

# Modelisation numerique de phenomenes d'upwelling cotiers par une methode d'elements finis non conformes

Lien Hua, François Thomasset

► To cite this version:

Lien Hua, François Thomasset. Modelisation numerique de phenomenes d'upwelling cotiers par une methode d'elements finis non conformes. [Rapport de recherche] IRIA-RR-366, IRIA. 1979, pp.48. hal-04716624

HAL Id: hal-04716624

https://inria.hal.science/hal-04716624v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Inv: 07382 RR 366

IRIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

<figure><img src="image_3.png" /></figure>

MODÉLISATION NUMÉRIQUE
DE PHÉNOMÈNES
D'"UPWELLING" CÔTIERS
PAR UNE MÉTHODE
D'ÉLÉMENTS FINIS
NON CONFORMES

Lien HUA
François THOMASSET

RR 366
Rapport de Recherche N° 366

Octobre 1979
MODELISATION NUMERIQUE DE PHENOMÈNES D' "UPWELLING" COTIERS
PAR UNE METHODE D'ELEMENTS FINIS NON CONFORMES.

Lien Hua* - François THOMASSET

Résumé :

Le terme d' "upwelling" désigne en Océanographie la remontée d'eaux froides plus profondes vers la surface. On s'intéresse ici à l'upwelling côtier induit par le vent.

On utilise un modèle d'océan à deux couches, chacune de densité constante ; les équations sont du type "eaux peu profondes" (équations de Navier Stokes intégrées sur la verticale dans chaque couche).

La méthode d'éléments finis présentés ici utilise des éléments linéaires non conformes et permet de prendre en compte aisément la ligne de côte et la topographie du Golfe du Lion ainsi qu'un champ de vent arbitraire. Un schéma semi implicite pour les termes de gravité (schéma de KWIZAK et ROBERT) permet de prendre des pas de temps de l'ordre de 5 minutes à 1 heure, et donc d'avoir des temps de calcul raisonnables. La méthode a été testée sur un cas académique et montre une bonne concordance avec la solution analytique. Finalement, nous présentons des résultats sur le Golfe du Lion.

Abstract :

In Ocean Dynamics "upwelling" is defined as an ascending motion by which deep, cool water is brought to the surface. We are concerned here with coastal upwelling induced by winds. We use a two layer ocean model ; in each layer shallow water equations are discretized with piecewise linear, non conforming finite elements. Thus we are able to deal easily with the coastal line and topography of the Gulf of Lions and with an arbitrary wind stress. A semi implicit time scheme (with respect to gravity termes : KWIZAK and ROBERT's scheme) allows time steps of order 5 minutes to 1 hour and therefore reasonnable CPU times. The method has been tested on an academic problem and shows good agreement with the analytic solution. Finally we show some results on the Gulf of Lion.

<footnote>*Laboratoire d'Océanographie Physique du Muséum d'Histoire Naturelle, Paris.</footnote>
# INTRODUCTION

## I. - HYPOTHESES DU MODELE

A). Equations du modèle

B). Conditions aux limites et initiales

C). Ordres de grandeur

D). Modes propres du système

## II. - METHODE NUMERIQUE

A). Formulation faible des équations

B). Choix de la méthode de discrétisation

C). Schéma de la discrétisation en temps

1. Schéma explicite

2. Stabilité du schéma

D). Schéma semi-implicite

1. Traitement implicite des ondes de gravité

2. Structure matricielle des systèmes obtenus

E). Avantages des éléments finis linéaires non conformes

F). "Coût" de la méthode.

## III - MODELE ACADEMIQUE

A). "Contexte physique"

B). Résultats numériques

## IV. - GOLFE DU LION

A). Domaine

B). Champ du vent

C). Résultats numériques

# CONCLUSION

# APPENDICE

# BIBLIOGRAPHIE
<figure><img src="image_4.png" /><figcaption>Figure 1</figcaption></figure>

Vue moyennée sur plusieurs coups de vent, des isolignes de la température de surface prise par le satellite NOAA IV.

Les pointillés désignent la présence de points-sources d'upwelling intenses.

Pour situer géographiquement la région, on a noté 3 points géographiques :

A : Marseille B = Pointe de l'Espiguette C = Cap d'Agde

Millot (1979)
# INTRODUCTION

Le terme d'upwelling désigne toute remontée d'eau de couches plus profondes, généralement plus froides (pour une stratification stable) et plus riches en éléments nutritifs, vers la surface.

L'origine de ces phénomènes repose dans l'existence d'une divergence horizontale du transport dans la couche superficielle, divergence qui doit être compensée par un mouvement vertical à cause de la continuité de la masse. (1)

Dans le cas de l'upwelling côtier causé par la composante parallèle à la côte, du vent, un schéma du phénomène serait le suivant : le transport d'Ekman (2) vers le large, causé par le vent en surface doit être compensé par un mouvement vertical d'eaux plus profondes au niveau de la côte, ce qui induit un afflux vers la côte des eaux de subsurface provenant du large. Dans la réalité, le phénomène est plus complexe : il peut devenir tridimensionnel aux points où le mélange est important et surtout il est très fortement transitoire. Du point de vue de la dynamique des océans, ces processus côtiers ne concernent que des régions relativement restreintes et conduit sur des échelles d'espace et de temps relativement petites (d'un point de vus géophysique) à d'importants gradients des paramètres physiques : contrastes très marqués de température et grande variabilité des champs de vitesses. C'est donc un phénomène relativement énergétique à petite échelle ("small scale-high magnitude process") qui occupe une part non négligeable dans le spectre d'énergie de la circulation océanique et qu'il est important de modéliser. Ces processus aboutissent à un enrichissement en éléments nutritifs des eaux plus oxygénées de la zone euphotique et il en résulte un fort accroissement de la productivité biologique. Leur importance écologique et économique est donc considérable pour la région du plateau continental qui est une des zones les plus fertiles de l'océan mondial.

<footnote>(1) Si $\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \neq 0$, alors $\frac{\partial w}{\partial z} \neq 0$ dans l'équation de continuité.</footnote>

<footnote>(2) Mouvement dans la couche de surface qui est dirigé sur la droite du vent dans l'hémisphère Nord et qui résulte d'un équilibre entre le cisaillement vertical induit par le vent et la force de Coriolis.</footnote>
Les travaux théoriques et numériques antérieurs concernent des échelles dites "meso-scale" (échelles de l'ordre de 100 km et de plusieurs jours au mois), qui sont importantes pour les upwellings des bords Est des océans (direction favorable des vents) tels ceux de l'Orégon ou du Pérou (cf. O'Brien et Peffley [9]). Nous essayons ici de simuler numériquement des upwellings d'échelle plus petite (variabilité de l'ordre de 10 km et de quelques heures à quelques jours) comme ceux qui se produisent l'été dans le Golfe du Lion (Fig. 1).

Ils sont induits par des coups de tramontane et de mistral, vents de N-NW, qui peuvent se lever très brusquement, puis rester stables en intensité et direction pendant plusieurs jours et couvrir tout le golfe.

Les observations in situ (cf. Millot (1979)) jointes aux images satellites de la température de surface ("signature" la plus immédiate de l'upwelling) montrent l'existence de véritables cellules d'upwelling très précisément localisées et ceci de manière répétitive pour les divers coups de vent qui se produisent sur un été. Elles semblent fortement liées à la géométrie de la côte, ce qui suggère une méthode d'éléments finis pour discrétiser le problème.

L'exposé comprend quatre parties :

1) - On essaie de re-situer les hypothèses du modèle dans le contexte des travaux antérieurs et en particulier, les limitations spatiales et temporelles des équations du problème.

2) - Exposé de la méthode numérique, du choix d'éléments finis non conformes et de la méthode semi-implicite mise au point.

3) - Présentation d'un test de l'algorithme sur un cas académique où l'on possède une solution analytique pour le problème linéarisé transitoire.

4) - Résultats numériques de la simulation pour le Golfe du Lion.

-------
<figure><img src="image_5.png" /><figcaption>Figure 2. Profil de densité selon la verticale</figcaption></figure>

La densité est définie par $\rho = 10^3 (1 = \sigma_t \times 10^{-3}) \text{ kg/m}^3$.

<figure><img src="image_6.png" /><figcaption>Figure 3. Géométrie du modèle à deux couches.</figcaption></figure>

$D(x,y)$ est la hauteur du fond au-dessus d'un niveau de référence arbitraire.

(Remarque : l'interface entre les deux couches ne peut pas intersecter la pente du fond).
# I. - HYPOTHESES DU MODELE

Le modèle n'est pour l'instant que dynamique et ne prend pas en compte les processus thermodynamiques.

La structure verticale est celle d'un océan à 2 couches, chacune de densité constante avec un saut de densité au niveau de la thermocline$^{(1)}$. La pression est hydrostatique. Les vitesses dans chaque couche sont indépendantes de la profondeur. Ceci constitue une hypothèse assez forte mais qui semble raisonnable pour une première approche du problème :

- en effet, les mesures in situ de profil de température montrent la présence d'une thermocline bien marquée vers 20-30m selon la position du point par rapport à la côte et le profil de densité montre que l'hypothèse de deux couches de densité constante semble justifiée (cf. Figure 2).

- les observations de courantométrie montrent d'autre part, un découplage assez net entre les courants de part et d'autre de la thermocline : la couche superficielle semble "glisser" presque sans frottement sur l'autre, en un mouvement de "dalle" ("slab-motion"), où les vitesses semblent indépendantes de la profondeur. Les échelles de temps (de quelques heures à quelques jours) et d'espace (de quelques kilomètres), justifient ces hypothèses qui sont celles des ondes longues. L'essentiel de la dynamique est donc pris en compte puisque le modèle permet la reproduction du mode barotrope et du premier mode barocline$^{(2)}$; en particulier, on ne fait aucune hypothèse a priori sur les mouvements de la surface libre telle que l'hypothèse du toit rigide appliqué à la surface, ce qui permettrait les plus grands pas d'intégration en filtrant les ondes de gravité de surface$^{(3)}$, mais qui peut fausser la possibilité d'ajustement géostrophique$^{(4)}$ et filtre également toute la classe des mouvements globalement divergents sur la colonne d'eau.

Le couplage entre les mouvements des deux couches est essentiellement dû au terme de pression, bien que l'on inclue une légère tension d'entraînement entre les deux couches mais son rôle reste faible.

<footnote>(1) zone de fort gradient de température selon la verticale.</footnote>

<footnote>(2) une définition grossière serait que le mode barotrope concerne le mouvement global de la colonne d'eau alors que le premier mode barocline permet de prendre en compte une première dépendance verticale du profil des vitesses horizontales par le cisaillement possible entre les deux couches.</footnote>

<footnote>(3) voir plus loin.</footnote>

<footnote>(4) équilibre entre les forces de pression et les forces de Coriolis.</footnote>
Le seul "forcing" extérieur est dû au vent qui agit en surface par sa tension
(loi du $C_{10}$) et une loi de frottement quadratique est prise sur le fond
($C_B \sim C_{10}$).

La turbulence horizontale est paramétrisée par un terme de diffusion turbulente.
Les déplacements verticaux liés à l'upwelling à la côte sont donc représentés
par le relèvement de l'interface, ce qui revient à dire que les vitesses verti-
cales peuvent être obtenues à partir de la vitesse lagrangienne de l'interface.

En résumé, la thermocline constitue un écran à la fois aux transferts de masse,
de chaleur et de quantité de mouvement par entraînement turbulent. Ceci constitue
une forte limitation du modèle qui devra être levée dans la poursuite ultérieure
des travaux. Il en résulte en particulier qu'il faut arrêter l'intégration dès
que l'interface vient crever la surface libre. Dans la réalité, aux points où
l'épaisseur de la couche supérieure devient mince, les phénomènes de brassage
et de mélange interviennent pour l'approfondir de nouveau.

A) - Equations du Modèle

Les équations sont celles de Hurlburt et O'Brien [6] . Elles sont obtenues par
intégration des équations de Navier-Stokes sur chaque couche ; les orientations
des axes et la signification des différentes quantités sont données sur la Fig.3.

$$
\begin{align*}
\text{a) } \frac{\partial \vec{u}_i}{\partial t} + (\vec{u}_i \cdot \nabla) \vec{u}_i + f_k \times \vec{u}_i &= -g \nabla (h_1 + h_2 + D) - g \delta_{i2} \frac{\rho_2 - \rho_1}{\rho_2} \nabla h_1 \\
(1) \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \partial h_i \partial t \\
\text{b) } \frac{\partial \vec{h}_i}{\partial t} + \text{div} \left( (\vec{h}_i \cdot \nabla) \vec{u}_i \right) &= 0 \quad \text{i=1,2.}
\end{align*}
$$

avec

• saut de densité relatif $\frac{\rho_2 - \rho_1}{\rho_2} = 2.10^{-3}$

. tension en surface dûe au vent $\vec{\tau}_s = \rho_a C_{10} |\vec{w}| \vec{w} = 1.6 \ 10^{-3} |\vec{w}| \vec{w}$
($\vec{w}$ est le vent à 10 m au-dessus du niveau de la mer)

. tension à l'interface entre les deux couches

$$
\vec{\tau}_{B1} = \tau_{S2} \rho C_1 |\vec{u}_1 - \vec{u}_2| (\vec{u}_1 - \vec{u}_2) \quad (C_I = 10^{-5})
$$

. tension de frottement au fond

$$
\tau_{B2} = \rho_2 C_B |\vec{u}_2| \vec{u}_2 \quad (C_B = 10^{-3})
$$

<footnote>(1) L'étude de l'existence et de l'unicité d'une solution, ainsi que l'analyse numérique d'un modèle analogue, est faite dans M. Benahmed [13].</footnote>
paramètre de Coriolis : $f = 10^{-4} s^{-1}$

accélération de pesanteur : $g = 10 \text{ m.s}^{-2}$

coefficient de viscosité turbulente horizontale $v = 10^2 \text{ m}^2 \text{ s}^{-1}$

symbole de Kronecker : $\delta_{i2}$ avec $\delta_{12} = 0, \delta_{22} = 1$.

## B) - Conditions aux limites et initiales.

Les équations sont intégrées à partir de l'état du repos et le vent agit à partie de $t=0$ ($w=10 \text{ m/s}$) et reste ensuite localement constant en direction, ceci n'exclut pas une variation spatiale et temporelle de l'intensité.

La condition aux limites sur toutes les parois latérales (la côte et les frontières océaniques du domaine) est le *glissement sans frottement* ("free-slip") ; elle permet aux parois solides de ne pas faire rejaillir trop loin la couche-limite visqueuse, lorsque la taille de la maille ne permet pas de résoudre spatialement la couche-limite (qui serait de l'ordre du Km.)$^{(1)}$.

En raison de ce choix de condition aux limites, les frontières océaniques du domaine ont été repoussées le plus loin possible de la zone concernée par le "forcing" du vent, afin d'éviter le problème de fausses réflexions d'ondes aux frontières océaniques ouvertes qui peuvent perturber la solution dans la zone à laquelle on s'intéresse$^{(2)}$.

D'autre part, bien que la topographie puisse être prise en compte dans le modèle, on ne peut considérer que le cas où l'interface et la surface libre intersectent un mur vertical à la côte et non pas une frontière pouvant présenter une pente, sinon le domaine d'intégration serait variable dans le temps.

## C) - Ordres de grandeur.

Pendant le temps d'intégration du modèle (jusqu'à 2 ou 3 jours), les termes prépondérants dans les équations du mouvement sont l'accélération temporelle, celle de Coriolis, le terme d'entraînement dû au vent et les termes de pression.

La valeur du coefficient de frottement à l'interface ($C_I = 10^{-5}$) a été choisi

<footnote>$^{(1)}$ $\vec{u} \cdot \vec{n} = 0$ et $\frac{\partial}{\partial u}(u \cdot v) = 0$ à la côte.</footnote>

<footnote>$^{(2)}$ On pourrait envisager des conditions aux frontières océaniques de type "radiation" qui permettent aux ondes de sortir, mais non au courant.</footnote>

<footnote>$^{(3)}$ Il sera donc important de les intégrer *exactement* de manière numérique.</footnote>
de manière à conserver la cohérence de l'ensemble des hypothèses du modèle
avec l'hypothèse que les mélanges de type thermohalin n'ont pas le temps de
se produire (ils ne deviennent importants qu'au-delà de quelques jours).
Bien que les équations soient proches de celles d'Euler, la vitesse des ondes
les plus rapides (ondes de gravité) est de l'ordre de 100m/s et est nettement
supérieure aux vitesses maximales du fluide, qui sont de l'ordre de 50cm/s ;
la solution physique ne contient pas de choc.

D) - Modes propres du système.

Il est important de pouvoir reproduire correctement la "physique" du système,
que la discrétisation numérique soit cohérente avec les diverses oscillations
propres possibles du système. Dans le cas de l'upwelling côtier en présence
de topographie et de stratification et en milieu tournant, celles-ci sont
multiples :

- oscillations d'inertie de Coriolis,

- ondes de gravité longues externes (à la surface) et internes (à l'inter-
face entre les deux couches), modifiées par la rotation de la Terre.

Une sous-classe importante de ces ondes est l'onde de Kelvin qui ne peut exister
qu'en présence d'une paroi latérale et en milieu tournant : elle est caractérisée
par une décroissance exponentielle vers le large de la dénivellation et se
propage en laissant la côte sur la droite.

- ondes de plateau continental induites par l'existence d'une pente du fond.

Il faut notamment que la vitesse "numérique" du modèle qui est $\frac{\Delta x}{\Delta t}$ où $\Delta x$ est la taille caractéristique d'une maille et $\Delta t$ le pas de temps d'un schéma explicite, soit plus grande que la vitesse de propagation la plus rapide dans le fluide, c'est-à-dire celles des ondes de gravité$^{(1)}$ externes.

Nous rappelons ici une dérivation de l'équation de propagation des ondes de
gravité qui sera la base du schéma semi-implicite de discrétisation en temps
proposé plus loin (cf. point II-D) 1., p. 18).

<footnote>(1) Les ondes de gravité sont des oscillations d'une particule autour d'une
position d'équilibre statique, (elles sont le mieux visualisées par les
mouvements d'une surface libre) où la force de gravité joue le rôle de
force de rappel.
Ces oscillations induisent une variation locale de la *divergence
horizontale* du fluide (à cause de l'équation de continuité).</footnote>
Dans le cas de notre modèle à deux couches, en négligeant les effets de la rotation terrestre, elles résultent de l'équilibre des seuls termes suivants :

$$ (2) \quad \left\{ \begin{array}{l} \displaystyle \frac{\partial \vec{u}_i}{\partial t} = - g \text{ grad } (\alpha_i h_i' + h_2') \\[1em] \displaystyle \frac{\partial h_i'}{\partial t} + \text{div}(H_i \vec{u}_i) = 0 \end{array} \right. $$

avec

$$ h_i = H_i + h_i'(x,y,t) $$

$$ \alpha_i = 1 - \delta_{i2} \frac{\rho_2 - \rho_1}{\rho_2} $$

par dérivation croisée, on obtient l'équation de propagation :

$$ \frac{\partial^2}{\partial t^2} h_i' - g H_i \Delta (\alpha_i h_i' + h_2') = 0 . $$

Elles admettent deux modes de propagation dits externe et interne

$c_1 = (g(h_1 + h_2))^\frac{1}{2}$ correspondant à la célérité des ondes de gravité de surface,

$c_2 = (g \frac{h_1 h_2}{h_1 + h_2})^{\frac{1}{2}}$ correspondant à la célérité des ondes de gravité à l'interface.

$c_1 = 100 \text{ m/s}$, $c_2 = 1 \text{ m/s}$, pour les ordres de grandeur du problème.

Avec la maille choisie, ceci conduirait à des pas de temps inférieurs à la minute, pour des temps d'intégration en temps réel de quelques jours, ce qui est donc extrêmement pénalisant.

En traitant de manière *implicite* les seuls termes qui interviennent dans le mécanisme des ondes de gravité externes et internes, on pourra s'affranchir de cette contrainte. Les résultats analytiques montrent, d'autre part, qu'une échelle spatiale caractéristique des ondes de gravité est le *rayon de déformation* de Rosby défini par $\frac{c_i}{f}$, $i=1, 2$, où $c_i$ est la vitesse des ondes de gravité et où $f$ est le paramètre de Coriolis. Pour l'upwelling, ce sont surtout les ondes internes de gravité qui interviennent pour propager l'énergie *en dehors* de la zone de forcing du vent et permettent ainsi l'ajustement géostrophique local, c'est-à-dire un retour du système à un état d'équilibre stationnaire non divergent horizontalement propre aux systèmes en milieu tournant. Etant donné les valeurs des paramètres du système, le rayon interne de déformation est de l'ordre de 10 kms. et c'est donc l'échelle spatiale que la maille du domaine devra au moins résoudre.
## II. - METHODE NUMERIQUE

### A) - Formulation faible des équations (Méthode de Galerkin).

Ainsi qu'il est habituel pour la méthode des éléments finis, on remplace d'abord les équations (1) par une formulation variationnelle. Pour cela, on considère une fonction vectorielle $v(x,y)$ "assez régulière" et satisfaisant aux conditions aux limites homogènes

$$ (3) \qquad \vec{v}.n = 0 \text{ sur } \Gamma = \text{frontière du domaine.} $$

Multiplions par $\vec{v}$ et intégrons le résultat sur le domaine $\Omega$.
En utilisant les notations classiques, on a

$$ (\rho_1, \rho_2) = \int_{\Omega} \varphi_1(x,y) \varphi_2(x,y) dxdy, $$

$$ (\vec{v}_1, \vec{v}_2) = \int_{\Omega} v_1(x,y) \cdot \vec{v}_2(x,y) dxdy, $$

$$ (\operatorname{grad} \vec{v}_1, \operatorname{grad} \vec{v}_2) = \int_{\Omega} \left( \frac{\partial \vec{v}_1}{\partial x} \frac{\partial \vec{v}_2}{\partial x} + \frac{\partial \vec{v}_1}{\partial y} \frac{\partial \vec{v}_2}{\partial y} \right) dxdy, $$

A l'aide de la formule de Green et des conditions aux limites, on obtient :

$$ \left(\frac{\partial \vec{u}_i}{\partial t} + (\vec{u}_i \cdot \nabla) \vec{u}_i + g \operatorname{grad}(\alpha_i h_i + h_2 + D) + C_F \left(\frac{|\vec{u}_i|}{h_i} \frac{\partial \vec{u}_i}{\partial x} + f \vec{k} \cdot \vec{u}_i, \vec{v}\right) + \right. \\ \left. + v (\operatorname{grad} \vec{u}_i, \operatorname{grad} \vec{v}) = \delta_{il} \left(\frac{|\vec{w}|}{\rho_1 h_1}, \vec{v}\right) . \right) $$

De même, en multipliant (1-a) par une fonction scalaire $\varphi$ et en intégrant sur $\Omega$, on obtient :

$$ (5) \qquad \left(\frac{\partial h_i}{\partial t}, \varphi\right) - \left(h_i \cdot \vec{u}_i, \operatorname{grad} \varphi\right) = 0. $$

Les équations (4) et (5) doivent être satisfaites pour $i = 1$ et $2$, et pour toutes fonctions $\vec{v}, \varphi$ assez régulières.

La méthode des éléments finis consiste à prendre pour fonctions "test" $\vec{v}$ et $\varphi$ des fonctions moins régulières, suffisamment simples à calculer, polynômiales par morceaux.
**Remarque** : On peut montrer que pour des éléments finis non conformes, la formulation précédente permet de conserver *exactement* la masse :

$$ \int_{\Omega} \frac{\partial h_i}{\partial t} dxdy \equiv 0 \quad \text{si} \quad \vec{u}_i \cdot \vec{n}|_{\Gamma} = 0. $$

## B) - Choix de la méthode de discrétisation.

On a vu que les échelles spatiales qu'il faut résoudre (rayon interne de déformation) imposent à la maille spatiale d'être au plus égale à environ 3 km., ce qui conduit par exemple pour le cas du Golfe du Lion à 1.000-2.000 triangles. Par conséquent, il nous apparaît préférable d'utiliser les éléments finis les plus simples possibles quitte à utiliser un grand nombre de triangles. C'est pourquoi nous avons choisi de discrétiser $\vec{u}_i$ et $h_i$ par éléments finis linéaires non conformes décrits au paragraphe suivant.

### Les éléments finis linéaires non conformes.

Cet élément fini a été introduit par Crouzeix et Raviart [5] et utilisé par Thomasset [10],[11], pour la résolution des équations de Navier-Stokes en deux dimensions d'espace.

Pour ces équations, on montre que l'erreur sur la vitesse en norme $H_1$ est en $O(\delta x)$ où $\delta x$ est l'ordre de grandeur de la maille.

Décrivons maintenant la méthode :

Nous recouvrons le domaine étudié avec un nombre fini de triangles, avec les règles habituelles aux éléments finis :

- les triangles ne se recouvrent pas ;

- un sommet de triangle ne peut appartenir à une arête d'un autre triangle, à moins d'être aussi sommet de ce triangle.

La figure 10 montre la triangulation que nous avons utilisée pour le Golfe du Lion, soit 524 triangles ($\delta x$ = 10 km.).

Nous utiliserons les notations suivantes :

$\varphi$ = ensemble des milieux de côtés de tous les triangles,

p = un milieu de côté,

$w_p$ = la fonction de base associée au point p :
• linéaire affine sur chaque triangle,

• égale à 1 au point p

• égale à 0 en tous les autres milieux de côtés.

Sur le triangle (abc) de la figure ci-contre, $w_p$ prend les valeurs suivantes:

<figure><img src="image_7.png" /></figure>

$$w_p(a) = w_p(b) = 1$$

$$w_p(c) = -1$$

$$w_p(q) = w_p(r) = 0$$

Il est clair que sur l'ensemble du domaine, la fonction $w_p$ est *continue* aux milieux des côtés, mais *discontinue* aux autres points des arêtes : c'est le sens de l'expression "non conforme".

D'autre part, $w_p$ est non nulle seulement sur les deux triangles adjacents au noeud $p$, soit (abc) et (abc') sur la figure.

On notera

$$S_p = \frac{1}{3} \text{ (aire en support de } w_p) = \frac{1}{3} \text{ aire (abc)} + \frac{1}{3} \text{ aire (abc')}.$$

L'une des propriétés fondamentales est l'orthogonalité des fonctions de base

$$(6) \qquad (w_p, w_q) = \int_{\Omega} w_p(x,y) w_q(x,y) dx dy = \begin{cases} S_p & \text{si } p = q \\ 0 & \text{si } p \neq q \end{cases}.$$

On représente les inconnues $h_i$ et $u_i$ (dans la couche $i$, $i=1,2$) comme des combinaisons linéaires des fonctions de base $w_p$ :

$$(7) \qquad h_i(x,y) = \sum_{p \in S} h_{i,p} w_p(x,y)$$

$$(8) \qquad u_i(x,y) = \sum_{p \in S} \begin{pmatrix} u_{x,i,p} \\ u_{y,i,p} \end{pmatrix} w_p(x,y),$$

avec $$h_{i,p} = h_i(p)$$

$$\begin{pmatrix} u_{x,i,p} \\ u_{y,i,p} \end{pmatrix} = \begin{pmatrix} u_i(p) \end{pmatrix}$$
(on rappelle que $\mathfrak{S}$ est l'ensemble de tous les noeuds (milieux de côtés) du domaine).

Compte tenu des conditions aux limites (1) sur $\vec{u}_2$, il convient de distinguer entre les noeuds intérieurs et frontaliers.

On note donc :

$\mathfrak{S}_0$ = l'ensemble des noeuds intérieurs

$\mathfrak{S}_1$ = l'ensemble des noeuds frontaliers.

Si $p$ appartient à $\wp_1$, on pose :

$\vec{\tau}_p = \begin{pmatrix} t_{x,p} \\ t_{y,p} \end{pmatrix} = \text{vecteur unitaire tangent à la frontière,}$

$\vec{u}_i(p) = \begin{pmatrix} u_{x,i,p} \\ u_{y,i,p} \end{pmatrix} \equiv u_{i,p} \begin{pmatrix} t_{x,p} \\ t_{y,p} \end{pmatrix}$

c'est-à-dire $u_{i,p} = \vec{u}_i(p) \cdot \vec{\tau}_p$.

Finalement

$$ (9) \qquad \vec{u}_i(x,y) = \sum_{p \in \mathfrak{S}_0} \begin{pmatrix} u_{x,i,p} \\ u_{y,i,p} \end{pmatrix} w_p(x,y) + \sum_{p \in \mathfrak{S}_1} u_{i,p} \begin{pmatrix} t_{x,p} \\ t_{y,p} \end{pmatrix} \times w_p(x,y). $$

Les inconnues sont maintenant les nombres :

$u_{x,i,p}(t)$, $u_{y,i,p}(t)$, pour $p \in \mathfrak{S}_0$

$u_{i,p}(t)$, pour $p \in \mathfrak{S}_1$

$h_{i,p}(t)$, pour $p \in \mathfrak{S}$

Pour établir les équations du problème discret, il reste à remplacer dans (4) et (5) :

$\vec{u}_i$ et $h_i$ par (7) et (9), $\vec{v}$ par $\begin{pmatrix} w_p \\ 0 \end{pmatrix}$ ou $\begin{pmatrix} 0 \\ w_p \end{pmatrix}$, ($p \in \mathfrak{S}_0$) par $\vec{\tau}_p w_p$ ($p \in \mathfrak{S}_1$)

<footnote>(1) Il serait possible de considérer d'autres conditions aux limites, par exemple $\vec{u} \cdot \vec{n} \neq 0$.</footnote>
et $\varphi$ par $w_p(r \in S)$.

On obtient un système d'équations différentielles ordinaires (en utilisant la relation d'orthogonalité (6)):

$p \in S_0$

$$
(10) \quad S_p \frac{d}{dt} u_{x,i,p} + ((\vec{u}_i \cdot \vec{\nabla}) u_{x,i}) + g \frac{\partial}{\partial x} (\alpha_i h_1 + h_2 + D) - f u_{y,i} + C_F \frac{|u_i|^2 u_{x,i}}{h_i^2}, w_p) + \\
\qquad + v(\operatorname{grad} u_{x,i}, \operatorname{grad} w_p) = \delta_{i1} \frac{|w|^2 w_x}{h_1^2}, w_p)
$$

$$
(11) \quad S_p \frac{d}{dt} u_{x,i,p} + ((\vec{u}_i \cdot \vec{\nabla}) u_{x,i}) + g \frac{\partial}{\partial y} (\alpha_i h_1 + h_2 + D) + f u_{x,i} + C_F \frac{|u_i|^2 u_{y,i}}{h_i^2} w_p) + \\
\qquad + v(\operatorname{grad} u_{y,i}, \operatorname{grad} w_p) = \delta_{i2} \frac{|w|^2 w_y}{h_1^2} w_p.
$$

$p \in S_1$

$$
(12) \quad S_p \frac{d}{dt} u_{i,p} + ((\vec{u}_i \cdot \vec{\nabla}) \vec{u}_i + g \operatorname{grad}(\alpha_i h_1 + h_2) + D) + \vec{f}_k \times \vec{u}_i + C_F \frac{|u_i|^2 \vec{u}_i}{h_i^2}, \vec{\tau}_p w_p) + \\
+ v (\tau_{x,p} \operatorname{grad} u_{x,i} + \tau_{y,p} \operatorname{grad} u_{y,i} \operatorname{grad} w_p) = \\
= \delta_{i1} \left( \frac{|w|^2 \vec{w}_i}{\rho_i h_1}, \vec{\tau}_p w_p \right).
$$

$r \in S$

$$
(13) \quad S_r \frac{d}{dt} h_{i,r} - (h_i \vec{u}_i, \operatorname{grad} w_r) = 0 .
$$

Remarque : Les *dérivées spatiales* dans les expressions ci-dessus sont prises
bien sûr au sens ordinaire, triangle par triangle. Ainsi, par exemple :

$$
(h_i \vec{u}_i, \operatorname{grad} w_r) = \sum_{T/P} \operatorname{grad} w_r|_{T} \cdot \int_{T} h_i \vec{u}_i \, dxdy,
$$
les termes de convection se calculent aisément :

$$ ((u_i \vec{\nabla}) u_{x,i} w_p) = \sum_{T/P} \frac{\text{aire}(T)}{3} \vec{u}_i(p) \cdot \left[ \vec{\nabla}(u_{x,i}) \right]_T $$

les termes non linéaires de frottement sont calculés par la formule d'intégration d'ordre 2 utilisant les noeuds comme points d'intégration :

$$ \frac{|\vec{u}_i| \vec{u}_i}{h_i} w_p \approx \sum_{T/P} \frac{\text{aire}(T)}{3} \frac{\vec{u}_i(p)}{h_{i,p}} $$

### C) - Schémas de discrétisation en temps.

#### 1) - Schéma explicite :

L'indice m désignera les valeurs des inconnues au temps (m $\delta t$). En utilisant des dérivées centrées en temps (schéma "leap-frog"), on obtient à la place de la première équation (10) :

$(\forall p \in S_0)$ :

$$ (14) \quad S_p \frac{u_{x,i,p}^{m+1} - u_{x,i,p}^{m-1}}{2 \delta t} + \left( (\vec{u}_i^m \cdot \vec{\nabla}) u_{x,i}^m + g \frac{\partial}{\partial x} (\alpha_i h_1^m + h_2^m + D) - f u_{y,i}^m + C_F \frac{|u_i^m|}{h_i^m} \right) w_p + v (\text{grad } u_{x,i}^m, \text{grad } w_p) = \\ = \delta i_1 \left( \frac{w}{\rho_1^m} \frac{w_x}{u_1^m}, w_p \right). $$

Les autres équations se transforment de manière évidente.

### Remarques :

Dans le premier pas de temps (m = 0), on utilise un schéma prédicteur-correcteur de Matsuno :

$$ (15) \quad p \in S_0 : S_p \frac{u_{x,i,p}^1 - u_{x,i,p}^0}{\delta t} + ((\vec{u}_i^0 \cdot \vec{\nabla}) u_{x,i}^0 + \dots + w_p) + v (\text{grad } u_{x,i}^0, \text{grad } w_p) = \delta i_1 \frac{w}{\rho_1^0} \frac{w_x}{h_1^0}, w_p. $$
$$
(16) \quad S_p \frac{u_i^{\frac{1}{2}}}{\delta t} \stackrel{u_i^{0}}{=} \frac{u_{i,p}^{0}}{\delta t} + \left( (\vec{u}_i^0 \cdot \nabla) \vec{u}_{y,i} + \dots + w_p \right) + \dots
$$
$$
S_r \frac{h_i^{\frac{1}{2}}}{\delta t} \stackrel{h_i^{0}}{=} \frac{h_{i,r}^{0}}{\delta t} - \left( (\vec{h}_i^0 \cdot \nabla) \vec{w}_r \right) w_r = 0
$$

$$
(17) \quad S_p \frac{u_{x,i,p}^1 - u_{x,i,p}^0}{\delta t} + \left( \left( \vec{u}_i^{\frac{1}{2}} \cdot \vec{\nabla} \right) u_{x,i}^{\frac{1}{2}} + \dots + w_p \right) \\
\qquad + v \left( \vec{\nabla} u_{x,i}^{\frac{1}{2}} \cdot \vec{\nabla} w_p \right) \\
\qquad = \delta_{i1} \left( \frac{|w|_x}{h_1^{\frac{1}{2}}} , w_p \right)
$$

$$
(18) \quad S_p \frac{u_{y,i,p}^1 - u_{y,i,p}^0}{\delta t} + \left( \left( \vec{u}_i^{\frac{1}{2}} \cdot \vec{\nabla} \right) u_{y,i}^{\frac{1}{2}} + \dots + w_p \right) \\
\qquad + v \left( \vec{\nabla} u_{y,i}^{\frac{1}{2}} \cdot \vec{\nabla} w_p \right) = \delta_{i1} \left( \frac{|w|_y}{h_1 h_1^{\frac{1}{2}}} , w_p \right)
$$

$$
(19) S_r \frac{h_{i,r}^1 - h_{i,r}^0}{\delta t} - \left( h_i^{\frac{1}{2}} \vec{u}_i^{\frac{1}{2}} \cdot \vec{\nabla} w_r \right)
$$

Il est connu que le schéma "saute-mouton" induit un mode "parasite" en plus
de la solution "physique" : on peut éliminer le mode parasite en intercalant
un pas de Matsuno (équations (15) à (19)), par exemple tous les 5 pas de temps.

2) - *Stabilité du schéma semi-explicite*

Ainsi qu'on l'a vu plus haut, la vitesse de propagation des ondes de gravité
est donnée par

$$
\begin{align*}
c_1 &\approx \sqrt{g(h_1 + h_2)} \\
c_2 &\approx \sqrt{\frac{\epsilon g h_1 h_2}{h_1 + h_2}} \qquad \epsilon = \frac{\rho_2 - \rho_1}{\rho_2}
\end{align*}
$$

Par exemple :

$$
\begin{align*}
h_1 &= 100\text{m} \quad , h_2 = 1.000\text{m} = c_1 \approx 100\text{m/s} \\
c_2 &\approx 1\text{m/s} \text{ .}
\end{align*}
$$
En appliquant la condition classique C.F.L., on obtient la condition nécessaire

$$ (20) \qquad c_1 \delta t < \delta x, $$

Avec $\delta x \approx 3$ km et les données numériques précédentes, on obtient une estimation de l'ordre de grandeur du pas de temps :

$$ \delta t < 30 \text{ s.} $$

(en pratique, il faut prendre un pas de temps trois fois plus petit).

Ainsi, le schéma semi-sxpicite nécessite l'intégration d'au moins 3.000 pas de temps pour représenter l'évolution sur un jour. Comme le temps-machine nécessaire pour exécuter un pas de temps est de l'ordre de quelques secondes sur IBM 370-168, il apparaît nécessaire d'utiliser un autre schéma.

Un schéma totalement implicite serait très coûteux ; il faudrait en effet résoudre 6 équations non linéaires couplées.

### D) - Schéma semi-implicite

#### 1) - Traitement implicite des termes de gravité :

On a vu que la contrainte la plus importante sur le pas de temps provenait de la vitesse de propagation des ondes de gravité. On peut donc penser qu'un traitement implicite des seuls termes qui interviennent dans leur mécanisme permettra de prendre des pas de temps plus grands.

Le principe de la méthode numérique est le même que celui de Kwizak et Robert (1971). Pour ne pas alourdir l'exposé, le détail des calculs et des équations obtenues a été renvoyé en Appendice.

Par rapport au schéma temporel explicite précédent, les modifications sont les suivantes :

- dans les équations du mouvement pour chaque couche, le terme de gradient qui était pris au temps $m$, est remplacé par la demi-somme du même terme pris aux temps $m-1$ et $m+1$.

- dans les équations de continuité, le terme de divergence est d'abord séparé en parties linéaires et non-linéaires, et seule la partie linéaire initialement prise au temps $m$ dans le schéma explicite est remplacée par la demi-somme du même terme pris aux temps $m-1$ et $m+1$.
L'élimination des degrés de liberté sur les vitesses conduit alors à un système d'équations sur les hauteurs d'eau qui sont les analogues, discrétisées, des équations de propagation des ondes de gravité exposées plus haut (équations de Helmholtz discrétisées).

L'inversion du système des équations de propagation donne les nouvelles valeurs au temps (m+1) pour les hauteurs d'eau dans chaque couche, à partir desquelles il est immédiat de calculer de manière explicite les nouvelles valeurs au temps (m+1) des vitesses dans chaque couche.

## 2) - *Structure matricielle du système obtenu.*

On note

$H_i$ : le vecteur de degrés de liberté de $h_i^1$ (fluctuations de $h_i$),

$u_i$ : le vecteur de degrés de liberté de $\vec{u}_i$.

Pour chaque couche $i = 1,2$, les équations peuvent se mettre sous la forme matricielle :

$$ \eta u_i^{m+1} + g \delta t \beta(\alpha_i H_1^{m+1} + H_2^{m+1}) = v_i^m . $$

$$ \eta H_i^{m+1} - H_i \delta t C u_i^{m+1} = K_i^m, \quad H_i \text{ = constante.} $$

Les matrices $\beta$ et $C$ discrétisent respectivement les opérateurs grad. et (-div). $-v_i^m$ et $K_i^m$ désignent les seconds membres des deux équations.

L'élimination des vitesses donne :

$$ (22) \qquad \eta H_i^{m+1} + g H_i \delta t^2 C \eta^{-1} \beta(\alpha_i H_1^{m+1} + H_2^{m+1}) = K_i^m + H_i \delta t C \eta^{-1} v_i^m . $$

Les matrices de masse $\eta$ et $\eta$ sont diagonales avec les éléments finis linéaires non conformes : ceci n'est plus vrai avec d'autres éléments (1) et dans ce cas la matrice $C\eta^{-1}$. $\beta$ n'est plus une matrice creuse, ce qui pose des difficultés pour la résolution numérique.

<footnote>(1) A moins de faire du "mass-lumping" ; on a vu que cela présente des inconvénients. En outre, la structure des matrices $\eta$ et $\eta$ serait quand même moins simple.</footnote>
Pour résoudre les équations (22) qui sont une forme possible des équations de Helmholtz, on utilise la surrelaxation par blocs, en groupant les dérivations des deux couches en un même noeud. Le paramètre optimal de surrelaxation dépend fortement de $\delta t$, mais il peut être choisi une fois pour toutes au début du calcul, car la matrice du système linéaire (22) $a_{qr}^{(i)}$, $i=1,2$ pour chaque couche, ne dépend pas du temps.

### E) - Avantages de la classe d'éléments finis choisie.

On remarque que les inconnues étant définies aux milieux des côtés, les vecteurs normaux et tangents à la frontière sont bien définis et il est très facile de traiter les conditions aux limites de glissement (ou d'autres conditions sur $\vec{u}_i, \vec{n}$).

Il existe, avec l'approximation choisie, des champs de vitesse à divergence nulle. En effet, on sait que pour des écoulements bidimensionnels satisfaisant $\overline{\vec{u} \equiv 0}$, une approximation par des éléments conformes linéaires ne donne que la solution triviale $u \equiv 0$, car le nombre de degrés de liberté est plus faible.

Cette remarque est importante pour la simulation d'écoulements océanographiques, où l'on sait qu'un état d'équilibre stationnaire possible est l'équilibre géostrophique qui est à divergence horizontale nulle.

Dans le schéma explicite, la propriété d'orthogonalité des fonctions de base fait que le schéma reste totalement explicite (pas de système linéaire à résoudre). Dans un travail précédent, Bégis et Crépon [3] ont utilisé des éléments finis conformes (vitesses et hauteurs d'eau sont alors continues sur le domaine), mais dans ce cas, il faut utiliser une formule d'intégration approchée dans le calcul de $(\frac{\partial \vec{u}}{\partial t} \cdot \vec{v})$ pour conserver le caractère explicite du schéma :

$$ \int_{\Gamma} \frac{\partial u_i}{\partial t} \vec{v} \, dxdy \approx \frac{\text{aire (T)}}{3} \sum_{S \in T} \frac{\partial \vec{u}_i}{\partial t}(s) \cdot \vec{v}(s) $$

(la sommation s'étend aux 3 sommets du triangle).
Un procédé, dit de "mass-lumping" introduit une erreur, qui peut ne pas être négligeable pour un phénomène fortement transitoire. Toutefois, l'influence exacte de l'intégration numérique approchée ne semble pas avoir été étudiée à fond jusqu'à présent (voir Thacker [12]).

## F) - Coût de la méthode.

L'encombrement en mémoire (pour le schéma semi-implicite) est le suivant :
le stockage de la géométrie du domaine, de la matrice du système, et des variables nécessite environ :

65 × NT mots (¹) en simple précision (avec NT = nombre de triangles)
38 × NT doubles mots

soit au total : 141 × NT mots.

Lorsque la topographie D est variable, il faut encore ajouter environ (21 × NT) mots mais seulement 2 à 3 × NT mots si la zone de variation de D est peu étendue.

Afin de comparer les performances des schémas explicite et semi-implicite, on donne ci-dessous les temps de calcul (sur IBM 370/168 du CIRCE) nécessaires pour faire évoluer de 0 à 5 heures le modèle académique décrit en II.

<table><thead><tr><th rowspan="2">t</th><th>explicite</th><th colspan="4">semi-implicite</th></tr><tr><th>1'</th><th>5'</th><th>10'</th><th>30'</th><th>60'</th></tr></thead><tbody><tr><td>.</td><td>--</td><td>1.2</td><td>1.4</td><td>1.63</td><td>1.68</td></tr><tr><td>Temps machine<br>pour intégrer<br>de 0 à 5 heures</td><td>7'</td><td>2'50"</td><td>1'26"</td><td>46"</td><td>33"</td></tr></tbody></table>

<footnote>(¹) Ce chiffre peut être réduit de moitié en utilisant pour certains entiers des mots de 2 octets, facilité offerte sur IBM.</footnote>
### III. - MODELE ACADEMIQUE

Nous avons voulu dans un premier temps tester la validité des algorithmes mis au point sur un cas académique, pour lequel il existe une solution *analytique* au problème linéaire transitoire.

Le cas "linéaire" correspond au cas où l'on néglige les termes d'advection et où les épaisseurs des deux couches ne sont supposées constantes que dans les termes suivants :

$$ \frac{\tau}{\rho h_1} \quad \text{et} \quad \text{div}(h_i u_i). $$

#### A). - Contexte physique

Le contexte "physique" de cette étude académique est l'étude de la réponse du système à un vent favorable à l'upwelling, c'est-à-dire à la côte dans le cas d'une côte rectiligne, et tel que le champ de vent présente des variations dans la direction *parallèle* à la côte et soit transitoire.

Le cas idéalisé le plus simple est un vent impulsionnel ne soufflant que dans une bande limitée en largeur dans la direction de la côte (cf. Allen [2], Crépon [4]). On a simulé le cas où la tension dûe au vent est une fonction "chapeau" dans cette bande (Fig. 4).

L'essentiel de la réponse du système est l'excitation d'une onde de Kelvin interne, se propageant à la célérité des ondes de gravité internes vers le sud (en laissant la côte à sa droite), c'est-à-dire dans la direction opposée au vent et présentant une décroissance exponentielle vers le large sur une échelle spatiale de l'ordre du rayon interne de déformation (~10km.).

Les caractéristiques de l'onde sont résumées sur la Figure 4. Cette onde, qui cause donc un relèvement important de l'interface aux points de la côte où elle masse, peut donc engendrer un upwelling en dehors de la zone du "forcing" atmosphérique, autrement dit, en un point donné de la côte, l'upwelling peut ne pas être corrélat au champ de vent local.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 4 : (a) montre la variation de la tension exercée ar le vent selon l'axe y : le "forcing" est limité à la bande [-b, +b]</figcaption></figure>

<caption>(b) ... (e) montrent l'allure de la dénivellation de l'interface à la côte à différents instants ; noter la dissymétrie de la réponse dûe à la présence d'une onde de Kelvin se propageant vers y < 0, dans la direction opposée au vent.
Cf. Crépon (1979).</caption>
On peut remarquer également que l'amplitude de l'onde atteint une valeur asymptotique lorsque le temps croît ; donc, dans la zone "forte", bien qu'il y ait un apport continu d'énergie par le vent, la dénivellation de l'interface reste bornée dans le temps, et l'excès d'énergie apporté par le forcing du vent est propagé en dehors de la zone forcée par l'onde de Kelvin interne.

<figure><img src="image_13.png" /><figcaption>Figure 5</figcaption></figure>

Transport dans la couche de surface

A la côte, le transport parallèle lié à l'existence de l'onde de Kelvin, compense directement le transport d'Eckman dans la bande forcée par le vent. Crépon (1979)
La Figure 5 montre qu'à l'onde de Kelvin interne est associé en surface un transport longitudinal, parallèle à la côte, décroissant exponentiellement dans la direction du large et que ce transport compense directement le transport induit par le vent vers le large dans la bande où le vent souffle (transport d'Ekman). Ceci exclut en particulier, une condition aux limites d'adhérence ($u.n \equiv 0$ et $\vec{u} \cdot \vec{t} \equiv 0$) qui empêcherait l'existence d'un tel cou- rant longitudinal collé à la côte. Ce sont les principaux traits de la physique du problème que l'on essaie de retrouver dans les résultats de la simulation numérique.

## B) - Résultats numériques.

Le maillage du domaine académique est donné sur la Figure 6. L'extension N.-S. du domaine (direction parallèle à la côte) est de 125km, l'extension N.-W. est de 100 km.

Le vent souffle dans une bande de largeur N.-S. de 50 km. Dans son extension E.-W. cette bande de vent est limitée à une largeur de 37km., ce qui est nettement supérieur au rayon interne de déformation et ne perturbera pas trop la solution à la côte.

La distance minimale de résolution dans la direction perpendiculaire à la côte est de 3km, ce qui permet donc bien d'atteindre le rayon interne de déformation (10 km.).

On voit nettement sur les Figures 7A, 7B, 7C, qui représentent les isolignes de la dénivellation de l'interface, à différents instants, l'existence d'une onde de Kelvin se propageant vers le sud. L'intégration du modèle est arrêtée lorsque l'onde atteint le bord sud du domaine pour éviter des problèmes de fausses réflexions d'onde.

On remarque la présence d'un downwelling à la limite Est de la zone de vent, ce qui est consistant avec la présence d'un rotationnel du champ de vent induit par la limitation du champ de vent dans la direction pendiculaire à la côte.
<figure><img src="image_14.png" /><figcaption>Fig. 6</figcaption></figure>
<figure><img src="image_15.png" /></figure>

T =16 H. 0 MN. 0 S. LIGNES ISO-H (THERMOCLINE)

AAA7905820150146 BLOC NO 65

Figure 7 A. : Isolignes de la dénivellation de l'interface, 16 heures après le début du vent.

On commence à voir une onde de Kelvin (décroissance exponentielle de la dénivellation vers le large), qui commence à se propager vers le sud, dans la direction opposée au vent.

A la limite Est de la zone de vent, un léger downwelling est dû à la décroissance du champ du vent.
<figure><img src="image_16.png" /><figcaption>T =30 H. 0 MN. 0 S. LIGNES ISO-H (THERMOCLINE)</figcaption></figure>

ΔΔΔ7905920265896 BLOC NO 19

Figure 7 B. : Isolignes à 30 heures.
<figure><img src="image_17.png" /><figcaption>Figure 7 C. : Isolignes à 40 h.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

T = 30 H. 0 MN. 0 S. VITESSES DANS LA COUCHET

VITESSE MAX.=0.107E00

AAA7905920265896 BLOC NO 16

<caption>Fig. 8. A : Vitesses dans la 1e couche à 30 h.</caption>

Les vitesses à la côte longitudinales sont induites par le passage
de l'onde de Kelvin.
<figure><img src="image_19.png" /></figure>

T. =30 H. O MN. O S. VITESSES DANS LA COUCHE 2
VITESSE MAX.=0.019E00
AAA7905920265896 BLOC NO 17

<caption>Fig. 8. B : Vitesses dans la 2e couche à 30 h.</caption>

Les vitesses parallèles à la côte coincident avec le passage
de l'onde de Kelvin.
<figure><img src="image_20.png" /><figcaption>Figure 9</figcaption></figure>
Le transport longitudinal (parallèle à la côte) dans le première couche, est nettement visible sur les champs de vitesse des figures 8, 9.
Le transport longitudinal important coincide avec les points atteints par l'onde de Kelvin.

Au début de l'intégration, les effets non-linéaires sont relativement faibles et la Figure 9 montre le bon accord entre la solution linéaire analytique et la solution numérique ($t \approx 16h. = \frac{2b}{c}$) sur l'allure spatiale de l'interface à la côte.

On peut en conclure que les effets non linéaires pris en compte dans le modèle numérique restent faibles pendant le temps d'intégration concerné. Leurs principaux effets sont :

* propagation plus lente de l'onde de Kelvin qui reste au début plus longtemps dans le zone de forcing, ce qui explique l'allure plus confinée spatialement et une dispersion *non linéaire* de l'onde :
  le front à l'avant de l'onde se propage plus lentement d'où il résulte l'aspect pentu et dissymétrique de la "crête" de l'onde.

* État asymptotique : la dénivellation calculée de l'interface est importante par rapport au résultat donné par la solution linéaire :
  le terme du forcing est en $\frac{\tau}{\rho h_1}$, morsqu'il y a upwelling, c'est-à-dire lorsque $h_1$ diminue, le terme du forcing croît et donc l'apport d'énergie croît. Ces résultats qualitatifs sur les effets non-linéaires sont en accord avec les conclusions de Adamec et O'Brien [1].

L'évolution au cours du temps des dénivellations et des champs de vitesse montrent que les oscillations d'inertie de Coriolis décroissent lorsque le temps croît, ce qui est aussi en accord avec la solution linéaire transitoire analytique, (cf. Crépon [4]).
<figure><img src="image_21.png" /><figcaption>Fig. 10</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 11 :</figcaption></figure>

Champ de la tension exercée
par le vent : W ≈ 10m/s. Le vent
est nul au-delà de 200 km du cap
d'Agde. A ≡ Marseille ;
B ≡ Pointe de l'Espiguette ;
C ≡ Cap d'Agde.
T =32 H. O TN. O S. LIGNES ISO-H (THERMOCLINE)
ADDY909522330800 BLOC NO 5

Fig. 12. A : Isolignes du relèvement de l'interface 32 heures après le début du vent (les hauteurs sont exprimées en mètres).

A ≡ Marseille

B ≡ Pointe de l'Espiguette

C ≡ Cap d'Agde

Les points désignent l'emplacement de maximas relatifs à la dénivellation de l'interface.

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 12. B : Isolines du relèvement de l'interface, 45 heures après le début du vent.
Même légende que Fig. 11. A.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 12. C : Isolignes de la dénivellation de l'interface, 60 heures après le début du vent.
La dénivellation est de 28 m à la pointe de l'Espiguette et l'interface est près
de crever la surface libre.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. 13. A : Champs de vitesse de la 1e couche.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. 13. B : Champs de vitesses dans la 2e couche.
La vitesse maximale est de 6,7 cm/s.</figcaption></figure>
## IV. - CAS DU GOLFE DU LION.

### A) . Maillage du domaine.

La zone correspondant au Golfe du Lion est raffinée en respectant la direction privilégiée de la rupture de pente du plateau continental.

Les frontières océaniques au large du Golfe ont été repoussées le plus loin possible jusqu'à des limites "naturelles" : côtes de Corse, Sardaigne, limite Est de la Mer Ligure, au sud, limite de la latitude 40° N., Baléares au S.-W., etc., ceci pour éviter qu'il y ait des problèmes de réflexions parasites d'ondes à ces frontières. Il faut environ 6 jours pour qu'une onde interne engendrée dans le Golfe du Lion à la côte, puisse atteindre la frontière "océanique" la plus proche : en fait, le temps d'intégration du modèle est au maximum de 3 jours, à cause du problème de l'interface venant crever la surface libre.

### B) . Champ de vent.

L'étude des données météorologiques montrent que les effets orographiques sont assez importants sur le Golfe, en ce qui concerne les variations spatiales et temporelles du vent.

Bien que le modèle puisse prendre en compte un vent variable dans le temps et l'espace, celui-ci reste encore mal mesuré et constitue une forte limitation pour des simulations "réalistes".

L'extension spatiale du champ de vent (cas d'une tramontane généralisée sur le Golfe du Lion) est donnée sur la Figure 11. L'intensité maximale est de 10m/s.

### C). Résultats numériques.

*Champs de dénivellation de l'interface.*

On sait que la variation de la température des eaux de surface est une des signatures les plus immédiates d'un phénomène d'upwelling : pour un temps inférieur à 3 jours, les phénomènes d'advection qui pourraient transporter le champ de la température peuvent être négligés et les endroits où le relèvement de l'interface sera le plus important correspondant aux eaux de surface les plus froides, car le brassage vertical y est le plus fort.
Ceci semble confirmé par la comparaison des Figures 1 et 12 (A,B,C) qui montrent respectivement les isolignes de la température de surface obtenues à partir de satellite, moyenne sur plusieurs coups de vent et les isolignes du relèvement de l'interface à différents instants après le début de la tramontane. On peut remarquer la coincidence de la localisation des maxima de relèvement de l'interface avec les points sources les plus froids et également au large de la pointe de l'Espiguette d'allure en faisceau des isolignes sur les deux figures.

Ceci semble confirmer que la variabilité à l'échelle de 10km de la direction et des décrochements de la côte soit essentielle pour expliquer l'aspect en taches ("patchines") des upwellings observés l'été dans le Golfe du Lion. Ces taches sont, d'autre part, bien confinées dans une zone de largeur de l'ordre de 10 km de la côte.

Les observations des champs de vitesse (fig. 13) dans les deux couches montrent :

- le transport d'Ekman dans la couche au large, à la droite du vent en opposition de phase avec le courant dans le deuxième couche,

- une structure compliquée le long de la côte dans les deux couches et correspondant aux phénomènes de re-circulation entre les diverses cellules d'upwelling imposés par la conservation de la masse. Des courants importants sont situés dans une zone intéressant une largeur de l'ordre de 10km à partir de la côte.

- dans la zone intermédiaire entre le large et cette couche limite côtière, la variabilité observée dans les courants obtenus par la simulation numérique n'est pas aussi grande que celle qui est observée en réalité dans le Golfe du Lion, et qui serait plutôt induite par des variations locales dans le cham de vent, paramètre actuellement non accessible.

## CONCLUSION

Au vu des résultats ci-dessus, il semble que l'emploi d'une méthode d'éléments finis soit essentielle pour :

- reproduire certaines "variabilités" des écoulements observés dans l'océan : "patchines" des écoulements observés à partir de données satellites, etc.)

- permettre de prendre en compte des géométries compliquées de la ligne de côte.

Les étapes ultérieures de ce travail seront :

- le passage à une maille plus fine ($\delta x \sim 3$km) sur le Golfe du Lion ;

- la prise en compte de processus physiques plus complexes dans le modèle (par exemple, la possibilité de mélanges entre les deux couches).
# APPENDICE

## Schéma semi-implicite

Avant de présenter le schéma, décomposons les hauteurs $h_i$ selon :

$$h_1(x,y,t) = h'_1(x,y,t) + H_1$$

$$h_2(x,y,t) = h'_2(x,y,t) + H_2 - \mathcal{D}(x,y)$$

<figure><img src="image_28.png" /></figure>

$H_1$ et $H_2$ sont des hauteurs moyennes, constantes dans le temps et dans l'espace. On pose :

$$h'_{i,r}(t) = h'_{i}(r,t) = \text{valeur de } h'_{i} \text{ au noeud } r$$

$$D_{i,r} = D_{i}(r)$$

$$h'_{i}(x,y;t) = \sum_{r \in g} h'_{i,r}(t) \omega_{r}(x,y)$$

$$h'_{i}^{m}(x,y) = h'_{i}(x,y, m\delta t)$$

$$h'_{i,r}^{m} = h'_{i,r}(m\delta t)$$

$$\alpha_1 = 1 \qquad \alpha_2 = 1 - \epsilon$$

On utilise $H_1$ et $H_2$ pour linéariser les équations ; on prend un schéma "saute-mouton", où les termes de gravité sont traités implicitement, et qui s'écrit avec les notations précédentes :

$(\forall p \in g_0)$

$$S_p \frac{u_{x,i,p}^{m+1} - u_{x,i,p}^{m-1}}{2 \delta t} + g\left(\frac{\partial}{\partial x}\left(\alpha_i \frac{h_1'^{m+1} + h_1'^{m-1}}{2} + \frac{h_2'^{m+1} + h_2'^{m-1}}{2}\right), w_p\right) \\ + \left(\left(u_{i,\vec{v}}^{m-1}\right)u_{i,x}^{m} + C_{F,i} \frac{|u_{i,x}^{m}|^m}{h_i^m} - \delta_{i1} \frac{|w|w_x}{\rho_1 h_1^m} - f \left(u_{y,i,p}^{m-1}, w_p\right)\right)$$
$$ + v(\vec{\text{grad}} u_{x,i}^m, \vec{\text{grad}} w_p) = 0 $$

$$
\begin{align*}
& S_p \times \frac{u_{V,i,p}^{m+1} - u_{V,i,p}^{m-1}}{2\delta t} + g\left(\frac{\partial}{\partial y}\left(\alpha_i \frac{h_1^{m+1} + h_2^{m-1}}{2} + \frac{h_2^{m+1} + h_1^{m-1}}{2}\right), w_p\right) \\
& \quad + ((u_i^m \cdot \nabla) u_{i,y}^m) C_{F_i} \frac{|u_i^m|_{V,i,y}^m}{h_i^m} - \delta_{i1} \frac{|w_i|_{V,i}^m}{\rho_1 h_1^m} + f_i u_{x,i,p}^m w_p \\
& \quad + v(\vec{\text{grad}} u_{y,i}^m, \vec{\text{grad}} w_p) = 0
\end{align*}
$$

$$ (\forall p \in g_1) : $$

$$
\begin{align*}
& S_p \frac{u_{i,p}^{m+1} - u_{i,p}^{m-1}}{2 \delta t} + g(\vec{\text{grad}}(\alpha_i \frac{h_1^{m+1} + h_2^{m-1}}{2} + \frac{h_2^{m+1} + h_1^{m-1}}{2}), \vec{\tau}_p w_p) \\
& \quad + ((u_i^m \cdot \nabla) u_i^m + C_{F_i} \frac{|u_i^m|^2}{h_i^m} - \delta_{i1} \frac{|w_i|^2}{\rho_1 h_1^m}, \vec{\tau}_p w_p) \\
& \quad + v(\vec{\text{grad}} u_{x,i}^m + \vec{\text{grad}} u_{y,i}^m, \vec{\text{grad}} w_p) \\
& = 0
\end{align*}
$$

$$ \forall r \in g : $$

$$
\begin{align*}
S_r \frac{h_1^{m+1} - h_2^{m-1}}{2\delta t} - H_1 \left( \frac{u_1^{m+1} + u_2^{m-1}}{2} \right), & \quad \text{grad } w_r \\
& = (h_1^m \vec{u}_1^m, \text{grad } w_r) = 0
\end{align*}
$$

$$
\begin{align*}
S_r \frac{h_2^{m+1} - h_1^{m-1}}{2 \delta t} - ((H_2 - D) \frac{u_2^{m+1} + u_1^{m-1}}{2}), & \quad \text{grad } w_r \\
& = (h_2^m \vec{u}_2^m, \text{grad } w_r) = 0
\end{align*}
$$

Il s'agit maintenant de résoudre ces six systèmes d'équations couplées. En fait grâce à l'orthogonalité des fonctions de base, on voit qu'il est possible d'éliminer les vitesses (de la même façon qu'au paragraphe I-D pour former les
équations de propagation des ondes).

On peut réécrire les équations sous la forme : (1)

$$
\forall p \in g : S_p u_{x,i,p}^{m+1} + g \delta t \left( \frac{\partial}{\partial x} (\alpha_i h_1^{m+1} + h_2^{m+1}) \cdot w_p \right) = S_p v_{x,i,p}^m
$$

$$
S_p u_{y,i,p}^{m+1} + g \delta t \left( \frac{\partial}{\partial y} (\alpha_i h_1^{m+1} + h_2^{m+1}) \cdot w_p \right) = S_p v_{y,i,p}^m
$$

$$
\forall r \in g : S_r h_{1,r}^{m+1} - H_1 \delta t (\vec{u}_1^{m+1}, \text{grad } w_r) = S_r k_1^m
$$

$$
S_r h_{2,r}^{m+1} - \delta t ((H_2 - D)u_{2,q}^{m+1}) \vec{\text{grad}} w_r = S_r k_2^m
$$

L'élimination des degrés de liberté des vitesses dans les équations précédentes
donne :

Vr ∈ g;

$$
S_r h_{1,r}^{m+1} + g H_1 \delta t^2 \sum_{q \in \phi} a_{qr}^{(1)} \times (h_{1,q}^{m+1} + h_{2,q}^{m+1}) \\
= S_r k_{1,r}^m + E_1 \delta t \sum_{p \in \phi_0} v_{x,1,p} (w_p, \frac{\partial w}{\partial x}) + v_{y,1,p} (w_p, \frac{\partial w}{\partial y})
$$

$$
S_r h_{2,r}^{m+1} + g \delta t^2 \sum_{q \in \phi} a_{qr}^{(2)} \times ((1-\epsilon)h_{2,q}^{m+1} + h_{2,q}^{m+1}) \\
= S_r k_{1,r}^m + \delta t \sum_{p \in \phi_0} v_{x,2,p}((H_2-D)w_p, \frac{\partial w}{\partial x}) + v_{y,2,p}((H_2-D)w_p, \frac{\partial w}{\partial x})
$$

avec :

$$
a_{qr}^{(1)} = \sum_{p \in g_0} S_p^{-1} \left( \left( w_p, \frac{\partial w}{\partial x} \right) \left( w_p, \frac{\partial w}{\partial x} \right) + \left( w_p, \frac{\partial w}{\partial y} \right) \left( w_p, \frac{\partial w}{\partial y} \right) \right)
$$

$$
a_{qr}^{(2)} = \sum_{p \in g_0} S_p^{-1} \left( \left( (H_2 - D) w_p, \frac{\partial w}{\partial x} \right) \left( (H_2 - D) w_p, \frac{\partial w}{\partial x} \right) + \left( (H_2 - D) w_p, \frac{\partial w}{\partial y} \right) \left( (H_2 - D) w_p, \frac{\partial w}{\partial y} \right) \right)
$$

<footnote>(1) En supposant que $u=0$ sur $\Gamma$ ; on laisse au lecteur le soin d'établir les formules dans le cas général</footnote>
La structure des matrices $a_{qr}^{(1)}$ et $a_{qr}^{(2)}$ à inverser est assez simple :

<figure><img src="image_29.png" /></figure>

$a_{qr}^{(i)}$ est non nul si :

- ou bien $q$ est voisin de $r$

- ou tien $q$ est voisin de $p$, et $p$ voisin de $r$.

Ces deux matrices sont donc très *creuses* : le stockage nécessaire pour chacune d'elle (en vue de résoudre par surrelaxation) est de l'ordre de :

(13' x NM) avec NM = nombre de côtés. ≈ 3/2 nombre de triangles.

Si la profondeur D est constante sur une partie du domaine, $a_{qr}^{(2)} = a_{qr}^{(1)} \times (H_2-D_r)$ pour les points de cette partie ; si de telles régions à fond plat sont importantes il est avantageux de ne stocker que les éléments $a_{qr}^{(2)}$ qui différent de $a_{qr}^{(1)} \times (H_2-D_r)$ (ceux qui correspondent à des régions à profondeur variable).

Ceci peut se faire au prix d'un tableau de longueur NM.
BIBLIOGRAPHIE

[1] Adamec et J.J. O'Brien "The seasonal upwelling in the Gulf of Guinea to remote forcing" J. of Physical Oceanography, Vol. 8, n° 6, Novembre 1978, pp. 1050-1060.

[2] J.S. Allen "Some aspects of the forced-wave response of stratified coastal regions", J. of Physical Oceanography, vol. 6 n° 1 (1976), pp. 113-119.

[3] D. Bégis, M. Crépon, C. Millot "Un problème d'environnement : modélisation de la dynamique océanique du plateau continental" Rapport de Recherche IRIA-LABORIA, n° 236, Mai 1977.

[4] M. Crépon "Transient upwelling generated by 2-dimensional atmospheric forcing" (1979), to appear.

[5] M. Crouzeix, P.A. Raviart "Conforming and non-conforming finite element methods for solving the stationary Stokes equations", R.A.I.R.O., (R-3), Décembre 1973, pp. 33-75.

[6] H.E. Hurlburt, J.J. O'Brien "A numerical method of coastal upwelling", J. of Physical Oceanography, vol. 2, N° 1, Janvier 1972, pp. 11-26.

[7] M. Kwizak, A.J. Robert, "A semi-implicit scheme for grid-point atmospheric models of the primitive equations" Mon. Wes. Rev. 99, (1971), pp. 32-36.

[8] C. Millot "Wind induced upwellings in the Gulf of Lions. Oceanologica Acta, (1979), vol. 2 n° 3.

[9] Peffley, J.J.O'Brien, "A three-dimensional simulation of coastal upwelling of Oregon", J. of Physical Oceanography, vol 6, n° 2, (1976), pp. 164-180.

[10] F. Thomasset "Numerical solution of the Navier-Stokes equations by finite element methods", Von Karman Institute, Rhode St. Genèse, Belgium, Lecture Series 86, (1977).

[11] R. Temam et F. Thomasset "Numerical solution of the Navier-Stokes equations by a finite element method", Proceedings of the "IInd International Symposium, Finite Element Methods in Flow Problems", Santa Margherita Ligure, Italy (1976), Springer ed.

[12] W.C. Thacker, "Comparison of Finite Element and Finite Differences Schemes", Journal of Physical Oceanography, vol. 8 (1978).

[13] M. Benahmed "Identification de non-linéarités ou de paramètres répartis dans deux équations aux dérivées partielles non linéaires modélisant un gisement d'hydrocarbures et un écoulement dans un canal ouvert". Thèse d'Etat, Université de Paris IX, (1978).