{'fraction_non_alphanumeric': 0.0364541122313353, 'fraction_numerical': 0.013210360813045028, 'mean_word_length': 4.738002232142857, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transport and emergent stratification in the equilibrated Eady model: the vortex gas scaling regime

We numerically and theoretically investigate the Boussinesq Eady model, where a rapidly rotating density-stratified layer of fluid is subject to a meridional temperature gradient in thermal wind balance with a uniform vertically sheared zonal flow. Through a suite of numerical simulations, we show that the transport properties of the resulting turbulent flow are governed by quasi-geostrophic (QG) dynamics in the rapidly rotating strongly stratified regime. The 'vortex gas' scaling predictions put forward in the context of the two-layer QG model carry over to this fully 3D system: the functional dependence of the meridional flux on the control parameters is the same, the two ajustable parameters entering the theory taking slightly different values. In line with the QG prediction, the meridional buoyancy flux is depth-independent. The vertical buoyancy flux is such that turbulence transports buoyancy along isopycnals, except in narrow layers near the the top and bottom boundaries, the thickness of which decreases as the diffusivities go to zero. The emergent (re)stratification is set by a simple balance between the vertical buoyancy flux and diffusion along the vertical direction. Overall, this study demonstrates how the vortex-gas scaling theory can be adapted to quantitatively predict the magnitude and vertical structure of the meridional and vertical buoyancy fluxes, and of the emergent stratification, without additional fitting parameters.


Introduction
The oceans and atmospheres of planets and their satellites are shallow fluid layers subject to the combined effects of rapid global rotation and strong density stratification. Planetary atmospheres are set into motion by differential heating in the meridional (North-South) direction -either as a result of radiative heating by a nearby star or because of heat coming from the deep planetary interior -while oceans are also subject to mechanical forcing by predominantly zonal (East-West) atmospheric winds.
In both cases a meridional temperature gradient emerges and coexists with a vertically sheared zonal flow as a result of thermal wind balance. This base state is unstable, however, and turbulent motion rapidly arises as a consequence of baroclinic instability (Pedlosky 1979;Salmon 1998;Vallis 2006). The resulting 'baroclinic' turbulence enhances heat transport in the meridional direction, which in turn affects the equilibrated meridional temperature profile, but also the local density stratification of the fluid layer. Baroclinic turbulence strongly contributes to meridional heat transport at midlatitudes in Earth's atmosphere and other planetary atmospheres, like Jupiter's, and plays a key role in setting the planets' climate (Liu & Schneider 2010;Read et al. 2020). It is also a dominant feature of ocean currents, most notably in the Southern Ocean, where baroclinic instability of the Antarctic Circumpolar Current (ACC) flowing around Antarctica leads to enhanced meridional transport (Nowlin Jr & Klinck 1986;Marshall & Radko 2003;Volkov et al. 2010). The resulting heat and salt transport sets the stratification of the Southern Ocean and to some extent that of all ocean basins (Wolfe & Cessi 2010;Nikurashin & Vallis 2011, 2012.
One crucial characteristic of baroclinic turbulence in the Southern Ocean is that its horizontal integral scale is much shorter than the transverse extent of the ACC (a similar scale separation arises in Jupiter's atmosphere but not in the Earth's atmosphere). The scale separation can be leveraged to parameterise the meridional transport in the form of a diffusive closure, where the diffusivity is inferred from a 'local' model consisting in an isolated patch of fluid much smaller than the size of the ocean basin but much larger than the integral scale of the flow (Held 1999). A hierarchy of local models can be considered: in the simplest instances the patch of fluid is modelled as a superposition of two vertically invariant layers of fluid stacked in the vertical direction. The equations of motion are further simplified by considering the quasi-geostrophic (QG) approximation to the dynamics of rapidly rotating strongly stratified fluid layers. This leads to the two-layer quasi-geostrophic (2LQG) model, where the governing equations reduce to a conservation equation for the potential vorticity (a scalar field) inside each fluid layer (Phillips 1954). Even for this canonical model, however, the transport properties of the equilibrated turbulent flow in the moderate-to-low drag regime have been captured only recently by a scaling theory, which we coined the 'vortex gas' scaling regime and extended to the β-plane (Gallet & Ferrari 2020, 2021. The 2LQG model is a rather crude description of the flow, and one may wonder whether any of these scaling-laws carry over to a fully 3D model of a patch of ocean or atmosphere experiencing weak bottom drag (Rivière et al. 2004). A 3D model also leads to key additional questions that must be addressed to design a skillful parameterisation of turbulent transport: beyond the meridional buoyancy flux, what is the magnitude of the vertical buoyancy flux, i.e., what is the direction of the eddy-induced buoyancy current vector? In a model with a continuous vertical direction, can we predict the vertical structure of the meridional and vertical buoyancy fluxes? Finally, how does this vertical buoyancy flux feed back onto the vertical density stratification?
In the present study we thus relax both the finite-layer and the quasi-geostrophic approximations, i.e., we consider a three-dimensional patch of rapidly rotating densitystratified fluid described by the Boussinesq equations. More precisely, we focus on the Boussinesq Eady model, where a uniform meridional temperature gradient coexists with a uniform vertically sheared zonal flow (Eady 1949). We restrict attention to the large-Richardson-number regime (buoyancy frequency much greater than the vertical shear) where perturbations amplify through a 'geostrophic instability' mechanism, see Stone (1971). The departures from the base state thus rapidly evolve into 3D turbulence that we solve using periodic boundary conditions in the horizontal directions. This Boussinesq Eady model is introduced in section 2. We present a suite of numerical simulations of the Eady model in section 3. We show that the numerical results are consistent with QG dynamics. We then recall the prediction of the vortex gas scaling theory for the magnitude of meridional transport and we validate this prediction against the numerical data. In section 4, we extend these predictions to the vertical buoyancy flux and emergent vertical stratification, finding good agreement with the numerical data. In section 5, we turn to the vertical structure of these quantities, which we predict to be depth-invariant within the low-diffusivity QG framework. Based on the Boussinesq data, we show that they are indeed all depth-independent in the bulk of the layer within the low-diffusivity QG regime, but that extremely low diffusivities are needed for that depth-independence to hold near the top and bottom boundaries. The turbulent kinetic energy profile gradually becomes depth-independent as one enters the low-friction vortex-gas regime, in line with a barotropization of the flow. We summarize the results and conclude in section 6.

Base state
We consider the Boussinesq equations for a rotating density-stratified layer of fluid inside a 3D domain (x, y, z) ∈ [0, L] × [0, L] × [0, H], with gravity and rotation along the vertical axis z: where u(x, y, z, t) and θ(x, y, z, t) denote the velocity and (potential) temperature fields, p(x, y, z, t) is the kinematic pressure, α denotes the thermal expansion coefficient, g denotes gravity. The problem having very different characteristic length scales in the horizontal and vertical directions, we use two different viscosities ν z and ν ⊥ for the vertical and horizontal derivatives arising in the viscous term (i.e., for the ∂ zz and ∆ ⊥ = ∂ xx +∂ yy terms, respectively). Similarly, we use two different thermal diffusivities for horizontal (ν b;⊥ ) and vertical (ν b;z ) thermal diffusion. For simplicity, temperature θ is considered to be the only stratifying agent in the equations above. We focus on a steady base state with a vertically sheared zonal flow u = Sz e x in thermal-wind balance with a temperature gradient along the meridional direction y. Substitution into the equations above leads to the following base state: The base state is sketched in figure 1. The meridional temperature gradient in (2.3) is directly proportional to the vertical shear S as a result of thermal wind balance, and we have included a uniform background vertical temperature gradient (uniform background buoyancy stratification N 2 ). In the following, we thus consider base states both with and without background vertical stratification (N = 0 and N = 0, respectively). Strictly speaking, the base state (2.3) is a solution to equations (2.1-2.2) only if supplemented with the right boundary conditions. First, we neglect the tiny Ekman layer connecting the shear flow (2.3) to the stress-free upper boundary. Secondly, a weak vertical heat flux needs to be imposed at the boundaries to maintain a base state with N = 0, the magnitude of which vanishes as the vertical diffusivities tend to zero. We validate this approach a posteriori by showing some solutions with N = 0. For these solutions a strong density stratification compatible with the top and bottom insulating boundary conditions emerges, and the resulting equilibrated state follows the same scaling behavior as solutions with strong imposed background stratification N = 0. In line with the theoretical predictions to come, the emergent stratification is uniform, which justifies the compatibility between simulations performed with and without an imposed uniform background stratification N 2 .

Departure from the base state: governing equations
Consider arbitrary departures from the base state: u(x, y, z, t) = Sz e x + v(x, y, z, t) , (2.4) αgθ(x, y, z, t) = −f S y + N 2 z + b(x, y, z, t) , (2.5) where b denotes the buoyancy departure. The velocity departure v = (u, v, w) is divergence-free in the Boussinesq framework, ∇ · v = 0. Substituting (2.4-2.5) into equations (2.1-2.2) leads to the evolution equations for v and b: We consider periodic boundary conditions for the departure fields v and b in the horizontal directions. The top and bottom boundaries are no-flux: and v satisfies free-slip boundary conditions at the top: where v ⊥ denotes the horizontal components of v. Strictly speaking, the bottom boundary condition should be no-slip, v| z=0 = 0. Such a no-slip boundary condition induces strong friction and pumping associated with the bottom Ekman boundary layer, of typical thickness ν z /f . However, the picture of a laminar Ekman layer acting on a perfectly flat ocean bottom is probably a very naive representation of the turbulent boundary layer near the rough ocean floor. We thus turn to a bottom boundary condition that allows us to vary the magnitude of bottom friction independently of the vertical viscosity. We replace the no-slip bottom boundary condition by the following frictional boundary condition: This boundary condition corresponds to a linear drag law, somewhat similar to but weaker than pure Ekman friction. The strength of bottom friction can be tuned through the value of the friction coefficient κ (not to be confused with the buoyancy diffusivities ν b;⊥ and ν b;z ). As compared to a standard no-slip boundary condition, (2.10) allows us to study the regime of weak bottom drag while keeping vertical viscosity values compatible with Direct Numerical Simulation (DNS). The boundary condition (2.10) reduces to a standard no-slip boundary condition for κ → ∞, while it reduces to a stress-free boundary condition for κ → 0. For arbitrary κ, the flow near the bottom of the fluid domain resembles a truncated Ekman spiral, with reduced shear and reduced pumping as compared to the standard Ekman boundary layer over a no-slip boundary. In appendix A, we compute the total damping and pumping induced by the boundary condition (2.10) on the bulk flow.

Non-dimensionalization
We non-dimensionalize the equations using the length scale H and the time scale f −1 : Dropping the ♯ symbols in the following to alleviate notations, the dimensionless governing equations are: where the Rossby number is defined as Ro = S/f . The dimensionless boundary conditions are: (2.17) 2.4. Numerical implementation We wish to characterize the buoyancy transport achieved by statistically steady solutions to equations (2.13-2.14) with the boundary conditions (2.15-2.17) (see Bachman & Fox-Kemper (2013) for a study of the spin-down problem). Denoting as · a time and volume average, we will consider the time-and volume-averaged meridional and vertical buoyancy fluxes, vb and wb , respectively. We also consider the vertical profiles of the meridional and vertical buoyancy fluxes, vb(z) and wb(z), where the overbar denotes a horizontal area average together with a time average. Finally, a quantity of interest is the emergent stratification b(z) that arises in the equilibrated state. The total vertical buoyancy stratification is (N/f ) 2 z + b(z), and, assuming that the emergent stratification is approximately uniform, we define the Rossby deformation radius λ as: where λ is non-dimensionalized with H. We have performed a suite of numerical simulations of equations (2.13-2.14) using Coral, a pseudo-spectral, scalable, time-stepping solver for differential equations (Miquel 2021). In Coral, the top and bottom boundary conditions are imposed through basis recombination. The variables are expanded on bases of functions obtained as tensor products of Fourier modes along the horizontal and suitable linear combinations of Chebyshev polynomials, each of which obeys the boundary conditions along z. The divergence-free constraint is readily implemented in Coral by introducing the standard toroidal and poloidal velocity potentials ψ(x, y, z) and φ(x, y, z), respectively. The frictional bottom boundary condition is dealt with by defining modified velocity potentials ψ and φ as ψ(x, y, z) = P ψ (z) ψ(x, y, z) and φ(x, y, z) = P φ (z) φ(x, y, z), where P ψ (z) and P φ (z) are carefully chosen quadratic polynomials, the coefficients of which depend on the drag coefficient. In particular, these polynomials are tailored so that imposing the Robin boundary conditions (2.16) on the velocity field amounts to imposing standard stressfree type boundary conditions for the modified potentials: ∂ z ψ = φ = ∂ zz φ = 0. We can thus solve for ψ and φ using standard Galerkin basis recombination. The code has been benchmarked by careful comparison with results from linear stability analysis, checks that the various power integrals are exactly satisfied in the nonlinear regime, and comparison of the overall buoyancy transport against a run of the same problem performed with the Dedalus solver (Burns et al. 2020). The numerical runs are initialized with either low-amplitude noise or a checkpoint from a previous simulation. We average the relevant quantities over the statistically steady regime, and make sure that the statistics are correctly converged by comparing with averages performed over the first half of the statistically steady signal only.

The low-diffusivity regime
We focus on the regime where the horizontal extent of the domain is large compared to the energy-containing scale and the mixing-length of the flow (estimated as the ratio of root-mean-square horizontal buoyancy fluctuations divided by the meridional buoyancy gradient, see e.g. Thompson & Young (2006)). Through variations of the domain size, we checked that the transport properties of the flow are then independent of the horizontal domain size (for brevity, we do not report these validation runs in the Tables below). The diffusivities are chosen to be equal for momentum and buoyancy, i.e., E b;⊥ = E ⊥ and E b;z = E z in all of the numerical runs. This corresponds to setting the Prandtl number to unity, an appropriate choice in the present context where the viscosities and diffusivities represent mixing by small-scale turbulence, and not molecular processes. The horizontal diffusivities are small enough not to affect the transport properties of the equilibrated flow. The role of the vertical diffusivities is more subtle: anticipating the results in sections 3 and 4, we choose small enough vertical diffusivities for the transport properties of the flow to be well-described by bulk diffusivity-free QG dynamics with bottom friction. In that regime, however, the emergent stratification depends on the vertical buoyancy diffusivity, as described in section 4.2. Table 1 lists the various control parameters for the numerical runs considered in the present study: we vary the background stratification, the Rossby number, the diffusivities and the friction coefficient.
The qualitative aspect of the equilibrated flow is illustrated in Figure 2, where we show a 3D rendering of the total buoyancy field together with horizontal slices of the dimensionless vertical vorticity ζ at mid-depth and at the top surface. We observe b(x, y, z, t) . Equilibrated baroclinic turbulence illustrated with snapshots from runs 21 (top, low drag) and 3 (bottom, moderate drag). a: 3D rendering of the total buoyancy field. Notice the very different vertical and horizontal scales, and thus the very weak isopycnal slopes. b: the vertical vorticity at mid-depth shows coherent vortices, with slight cylcone/anticyclone asymmetry in this marginally QG run. c: The surface vertical vorticity exhibits sharp frontal structures where relative vorticity greatly exceeds planetary vorticity. d and e: same as panels b and c, but for larger dimensionless drag κ * in a strongly QG regime. The vorticity at mid-depth shows a less dilute vortex gas than in panel b. QG cyclone/anticyclone symmetry is well satisfied, including in the sharp SQG frontal structures visible in panel e. that the flow remains strongly stratified in the equilibrated state, with weak isopycnal slopes in the bulk of the fluid domain. Panels a, b and c correspond to a run that is marginally QG and has very low dimensionless drag. The bulk vertical vorticity exhibits a relatively dilute gas of vortices, although there is some asymmetry between cyclones and anticyclones. This asymmetry probably stems from the sharp non-QG frontal structures that develop at the surface (Hoskins & Bretherton 1971;Klein et al. 2008;Ragone & Badin 2016;Siegelman et al. 2020), see panel c. We will see in the following that the transport properties of this flow are correctly described by QG theory despite the slight cyclone/anticyclone asymmetry of the interior flow. Panels d and e illustrate a run in the strongly QG regime with moderately low drag. The vortex gas is less dilute as a result of the larger drag. QG cyclone-anticyclone symmetry is wellsatisfied both at mid-depth (panel d) and in the surface quasi-geostrophic (SQG) fronts visible in panel e.
In Figure 3, we show the overall meridional buoyancy flux vb as a function of the friction coefficient κ. The buoyancy flux varies over more than four orders of magnitude over the entire data set. Meridional transport increases as the friction coefficient decreases (for otherwise constant parameters), and it decreases rapidly with decreasing Rossby number.

Bulk QG dynamics: further reducing the number of control parameters
The qualitative features of the equilibrated flows in Figure 2 point to QG dynamics. One can derive the QG approximation to the Boussinesq Eady system following the standard procedure reviewed in appendix B (see also Salmon (1998); Vallis (2006)). Decompose the pressure field into a time and horizontal mean p(z) plus fluctuations p(x, y, z, t): p(x, y, z, t) = p(z) +p(x, y, z, t), Notice that the horizontal and time derivatives of p andp are equal, so that p andp can be used interchangeably in many (but not all) of the expressions to come. Assuming that the emergent stratification is approximately uniform and neglecting the diffusive terms in the bulk of the flow, the potential vorticity (PV) conservation equation reads (see appendix B): The boundary conditions are the buoyancy equation written at the boundaries, remembering the hydrostatic relation b = ∂ z p. At the top free-surface one obtains: The second boundary condition stems from the buoyancy equation written just above the partial Ekman layer connecting the bulk flow to the frictional bottom boundary condition (a height denoted as z = 0 + ): where w| 0 + denotes the vertical pumping velocity induced by the partial Ekman layer. As shown in appendix A, this pumping velocity takes the form: where the effective friction coefficient acting on the interior flow is: both κ and κ eff being non-dimensionalized with f . Expression (3.7) reduces to the standard Ekman friction coefficient in the no-slip limit κ → ∞ (Pedlosky 1979). The buoyancy equation at z = 0 + finally reads: QG dynamics is governed by the conservation equation (3.2) for the bulk PV (3.3), with the top and bottom boundary conditions (3.4) and (3.8). The fluid communicates in the vertical direction through slight vertical displacements of the isopycnals, which induce modifications in vertical vorticity ζ = ∆ ⊥ p through the conservation of PV (for the specific case of a cylindrical fluid column contained between two isopycnals, this process can be understood as the conservation of angular momentum as the column gets compressed or stretched). Diffusion plays no roles in this process, the consequence being that κ and E z enter the QG dynamics only through the effective friction coefficient κ eff associated with pumping in the bottom (partial) Ekman layer. The number of dimensionless parameters relevant to the QG regime can be further reduced through the introduction of the QG scalings: Through this change of variables the QG equations become: An instructive feature of this set of equations is that it involves a single control parameter: which is the Boussinesq equivalent of the relevant dimensionless drag coefficient in the QG framework (Thompson & Young 2006;Gallet & Ferrari 2020). The equilibrated solutions to the set of equations (3.10-3.13) are characterized by a meridional buoyancy transport D * = ∂ X P ∂ z P = vb /(Ro 2 λ) that depends only on κ * . We conclude that: where the unknown function F denotes the dependence of the QG buoyancy eddy diffusivity D * on the QG dimensionless drag κ * . As discussed in Thompson & Young (2006) and Gallet & Ferrari (2020), the 2LQG equivalent of D * readily quantifies the ratio of the meridional buoyancy flux over the meridional buoyancy gradient, hence the name 'diffusivity'. Indeed, while the change of variables above provides a rigorous way to identify D * and κ * as the 'order' and 'control' parameters of the QG regime, the form of D * and κ * could have been guessed from studies of the strongly idealized 2LQG model: as shown in Thompson & Young (2006), the 2LQG model on the f -plane reduces to the study of D * versus κ * , where D * is the meridional flux divided by the product of the deformation radius with the squared shearing velocity, while κ * is the bottom drag coefficient multiplied by the deformation radius and divided by the shearing velocity. Equation (3.15) suggests a quantitative way to assess whether meridional buoyancy transport is governed by QG dynamics in the present simulations of the Boussinesq system: in Figure 4, we plot D * as a function of κ * for the entire dataset. This representation leads to a collapse onto a master curve of the entire data set, which includes various values of the Rossby number, background stratification and vertical Ekman numbers. This confirms the validity of the QG scalings for bulk transport in the present numerical data, despite the emergence of potentially non-QG frontal dynamics near the boundaries.
We further compare the results of the Boussinesq system to QG theory by solving the QG Eady problem explicitly. A trivial solution to the QG PV conservation equation (3.10) is Q = 0. In appendix B, we recall how this constraint leads to an inversion relation providing P | 1 and P | 0 + in terms of ∂ z P | 1 and ∂ z P | 0 + at each time step. Solving the 3D QG problem then reduces to solving the coupled equations (3.12) and (3.13) for the 2D fields ∂ z P | 1 and ∂ z P | 0 + , using the inversion relation to infer P | 1 and P | 0 + at each time step. The approach is very similar to the surface quasi-geostrophy (SQG) (Blumen 1978;Held et al. 1995;Callies et al. 2016;Lapeyre 2017), the difference being that there are two horizontal boundaries in the present system (at top and bottom), as opposed to a single horizontal boundary in the standard SQG system. The second boundary crucially allows for the emergence of depth-invariant barotropic eddies in the present system. The effectively 2D equations are solved on a GPU in a large-enough domain for finite-size effects to be negligible. To damp the small-scale filaments, we also add hyperviscous terms −µ∆ 4 X,⊥ ∂ z P | 1 and −µ∆ 4 X,⊥ ∂ z P | 0 + to the right-hand sides of (3.12) and (3.13), respectively, where the dimensionless hyperdiffusivity µ is chosen small enough to not affect the meridional heat flux in the equilibrated state. The resulting data points are plotted using crosses in Figure 4. They agree very well with the 3D Boussinesq data, indicating that PV vanishes in the interior of the 3D domain as a result of the dissipative terms, with negligible PV injection by frontal structures at the boundaries (in contrast  Figure 4. Dimensionless eddy-induced diffusivity as a function of the dimensionless effective friction coefficient, using the QG non-dimensionalization. The entire data set falls onto a master curve, indicating that QG dynamics governs the overall buoyancy transport. Same symbols and colours as in Figure 3. The thick vertical dashed lines indicate the predicted values of κ * according to equation (4.9) for the runs without background stratification (black circles). The predicted values of κ * agree with the DNS ones within 5%. Black crosses are the solutions to the QG system. The solid line is the vortex-gas prediction (3.16) with c1 = 0.32 and c2 = 0.61.
To summarize, we have established that meridional transport in the equilibrated Eady Boussinesq system is governed by QG dynamics, as illustrated by the collapse of the Boussinesq data in the representation of Figure 4 and their good agreement with solutions of the QG system. Developing a theory for the magnitude of the meridional transport in the equilibrated Eady problem then reduces to the determination of the function F in the scaling relation (3.15).

Vortex gas scaling theory
We recently introduced such a scaling theory, coined the vortex gas scaling regime (Gallet & Ferrari 2020), within the framework of the 2LQG model. The central idea is that the barotropic flow consists in a dilute gas of isolated vortices wandering around through mutual advection (Thompson & Young 2006;Carnevale et al. 1991). The vortex cores are small compared to the inter-vortex distance, with a radius comparable to the Rossby deformation radius. Buoyancy fluctuations arise through the distortion of the background meridional gradient by the vortical flow. Focussing on a single barotropic vortex dipole shows that the mixing length is comparable to the inter-vortex distance, while the diffusivity is given by the product of the inter-vortex distance with the velocity at which the vortex cores move. Two energetic arguments complete the scaling theory: first, a typical velocity scale is estimated by equating the potential energy drop with the final kinetic energy as a fluid column travels over one mixing-length (seen as the mean-free path of the turbulent fluid motion). This 'slantwise free-fall' argument yields the typical velocity of a fluid column located between two vortices, which is also the velocity scale at which the vortex cores move. Second, at equilibrium the energy released through baroclinic instability must be balanced by frictional dissipation. Significantly larger velocities arise at the immediate periphery of the vortex cores, where fluid particles circle rapidly without transporting much buoyancy. These large velocities are associated with large frictional dissipation, with important consequences for the energy power integral, which is the last scaling relation of the theory. This line of arguments provides a theoretical expression for the function F in (3.15). Specifically, for linear bottom friction we obtain: ≃ 0.72 (Gallet & Ferrari 2021). For the present Boussinesq Eady data, we show in Figure 4 that the vortex gas prediction (3.16) accurately captures the master curve, with c 1 = 0.32 and c 2 = 0.61. This indicates that the physical intuition gathered from the 2LQG model carries over to the fully 3D Eady model.

Vertical transport and emergent stratification
As compared to the 2LQG system, the Eady model includes a continuous vertical direction. The eddy-induced buoyancy current is a vector with both a meridional and a vertical component. The vertical buoyancy flux is crucial because it plays a central role in setting the emergent vertical stratification. In the following we propose simple scaling arguments to quantitatively predict the vertical buoyancy flux before computing the emergent stratification. We will show that the latter can be negligible, comparable, or much greater than the background stratification depending on the parameter regime. In particular, we will show that restratification can be strong enough to induce QG dynamics even in the absence of an imposed background vertical stratification.

Bulk eddy transport is along isopycnals
In a similar fashion to (3.1), decompose the buoyancy field into a time and horizontal mean b(z) plus fluctuationsb(x, y, z, t): b(x, y, z, t) = b(z) +b(x, y, z, t), withb = 0 . (4.1) Upon multiplying the buoyancy equation (2.14) with b before averaging temporally and horizontally, assuming that the fluctuations are much weaker than the mean,b ≪ b, one obtains: If the diffusive terms on the right-hand side are negligible, we recover the standard conclusion that the eddy-induced buoyancy current is along isopycnals (Gent & Mcwilliams 1990;Young 2012). In terms of the overall magnitude of the vertical buoyancy flux, this  Figure 5. Inset: Ratio r between the two sides of the approximate relation (4.3). This ratio is always close to unity, with a typical value r ≃ 0.85 over the entire data set. This ratio approaches unity as the diffusivities are lowered, indicating along-isopycnal transport. Main figure: Dimensionless vertical buoyancy flux as a function of the dimensionless effective friction coefficient, both being non-dimensionalized following the QG scalings. Combining the vortex-gas theory with pure along-isopycnal transport leads to the prediction on the right-hand side of (4.4), shown as a solid line. Same symbols and colours as in Figure 3.
leads to the simple relation: the slope of the isopycnals being Ro/λ 2 . In Figure 5 we plot the ratio of the left-hand side over the right-hand side of equation (4.3), denoted as r, for solutions of the Boussinesq system. This ratio is indeed approximately constant at low drag, with a value r ≃ 0.85 close to unity, the slight departure from one being due to the small diffusive contributions in (4.2). The ratio r approaches unity if the diffusive terms are small enough in a regime where the assumptions of QG are well satisfied. Solutions to the QG system -equations (3.12-3.13) together with inversion relations (B 17-B 18) -satisfy r = 1 exactly. Overall, the combination of the vortex-gas scaling prediction with the Gent-McWilliams argument leads to a good theoretical prediction for the vertical buoyancy flux: as illustrated in Figure 5.

Strength of the emergent stratification
One then estimates the emergent stratification from the vertical buoyancy flux through the potential energy evolution equation. Multiplying the buoyancy equation (2.14) by z before averaging over space and time yields, after using w = 0 and performing some integrations by parts using the boundary conditions: (4.5) We wish to show that the term involving the mean meridional flow v(z) is negligible as compared to the vertical buoyancy flux on the left-hand side. This mean meridional velocity corresponds to an ageostrophic flow, because a mean zonal buoyancy gradient maintaining v(z) through thermal wind balance would be incompatible with the periodic boundary conditions. The profile v(z) can be estimated by time-and area-averaging the zonal component of the momentum equation (2.6): The second term on the right-hand side is negligible in the bulk of the domain because E z ≪ 1. The first term on the right-hand side involves wu. It leads to a contribution −Ro wu to the left-hand side of (4.5), to be compared with the second term, wb . A crude estimate of the ratio of these two contributions is Ro u 2 1/2 / b 2 1/2 , which using thermal-wind balance for a flow of horizontal scale λ is of order Ro/λ ≪ 1. This simple estimate is a first indication that the v term in equation (4.5) is negligible as compared to the vertical buoyancy flux.
We stress the fact, however, that wu and v are even smaller than the crude estimate above. Indeed, to lowest order the magnitude of wu can be estimated using QG theory, the diagnostic vertical velocity being inferred from equation (B 7). One can then leverage the cyclone-anticyclone symmetry of the QG system to show that the resulting wu vanishes, as explained in appendix B.3. In other words, while QG dynamics induces a nonzero vertical buoyancy flux wb in the bulk of the domain, it leads to a vanishing vertical flux of zonal momentum wu for symmetry reasons. Nonzero wu in the bulk may only arise at next order in the QG expansion, where cyclone-anticyclone symmetry is broken (Muraki et al. 1999;Hakim et al. 2002;Gallet et al. 2014). We conclude that wu is smaller than the estimate of the previous paragraph by (at least) a factor corresponding to the isopycnal slope, Ro/λ 2 ≪ 1. The v term is thus fully negligible as compared to wb in equation (4.5).
Having shown that the dominant balance in equation (4.5) is between the vertical buoyancy flux and the diffusive flux associated with the emergent stratification, we obtain the following expression for the emergent stratification: One should notice that the vertical buoyancy diffusivity comes back into play through the coefficient E b;z in the denominator. While the vertical momentum diffusivity (the viscosity coefficient) affects the scaling behavior of the system only marginally, through a modification of the effective friction coefficient (3.7), the vertical buoyancy diffusivity greatly impacts the emergent stratification (4.7). The transport properties of the system are independent of the vertical buoyancy diffusivity only in the regime where the emergent stratification (4.7) is negligible as compared to the imposed background stratification. Otherwise, changes in E b;z induce changes in the emergent stratification (4.7) and thus in the Rossby deformation radius λ, the consequence being that the transport properties of the equilibrated state are shifted following the QG master curves in Figures 4 and 5. From equation (4.7), we can estimate the emergent stratification directly in terms of the control parameters of the problem. While this can be done for an arbitrary background stratification, we focus on two limiting situations: the case without imposed background  (0), over the background one, (N/f ) 2 , versus the ratio of the right-hand side of (4.10) over E b;z . In agreement with the criterion (4.10), the emergent stratification is negligible as compared to the background one when the abscissa is much less than unity. In this regime the abscissa readily provides the theoretical prediction for the emergent stratification divided by the background one, and indeed the data points fall onto the diagonal. The emergent stratification becomes comparable to the background one as the abscissa approaches unity. Same symbols and colours as in Figure 3.
stratification, and the case where the emergent stratification is negligible as compared to the imposed background one.
In the absence of imposed background stratification, where b(1) − b(0) = λ 2 , we insert expression (4.4) for the vertical flux into (4.7), before writing the resulting relation as: (4.8) We invert this relation to obtain c 2 /κ * , which leads to the following expression for the emergent Rossby deformation radius (emergent stratification): where W denotes the Lambert function. These equations relate emergent quantities (on the left-hand side) to control parameters only (on the right-hand side). The first expression in (4.9) provides the emergent stratification, while the second one readily gives the dimensionless friction κ * arising in the QG framework (the abscissa in Figures 4 and  5). We have performed two numerical runs without background stratification to validate these predictions. In Figure 4, we show that the resulting meridional flux falls onto the master curve, at an abscissa given by the theoretical expression (4.9) within 5% accuracy. We now consider the case of a nonzero imposed background stratification (N/f ) 2 . We wish to determine when the emergent stratification is negligible as compared to the imposed one. To wit, we estimate the emergent stratification by inserting into (4.7) expression (4.4) for the vertical flux, in which we replace λ with N/f . Demanding that

Buoyancy fluxes
Kinetic energy Mixing for decreasing diffusivities a.
b. this emergent stratification be negligible as compared to the imposed one, (N/f ) 2 , then amounts to satisfying the inequality: The criterion (4.10) is expressed in terms of the control parameters of the system only. It allows one to assess a priori whether or not the emergent stratification is negligible as compared to the background one, as illustrated in Figure 6. When the inequality (4.10) is satisfied, the emergent stratification is negligible and the meridional and vertical buoyancy fluxes are readily obtained by replacing λ with (N/f ) in expressions (3.16) and (4.4).
To summarize, combining the vortex-gas scaling theory with both the along-isopycnaltransport argument and a simple advective-diffusive buoyancy-flux balance along the vertical direction, we obtain theoretical expressions for the meridional buoyancy flux, the vertical buoyancy flux and the emergent stratification in terms of the control parameters of the problem.

Vertical structure
Having established that the 2LQG model has skill to predict the functional dependence of meridional transport in the Eady problem, we turn to the vertical structure of the buoyancy fluxes, a question that can be addressed only by going beyond 2LQG. An accurate description of the vertical structure of the transport tensor is crucial step in the development of a parameterization of baroclinic turbulence to be used in a fully 3D global ocean model.

Depth-invariant meridional flux
The QG framework readily provides a prediction for the vertical structure of the meridional buoyancy flux. According to equation (3.2), the potential vorticity (3.3) is a material invariant of the undamped system. If the initial condition corresponds to vanishing PV, the PV will remain zero at subsequent times and the meridional PV flux vanishes at any depth. Within the QG framework this PV flux reads (Vallis 2006): which indicates that the meridional buoyancy flux is depth-independent. This QG prediction turns out to be in very good agreement with the numerical data, as illustrated by the meridional buoyancy flux profiles in Figure 7: the low-drag profile in the topleft panel of Figure 7a is depth-independent to very good accuracy. The moderate-drag profile in the bottom-left panel of Figure 7a is also depth-independent in the interior, with slight departures restricted to a boundary-layer near the surface. Finally, solutions to the QG system -equations (3.12-3.13) together with inversion relations (B 17-B 18)are characterized by depth-independent meridional flux profiles by construction.

Vertical flux and bulk stratification
Having determined the vertical structure of the meridional buoyancy flux, we turn to the vertical buoyancy flux using relation (4.2). We have seen that the overall contribution from dissipative terms on the right-hand side of that relation is sometimes not entirely negligible, with r ≃ 0.85 instead of r = 1, but how is that contribution distributed with depth?
We answer that question in Figure 7, where we compare the two flux terms on the left-hand side of (4.2). Within the central half of the domain, these two terms balance each other to a good accuracy: the eddy transport current is directed along isopycnals in the bulk of the domain. However, for several runs this balance becomes unsatisfactory as we move closer to the top and bottom boundaries, where the vertical flux needs to be compatible with the boundary conditions: w vanishes at the top free surface, while it matches the pumping velocity near the bottom boundary. In other words, the diffusive terms in (4.2) are not entirely negligible near the boundaries. In particular, we expect a strong contribution from the horizontal diffusive term in (4.2) at the surface, where surface quasi-geostrophic (SQG) dynamics predicts a forward cascade of buoyancy (towards small horizontal scales) Lapeyre 2017). The depthindependent meridional buoyancy flux in equation (4.2) is thus balanced by the vertical buoyancy flux term in the bulk of the domain and by the diffusive terms at the boundaries. The thickness of the surface-influenced region decreases with the diffusivities within the SQG framework. However, when the diffusivities reach very low values, sharp fronts develop, the latter being better described within the semi-geostrophic framework. Close inspection of our numerical runs indicates that the surface-influenced region is small only if the diffusivities are very weak, with large Reynolds (and Péclet) numbers, in a regime with weak isopycnal slopes and Ro ≪ 1. Figure 7 provides profiles from a numerical run in such a QG regime with extremely weak diffusivities. Both the meridional and the vertical buoyancy flux are depth-independent to a good approximation. The scaled buoyancy fluxes in figure 7 are approximately equal in the bulk of the domain, which shows that the buoyancy flux is directed along isopycnals. We also provide the profile of the mixing term ǫ b = E b;z (∂ z b) 2 +E b;⊥ |∇ ⊥ b| 2 for a series of runs where we lower both the vertical and horizontal diffusivities for otherwise constant parameters. We observe that there is significant diapycnal mixing near the top boundary, over a depth that decreases (albeit slowly) as the diffusivities are lowered. As the latter become asymptotically small, we expect the diapycnal mixing to be confined to the immediate vicinity of the upper boundary, the eddy-induced buoyancy current being then along isopycnals in most of the fluid domain.
The vertical structure of the emergent stratification is deduced from the horizontaland time-average of the buoyancy equation (2.14): where we have neglected the very weak contribution from v (for the two runs displayed in Figure 7, this contribution is smaller than the contribution from wb by more than two orders of magnitude). The emergent stratification profile is thus readily given by that of the vertical buoyancy flux. In the bulk of the domain, the emergent stratification is thus uniform in the low-diffusivity strongly-QG regime where wb is depth-independent.

Gradual barotropization
It is worth stressing the fact that the depth-invariance of the fluxes arises even when the turbulent kinetic energy profile (TKE) retains some significant depth dependence, as dictated by the QG result and illustrated by the profiles for a κ * = O(1) run in the bottom panels in Figure 7a (κ * ≃ 0.6). For lower drag, κ * = O(10 −1 ), the flow barotropizes with a TKE profile that gradually becomes depth-invariant, as seen in the upper panels in Figure 7a (κ * ≃ 0.14). The approach to the asymptotic low-drag predominantly barotropic state is thus rather slow, and this slow approach can be understood from QG vortex gas dynamics: multiply the uniform bulk QG PV (3.3) withb = ∂ zp , where the tilde denotes departure from the horizontal and temporal mean, before averaging horizontally and over time. After a few integrations by parts one is led to: The vortex gas scaling regime predicts that v 2 is asymptotically greater thanb 2 at low drag. The second term on the left-hand side of (5.3) is thus negligible at low drag and one concludes that the flow is barotropic. However, the approach to this asymptotically barotropic state is rather slow as κ * decreases. Indeed, the vortex gas scalings state that the (appropriately non-dimensionalized) buoyancy variance scales as D * , while the horizontal kinetic energy scales as D * log D * (see Gallet & Ferrari (2020) for the derivation of these scaling estimates). While D * increases very rapidly with decreasing drag -as exp(c 2 /κ * ), see equation (3.16) -the ratio of the first over the second term on the left-hand side of (5.3) increases only as log D * ∼ 1/κ * , hence the slow approach to the asymptotic barotropic state. The vortex gas scaling theory performs surprisingly well even in the absence of complete barotropization. To some extent, this can be understood from 2LQG dynamics: the vortex gas theory is based on the idea that a barotropic vortex gas distorts the background buoyancy gradient, inducing buoyancy fluctuations at the inter-vortex scale. Within 2LQG, this regime corresponds to the large-scale-dynamics approximation given by equations (32) and (33) in Thompson & Young (2006) (see also Salmon (1980) and Larichev & Held (1995)), where the buoyancy field is stirred by the barotropic flow only. An interesting point is that complete barotropization is not required for this large-scaledynamics approximation to hold. Instead, only scale separation between the inter-vortex distance and the vortex core radius is necessary.
Of course, partial barotropization also affects the main power integral entering the theory, that is, the balance between the rate of release of available potential energy -proportional to the dimensionless meridional buoyancy flux D * -and the frictional dissipation term. Assuming that the latter corresponds mainly to dissipation by the barotropic flow, the rate of frictional dissipation scales as κ * D * log D * (where the logarithmic correction stems from the fast velocities arising in the vicinity of the vortex cores, see Gallet & Ferrari (2020)). Balancing this estimate for the frictional dissipation rate with the rate D * of release of APE leads to the scaling prediction (3.16) for D * .
The assumption that frictional dissipation is associated with the barotropic flow is questionable when barotropization is only partial. In the present 3D model, the frictional dissipation rate is given by the product of the effective friction coefficient with the squared bottom velocity. According to the relation (5.3), the bottom squared velocity may depart from the squared barotropic velocity by a correction proportional to b 2 (0), the latter scaling as D * in the vortex gas theory. Including this correction leads to the frictional dissipation rate being proportional to κ * D * [log(D * ) + const.] = κ * D * log(D * /const.). Balancing this last expression for the frictional dissipation rate with the dimensionless rate of release of APE, D * , leads again to the scaling prediction (3.16) for the meridional diffusivity. We conclude that partial barotropization does not affect the scaling prediction for the meridional buoyancy flux.
In summary, it seems that the vortex gas theory only requires that the gas be dilute: the mixing length must be much greater than λ, the ratio of the two scaling as √ D * . Using the scaling expression (3.16) for D * , we conclude that the vortex gas scaling predictions hold provided exp(const./κ * ) ≫ 1, whereas strong barotropization of the flow requires the much more stringent criterion log D * ∼ 1/κ * ≫ 1.

Conclusion
Through a suite of numerical runs of the Boussinesq Eady problem, we have shown that the strongly stratified rapidly rotating regime is governed by quasi-geostrophic dynamics. One consequence is that the dependence of the transport properties of the turbulent flow on all control parameters is encoded in a scaling relation between the dimensionless diffusivity and the dimensionless drag. That scaling relation is correctly captured by the vortex gas scaling theory, see equation (3.16), initially put forward in the context of the 2LQG model and validated against numerical simulations with both linear and quadratic drag. In the context of the Eady problem, the two adjustable parameters arising in that scaling relation take slightly different values than in the 2LQG problem. We stress the fact that adjusting these two parameters leads to a quantitative theory for the magnitude and vertical structure of the meridional buoyancy flux, the vertical buoyancy flux and the emergent stratification, in terms of the background shear, the Coriolis parameter, the background stratification, the bottom friction and the diffusivities.
Indeed, the estimate (3.16) for the magnitude of the meridional buoyancy flux readily translates into a quantitative prediction for the vertical buoyancy flux through the standard along-isopycnal-transport argument. This argument is validated by comparing the numerical value of r with the theoretical value r = 1: we obtain numerically r ≃ 0.85 over the suite of simulations, with r → 1 as we enter the very-low-diffusivity QG regime.
The vertical buoyancy flux in turn leads to a quantitative prediction for the magnitude of the emergent stratification, set by a balance between eddy-induced vertical buoyancy transport and diffusion. This prediction for the emergent stratification allows us to determine the region of parameter space where the emergent stratification is negligible as compared to the imposed background stratification, and the region of parameter space where the sole emergent stratification is sufficient to induce QG dynamics in the absence of a background stratification.
Beyond scaling analysis, the 3D Eady model allows one to address the vertical structure of the eddy-induced transport. QG theory predicts that the meridional buoyancy flux is depth-invariant. This prediction is very specific to the Eady model, however, as it stems from the absence of interior PV gradients. One expects some vertical structure for the meridional buoyancy flux in a more realistic model with interior PV gradients. The prediction of a depth-invariant meridional flux is well satisfied by the numerical data (more so in the low-drag regime). The along-isopycnal-transport argument then implies that the vertical buoyancy flux is also depth-independent. The latter prediction is indeed validated by the strongly QG low-diffusivity data, although some depth-dependence remains within small boundary layers near the top and bottom boundaries (the thickness of which decreases, albeit rather slowly, with the diffusivities). A depth-independent vertical buoyancy flux in turn induces a uniform emergent stratification. The turbulent kinetic energy profile retains some depth dependence for moderate drag and gradually becomes more depth-invariant as the drag is decreased, in line with a barotropization of the flow. This slow barotropization of the flow as the drag coefficient decreases can be explained within the vortex gas scaling theory, and we have provided arguments indicating that some baroclinicity does not invalidate the theory. In particular, the theory could have skill to describe the transport induced by the mid-ocean eddies considered in Arbic & Flierl (2004) based on current-meter data (Wunsch (1997) and references therein), which retain some significant baroclinicity. For instance, the bottom panels in Figure 7a illustrate an equilibrated flow with comparable baroclinic and barotropic turbulent kinetic energies. Such near equipartition of barotropic and baroclinic energies is consistent with ACC dynamics. The dimensionless drag coefficient is κ * ≃ 0.6 for this run. Assuming the ocean-like parameter values Ro = 0.25 and N/f = 30, this value for κ * translates into a dimensionless linear drag coefficient κ eff = 0.005. Using a Coriolis parameter f ≃ 10 −4 rad.s −1 the corresponding dimensional magnitude of the linear drag coefficient is (25 days) −1 . The order of magnitude is correct, this value being greater by a factor of four than the typical estimates reported in Arbic & Flierl (2004). Vorticity snapshots associated with this run are provided in Figure 2d,e. We observe that the vortex gas is moderately dilute, with an inter-vortex distance comparable to the core radius. The inter-vortex distance increases to 3 to 5 core radii for the lower-drag run illustrated in Figure 2b,c. These values for the inter-vortex distance are within the range of typical oceanic observations. Beyond the 2LQG approach of Arbic & Flierl (2004), a more quantitative comparison with ACC dynamics and mid-ocean eddies in general would require the inclusion of the vertical structure of the current profile and density stratification into the 3D model, together with quadratic bottom drag.
Of course, whether the equilibrated Eady flow is governed by QG vortex gas dynamics depends on the values of the control parameters, and the runs in Table 1 represent only a subset of our entire suite of numerical runs. Satisfactory collapse onto the QG vortex-gas master curve in Figure 4 arises provided (i) the diffusivities are sufficiently low and (ii) the assumptions of QG dynamics are well satisfied. Specifically, for a run to be retained in Table 1 we demand that the Reynolds numbers be large enough, and that the isopycnal slope and vortex core vorticity be low enough. Regarding the low  Figure 4 for the entire dataset. The blue stars denote the runs retained throughout the manuscript and reported in Table 1. They are characterized by large enough Reynolds numbers, Rev > 4000 and Re h > 350, together with low enough core vorticity to ensure QG dynamics, vb /(Roλ 2 ) < 0.5. The runs that violate any of these inequalities are represented with black crosses and reported in Table 2. They are characterized by a lower eddy diffusivity as compared to the QG vortex-gas master curve in Figure 4. diffusivities, we used a threshold value for the vertical Reynolds (and Péclet) number Re v = Ro/E z of the base shear flow, Re v > 4000, together with a threshold value for the emergent horizontal Reynolds (and Péclet) number Re h = v 2 λ/E ⊥ > 350. Regarding the assumptions of QG dynamics, we find that the isopycnal slope is small in all the numerical runs. However, the bulk vertical vorticity ζ reaches large values inside the vortex cores, especially in the low-drag regime, see Figure 2b. Within the vortex gas theory, the circulation of the vortices scales like the meridional eddy diffusivity vb /Ro, while the vortex core radius scales like the Rossby deformation radius λ. Demanding that the core vorticity be at most comparable to f then amounts to satisfying the inequality vb /(Roλ 2 ) < 0.5, where the threshold value 0.5 has been chosen somewhat arbitrarily. The runs in Table 1 satisfy the three inequalities above. The emergent Rossby number, measured as the rms vertical vorticity at mid-depth divided by f , is at most of the order of 0.1 for these runs. The remaining runs are reported in Table 2. They violate at least one of the three inequalities above, which results in an eddy diffusivity that lies below the QG vortex-gas master curve, as illustrated in Figure 8. An extreme example of breakdown of the QG model when the assumptions of QG dynamics are not satisfied is reported in Molemaker et al. (2010). They perform numerical simulations without bottom friction. In the absence of a large-scale energy sink for the QG inverse cascade, the flow intensifies until the emergent Rossby number is of the order of unity. The QG expansion then breaks down and, consistently, their Boussinesq simulations greatly depart from QG dynamics, with a forward energy cascade that halts the intensification of the flow.
In summary, our study demonstrates how the vortex-gas scaling theory can be extended to describe the structure and magnitude of eddy-induced transport by fully 3D baroclinic turbulence. The skill of the vortex-gas scaling theory for both the 2LQG model and the Eady model shows that the theoretical predictions hold regardless of whether the flow is driven by boundary dynamics or interior PV gradients. The model is strongly idealized, however, and should be complexified if the resulting parameterization is to describe real atmospheres or oceans. In the atmospheric context, one may want to include β and address the eddy-induced transport in the equilibrated Charney model. The simple advective-diffusive balance that sets the emergent vertical stratification could be modified as well, with the goal of better capturing the feedback of baroclinic turbulence on the vertical structure of the atmosphere (Held & Suarez 1994).
β-plane dynamics induced by planetary curvature or sloping topography are also important ingredients of ocean dynamics. In the context of the ACC, an accurate description of mesoscale transport probably requires a base state with depth-dependent zonal shear and an inversion of the meridional potential vorticity gradient at some finite depth (Charney & Stern 1962;Pedlosky 1979;Smith & Marshall 2009;Abernathey et al. 2010). The vortex gas approach has been extended to the β-plane within the 2LQG framework, but the skill of the resulting theory remains to be assessed for such a vertically-structured fully 3D model. while κ 2 is given by: An alternative approach to determining κ eff consists in computing the friction term induced by Ekman pumping associated with a bulk flow that varies slowly with x and y. Denoting the bulk vertical vorticity as ζ, the pumping velocity just above the partial Ekman spiral reads w| 0 + = κ 3 ζ| 0 + . The goal here is to show that κ 3 is equal to κ eff above.
For simplicity, consider a purely barotropic flow U (y) that varies slowly in y. The bulk flow has vertical vorticity ζ = −U ′ (y). The total flow, including the partial Ekman spiral, is: and the incompressibility condition imposes a vertical velocity such that: Integrating over the Ekman layer thickness yields the pumping velocity: and we conclude that: One can check that the two approaches yield the same effective friction coefficient κ eff , that is κ eff = κ 3 = κ 1 + κ 2 .
where k = |k|. The solution to this equation whose vertical derivative matches {∂ z P | 0 } k and {∂ z P | 1 } k at the boundaries is: P k (z, T ) = −{∂ z P | 0 } k cosh[k(z − 1)] + {∂ z P | 1 } k cosh(kz) k sinh k , (B 16) Evaluating this quantity at the boundaries z = 0 and z = 1 yields the two inversion relations providing {P | 0 } k and {P | 1 } k in terms of {∂ z P | 0 } k and {∂ z P | 1 } k : A time step of the QG Eady problem consists in (i) marching the 2D equations (3.12) and (3.13) for one timestep, before (ii) inferring the new fields P | 0 + and P | 1 using the inversion relations (B 17) and (B 18).
B.3. Symmetry argument for wu Consider a solutionp = P(x, y, z, t), q = Q(x, y, z, t) to the QG system: the functions P and Q satisfy equations (3.2-3.3) with the boundary conditions (3.4) and (3.8). The vertical velocity associated with this solution is inferred from equation (B 7), which reads in terms of the original variables: −λ 2 w = ∂ tzp − Ro ∂ xp + Ro z∂ xzp + J(p, ∂ zp ) .

(B 19)
Introduce the following transformed solution: p(x,ŷ, z, t) = −P(x, y = −ŷ, z, t) ,q(x,ŷ, z, t) = −Q(x, y = −ŷ, z, t) . (B 20) The fieldsp(x,ŷ, z, t) andq(x,ŷ, z, t) also satisfy the QG system. To check this, the operator ∂ y should be understood as a derivative with respect to the second variable, and thus replaced by ∂ŷ in equations (3.2-3.3) and boundary conditions (3.4) and (3.8) wherẽ p and q have been replaced byp andq, respectively. There is a direct correspondence between the fields associated with the original solution and the fields associated with the transformed one (the latter being evaluated atŷ = −y). Denoting the former with the standard notations and the latter with hat variables, we obtain: where the hat variables are evaluated atŷ = −y, while the variables associated with the original solution are evaluated at y. The last relation has been obtained by substitution into equation (B 19) for the vertical velocity. This correspondence indicates, for instance, that the meridional and vertical buoyancy fluxes are equal for the original and for the transformed solutions:vb = vb andŵb = wb. However, the vertical flux of zonal momentum is opposite between the two solutions:ŵû = −wu. Because the original solution and the transformed ones are equally probable, we conclude that wu = 0 within the QG framework, while vb = 0 and wb = 0.