Laurent Duchemin

Christophe Josserand

Rarefied gas correction for bubble entrapment in drop impacts

Keywords: Drop impact, rarefied gas impact de gouttes, gaz raréfiés

Version française abrégée 1 Introduction

Drop impact is crucial in many multiphase flows ranging from raindrops to combustion chambers or ink-jet printing and it has become an emblematic problem of surface flows [START_REF] Edgerton | Flash[END_REF][START_REF] Rein | Phenomena of liquid drop impact on solid and liquid surfaces[END_REF]. Depending on the impact parameters (drop diameter, velocity), fluid properties (viscosity, density, surface tension) and impacted surface (liquid deep or thin film, solid substrate), it can lead to many different outcomes: spreading, rebound, prompt splash, crown splash, cavity and jet formation to cite the most famous ones [START_REF] Rioboo | Outcomes from a drop impact on solid surfaces[END_REF]. Often, the influence of the surrounding gas is neglected in the analysis because of the high density and viscosity ratios between the gas and drop liquid. Indeed, beside the entrapment of a gas bubble at impact due to lubrication effect underneath the drop [START_REF] Thoroddsen | Air entrapment under an impacting drop[END_REF][START_REF] Thoroddsen | The air bubble entrapped under a drop impacting on a solid surface[END_REF] and some aerodynamic corrections to the corolla dynamics, no significant effects of the gas was noticed so far. However, the situation has suddenly changed recently in a striking experiments on drop impacts on a smooth solid substrate [START_REF] Xu | Drop splashing on a dry smooth surface[END_REF]: there, by changing only the operating pressure, they observe that splashes were suppressed as the pressure was lower below a critical level, emphasizing thus the crucial role played by the gas in the splashing dynamics. Although different scenari has been proposed to explain such effect, involving in particular gas compressibility [START_REF] Xu | Drop splashing on a dry smooth surface[END_REF], singular bubble entrapment dynamics [START_REF] Mandre | Precursors to splashing of liquid droplets on a solid surface[END_REF][START_REF] Mani | Events before droplet splashing on a solid surface[END_REF], thin film skating [START_REF] Duchemin | Curvature singularity and film-skting during drop impact[END_REF] and film wetting dynamics [START_REF] Kolinski | Skating on a film of air: drops impacting on a surface[END_REF][START_REF] Driscoll | Ultrafast interference imaging of air in splashing dynamics[END_REF], the surrounding gas influence remains yet a vibrant question of scientific debates.

Of particular interest is the coupled dynamics between the drop and the gas underneath it just before the impact. In this case, it can be shown that the thin film air dynamics can be considered within the lubrication approximation while the liquid viscosity can be neglected as far as thin liquid jets are not formed [START_REF] Mandre | Precursors to splashing of liquid droplets on a solid surface[END_REF][START_REF] Duchemin | Curvature singularity and film-skting during drop impact[END_REF][START_REF] Korobkin | Trapping of air in impact between a body and shallow water[END_REF]. Then neglecting the surface tension one can see that a finite time singularity arises as a gas bubble is entrapped by the dynamics. This singularity behavior has to be regularized physically at least by the surface tension and the liquid viscosity but is has been argued that the resulting violent dynamics might be relevant in the splashing dynamics. Within the lubrication approximation for the gas film, a liquid jet is then formed that skates on the very thin (but non zero) gas layer. Eventually it has been shown experimentally that the liquid wets the solid substrate [START_REF] Kolinski | Skating on a film of air: drops impacting on a surface[END_REF][START_REF] Driscoll | Ultrafast interference imaging of air in splashing dynamics[END_REF], something that cannot be explained within such lubrication approximation when surface tension is present. In fact, different effects can be proposed to explain the liquid contact with the substrate when the gas layer is very thin: notably surface (Van der Waals for instance) interaction between the liquid interface and the substrate, interface fluctuations and/or surface roughness, and finite size (or rarefied gas) effect in the gas layer. In this paper, we focus on this latter case, that is when the gas layer thickness becomes of the order of the mean free path of the gas, leading firstly to a corrected lubrication equation for the thin film. In the next section, we recall the general scaling argument obtained for drop impact on a solid substrate using the classical lubrication equation for incompressible fluids. Then in section 3 we introduce the correction when the gas thickness is of the order of the mean-free path. Finally we discuss in section 4 the properties of the singularity in this case, in the absence of surface tension and we draw some perspectives for this work.

Scaling analysis

We consider the impact of a liquid drop of diameter D with a vertical velocity U on a smooth solid substrate. The liquid and surrounding gas densities are noted ρ l and ρ g , their dynamical viscosity µ l and µ g respectively, the surface tension γ. Considering the high Reynolds and Froude numbers of the drop

Re = ρ l U D µ l Bo = U √ gD
we assume that we can neglect the liquid viscosity and the gravity in the dynamics of the cushioning air film beneath the drop at impact. Therefore, considering incompressible fluids we obtain the following set of equations describing the evolution of the drop surface h(r, t) in axisymmetric geometry when it is approaching the substrate:

(∂Ω) ∂ t ϕ + 1 2 ∇ϕ 2 + p ρ l + γ ρ l κ = C(t), (1)
(∂Ω)

∂ t h = 1 12rµ g ∂ r (rh 3 ∂ r p), (2)
(∂Ω)

∂ t h = ∂ z ϕ -∂ r ϕ∂ r h, (3)
(Ω) ∆ϕ = 0, ( 4
)
where ϕ is the liquid velocity potential (the velocity in the liquid is u = ∇ϕ), and Ω is the drop volume, ∂Ω its interface. The first equation ( 1) is the Bernoulli equation valid in the liquid and written at the interface, while the last one (4) is the incompressible condition for such potential flow. The interface motion is described by the advection equation (3) while the lubrication equation ( 2) allows the determination of the pressure p(r, t) in the gas film. The 1/12 numerical prefactor in this equation was obtained by considering no slip condition for the gas velocity both on the solid substrate and on the fluid interface. This system of equation can be solved numerically using boundary integral method and we resort with numerical simulation on the interface only [START_REF] Perriot | Shape and stability of axisymmetric levitated viscous drops[END_REF][START_REF] Duchemin | Curvature singularity and film-skting during drop impact[END_REF]. Notice that lubrication is only valid where the slope of the interface is small enough and additional terms should be considered otherwise.

As the drop is approaching the solid substrate, cushioning of the gas leads to high pressure gradient beneath the drop. The high pressure created at the bottom of the drop deforms it shape so that a gas pocket is entrapped around the impact center. Simple scaling arguments can help to estimate the typical sizes of this entrapped gas bubble. Indeed, considering the time t = 0 as the time of impact in the absence of air, one can estimate the typical vertical H and horizontal R scales of the drop deformation as H(t) ∼ U t and R(t) ∼ √ DH thanks to a geometrical argument based on the intersection between the falling drop and the substrate. Introducing these scaling in the lubrication equation yields the following scaling for the lubrication pressure in the gas layer:

P l ∼ µ g R 2 U H 3 .
On the other hand, considering the impact pressure in the liquid needed to deviate horizontally a volume of liquid of typical size R (this argument was first given in [START_REF] Josserand | Droplet splashing on a thin liquid film[END_REF]), one obtain:

P i ∼ 1 R 2 d dt ρ l R 3 U ∼ ρ l U 2 D R .
Thus, the lubrication pressure P l is strong enough to deviate the liquid until a critical height H * where the liquid somehow has to touches the solid substrate, that is when P i = P l , yielding:

H * ∼ D µ g ρ l U D 2/3 = DSt 2/3
introducing the Stokes number as the inverse of a Reynolds number balancing the liquid inertia with the gas viscous effects:

St = µ g ρ l U D .
Therefore, one expects the liquid to contact the substrate for t * ∼ St 2/3 D/U entrapping a gas bubble of height H * and radius R * ∼ DSt 1/3 . Notice that in practical (experimental) situations, the Stokes number is very small: for instance for a 2 mm diameter droplet of water impacting at a velocity of 1 meter per second, we obtain St ∼ 10 -8 , so that the entrapped bubble radius is only few thousands of the drop diameter.

However, it has been observed in numerical simulations of the set of equations (1,2,3,4) that this contact would arise as a finite time singularity in the absence of surface tension (taking γ = 0 in the equation 1) [START_REF] Korobkin | Trapping of air in impact between a body and shallow water[END_REF][START_REF] Mandre | Precursors to splashing of liquid droplets on a solid surface[END_REF][START_REF] Duchemin | Curvature singularity and film-skting during drop impact[END_REF]. Such singularity exhibits a divergence of the pressure following

P ∼ h -1/2
min and of the interface curvature κ ∼ h -2 min which can be explained following [START_REF] Duchemin | Curvature singularity and film-skting during drop impact[END_REF]. Writing the set of governing equation in the frame moving radially with the geometrical intersection between the falling sphere and the substrate and using the dimensionless variables defined by

x = r -R(t) D , z = z D , t = U t D , h(x, t) = h(r, t) D , p = p ρ l U 2 , and φ = ϕ U D ,
the following system of equation is obtained:

∂ t φ -˙R∂ x φ + 1 2 ∇ φ2 + p = C(t), (5)
∂ th -˙R∂ xh = 1 12St(x + St 1/3 ) ∂ x (x + St 1/3 ) h3 ∂ x p , (6)
∂ th -˙R∂ xh = ∂ z ϕ -∂ r ϕ∂ r h, (7)
∆ φ = 0. ( 8
)
The dimensionless radius R = H/D gives for the singularity dimensionles radial velocity ˙R ∼ D/H. Thus developing the former set of equation near the singularity where Ṙ ∼ St -1/3 (we drop the ˜thereafter for the sake of simplicity) and assuming that the time derivatives are subdominant compared with the Ṙ∂ x terms, in good agreement with the numerics [START_REF] Duchemin | Curvature singularity and film-skting during drop impact[END_REF]. Then seeking for a self-similar structure for the interface near the singularity in the form:

h(x, t) = h min f ( x l(t)
), p(x, t) = P 0 (t)g( x l(t)

) and ϕ(x, t) = ϕ 0 (t)Φ( x l(t) , z l(t) )

we obtain the following system of equation at the interface at the dominant order in the self similar variable ξ = x/l(t) and χ = z/l(t):

-St -1/3 ϕ 0 (t) l(t) ∂ ξ Φ + ϕ 0 (t) 2 2l(t) 2 (∇Φ) 2 + P 0 (t)g = C(t) (9)
-St -1/3 h min (t) l(t) f = 1 12St h min (t) 3 P 0 (t) l(t) 2 (f 3 g ) ( 10
) -St -1/3 h min (t) l(t) f = ϕ 0 (t) l(t) ∂ χ Φ - h min (t) l(t) ∂ ξ Φf ( 11
)
where the prime stands for the derivative of the function to the variable ξ.

The lubrication equation [START_REF] Kolinski | Skating on a film of air: drops impacting on a surface[END_REF] gives the followinf relation:

P 0 (t) ∼ St 2/3 l(t) h min (t) 2 ,
and from the interface dynamics eq. ( 11) one can see that two situations have to be considered: -h min (t) l(t) leading to the observed numerical scalings since ϕ 0 (t) ∼ St -1/3 h min (t) that gives P 0 (t) ∼ St -2/3 h min (t)/l(t) and thus:

l(t) ∼ St -2/3 h min (t) 3/2 , P 0 (t) ∼ h min (t) -1/2 and κ ∼ St 4/3 h min (t) -2 .
and this regime is found to be valid for thick films, h min St 4/3 . -h min (t) l(t), which gives ϕ 0 (t) ∼ St -1/3 l(t) and yielding:

P 0 (t) ∼ St -2/3 , l(t) ∼ St -4/3 h min (t) 2 andκ ∼ St 8/3 h min (t) -3 .
Similarly, this regime is valid for thin films, h min St 4/3 . Finally, notice that this second regime has never been reached numerically because of the small values of the Stokes number considered.

Influence on the mean free path

It is interesting to observe that under typical experimental condition the Stokes number is very small so that the system of equations used has to be questionned. In particular, when the typical air layer becomes of the order of the mean free path λ one has to considered rarefied gas correction to the gas dynamics. This can be quantified by the Knudsen number Kn defined as the ratio between the mean free path and the typical air layer thickness which gives here: Kn = λ DSt 2/3 , and one expect rarefied gas effects to enter into account for Kn > 10 -3 . From kinetic theory, we have that the mean free path is related to the gas pressure P g through:

λ = k B T √ 2πd 2 P g ,
where T is the temperature, k B the Boltzmann constant and d the typical size of the atoms of the gas. For ambiant temperature T = 300 K and ambient pressure P g0 = 10 5 Pa on obtain the typical mean free path in the air λ 0 = 70 nm so that one can write the simple relation:

λ λ 0 = P g0 P g .
Considering typical Stokes number of the order of 10 -8 , we obtain for millimetric droplet Knudsen numbers of the order one so that correction have to be introduced in the gas dynamics. Finally, let us remark that decreasing the external gas pressure enhance this effect since it increases λ. For such Knudsen numbers, a simple way to account for the correction due to the rarefied gas context is to introduce a slipping velocity for the gas so that the no-slip boundary condition at a solid interface transforms into the Navier-slip condition:

u t = λ ∂u t ∂n ,
where n is the normal direction at the interface and u t the tangiantial velocity. This condition comes from the fact that at the level of the mean free path the no-slip condition is meaningless and cannot be imposed. Solving the Stokes equation between the solid substrate and the drop interface located at z = h(r, t) and assuming that the Navier-slip condition applies on both side we obtain the following relation for the radial velocity u in the gas layer, under the thin film approximation where the interface slope is supposed to be small:

u(z, t) = - 1 2µ g ∂ r p hz -z 2 + λh ,
so that the lubrication equation ( 2) becomes:

∂ t h = 1 12rµ g ∂ r (rh 3 ∂ r p) + 6λ∂ r (rh 2 ∂ r p) .
Therefore, for film thickness h λ one expects that the second term in the lubrication dominates so that the singularity features should be changed.

Discussion

First of all, we will assume that the position of the singularity is not strongly affected by this correction, that is we can consider the gas bubble entrapped to have the same features than before. This assumption is reasonable since the Knudsen number is of order one but since the Knudsen number increases when the pressure is lowered, the other situation where the bubble size itself is selected by the rarefied gas lubrication term should also be considered in further work. Then writing the set of equation again in the moving frame and in dimensionless units and then seeking for similarity solution using the same change of variables, we obtain at the dominant order near the singularity:

-St -1/3 ϕ 0 (t) l(t) ∂ ξ Φ + ϕ 0 (t) 2 2l(t) 2 (∇Φ) 2 + P 0 (t)g = C(t) (12)
-St -1/3 h min (t)

l(t) f = 1 12St h min (t) 3 P 0 (t) l(t) 2 f 3 g + 6 λ h min (t) f 2 g (13) -St -1/3 h min (t) l(t) f = ϕ 0 (t) l(t) ∂ χ Φ - h min (t) l(t) ∂ ξ Φf (14)
Notice that λ is now the mean free path made dimensionless using the drop diameter. Now, two regimes can also be identified due to the two terms in the right hand side of the lubrication equation ( 13).

-h min (t) λ, the pressure relation due to the dominant term in the lubrication equation remains:

P 0 (t) ∼ St 2/3 l(t) h min (t) 2 .
-h min (t) λ, the pressure relation is determined by the other term in the lubrication equation, yielding

P 0 (t) ∼ St 2/3 l(t)
λh min (t) .

Finally, we resort with two dynamical scenari for the singularity dynamics, depending on the ratio between the mean free path λ and the critical thickness St 4/3 separating the two regimes of the singularity without slip condition, namely:

-λ St 4/3 , two dynamical regimes follow. First, when h min λ, then the "usual" scaling are valide:

l(t) ∼ St -2/3 h min (t) 3/2 , P 0 (t) ∼ h min (t) -1/2 and κ ∼ St 4/3 h min (t) -2 .
then it is followed by an other regime as h min decreases. When h min λ, the following scalings are obtained: λ the "usual" scaling holds: l(t) ∼ St -2/3 h min (t) 3/2 , P 0 (t) ∼ h min (t) -1/2 and κ ∼ St 4/3 h min (t) -2 .

l(t) ∼ √ λh min (t)St -2/
It is followed by a regime where St 4/3 h min λ which gives the second scaling obtained in the beginning P 0 (t) ∼ St -2/3 , l(t) ∼ St -4/3 h min (t) 2 and κ ∼ St 8/3 h min (t) -3 .

Such regime is also folowed by a new regime when St 4/3 λ h min yielding: l(t) ∼ λh min (t)St -4/3 P 0 (t) ∼ St -2/3 and κ ∼ St 8/3 λ 2 h min (t) .

stressed by Paul in our collaboration with him on the Rayleigh-Taylor instability nonlinear dynamics [START_REF] Duchemin | Asymptotic behavior of the rayleightaylor instability[END_REF], allowing the safe used of the method improved there to the case of drop impact. Dis-moi si je dis des conneries la, c'est pour donner une note sentimentale... Question: remercie-t-on Cedric Croizet car il m'a indique comment prendre en compte la limite de gaz rarefie ou lui propose-t-on d'etre co-auteur? A discuter.

C.J. wants to acknowledge the support of the Agence Nationale de la Recherche through its Grant DEFORMATION ANR-09-JCJC-0022-01 and of the Programme Émergence(s) of the Ville de Paris.

  3 , P 0 (t) ∼ 1 √ λ and κ ∼ St 4/3 λh min (t) . and we remains thereafter within the configuration where h min (t)/l(t) ∼ St 2/3 / √ λ 1. In this case, the regime where h min (t) l(t) is not present. -λ St 4/3 , then for h min St 4/3

Acknowledgements

It is our pleasure to thank Paul Clavin for stimulating discussions on this subject initiated in a workshop organized by him in Peyresq. Moreover, we would like to emphasize that a similar boundary integral method was hardly


Carine Lucas
email: carine.lucas@univ-orleans.fr

James C Mcwilliams

Antoine Rousseau
email: antoine.rousseau@inria.fr

Modélisation Mathématique et Analyse Numérique ON NONTRADITIONAL QUASI-GEOSTROPHIC EQUATIONS

Keywords: 1991 Mathematics Subject Classification. (2000): 35Q35, 76U05, 76M45, 35B40. December 20, 2016 Ocean modeling, Coriolis force, traditional approximation, tilted quasi-geostrophic equations, slanted rotation

In this article, we work on nontraditional models where the so-called traditional approximation on the Coriolis force is removed. In the derivation of the quasi-geostrophic equations, we carefully consider terms in δ/ε, where δ (aspect ratio) and ε (Rossby number) are both small numbers. We provide here some rigorous crossed-asymptotics with regards to these parameters, prove some mathematical results and compare QHQG and QG models.

Résumé. Dans cet article, nous traitons de modèles dans lesquels l'approximation traditionnelle sur la force de Coriolis est relaxée. Dans la dérivation des équations quasi-géostrophiques, nous considérons tout particulièrement les termes en δ/ε, où δ (rapport d'aspect) et ε (nombre de Rossby) sont deux petits paramètres. Nous donnons ici des limites croisées obtenues rigoureusement en fonction de ces deux paramètres, démontrons des résultats mathématiques et comparons modèles QHQG et QG.

Introduction

Equations used to model oceans rely on several hypotheses, depending on the underlying accuracy requirements in numerical simulations. One of these hypotheses is the so-called traditional approximation, which consists in neglecting the cosine part of the Coriolis force, modifying the zonal and vertical components of the momentum equations (see f * terms in Equations (1a) and (1c)). This approximation has been widely discussed in the literature (see [START_REF] Eckart | Hydrodynamics of oceans and atmospheres[END_REF] and the correspondence in [START_REF] Phillips | The equations of motion for a shallow rotating atmosphere and the "traditional approximation[END_REF][START_REF] Phillips | Reply (to George Veronis)[END_REF][START_REF] Veronis | Comments on Phillips' proposed simplification of the equations of motion for a shallow rotating atmosphere[END_REF][START_REF] Wangsness | Comments on "the equations of motion for a shallow rotating atmosphere and the 'traditionnal approximation[END_REF]).

Discussions on the traditional approximation usually appear when studying oceans near the equator. Indeed, the Coriolis force expresses as -→ Ω × U where U is the fluid velocity, and -→ Ω = Ω t (0, cos θ, sin θ), θ being the latitude. When modeling ocean in the equatorial zone, the latitude becomes smaller and the cosine part of the Coriolis force is not negligible compared to the sine part. This is the approach developed for example in [START_REF] De Verdière | Flows in a rotating spherical shell: the equatorial case[END_REF], based on [START_REF] Veronis | Large scale ocean circulation[END_REF], or in [START_REF] Raymond | Equatorial Meridional Flows: Rotationally Induced Circulations[END_REF], where the author refers to [START_REF] White | Dynamically consistent quasi-hydrostatic equations for global models with a complete representation of the Coriolis force[END_REF] to assert that the cos θ Coriolis terms cannot be omitted. These cosine terms are traditionally not considered in the quasi-geostrophic equations that are very familiar to

Quasi-geostrophic equations

Since the QG equations are obtained in the zero-limit of the Rossby number (large rotational effects), we choose to retain all rotating terms in the original equations, before performing the asymptotic analysis. We will see in Equation ( 16) that the modified QG equations, that we will call quasi-hydrostatic quasi-geostrophic (QHQG) equations, are very similar to the traditional QG ones, except that they raise a new vertical direction (denoted Z hereafter), which differs from the traditional vertical direction z, see also [START_REF] Julien | Generalized quasi-geostrophy for spatially anisotropic rotationally constrained flows[END_REF]. The tilt between z and Z is proportional to the nondimensional parameter λ introduced in Section 1.1.1, which mesures the ratio between traditional and nontraditional Coriolis terms. Experimental and numerical evidences of this tilted vertical direction can be found in [START_REF] Gerkema | Geophysical and astrophysical fluid dynamics beyond the traditional approximation[END_REF][START_REF] Sheremet | Laboratory experiments with tilted convective plumes on a centrifuge: a finite angle between the buoyancy force and the axis of rotation[END_REF][START_REF] Straneo | Idealized models of slantwise convection in a baroclinic flow[END_REF][START_REF] Wirth | Tilted convective plumes in numerical experiments[END_REF].

Derivation of the QHQG model

In this section we present the derivation of the quasi-hydrostatic quasi-geostrophic (QHQG) equations. The derivation follows classical principles (as in [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF]): scaling, asymptotic expansion with respect to a small parameter, equations satisfied at order zero and one. Here, the small parameter (denoted ε in the sequel) is the Rossby number, so that we underline the effect of rotating terms. In order to account for the complete Coriolis force (see e.g., [START_REF] Lucas | Quasi-hydrostatic primitive equations for ocean global circulation models[END_REF] and references therein), we retain all the rotating terms in the original equations, including the terms that are usually neglected in the so-called traditional approximation. This denomination was introduced by Carl Eckart [START_REF] Eckart | Hydrodynamics of oceans and atmospheres[END_REF].

Scaling Parameters and Scaled Equations

We consider a three-dimensional domain with periodic boundary conditions in the horizontal directions, rigid lid and flat bottom in the vertical. The governing non-hydrostatic equations, including the complete Coriolis force, read:

Du Dt -f v + f * w = -∂ x ϕ + F u , (1a)
Dv Dt + f u = -∂ y ϕ + F v , (1b)
Dw Dt -f * u + g 0 = -∂ z ϕ + F w , (1c)
∂ x u + ∂ y v + ∂ z w = 0, (1d) D Dt = F . (1e)
Here (u, v, w) and are respectively the three-dimensional velocity and density of the fluid, 0 stands for the averaged density of the fluid, and ϕ is the renormalized pressure, ϕ = p/ 0 . The scalars f = 2Ω sin(θ) and f * = 2Ω cos(θ) are the Coriolis parameters where Ω stands for the angular velocity of the earth and θ is the latitude; g is the universal gravity constant. The terms F • correspond to external forces (including diffusion terms). Equations (1a) to (1c) describe the conservation of momentum, where D/Dt is the material derivative D/Dt = ∂ t +u∂ x +v∂ y +w∂ z , and Equation (1d) corresponds to the conservation of mass. Finally, Equation (1e) describes the advection of tracers (here the density ). The density and the pressure may be classically decomposed as (x, y, z, t) = ρ(z) + ρ(x, y, z, t) and ϕ(x, y, z, t) = ϕ(z) + φ(x, y, z, t),

where ρ and ϕ are the (known) background density and potential, depending only on the vertical variable. The functions ρ and ϕ are in hydrostatic balance (∂ z ϕ = -ρg). We also denote by N 2 (z) = -∂ z ρ(z) the buoyancy frequency, assuming that ∂ z ρ is negative and bounded away from zero.

Before going further in the derivation of the corresponding QG model, let us insist on the fact that we keep in Equations (1a) and (1c) the Coriolis terms f * w and f * u: we do not use the traditional approximation. We will finally lead to a slightly modified QG model (see Equation ( 16)). We think that it is a relevant modification, since the QG approximation aims at underlying the Earth's rotation effects: one should thus include every rotation term in the primitive equations prior to an asymptotic expansion with respect to the Rossby number. In the context of the β-plane approximation (see [START_REF] Dellar | Variations on a beta-plane: derivation of non-traditional beta-plane equations from Hamilton's principle on a sphere[END_REF] for example), as suggested in [START_REF] Gerkema | Geophysical and astrophysical fluid dynamics beyond the traditional approximation[END_REF] we have, with θ 0 the average latitude and R e the mean radius of the Earth:

f = f 0 + 1 R e 2Ω cos θ 0 y, with f 0 = 2Ω sin θ 0 , f * = f * 0 = 2Ω cos θ 0 .
We now introduce the following dimensionless variables, as it is classically done in QG modeling:

(x, y) = L (x , y ), z = H z , t = L U t , u = U u , v = U v , w = U H L w , ρ = P ρ and ϕ = HP gϕ , ρ = 0f0U L gH ρ , φ = f 0 U Lφ .
The Rossby number ε = U/f 0 L is the fundamental ordering parameter in the following asymptotic expansion. A secondary ordering parameter is the scale ratio of ρ to : this ratio is assumed to be ε as usual for large scale oceanography, that is1 , we assume:

0 f 0 U L gH = P ε.
Finally, the density may be expressed in terms of nondimensionnal quantities:

= P ρ (z) + ερ ,
as well as the Coriolis parameter:

f = f 0 (1 + εβ 0 y ) with β 0 = 2ΩL 2 cos θ 0 U R e .
We end this section with the non-hydrostatic scaled equations (we naturally drop the primes):

Du Dt - 1 ε + β 0 y v + λ ε w = - 1 ε ∂ x φ + ν 1h ∆ h u + ν 1z ∂ 2 zz u, (2a)
Dv Dt + 1 ε + β 0 y u = - 1 ε ∂ y φ + ν 1h ∆ h v + ν 1z ∂ 2 zz v, (2b)
δ 2 Dw Dt - λ ε u + ρ ε = - 1 ε ∂ z φ + ν 2h ∆ h w + ν 2z ∂ 2 zz w, (2c)
∂ x u + ∂ y v + ∂ z w = 0, (2d) Dρ Dt + w ε ∂ z ρ = ν 3h ∆ h ρ + ν 3z ∂ 2 zz ρ, (2e)
where we recall that ε is the Rossby number (meant to go to zero), δ = H/L is the domain aspect ratio, and λ = δ cot θ 0 . The operator ∆ h is the horizontal Laplacian. In those equations, we do not consider any external forcing, except diffusion terms with kinematic viscosities (ν 1• , ν 2• > 0) and with eddy viscosities (ν 3• > 0) that are required for geophysical flows (see e.g. [7, chapter 3, p. 39]). When we considered the scaling numbers introduced above, we have implicitly assumed that the leading term at the left-hand-side of Equation (1a) was f v, which means (when comparing the two Coriolis terms in this equation) that δ cot θ 0 should not be too large: λ 1.

(3) Fortunately, because the aspect ratio δ = H/L is rather small in large ocean models, the condition (3) is easily satisfied.

Remark 1. Actually, the term δ 2 Dw/Dt could be set to zero in (2c) above with no modification in the sequel: indeed the reader will see below that w vanishes at the main order, hence δ 2 Dw/Dt = O(ε) can be neglected prior to the QG approximation.

The final model is thus called QHQG since the differences between the new QG model and the traditional one rely only on the terms related to the Coriolis force.

The key point of the present work is to draw the reader's attention on the fact that the ratio λ/ε = cot θ 0 δ/ε is not necessarily small compared to 1/ε, and that it may have some physical impact (see Section 1.4).

Geostrophic Balance

We now consider an asymptotic expansion of all variables with respect to the Rossby number: for every unknown function γ, we write the formal asymptotic expansion γ = γ (0) + εγ (1) + ε 2 γ (2) + ... where (γ (j) ) j≥0 behave as O(1) as ε goes to zero. Equations (2a)-(2c) give, keeping only the leading order terms in ε:

-v (0) + λw (0) = -∂ x φ (0) , (4a)
u (0) = -∂ y φ (0) , (4b)
-λu (0) + ρ (0) = -∂ z φ (0) . (4c)
We complement this system with Equation (2e) (written to the order zero and since ∂ z ρ never vanishes) which leads to w (0) = 0.

Remark 2. An alternative way to obtain the main order of w is based on the incompressibility contition. It reads 0) and this traditionally leads to w (0) = 0, thanks to Equations (4a), (4b) and boundary conditions on w (see [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF]). Here, the incompressibility condition does not provide ∂ z w (0) = 0, but we have, denoting

∂ z w (0) = -∂ x u (0) -∂ y v (
∂ Z = ∂ z + λ∂ y : ∂ Z w (0) = ∂ z w (0) + λ∂ y w (0) = -∂ x u (0) -∂ y v (0) + λ∂ y w (0) = curl(∂ y φ, ∂ x φ) (5)
= 0.

Thanks to [START_REF] Charve | Convergence of weak solutions for the primitive system of the quasigeostrophic equations[END_REF] and to homogeneous boundary conditions on w (0) , we obtain again that w (0) = 0.

The geostrophic equations read:

-v (0) = -∂ x φ (0) , (6a)
u (0) = -∂ y φ (0) , (6b)
ρ (0) = -∂ z φ (0) -λ∂ y φ (0) = -∂ Z φ (0) , (6c)
w (0) = 0. (6d) 1.1.3. Quasi-Geostrophic Equations
Now we need the first order equations in order to determine the evolution of φ (0) . We denote by d g the zero-order material derivative:

d g = ∂ t + u (0) ∂ x + v (0) ∂ y .
The first order equations are:

d g u (0) -β 0 yv (0) -v (1) + λw (1) = -∂ x φ (1) + ν 1h ∆ h u (0) + ν 1z ∂ 2 zz u (0) , (7a)
d g v (0) + β 0 yu (0) + u (1) = -∂ y φ (1) + ν 1h ∆ h v (0) + ν 1z ∂ 2 zz v (0) , (7b)
-λu (1) + ρ (1) = -∂ z φ (1) , (7c)

∂ x u (1) + ∂ y v (1) + ∂ z w (1) = 0, ( 7d
)
d g ρ (0) + w (1) ∂ z ρ = ν 3h ∆ h ρ (0) + ν 3z ∂ 2 zz ρ (0) . ( 7e
)
We now take the curl of Equations ( 7a)-(7b) to obtain, thanks to Equation (7d)

d g ∂ x v (0) -∂ y u (0) -∂ z w (1) -λ∂ y w (1) + β 0 v (0) = ν 1h ∆ h ∂ x v (0) -∂ y u (0) + ν 1z ∂ 2 zz ∂ x v (0) -∂ y u (0) . ( 8
)
We notice, as for the traditional QG equations, that β 0 v (0) = d g (β 0 y). We thus try to express the quantity -∂ z w (1) -λ∂ y w (1) = -∂ Z w (1) as d g (Γ) plus viscous terms, where Γ is a function to be defined. To this aim, we will extensively make use of Equation (7e) that we reformulate:

w (1) = N -2 d g ρ (0) -ν 3h ∆ h ρ (0) -ν 3z ∂ 2 zz ρ (0) = d g N -2 ρ (0) -ν 3h N -2 ∆ h ρ (0) -ν 3z N -2 ∂ 2 zz ρ (0) . (9)
Given ( 9), we may compute the required quantity

∂ z w (1) + λ∂ y w (1) = ∂ z d g N -2 ρ (0) + λ∂ y d g N -2 ρ (0) -ν 3h (∂ z + λ∂ y ) N -2 ∆ h ρ (0) -ν 3z (∂ z + λ∂ y ) N -2 ∂ 2 zz ρ (0) . ( 10
)
We remark that for any function γ and any variable * we have the identity

(d g (γ)) * = d g (γ * ) + u (0) * ∂ x γ + v (0) * ∂ y γ,
so that we can write Equation ( 10)

∂ z w (1) + λ∂ y w (1) = d g ∂ z N -2 ρ (0) + λd g ∂ y N -2 ρ (0) -ν 3h (∂ z + λ∂ y ) N -2 ∆ h ρ (0) -ν 3z (∂ z + λ∂ y ) N -2 ∂ 2 zz ρ (0) + R, (11)
where the remainder R, according to the remark above, is

R = ∂ z u (0) ∂ x N -2 ρ (0) + ∂ z v (0) ∂ y N -2 ρ (0) + λ∂ y u (0) ∂ x N -2 ρ (0) + λ∂ y v (0) ∂ y N -2 ρ (0) = ∂ x N -2 ρ (0) ∂ z u (0) + λ∂ y u (0) + ∂ y N -2 ρ (0) ∂ z v (0) + λ∂ y v (0) .
Using (6a)-(6c) again, we have

R = N -2 ∂ x ρ (0) ∂ y ρ (0) -N -2 ∂ y ρ (0) ∂ x ρ (0) = 0,
which simplifies Equation [START_REF] Embid | Low Froude number limiting dynamics for stably stratified flow with small or finite Rossby numbers[END_REF] as follows:

∂ Z w (1) = ∂ z w (1) + λ∂ y w (1) = d g ∂ z N -2 ρ (0) + λ∂ y N -2 ρ (0) -ν 3h (∂ z + λ∂ y ) N -2 ∆ h ρ (0) -ν 3z (∂ z + λ∂ y ) N -2 ∂ 2 zz ρ (0) . ( 12
)
Back to Equation ( 8), we replace expression [START_REF] Gerkema | Geophysical and astrophysical fluid dynamics beyond the traditional approximation[END_REF] and get :

d g ∂ x v (0) -∂ y u (0) -∂ z N -2 ρ (0) -λ∂ y N -2 ρ (0) + β 0 y = ν 1h ∆ h ∂ x v (0) -∂ y u (0) + ν 1z ∂ 2 zz ∂ x v (0) -∂ y u (0) -ν 3h (∂ z + λ∂ y ) N -2 ∆ h ρ (0) -ν 3z (∂ z + λ∂ y ) N -2 ∂ 2 zz ρ (0) . ( 13
)
Let us now rewrite Equation ( 13), expressing everything in terms of φ (0) . We have

∂ t -∂ y φ (0) ∂ x + ∂ x φ (0) ∂ y ζ = ν 1h ∆ h ∆ h φ (0) + ν 1z ∂ 2 zz ∆ h φ (0) +ν 3h (∂ z + λ∂ y ) N -2 ∆ h (∂ z φ (0) + λ∂ y φ (0) ) + ν 3z (∂ z + λ∂ y ) N -2 ∂ 2 zz (∂ z φ (0) + λ∂ y φ (0) ) , (14)
ζ = ∆ h φ (0) + N -2 (∂ z + λ∂ y ) 2 φ (0) + ∂ z N -2 (∂ z + λ∂ y ) φ (0) + β 0 y. ( 15
)
We could write this more compactly for Φ = φ (0) with ∂ Z = ∂ z + λ∂ y and get the following viscous quasihydrostatic quasi-geostrophic (QHQG) equation:

(∂ t -∂ y Φ∂ x + ∂ x Φ∂ y ) ∆ h Φ + ∂ Z N -2 ∂ Z Φ + β 0 y = ν 1h ∆ h ∆ h Φ + ν 1z ∂ 2 zz ∆ h Φ + ν 3h ∂ Z N -2 ∆ h ∂ Z Φ + ν 3z ∂ Z N -2 ∂ 2 zz ∂ Z Φ . ( 16
)
One can thus easily recognize the traditional QG equation (see Equation (2.23) in [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF] without viscosity or Equation (2.7) in [START_REF] Bennett | The dissipative quasigeostrophic equations[END_REF]), except that the differential operator ∂ z is replaced by the tilted one ∂ Z = ∂ z + λ∂ y , see [START_REF] Julien | Generalized quasi-geostrophy for spatially anisotropic rotationally constrained flows[END_REF]. We recall here that λ = δ cot θ 0 is proportional to the domain aspect ratio. In particular, we recover the traditional QG equation when setting λ = 0 in Equations ( 16). These equations are also called quasigeostrophic equations for slanted rotation in [START_REF] Embid | Low Froude number limiting dynamics for stably stratified flow with small or finite Rossby numbers[END_REF]. Equation ( 16) complemented by adapted initial conditions and, eventually, boundary conditions, constitutes the QHQG model that we will study in the following.

Existence of solutions to the inviscid QHQG

In this section, we are interested in existence results for the inviscid quasi-hydrostatic quasi-geostrophic equation (Equation [START_REF] Lucas | New developments and cosine effect in the viscous Shallow-Water and quasi-geostrophic equations[END_REF] with vanishing viscosities ν 1h , ν 1z and ν 3h , ν 3z ):

(∂ t -∂ y Φ∂ x + ∂ x Φ∂ y ) ∆ h Φ + ∂ Z N -2 ∂ Z Φ + β 0 y = 0. ( 17
)
We consider a periodic domain in the horizontal variables with rigid boundaries at the top and bottom. More precisely, we denote by

B = Σ × (0, h) = -1 2 , 1 2 × -1 2 , 1 2
× (0, h) the spatial domain where we establish existence of solutions, with t ∈ [0, T ] (where T > 0) the time interval.

We assume the solutions to be horizontally periodic with period 1. The pressure potential Φ is defined by the QHQG equation ( 17) up to a constant at each time; then we add the following condition: B Φ = 0. We define u u u = (u (0) , v (0) ) the horizontal velocity at the leading order, we omit the (0) superscripts and rewrite the QHQG equation together with the geostrophic equations as:

v = ∂ x Φ, (18)
u = -∂ y Φ, (19)
ρ = -∂ z Φ -λ∂ y Φ, (20)
ω := ∂ x v -∂ y u -∂ z N -2 ρ -λ∂ y N -2 ρ in B, (21)
∂ t ω + u u u • ∇ω = -β 0 v in B × [0, T ], (22)
with the initial condition ω(x, y, z, 0) = ω 0 (x, y, z),

and the boundary conditions ρ = 0 at z = 0 and z = h.

Existence results are based on the following property:

Lemma 1. The partial differential operator L defined by

LΦ := ∆ h Φ + ∂ z N -2 (∂ z Φ + λ∂ y Φ) + λ∂ y N -2 (∂ z Φ + λ∂ y Φ) in B
is strictly elliptic.

Proof. The function N -2 (z) is assumed to be bounded away from 0 for z ∈ [0, h].

We expand L as:

LΦ = ∂ 2 xx Φ + 1 + λ 2 N -2 ∂ 2 yy Φ + N -2 ∂ 2 zz Φ + 2λN -2 ∂ 2 yz Φ + λ∂ z N -2 ∂ y Φ + ∂ z N -2 ∂ z Φ = 3 i=1 3 j=1 a ij ∂ 2 i,j Φ + 3 i=1 b i ∂ i Φ,
where:

A = (a ij ) i,j=1,...,3 =   1 0 0 0 1 + λ 2 N -2 λN -2 0 λN -2 N -2   , and b = (b i ) i=1,...,3 = 0, λ∂ z (N -2 ), ∂ z (N -2 ) .
For all ξ = (ξ 1 , ξ 2 , ξ 3 ) , we have, for almost every z in [0, h]:

3 i=1 3 j=1 a ij ξ i ξ j = ξ 2 1 + 1 + λ 2 N -2 ξ 2 2 + N -2 ξ 2 3 + 2λN -2 ξ 2 ξ 3 ≥ ξ 2 1 + 1 2 ξ 2 2 + 1 2λ 2 + N 2 ξ 2 3 (thanks to Young inequality) ≥ θ|ξ| 2 with θ := min 1 2 , 1 2λ 2 + N 2 > 0.
The condition θ > 0 is satisfied as soon as λ > 0. In the case where λ = 0 (traditional equation), one needs to assume that N -2 is bounded away from zero (see [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF]). Lemma 1 generalizes the ellipticity result of [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF] for any λ ∈ R and we can now obtain the same existence results. Indeed, System (18)-( 24) only differs from the corresponding system in [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF] in Equations ( 20) and [START_REF] Phillips | Reply (to George Veronis)[END_REF].

We first obtain the local existence of solutions: Proof. The proof follows the lines of [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF], considering an iterative process:

• starting from ξ 0 (x, y, z, t) = ω 0 (x, y, z) given,

• for k ≥ 0 compute Φ k thanks to the relation:

           ∂ 2 xx Φ k + ∂ 2 yy Φ k + ∂ z ∂ z Φ k + λ∂ y Φ k N 2 + λ N 2 ∂ y ∂ z Φ k + λ∂ y Φ k = ξ k in B, ∂ z Φ k + λ∂ y Φ k = 0 for z = 0 and z = h, B Φ k = 0, (25)
• then set

u k = -∂ y Φ k and v k = ∂ x Φ k , • for every z ∈ [0, h], compute ξ k+1 solution of ∂ t ξ k+1 + u k , v k • ∇ξ k+1 = -β 0 v k ,
with ξ k+1 (t = 0) = ω 0 , with periodic boundary conditions in x and y.

The elliptic result (Lemma 1) still gives the estimate: 18), ( 19) and ( 22) being unchanged, the short-time existence of the k-th vorticity iterate ξ k follows, as well as the upper bound on its H s norm.

|Φ k | s+2 ≤ C 0 |ξ k | s with C 0 = C 0 (B, N (z), λ) > 0 such that, Equations (
The end of the proof, for the convergence of the iterates {ξ k (t)} k≥0 , is not modified by the λ coefficient.

We are also able to adapt the proof of global solutions for the QG model by [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF] to get: Proof. From Lemma 1, relations ( 18), ( 19), [START_REF] Phillips | The equations of motion for a shallow rotating atmosphere and the "traditional approximation[END_REF] and from

           ∂ 2 xx Φ + ∂ 2 yy Φ + ∂ z ∂ z Φ + λ∂ y Φ N 2 + λ N 2 ∂ y (∂ z Φ + λ∂ y Φ) = ω in B, ∂ z Φ + λ∂ y Φ = 0 for z = 0 and z = h, B Φ = 0,
we have the elliptic estimate:

|u u u| s+1 + |ρ| s+1 ≤ C|ω| s ,
such that we also have a global estimate on the QHQG vorticity ω of the form

|ω(t)| s ≤ K(t) ∀t ∈ [0, T ],
where the function K only depends on the H s norm of the initial condition ω 0 . Thanks to the short-time existence, the global existence can be proved, using the global bound

|ω(t)| s ≤ max t∈[0,T ] K(t) ∀t ∈ [0, T ]
to reach any time T thanks to an iterative process.

Existence of solutions to the viscous QHQG

As in Section 1.2, we are interested here in adapting existence results of [START_REF] Bennett | The dissipative quasigeostrophic equations[END_REF] to the QHQG equation ( 16). Let us consider the QHQG equation with vanishing vertical viscosities ν 1z and ν 3z :

(∂ t -∂ y Φ∂ x + ∂ x Φ∂ y ) ∆ h Φ + ∂ Z N -2 ∂ Z Φ + β 0 y = ν 1h ∆ h ∆ h Φ + ν 3h ∂ Z N -2 ∆ h ∂ Z Φ , (26)
in the spatial domain B = (0, 1) × (0, 1) × (0, h), for t ∈ [0, T ] with the following boundary conditions at the top and at the bottom:

∂ z Φ = 0 for z = 0 and z = h, (27)
and with periodic boundary conditions (with period 1) for the horizontal variables. The initial condition on the vorticity ω given by ω

= ∆ h Φ + ∂ Z N -2 ∂ Z Φ is: ω(x, y, z, 0) = ω I (x, y, z), (28)
where ω I is periodic with period 1 in the horizontal directions.

In the following we denote by f ∈ C n (D) a function that is 1-periodic in x and y and with all its derivatives (with respect to variables in D) ∂ • f continuous on D for 0 ≤ ≤ n. Proof. The proof is exactly the one of [START_REF] Bennett | The dissipative quasigeostrophic equations[END_REF], based on a Schauder fixed point for the existence and energy estimates for the uniqueness. As in Section 1.2, the only novelty for the QHQG model ( 26)-( 28) is the ∂ Z term in the vorticity instead of ∂ z : Lemma 1 ensures the ellipticity of the corresponding operator and the rest of the proof remains unchanged. We refer to [START_REF] Bennett | The dissipative quasigeostrophic equations[END_REF] for details and for the expression of T * in terms of β.

Simple physical properties of the QHQG model

We detail here some basic physical properties that can be obtained for the QHQG model and we perform comparisons with the well-known QG system.

We recall the QHQG equation ( 16):

(∂ t -∂ y Φ∂ x + ∂ x Φ∂ y ) ∆ h Φ + ∂ Z N -2 ∂ Z Φ + β 0 y = ν 1h ∆ h ∆ h Φ + ν 1z ∂ 2 zz ∆ h Φ + ν 3h ∂ Z N -2 ∆ h ∂ Z Φ + ν 3z ∂ Z N -2 ∂ 2 zz ∂ Z Φ ,
where ∂ Z = ∂ z + λ∂ y and λ = δ cot θ 0 .

Coordinate transformation

We define a coordinate transform from (x, t) = (x, y, z, t) to (X, T ) = (X, Y, Z, T ) by

X = x, Y = y -λz, Z = z, T = t. ( 29
)
The transformed QHQG system is isomorphic to ( 16) because (∂ t -∂ y Φ∂ x + ∂ x Φ∂ y ) and ∆ h , have identical functional forms in the transformed coordinates; N 2 (z) = N 2 (Z); and βy = β(Y + λZ) can be replaced by βY because (∂ t -∂ y Φ∂ x + ∂ x Φ∂ y ) (λZ) = 0. Thus, all QG solutions Φ(x, t) are also QHQG solutions Φ(X, T ) if the initial and boundary conditions are consistent. Initial conditions are equivalently specified at t = 0 or T = 0. Planar vertical boundary surfaces at z = z 0 are also planar at Z = z 0 . The same is true in x and X.

Only in y is a planar surface in y no longer planar in Y ; therefore, the simple equivalence is for solutions that have an unbounded y domain, or are y-periodic, or decay away in y before the boundary.

Separable solutions and vertical modes

In the inviscid case (see Equation ( 17)), QHQG has vertically separable solutions when linearized by ne-

glecting (∂ x Φ∂ y -∂ y Φ∂ x ) ∆ h Φ + ∂ Z N -2 ∂ Z Φ , even for general N (z), Φ = Π(X, Y, T )F (Z)
. (QG has an analogous property.) Such QHQG solutions have a vertically upward phase tilt in the (y, z) plane with a slope of dz/dy = λ -1 relative to the z axis. This direction is aligned with the full rotation vector. Many phenomena influenced by the non-traditional approximation (e.g., convection and centrifugal instability as reviewed in [START_REF] Gerkema | Geophysical and astrophysical fluid dynamics beyond the traditional approximation[END_REF]) are known to exhibit this type of phase tilt. If we assume constant density at the vertical boundaries, F is determined by the 1D eigenvalue problem,

∂ Z (N -2 ∂ Z F ) + R -2 F = 0 , ∂ Z F = 0 at Z = Z 0 , Z 1 . ( 30
)
The eigenvalue R is the horizontal radius of deformation. With N constant, the eigenmodes F are cosine functions in Z.

Thermal wind balance

The geostrophic equations (6a)-(6c) have a steady solution for a zonal flow. This can be expressed as

u = U (Y, Z), = - Z 0 N 2 (Z ) dZ + B(Y, Z) with ∂ Z U = -∂ Y B,
∂ t ∆ h Φ + (∂ z + λ∂ y ) N -2 (∂ z + λ∂ y )Φ + β∂ x Φ = 0. (31)
With constant N in either a vertically bounded or unbounded domain, eigenmodes are

Φ ∝ e i(kx+ly+mz-ωt) , (32)
with a dispersion relation,

ω = -βk K 2 , K 2 = k 2 + l 2 + 1 N 2 (m + λl) 2 . ( 33
)
We can equivalently write this as Φ ∝ e i(kX+lY +M Z-ωT ) , (34) with M = m + λl.

Let us define the function

ρ k,l,m (λ) = ω(λ) ω (0)
, that will evidence the difference between QG and QHQG models. We plot in Figure 1 below the behavior of ρ k,l,m (λ) as a function of λ for various values of k, l, m.

Naturally, for θ 0 = 0 we have a singularity (λ = +∞) and ρ(λ) goes to zero. This advocates for the fact that QHQG and QG have drastically different behaviors "at the equator", where QG is known to be irrelevant. The results in Figure 1 show that the barotropic mode m = 0 has a significantly lower frequency even for non-infinite (reasonable) values of λ. Indeed, assuming l > 0 (otherwise ρ = 1) a difference of 10% can be seen with λ = 4.10 -3 with (k, l, m) = (0, 1, 0), which corresponds (with δ = 10 -3 which is a good typical value in the North Atlantic ocean) to θ 0 = 14 o . This is not that close from the equator. For the mode (k, l, m) = (1, 1, 0) the 10% margin leads to θ 0 = 11 o (λ = 5.10 -3 ). Even if this 10% difference cannot be seen for m > 0 (see Figure 1), the barotropic mode is the more energetic one and this is another reason to claim that QH modeling should be used as much as possible.

Vortex solutions

A simple vortex solution is a nonlinear stationary state when β = 0 and N is constant. In QG this occurs for any axisymmetric profile, Φ(r, z), where r = x 2 + y 2 is the radial coordinate, which is sufficient to make (∂ x Φ∂ y -∂ y Φ∂ x ) ∆ h Φ + ∂ z N -2 ∂ z Φ = 0 even when the arguments have large amplitude. Physical interest usually lies in profiles that are localized in r, e.g., a monopole with Φ(r) decaying away from a central extremum. In QHQG nonlinear stationary solutions exist for any profile Φ(R, Z), where R 2 = X 2 + Y 2 . Thus, QH vortices are meridionally tilted rather than verically aligned. This tilted structure was previously proposed for Meddy vortices by [START_REF] Semenova | Dynamically equilibrium shape of intrusive vortex formations in the ocean[END_REF].

1.4.6. Fourier-space estimates for non-traditional effects Assuming a wavenumber-space characterization of the solution, as commonly done for turbulent flows, in terms of vertical and horizontal wavenumbers, kv and kh, or equivalently local values of H and L, we can ask when the λ value is not small. In QG the common view (sometimes called Charney's stretched isotropy) is that the Burger number, Bu = N H/2ΩL = kh/(2Ωkv/N ), is order one while 2Ω/N is small. If H/L ∼ 2Ω/N 1, then λ ∼ (2Ω/N ) cot θ 0 , which will be small except when θ 0 is very close to the Equator. Alternatively, if we consider flow patterns with 2Ωkv/N kh ∼ r 1, then λ ∼ (2Ω/N r) cot θ 0 . For small enough r, λ need not be small, and the QH correction to QG will be important. This happens when the aspect ratio H/L is large, the stratification is weak, and/or θ 0 is small. Large aspect ratio can be described from the QG perspective as atypically "tall" flows.

Place of the nontraditional models

The objective of this section is to show that the quasi-hydrostatic quasi-geostrophic equation ( 16) is an intermediate model between the non-hydrostatic equations and the usual (traditional) quasi-geostrophic equation, namely Equation ( 16) with λ = 0 (or equivalently Z = z).

Hierarchy of models

We can represent these models on the diagram of Figure 2.

3 λ = 0 ε → 0 ε → 0 λ = 0 1 2 0 Figure 2.
Non-hydrostatic equations and quasi-geostrophic equations: models with and without nontraditional terms. NH: non-hydrostatic equations, see [START_REF] Bennett | The dissipative quasigeostrophic equations[END_REF]; NHT: non-hydrostatic traditional equations, see [START_REF] Bennett | The dissipative quasigeostrophic equations[END_REF] without the f * terms; QHQG: quasi-hydrostatic quasi-geostrophic equation, see [START_REF] Lucas | New developments and cosine effect in the viscous Shallow-Water and quasi-geostrophic equations[END_REF]; QG: (traditional) quasi-geostrophic equation, see [START_REF] Lucas | New developments and cosine effect in the viscous Shallow-Water and quasi-geostrophic equations[END_REF] with λ(= δ cot θ 0 ) = 0.

The convergence of NH towards QHQG (see 1 in Figure 2) is proved in [START_REF] Embid | Low Froude number limiting dynamics for stably stratified flow with small or finite Rossby numbers[END_REF]: based on techniques developed in [START_REF] Embid | Averaging over fast gravity waves for geophysical flows with arbitrary potential vorticity Communications in Partial Differential Equations[END_REF], the authors prove that the solution of NH weakly converges, when ε tends to zero, to a limit that satisfies QHQG. The convergence cannot be in the strong sense due to fast oscillations that can appear in NH and that disappear in QHQG.

The convergence of NHT (that are exactly NH with λ = 0) to QG (that is 2 in Figure 2) have been studied in several papers; let us mention [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF] in a periodic domain in the horizontal variables, with regular initial data, or in [START_REF] Charve | Convergence of weak solutions for the primitive system of the quasigeostrophic equations[END_REF] in a three dimensional periodic domain, for ill prepared data.

In the following of this section, we are interested in the right part 3 of the diagram 2: we prove that a solution of QHQG converges to a solution of QG when δ (and consequently λ) tends to zero.

Convergence as δ (or λ) tends to zero

We are interested in the limit when the aspect ratio δ tends to zero, which means that λ also tends to zero. More precisely, we prove that the solution of the quasi-hydrostatic quasi-geostrophic model ( 16) converges to its corresponding hydrostatic model when the aspect ratio δ tends to 0, that is, as λ = δ cot θ 0 , when the nontraditional Coriolis parameter λ tends to 0. We assume ν 1z = 0 and ν 3z of order ε to be consistent with boundary conditions (see [START_REF] Macwilliams | A note on a consistent quasigeostrophic model in a multiply connected domain[END_REF]).

Let us denote by φ λ the solution of the nontraditional (quasi-hydrostatic) quasi-geostrophic model for a fixed λ:

∂ t -∂ y φ λ ∂ x + ∂ x φ λ ∂ y ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ + β 0 y = ν 1h ∆ h ∆ h φ λ + ν 3h ∂ Z N -2 ∆ h ∂ Z φ λ + ν 3z ∂ Z N -2 ∂ 2 zz ∂ Z φ λ . ( 35
)
in the domain D (where ∂ Z = ∂ z + λ∂ y ), with periodic conditions in x and y, and verifying ∂ Z φ λ = 0 on the top (z = h) and bottom (z = 0). At t = 0, the function is given by φ λ 0 . We assume the data to be well-prepared, that is φ λ 0 → φ 0 when λ → 0. We want to prove that φ λ converges when λ → 0 towards a function φ, solution of the following traditional (hydrostatic) quasi-geostrophic model:

∂ t -∂ y φ (G) ∂ x + ∂ x φ (G) ∂ y ∆ h φ (G) + N -2 ∂ 2 zz φ (G) + ∂ z N -2 ∂ z φ (G) + β 0 y = ν 1h ∆ h ∆ h φ (G) + ν 3h ∂ z N -2 ∆ h ∂ z φ (G) + ν 3z ∂ z N -2 ∂ 2 zz ∂ z φ (G) , (36)
with ∂ z φ (G) = 0 on the top and bottom, and φ (G) (t = 0) = φ 0 .

2.2.1.

A priori estimates for φ λ .

Let us multiply equation ( 35) by φ λ and integrate over the periodic domain (in x and y) D:

D ∂ t ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ φ λ = ν 1h D ∆ h ∆ h φ λ φ λ +ν 3h D ∂ Z N -2 ∆ h ∂ Z φ λ φ λ +ν 3z D ∂ Z N -2 ∂ 2 zz ∂ Z φ λ φ λ
(due to the properties of the jacobian), which leads to:

- 1 2 d dt ∇ h φ λ 2 2 - 1 2 d dt N -1 ∂ Z φ λ 2 2 = ν 1h ∆ h φ λ 2 2 + ν 3h N -1 ∇ h ∂ Z φ λ 2 2 + ν 3z N -1 ∂ z ∂ Z φ λ 2 2 + ν 3z D ∂ z N -2 ∂ z ∂ Z φ λ ∂ Z φ λ ,
that reads:

∇ h φ λ is bounded in L ∞ (0, T, L 2 (D)) if, ∀λ > 0, ∇ h φ λ (t = 0) is in L 2 (D), ∂ Z φ λ is bounded in L ∞ (0, T, L 2 (D)) if, ∀λ > 0, ∂ Z φ λ (t = 0) is in L 2 (D), and
∆ h φ λ is bounded in L 2 (0, T, L 2 (D)), ∇∂ Z φ λ is bounded in L 2 (0, T, L 2 (D)),
as the Brunt-Väisälä frequency N (z) is bounded (from above and away from zero). Let us also multiply equation (35) by ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ and integrate over D. We obtain:

D ∂ t ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ + β 0 ∂ x φ λ ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ = ν 1h D ∆ h ∆ h φ λ ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ + ν 3h D ∂ Z N -2 ∆ h ∂ Z φ λ ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ + ν 3z D ∂ Z N -2 ∂ 2 zz ∂ Z φ λ ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ , that gives: 1 2 d dt ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ 2 2 + ν 1h ∇ h ∆ h φ λ 2 2 + (ν 1h + ν 3h ) N -1 ∆ h ∂ Z φ λ 2 2 + ν 3z N -1 ∇ h ∂ 2 zZ φ λ 2 2 + ν 3h ∂ Z (N -2 ∇ h ∂ Z φ λ ) 2 2 + ν 3z ∂ Z (N -2 ∂ z ∂ Z φ λ ) 2 2 = -β 0 D ∂ x φ λ ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ + ν 3z D ∂ z N -2 ∂ 2 zZ φ λ ∆ h ∂ Z φ λ .
The r.h.s. can be bounded by

C 1 + C 2 ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ 2 2 + C 3 ∆ h ∂ Z φ λ 2
2 using the previous results, where the C i are constants, leading to:

∆ h φ λ is bounded in L ∞ (0, T, L 2 (D))
if, ∀λ > 0, ∆ h φ λ (t = 0) is in L 2 (D),

∂ Z N -2 ∂ Z φ λ is bounded in L ∞ (0, T, L 2 (D)) if, ∀λ > 0, ∂ Z N -2 ∂ Z φ λ (t = 0) is in L 2 (D), N -1 ∂ Z ∇ h φ λ is bounded in L ∞ (0, T, L 2 (D)) if, ∀λ > 0, N -1 ∂ Z ∇ h φ λ (t = 0) is in L 2 (D)
, and

∇ h ∆ h φ λ is bounded in L 2 (0, T, L 2 (D)), ∂ Z ∆ h φ λ is bounded in L 2 (0, T, L 2 (D)), ∂ Z N -2 ∇∂ Z φ λ is bounded in L 2 (0, T, L 2 (D)), ∂ 2 zZ ∇ h φ λ is bounded in L 2 (0, T, L 2 (D)).
From equality (35), thanks to the previous properties, we also have:

∂ t ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ is bounded in L 2 (0, T, H -1 (D)).
2.2.2. Convergence when λ → 0.

From the previous results, the gradient of φ λ is in L ∞ (0, T, L 2 (D)), so with Sobolev embeddings, φ λ strongly converges to a limit denoted φ in L ∞ (0, T, L p (D)), with p < 6. Let us prove that φ satisfies the traditional quasi-geostrophic equation (36).

We just proved that ∂ t ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ is bounded in L 2 (0, T, H -1 (D)). With the previous results, we also have ∇ ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ in L 2 (0, T, L 2 (D)): thanks to Aubin-Lions theorem ∆ h φ λ + ∂ Z N -2 ∂ Z φ λ strongly converges towards ∆ h φ+∂ z N -2 ∂ z φ in C 0 (0, T, L 2 (D)) (we recall that ∂ Z = ∂ z +λ∂ y ). So we can pass to the limit in the jacobian thanks to the strong convergence, the other terms do not bring any difficulty. Then, up to a subsequence, φ λ strongly converges in L ∞ (0, T, L p (D)), with p < 6 towards φ, weak solution of the quasi-geostrophic equation (36).

Conclusion

In this work we have shown that results of the reference paper [START_REF] Bourgeois | Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean[END_REF] could be extended to quasi-hydrostatic modeling, leading to the QHQG equations (see also slanted or tilted models in [START_REF] Embid | Low Froude number limiting dynamics for stably stratified flow with small or finite Rossby numbers[END_REF] and [START_REF] Julien | Generalized quasi-geostrophy for spatially anisotropic rotationally constrained flows[END_REF]) with few adjustments. To the authors' opinion, this should advocate for the use of QHQG modeling rather than QG when the nontraditional parameter λ is not small (in particular close to the equator). However, the convergence studies confirm that QG is a relevant approximation of QHQG for small λ configurations. A numerical study to quantify more precisely the effect of the traditional approximation could be performed but it is beyond the scope of this work.

Theorem 1 . 1 (

 11 Short-time existence of solutions to the QHQG model). If the initial vorticity ω 0 is in H s (B) for some s ≥ 3 with |ω 0 | s ≤ M , then there exists a time T * > 0 and a solution ω in C ([0, T * ], H s (B)) to the QHQG model (18)-(24), where T * depends only on M, B, λ and β 0 . The vorticity ω satisfies the estimate ω s,T * ≤ 2M .

Theorem 1 . 2 (

 12 Global existence of solutions to the QHQG model). If ω 0 is in H s (B) for some s ≥ 3, then given any time T > 0, there exists a solution ω in C ([0, T ], H s (B)) to the QHQG model (18)-[START_REF] Sheremet | Laboratory experiments with tilted convective plumes on a centrifuge: a finite angle between the buoyancy force and the axis of rotation[END_REF].

Theorem 1 . 3 (

 13 Existence of solution to the dissipative QHQG model). If N (z) is continuous and continuously differentiable on [0, h], if ω I and its horizontal derivatives (up to the second order) are in C γ (B), then there exists a unique solution (Φ * , ω * ) to the QHQG problem (26)-(28) in the time interval [0, T * ], where T * is inversely proportional to the β 0 parameter, with:• Φ * has all its spatial partial derivatives ∂ • Φ * (for 0 ≤ ≤ 2) in C 0 (B × [0, T * ]) and its second order spatial derivatives∂ 2 • Φ * in C γ (B) uniformly in t ∈ [0, T * ],• ω * , its horizontal derivatives (up to the second order) and ∂ t ω * are in C 0 (B × [0, T * ]) and in C γ (B) uniformly in t ∈ [0, T * ].

  which is thermal wind balance. 1.4.4. Rossby wave modes One class of simple solutions is Rossby wave modes. They satisfy a linearized PDE usually justified by an assumption of small amplitude flow,

Figure 1 .

 1 Figure 1. Values of ρ k,l,m (λ) for several values of m. The main differences can be seen with m = 0. Other values of k, l, m (in particular with m > 0) do not bring any additional information. The computations are made with a typical value of the Brunt-Väisälä frequency N = 0.01.

Theorem 2 . 1 .

 21 Let φ λ be a sequence of weak solutions of the nontraditional quasi-geostrophic equation (35) with initial data φ λ (t = 0) = φ λ 0 satisfying φ λ 0 → φ 0 in L 1 (D) with, for all λ > 0,∇φ λ 0 ∈ L 2 (D), ∆φ λ 0 ∈ L 2 (D), and ∂ Z φ λ 0 ∈ L 2 (D).

This assumption correspond to a small Froude number, of order ε.

Acknowledgments

This work was supported by ANR through contract ANR-11-MONU-005 (COMODO). The authors are grateful to the anonymous referees for their suggestions and comments that widely improved the manuscript. The first author would also like to thank Philippe Grillot for helpful discussion related to this work.