A V E N A I L L E

S M G R I F F I

The catalytic role of beta effect in barotropization processes

The vertical structure of freely evolving, continuously stratified, quasi-geostrophic flow is investigated. We predict the final state organization, and in particular its vertical structure, using statistical mechanics and these predictions are tested against numerical simulations. The key role played by conservation laws in each layer, including the finegrained enstrophy, is discussed. In general, the conservation laws, and in particular that enstrophy is conserved layer-wise, prevent complete barotropization, i.e., the tendency to reach the gravest vertical mode. The peculiar role of the β-effect, i.e. of the existence of planetary vorticity gradients, is discussed. In particular, it is shown that increasing β increases the tendency toward barotropization through turbulent stirring. The effectiveness of barotropization may be partly parameterized using the Rhines scale 2πE 1/4 0 /β 1/2 . As this parameter decreases (beta increases) then barotropization can progress further, because the beta term provides enstrophy to each layer. However, if the beta effect is too large then the statistical mechanical predictions fail and wave dynamics prevent complete barotropization.

Introduction

The vertical structure of eddying flow in the oceanic mesoscale is a fundamental problem in geophysical fluid dynamics, one that has has been reinvigorated by the need to interpret altimetric observation of surface velocity fields [START_REF] Scott | Direct Evidence of an Oceanic Inverse Kinetic Energy Cascade from Satellite Altimetry[END_REF][START_REF] Lapeyre | What Vertical Mode Does the Altimeter Reflect? On the Decomposition in Baroclinic Modes and on a Surface-Trapped Mode[END_REF]). Energy cascades and self-organization processes are fairly well understood and characterized in two-dimensional turbulence, but how these results generalize to rotating, continuously stratified flows is still an open question. We ask in this paper what is the final organization of a freely evolving inviscid quasi-geostrophic flow with continuous stratification. In particular, we examine what are the precise conditions for the oft-cited barotropization process to occur. Barotropization refers to the tendency of a quasi-geostrophic flow to reach the gravest vertical mode [START_REF] Charney | Geostrophic Turbulence[END_REF][START_REF] Rhines | The dynamics of unsteady currents[END_REF]. Because the gravest vertical mode is the barotropic one, i.e. a depth independent flow, "barotropization" means a tendency toward the formation of depth independent flows. We study in particular the key role played by the beta effect (the existence of planetary vorticity gradients) in such barotropization processes. It has been previously noticed that the presence of these large scale planetary vorticity gradients favors barotropic process, see e.g. [START_REF] Smith | The Scales and Equilibration of Midocean Eddies: Freely Evolving Flow[END_REF], but this phenomenon remains to our knowledge unexplored and unexplained.

Continuously stratified quasi-geostrophic flows take place in three dimensions, but their dynamics are quasi-two-dimensional because the non-divergent advecting velocity field has only horizontal components. The dynamics admits similar conservation laws as the two-dimensional Euler equations, including among others the total energy and the enstrophy of each layer. These conservation laws imply an inverse energy cascade toward large scales, and weak energy dissipation in presence of weak viscosity [START_REF] Kraichnan | Two-dimensional turbulence[END_REF]. This situation contrasts with three dimensional turbulence where the energy cascades forward toward small scales, giving rise to finite energy dissipation no matter how weak the viscosity.

Just as for two-dimensional Euler flows, inviscid, freely evolving continuously stratified quasi-geostrophic flows are disordered and an unstable initial condition rapidly evolves into a strongly disordered fine-grained flow with filaments being stretched, folded, and consequently cascading toward smaller and smaller scales. However, the inverse energy cascade leads to the formation of robust, large scale coherent structures filling the domain in which the flow takes place, contrasting strongly with the analogous three dimensional flow with its energy transfer to small scales.

The aim of this paper is to provide a prediction for the vertical structure of the large scale flow organization resulting from inviscid, freely evolving continuously stratified quasi-geostrophic dynamics using a statistical mechanics approach. The goal of the statistical mechanics approach is to drastically reduce the complexity of the problem of determining the large scale flow organization to the study of a few key parameters given by the dynamical invariants, as for instance the total energy and the fine-grained enstrophy in each layer. Such a theory has been independently proposed by [START_REF] Robert | Etats d'equilibre statistique pour l'ecoulement bidimensionnel d'un fluide parfait[END_REF]; [START_REF] Robert | Statistical equilibrium states for two-dimensional flows[END_REF] and [START_REF] Miller | Statistical mechanics of euler equations in two dimensions[END_REF], and henceforth will be referred to as MRS theory. It is an equilibrium theory, developed in the case without forcing and dissipation. Our contribution will be to apply this theory to the explicit computation of a class of equilibrium states of the stratified quasi-geostrophic equations.

A crucial ingredient of the theory comes from the constraints given by dynamical invariants. In the present paper, we will focus on the role played by the energy and the enstrophy constraints. In the framework of the MRS theory, it is justified in a low energy limit, and it is only a first step before more comprehensive studies including the effects of additional invariants, see e.g [START_REF] Majda | Nonlinear Dynamics and Statistical Theories for Basic Geophysical Flows[END_REF]; [START_REF] Corvellec | A complete theory of low-energy phase diagrams for two-dimensional turbulence equilibria[END_REF].

The phenomenology of two-dimensional and geostrophic turbulence is often explained by considering the energy and enstrophy constraints only, see e.g [START_REF] Kraichnan | Two-dimensional turbulence[END_REF]; [START_REF] Salmon | Lectures on Geophysical Fluid Dynamics[END_REF]. Actually, earlier statistical mechanics approaches were developed for Galerkin-truncated Euler or quasi-geostrophic dynamics, which allowed for a simpler theoretical treatment than in the continuous dynamics case, see e.g. [START_REF] Kraichnan | Inertial Ranges in Two-Dimensional Turbulence[END_REF]; [START_REF] Salmon | The equilibrium statistical mechanics of simple quasi-geostrophic models[END_REF]. In the case of Galerkin-truncated models, only the energy and the enstrophy are conserved quantities, so the statistical mechanics of these systems were called energy-enstrophy theory.

A first reason to consider MRS theory rather than earlier energy-enstrophy statistical mechanics approach is that real flows are continuous, not Galerkin-truncated. Because real flows are also forced-dissipated, our first working hypothesis is that the computation of the MRS equilibrium states of the continuous dynamics is a necessary step before taking into account the effect of weak forcing and dissipation. This is for instance the approach followed by [START_REF] Bouchet | Random Changes of Flow Topology in Two-Dimensional and Geophysical Turbulence[END_REF], who found that in the presence of weak forcing and dissipation, the dynamics were close to MRS equilibrium states which were different from energy-enstrophy equilibrium states.

One could argue that numerical models often at best conserve energy and enstrophy, so there would be at first sight no reason to consider higher order invariants to predict the simulated self-organization. However, the truncated dynamics in numerical models may lead to the formation of quasi-stationary states having a life time that tends to infinity as the resolution is refined. These quasi-stationary states are distinct from the energyenstrophy equilibria, and thus require more invariants to describe using a statistical theory. The formation of such quasi-stationary states is a characteristic of long-range interacting systems (e.g., [START_REF] Campa | Statistical mechanics and dynamics of solvable models with long-range interactions[END_REF] and references therein). Given these results, our second working hypothesis is that the quasi-stationary states found in the truncated dynamics of numerical simulations may be close to a state predicted by MRS theory.

There are only a few studies of statistical equilibria of the continuously stratified quasigeostrophic dynamics. [START_REF] Merryfield | Effects of stratification on quasi-geostrophic inviscid equilibria[END_REF] obtained and discussed a class of solutions in the framework of energy-enstrophy statistical theory above topography, for which there is linear relation between potential vorticity q and streamfunction ψ. The solutions obtained were critical points of the theory, but the actual entropy maxima were not selected. Here we consider the particular class of MRS equilibrium states characterized by a linear qψ relation, so we describe the same class of equilibrium states, except that we do not take into account bottom topography. The effect of bottom topography will be addressed in a future work. The novelty of our work is to relate each critical point to the energy and fine-grained enstrophy constraints and to select the actual entropy maxima in different limit cases. [START_REF] Schecter | Maximum entropy theory and the rapid relaxation of three-dimensional quasi-geostrophic turbulence[END_REF] provided a direct comparison between numerical simulations and numerical computation of MRS statistical equilibria, focusing on the particular case of initial conditions given by a random field statistically homogeneous in space, as in [START_REF] Mcwilliams | Anisotropy and Coherent Vortex Structures in Planetary Turbulence[END_REF]. This corresponds to a particular case for which the enstrophy profile is only weakly varying on the vertical. We will outline in this paper the important consequence of the conservation of enstrophy at each depth. This important role played by the enstrophy constraint allows us to discuss how the beta effect can favor barotropization. For a given initial streamfunction field, the beta effect modifies the initial distribution of potential vorticity, and therefore modifies the vertical profile of fine-grained enstrophy. More precisely, increasing the beta effect tends to increase the depth independent part of the fine-grained enstrophy profile. We will show that statistical equilibria associated with depth independent fine-grained enstrophy profiles are barotropic, which means that increasing the beta favors barotropization, according to statistical mechanics predictions.

The paper is organized as follows: the continuously stratified quasi-geostrophic equations are introduced on section 2. The MRS equilibrium statistical mechanics theory of such systems is presented in section 3. Here it is shown that the class of MRS statistical equilibria characterized by a linear relation between potential vorticity and streamfunction are solutions of a minimum enstrophy problem, which is solved in various limit cases, first without beta effect and then with beta effect. The theory predicts an increase of barotropization with increasing value of beta, which is tested against numerical simulations in section 4. A summary of the main results and a discussion of their application to the ocean is given section 5.

Continuously stratified quasi-geostrophic model

Equations of motion

We consider an initial value problem for a freely evolving, unforced, undissipated, quasigeostrophic flow with continuous stratification N (z) (see e.g., [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF], section 5.4):

∂ t q + J(ψ, q) = 0 ,
(2.1)

q = ∆ψ + ∂ ∂z f 2 0 N 2 ∂ ∂z ψ + βy , (2.2)
with J(ψ, q) = ∂ x ψ∂ y q -∂ y ψ∂ x q , and f = f 0 + βy is the Coriolis parameter with the β-plane approximation. Neglecting buoyancy and topography variations, the boundary condition at the bottom (z = -H, where H is a constant) is given by

∂ z ψ| z=-H = 0. (2.
3)

The boundary condition at the surface (defined as z = 0, using the rigid lid approximation), is given by the advection of buoyancy

∂ t b s + J(ψ, b s ) = 0, f 2 0 N 2 ∂ z ψ| z=0 = b s .
(2.4)

The reason we focus on quasi-geostrophic dynamics is that for these flows, there is no direct cascade of energy and no dissipative anomaly, see e.g. [START_REF] Salmon | Lectures on Geophysical Fluid Dynamics[END_REF], which makes possible a straightforward generalization of statistical mechanics theories developed initially in the framework of two-dimensional Euler equations. As far as the ocean is concerned, quasi-geostrophic dynamics is relevant to describe "mesoscale turbulence", but does not account for "sub-mesoscale turbulence", which involve ageostrophic effects, with part of the energy that may cascade forward, see e.g. [START_REF] Thomas | Sub-mesoscale processes and dynamics[END_REF][START_REF] Ferrari | Ocean Circulation Kinetic Energy: Reservoirs, Sources, and Sinks[END_REF] and references therein.

Modal decomposition

Laplacian and baroclinic eigenmodes

We consider the case of a doubly periodic domain

D = [-π π] × [-π π].
The Fourier modes e k,l (x, y) are a complete orthonormal basis of the 2D Laplacian ∆e k,l = -K 2 e k,l where (k, l) are the wavenumbers in each directions and K 2 = k 2 + l 2 . In the vertical, the complete orthonormal basis of barotropic-baroclinic modes is defined by the solutions F m (z) of the Sturm-Liouville eigenvalue problem

∂ ∂z f 2 0 N 2 ∂ ∂z F m = -λ 2 m F m with F ′ (0) = F ′ (-H) = 0, (2.5)
where m ≥ 0 is an integer, λ 2 m is the eigenvalue associated with mode m, which defines the m th -baroclinic Rossby radius of deformation λ -1 m . The barotropic mode is the (depth independent) mode F 0 , associated with λ 0 = 0.

Surface quasi-geostrophic (SQG) modes

In the following, it will be convenient to distinguish the interior dynamics from the surface dynamics. For a given potential vorticity field q(x, y, z, t) and surface buoyancy field b s (x, y, t) , the streamfunction can be written ψ = ψ int + ψ surf , where ψ int is the solution of

∆ψ + ∂ ∂z f 2 0 N 2 ∂ ∂z ψ = q -βy with ∂ z ψ| z=0 = ∂ z ψ| z=-H = 0, (2.6)
and where ψ surf is the solution of

∆ψ + ∂ ∂z f 2 0 N 2 ∂ ∂z ψ = 0, with ∂ z ψ| z=0 = N (0) f 0 2 b s , ∂ z ψ| z=-H = 0 . (2.7)
Each projection of ψ surf on a (horizontal) Laplacian eigenmode (l, k) defines a surface quasi-geostrophic eigenmode (SQG mode hereafter). When the stratification N is constant and when f 0 /(N K) ≪ H, these SQG modes are decreasing exponential functions

ψ surf,k,l = A k,l e zN K/f0 + O h K H , (2.8)
where the coefficients A k,l are determined using the boundary condition at the surface. When the stratification is not constant, but when f 0 /N (0)K ≪ H, the SQG modes are still characterized by a decreasing exponential, with an e-folding depth h K = f 0 /KN (0). The distinction between surface and interior streamfunction is sometimes useful to describe the dynamics , see e.g. [START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF], but a shortcoming of such a decomposition is that ψ int and ψ surf are not orthogonal.

Delta-sheet approximation and SQG-like modes

The boundary condition (2.4) can be formally replaced by the condition of no buoyancy variation (∂ z ψ = 0 at z = 0), provided that surface buoyancy anomalies are interpreted as a thin sheet of potential vorticity just below the rigid lid [START_REF] Bretherton | Critical layer instability in baroclinic flows[END_REF]. For this reason, and without loss of generality, we will consider that b s = 0 in the remainder of this paper.

Let us now consider a case with no surface buoyancy, but with a surface intensified potential vorticity field defined as q = q surf (x, y)Θ(z + H 1 ) + q int (x, y, z), with H 1 q surf ≫ Hq int , H 1 ≪ H and Θ the Heaviside function. Using the linearity of (2.6), one can still decompose the streamfunction field into ψ = ψ surf + ψ int where ψ surf is the flow induced by the PV field q = q surf Θ(z + H 1 ), obtained by inverting (2.2), and ψ int the flow induced by q int . For z < -H 1 , the streamfunction ψ surf satisfies equation (2.7), just as SQG modes. This is why we call it an SQG-like mode.

Equilibrium statistical mechanics of stratified quasi-geostrophic flow

In this section we briefly introduce the MRS equilibrium statistical mechanics theory, and introduce a general method to compute the MRS equilibria. This allows for explicit computation of statistical equilibria on the f -plane in subsection 3.3, and on the β-plane in subsection 3.4.

Theory

Let us consider an initial (fine-grained) PV field q 0 (x, y, z). The dynamics (2.1) induced by this PV field stirs the field itself, which rapidly leads to its filamentation at smaller and smaller scale, with the concomitant effect of stretching and folding. There are two (related) consequences, one in phase space, one in physical space.

In physical space, the filamentation toward small scales implies that after a sufficiently long time evolution, the PV field can be described at a macroscopic level by a probability field ρ(x, y, z, σ, t)dσ to measure a given PV level between σ and σ + dσ at a given point (x, y, z). This probability distribution function (pdf hereafter) ρ is normalized locally

∀ x, y, z N [ρ] = Σ dσ ρ = 1, (3.1)
where the integral is performed over the range Σ of PV levels, which is prescribed by the initial condition. The field ρ is a macrostate because it is associated with a huge number of fine-grained PV configurations q(x, y, z) (the microscopic states), in the sense that many realizations of the (microscopic, or fine-grained) PV field q leads to the same macroscopic state ρ.

The phase space is defined by all configurations of the field q(x, y, z). The cornerstone of the theory is to assume that stirring allows the system to explore evenly each possible configuration of the phase space that satisfies the constraints of the dynamics, so that these configurations can be considered equiprobable. When this assumption fails, as in the linear regime discussed in section 4.4, the statistical theory fails.

The main idea of the theory is then to count the number of microstate configurations associated with a macrostate (the pdf field ρ), and to say that the equilibrium state is the most probable macrostate that satisfies the constraints of the dynamics, i.e. the one that is associated with the largest number of microstates satisfying the constraints. One can then show that an overwhelming number of microstates are associated with this most probable macrostate [START_REF] Michel | Large deviations for young measures and statistical mechanics of infinite dimensional dynamical systems with conservation law[END_REF]). An important physical consequence is that an observer who wants to see the actual equilibrium macrostate would neither need to perform ensemble average nor time average, but would simply have to wait a sufficiently long time.

Note that the output of the theory is the pdf field ρ, but the quantity of interest to describe the flow structure is the coarse-grained PV field defined as

q = Σ dσ σρ. (3.2)
In practice, counting the number of microstates q associated with a macrostate ρ is a difficult task, because q is a continuous field, and it requires the use of large deviation theory, see e.g. [START_REF] Bouchet | Invariant measures of the 2D Euler and Vlasov equations[END_REF]. However, it can be shown that classical counting arguments do apply here, and that the most probable state ρ maximizes the usual Boltzmann mixing entropy

S = - D dxdy 0 -H dz Σ dσ ρ ln ρ , (3.3)
see [START_REF] Bouchet | Invariant measures of the 2D Euler and Vlasov equations[END_REF] and references therein for further details.

Let us now collect the constraints that ρ must satisfy, in addition to the local normalization (3.1). An infinite number of constraints are given by conservation of the global distribution of PV levels P(σ, z), which is prescribed by the initial global distribution of PV levels P 0 (σ, z):

P(σ, z)[ρ] = D dxdy ρ = P 0 (σ, z).
(3.4)

These conservation laws include the conservation of the fine-grained enstrophy

D dxdy Σ dσ σ 2 ρ = Z 0 . (3.5)
Another constraint is given by the conservation of the total energy

E = 1 2 0 -H dz D dxdy (∇ψ) 2 + f 2 0 N 2 (∂ z ψ) 2 . (3.6)
Because an overwhelming number of microstates are close to the equilibrium state ρ, the energy of the fluctuations are negligible with respect to the energy due to the coarsegrained PV field q (Michel & Robert 1994). Integrating 3.6 by parts and replacing q by q gives the constraint

E 0 = E[ρ]
, where E 0 is the energy of the initial condition, and where

E = 1 2 0 -H dz D dxdy (βy -q) ψ = 1 2 0 -H dz D dxdy Σ dσ (βy -σ) ψρ . (3.7)
Finally, the MRS theory provides a variational problem

S (E 0 , P 0 (σ, z)) = max ρ,N [ρ]=1 {S [ρ] | E [ρ] = E 0 & P [ρ] = P 0 (σ, z)} , (3.8)
which means that the equilibrium state is the density probability field ρ rsm that satisfies the local normalization constraint (3.1), the energy constraint E [ρ] = E 0 with E given by (3.7), the incompressibility constraint P(σ, z) [ρ] = P 0 (σ, z) with P given by (3.4), and that maximizes the entropy functional S given by (3.3). We see that the parameters of this problem are the values of the constraints, which are prescribed by the initial condition q 0 (x, y, z).

3.2. MRS equilibrium states characterized by a linear qψ relation

Simplification of the variational problem

In practice, MRS equilibrium states are difficult to compute, because such computations involve the resolution of a variational problem with an infinite number of constraints. Computing the critical points of (3.8) is straightforward, see e.g. equation (A 1) of appendix A, but showing that second order variations of the entropy functional are negative is in general difficult. However, recent results have allowed considerable simplifications of the analytical computations of these equilibrium states [START_REF] Bouchet | Simpler variational problems for statistical equilibria of the 2d euler equation and other systems with long range interactions[END_REF]. The idea is first to prove equivalence (usually for a restricted range of parameters) between the complicated variational problem (3.8) and other variational problems more simple to handle, involving less constraints than the complicated problem, but characterized by the same critical points, and then to explicitly compute solutions of the simpler variational problem. The equivalence between the variational problems ensures that solution of the simple problem are also solutions of the more complicated one. We will follow this approach in the following.

Minimum enstrophy variational problem

Starting from the MRS variational (3.8), it is shown in appendix A that any MRS equilibrium state characterized by a linear qψ relation is a solution of the variational problem

Z tot cg (E 0 , Z 0 (z)) = min q Z tot cg [q] | E [q] = E 0 (3.9) with Z tot cg [q] = 0 -H dz Z cg Z 0 , Z cg = 1 2 D dxdy q 2 , (3.10)
where Z 0 (z) is the fine-grained enstrophy profile, E 0 the energy. This variational problem amounts to find the minimizer q m (a coarse grained PV field) of the coarse-grained enstrophy Z tot cg (with

Z tot cg [q m ] = Z tot cg (E 0 , Z 0 (z))
), among all the fields q satisfying the constraint E[q] = E 0 given by (3.6).

Critical points of the variational problem (3.9) are computed by introducing a Lagrange multiplier β t associated with the energy constraint, and by solving δZ tot cg + β t δE = 0, where first variations of the functionals are taken with respect to the field q. Because β t is the Lagrange parameter associated with the energy constraint, it is called inverse temperature †. Using δZ tot cg = D dxdy 0 -H dz qδq/Z 0 and δE =dxdydz ψδq, one finds that critical points of this problem are q = β t Z 0 ψ.

(3.11) † The notation β is traditionally used in thermodynamics for inverse temperature. We use here the convention βt for this inverse temperature, since β = ∂yf refers to the beta effect, i.e. to the variation of the Coriolis parameter with latitude.

The flow structure of these critical points can then be computed by solving

∆ψ + ∂ ∂z f 2 0 N 2 ∂ ∂z ψ + βy = Z 0 β t ψ with ∂ z ψ| z=0,z=-H = 0.
(3.12)

Critical points can then be classified according to the value of their inverse temperature β t . All the problem is then to find which of these critical point are actual minima of the coarse-grained enstrophy Z tot cg [q] for a given energy E[q] = E 0 . MRS equilibrium states characterized by linear qψ relations are expected either in a low energy limit or in a strong mixing limit, see [START_REF] Bouchet | Statistical mechanics of two dimensional and geophysical turbulent flows[END_REF] and references therein. Let us estimate when the low energy limit is valid. For a given distribution P 0 (σ, z) of PV levels σ, one can compute the maximum energy E max among all the energies of the PV fields characterized by the distribution P 0 (σ, z). Low energy states characterized by the same distribution P 0 (σ, z) are therefore those characterized by E ≪ E max . The strong mixing limit, which can in some case overlap the low energy limit, is obtained by assuming β t σψ ≪ 1 in the expression (A 1) of the critical points of the MRS variational problem (3.8), see [START_REF] Chavanis | Classification of self-organized vortices in twodimensional turbulence: the case of a bounded domain[END_REF]. Note that these results would remain valid in the presence of bottom topography. The only change would be the expression of the energy, and the bottom boundary condition.

Link with the phenomenological principle of enstrophy minimization and other equilibrium states

At a phenomenological level, the variational problem (3.9) can be interpreted as a generalization of the minimum enstrophy principle of [START_REF] Bretherton | Two-dimensional turbulence above topography[END_REF] to the continuous stratified flow: the equilibrium state minimizes the vertical integral of the coarse-grained enstrophy normalized by the fine-grained enstrophy in each layer, with a global constraint given by energy conservation. This normalization is necessary to give a measure of the degree of mixing comparable from one depth z to another.

It has been shown by [START_REF] Carnevale | Nonlinear stability and statistical mechanics of flow over topography[END_REF] that solutions of the minimum enstrophy problem of [START_REF] Bretherton | Two-dimensional turbulence above topography[END_REF] for one layer quasi-geostrophic flows, are energy-enstrophy equilibrium states of the Galerkin-truncated dynamics of these models, in the limit when the high-wave number cut-off goes to infinity. These energy-enstrophy equilibrium states were first introduced in the framework of 2D truncated Euler flows by [START_REF] Kraichnan | Inertial Ranges in Two-Dimensional Turbulence[END_REF]. [START_REF] Salmon | The equilibrium statistical mechanics of simple quasi-geostrophic models[END_REF] generalized the theory to one or few layers truncated quasi-geostrophic flows. [START_REF] Frederiksen | Statistical Dynamics of Two-Dimensional Inviscid Flow on a Sphere[END_REF] applied the theory for a barotropic flow on a sphere, and Frederiksen (1991b) generalized these results for two-layers quasi-geostrophic flow on a sphere. In this latter case, interesting stability results were obtained by Frederiksen (1991a).

The case of continuously stratified flow above topography was addressed by Merryfield (1998) for truncated dynamics. [START_REF] Merryfield | Effects of stratification on quasi-geostrophic inviscid equilibria[END_REF] computed critical points of the statistical mechanics variational problem, found q = µ t (z)ψ and discussed the vertical structures of these flows depending on µ t . We showed in the previous subsection that the parameter µ t is related to the fine-grained enstrophy profile through µ t = β t Z 0 (z). In addition, the critical states described by [START_REF] Merryfield | Effects of stratification on quasi-geostrophic inviscid equilibria[END_REF] could either be saddle, minima or maxima of the entropy functional. We show in the next section that we can select which of them are actual entropy maxima (or equivalently coarse-grained enstrophy minima), in different limit cases.

Transients, mean and sharp equilibrium states

In the framework of energy-enstrophy theory, a distinction is often made between "transients" and "mean" energy, see e.g. [START_REF] Salmon | The equilibrium statistical mechanics of simple quasi-geostrophic models[END_REF]; Frederiksen (1991b); [START_REF] Majda | Nonlinear Dynamics and Statistical Theories for Basic Geophysical Flows[END_REF], among others. Statistical equilibria are called "sharp" when the "transients" vanish. In earlier works, the "mean" quantities were often computed by summing over all the microstates distributed according to the canonical measure, which were implicitly assumed to be equivalent to the microcanonical measure. In the absence of topography and beta plane, the "mean" states defined in this way are always zero in a doubly periodic domain, whatever the resolution of the truncated model, and all the energy is in the "transients".

By contrast, the equilibrium states computed in the MRS framework in the microcanonical ensemble are always sharp: the energy of the fluctuations necessarily vanish, which can be shown with large deviation theory, see [START_REF] Bouchet | Invariant measures of the 2D Euler and Vlasov equations[END_REF] and references therein. This apparent paradox is solved when one realizes that microcanonical and canonical ensembles are, in general, not equivalent for long-ranged interacting systems, see e.g. [START_REF] Ellis | Large deviation principles and complete equivalence and nonequivalence results for pure and mixed ensembles[END_REF] for general considerations and Venaille & Bouchet (2011b) for application to two-dimensional and geophysical flows.

In the case of a doubly periodic domain without bottom topography and beta effect, the MRS equilibrium states (or the energy-enstrophy equilibria) are degenerate, because of the symmetries of the domain, see [START_REF] Bouchet | Statistical mechanics of two dimensional and geophysical turbulent flows[END_REF] for more details and for a discussion on how to remove the degeneracy, see also [START_REF] Herbert | Phase transitions and marginal ensemble equivalence for freely evolving flows on a rotating sphere[END_REF] for the case of a spherical domain. In practice, the dynamics breaks the system symmetry by selecting one of the degenerate equilibrium states. Of course, if one considers the average potential vorticity field over all the degenerate equilibria, then one recovers zero, and this is what happens when computing the "mean" state in the canonical ensemble. We stress that for either quasi-geostrophic flows or truncated quasi-geostrophic flows in the limit of infinite resolution, the freely evolving dynamics itself can not jump from one state to another once the symmetry is broken: in the microcanonical ensemble, there are no "transients". We also stress that the physically relevant statistical ensemble for an isolated system such as a freely evolving flow is the microcanonical ensemble.

Enstrophy minima on a f -plane

In the previous subsection, we have shown that MRS equilibrium states characterized by a linear qψ relationship are solution of a minimum coarse-grained enstrophy variational problem (3.9). Such MRS equilibrium states will be referred to as "coarse-grained enstrophy minima" in the following. We found in the previous section that the affine relation of these coarse-grained enstrophy minima are of the form q = β t Z 0 ψ, with a coefficient proportional to the fine-grained enstrophy Z 0 , and are solutions of (3.12). The critical points of the variational problem (3.9) are any states q satisfying these conditions. The aim of this subsection is to find which of these critical point q are the actual coarse-grained enstrophy minima, solutions of variational problem (3.9), when there is no β-effect, but for arbitrary stratification N .

Injecting (3.11) in (3.10), and using the expression (3.7) with the constraint E 0 = E[ρ], one finds that the coarse-grained enstrophy of each critical point q is given by

Z tot cg = -β t E 0 . (3.13)
Remarkably, the coarse grained enstrophy of a critical point q depends only on the inverse temperature β t . We conclude that coarse-grained enstrophy minimum states are the solutions of (3.12) associated with the largest value β t . Projecting (3.12) on the Laplacian eigenmode e l,k (x, y) gives with

∂ ∂z f 2 0 N 2 ∂ ∂z ψ k,l = β t Z 0 + K 2 ψ k,l , (3.14) 0 0.5 1 -1 -0.8 -0.6 -0.4 -0.2 0 z Z 0 0 0.5 1 -1 -0.8 -0.6 -0.4 -0.2 0 z ψ(z)/ψ(0) -H 1 / Z surf h A A B C B C
∂ z ψ k,l | z=0,-H = 0, K 2 = k 2 + l 2 , ψ = k,l ψ k,l e kl .
We see that each critical point is characterized by a given wavenumber modulus K. Its vertical structure and the corresponding value of β t must be computed numerically in the case of arbitrary profiles Z 0 (z). Let us consider the example shown in Fig. 1, for a two-step fine-grained enstrophy profile

Z 0 = Z surf Θ (z + H 1 ) + Z int Θ (-z -H 1 ) , H 1 ≪ H, (3.15)
where Θ is the Heaviside function, and for Z int /Z surf varying between 0 and 1. We find that the minimum coarse-grained enstrophy states are always characterized by the gravest horizontal mode on the horizontal (K = 1). As for the vertical structure, we observe Fig. 1 a tendency toward more barotropic flows when the ratio Z int /Z surf tends to one.

It is instructive to consider the limit cases Z surf = Z int and Z int = 0, corresponding respectively to case A and C of Fig. 1, for which minimum enstrophy can be explicitly solved. When the enstrophy profile is depth independent (Z 0 (z) = Z int , case C of figure 1), solutions of (3.14) are given by the modes F m (z)e K (x, y) defined in subsection (2.2), and are associated with inverse temperatures β t =λ 2 m + K 2 /Z int . We see that the maximum value of β t is reached for the gravest horizontal mode (K = 1) and the gravest vertical mode (λ 0 = 0). It means that the coarse-grained enstrophy minimum state is barotropic.

When Z int = 0 (case A of Fig. 1), it is straightforward to show that critical points, the solutions of (3.14), are SQG-like modes (see subsection 2.3) of e-folding depth h = f 0 /N (0)K, associated with β t = -N (0)K/(H 1 Z surf f 0 ) with corrections of order H 1 /H and of order h/H. The coarse-grained enstrophy minimum is therefore the SQG-like mode associated with K = 1 and with the largest e-folding depth h = Lf /2πN , where L is the domain length scale.

These examples show the importance of the conservation of fine-grained enstrophy to the vertical structure of the equilibrium state. The main result is that statistical mechanics predicts a tendency for the flow to reach the gravest Laplacian mode on the horizontal. The vertical structure associated with this state is fully prescribed by solving (3.14) with K = 1. Because the barotropic component of such flows are larger than solutions of (3.14) with K > 1, we can say that the inverse cascade on the horizontal is associated with a tendency to reach the gravest vertical mode compatible with the vertical fine-grained enstrophy profile Z 0 . This means a tendency toward barotropization, although in general, the fact that the profile Z 0 is non constant prevents complete barotropization.

Note that the previous results are presented in the case of a doubly periodic domain, but it would be straightforward to generalize them to channel geometries, or to any other domains with boundaries. The only change would be the spatial structure of the Laplacian eigenmodes. For these domain geometries, one would also need to discuss the effect of the circulation values at the boundaries.

3.4. Including the β-effect .

We now discuss how the results from the previous section generalize in the presence of beta. For a given initial condition ψ 0 (x, y, z), increasing β increases the contribution of the (depth independent) available potential enstrophy defined as (3.16) to the total fine-grained enstrophy profile Z 0 (z) = D dxdy q 2 0 , where q 0 is the initial PV field that can be computed by injecting ψ 0 in (2.2). For sufficiently large values of β, the PV field is dominated by the beta effect (q 0 ≈ βy), Z 0 therefore tends to Z p and becomes depth independent. Because statistical equilibria computed in the previous subsection were fully barotropic when the fine-grained enstrophy Z 0 was depth-independent, we expect a tendency toward barotropization by increasing β.

Z p = β 2 D dxdy y 2 ,
Consider for instance the case where ψ 0 is an SQG-like mode such that q 0 = βy for z < -H 1 , with H 1 ≪ H and associated with a surface fine-grained enstrophy Z 0 = Z surf for z > -H 1 . The initial profile of fine-grained enstrophy Z 0 is thus given by (3.15), with the interior enstrophy Z int = Z p prescribed by the available potential enstrophy (3.16). The question is then to determine whether the results obtained previously without beta effect still hold in presence of it.

In oder to apply properly MRS theory in the presence of beta, one must consider southern and northern boundaries. Indeed, MRS theory relies on the conservation of the global fine-grained PV distribution. Strictly speaking, this distribution is not conserved in the doubly periodic configuration on a beta plane, since any fluid particle that travels 2π in the y direction gains a value -2πβ, and one can not know a priori how many times a fluid particle will turn over the domain in the y direction before the flow reaches equilibrium. In practice, the theory can still be used in the doubly periodic case to make qualitative predictions, as explained in the next section.

In the channel configuration, the streamfunction is a constant at y = ±π, and we can compute the statistical equilibria in two limit cases. First, when β = 0, we recover the configuration with surface intensification of the fine grained enstrophy profile Z 0 , for which we found that statistical equilibria were the SQG-like modes associated with the gravest Laplacian horizontal modes, see subsection 3.3. Second, with β sufficiently large, such that Z p ≫ Z surf , the statistical equilibria can be computed by considering the case of a constant fine-grained enstrophy profile Z 0 . In the general case, the values of the streamfunction at the northern and southern boundaries are different, and are determined using constraints given by mass conservation and circulation conservation along each boundary, see e.g. [START_REF] Pedlosky | Geophysical fluid dynamics[END_REF]). Here we consider the simple case ψ(x, y = ±π, z) = 0, with zero total circulation ( D dxdyq = 0) at each depth. It is necessary to take into account the conservation of linear momentum (associated with translational invariance of the problem in the x direction), which provides an additional constraint L[q] = L 0 , with

L[q] = 0 -H dz D dxdy qy.
(3.17)

Critical points of the variational problem (3.9) with this additional constraint satisfy q = (β t ψµy) Z 0 , (3.18)

where µ and β t are respectively the Lagrange multipliers associated with the linear momentum conservation and the energy conservation. Equilibrium states therefore satisfy

Z 0 (β t ψ -µy) = ∆ψ + ∂ ∂z f 2 0 N 2 ∂ ∂z ψ + βy . (3.19)
Taking µ = -β/Z 0 , we recover the f -plane case of the previous subsection: solutions of Eq. (3.19) are barotropic-baroclinic modes on the vertical with eigenvalue -λ 2 m , and are Laplacian eigenmodes in the channel (with eigenvalue denoted by -K 2 ), with

β t Z 0 = -K 2 + λ 2
m . These states characterized by β t = -K 2 + λ 2 m /Z 0 , µ = -β/Z 0 are drifting towards the west at speed of Rossby waves:

U drif t = -β/ λ 2 m + K 2 .
The enstrophy minimizer among these states is the one associated with the smallest value of β t , namely the barotropic mode associated with the gravest Laplacian eigenmode consistent with the circulation and momentum constraints. When µ = β/Z 0 , Eq. (3.19) can be inverted, and ψ is barotropic since (β -µZ 0 )y is depth-independent. The computation of the minimum enstrophy state is therefore equivalent to computing the minimum enstrophy state of a barotropic flow in a channel. A complete discussion of minimum enstrophy states in a barotropic channel are presented in [START_REF] Corvellec | A complete theory of low-energy phase diagrams for two-dimensional turbulence equilibria[END_REF], see also [START_REF] Majda | Nonlinear Dynamics and Statistical Theories for Basic Geophysical Flows[END_REF], with detailed computations including the effect of bottom topography and of non-zero circulations at the boundaries. The important point here is that in the large beta limit, the statistical equilibrium states are barotropic.

Numerical simulations

Experimental set-up

We consider in this section the final state organization of an initial SQG-like mode (defined in subsection 2.3), varying the values of β in a doubly periodic square domain of size L = 2π. More precisely, we consider a vertical discretization and an horizontal Galerkin truncation of the dynamics (2.1), for an initial potential vorticity field q 0 = q 0,surf (x, y)Θ (z + H 1 ) + βy, such that q 0 = βy in the interior (-H < z < -H 1 ) and q 0 ≈ q 0surf in a surface layer z > -H 1 . The surface PV q surf (x, y) is a random field with random phases in spectral space, and a Gaussian power spectrum peaked at wavenumber K 0 = 5, with variance δK 0 = 2, and normalized such that the total energy is equal to one (E 0 = 1). This corresponds to the case discussed subsection (3.4), with a vertical profile of fine-grained enstrophy given by (3.15), where the interior enstrophy is given by the available potential enstrophy (3.16): for z < -H 1 , Z 0 = 4π 4 β 2 /3.

We perform simulations of the dynamics by considering a vertical discretization with 10 layers of equal depth, horizontal discretization of 512 2 , H = 1, and F = (Lf 0 /HN ) 2 = 1, using a pseudo-spectral quasi-geostrophic model for the horizontal dynamics without small scale dissipation [START_REF] Smith | The Scales and Equilibration of Midocean Eddies: Freely Evolving Flow[END_REF]. We choose H 1 = H/10 for the initial condition, so that there is non zero enstrophy only in the upper layer in the absence of a beta effect. We also perform experiments in the presence of small scale dissipation (Laplacian viscosity or hyperviscosity), and found no qualitative differences as far as the large scale horizontal structure and the vertical structure of the flow were concerned. The only difference in presence of small scale dissipation is the decay of the enstrophy in each layer, because small scale dissipation smooths out PV filaments. Because of the energy cascade, the energy remains constant at lowest order even in presence of small scale dissipation.

Time integration of the freely evolving dynamics proceeds for about 250 eddy turnover times. Typically, i) the unstable initial condition leads to a strong turbulent stirring and the flow self-organizes rapidly into a few vortices, which takes a few eddy turnover times ii) same sign vortices eventually merge together on a longer time scale (a few dozen of eddy turnover times), and iii) the remaining dipole evolves very slowly on a slow time scale (more than hundreds of eddy turnover times). In absence of β-effect, a good indicator of the convergence is given by the convergence of the qψ relation to a single functional relation. As we shall see in the following, this phenomenology is complicated by the presence of beta.

We saw in the previous section that a proper application of the MRS theory in presence of a beta plane would require considering northern or southern boundaries, since the finegrained PV distribution is not conserved in the doubly-periodic domain case. However, we still can use the theory to interpret qualitatively the simulations in the doubly periodic case. For small values of β, planetary vorticity gradients are negligible with respect to relative vorticity and stretching terms in the PV expression, in which case one can check that the fine grained PV distribution remains nearly the same during the flow evolution, and MRS can be applied. For large values of β, the fluid particles can not escape in the meridional direction, they are confined in zonal bands having a width of the order of the Rhines scale † L R = 2π ∼ E 1/4 0 /β 1/2 , and thus the fine-grained PV distribution remains close to the initial one, and MRS theory can be applied.

SQG-like dynamics (β = 0)

The case with β = 0 is presented in Fig. 2. The initial PV and streamfunction fields have horizontal structures around wave number K 0 = 5, which is associated with an SQG-like mode of typical e-folding depth f 0 /(N K 0 ) ≈ 0.1.

An inverse cascade in the horizontal leads to flow structures with an horizontal wavenumber K decreasing with time, associated with a tendency toward barotropization: the e-folding depth of the SQG-like mode increases as f 0 /2N K. According to statistical mechanics arguments of the previous section, the flow should tend to the "gravest" horizontal SQG-like mode (K = 1), which is also the "gravest" vertical SQG mode, i.e. the one associated with the largest e-folding depth. The concomitant horizontal inverse cascade (most of the kinetic energy is in the gravest horizontal mode K = 1 at the end of the simulation) and the increase of the e-folding depth are observed on Fig. 2, showing good qualitative agreement between statistical mechanics and numerical simulations.

There are however quantitative differences between the observed final state organization and the prediction of MRS equilibria characterized by a linear qψ relation. As seen on Fig. 3, the qψ relation is well defined, but it is not linear. Taking into account higher order invariants in the MRS theory would allow to better describe this sinh-like relation, see e.g. [START_REF] Bouchet | Random Changes of Flow Topology in Two-Dimensional and Geophysical Turbulence[END_REF], but it may not be sufficient to describe the "unmixed" cores of the remaining vortex in the final dipole [START_REF] Schecter | Maximum entropy theory and the rapid relaxation of three-dimensional quasi-geostrophic turbulence[END_REF].

The presence of these "unmixed" PV blobs implies the existence of structures smaller than wavenumber K = 1. Because each projection of the PV field on a different wavenum- † Here E0 = 1 is the initial energy, so that E 1/2 0 is a good metric for typical velocities of the flow (assuming it becomes barotropic). ber is associated with an SQG-like mode associated with a different e-folding depth, the vertical structure of the flow is actually a linear combination of SQG-like modes associated with different wavenumbers. A consequence is that the effective e-folding depth of the kinetic energy is f 0 /(2N K eff ). We estimated the coefficient K eff ≈ 3.5 by considering the linear q-ψ relation passing through the extremal points of the observed q-ψ relation of Fig. 3. Let us conclude on this case by discussing the relaxation towards equilibrium: we observed that the qψ relation did not change at all between 150 and 250 eddy turnover times, proving that the system reached a stationary state. In addition, we see in Fig. 8-a that the kinetic energy of each vertical mode (and therefore of each layer) reaches a plateau after a few eddy turnover times. As expected from the MRS theory, the statistical equilibrium is sharp, there is no transient.

Effect of beta (β ∼ 1)

We now consider the free evolution of the same initial condition ψ 0 (x, y, z) as in the previous subsection, but on a beta plane. Because β = 0, the initial condition for potential vorticity is different than in the previous section: the total energy and initial velocity fields are the same as before, but the vertical profile of fine grained enstrophy is different, since there is now available potential enstrophy in the interior. As explained in subsection 3.3, this contribution from interior fine-grained enstrophy shall increase the relative barotropic component of the statistical equilibrium state. This is what is actually observed in the final state organization of Fig. 4. We conclude that in this regime, the beta effect is a catalyst † of barotropization, as predicted by statistical mechanics.

Additionally, the qψ relation presented in Fig. 5 presents a very good agreement with the prediction of our computation of barotropic MRS equilibria carried in subsection 3.4: we observe a linear relation between (qβy) and ψ with a slope β t Z 0 = -1. This means that the dipole structure observed in the streamfunction field is the actual gravest horizontal Laplacian eigenmode K = 1, and that this eigenmode is drifting westward at a speed β/K 2 . In addition, we see in Fig. 8 that after a few eddy turnover times, almost all the kinetic energy is in the barotropic mode, and that the energy levels in each vertical mode reach a plateau: there are no transients.

We conjecture that the good agreement with MRS theory occurs since the initial distribution of PV levels and the initial energy place the flow in the low energy regime considered in the Appendix. In addition, for intermediate values of β, the Rhines scale

L R = 2πE 1/4 O /β 1/2
is of the order of the domain scale. The beta effect limit the meridional displacements of fluid particles at this scale; this is why the computations performed Appendix 2 in a channel geometry are qualitatively relevant to this case.

Effect of strong Rossby waves (β ≫ 1)

For large values of beta there is still an early stage, lasting a few eddy turnover times, of turbulent stirring of the unstable initial condition. But this initial stirring is limited to scales smaller than the Rhines scale L R = 2πE 1/4 0 /β 1/2 , which eventually leads to a zonal PV field perturbed by Rossby waves around this scale, see Fig. 6. The dynamics is then dominated by Rossby waves, just as in a single layer model [START_REF] Rhines | The dynamics of unsteady currents[END_REF][START_REF] Rhines | Waves and turbulence on a β-plane[END_REF].

Because the perturbation is confined in the upper layer, irreversible turbulent stirring occurs in the upper layer only, but not at the bottom. We did check that the global distribution of coarse-grained PV levels did not change in the lower layer through time, while this distribution was changed in the upper layer. In this case, the fundamental assumption that the system explores evenly the available phase space through turbulent stirring is not valid. One might reasonably argue that the statistical mechanic theory is incomplete, for it cannot identify the phase space that can be explored. In any case, we conclude that large values of beta prevent the convergence toward this equilibrium, leaving the question of the convergence itself to future work.

Turbulent stirring vs wave dynamics

We summarize the effect of planetary vorticity gradients on Fig. 7 by considering the ratio E baro /E 0 of barotropic energy to total energy, after 250 eddy turnover times, for varying values of β. The plain red line represents statistical mechanics predictions. In order to allow comparison with the numerical experiments, for which the horizontal PV field did present structures smaller than the gravest mode K = 1 for low values of β (see subsection 4.2), we computed the vertical structure with equation (3.14), using K eff = 3.5, and using the enstrophy profile Z 0 (β) of the initial condition ψ 0 . Strictly speaking, this effective horizontal wavenumber K eff should vary with beta between K eff = 3.5 for β = 0 and K eff = 1 for β > 1, but taking into account this variations would not change much the shape of the red plain curve. . Relation between qβy and streamfunction in the upper layer, at t = 100 eddy turnover times, for β = 2.1. The fact that this relation is close to be linear shows that the small energy limit considered in this paper is relevant to this case. The coefficient of the slope is approximatively one, meaning that the energy is condensed in the gravest horizontal Laplacian eigenmode, which is propagating westward at the speed of barotropic waves. This is consistent with the computation of minimum enstrophy states on a beta plane, see subsection 3.4.

The critical value of β between the turbulent stirring regime for which statistical mechanics predictions are useful and the wave regime can be estimated by considering the case in which the Rhines scale

L Rh ∼ 2π E 1/2 0 /β 1/2
is of the order of the domain scale L = 2π. Here the total energy is E 0 = 1, which renders the critical value β ∼ 1 in our simulations.

Because the ratio E baro /E 0 gives incomplete information on the flow structure, and because dynamical information is useful to understand where statistical mechanics fails to predict the vertical structure, we represent in Fig. 8 the temporal evolution of the energies E m of the baroclinic modes † defined in (2.5).

For small values of beta (see the first panel of Fig. 8), the dynamics reaches the equilibrium after stirring during few eddy turnover times, and the contribution of each mode E m simply reflects the SQG-like structure of the final state. For values of beta of order one (see the second panel of Fig. 8), the same stirring mechanism leads to a barotropic flow after a few eddy turnover times. For values of β larger than one, there is still an initial stirring regime, but the duration of this regime tends to diminish with larger values of beta. We observe a slow transfer of energy from first baroclinic to the barotropic mode in the case β = 10 (see the third panel of Fig. 8), but it is unclear if it will eventually lead to a condensation of the energy in the barotropic mode. In the case β = 100 (see the fourth panel of Fig. 8), there is in average almost no energy transfers between the vertical modes, after the initial stirring regime. Again, it is unclear if interactions between waves would eventually condense their energy into the barotropic mode in a longer numerical run. The strong temporal fluctuations of the contribution of each vertical mode visible on the wave regime of figure Fig. 8 may explain why they are large variations of E baro /E tot (in addition to the mean tendency of a decay with β) in the wave regime figure 7.

To conclude, in the turbulent stirring regime (i.e. L R ∼ L or L R > L) MRS statistical mechanics predictions are qualitatively correct, and there are no "transients" (no temporal fluctuations) once the flow is self-organized. By contrast, in the wave regime (i.e. L R ≪ L) , there are strong temporal fluctuations in the mode amplitudes, and we believe that statistical mechanics (either MRS of energy-enstrophy equilibria for the truncated dynamics) can not account for these observations.

Discussions and conclusions

We have used statistical mechanics to predict the flow structure resulting from the self-organization of freely evolving inviscid stratified quasi-geostrophic flows. The only assumption of the theory is that the flow explores evenly the available phase space. In order to compute explicitly statistical equilibria, and to discuss basic physical properties of these equilibria, we have focused on MRS equilibria characterized by linear relation between potential vorticity and streamfunction. For such MRS equilibrium states, we need the total energy and the fine-grained enstrophy in each layer, although all the invariants are implicitly taken into account.

We explained that these states are expected in a low energy limit or in a strong mixing limit. By applying a method proposed by [START_REF] Bouchet | Simpler variational problems for statistical equilibria of the 2d euler equation and other systems with long range interactions[END_REF], we have shown that MRS equi- L/2π, at the end of each run, for the same initial condition ψ0 (here L = 2π, E0 = 1, so only β varies). Red continuous line: predictions from statistical mechanics. The vertical structure of the equilibrium state is computed with equation (3.14), using K eff = 3.5, and using the enstrophy profile Z0(β) of the initial condition ψ0. The dashed-dotted line separates the wave regime from the turbulent stirring regime, in which we expect the statistical mechanics predictions to be valid. libria characterized by a linear qψ relation are solutions of a minimum coarse-grained enstrophy variational problem, which generalizes the Bretherton-Haidvogel minimum enstrophy principle to the stratified case.

The central result of the paper is the elucidation of the physical consequences of the conservation of the fine-grained enstrophy Z 0 (z) on the structure of the equilibrium states. We showed first that the affine coefficients between the coarse-grained PV field and the streamfunction are proportional to the fine-grained enstrophy: q = β t Z 0 (z)ψ. This relation allowed to compute explicitly the statistical equilibria, and to discuss their structure depending on the profile Z 0 (z). When this profile Z 0 is depth independent, the equilibrium state is fully barotropic, whatever the stratification profile N . When this profile Z 0 is surface intensified, the equilibrium state is an SQG-like mode, characterized by an e-folding depth ∼ f 0 /N K, where 1/K is bounded by the domain scale. Since the e-folding depth increases with K -1 , larger horizontal flow structures imply "more barotropic" flows. Since statistical equilibria are associated with the gravest Laplacian horizontal modes (K = 1), the equilibrium state is the "most barotropic one" given an enstrophy profile Z 0 .

We conclude that the dynamics leads to the gravest horizontal mode on the horizontal (inverse cascade), associated with the gravest vertical structure consistent with the conservation of the fine grained enstrophy. The flow becomes fully barotropic when the fine-grained enstrophy profile is depth independent.

These results can be used to understand the role of beta in barotropization in the following way: consider an initial SQG-like flow, with β = 0, prescribed by its streamfunction; its enstrophy profile is surface intensified, and so will be the equilibrium sate. But when one switches on the beta effect, with the same initial streamfunction, the contribution of the depth independent part of the fine-grained enstrophy increases, which means a tendency toward a more barotropic equilibria, according to the previous conclusion. This result reflects the fact that in physical space, the initial SQG-like mode stirs the interior PV field (initially a beta plane), which in turn induces an interior flow, which stirs even more the interior PV field, and so on.

When the beta effect becomes large enough, the PV distribution at each depth z becomes prescribed by the initial beta plane, so the most probable state should be also characterized by a depth independent PV field, i.e. by a barotropic flow.

To reach the equilibrium state, the PV field must be stirred enough to explore the phase space. Yet large values of beta stabilize the initial condition, the system is trapped in a stable state different from the equilibrium state, and the flow dynamics becomes dominated by the interaction between Rossby waves rather than by turbulent stirring. The reason for this trapping is that the dynamics can not provide by itself sufficient perturbations to escape from the stable state and to explore other regions of the phase space, which would attract the system toward the equilibrium state. We observed in numerical simulations that this wave regime did not lead to a barotropization of the flow after a few hundred turnover times. We do not know if still longer numerical runs would lead to more barotropization. We estimated the critical value of beta between the turbulent stirring regime and the Rossby wave regime as the value such that the Rhines scale is of the order of the domain scale, which gives β = 4π 2 E 1/2 0 /L 2 . In this paper, we focused on the effect of β on barotropization process. More generally, our results show that vertical transfers of energy and momentum are favored by the presence of any lateral potential vorticity gradients, since these gradients provide a source of available perturbation enstrophy. These lateral PV gradients may for instance be due to large scale mean flows set by external process (large scale wind pattern for the oceans, and large scale temperature gradients for the atmosphere). A similar situation occurs in the framework of interactions between waves (or eddies) and mean flows: horizontal momentum can be transfered vertically by inviscid form drag if there are horizontal fluxes of potential vorticity, which require themselves the presence of large scale potential vorticity gradients, see e.g. [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF].

Another source of fine-grained potential vorticity would be provided by the addition of bottom topography. This should in this case play against barotropization, since the topography induces potential fine-grained enstrophy in the lower layer only. In fact, an initially surface intensified flow may evolve towards a bottom trapped current above the topographic anomaly, which can be explained by the statistical mechanics arguments presented in this paper, see [START_REF] Venaille | Bottom trapped currents as statistical equilibrium states above topographic anomalies[END_REF].

We now discuss how these results may apply to the mesoscale ocean, and in what respect they may provide interesting eddy parametrization. First, we have considered the relaxation of an initial condition by an inviscid flow. Real oceanic flows are forced dissipated. Statistical mechanics predictions should apply if a typical time scale for selforganization is smaller than typical time scale of forcing and dissipation. Issues of forcing and dissipation lead to another difficulty: what is the domain in which the flow selforganizes ? Strictly speaking, it should be the oceanic basin. But forcing and dissipation become dominant at basin scale (with, for instance, the Sverdrup balance setting the gyre structures). In practice, on a phenomenological basis, one could still apply the statistical mechanics predictions by considering an artificial domain scale with prescribed scale of arrest for the inverse cascade, governed say by dissipation processes, but which could not be predicted by statistical mechanics itself.

We have also neglected bottom friction, because equilibrium statistical mechanics apply for non-dissipative systems. It is generally believed that bottom friction plays a key role in the vertical structure of oceanic flows [START_REF] Arbic | Baroclinically Unstable Geostrophic Turbulence in the Limits of Strong and Weak Bottom Ekman Friction: Application to Midocean Eddies[END_REF]; [START_REF] Thompson | Scaling Baroclinic Eddy Fluxes: Vortices and Energy Balance[END_REF]. On a phenomenological basis, quasi-statistical equilibria could be computed in the high bottom friction limit by adding a constraint of zero velocity at the bottom, for instance by considering ψ(-H) = 0 in equation (3.14).

Third, equilibrium statistical mechanics predict the final state organization of the freely evolving inviscid dynamics, but neither predicts the route toward this state nor the corresponding typical time scales for the convergence toward this state. Various parametrizations that relax the flow toward the equilibrium state, following a path of maximum entropy production have been proposed, see [START_REF] Kazantsev | Subgrid-Scale Eddy Parameterization by Statistical Mechanics in a Barotropic Ocean Model[END_REF]; [START_REF] Robert | Relaxation towards a statistical equilibrium state in twodimensional perfect fluid dynamics[END_REF]. These parametrization satisfy basic physical constraints satisfied by the dynamics (PV distribution conservation and energy conservation), but the actual dynamics may in some cases follow a different path than one of maximum entropy production.

For instance, [START_REF] Smith | The Scales and Equilibration of Midocean Eddies: Freely Evolving Flow[END_REF]; [START_REF] Fu | Nonlinear energy and enstrophy transfers in a realistically stratified ocean[END_REF] reported the existence of two time scales for the flow evolution in a freely evolving quasi-geostrophic turbulent flow with surface intensified stratification. On a fast time scale, the dynamics lead to a inverse cascade of vertical modes toward the first baroclinic modes, which is surface intensified, and observed the tendency toward barotropization at much larger time scale. The existence of these two time scales is a dynamical effect that can not be accounted by statistical mechanics.

The results presented in this paper have highlighted the important role of the finegrained enstrophy profile, for this imposes strong constraints for the eddy structures. On a practical point of view, from the perspective of eddy parametrizations, our result suggest that rather than assuming a vertical structure of eddies given by SQG modes, or given by combination of barotropic and first baroclinic modes, one could compute the eddy structure with equation (3.14), assuming that K is the scale of arrest, and Z 0 an enstrophy profile that could be deduced from the resolved flow. For instance, Z 0 could be set by the structure of the most unstable mode of local baroclinic instabilities.

Observations indicate the presence of long-lived surface intensified eddies in the ocean [START_REF] Chelton | Global observations of large oceanic eddies[END_REF][START_REF] Johnson | Equatorial Pacific 13C Water Eddies in the Eastern Subtropical South Pacific Ocean[END_REF][START_REF] Chaigneau | Vertical structure of mesoscale eddies in the eastern South Pacific Ocean: A composite analysis from altimetry and Argo profiling floats[END_REF], and we here speculate as to how these eddies may be interpreted according to a statistical mechanical theory. First, the eddies may be interpreted as local statistical equilibrium states of the continuously stratified dynamics associated with a surface intensified fine-grained enstrophy profile, much as we considered in the present paper. Second, they could be interpreted as statistical equilibrium states of a 1.5 layer quasi-geostrophic model assuming that only the upper layer of the ocean is active, see e.g. Venaille & Bouchet (2011a). Third, they could be far from equilibrium states driven by physical mechanisms preventing convergence towards statistical equilibrium (dynamical effects due to non-uniform stratification, bottom friction, large scale forcing). Understanding which hypothesis is relevant to explain the formation of long-lived surface intensified ocean vortices will require further research.

where E is given by (3.7), where the pdf field is locally normalized using (3.1). Equations (A 5) and (A 4) give α(σ, z) = σ 2 4(Z 0 -Z cg )

and ρ = exp -(σ-q) 2 4(Z0-Zcg)

(4π(Z 0 -Z cg ))

1/2

Injecting these expressions into (A 7) yields

G α [ρ] = 1 2 0 -H dz w Z cg Z 0 + cst, with w(x) = -ln(1 -x) + 1 1 -x
In the limit x ≪ 1, we have w(x) ≈ x, so finding the minimizer ρ m of G α with the (global) constraint on the energy and the enstrophy is equivalent of finding the minimizer of the functional 0 -H dz Z cg [q]/Z 0 with the same constraints. The hypothesis that Z cg ≪ Z 0 is consistent with the "strong mixing limit" (β t ψσ ≪ 1) or the low energy limit for which MRS equilibria are characterized by a linear relation between potential vorticity and streamfunction. We conclude that MRS equilibrium states characterized by a linear qψ relation are solutions of the variational problem (3.9).

Figure 1 .

 1 Figure 1. Left panel: three different fine grained enstrophy profiles. Right panel: corresponding vertical structure of statistical equilibrium states (ψ(z)/ψ(0) on the left panel), in the case of constant stratification ( f 2 0 /N 2 = 0.1). The e-folding depth in case A is h = f0/N K (with here K = 1 for the statistical equilibrium state).

Figure 2 .

 2 Figure 2. case β = 0 (SQG-like). Only the fields in upper, middle and lower layer are shown.

Figure 3 .

 3 Figure 3. qψ relation in the upper layer for t = 250 eddy turnover time.

Figure 4 .

 4 Figure 4. case β = 2.1. See Fig. 2 for legend. Note that the colors scales are changed at each time; this is why the interior " beta plane" is not visible in the upper panel of potential vorticity.

Figure 5

 5 Figure5. Relation between qβy and streamfunction in the upper layer, at t = 100 eddy turnover times, for β = 2.1. The fact that this relation is close to be linear shows that the small energy limit considered in this paper is relevant to this case. The coefficient of the slope is approximatively one, meaning that the energy is condensed in the gravest horizontal Laplacian eigenmode, which is propagating westward at the speed of barotropic waves. This is consistent with the computation of minimum enstrophy states on a beta plane, see subsection 3.4.

Figure 6 .

 6 Figure 6. case β = 10, see Fig. 2 for legend.

Figure 7 .

 7 Figure 7. Variation of the Barotropic energy E baro normalized by the total energy E0, varying the ratio of the domain scale with the Rhines scale L/LR = β 1/2 E -1/4 0

Figure 8 .

 8 Figure 8. Temporal evolution of the contribution of the energy of each baroclinic mode to the total energy, for different values of β. The modes are indexed by m, with m = 0 the barotropic mode and m = 1 the first baroclinic mode.

† By catalyst we do not mean a dynamical effect, as in chemistry, but simply the fact that barotropization in the final state is enhanced in presence of a beta plane.

† Note that with this notation, Em=0 = E baro is the energy of the barotropic mode, which is in general different from the total energy of the initial condition E0.

Acknowledgments This work was supported by DoE grant DE-SC0005189 and NOAA grant NA08OAR4320752. The authors sincerely thank Peter Rhines and two other reviewers for their very helpful comments provided during the review process. The authors also warmly thank F. Bouchet and J. Sommeria for fruitful discussions, and K.S. Smith for sharing his QG code.

Appendix A. Minimum enstrophy states

The aim of this appendix is to show that the class of MRS statistical associated with a linear qψ relation can be computed by solving a minimum enstrophy variational problem satisfied by the coarse-grained state q. This is a direct consequence of a more general result of [START_REF] Bouchet | Simpler variational problems for statistical equilibria of the 2d euler equation and other systems with long range interactions[END_REF], see also [START_REF] Majda | Nonlinear Dynamics and Statistical Theories for Basic Geophysical Flows[END_REF]; [START_REF] Naso | Statistical mechanics of two-dimensional Euler flows and minimum enstrophy states[END_REF].

Here bottom topography is omitted to simplify the presentation, but taking into account this additional parameter is straightforward. In addition, the results obtained in this appendix are valid whatever the domain geometry. The only complication added by the presence of boundaries is that the circulation at these boundaries becomes a parameter that must be discussed.

A.1. Variational problem of the MRS statistical theory

We consider the variational problem (3.8) given by MRS theory. To compute its critical points, we introduce the Lagrange multipliers β t ,α(σ, z), ζ(x, y) associated respectively with the energy E (3.7), with the global PV distribution P (3.4), and with the local normalization constraint N given by (3.1). Critical points are solutions of

where first variations are taken with respect to the field ρ. This equation and the normalization constraint lead to

The coarse-grained potential vorticity field associated with this state is

and the local fluctuations of potential vorticity are

where f ′ α is the derivative of f α .

A.2. Case of a linear qψ relation

Here we focus on a particular class of solutions, assuming

where the coefficient Ω(z) will be expressed in term of the constraints of the problem. Equations (A 2) and (A 3) give the qψ relation and the fluctuations

Integrating the fluctuations on the horizontal leads to

Using the relations (A 1-A 2), we find that

A.3. Link with a minimum enstrophy problem Solution of dual, less constrained variational problem is always a solution of the more constrained problem for some values of the constraints [START_REF] Ellis | Large deviation principles and complete equivalence and nonequivalence results for pure and mixed ensembles[END_REF]. We relax the constraints on the PV distribution by introducing the Legendre Transform G α of the entropy associated with this constraint, taking into account the conservation of the total energy E 0 and of the fine-grained enstrophy Z 0 (z) in each layer:

G (E 0 , Z 0 (z), α) = min