J Demange

L Debreu

P Marchesiello

F Lemarié

E Blayo

On the use of a depth-dependent barotropic mode in ocean models: impact on the stability of the coupled barotropic/baroclinic system

Keywords: Mode splitting, Barotropic mode, Stability analysis

Evolution of the oceanic free-surface is responsible for the propagation of fast surface gravity waves which roughly propagates at speed gH (with g the gravity and H the local water depth). In the deep ocean, this phase speed is roughly two orders of magnitude faster than the fastest internal gravity waves. The steep stability constraint imposed by those fast surface waves on the time-step of numerical models is handled using a splitting between slow (internal / baroclinic) and fast (external / barotropic) motions to allow the possibility to adopt specific numerical treatments in each component. The barotropic mode is traditionally approximated by the vertically integrated flow because it has only slight vertical variations. However, the implications of this assumption on the stability of the splitting are not well documented. In this paper, we describe a stability analysis of the mode-splitting technique based on an eigenvector decomposition using the true (depth-dependent) barotropic mode. We show that the use of such a depth-dependent barotropic mode allows a much stable integration of the mode-split equations. As a consequence, the amount of dissipation required to achieve stable integrations, usually applied through averaging filters, can be drastically reduced. It results in a much improved effective resolution even for complex flows. In addition, the formulation of a new mode splitting algorithm using the depthdependent barotropic mode is introduced. The benefits of this new formulation are illustrated by idealized numerical experiments.

Sur l'utilisation d'un mode barotrope dépendant de la profondeur dans les modèles d'océan: impact sur la stabilité du système couplé barotrope/barocline Résumé : Des déplacements de la surface libre de l'océan résulte la propagation d'ondes de gravité rapides qui se propagent à une vitesse de l'ordre de gH (avec g la gravité et H la hauteur d'eau locale). Là où l'océan est profond, cette vitesse de phase est environ 2 ordres de grandeur plus rapide que l'onde interne de gravité la plus rapide. La forte contrainte de stabilité imposée par ces ondes rapides de surface sur le pas de temps des modèles est traitée en exploitant une séparation d'échelles entre les mouvements lents (i.e. le mode barocline/interne) et rapides (i.e. le mode barotrope/externe) afin de permettre l'utilisation de schémas numériques adaptés à chaque composante. Le mode barotrope est traditionnellement approché par l'intégrale sur la profondeur de l'écoulement car ses variations sur la verticale sont généralement faibles. Cependant, les conséquences de cette hypothèse sur la stabilité de la séparation d'échelles ont été très peu étudiées jusqu'à maintenant. Dans ce papier, nous développons une analyse de stabilité de la technique de séparation d'échelles en nous basant sur une décomposition en vecteurs propres exhibant le "vrai" mode barotrope (i.e. qui varie avec la profondeur). Nous montrons que l'utilisation de ce mode barotrope dépendant de la profondeur permet une intégration beaucoup plus stable du système couplé barotrope/barocline. Par conséquent, la quantité de dissipation habituellement requise pour stabiliser cette intégration (via des filtres de moyennage) peut être drastiquement réduite. De plus, la formulation d'un nouvel algorithme de séparation d'échelles faisant intervenir un mode barotrope dépendant de la profondeur est introduit. Les avantages de cette nouvelle formulation sont illustrés à l'aide d'expériences numériques.

Mots-clés : Séparation d'échelles, mode barotrope, analyse de stabilité

Introduction

Most of the current ocean circulation models have now relaxed the rigid lid assumption and directly integrate a prognostic equation for the free surface evolution. In this case, external gravity waves which roughly propagate at speed gH (with g the gravity and H the water depth) are explicitly resolved. Because a tridimensional implicit time stepping algorithm seems computationally impractical, the numerical integration of the gravest external waves introduces a strong stability constraint on the model time step. Under several assumptions (including flat bottom), the linear stratified primitive equations can be projected onto a set of normal modes which leads to an orthogonal separation of the modes. This highlights a fast (barotropic) mode and slow (baroclinic) modes.

The usual approach in state-of-the-art oceanic models is to introduce a mode splitting where the barotropic mode is integrated separately from the baroclinic modes. In order to satisfy the stability condition, the barotropic part can be integrated using either a 2D implicit time stepping algorithm [START_REF] Dukowicz | Implicit free-surface method for the bryan-coxsemtner ocean model[END_REF] or a time splitting approach [START_REF] Peter | The development of a free-surface bryan-cox-semtner ocean model[END_REF]. The use of an implicit algorithm has several drawbacks in particular in term of accuracy (large dispersive errors) and performance on parallel computers, especially at high resolution (poor scaling properties). As a consequence, most of the current ocean models use a time splitting approach where the barotropic component is integrated with several small time steps within a larger baroclinic time step ( [START_REF] Berntsen | Efficient numerical simulation of ocean hydrodynamics by a splitting procedure[END_REF], [START_REF] Peter | The development of a free-surface bryan-cox-semtner ocean model[END_REF]). In this context, the usual assumption is to consider that the external mode is vertically constant (i.e. depth independent). However, as pointed out in [START_REF] Peter | The development of a free-surface bryan-cox-semtner ocean model[END_REF], the density field has a nonzero barotropic component which invalidates the depth independent assumption.

Although depth-averaging provides an approximation to the barotropic mode , this approximation leads to a non orthogonal separation of slow and fast modes, even in the linear case. A reconciliation of the estimates of the barotropic mode coming from the 2D and the 3D parts is thus needed, and stability of the split between fast and slow modes is not guaranteed (indeed, vertical averaging does not exactly separate out the fast and slow dynamics).

To be more specific, let us introduce the 2D barotropic system linearized around a state of rest

       ∂u ∂t + g ∂η ∂x = - 1 ρ 0 ∂ ∂x 1 H 0 -H p h dz ∂η ∂t + ∂Hu ∂x = 0 ( 1
)
where ū is the depth average flow, η the free surface, and the right hand side is the depth integrated internal pressure gradient. This last term is held constant during barotropic integration because the associated purely barotropic component is expected to be rather small. However, we argue that this term is the main source of instability of the splitting under the depth-independent assumption. To prevent instabilities associated to these splitting errors, a time filtering of the barotropic variables has to be applied. This can be achieved either by using a diffusive time stepping algorithm in the barotropic integration itself [START_REF] Peter | The development of a free-surface bryan-cox-semtner ocean model[END_REF] or in the baroclinic integration ([8]). This time-filtering is usually done using explicit averaging filters ( [START_REF] Nadiga | On simulating flows with multiple time scales using a method of averages[END_REF], [START_REF] Shchepetkin | Computational kernel algorithms for fine-scale, multiprocess, longtime oceanic simulations[END_REF]) applied to the barotropic solution at all subtime steps. In addition to the inexact splitting above, several other reasons motivate the need for some form of time filtering [START_REF] Shchepetkin | The regional oceanic modeling system (roms): a split-explicit, free-surface, topography-following-coordinate oceanic model[END_REF][START_REF] Shchepetkin | Computational kernel algorithms for fine-scale, multiprocess, longtime oceanic simulations[END_REF]. First, in the nonlinear case, the right hand side of (1) is integrated from -H to η. Because the free surface η itself evolves during the barotropic subtime steps, the consistency is not maintained. An efficient remedy based on a redefinition of the barotropic pressure-gradient terms to account for the local variations in density field is proposed in [START_REF] Shchepetkin | The regional oceanic modeling system (roms): a split-explicit, free-surface, topography-following-coordinate oceanic model[END_REF]. Second, aliasing errors due to nonlinearities (when an integral of the advection terms u • ∇u is also included in (1)) are an additional source of instability controlled by time averaging . The benefit effect of recomputing at least the fast part of this term has been studied in [START_REF] Morel | Time splitting and linear stability of the slow part of the barotropic component[END_REF].

The paper deals with the impact of the aforementioned inexact splitting on the stability of numerical models. In [START_REF] Hallberg | Stable split time stepping schemes for large-scale ocean modeling[END_REF][START_REF] Higdon | Barotropic-baroclinic time splitting for ocean circulation modeling[END_REF][START_REF] Kamenkovich | On the time-splitting scheme used in the princeton ocean model[END_REF] , detailed numerical analysis of various 2D/3D time stepping algorithms are presented. Those studies specifically looked at the impact of adding diffusion either in the 2D or in the 3D part of a numerical model. In the present work, the focus is also on the quantification of the impact of the depth-independent assumption on stability of the modesplitting procedure. Indeed, among aliasing errors, enforcement of 2D/3D consistency, and the inexact splitting, it is yet not well understood which source of error is the most damaging for stability (and consequently which one of those errors justifies the amount of extra diffusion put through time-filtering). An other objective is to provide a general framework for the stability analysis of the mode splitting approach. Unlike previous studies, our analysis is based on a projection of the solution on the basis formed by unapproximated (depth dependent) barotropic and baroclinic modes.

The analysis done here leads to a simple expression of the minimum amount of dissipation required to achieve a stable integration of the mode-split equations. We show that usual explicit averaging filtering techniques are responsible for a much larger dissipation than needed for stability reasons.

The paper is organized as follows. In a first section, we briefly recall the classical theory of vertical mode decomposition and introduce important parameters for the stability analysis. We then proceed to the stability analysis focusing on the amount of diffusion required to control the instability arising from the depth independent barotropic mode assumption. Idealized numerical experiments are performed and the implementation of an alternative splitting technique based on the depth-dependent barotropic mode is introduced.

Normal mode decomposition

First, we briefly recall the normal mode decomposition theory (e.g [START_REF] Kundu | Fluid Mechanics, Second Edition[END_REF][START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF]). We consider a linearization of the 2D (x-z) primitive (i.e. under the hydrostatic and Boussinesq assumptions) equations around a motionless state in hydrostatic equilibrium

∂u ∂t + 1 ρ 0 ∂p ∂x = 0 (2) ∂p ∂z = -gρ (3)
∂u ∂x + ∂w ∂z = 0 (4) ∂ρ ∂t + w dρ dz = 0 (5)
In the vertical direction the model extends from the flat bottom z = -H to the top given by the free surface elevation z = η(x, t). Here u(x, z, t) and w(x, z, t) denote the perturbation components of horizontal and vertical fluid velocities (primes are omitted for clarity). p(x, z, t) and ρ(x, z, t) denotes pressure and density perturbations around a state (p(z), ρ(z)) satisfying the hydrostatic relation dp(z) dz = -ρ(z)g where ρ(z) is a reference density profile.

Inria

The (linearized) surface and bottom boundary conditions read

∂η ∂t = w(z = 0) at z = 0 (6) w = 0 at z = -H (7)
On a vertical grid with n layers, we assume that the discrete solution of (2-5) can be decomposed using vertical modes M q (z) ( [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF])

u(x, z, t) = n-1 q=0 u q (x, t)M q (z) (8) p(x, z, t) = ρ 0 g n-1 q=0 h q (x, t)M q (z) (9) ρ(x, z, t) = -ρ 0 n-1 q=0 h q (x, t) dM q (z) dz ( 10
)
where the modes M q (z) are the eigenvectors of the following Sturm Liouville problem

ΛM q = λ q M q (11) dM q dz z=-H = 0 (12) dM q dz z=0 = - N 2 (0) g M q (0) (13)
where Λ = -d dz N -2 d dz , with N 2 (z) = -g ρ 0 dρ dz the Brünt-Vaisala frequency assumed to be positive. Because Λ is a compact symmetric operator, it admits a basis of orthonormal eigenvectors M q (z) with positive eigenvalues λ q . The vertical modes are orthonormal with respect to the dot product f, g = 1 H 0 -H f (z)g(z)dz (the vertical integration extends from -H to 0 instead of η in this linear framework), so that we can write u q =< u, M q > , and h q = 1 ρ 0 g < p, M q >. The time evolution of u q and h q can be obtained as follows:

• (2) is multiplied by M q (z) and vertically integrated over [-H, 0] (i.e. by taking the dot product) to obtain the time evolution of u q ∂u q ∂t + g ∂h q ∂x = 0

• (3) is multiplied by N -2 (z) dM q (z) dz and vertically integrated over [-H, 0]. The result is integrated by part using [START_REF] Kamenkovich | On the time-splitting scheme used in the princeton ocean model[END_REF] and the continuity equation ( 4) is used to replace the vertical derivatives of w by the horizontal derivative of u. This leads to the time evolution of h q ∂h q ∂t + 1 gλ q ∂u q ∂x = 0

RR n°8589

We thus get a set of n uncoupled systems for 0 ≤ q ≤ n -1:

∂u q ∂t + g ∂h q ∂x = 0 ( 14
)
∂h q ∂t + 1 gλ q ∂u q ∂x = 0 (15)
This can also be expressed in terms of the characteristic variables y ± q = u q ± g c q h q , with c q = 1 λ q , which leads to n transport equations ∂y q ∂t ± c q ∂y q ∂x = 0

The c q term corresponds to the phase speed associated to vertical mode M q (in particular, c 0 is the speed of external gravity waves).

2.1 Linear stratification case (i.e. N = cste)

In this paragraph, we consider the particular case of a constant background stratification N . This simplified case will be useful for the stability analysis presented in Sec. (3.2.1). Let us introduce the dimensionless parameter ǫ = N 2 H/g whose typical value is such that ǫ ≪ 1 (Boussinesq approximation). Indeed , for N = 10 -3 s -1 and H = 4000 m, we get ǫ ≈ 4.10 -4 . The normal mode definition [START_REF] Kamenkovich | On the time-splitting scheme used in the princeton ocean model[END_REF][START_REF] Peter | The development of a free-surface bryan-cox-semtner ocean model[END_REF][START_REF] Kundu | Fluid Mechanics, Second Edition[END_REF] can be expressed in terms of ǫ (considering

λ q = c -2 q )                - d 2 M q dz 2 = ǫ g c 2 q H M q dM q dz -H = 0 dM q dz 0 = - ǫ H M q (0) (16)
For a constant background stratification, the solution of Sturm-Liouville problem (16) can be found analytically :

M q (z) = M q (0)   cos ǫg c 2 q H z - ǫc 2 q gH sin ǫg c 2 q H z   (17) tan ǫgH c 2 q = ǫc 2 q gH (18)
The eigenvalues λ q = c -2 q are supposed to be sorted in increasing order. A first approximation of the values of c q can be obtained by solving [START_REF] Nadiga | On simulating flows with multiple time scales using a method of averages[END_REF] with a rigid lid assumption dM q dz 0 = 0 for which λ 0 = 0 (c 0 tends to infinity) and c q = N H qπ are small. The knowledge of these orders of magnitude allows now to refine the solutions. The gravest mode (i.e. q = 0 ) is the barotropic Inria mode for which ǫgH c 2 0 ≪ 1 since c 0 is large. A Taylor expansion of system [START_REF] Shchepetkin | Computational kernel algorithms for fine-scale, multiprocess, longtime oceanic simulations[END_REF] in term of ǫ leads at second order to

       c 0 = α 0 gH with α 0 = 1 + ǫ 6 - ǫ 2 360 + O(ǫ 3 ) M 0 (z) = 1 -ǫ 1
M 0 (0) = 1 - ǫ 3 + 7 ǫ 2 90
The other eigenvalues (c q ) are obtained by using the fact that c q are small enough to consider that equation in ( 18) can be replaced by tan ǫgH c 2 q = 0. We can find that the baroclinic modes (i.e. for q ≥ 1) are expressed as:

q ≥ 1,        c q = α q gH with α q = √ ǫ qπ + O(ǫ 2 ) M q (z) = √ 2[cos qπ H z - ǫ qπ sin qπ H z ] + O(ǫ 2 ) (20)
Figure [START_REF] Arakawa | Vertical differencing of the primitive equations based on the charney-phillips grid in hybrid sigma-p vertical coordinates[END_REF] shows the ratio c 1 /c 0 computed from the outputs of a global 1/2 • simulation using the NEMO ocean model [START_REF] Molines | Definition of the global 1/2 • experiment with CORE forcing[END_REF]. In agreement with the theory (from ( 19) and (20) we get c 1 /c 0 ≈ √ N 2 H/(gπ)) its value is small everywhere with the exception of shallow areas where H (hence c 0 ) decreases. Note that the stability analysis will essentially use the fact that α 0 is close to 1 but the knowledge of its exact value, obtained here in the case of constant N , is not really required. In the case of non constant N , a good approximation of c 2 0 /(gH) has been given in [START_REF] Dukowicz | Structure of the barotropic mode in layered ocean models[END_REF] (eq. ( 29)) for layered ocean models. Under the hydrostatic and Boussinesq assumptions, applying the same technique as in [START_REF] Dukowicz | Structure of the barotropic mode in layered ocean models[END_REF] leads to

α 0 = c 0 √ gH = 1 + 1 gH 2 0 -H 0 z N 2 (H + z ′ )dz ′ dz
and exactly matches the first order term of our development (ǫ/6) in the case of constant N . In the following we will also use the following equalities (which hold for non constant N ):

1 H 0 -H M q (z)dz = M q (0) λ q gH (21) n-1 q=0 1 H 0 -H M q (z)dz 2 = 1 (22)
The first equality is obtained by vertically integrating the normal mode definition [START_REF] Kamenkovich | On the time-splitting scheme used in the princeton ocean model[END_REF] while the second equality is obtained by writing the decomposition of a constant velocity field u = n-1 q=0

u q M q = 1 which leads to u q =< 1, M q >= 1 H 0 -H M q (z)dz and 1 H 0 -H udz = 1 = n-1 q=0 1 H 0 -H M q (z)dz 2 . RR n°8589

Rigid lid limit

Considering that ǫ is small, the barotropic mode given by ( 19) can be approximated by M rl 0 (z) = 1 where the exponent rl stands for rigid lid. It is indeed equivalent to compute the vertical mode assuming a rigid-lid ∂η ∂t = 0 which results in a modification of the surface boundary condition in dM rl 0 (z) dz 0 = 0. Having done this assumption implies a depth independent barotropic mode for which dM rl 0 (z) dz = 0. It implies that equation ( 15) is not valid for q = 0 and in addition, instead of (21) we have for q ≥ 1, we have

0 -H M rl q (z) = 0.
This leads to values of u rl 0 (x, t) and h rl 0 (x, t) obtained by simple vertical averaging:

u rl 0 (x, t) = u(x, t) = 1 H 0 -H u(x, z, t), h rl 0 (x, t) = 1 ρ 0 g p(x, t) = 1 ρ 0 gH 0 -H p(x, z, t)
A direct integration of h rl 0 (x, t) is not possible but the pressure can be decomposed in surface and internal values using

p(x, z, t) = ρ 0 gη + p h (x, z, t) with p h (x, z, t) = 0 (23) ∂p h ∂z = -ρg (24)
so that h rl 0 = η + 1 ρ 0 gH 0 -H p h . Time evolution of the free surface η is given by the surface boundary condition [START_REF] Dukowicz | Implicit free-surface method for the bryan-coxsemtner ocean model[END_REF] :

∂η ∂t = w(z = 0)
and w at the surface is deduced from the vertical integration of the continuity equation ( 4) to get w(x, z = 0) = -H ∂u ∂x . The corresponding barotropic system is thus given by:

       ∂u ∂t + g ∂η ∂x = - 1 ρ 0 ∂ ∂x 1 H 0 -H p h dz ∂η ∂t + ∂Hu ∂x = 0 (25)
The vertically integrated internal pressure gradient evolves slowly in comparison with barotopic variables. While the use of the top boundary condition in the normal mode decomposition leads to n independent shallow water systems of type [START_REF] Molines | Definition of the global 1/2 • experiment with CORE forcing[END_REF][START_REF] Morel | Time splitting and linear stability of the slow part of the barotropic component[END_REF], using the rigid lid approximation leads to n-1 independent shallow water system of type [START_REF] Molines | Definition of the global 1/2 • experiment with CORE forcing[END_REF][START_REF] Morel | Time splitting and linear stability of the slow part of the barotropic component[END_REF] for q ≥ 1 and the system (25) which actually includes contribution from all the modes. More precisely, the vertical average of internal pressure can be decomposed as

1 H 0 -H p h dz = 1 H 0 -H q ρ 0 gh q [M q (z) -M q (0)] dz = ρ 0 g q 1 H 0 -H M q (z) -M q (0) h q (26) RR n°8589
The contribution of the barotropic mode is not zero but proportional to 1

H 0 -H M 0 (z) -M 0 (0)
which can be expressed at first order in ǫ using (19,21) as

1 H 0 -H M 0 (z) -M 0 (0) = ǫ 3
Maintaining this term constant during the barotropic integration is the main source of instability of the usual mode splitting algorithm.

3 Formulation of the stability analysis using normal mode decomposition

In this section, we explore the stability of the time splitting approach. In order to proceed to the stability analysis of the barotropic/baroclinic mode splitting approach, we will use the basis composed of the rigorous vertical modes M q (z). Note that expanding the velocity and density fields using the approximated vertical modes M rl q (z) is not possible since they do not satisfy the correct surface boundary conditions.

2D-3D correction

The idea of the mode splitting approach is to integrate separately the 2D and 3D components. When using the rigorous normal mode decomposition, the barotropic component evolves independently of the baroclinic ones so that no correction is needed after the integration of the barotropic system. Contrarily if the depth-independent approximation is made, the 3D system is corrected in order to maintain consistency with the 2D counterpart. An example of a simple time integration algorithm is depicted in figure [START_REF] Berntsen | Efficient numerical simulation of ocean hydrodynamics by a splitting procedure[END_REF]. These corrections are twofold. The first one is to modify the 3D velocities so that their vertical average matches the barotropic velocity. If u n+1 and u n+1,c denotes the 3D velocities at time n + 1 before and after correction , the correction step reads

u n ,h n ūn,2d ,η n u n+1,3d ,h n+1,3d ūn+1,2d ,η n+1 1 ρ0 @ @x 1 H Z 0 -H phdz ∆t δt = ∆t Nsplit
u n+1,c = u n+1 + ūn+1,2d - 1 H 0 z u n+1 dz

Inria

In term of projected variables u q , this correction can be written as

u n+1,c q =< u n+1,c , M q >= u n+1 q + 1 H 0 -H M q (z)dz ūn+1,2d - 1 H 0 z u n+1 dz (27)
The second correction is for consistency reason. In the 2D integration, the free surface evolves and this leads to a modification of the surface pressure gradient p(z = 0) = ρgη. But this evolution has to be in agreement with the modal decomposition which says that p(z = 0) = ρ 0 g q h q M q (0).

This second correction has to ensure that

p n+1 (z = 0) = ρgη n+1 = ρ 0 g n-1 q=0 h n+1 q M q (0) ( 28
)
where h n+1 q has evolved through the 3D integration. In the literature this correction is also known as constancy preservation (e.g. [START_REF] Griffies | Tracer conservation with an explicit free surface method for z-coordinate ocean models[END_REF][START_REF] Shchepetkin | Computational kernel algorithms for fine-scale, multiprocess, longtime oceanic simulations[END_REF]). Indeed (28) shows that the density field (hence the h q components) has to be integrated with a velocity field which has a vertical average that matches the one needed to make the free surface η evolve from n to n+1. In practice two choices are often made. The first option [START_REF] Shchepetkin | Computational kernel algorithms for fine-scale, multiprocess, longtime oceanic simulations[END_REF] is to integrate the barotropic variables and to deduce a transport which divergence matches the evolution of η from n to n + 1. This transport is then used to correct (like in eq (27)) the depth average of the 3D velocity that has been used in the integration of the density equation. The second option is to forget the 2D free surface and to recompute it from the 3D velocity fields used to advance the density field. In the next, we will follow the first option. In that case the correction can be written using ( 22)

h n+1,c q = h n+1 q + α 4 q M q (0) η n+1,2d -η n+1,3d , where η n+1,3d = n-1 j=0 h n+1 j M j (0) (29)
Indeed it leads to h n+1,c q M q (0) = h n+1 q M q (0) + η n+1,2d -η n+1,3d α 4 q M q (0) 2 . Since α 4 q M q (0) 2 = 1 using (21,22), it shows that h n+1,c q M q (0) = η n+1,3d . It is also easy to verify that it exactly corresponds to a correction of the depth average of the 3D velocity that are used to integrate the density field by a two dimensional field ū that satisfies η n+1,2d -η n,2d ∆t = -∂H ū ∂x . Note that in previous stability analysis found in the literature this correction is not taken into account even if its impact on the stability of the time splitting algorithm is important since the density field has immediately the feedback from the 2D integration.

Stability analysis using normal mode decomposition

Let X q = u q h q be the vector of velocity and pressure projections. The system (14,15) can be written as:

∂X q ∂t + A q ∂X q ∂x = 0 with A q =   0 g c 2 q g 0   RR n°8589
The correction (27,29) is written has:

X n+1,c q = X n+1 q + C q ūn+1,2d -ūn+1,3d η n+1,2d -η n+1,3d (30)
where

C q = α 2 q M q (0) 0 0 α 4 q M q (0)
The objective is now to express the vector ūn+1,2d -ūn+1,3d

η n+1,2d -η n+1,3d
according to the 2D barotropic integration. The barotropic system (25) can be rewritten as

     ∂u ∂t + g ∂ζ ∂x = 0 ∂ζ ∂t + ∂Hu ∂x = 0 (31)
where

ζ = η+ρ 0 g 1 H 0 -H
p h dz and where we have used the fact that the internal pressure gradient is maintained constant during the barotropic integration. Using A 2d the matrix corresponding to the discrete 2D time stepping algorithm, we write ūn+1,2d

ζ n+1,2d = A 2d N split ūn ζ n (32)
where N split is the number of barotropic sub time steps. The expression in term of evolution of the free surface is thus given by: ūn+1,2d

η n+1,2d = A 2d N split ūn η n + A 2d N split -I   0 1 ρ 0 gH 0 -H p h dz  
The 2D variables ū and η can be expressed as

ūn,3d = q 1 H 0 -H M q (0)u n q , η n,3d = q h n q M q (0)
As shown before (26), the decomposition of the internal pressure gradient is:

1 H 0 -H p h dz = ρ 0 g q 1 H 0 -H M q (z) -M q (0) h q (33)
Let us first assume that the internal pressure gradient is computed at time n in (33). In that case we deduce the following expression:

ūn+1,2d η n+1,2d = i V i A 2d i u n i h n i with A 2d i = I + W i A 2d N split -I , V i = α 2 i M i (0) 0 0 M i (0) , W i = 1 0 0 α 2 i

Inria

A 2d i is the matrix that makes the projections evolve during the barotropic time step. It is close to A 2d N split for the barotropic component (i = 0). In the 3D baroclinic time step, we denote by A 3d i the amplification matrix such that

u n+1,3d i h n+1,3d i = A 3d i u n i h n i so that we deduce that ūn+1,3d η n+1,3d = i V i A 3d i u n i h n i we obtain the relation ūn+1,2d η n+1,2d - ūn+1,3d η n+1,3d = i V i A 2d i -A 3d i u n i h n i (34)
Combining ( 30) and (34), we finally get:

X n+1,corrected q = X n+1,3d q + C q i W 1 i A 2d i -A 3d i X n,3d i = A 3d i X n,3d q + C q i W 1 i A 2d i -A 3d i X n,3d i (35)
The main source of instability here is in the evolution of the barotropic components X 0 = u 0 h 0 and thus we pay a particular attention on these components in the next paragraph.

Stability via diffusion at the barotropic level

We first study the stability of the 2D integration alone with respect to the barotropic variables (u 0 , h 0 ) and so look at A 2d 0 . The question is what is the minimum amount of diffusion to add to the 2D integration to compensate the fact that the depth integrated vertical pressure gradient is held constant. Let's suppose that the 2d system (32) has been integrated exactly and integrate a damping factor d 0 (a function of µ 0 ). In the following, this damping will be put either with an averaging filter or by using a diffusive 2D time stepping. A 2d corresponds to a shallow water system with a propagation speed of gH (or c 0 /α 0 ) leading to:

[A 2d ] N split = P rl 0 d 0 e -iµ0/α0 0 0 d 0 e iµ0/α0 (P rl 0 ) -1 , with (P rl 0 ) -1 =   1 α 0 g c 0 1 -α 0 g c 0  
where µ 0 = kc 0 ∆t. So that

A 2d 0 = I + W 0 P rl 0 d 0 e -iµ0/α0 -1 0 0 d 0 e iµ0/α0 -1 (P rl 0 ) -1 , with W 0 = 1 0 0 α 2 0
The matrix A 2d 0 has two complex conjugate eigenvalues and its determinant has the following simple expression :

|λ 0 | 2 = det(A 2d 0 ) = d 0 α 2 0 d 0 -(α 2 0 -1) cos µ 0 /α 0 (36)
Let us first study low frequencies (or large horizontal scales) µ 0 ≪ 1 and assume that at those scales d 0 can be developed as d 0 = 1 -γµ 2 0 . A second order Taylor expansion leads to

|λ 0 | 2 = 1 + 1 2 α 2 0 -1 α 2 0 -γ(1 + α 2 0 ) µ 2 0 ( 37
)
Since α 0 > 1, this shows that if γ = 0 (no filtering) the model is unstable at large scales. The minimum value of γ is given by

γ ≥ α 2 0 -1 2α 2 0 (α 2 0 + 1)
and (assuming constant N ) using expression (19) of α 0 and a first order expansion in ǫ we get:

γ ≥ ǫ 12 (38)
(38) is a necessary condition for the stability of the mode splitting approach. It implies that at large scales the filter has to correspond to a second order filter (and not higher order) with a minimum value given by (38).

As can be seen from expression (36), the maximum value of λ 0 are attained at each point where cos µ 0 /α 0 = -1 for which we have λ 0 = d 0 α 2 0 (d 0 + 1) -1 .

Discussion on usual filters and minimal diffusion

The preceding results gave us insight on the minimal amount of diffusion that should be present in the 2D step. The objective is now to see if this criterion is full filed by some usual filters. We consider three kinds of filters : two averaging filters with different weights (Flat and Cosine) and a diffusive 2D time stepping. The averaging filters are implemented as ūn+1,2d =

N filter m=0 a m u m,2d
where the weights a m are normalized and such that the average is centered at time n + 1 (see [START_REF] Shchepetkin | The regional oceanic modeling system (roms): a split-explicit, free-surface, topography-following-coordinate oceanic model[END_REF]):

N filter m=0 a m = 1, N filter m=0 a m m N split = 1
Note that when using averaging filters, in order to be able to center the average at time n+1, the window of integration of the barotropic equations has to extend beyond n + 1. The expression of the corresponding damping factor d 0 can easily be derived by Fourier transform:

• Flat filter over n + 1 - N filter N split , n + 1 + N filter N split , with N filter ≤ N split d Flat 0 e iµ0 = 1 2N filter + 1 m=N split +N filter m=N split -N filter e im µ 0 N split = 1 2N filter + 1 cos µ 0 N filter N split + sin µ 0 N filter N split tan µ0 2N split e iµ0 ≈ sin µ 0 N filter N split µ 0 N filter N split e iµ0 for large N split

Inria

The development of d Flat 0 at low frequencies leads to

d Flat 0 (N filter = N split ) = 1 - 1 6 µ 2 0 ≈ 1 -0.167µ 2 0 d Flat 0 (N filter = N split /2) = 1 - 1 24 µ 2 0 ≈ 1 -0.042µ 2 0 ( 39
)
• Cosine filter

d cos 0 e iµ0 = 1 N split m=3N split /2 m=N split /2 1 + cos 2π j N split e ij µ 0 N split = 2 N split (sin π N split ) 2 sin µ0 2 tan µ0 2N split cos µ0 N split -cos 2π N split e iµ0 ≈ sin(µ 0 /2) µ0 2 (1 - µ 2 0 4π 2 )
for large N split

The development of d cos 0 at low frequencies leads to

d cos 0 = 1 - 1 24 - 1 4π 2 µ 2 0 ≈ 1 -0.016µ 2 0 ( 40
)
• A dissipative Forward backward scheme. Here we don't use an averaging filter but add diffusion inside the 2D time stepping itself.

u n+1 = u n -δtg∂ x η n η n+1 = η n -δt∂ x [(1 + θ)u n+1 -θu n ]H (41)
where δt = ∆t N split . The corresponding amplification factor over N split time steps writes:

d FB 0 =   1 -θ µ 0 N split 2   N split
The development of d 0 at low frequencies leads to

d FB 0 = 1 - θ 2N split µ 2 0
So that stability condition (38) requires that

θ ≥ 2ǫN split 12 = ǫN split 6 ( 42
)
This shows that the necessary damping increases both when ǫ increases and when the splitting ratio N split increases (knowing that N split is limited by c 0 c 1 ).

Note that the stability condition of the Forward Backward scheme (41) requires (1+θ)µ 0 ≤ 2 on a A-grid and (1 + θ)µ 0 ≤ 1 on a C-grid.

To illustrate the damping effect of these filters, we fixed some parameters. First we suppose that the 2D integration is integrated with a Courant number equal to CN 0 which leads to kc 0 δt = µ 0 N split = CN 0 . N split has a maximum value which is given by the stability of the first baroclinic mode. Assuming that the maximum Courant number of the 3D time stepping is 1 (kc 1 ∆t ≤ 1), this leads to

kc 1 ∆t = µ 0 c 1 c 0 ≤ 1 ←→ N split ≤ 1 CN 0 c 0 c 1
At first order, the ratio c 0 c 1 (cf 19,20) is given by π N g H . Assuming CN 0 = 0.75, for a value of H = 4000m, we study the cases of table [START_REF] Arakawa | Vertical differencing of the primitive equations based on the charney-phillips grid in hybrid sigma-p vertical coordinates[END_REF]. As expected, the dissipative Forward Backward scheme has much less damping than the averaging filters. Indeed, comparing the developments at large scales of d 0 for the averaging filters (39,40) shows that it is much more than needed by condition (38) which leads to d 0 = 1 -ǫ 6 µ 2 0 = 1 -0.00683µ 2 0 (for the maximum value of N = 10 -2 ).

N µ 0 = c 0 c 1 ǫ N split,max

Stability with right hand side extrapolation

In the preceding paragraph, the depth integrated internal pressure gradient was the one computed at time n. There are obviously a lot of other choices (e.g. [START_REF] Shchepetkin | The regional oceanic modeling system (roms): a split-explicit, free-surface, topography-following-coordinate oceanic model[END_REF]) than can be studied in the Inria same framework. Here we look at the possibility of a time extrapolation of this term by using (1 + β)h n q -βh n-1 q in (33) (Adams Bashforth like extrapolation). In that case we get, we get ūn+1,2d

η n+1,2d = i V i A 2d i u n i h n i + β V -1 i A fb,2d N split -I W 2 i u n i -u n-1 i h n i -h n-1 i ( 43
)
with

W 2 i = 0 0 0 (α 2 i -1)M i (0)
Looking at the eigenvalues of the A 2d 0 matrix and doing a second order Taylor development of its module around µ 0 = 0 leads to

|λ 0 | = 1 - (3 -2β)(1 -α 2 0 ) + 2α 2 0 (1 + α 2 0 )γ 4α 2 0 µ 2 0
For β = 1, we recover the expression (37) and for β = 3 2 the scheme is second order accurate at low frequencies. Now we look at the global spectrum when there is no damping to see the effect of the potential extrapolation (β > 1). Here we use a perturbation analysis. We inject in the characteristic polynomial of (43) the expression λ 0 = e iµ0 (1 + ξ) and solve at first order for ξ (and at first order in ǫ) to find

ξ = 1 6 (1 -e iµ0 )(1 + β(-1 + e iµ0 ))ǫ
This leads to a module of the amplification factor equal to:

|λ 0 | = 1 - 1 3 (-1 + 2(-1 + β) cos µ 0 ) sin µ 0 2 2 ǫ
For β ≥ 1, maximum values of |λ 0 | are attained at each point where cos µ 0 = -1 (µ 0 = (2p + 1)π with p an integer). At these points we have:

|λ 0 | max = 1 + 1 3 (-1 + 2β) ǫ
and the apparent problem is that this maximum value is increased when β is increased. Figure [START_REF] Debreu | Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation[END_REF] shows this amplification factor for the case ǫ = 0.041, corresponding to N = 10 -2 (the exact amplification is not shown since it matches perfectly the approximation).

If filtering is added, the first objective will be to counter the two instabilities seen at large scales

(present if β < 3 2
) and at each aliasing frequency µ 0 = (2p + 1)π. As before, we now add a damping factor d 0 and compute the module of the eigenvalue at µ 0 = π to get

|λ 0 | max (µ 0 = π) = d 0 + 1 3 (-1 + 2β) ǫ ( 44
)
We now consider the possibility of adding second or fourth order diffusion to the 2D integration.

• Second order diffusion. Let's suppose we use a second order filter and express d 0 at second order in µ 0 1 as d 0 = 1 -γµ 2 0 . The condition at |λ 0 | max (µ 0 = π) ≤ 1 writes: For β = 3 2 , there is no instability at large scales and the previous condition leads to:

γ ≥ 1 π 2 1 3 (-1 + 2β) ǫ
γ ≥ 2 3π 2 ǫ
which is (only) slightly less restrictive that condition (38). The right hand side extrapolation allows to be second order accurate at large scales but does not really reduce the amount of required damping.

• Fourth order diffusion. Let's assume now that right hand side extrapolation has been done

(β ≥ 3 2
) so that second order diffusion at large scales is not required. Thus we can also try to apply fourth order diffusion with d 0 = 1 -γ 4 µ 4 0 and the condition |λ 0 | max (µ 0 = π) ≤ 1 leads to:

γ 4 ≥ 1 π 4 1 3 (-1 + 2β) ǫ which give for β = 3/2 γ 4 ≥ 2 3π 4 ǫ (45)
This is however a rather large biharmonic diffusion coefficient. Indeed if view in term of added diffusion in the 2D time stepping:

d 0 = 1 -γ 2d µ 0 N split α 0 4 N split
Developing for µ 0 small, leads to

d 0 = 1 - γ 4 N 3 split µ 4 0 (46)
Inria so that (45) is satisfy if

γ 2d ≥ 2 3π 4 ǫN 3 split
For the extreme case used abose of ǫ = 0.041, N split = 20 this gives γ 2d ≥ 2.2. This large value of biharmonic diffusion would severely constrain the time step of the barotropic integration.

Addition of the first baroclinic mode

We here add one baroclinic mode in the analysis, mainly to look at the damping of the first baroclinic mode associated with the 2D filtering. We assume no right hand side extrapolation β = 1. The system coming from (35) can be written as :

(X q ) n+1,corrected q=1,2 = G ((X q ) n ) q=1,2 with G =   C 0 W 1 0 A 2d 0 + I -C 0 W 1 0 )A 3d 0 C 0 W 1 1 A 2d 1 -A 3d 1 C 1 W 1 0 A 2d 0 -A 3d 0 A 3d 1 + C 1 W 1 1 A 2d 1 -A 3d 1  
The amplifications of the barotropic and baroclinic mode are plotted on figures [START_REF] Dukowicz | Structure of the barotropic mode in layered ocean models[END_REF][START_REF] Dukowicz | Implicit free-surface method for the bryan-coxsemtner ocean model[END_REF][START_REF] Griffies | Tracer conservation with an explicit free surface method for z-coordinate ocean models[END_REF]. Note that at the difference to the preceding plots, here the amplification is plotted against an inverse horizontale scale k∆x = µ 0 /(N split CN 0 ). Figure [START_REF] Dukowicz | Structure of the barotropic mode in layered ocean models[END_REF] shows the amplifications without any filtering. The barotropic mode is unstable at all scales while the baroclinic modes have small instabilities at small scales 2 . 7)) shows the same amplifications factor with added filtering for the case N = 10 -2 (resp. (N = 10 -3 )). The filtering is added either using a flat weight filter (over t, t + 2∆t) or using a 2D Forward backward diffusive time stepping (with θ = 0.13 as above).

All filters are able to stabilize the integration and as expected the dissipative Forward Backward scheme is less diffusive than the averaging filters.

2 This last point can also easily be derived by building the matrix G char of evolution of the characteristics variables Yq = y + q y - q and computing the eigenvalues of G char as perturbations of a diagonal matrix with exact amplification factors Diag[(e iµ 0 , e -iµ 0 , e iµ 1 , e -iµ 1 )] with µ 1 = kc 1 ∆t RR n°8589

Partial conclusions

The stability analysis has first confirm that the source of instability of the inexact splitting comes from the fast part of the depth integrated internal pressure gradient that is held constant during the barotropic integration. The time extrapolation of the right hand side enables to be second order at large scales for the barotropic mode but does not improve the stability. In this stability analysis framework, it is naturally sufficient to study the evolution of the (true) barotropic mode since the other baroclinic mode require much less damping.

4 Numerical implementation of a barotropic/baroclinic mode splitting and idealized test case

In this section, we will illustrate the theoretical study on a simple idealized test case that rely on the linearized (around a state at rest) primitive equations [START_REF] Berntsen | Efficient numerical simulation of ocean hydrodynamics by a splitting procedure[END_REF][START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF][START_REF] Debreu | Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation[END_REF][START_REF] Dukowicz | Structure of the barotropic mode in layered ocean models[END_REF]. We also formulate a splitting based on the use of the true depth-dependent barotropic mode. We begin by its practical implementation then show how the barotropic corrections are handled and finally show the results

Inria

of the numerical experiments.

4.1 Practical implemtation of the (depth-dependent) barotropic/baroclinic mode splitting

The integration of the barotropic part is given by eqs [START_REF] Molines | Definition of the global 1/2 • experiment with CORE forcing[END_REF][START_REF] Morel | Time splitting and linear stability of the slow part of the barotropic component[END_REF] for q = 0

∂u 0 ∂t + g ∂h 0 ∂x = 0, ∂h 0 ∂t + c 2 0 g ∂u 0 ∂x = 0 (47)
For the practical implementation, we reformulate this equation in term of u 0 , η, ρ b where

ρ b = 1 H 0 -H ρ(x, z, t)N -2 (z) dM 0 (z) dz dz
ρ b is the barotropic component of the density field and is zero when the barotropic mode is assumed to be depth-independent. One of the reason of introducing ρ b and η is to able to have access to the free surface elevation during the barotropic integration. It may be useful for example if the free surface is used in the boundary conditions. But there is a more fundamental reason. If h 0 is computed inside the barotropic integration, then at the end of this integration the correction naturally apply to the barotropic component of the pressure field. However when discretized on a vertical Lorenz grid, the presence of a computational mode ( [START_REF] Arakawa | Vertical differencing of the primitive equations based on the charney-phillips grid in hybrid sigma-p vertical coordinates[END_REF], [START_REF] Holdaway | Comparison of lorenz and charney-phillips vertical discretisations for dynamics-boundary layer coupling. part i: Steady states[END_REF]) prevents from uniquely deducing from a correction on the pressure field a correction on the density field (that is what is needed at the end).

From the decomposition of ρ (10), we can prove that:

h 0 = M 0 (0) λ 0 gH η - 1 λ 0 ρ b ρ 0 ( 48
)
Note that is this expression M 0 (0) λ 0 gH is actually equals to the vertical integral of the barotropic mode 0 -H M 0 (z)dz. Integrating (48) in (47) leads to

∂u 0 ∂t + ∂ ∂x M 0 (0) λ 0 H η - g λ 0 ρ 0 ρ b = 0, ∂ρ b ∂t - ρ 0 M 0 (0) gH ∂η ∂t - ρ 0 g ∂u 0 ∂x = 0
The time evolution of η is as usual given by the free surface condition:

∂η ∂t + ∂H ū ∂x = 0
where ū is expressed as a function of u 0 and the truly baroclinic part:

ū(x, t) = 1 H 0 -H u(x, z, t) = 1 H 0 -H [u 0 (x, t)M 0 (z) + u ′ (x, z, t)] = M 0 (0) λ 0 gH u 0 (x, t)+ 1 H 0 -H u ′ (x, z, t) 1 H 0 -H
u ′ (x, z, t) contains only baroclinic modes and thus is the slow part of the depth averaged velocity. In practice, this term is computed at the beginning of the barotropic integration and is held constant over the barotropic time steps (even if it can also be extrapolated). The final form of the barotropic equations is:

                         ∂u 0 ∂t + ∂ ∂x M 0 (0) λ 0 H η - g λ 0 ρ 0 ρ b = 0 ∂ρ b ∂t + ρ 0 M 0 (0) gH ∂H ū ∂x - ρ 0 g ∂u 0 ∂x = 0 ū = M 0 (0) λ 0 gH u 0 + 1 H 0 -H u ′ (x, z, t) ∂η ∂t + ∂H ū ∂x = 0
and we recall that when a vertically constant barotropic mode is assumed, the barotropic system writes:

       ∂ ū ∂t + g ∂η ∂x = - 1 ρ 0 ∂ ∂x 1 H 0 -H p h dz ∂ η ∂t + ∂H ū ∂x = 0

Correction of 3d variables

The correction step, that makes the barotropic mode coincides with the 3D part has been explained in (3.1) for the case of depth-independent barotropic mode. In the depth-dependent case, it includes in addition a correction of the density field and writes:

           1 H 0 -H u n+1,c M 0 = u n+1 0 1 H 0 -H ρ n+1,c N -2 dM 0 (z) dz = ρ n+1 b p n+1 (0) = ρ 0 gη n+1
The velocity correction is done as follows

u n+1,c = u n+1 + u n+1 0 - 1 H 0 -H u n+1 M 0 (z)dz M 0 (z)
while the density correction is done as follows

ρ n+1,c = ρ n+1 + 1 λ 0 -M0(0) 2 gH ρ n+1 b - 1 H 0 -H ρ n+1 N -2 dM 0 (z) dz dM 0 (z) dz

Test case

Configuration and model initialization The domain is periodic of size L x = 10000km with a depth of H = 4km. The Brünt Vaisala frequency N is taken equal to 10 -3 s -1 . The model is initialized with a barotropic solution with a right travelling wave (y - 0 = u 0 -

g c 0 h 0 = 0). h 0 (x, t = 0) = 10 sin 2π L x x e -(x-Lx/2) 2 /∆ 2 , u 0 (x, t = 0) = g c 0 h 0 (x, 0)
Inria with ∆ = 2000km. The velocity, density and free surface are then deduced:

u(x, z, t = 0) = u 0 (x, t = 0)M 0 (z), ρ(x, z, t = 0) = -ρ 0 h 0 (x, t = 0) dM 0 (z) dz , η(x, t = 0) = h 0 (x, t = 0)M 0 (0)
The initial condition ρ(x, z, t = 0) is plotted on figure (8): Results As a diagnostic of the induced damping of the different filtering techniques we look at the evolution of the maximum of the free surface elevation. Figure [START_REF] Higdon | Barotropic-baroclinic time splitting for ocean circulation modeling[END_REF] shows the results obtained with filtering techniques presented in (3.2.1). In addition to the flat and cosine weights described above, we add results obtained with the second order power law averaging filter implemented in the ROMS ( [START_REF] Shchepetkin | The regional oceanic modeling system (roms): a split-explicit, free-surface, topography-following-coordinate oceanic model[END_REF]). The reference solution corresponds to a simulation without splitting (the time step is constrained by the barotropic mode). The amount of dissipation vary quite strongly between the different filters and as expected the use of the second order power law filter reduces the damping of the free surface elevation. Of course, without any filtering, the solutions are blowing up quite rapidly. Figure [START_REF] Holdaway | Comparison of lorenz and charney-phillips vertical discretisations for dynamics-boundary layer coupling. part i: Steady states[END_REF] Barotropic/baroclinic mode splitting for free surface ocean models remains an issue to consider since it can require a large amount of unphysical diffusion to achieve a stable integration of the mode-split equations. In this paper, we introduce a framework for the stability analysis of the splitting technique. This is based on a decomposition that uses the true (depth-dependent) barotropic mode as opposed to the traditional depth-independent assumption adopted in realistic oceanic models which require extra sources of diffusion to maintain stability (either through time-filters or dissipative integration schemes). Our study reveals that the amount of diffusion induced by classical averaging filters is much larger than needed to compensante the inexact mode splitting (i.e. under the depth-independent assumption). We thus favor the use of slightly dissipative 2D time stepping algorithms. Moreover, our stability analysis allows to quantify the minimum amount of required diffusion necessary to counteract mode splitting instabilities arising linked with the depth-dependant barotropic mode approximation. The numerical experiments were done here in a very simplified model where all the assumptions of the normal mode decomposition are valid. We have however also run several realistic experiments using the ROMS model. Adding to the 2D barotropic time stepping algorithm an amount of laplacian diffusion corresponding to the theoretical study presented in this paper (and removing the existing averaging filter) has led to long terms stable runs. Using 2D dissipative time stepping algorithms has several additional advantages over the averaging filters. First the barotropic integration stops at time n + 1 and does not require additional time steps like in averaging filters. This lowers the computational cost, especially on parallel computers where the 2D integration is the less scalable part of the numerical model. It also allows to obtain a continuous free surface elevation. Finally it strongly simplifies the implementation of grid nesting with coupling at the barotropic level ( [START_REF] Debreu | Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation[END_REF]).

The formulation of a mode splitting technique that uses the depth-dependent barotropic mode has been introduced. In this formulation, the barotropic component of both velocities and density are integrated with small time steps. Several developments are still needed to be able to use this formulation in a realistic context (e.g. non flat bottom) where the normal mode decomposition in not valid and possibly where the density itself is not integrated but diagnosed from temperature and salinity. But the hope is that, even if approximated in this case, this formulation will allow the use of less diffusion than in current ocean models to maintain good stability properties.

Figure 1 :

 1 Figure 1: Ratio c 1 /c 0 computed from the outputs of a global ORCA 1/2 • simulation[START_REF] Molines | Definition of the global 1/2 • experiment with CORE forcing[END_REF]

Figure 2 :

 2 Figure 2: Time splitting algorithm

Figure 3 :

 3 Figure 3: Damping factor as the function of µ 0 for the flat and cosine averaging filters and for the dissipative Forward Backward scheme (with θ = 0.13) for N = 10 -2 (left), N = 10 -3 (right)

Figure 4 :

 4 Figure 4: Amplification factor of the barotropic mode |λ 0 | without (β = 1) and with (β = 3 2 ) right hand side extrapolation

Figure 5 :

 5 Figure 5: Barotropic and first baroclinic mode amplification without filtering for different splitting ratio N split . N = 10 -2 left, N = 10 -3 right

Flat weights, Nsplit = 10 Flat 10 k∆xFigure 6 : 2 0k∆xk∆x

 101062 Figure 6: Barotropic (left) and first baroclinic mode (right) amplification with different filters and different splitting ratio N split for N = 10 -2

Figure 7 :

 7 Figure 7: Barotropic (left) and first baroclinic mode (right) amplification with different filters and different splitting ratio N split for N = 10 -3

Figure 8 :

 8 Figure 8: Initial value of the density field

Figure 10 :

 10 Figure 10: Time evolution of the maximum of the free surface elevation. Comparaison of usual filters against a reference solution without splitting

1 2 . 2 9 3

 1229 Linear stratification case (i.e. N = cste) . . . . . . . . . . . . . . . . . . . . . . . 6 Rigid lid limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Formulation of the stability analysis using normal mode decomposition 10 3.1 2D-3D correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2 Stability analysis using normal mode decomposition . . . . . . . . . . . . . . . . 11 3.2.1 Stability via diffusion at the barotropic level . . . . . . . . . . . . . . . . 13 3.2.2 Stability with right hand side extrapolation . . . . . . . . . . . . . . . . . 16 3.2.3 Addition of the first baroclinic mode . . . . . . . . . . . . . . . . . . . . . 19 3.3 Partial conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4 Numerical implementation of a barotropic/baroclinic mode splitting and idealized test case 20 4.1 Practical implemtation of the (depth-dependent) barotropic/baroclinic mode splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.2 Correction of 3d variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.3 Test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 5 Conclusion 25 RR n°8589 RESEARCH CENTRE GRENOBLE -RHÔNE-ALPES Inovallée 655 avenue de l'Europe Montbonnot 38334 Saint Ismier Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Table 1 :

 1 Set of parameters N and corresponding values of µ

	10 -2	15	0.041	20
	10 -3	155	0.0041	207

0 , ǫ, N split,max for H = 4000m Figure

[START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF]

shows the damping factor d 0 for the two averaging filters and for the dissipative Forward Backward scheme. For this last scheme, θ is fixed by the condition (42) for the extreme case of N = 10 -2 (ǫ = 0.041, N split,max = 20) which leads to a value of θ = ǫN split,max 6 ≈ 0.13.

Two different lengths of the flat weight filter are considered with either N filter = N split (in this case the barotropic equations are integrated from t to t + 2∆t) or N filter = N split 2 .

  add two others numerical solutions. The first one corresponds to the use of a diffusive Time evolution of the maximum of the free surface elevation. Comparaison of usual filters against a reference solution without splitting 2D time stepping using a forward backward scheme with θ = 0.13 and the second one is based on the use of the true depth-dependent barotropic mode (see 4.1) without any filtering. The use of a dissipative 2D time stepping is stable and introduces much less damping than the averaging filters. Since there is no damping, the solution obtained with the true barotropic mode, matches exactly the reference solution.

	5
				Reference (no splitting)
	4			Flat filter over [t, t + 2∆t]
				Flat filter over [t +∆t/2,t+ 3∆t/2]
				Cosine Filter
				Power law Filter
	3
	2
	1
	0
	0	200	400	600	800	1000
			Time in hours
	0 1 2 4 5 Figure 9: 0 3	200	400	600 Reference (no splitting) Flat filter over [t, t + 2∆t] Flat filter over [t +∆t/2,t+ 3∆t/2] 800 Cosine Filter Power law Filter Forward Backward θ =0.13 Modal decomposition	1000
			Time in hours

RR n°8589

Here we assume that the additional filtering is achieved through a

2D dissipative time filter. Obviously µ 0 is not small at the value we are interested in (µ 0 = π) but the Taylor development is actually done in term of µ 0 /N split RR n°8589

Acknowledgements

E. Blayo, J. Demange, L. Debreu, and P. Marchesiello were funded by the ANR through contract ANR-11-MONU-005 (COMODO). F. Lemarié acknowledges the support of the French LEFE MANU and GMMC programs through project CHRONOS.


J Demange

L Debreu

P Marchesiello

F Lemarié

E Blayo

Numerical representation of internal waves propagation

Keywords: Internal gravity waves, Grid staggering, Direct space time algorithms

Similar to surface waves propagating at the interface of two fluid of different densities (like air and water), internal waves in the oceanic interior travel along surfaces separating waters of different densities (e.g. at the thermocline). Due to their key role in the global distribution of (physical) diapycnal mixing and mass transport, proper representation of internal wave dynamics in numerical models should be considered a priority since global climate models are now configured with increasingly higher horizontal/vertical resolution. However, in most state-of-the-art oceanic models, important terms involved in the propagation of internal waves (namely the horizontal pressure gradient and horizontal divergence in the continuity equation) are generally discretized using very basic numerics (i.e. second-order approximations) in space and time. In this paper, we investigate the benefits of higher-order approximations in terms of the discrete dispersion relation (in the linear theory) on staggered and nonstaggered computational grids. A fourth-order scheme discretized on a C-grid to approximate both pressure gradient and horizontal divergence terms provides clear improvements but, unlike nonstaggered grids, prevents the use of monotonic or nonoscillatory schemes. Since our study suggests that better numerics is required, second and fourth order direct space-time algorithms are designed, thus paving the way toward the use of efficient high-order discretizations of internal gravity waves in oceanic models, while maintaining good stability properties (those schemes are stable for Courant numbers smaller than 1). Finally, important results obtained at a theoretical level are illustrated at a discrete level using two-dimensional (x,z) idealized experiments.

Représentation numérique de la propagation des ondes internes

Résumé : De manière analogue aux ondes de surface se propageant à l'interface entre deux fluides de densité différente (tels que l'air et l'eau), les ondes internes dans l'océan intérieur se propagent le long des surfaces séparant des masses d'eau de différentes densités (e.g. au niveau de la thermocline). Du fait de leur rôle clé dans la distribution globale du mélange diapycnal et du transport de masses d'eau, la bonne représentation de la dynamique des ondes internes dans les modèles numérique doit devenir une priorité car les modèles globaux de climat sont maintenant intégrés à des résolutions permettant de marginalement la résoudre. Cependant, dans la plupart des modèles d'océan, les termes importants intervenant dans la propagation des ondes internes (à savoir le gradient de pression horizontal et le calcul de divergence dans l'équation de continuité) sont généralement discrétisés avec des schémas d'ordre peu élevé en espace et en temps (typiquement à l'ordre 2). Dans ce papier, nous illustrons les bénéfices d'une augmentation de la précision pour différents arrangements des variables sur la grille de calcul (i.e. grilles décalées ou colocalisées) en se basant sur la relation de dispersion au niveau discret (dans le cas linéaire). L'utilisation d'un schéma centré d'ordre 4 sur grille C pour approcher les termes de gradient de pression et de divergence horizontale apporte des améliorations significatives mais, à l'inverse des grilles colocalisées, ne permet pas l'utilisation de schémas monotones ou nonoscillants. Notre étude suggère qu'une amélioration des schémas numériques est nécessaire, nous proposons donc le développement de schémas couplés espace-temps spécifiquement destinés à l'intégration des ondes internes. Des schémas d'ordre 2 et 4 sont mis au point, ceux-ci allient de bonnes propriétés de précision et de stabilité (ils sont stables pour des nombres de Courant inférieurs à 1). Finalement, les résultats obtenus au niveau théorique sont illustrés au niveau discret à l'aide d'expériences effectuées dans un cadre semi-idéalisé.

Mots-clés : Ondes internes de gravité, Arrangement des variables, Schémas couplés espacetemps

Introduction

Internal waves are ubiquitous in the oceanic interior and have received much attention by oceanographers due to the pivotal role of internal-wave driven mixing in maintaining the ocean stratification and meridional overturning circulation. Internal waves can arise in a variety of ways, by the action of atmospheric fronts, tides, wind-stress, current-topography interactions, etc. Numerous studies have been dedicated to a better understanding of internal waves generation, propagation, dissipation, and their role on Earth's climate, either from observations or numerical simulations. In a numerical model, proper representation of internal waves dynamics is directly linked to horizontal and vertical resolution, but also to the numerical schemes used to discretize horizontal pressure gradient and horizontal divergence in the continuity equation.

Most oceanic numerical models employ low-order space-time discretization schemes to integrate internal waves dynamics. This choice has been historically made to offer the best efficiency factor (i.e. the CFL number divided by the number of evaluation of the right hand side) with less emphasis on accuracy. Indeed, for low-resolution global climate models, most of the internal wave spectrum is not resolved, hence the numerical algorithms were only chosen to maintain good stability properties (e.g. [START_REF] Brown | An economical time-differencing system for numerical weather prediction[END_REF]). However, thanks to advances in computational power, global climate models are now configured with increasingly higher horizontal/vertical resolution (it must be clear that both horizontal and vertical wavenumbers are involved in the dispersion relation). Since the typical grid-spacing in the new generation of global models is around 1/10 • (at midlatitudes), a wide part of the internal waves spectrum could be potentially represented by the corresponding computational grid. It is, thus, legitimate to reassess the historical numerical choices for the space-time integration of internal gravity waves. In particular, current low-order schemes are prone to dispersion errors and may then induce erroneous energy transfers through nonlinearities and aliasing [START_REF] Gouillon | Internal wave propagation and numerically induced diapycnal mixing in oceanic general circulation models[END_REF].

The objective of this paper is to check numerical schemes against the linear theory for which the characteristics of internal waves are well known. Many studies have emphasized the role of the grid staggering (cf. figure 1) for the proper representation of internal and external inertiagravity waves [START_REF] Dukowicz | Mesh effects for rossby waves[END_REF][START_REF] Beckers | Selection of a staggered grid for inertia-gravity waves in shallow water[END_REF]. More specifically, internal gravity waves representation can be studied using the following system of equations

∂u ∂t -f v + g ∂h ∂x = 0 ∂v ∂t + f u + g ∂h ∂y = 0 ∂h ∂t + H ∂u ∂x + ∂v ∂y = 0 (1)
where u,v are the horizontal velocity components, f is the Coriolis parameter, H the constant water depth and h the displacement. To solve (1), the criteria to discriminate the available choice of computational grids (structured or unstructured) and staggerings are

• the treatment of the Coriolis force

• the proper representation of Rossby waves

• the computation of the pressure gradient term Nowadays, the Arakawa C-grid is the preferred staggering for high resolution simulations because it is assumed that the computational mode in the averaging of the Coriolis force is controlled by viscous or diffusive terms. On the contrary, the Arakawa B-grid should remain the preferred approach for coarse resolution models because it allows a natural discretization of the Coriolis term without averaging. The A-grid has the same advantage but requires averaging to compute the pressure gradient term and the continuity equation. For structured grids, most comparisons between A, B and C-grids have been done using low order (second order) schemes. In this paper, we first show that higher order schemes show clear improvements whatever the grid staggering.

For internal waves, the discretization order is directly linked to the order of approximation of both the pressure gradient term and of the horizontal divergence term of the continuity equation.

In the following, we neglect the Coriolis term and focus on the pressure gradient term in the 2D (x, z) framework. In this (x,z) case, the C-grid will be referred to as the staggered grid and the A-grid to as the non staggered grid. The A-grid (like the B-grid in 3D) is prone to a computational mode when computing the pressure gradient term. This problem has been extensively studied in several papers in the shallow water case both for structured and unstructured grids [START_REF] Roux | Spurious inertial oscillations in shallow-water models[END_REF]. In the shallow-water case, a solution is to add a velocity-dependent dissipation (e.g. using a discontinuous Galerkin method) in order to introduce enough damping of the computational mode. In the 3D case, in the presence of internal gravity waves, this is more difficult to handle since this diffusion term should be proportional to the speed of internal gravity waves which are usually not explicitly computed (indeed, it would require the solution of a boundary value problem every time-step at each grid point). Here we propose to compute these phase speeds via a vertical normal mode decomposition, thus enabling an efficient damping of the computational mode without affecting the main signal. The computation of these speeds also enables the setup of direct space-time algorithms. It thus opens the path toward fully space time algorithms for numerical ocean models.

The paper is organized as follows. Section (2) contains a brief remainder of the normal mode decomposition and the derivation of internal gravity waves characteristics. Section (3) recalls the usual discretization properties, with a focus on higher-order approximations of the pressure gradient term and horizontal divergence, both on unstaggered and staggered grids. We also introduce an internal wave speed dependent dissipation for the damping of the computational mode on the non staggered grids. In section (4), idealized numerical experiments illustrate the theoretical results obtained in section (3). In the last section, we focus on the derivation of direct space time algorithms both on unstaggered and staggered grids, resulting in highly efficient algorithms for the numerical propagation of internal gravity waves. These schemes are also reliable for integration of external gravity waves of the barotropic part of a 3D ocean model (where the propagation speed is know in advance). Conclusions and perspectives are then drawn.

Inria 2 Normal mode decomposition

In this section, we briefly summarize salient points regarding the normal mode decomposition of the 2D (x-z) linearized (around u = u 0 ) flat bottom inviscid primitive equations (i.e. under the hydrostatic and Boussinesq assumptions). See [START_REF] Kundu | Fluid Mechanics[END_REF][START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF][START_REF] Demange | Depth-dependent barotropic integration and stability through vertical modes in free-surface ocean models[END_REF] for more details about the vertical mode decomposition. In the following, u(x, z, t) and w(x, z, t) denote the perturbation components of fluid velocities, p(x, z, t) and ρ(x, z, t) denote pressure and density perturbations around a state (p(z), ρ(z)) satisfying the hydrostatic balance dp(z) dz = -ρ(z)g where ρ(z) is a reference density profile. Using those notations, the system of interest reads

∂u ∂t + u 0 ∂u ∂x + 1 ρ 0 ∂p ∂x = 0 Momentum conservation ( 2
)
∂p ∂z = -gρ Hydrostatic equilibrium (3) ∂u ∂x + ∂w ∂z = 0 Continuity / Incompressibility (4) ∂ρ ∂t + u 0 ∂ρ ∂x + w dρ dz = 0 Energy conservation (5)
In the vertical direction the model extends from the flat bottom z = -H to the free surface z = η(x, t). Since we are interested in internal gravity waves dynamics we assume a rigid lid i.e. ∂η ∂t = 0 1 . Under this assumption, the surface and bottom boundary conditions read

∂η ∂t = 0 (6) w(z = -H) = 0 (7)
On a discrete grid, with n vertical levels, the solutions of (2,3,4,5) can be decomposed using vertical modes M q (z):

                       u(x, z, t) = n-1 q=0 u q (x, t)M q (z) p(x, z, t) = ρ 0 g n-1 q=0 h q (x, t)M q (z) ρ(x, z, t) = -ρ 0 n-1 q=0 h q (x, t) dM q (z) dz (8)
where the modes M q (z) are the eigenvectors of the following Sturm-Liouville problem 2

   ΛM q = λ q M q ∂ z M q | 0 = 0 ∂ z M q | -H = 0 (9)
1 This rigid lid assumption classically leads to a surface pressure that is diagnosed in order to get a non divergent depth integrated flow.

2 note that because of the rigid-lid assumption the top boundary condition for the Sturm-Liouville problem [START_REF] Gouillon | Internal wave propagation and numerically induced diapycnal mixing in oceanic general circulation models[END_REF] is no longer dMq dz z=0 = -N 2 (0) g Mq(0) as in [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF][START_REF] Demange | Depth-dependent barotropic integration and stability through vertical modes in free-surface ocean models[END_REF], but simply dMq dz z=0 = 0

RR n°8590

where Λ = -

d dz N -2 d dz , with N 2 (z) = - g ρ 0
dρ dz the Brünt-Vaisala frequency assumed to be positive. Each mode M q (z) is associated to a positive eigenvalue λ q = c -2 q (sorted in increasing order). The vertical modes are orthonormal with respect to the dot product f, g = 1 H 0 -H f (z)g(z)dz. The first mode (q = 0) is called the barotropic mode and is depth-independent (i.e. λ 0 = 0) under the rigid lid assumption. For constant N , an analytical expression of the modes and their associated eigenvalues can easily be found [5] M q (z)

M q (z) = √ 2 cos qπ H z , c q = N H qπ , q ≥ 1 (10)
Figure 2: Baroclinic modes M q (z) (defined in (10)) for 1 ≤ q ≤ 4 with respect to the depth z.

of equations (2,3,4,5) onto the vertical modes M q leads to the following system for q ≥ 1 :

∂ t u q + u 0 ∂u q ∂x + g∂ x h q = 0 (11)
∂ t h q + u 0 ∂h q ∂x + c 2 q g ∂ x u q = 0 ( 12
)
At this point, let us draw two important remarks

• (11) is obtained by multiplying the momentum equation ( 2) by M q and integrating over [-H, 0]. The order of accuracy for the discretization of ∂ x h q is thus directly linked to the order of accuracy of the pressure gradient discretization.

• (12) is obtained by multiplying the density equation ( 5) by N -2 (z) dM q dz , integrating over [-H, 0] and using the continuity equation. The order of accuracy for the discretization of the term ∂ x u q is linked to the order of approximation of the horizontal divergence in the continuity equation.

Inria

As mentioned earlier, the pressure gradient term and the divergence in the continuity equation are generally discretized using low-order (second-order) approximations.

In [START_REF] Demange | Depth-dependent barotropic integration and stability through vertical modes in free-surface ocean models[END_REF][START_REF] Dukowicz | Mesh effects for rossby waves[END_REF], each modal projection leads to a shallow water system which can also be expressed in terms of characteristic variables. The corresponding characteristic variables y q are y ± q = u q ± g c q h q which propagate at speed u 0 ± c q , they indeed satisfy the transport equation

∂y ± q ∂t + (u 0 ± c q ) ∂y q ∂x = 0 (13)
In practice, internal waves can be forced by barotropic tides at a given frequency ω, so that the wavenumber associated to the baroclinic mode q is given by

k q = ω c q
The main objective of this paper is to look at the numerical accuracy of the discrete propagation of these internal waves in an ocean model.

3 Discretization of internal waves in ocean models

Accuracy and grid staggering

In most ocean models, the approximation of [START_REF] Demange | Depth-dependent barotropic integration and stability through vertical modes in free-surface ocean models[END_REF][START_REF] Dukowicz | Mesh effects for rossby waves[END_REF] is done with second-order accuracy (in space and time). Indeed, both horizontal pressure gradient and horizontal divergence are computed at second order. Numerical models are generally discretized on a grid where velocity and pressure are staggered (Arakawa C-grid in the horizontal). The C-grid is superior to the non staggered Agrid at high resolution whereas the non staggered grid is advantageous for the computation of the Coriolis term, prominent at low resolution [START_REF] Griffies | Developments in ocean climate modelling[END_REF]. In the following, we consider only the pressure gradient term in a 2D x-z domain without rotation so that this problem is absent. However, we keep in mind that staggered grids are not fundamentally better especially for applications at medium and coarse resolution (e.g. climate simulations). Let us rewrite baroclinic mode q around a state at rest (u 0 = 0):

       ∂u q ∂t + g ∂h q ∂x = 0 ∂h q ∂t + c 2 q g ∂u q ∂x = 0
As mentioned earlier, we consider two staggerings : a first one analogous to A-GRID and a second one analogous to a C-GRID in 2D. The corresponding variable arrangement is represented in figure [START_REF] Brown | An economical time-differencing system for numerical weather prediction[END_REF]. Assuming that the space increment ∆x is constant, the second and fourth order approximation of spatial derivatives on the A-and C-grids are

u i u i+1 u i-1 h i+ 1 2 h i-1 2 u i u i+1 u i-1 h i-1 h i h i+1
A-GRID: ∂f (x) ∂x 2 nd order = 1 2∆x (f (x + ∆x) -f (x -∆x)) + O(∆x) 2 ∂f (x) ∂x 4 th order = 1 2∆x (f (x + ∆x) -f (x -∆x)) - 1 12∆x (f (x + 2∆x) -2f (x + ∆x) + 2f (x -∆x) -f (x -2∆x)) + O(∆x) 4 C-GRID: ∂f (x) ∂x 2 nd order = 1 ∆x (f (x + ∆x/2) -f (x -∆x/2)) + O(∆x) 2 ∂f (x) ∂x 4 th order = 1 ∆x (f (x + ∆x/2) -f (x -∆x/2)) - 1 24∆x (f (x + 3∆x/2) -3f (x + ∆x/2) + 3f (x -∆x/2) -f (x -3∆x/2)) + O(∆x) 4
The exact phase velocity is given by v exact q = ω exact k = c q while the numerical phase velocity is given by v num q = ω num k . A simple Fourier analysis (see [START_REF] Durran | Numerical Methods for Fluid Dynamics : With Application to Geophysics[END_REF]) gives the numerical phase error (4). As expected, higher order discretization of horizontal pressure gradient and horizontal divergence naturally lead to a better propagation of internal gravity waves. Indeed, we see that the phase error is significantly reduced when going from a second-order to a fourth-order approximation of ∂ x h q and ∂ x u q in (11-12), whatever the grid-staggering. In particular for the C-grid the improvement is over the whole wave spectrum while for the A grid the 2∆x wave remains unchanged (computational mode). On the A-grid, the Figure 4: Phase error with respect to the wavenumber k∆x for A and C grid for second and fourth order approximations error reduction is significant on a large part of the wave spectrum. The computational mode on the A-grid justifies the need for additional numerical diffusion which is the subject of subsection (3.2). This is less needed for the C-grid because all scales are well resolved, especially when using a fourth-order approximation.

v num /v exact indicated
Practical implementation of the fourth order discretizations As previously mentioned, better accuracy in the propagation of internal waves requires more accurate approximations of both the horizontal pressure gradient and the continuity equation. Both components have to be computed at high order since their role is equally important. The horizontal divergence should thus be computed using fourth order approximation and the result is used to deduce the vertical velocity from the continuity equation ( 4). Note that this high order computation of the horizontal divergence is also more in agreement with what is done in some numerical ocean models (e.g. NEMO, ROMS) for the nonlinear advection term (not present in our linear equations) : the velocity at mid points is computed using fourth order interpolation. In this case, using only a second order approximation in the continuity equation, introduces a problem of constancy preservation of momentum in a way identical to the constancy preservation for tracers [START_REF] Shchepetkin | The regional oceanic modeling system (roms): a split-explicit, free-surface, topography-following-coordinate oceanic model[END_REF].

Damping of the numerical mode of the nonstagerred grids

Let us note A h , A v the horizontal/vertical viscosity coefficients and K h , K v the horizontal/vertical diffusion coefficients. We allow the possibility for the A h , A v , K h , K v coefficients to be constant or dependent on the grid size and local velocity to represent second order implicit diffusion associated to odd-ordered advection schemes/time integration scheme. For vertical diffusion, we assume a simple parameterization based on the Brunt-Vaisala frequency [START_REF] Osborn | Estimates of the local rate of vertical diffusion from dissipation measurements[END_REF] which will enable us to project the equations onto the vertical modes, as follows

       ∂u ∂t + u 0 ∂u ∂x + 1 ρ 0 ∂p ∂x = A h ∂ 2 u ∂x 2 + ∂ ∂z A v N -2 (z) ∂u ∂z ∂ρ ∂t + u 0 ∂ρ ∂x + w ∂ρ ∂z = K h ∂ 2 ρ ∂x 2 + ∂ ∂z K v N -2 (z) ∂ρ ∂z (14)

RR n°8590

Using (8,9) we obtain

       ∂u q ∂t + u 0 ∂u q ∂x + g ∂h q ∂x = A h ∂ 2 u q ∂x 2 -λ q A v u q ∂h q ∂t + u 0 ∂h q ∂x + c 2 q g ∂u q ∂x = K h ∂ 2 h q ∂x 2 -λ q K v h q
which in term of characteristic variables y ± q = u q ± g c q h q leads to

∂y + q ∂t + (u 0 + c q ) ∂y + q ∂x = α h ∂ 2 y + q ∂x 2 + β h ∂ 2 y - q ∂x 2 -λ q α v y + q -λ q β v y - q ∂y - q ∂t + (u 0 -c q ) ∂y - q ∂x = α h ∂ 2 y - q ∂x 2 + β h ∂ 2 y + q ∂x 2 -λ q α v y - q -λ q β v y + q ( 15
)
where

α h = 1 2 (A h + K h ) , β h = 1 2 (A h -K h ) , α v = 1 2 (A v + K v ) , β v = 1 2 (A v -K v )
From ( 15) we can draw a few Comments:

1. We can notice that if the numerical Prandlt number A h K h is not equal to one (i.e. for β h = 0), the two characteristics are mixed together. The same holds for vertical viscosities/diffusivities (β v = 0). More generally if different advection schemes (diffusive or not) are used for density and momentum, we expect the two characteristics to become dependent one to the other in contradiction with the theory.

2. In terms of characteristics variables, Horizontal diffusion (through A h , K h ) induces a diffusion coefficient which is independent of the propagation speed3 . The exact dispersion relation for the characteristic variables is ω ± q = ±kc q . With A v , K v = 0, the dispersion relation associated to [START_REF] Griffies | Developments in ocean climate modelling[END_REF] is

ω ± q = ±k c 2 q -β 2 h k 2 -iα h k 2 q
For small values of β h such that

β 2 h k 2 ≤ c 2 q we get ω ± q = ±kc q 1 - β h k c q 2 -iα h k 2 q
which shows dispersion errors, while for larger values such that β 2 h k 2 ≥ c 2 q , it leads to purely imaginary frequency

ω ± q = -i   ±k q c q β h k c q 2 -1 + α h k 2 q  
Previous relations show that if the difference between horizontal viscosity and diffusivity is too large, some of the internal waves will stop propagating.

3. When projected onto the characteristic variables, vertical diffusion based on the Brunt-Vaisala frequency results in a friction term with a coefficient dependent on the internal wave speed (the λ q β v y ± q terms in ( 15) ). The coefficient behaves like

1 c 2 q
, it is therefore larger for high baroclinic modes. Moreover, this friction term is not a scale selective filter.

Diffusion based on internal waves phase speed The vertical mode decomposition [START_REF] Griffies | Developments in ocean climate modelling[END_REF] allows to compute an internal wave speed dependent diffusion. Since we are mainly interested in the amount of diffusion necessary to properly integrate internal waves, the idea is not to formulate the numerical model entirely in terms of the vertical modes (e.g. [START_REF] Griffiths | Internal tide generation at the continental shelf modeled using a modal decomposition: Two-dimensional results[END_REF]). Indeed, from our perspective, the projection onto the vertical modes will be done only to compute the diffusive term. This diffusive term will be then expressed back in terms of primitive variables. Any dissipative numerical scheme can be used for the propagation of internal gravity waves, like (high order) upwind schemes or monotone advection schemes. However, this choice has to be done carefully for the following reasons:

• The use of first order upwind schemes (whose diffusion coefficients, associated to the leading term in the truncation error, equal 1 2 c q ∆x) would lead to a too diffusive solution. But the vertical mode decomposition offers the possibility to build high order non oscillatory schemes acting directly on the characteristic variables that take into account the local regularity of the discrete solution. The design of non oscillatory schemes may also be required in case of internal waves breaking behaving like hydraulic jumps.

• The numerical advection schemes when formulated in terms of characteristic variables take into account the direction of propagation. This makes the problem of open boundary conditions specification well posed ( [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF]) with respect to internal waves.

• There are fundamental differences between formulation (15) based on characteristic variables and the traditional formulation [START_REF] Beckers | Selection of a staggered grid for inertia-gravity waves in shallow water[END_REF] based on primitive variables. In the diffusion based on characteristic variables, an irregularity in the density field is immediately "felt" by momentum through the added diffusion (and vice versa). This is not the case when the diffusive term are put independently on momentum and tracers, as in [START_REF] Beckers | Selection of a staggered grid for inertia-gravity waves in shallow water[END_REF].

Numerical experiment

In this section, we use idealized numerical experiments to look first at the impact of increasing the order of the discrete approximations of the pressure gradient and horizontal divergence terms (both on A and C-grid), and secondly at the use of internal wave phase speed dependent diffusion coefficients.

Design of the numerical test

The domain is a vertical section with L x = 2000km and H = 4km. The initial condition is given in term of modal decomposition and corresponds to right travelling only waves (i.e. y - q = u q -g c q h q = 0):

For 1 ≤ n ≤ n max , h n (x, t = 0) = 1 n cos(k n x)e -(x-x0) 2 /σ 2 , u n (x, t = 0) = g c n h n (x, t = 0) (16)

RR n°8590

where k n = ω c n , ω = ω M 2 = 1, 405.10 -4 s -1 (corresponding to the M2 barotropic tide), x 0 = L/5, σ 2 = 6000L x . The Brunt-Vaisala Frequency is constant and set to N = 2.10 -3 s -1 . Each mode has wavelength λ n = 2π/k n = 2πc n /ω. Let n max be the number of modes, the smallest wavelength is given by λ min = λ nmax = 2πc nmax /ω.

In the case of constant N , we have

λ min = 2πN H n max πω = 2N H n max ω
So that if we want the initial wave to be resolved by the grid we have to choose n max such that λ min ≥ 2∆x which leads to

n max ≤ N H ω∆x (17)
In the numerical experiments, we will use ∆x = 4km and ∆x = 2km. Using [START_REF] Osborn | Estimates of the local rate of vertical diffusion from dissipation measurements[END_REF] for the coarsest resolution case ∆x = 4km we get n max ≤ 2.10 -3 10 4 1, 405.10 -4 10 4 ≈ 14.3

In ( 16), n max is fixed to 14. The corresponding initial condition is plotted in figure (5). The time integration is done using a fourth-order Runge-Kutta method (referred to as RK4) to match the order of accuracy in space.

Low vs high order discretization of internal waves

In this first series of experiments, we investigate the impact of going from a second order (C2) scheme to a fourth order (C4) scheme to discretize both horizontal pressure gradient and horizontal divergence in the continuity equation. The impact of grid staggering is also investigated. Figures [START_REF] Leonard | Conservative explicit unrestricted time step multidimensional constancy preserving advection schemes[END_REF][START_REF] Griffiths | Internal tide generation at the continental shelf modeled using a modal decomposition: Two-dimensional results[END_REF] show the solutions at final time for the ∆x = 4km and ∆x = 2km resolution. Vertical lines are added to more easily see the dispersion error in comparison with the analytical solution. The corresponding relative RMS errors are indicated on figure [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF]. These simulations is constant with the theoretical dispersive errors given in figure (4). Figures [START_REF] Leonard | Conservative explicit unrestricted time step multidimensional constancy preserving advection schemes[END_REF][START_REF] Griffiths | Internal tide generation at the continental shelf modeled using a modal decomposition: Two-dimensional results[END_REF] show that fourth-order accuracy leads to a significant improvement for the A-grid at large scales. Concerning the C-grid the improvement is over the whole spectrum, as seen in figure [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF] the resulting RMS error at high resolution (∆x = 2km) is much lower than the one obtained with second order Inria schemes. This favors the use of the staggered grid. However we remind the computational mode associated with the treatment of the Coriolis term in the horizontal with this kind of staggering. Moreover, in these configurations, only the A-grid allows the design of monotone schemes for internal-gravity waves propagation since the velocity and density are collocated so that the discrete definition of the characteristic variables, which combine these two fields, is straightforward and does not require any averaging.

Diffusion scaled on internal wave phase speed

The preceding results clearly show the emergence of a numerical mode on the non staggered A-grid. The objective, here, is to add diffusion based on the internal wave phase speed to damp this computational mode. To be able to compare with the fourth order centered scheme, we have implemented a fourth order scheme with total variation diminishing (TVD) limiters. The limiters used in these simulations correspond to the minimum of the TVD region ( [START_REF] Durran | Numerical Methods for Fluid Dynamics : With Application to Geophysics[END_REF]). This Figure 7: Comparison of second (C2) and fourth (C4) order centered schemes against the analytical solution for the A and C-grid. Density field for ∆x = 2km naturally induces a velocity (c q ) dependent (nonlinear) diffusion for each internal waves. In addition, we compare with a monotony preserving (MP) scheme which is known to perform better around smooth extrema ( [4,[START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations[END_REF]. The results in figure [START_REF] Gouillon | Internal wave propagation and numerically induced diapycnal mixing in oceanic general circulation models[END_REF] show that this effectively damps the computational mode, as confirmed by the plot of the corresponding RMS error in figure [START_REF] Shchepetkin | The regional oceanic modeling system (roms): a split-explicit, free-surface, topography-following-coordinate oceanic model[END_REF]. Note also the much better behavior of the MP scheme against the TVD scheme. Even if a little of dispersion is visible for high baroclinic modes, the 2km fourth order MP solution is closed to the analytical solution. Unlike all other schemes introduced so far, it leads to a non-oscillatory scheme for the propagation of internal gravity waves.

Direct space time algorithms

Even if seldom implemented in ocean models, direct space time algorithms are sometimes used to discretize momentum advection terms (e.g. [START_REF] Adcroft | MITGCM USER MANUAL[END_REF]). The advantages of this approach are well known: it allows to construct high order schemes in space and time with minimal stencil through only one evaluation of the right hand side together with a stability condition which remains unity when the order is increased. For multidimensionnel problems, direct space time algorithms have direct extensions via a Multidimensional Advective-Conservative Hybrid Operator (MACHO) approach ( [START_REF] Leonard | Conservative explicit unrestricted time step multidimensional constancy preserving advection schemes[END_REF]) or via a directional splitting approach ( [START_REF] Bott | Improving the time-splitting errors of one-dimensional advection schemes in multidimensional applications[END_REF]). The computation of the eigenvalues of the normal mode decomposition allows the design of direct space time algorithms specifically designed for the internal waves propagation. The schemes derived in this section, both on staggered and nonstaggered grids, are basically constructed for a shallow water model, it can therefore be used in several different contexts. The objective is to build direct space time schemes for the following system of equations

∂u ∂t + g ∂h ∂x = 0 ∂h ∂t + c 2 g ∂u ∂x = 0 (18)

Non Staggered grid formulation

The formulation of space time algorithms in the scalar case is well known and can be directly extended to system (18) when the grid is not staggered. On the non staggered grid, the resulting scheme can either be formulated in terms of characteristic (as in ( 13)) or primitive variables (as in ( 18)) since there is no ambiguity in the discrete definition of the characteristic variables. We chose to formulate the scheme in term of primitive variables so that the extension to the non staggered grid will follow the same ideas. The time evolution of u i and h i is written as

u n+1 i = u n i -g ∆t ∆x F u i+1/2 -F u i-1/2 + G u i+1/2 -G u i-1/2 h n+1 i = h n i - c 2 g ∆t ∆x F h i+1/2 -F h i-1/2 + G h i+1/2 -G h i-1/2
where the interfacial fluxes are given by:

F u i+1/2 = α i-1 h n i-1 + α i h n i + α i+1 h n i+1 + α i+2 h n i+2 G u i+1/2 = β i-1 u n i-1 + β i u n i + β i+1 u n i+1 + β i+2 u n i+2 (19)
RR n°8590 For the first order hyperbolic shallow-water system (18), the obtained coefficients simply translates the fact that the even order temporal derivates of u (resp. h) is replaced by even order spatial derivatives of u (resp. h) and that odd order temporal derivates of u (resp. h) is replaced by odd order spatial derivatives of h (resp. u). As for other direct space time algorithms, both second and fourth order schemes are exact for µ = ±1.

Staggered grid formulation

For the staggered grid formulation, we start from a similar expression on a staggered grid: (h is located at i + 1/2 points, see figure (3)):

u n+1 i = u n i -g ∆t ∆x F u i+1/2 -F u i-1/2 + G u i+1/2 -G u i-1/2 h n+1 i+1/2 = h n i+1/2 - c 2 g ∆t ∆x F h i+1 -F h i-1 + G h i+1 -G h i RR n°8590

Analysis

Space time dispersion errors are plotted on figure [START_REF] Demange | Depth-dependent barotropic integration and stability through vertical modes in free-surface ocean models[END_REF]. In particular, it can be shown that for the A-grid the numerical mode is removed for 1 ≥ |µ| ≥ 1 2 4 -√ 10 ≈ 0.65.

Unlike the nonstaggered case, the direct space time schemes are not exact for µ = 1 on the staggered grid4 . But the CFL stability condition is still one. This leads to a very efficient algorithm. The stability limit of these schemes can still be increased by using classical forward backward approaches. The results obtained with these direct space time schemes (not shown) are almost identical to the ones obtained with the separated space time approach (i.e. RK4-C2 or RK4-C4) both for staggered and nonstaggered grids, but again with a far better efficiency. In the case of the A-grid, TVD or monotonicity preserving direct space time schemes can easily be derived.

Figure 13: Same as figure 11 for the fourth order direct space time algorithm on a C-grid

Conclusion and perspectives

The numerical approximation of internal waves propagation has been studied in this paper. We, first, emphasized the fact that the use of high order approximations allow much better discrete dispersion relation, even for staggered grids. In addition to the use of high order pressure gradient computation, it requires a high order computation of the horizontal divergence term in the continuity equation. Thanks to an explicit computation of internal gravity waves phase speeds, we introduced an efficient flow-dependent damping term to dissipate the non resolved modes, which is important when using non staggered grids. We also derived direct space-time algorithms that show great efficiency both in terms of stability and accuracy. Those schemes can be easily extended to satisfy monotonicity constraints for internal wave propagation on non staggered grid, however the design of such monotonic schemes remains challenging on staggered grids. The study of internal wave breaking, acting like hydraulic jumps, could probably benefit from the use of such schemes.

Figure 1 :

 1 Figure 1: Arakawa A, B and C-grids

Figure ( 2 )

 2 Figure[START_REF] Bott | Improving the time-splitting errors of one-dimensional advection schemes in multidimensional applications[END_REF] shows the first four baroclinic modes for H = 4000m and N = 2.10 -3 s -1 . Projection

Figure 3 :

 3 Figure 3: Grid staggering : A grid (top), C grid (bottom)

Figure 5 :

 5 Figure 5: Density field : initial condition

Figure 6 :

 6 Figure 6: Comparison of second (C2) and fourth (C4) order centered schemes against the analytical solution for the A and C-grid. Density field for ∆x = 4km

Figure 8 :

 8 Figure 8: Relative RMS error (in %) for ∆x = 4km (left) and ∆x = 2km (right)

Figure 10 :

 10 Figure 10: Relative RMS error with fourth order centered (C4) and fourth order TVD (C4-TVD) schemes on a A-grid

Figure 11 :

 11 Figure 11: Dissipation error (left) and dispersion error (right) with respect to the wavenumber k∆x for the fourth order direct space time algorithm on a A-grid for different values of the Courant number µ.

Table 1 :

 1 Numerical phase error for centered second and fourth order approximations on the A and C-grids responding phase errors are plotted in figure

	A-Grid	C-Grid

in table

[START_REF] Adcroft | MITGCM USER MANUAL[END_REF]

as a function of the normalized wavenumber k∆x. The cor-

RR n°8590

This is not true however if A h comes from a diffusive time stepping. As an example, a first order dissipative forward backward scheme will lead to a coefficient A h proportional to c 2 q .Inria

This essentially comes from the fact that the Euler time scheme has an error that expands in term of ∆t while the spatial derivatives have an error that expands in term of ∆x/2 instead of ∆x on a non staggered grid.RR n°8590

Acknowledgements

E. Blayo, J. Demange, L. Debreu, and P. Marchesiello were funded by the ANR through contract ANR-11-MONU-005 (COMODO). F. Lemarié acknowledges the support of the French LEFE MANU and GMMC programs through project CHRONOS.

Demange & Debreu & Marchesiello & Lemarié & Blayo

Figure 9: Numerical solutions on the A-GRID. Results at 4 and 2km with the fourth order centered scheme (C4) and with added TVD (C4-TVD) or MP limiters (C4-MP) Exchanging u and h in [START_REF] Roux | Spurious inertial oscillations in shallow-water models[END_REF] leads to the expressions of F h i+1/2 , G h i+1/2 . Fourier transform and Taylor expansion at large scales shows that second order in space and time is achieved by setting

where µ = c ∆t ∆x is the Courant number. For this particular choice of the α and β weights, the resulting scheme is equivalent to a Lax Wendroff scheme on the characteristic variables. Fourth order space time accuracy is achieved by setting

Inria

where

F h i and G h i are obtained by interchanging u and h in (20) and by adding 1/2 to the indices. Second order in space and time is achieved by setting

Setting α i-1/2 = α i+3/2 = 0 leads to an unstable scheme (the module of the amplification factor is given at large scales by |λ| = 1 + µ 2 µ 2 8 -α i+3/2 ). Here we choose a value of α i-1/2 which leads to a stability condition of |µ| ≤ 1.

The details of the computation are not given but it consists in looking at the module of the amplification factor at large (θ = k∆x ≪ 1) and small scales (θ = k∆x = π) and in finding the constraint such that it is less than one (and this can indeed be achieved only for |µ| ≤ 1 among all the choices for α i-1/2 . α i-1/2 given by ( 21) thus leads to the largest stability limit of this family). Fourth order accuracy in space and time is obtained for Figure 12: Same as figure 11 for the second order direct space time algorithm on a C-grid

Here trying to minimize the stencil by setting α i-3/2 to 0 also leads to an unstable scheme. The choice that leads to a stability condition of |µ| ≤ 1 is Again this choice of α i-3/2 comes from a study of module of the amplification factor at large (second term in the minimum above) and small scales (first term in the minimum above).

Inria