# Vertical structure of buoyancy transport by ocean baroclinic turbulence

Julie Meunier, Benjamin Miquel, Basile Gallet

► To cite this version:

Julie Meunier, Benjamin Miquel, Basile Gallet. Vertical structure of buoyancy transport by ocean baroclinic turbulence. Geophysical Research Letters, 2023, 50 (17), 10.1029/2023GL103948. hal-04301208

HAL Id: hal-04301208

https://hal.science/hal-04301208v1

Submitted on 22 Nov 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2023GL103948

**Key Points:**

* We derive a prediction for the depth-dependence of the buoyancy flux associated with ocean baroclinic turbulence in an idealized setup

* The prediction is validated quantitatively by simulations of an idealized patch of ocean meant to resemble Southern Ocean conditions

* The prediction can be readily implemented into global models as a vertical profile for the Gent-McWilliams coefficient

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

B. Gallet,
basile.gallet@cea.fr

**Citation:**

Meunier, J., Miquel, B., & Gallet, B. (2023). Vertical structure of buoyancy transport by ocean baroclinic turbulence. *Geophysical Research Letters*, *50*, e2023GL103948. https://doi.org/10.1029/2023GL103948

Received 3 APR 2023
Accepted 31 JUL 2023

© 2023. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

## Vertical Structure of Buoyancy Transport by Ocean Baroclinic Turbulence

Julie Meunier¹, Benjamin Miquel² ID, and Basile Gallet¹ ID

¹Service de Physique de l'Etat Condensé, CNRS, CEA, Université Paris-Saclay, Gif-sur-Yvette, France, ²LMFA, UMR5509, Ecole Centrale de Lyon, CNRS, Univ Lyon, INSA Lyon, Université Claude Bernard Lyon 1, Ecully, France

**Abstract** Ocean mesoscale eddies enhance meridional buoyancy transport, notably in the Antarctic Circumpolar Current where they contribute to setting the deep stratification of the neighboring ocean basins. The much-needed parameterization of this buoyancy transport in global climate models requires a theory for the overall flux, but also for its vertical structure inside the fluid column. Based on the quasi-geostrophic dynamics of an idealized patch of ocean hosting an arbitrary vertically sheared zonal flow, we provide a quantitative prediction for the vertical structure of the buoyancy flux without adjustable parameters. The prediction agrees quantitatively with meridional flux profiles obtained through numerical simulations of an idealized patch of ocean with realistic parameter values. This work empowers modelers with an explicit and physically based expression for the vertical profile of buoyancy transport by ocean baroclinic turbulence, as opposed to the common practice of using arbitrary prescriptions for the depth-dependence of the transport coefficients.

**Plain Language Summary** Ocean mesoscale vortices are turbulent structures tens of kilometers wide that play a central role in transporting tracers such as heat, salt, and carbon. In the Southern Ocean, the associated buoyancy transport crucially sets the deep stratification of neighboring ocean basins. Because mesoscale vortices are not resolved by most state-of-the-art climate models, modelers resort to rather crude parameterizations where—in the absence of a better theory—the transport properties of the eddies are often assumed to be depth-invariant in the ocean interior. In this contribution we derive a quantitative and parameter-free prediction for the vertical structure of the turbulent buoyancy flux, which can be readily implemented in global models at little computational cost.

## 1. Introduction

The baroclinic instability of large-scale ocean currents generates mesoscale eddies that strongly enhance heat and tracer transport. In the Antarctic Circumpolar Current (ACC), the resulting turbulent buoyancy transport contributes to setting the slope of the Southern Ocean density surfaces and therefore the deep stratification of the neighboring ocean basins (Nikurashin & Vallis, 2011, 2012; Wolfe & Cessi, 2010). Mesoscale eddies have a core size comparable to the Rossby deformation radius, a length scale of the order of 60 km at midlatitudes and 15 km in the Southern Ocean, smaller than the coarse resolution of most global climate models. Parameterizing the transport induced by mesoscale eddies in such global models is thus crucial to obtain realistic ocean states that quantitatively reproduce the sloping density surfaces of the Southern ocean and the deep stratification of ocean basins. Physically-based parameterizations are inferred from the study of an isolated patch of ocean, where baroclinic turbulence has homogeneous statistics in the horizontal directions. The parameterization problem then consists in determining the scaling behavior of the overall diffusivity in terms of the various control parameters (shear flow magnitude, background stratification, bottom friction coefficient, etc.) but also the vertical structure of the various fluxes within the water column. Far more studies have addressed the former task (Arbic & Flierl, 2004a, 2004b; Arbic & Scott, 2007; Chang & Held, 2019; Gallet & Ferrari, 2020, 2021; Held & Larichev, 1996; Larichev & Held, 1995; Phillips, 1954; Salmon, 1978, 1980; Thompson & Young, 2006, 2007) than the latter (Stanley et al., 2020; Yankovsky et al., 2022; Zhang & Wolfe, 2022) in the ocean context. In the absence of a better theory many global models assume that the transport coefficients are depth-invariant in the ocean interior (see, e.g., Griffies et al., 2005), while other models consider surface-intensified coefficients with arbitrary prescriptions for their vertical structure (such as, e.g., assuming that the coefficients are proportional to the local squared buoyancy frequency, Danabasoglu & Marshall, 2007; Ferreira et al., 2005; Gent, 2011). The latter assumption of surface-intensified transport coefficients is at odds with idealized eddy-resolving channel
<figure><img src="image_3.png" /><figcaption>**Figure 1.** An idealized patch of ocean. A layer of fluid is subject to global rotation at a rate that varies linearly with the meridional coordinate y. The fluid is density stratified with an arbitrary profile N(z) for the buoyancy frequency. The background zonal shear flow has an arbitrary profile U(z). This flow coexists with a background meridional buoyancy gradient. Friction damps kinetic energy on the ocean floor.</figcaption></figure>

simulations, which point to a bottom-enhanced buoyancy transport coefficient instead (Abernathey et al., 2013) (the so-called Gent-McWilliams coefficient, see below).

To improve upon this unsatisfactory state of the art, in this Letter we derive a parameter-free prediction for the vertical structure of the turbulent buoyancy flux within the water column. We consider an idealized patch of ocean with arbitrary background zonal shear flow and stratification, and $\beta \neq 0$, see Figure 1. Water occupies a volume $(x, y, z) \in [0, L]^2 \times [-H, 0]$ with a stress-free boundary at $z=0$ and a linear-friction boundary condition at $z=-H$, in a frame rotating around the vertical axis with a local Coriolis parameter $f_0 + \beta y$, where y denotes the meridional (North-South) coordinate. The fluid layer is density-stratified with an arbitrary profile N(z) for the buoyancy frequency, and we restrict attention to a single stratifying agent. We focus on the quasi-geostrophic (QG) regime arising for fast rotation and strong stratification (Salmon, 1998; Vallis, 2017; Venaille et al., 2011). The base flow consists of an arbitrary zonal velocity profile U(z) in thermal wind balance with a z-dependent meridional buoyancy gradient $\partial_z B = -f_0 U'(z)$, where the prime symbol denotes

a vertical derivative. We consider arbitrary departures from this base state with periodic boundary conditions in the horizontal directions. We denote as $p(x, y, z, t)$ the departure from the base pressure field, with $u = -p_y$ the departure zonal velocity, $v = p_x$ the departure meridional velocity, $b = f_0 p_z$ the departure buoyancy and $w$ the subdominant (ageostrophic) vertical velocity. Non-dimensionalizing time and space using $|f_0|^{-1}$ and $H$, the dimensionless base flow is written as $U/|f_0|H = Ro U(z)$, where $Ro = |U(0)|/f_0 H$ is the Rossby number associated with the surface speed of the base flow and $U(z)$ denotes the base-flow profile normalized at the surface ($|U'(0)| = 1$). For brevity we use the same symbols for the dimensionless variables, except for the notation $\tilde{\beta} = \beta H/f_0$.

Consider a tracer $\tau$ stirred by the 3D flow and subject to horizontally uniform gradients (at lowest order in $Ro$) $G_y^{(\tau)}(z)$ and $G_z^{(\tau)}(z)$ in the meridional and vertical directions, respectively. The QG evolution equation for $\tau$ reads:

$$ \partial_t \tau + Ro U(z) \tau_x + J(p, \tau) = -p_x G_y^{(\tau)}(z) - w G_z^{(\tau)}(z) + D_{\tau}, \quad (1) $$

where the Jacobian is $J(g, h) = g_x h_y - g_y h_x$ and $D_{\tau}$ denotes small-scale diffusion.

Denoting with an overbar $\tau$ a time average together with a horizontal area average, the eddy-induced meridional and vertical fluxes of $\tau$ are related to the background gradients by a Gent-McWilliams/Redi (GM/R) diffusion tensor (Gent, 2011; Gent & McWilliams, 1990; Griffies, 1998; McDougall & McIntosh, 2001; Redi, 1982):

$$ \left( \frac{\overline{\nu \tau}}{\overline{w \tau}} \right) = \begin{bmatrix} -K_R & (K_{GM} - K_R)S & \left[ G_y^{(\tau)} \right] \\ -(K_{GM} + K_R)S & -K_R S^2 & \left[ G_z^{(\tau)} \right] \end{bmatrix} \quad (2) $$

where the Redi diffusivity $K_R(z)$ encodes diffusion along the mean isopycnal direction, the GM coefficient $K_{GM}(z)$ encodes the advective (or skew-diffusive) transport, and we denote the isopycnal slope of the base state as $S(z) = Ro U'(z)/N^2$. While (2) is often introduced based on physical intuition and educated guesses, we have recently proposed a direct derivation of this diffusion tensor from the QG dynamics of the present system (Meunier et al., 2023). For completeness we briefly recall a few results from this recent study.

The QG potential vorticity (QGPV) $q = \Delta_{\perp} p + \partial_z [p_z/N^2(z)]$ is governed by Equation 1 with $\tau = q$, $G_y^{(q)} = 0$ and $G_y^{(q)} = \tilde{\beta} - S'(z)$, while buoyancy is governed by Equation 1 with $\tau = b$, $G_z^{(b)} = N^2$ and $G_y^{(b)} = -Ro U'(z)$. Substitution of these background gradients into the flux-gradient relation (2) indicates that $K_{GM}(z)$ and $K_R(z)$ can alternatively be thought of as the effective diffusivities associated with the meridional transport of $b$ and $q$, respectively:

$$ K_{GM} = - \frac{\overline{\nu b}}{G_y^{(b)}} = \frac{\overline{\nu b}}{Ro U'(z)}, \quad K_R = - \frac{\overline{\nu q}}{G_y^{(q)}} = \frac{\overline{\nu q}}{S'(z) - \tilde{\beta}}. \quad (3) $$

Because the vertical velocity vanishes at the surface, the governing equations for $q$ and $b$ admit the same limiting form as one approaches the top boundary. Both tracers are advected by the surface horizontal flow, fluctuations being induced by distortions of a horizontally homogeneous background meridional gradient. The associated
meridional diffusivity is thus equal for b and q at the surface: $\overline{vb}(0)/G_y^{(b)}(0) = \overline{vq}(0)/G_y^{(q)}(0)$. Provided the friction coefficient is small, the same holds near the bottom boundary, at a depth $z = -1^+$ located just above the bottom Ekman layer: $\overline{vb}(-1^+)/G_y^{(b)}(-1^+) \approx \overline{vq}(-1^+)/G_y^{(q)}(-1^+)$ (while the friction-induced vertical pumping velocity is crucial for damping kinetic energy through the stretching of planetary vorticity, it has a negligible direct contribution to buoyancy transport for low drag coefficient). Using Equation 3 we recast these equalities as:

$$K_{GM}(0) = K_{R}(0), \quad K_{GM}(-1^{+}) \approx K_{R}(-1^{+}). \qquad (4)$$

The two equalities in Equation 4 are illustrated numerically in Meunier et al. (2023). An additional constraint on $K_{GM}$ and $K_R$ is obtained by substituting the definition of $q$ into the meridional QGPV flux $\overline{vq}$. After a few integrations by parts one obtains the Taylor-Bretherton relation $\overline{vq} = d(\overline{vb}/N^2)/dz$ (Bretherton, 1966; Dritschel & McIntyre, 2008; Smith & Marshall, 2009; Taylor, 1915; Young, 2012), and expressing the meridional fluxes using Equation 3:

$$K_R(S' - \tilde{\beta}) = \frac{d}{dz}(K_{GM} S). \qquad (5)$$

In the following we show that the constraints (4–5) allow for a perturbative derivation of the vertical structure of the eddy-induced buoyancy flux within the water column in two situations of interest.

## 2. Case I: The Impact of Weak $\beta$ on Eady Turbulence

The QG Eady model corresponds to depth-independent stratification $N^2$ and shear $\mathcal{U}'$ (linear zonal velocity profile, $\mathcal{U}(z) = z + 1$), together with $\beta = 0$. As discussed in Gallet et al. (2022), there is no background PV gradient in this setup and therefore a solution can be obtained by assuming $q = 0$ in the bulk of the domain. The meridional QGPV flux then vanishes, and from relation (5) we conclude that the meridional buoyancy flux, and thus $K_{GM}$, are independent of $z$.

As established in Meunier et al. (2023), $K_R(z)$ is given by the Taylor-Kubo eddy diffusivity coefficient associated with the horizontal geostrophic flow. That is, at every depth $z$ the coefficient $K_R(z)$ is given by the integral of the Lagrangian correlation function of the horizontal geostrophic flow. Because in the low-drag limit the Eady flow barotropizes, we expect the horizontal geostrophic flow to be depth-invariant, which leads to $K_R$ being independent of $z$. Using the boundary relation (4) we conclude that the GM and Redi coefficients are depth-invariant and equal to one another. The low-drag Eady model thus represents one limiting situation for which the depth-invariance and equality of the GM and Redi coefficients can be established. We stress the fact that the equality of the GM and Redi coefficients has been established based on the properties of the low-drag equilibrated state, namely barotropization, and the theory presented below is really a theory for such a low-drag equilibrated—or “turbulent”—state. By contrast, the theory would not hold to predict the vertical structure of an eigenmode obtained using linear stability analysis, whose transport properties typically display strong depth-dependence (see Supporting Information S1).

Consider now the impact of a weak planetary vorticity gradient on Eady turbulence, that is, a Charney model with weak $\beta$. Within QG, the impact of $\beta$ is characterized by the product of $\beta$ with the squared deformation radius over the typical velocity of the background shear flow (Chang & Held, 2021; Charney, 1947; Gallet & Ferrari, 2021; Thompson & Young, 2007). We thus define the $(z$-invariant) parameter $\beta_* = \tilde{\beta}N^2/R_0$. In the perturbative regime $\beta_* \ll 1$ the correction to the $z$-invariant $\beta_* = 0$ situation is small, and a standard expansion leads to $K_R(z) = K_{GM}(z)[1 + \mathcal{O}(\beta_*)]$, where the $\mathcal{O}(\beta_*)$ correction vanishes both at the top and at the bottom boundary in the low-drag weakly diffusive regime, see Equation 4. Substitution into Equation 5 yields $K'_{GM}(z) = -\beta_* K_{GM}(z) + \mathcal{O}(\beta_*^2)$ and, neglecting the $\mathcal{O}(\beta_*^2)$ correction, $K_{GM}(z) = \text{const.} \times e^{-\beta_* z}$. Substituting this expression for $K_{GM}(z)$ into Equation 3 and denoting the overall buoyancy flux as $\langle \overline{vb}(z) \rangle = \int_{-1}^{0} \overline{vb}(\bar{z})d\bar{z}$, we obtain a parameter-free prediction for the vertical structure $\overline{vb}(z)/\langle \overline{vb}(z) \rangle$ of the meridional buoyancy flux:

$$\frac{\overline{vb}(z)}{\langle \overline{vb}(z) \rangle} = \frac{\beta_*}{e^{\beta_*} - 1} e^{-\beta_* z}. \qquad (6)$$

To test this perturbative prediction, we have performed numerical simulations of this setup in the QG regime with periodic boundary conditions in the horizontal directions. As detailed in Supporting Information S1 (see also
<figure><img src="image_4.png" /><figcaption>**Figure 2.** An illustrative example: the Charney model with weak $\beta$. (a) Snapshot of the departure buoyancy field $b$ for $\beta_*=0.5$ (large values in red, low values in blue). (b) Vertical structure of the meridional buoyancy flux in the equilibrated state for increasing $\beta_*$ (solid line: DNS, dashed line: perturbative prediction (6)). The agreement with the prediction is excellent in the perturbative regime $\beta_* \ll 1$ and deteriorates somewhat as $\beta_*$ reaches $\mathcal{O}(1)$ values. The perturbative prediction performs always better than the common practice of parameterizing turbulent transport using a depth-invariant $K_{GM}$, which corresponds to depth-invariant $\overline{vb}$ for the present setup. The prediction (6) also performs better than using the meridional flux associated with the most unstable eigenmode, computed perturbatively for weak $\beta_*$, and represented for $\beta_* = 0.1$ as a green dash-dotted line.</figcaption></figure>

Meunier et al. (2023)), our numerical approach consists in time-stepping a set of primitive-like equations with tailored $\beta$ terms that are compatible with the horizontal periodic boundary conditions. Importantly, these tailored terms reduce to the standard $\beta$ terms in the QG limit. Because we focus on parameter values that are strongly QG, this approach is equivalent to (but more convenient than) directly solving the QG system. In Figure 2 we plot the vertical structure of the meridional buoyancy flux, $\overline{vb}(z)/\langle vb \rangle$, for increasing $\beta_*$. The numerical profiles are in excellent agreement with the parameter-free prediction (6) for low $\beta_*$. As expected, the perturbative prediction deteriorates somewhat as $\beta_*$ increases up to $\beta_* = 1$. In the next section we show that the perturbative regime accurately captures the typical oceanic situation, characterized by surface-intensified baroclinic turbulence.

### 3. Case II: Surface-Intensified Shear and Stratification

The perturbative approach developed in the preceding section is based on the small value of the background meridional PV gradient: when low-drag baroclinic turbulence is subjected to a weak meridional PV gradient, the profile of $K_{GM}(z)$ can be inferred by inserting $K_R(z) \approx K_{GM}(z)$ into the Taylor-Bretherton relation (5).

The meridional PV gradient associated with $\beta$ is modest in a typical oceanic setting, which may suggest that one can again use the approximate relation $K_R(z) \approx K_{GM}(z)$ in the bulk of the domain. However, the PV gradient associated with the $z$-dependent shear profile is much greater (see Figure 6 of Smith and Marshall (2009)). Fortunately, ocean baroclinic turbulence is surface-intensified and the largest shear-induced meridional PV gradient arises in the upper region of the fluid column, where the approximate equality $K_R \approx K_{GM}$ holds by virtue of the near-surface relation (4). Once again, one can thus substitute the approximate relation $K_R(z) \approx K_{GM}(z)$ into the Taylor-Bretherton relation (5) to compute the profile of $K_{GM}$, this time perturbatively in distance from the upper boundary. We conclude that a useful approximation to the vertical dependence of $K_{GM}$ should be obtained by substituting $K_R(z) \approx K_{GM}(z)$ into Equation 5 throughout the entire water column. As for Case I, one way to derail this procedure would be to have a surprisingly large meridional PV flux arise in the interior of the domain despite the very weak PV gradient. This typically happens for a nearly marginal eigenmode, but not for the present low-drag equilibrated—or “turbulent”—states. In that respect the theory below is really a theory for such equilibrated baroclinic turbulence.

After re-arranging, the substitution of $K_R(z) \approx K_{GM}(z)$ into Equation 5 leads to the following ODE for the vertical structure of the GM coefficient:

$$ \frac{d}{dz} \ln K_{GM} = - \frac{\tilde{\beta}}{S(z)}. \qquad (7) $$
This relation points to the crucial role of $\beta$ in setting the vertical structure of the eddy-induced buoyancy flux: according to Equation 7 the common assumption of a depth-invariant GM coefficient is valid for $\beta = 0$ only. For arbitrary $\beta$ Equation 7 can be integrated into:

$$K_{GM}(z) = \text{const.} \times \exp\left[-\int_{0}^{z} \frac{\tilde{\beta}}{S(\tilde{z})} d\tilde{z}\right]. \qquad (8)$$

We have obtained an explicit expression for the vertical structure of the GM coefficient in terms of the vertical profiles of background stratification and shear. Using Equation 3, the expression (8) can be recast into a parameter-free prediction for the vertical structure of the meridional buoyancy flux:

$$\frac{\overline{vb}(z)}{\langle vb \rangle} = \frac{\mathcal{U}'(z) \exp\left[-\int_{0}^{z} \frac{\beta}{S(z)} d\tilde{z}\right]}{\int_{-1}^{0} \mathcal{U}'(z) \exp\left[-\int_{0}^{z} \frac{\tilde{\beta}}{S(\tilde{z})} d\tilde{z}\right] dz} \qquad (9)$$

To test the prediction (9) we have performed numerical simulations of surface-intensified baroclinic turbulence with parameter values typical of the ACC. The dimensionless stratification profile is linear in $z$ and surface intensified, $N^2(z) = a_0 + a_1 z$, with constant coefficients $a_0$ and $a_1$. The shear flow has an exponential profile $\mathcal{U}(z) = s e^{z/\ell}$ with an e-folding scale $\ell$ (in units of $H$) and a sign prefactor $s = +1$ for an eastward flow and $s = -1$ for a westward one. We perform a base run with dimensionless QG parameter values similar to the situation addressed by Smith and Marshall (2009): dimensional magnitude $|f_0| = 1.23 \times 10^{-4}$ s$^{-1}$ for the Coriolis parameter and $\beta = 1.23 \times 10^{-11}$ m$^{-1}$ s$^{-1}$ for the planetary vorticity gradient (corresponding to a latitude of 57.5°S), depth of fluid equal to $H = 4,000$ m, eastward shear flow with surface speed $U(0) = 0.15$ m s$^{-1}$ and e-folding scale of 2,000 m. The dimensional buoyancy frequency ranges from $8.7 \times 10^{-4}$ s$^{-1}$ at the bottom to $2.46 \times 10^{-3}$ s$^{-1}$ at the surface, which corresponds to a Rossby deformation radius $\lambda \approx 19$ km based on the rough WKB estimate $\lambda/H = \int_{-1}^{0} N(z) dz/\pi$ (recalling that $N(z)$ is non-dimensionalized with $|f_0|$). In terms of dimensionless parameters, these values translate into a Rossby number $Ro = 0.3$, a vertical scale $\ell = 0.5$ for the shear flow, a dimensionless planetary vorticity gradient $\tilde{\beta} = 4.0 \times 10^{-4}$ and stratification coefficients $a_0 = 400$ and $a_1 = 350$. To ensure that the base numerical run indeed corresponds to the fully QG regime, we have used values for $Ro$ and $\tilde{\beta}$ that are smaller by a factor of 10 (i.e., we use $Ro = 0.03$ and $\tilde{\beta} = 4.0 \times 10^{-5}$), which leaves invariant the dissipation-free QG dynamics.

Together with this base run we have performed a run without $\beta$ and a run with $\beta > 0$ and a westward base flow. These additional runs are performed with slightly larger stratification ($a_0 = 800$ and $a_1 = 700$) using the inferred values $Ro = 0.3$ and $\tilde{\beta} = 4.0 \times 10^{-4}$ for the dimensionless buoyancy and planetary vorticity gradients. Finally, we have repeated similar runs using the larger value $\ell = 1$ for the e-folding scale of the shear (see Supporting Information S1 for the values of the other parameters).

In Figure 3 we provide snapshots of the buoyancy and velocity fields in the equilibrated state of the base run. As expected the turbulence is surface intensified and so is the meridional buoyancy flux $\overline{vb}$, provided in the upper-right panel of Figure 4. Substituting the linear profile for $N^2(z)$ and the exponential profile for $\mathcal{U}(z)$ into expression (9), the theoretical prediction for the vertical structure of the meridional buoyancy flux becomes:

$$\frac{\overline{vb}(z)}{\langle vb \rangle} = \frac{\exp\left\{\frac{z}{\ell} + \frac{s\tilde{\beta}\ell^2}{Ro}[a_0 + a_1(z + \ell)]e^{-\frac{z}{\ell}}\right\}}{\int_{-1}^{0} \exp\left\{\frac{z}{\ell} + \frac{s\tilde{\beta}\ell^2}{Ro}[a_0 + a_1(z + \ell)]e^{-\frac{z}{\ell}}\right\} dz} \qquad (10)$$

We compare this prediction to the numerically determined meridional flux profiles in Figure 4. The agreement is very good for both values of $\ell$, both with and without $\beta$, and for both eastward and westward flows. For $\beta = 0$ the theoretical prediction is that of a depth-invariant GM coefficient, and thus a meridional buoyancy flux that inherits the vertical structure $\mathcal{U}'(z)$ of the background shear. The good agreement with the numerical profiles validates this prediction and indicates that a depth-invariant GM coefficient is indeed an excellent parameterization when $\beta = 0$. For $\beta \neq 0$, however, the prediction (10) departs from the common practice of using a depth-invariant GM coefficient. In all cases the prediction (10) better captures the vertical structure of $\overline{vb}$, without adjustable parameters (see Figure 4). The difference between the two predictions—Equation 10 versus uniform $K_{GM}$—is modest
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Surface-intensified baroclinic turbulence from the base run, meant to resemble a patch of the Antarctic Circumpolar Current (positive values in red and negative values in blue, background profiles in the upper-right panel).</figcaption></figure>

for $\ell = 0.5$ and greater for $\ell = 1.0$. In particular, using a uniform $K_{GM}$ would lead to the same vertical structure for the buoyancy flux regardless of whether the base flow is directed eastward or westward. By contrast, the numerical data indicate that the vertical structure strongly depends on the direction of the base flow: for $\ell = 0.5$ the bottom-to-top meridional flux ratio, evaluated as $\overline{vb}(-0.95)/\overline{vb}(-0.05)$, is 12% for a westward flow and 24% for an eastward flow. For $\ell = 1$ this ratio is 15% for a westward flow and 73% for an eastward flow. To illustrate these differences, we represent in Figure 4 the uniform-$K_{GM}$ prediction as an orange dotted-line for comparison with the present prediction, demanding that the two predictions be equal at the top surface $z = 0$.

**4. Conclusion**

The predictions (8) and (9) are based on a perturbative approach that holds in the near-surface and near-bottom regions of the fluid column for arbitrary meridional potential vorticity gradient, and throughout the entire fluid column when the meridional potential vorticity gradient is weak. The present perturbative framework is useful for baroclinic turbulence in the ocean, where the shear flow and meridional PV gradient are boundary-intensified, with weaker PV flux in the interior. One can then combine the Taylor-Bretherton relation between the buoyancy and PV fluxes with the near-equality of the GM and Redi coefficients in the vicinity of the boundaries. This leads to a prediction for the vertical structure of the buoyancy flux that agrees well with the profiles extracted from direct numerical simulations, see Figure 4. It would be interesting to further investigate the range of validity of the predictions (8) and (9) beyond the present oceanographically relevant situations. For instance, a system with a vanishing meridional buoyancy gradient at the bottom ($G_y^{(b)}(-1) = 0$) may emphasize the role of bottom friction and disrupt the relation $K_{GM}(-1^+) \approx K_R(-1^+)$. More generally, while surprisingly successful the present perturbative approach should probably be used with caution whenever the exponential factor in Equation 8 varies by much more than a factor of two within the water column.

For eastward shear flows (positive shear) the right-hand side of Equation 7 is negative: $K_{GM}(z)$ is greater at depth according to both the theory and the numerics, even though the turbulence is surface-intensified. This prediction
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Vertical structure of the meridional buoyancy flux for an eastward flow, for a westward flow and for the case $\beta = 0$, using either $\ell = 0.5$ or $\ell = 1$. The solid line is the profile extracted from the numerical runs. The dashed line is the theoretical prediction (10). The orange dotted-line corresponds to the uniform-$K_{GM}$ model that matches the surface value of the full prediction (10) (see text for details). The prediction (10) reduces to a depth-invariant $K_{GM}$ when $\beta = 0$, which agrees accurately with the $\beta = 0$ numerical profiles. For $\beta \neq 0$ the prediction (10) departs from a uniform $K_{GM}$ and agrees well with the numerical profiles.</figcaption></figure>

is fully compatible with the $K_{GM}$-profile reported by Abernathey et al. (2013) and challenges models where the $K_{GM}$-profile is assumed to be proportional to the profile of $N^2(z)$ (Ferreira et al., 2005). The present results also seem to invalidate the idea that the vertical structure of the flux could be governed by a single baroclinic mode (Stanley et al., 2020). Indeed, the modal decomposition (Flierl, 1978) is the same for all panels of Figure 4 and yet the buoyancy flux profiles differ strongly between panels.

Another idea put forward in the atmospheric context is that the flux profiles in the equilibrated state resemble those of the most unstable mode inferred from linear stability analysis (Chai & Vallis, 2014; Green, 1970; Held & O'Brien, 1992). An issue with this approach is that only the equilibrated state is governed by the diffusion tensor (2), see the derivation in Meunier et al. (2023). In particular, Equation 2 indicates that the ratio of the vertical to the meridional buoyancy flux is given by the mean isopycnal slope $S$ (adiabatic transport). As discussed in Eady (1949) and Vallis (2017), this constraint does not hold for an unstable eigenmode because of the non-stationary terms, the associated profiles $\overline{wb}(z)$ and $\overline{vb}(z)$ being therefore incompatible with Equation 2 (in other words, one would infer a different profile for $K_{GM}(z)$ based on $\overline{vb}(z)$ or $\overline{wb}(z)$). We have nevertheless computed the most unstable eigenmode of the present Charney model, perturbatively for weak $\beta_*$ (see Supporting Information S1). As shown in Figure 2, the associated meridional buoyancy flux overpredicts the variations of $\overline{vb}(z)$ with depth and compares unfavorably with the present prediction (6). The most-unstable-mode approach may be better-suited for weakly nonlinear atmospheric states characterized by a weak supercriticality $\xi = 1/\beta_*$, as opposed to the present large-supercriticality oceanic situations (Jansen & Ferrari, 2012).

The success of the perturbative approximation $K_R = K_{GM}$ throughout the entire water column for case II above may come as a surprise to the reader accustomed to channel simulations, where $K_R$ typically exceeds $K_{GM}$ in the interior (see e.g., Abernathey et al., 2013). The reason for this success is that the meridional QGPV gradient is small in the interior and around the so-called “steering levels” (Abernathey et al., 2010, 2013; Green, 1970; Smith
& Marshall, 2009; Tréguier, 1999), making the QGPV flux $\overline{\nu q}$ negligible there (see e.g., Figure 6 of Smith and Marshall (2009)). One thus makes a negligible error by inferring the buoyancy and QGPV flux profiles using the approximation $K_R = K_{GM}$ throughout the entire water column.

The perturbative prediction (8) for the vertical structure of the GM coefficient is simple to implement, it is easily extended to a patch of ocean subject both to zonal and meridional large-scale gradients and shear flows, it is free of adjustable parameters—except for the overall magnitude of the transport—and it compares very favorably with the common practice of using a depth-invariant GM coefficient. The implementation of Equation 8 in a global model should lead to a more accurate description of the stratification of the Southern Ocean, and therefore of neighboring ocean basins. Beyond this modeling application, the physically-based vertical structure (9) for the buoyancy flux could be of use to infer the buoyancy flux throughout the entire water column based on near-surface data. Indeed, Figure 4 shows that the prediction (9) allows one to propagate the value of the near-surface flux to the interior of the water column in a way that agrees closely with the full DNS profile. By contrast, propagating the near-surface information using a uniform GM coefficient would lead to the orange line in Figure 4, which at depth typically departs from the DNS profile by 40%–100% depending on the situation.

## Data Availability Statement

The open-source numerical code is archived at https://doi.org/10.5281/zenodo.5458888 and documented in Miquel (2021) (see also github.com/BenMql/coral). The plots can be reproduced using the data archived on figshare at https://doi.org/10.6084/m9.fishshare.23864517.

## Acknowledgments

This research is supported by the European Research Council under grant agreement FLAVE 757239. The numerical study was performed using HPC resources from GENCI-CINES and TGCC (Grants 2021-A0102A10803, 2022-A0122A12489, and 2023-A0142A12489).

## References

* Abernathey, R., Ferreira, D., & Klocker, A. (2013). Diagnostics of isopycnal mixing in a circumpolar channel. *Ocean Modelling*, 72, 1–16. https://doi.org/10.1016/j.ocemod.2013.07.004

* Abernathey, R., Marshall, J., Mazloff, M., & Shuckburgh, E. (2010). Enhancement of mesoscale eddy stirring at steering levels in the southern ocean. *Journal of Physical Oceanography*, 40(1), 170–184. https://doi.org/10.1175/2009jpo4201.1

* Arbic, B., & Flierl, G. (2004a). Baroclinically unstable geostrophic turbulence in the limits of strong and weak bottom Ekman friction: Application to midocean eddies. *Journal of Physical Oceanography*, 34(10), 2257–2273. https://doi.org/10.1175/1520-0485(2004)034<2257:bugtit>2.0.co;2

* Arbic, B., & Flierl, G. (2004b). Effects of mean flow direction on energy, isotropy, and coherence of baroclinically unstable beta-plane geostrophic turbulence. *Journal of Physical Oceanography*, 34(1), 77–93. Retrieved from https://journals.ametsoc.org/downloadpdf/journals/ phoc/34/1/1520-0485_2004_034_0077_eomfdo_2.0.co_2.xml

* Arbic, B., & Scott, R. (2007). On quadratic bottom drag, geostrophic turbulence, and oceanic mesoscale eddies. *Journal of Physical Oceanogra- phy*, 38(1), 84–103. https://doi.org/10.1175/2007jpo3653.1

* Bretherton, F. P. (1966). Critical layer instability in baroclinic flows. *Quarterly Journal of the Royal Meteorological Society*, 92(393), 325–334. https://doi.org/10.1002/qj.49709239302

* Chai, J., & Vallis, G. K. (2014). The role of criticality on the horizontal and vertical scales of extratropical eddies in a dry GCM. *Journal of the Atmospheric Sciences*, 71(7), 2300–2318. https://doi.org/10.1175/jas-d-13-0351.1

* Chang, C.-Y., & Held, I. (2019). The control of surface friction on the scales of baroclinic eddies in a homogeneous quasigeostrophic two-layer model. *Journal of the Atmospheric Sciences*, 76(6), 1627–1643. https://doi.org/10.1175/jas-d-18-0333.1

* Chang, C.-Y., & Held, I. M. (2021). The parameter dependence of eddy heat flux in a homogeneous quasigeostrophic two-layer model on a $\beta$ plane with quadratic friction. *Journal of the Atmospheric Sciences*, 78(1), 97–106. https://doi.org/10.1175/jas-d-20-0145.1

* Charney, J. (1947). The dynamics of long waves in a baroclinic westerly current. *Journal of Meteorology*, 4(5), 135–162. https://doi.org/10.1175/1520-0469(1947)004<0136:tdlwi>2.0.co;2

* Danabasoglu, G., & Marshall, J. (2007). Effects of vertical variations of thickness diffusivity in an ocean general circulation model. *Ocean Model- ling*, 18(2), 122–141. https://doi.org/10.1016/j.ocemod.2007.03.006

* Dritschel, D., & McIntyre, M. (2008). Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers. *Journal of the Atmospheric Sciences*, 65(3), 855–874. https://doi.org/10.1175/2007jas2227.1

* Eady, E. T. (1949). Long waves and cyclone waves. *Tellus*, 1(3), 33–52. https://doi.org/10.3402/tellusa.v1i3.8507

* Ferreira, D., Marshall, J., & Heimbach, P. (2005). Estimating eddy stresses by fitting dynamics to observations using a residual-mean ocean circulation model and its adjoint. *Journal of Physical Oceanography*, 35(10), 1891–1910. https://doi.org/10.1175/jpo2785.1

* Flierl, G. (1978). Models of vertical structure and the calibration of two-layer models. *Dynamics of Atmospheres and Oceans*, 2(4), 342–381. https://doi.org/10.1016/0377-0265(78)90002-7

* Gallet, B., & Ferrari, R. (2020). The vortex gas scaling regime of baroclinic turbulence. *Proceedings of the National Academy of Sciences of the United States of America*, 117(9), 4491–4497. https://doi.org/10.1073/pnas.1916272117

* Gallet, B., & Ferrari, R. (2021). A quantitative scaling theory for meridional heat transport in planetary atmospheres and oceans. *AGU Advances*, 2(3), e2020AV000362. https://doi.org/10.1029/2020AV000362

* Gallet, B., Miquel, B., Hadjerci, G., Burns, K. J., Flierl, G. R., & Ferrari, R. (2022). Transport and emergent stratification in the equilibrated Eady model: The vortex-gas scaling regime. *Journal of Fluid Mechanics*, 948, A31. https://doi.org/10.1017/jfm.2022.501

* Gent, P. R. (2011). The Gent-McWilliams parameterization: 20/20 hindsight. *Ocean Modelling*, 39(1–2), 2–9. https://doi.org/10.1016/j.ocemod.2010.08.002

* Gent, P. R., & McWilliams, J. (1990). Isopycnal mixing in ocean circulation models. *Journal of Physical Oceanography*, 20(1), 150–155. https:// doi.org/10.1175/1520-0485(1990)020<0150:imiocm>2.0.co;2
Green, J. (1970). Transfer properties of the large-scale eddies and the general circulation of the atmosphere. *Quarterly Journal of the Royal Meteorological Society*, **96**(408), 157–185. https://doi.org/10.1002/qj.49709640802

Griffies, S. M. (1998). The Gent-Mcwilliams skew flux. *Journal of Physical Oceanography*, **28**(5), 831–841. https://doi.org/10.1175/1520-0485(1998)028<0831:tgmsf>2.0.co;2

Griffies, S. M., Gnanadesikan, A., Dixon, K. W., Dunne, J., Gerdes, R., Harrison, M. J., et al. (2005). Formulation of an ocean model for global climate simulations. *Ocean Science*, **1**(1), 45–79. https://doi.org/10.5194/os-1-45-2005

Held, I., & Larichev, V. (1996). A scaling theory for horizontally homogeneous, baroclinically unstable flow on a beta plane. *Journal of the Atmospheric Sciences*, **53**(7), 946–952. https://doi.org/10.1175/1520-0469(1996)053<0946:astfhh>2.0.co;2

Held, I., & O'Brien, E. (1992). Quasigeostrophic turbulence in a three-layer model: Effects of vertical structure in the mean shear. *Journal of the Atmospheric Sciences*, **49**(19), 1861–1870. https://doi.org/10.1175/1520-0469(1992)049<1861:qtiatl>2.0.co;2

Jansen, M., & Ferrari, R. (2012). Macroturbulent equilibration in a thermally forced primitive equation system. *Journal of the Atmospheric Sciences*, **69**(2), 695–713. https://doi.org/10.1175/jas-d-11-041.1

Larichev, V., & Held, I. (1995). Eddy amplitudes and fluxes in a homogeneous model of fully developed baroclinic instability. *Journal of Physical Oceanography*, **25**(10), 2285–2297. https://doi.org/10.1175/1520-0485(1995)025<2285:eaafia>2.0.co;2

McDougall, T. J., & McIntosh, P. C. (2001). The temporal-residual-mean velocity. Part II: Isopycnal interpretation and the tracer and momentum equations. *Journal of Physical Oceanography*, **31**(5), 1222–1246. https://doi.org/10.1175/1520-0485(2001)031<1222:ttrmvp>2.0.co;2

Meunier, J., Miquel, B., & Gallet, B. (2023). A direct derivation of the Gent-Mcwilliams/Redi diffusion tensor from quasi-geostrophic dynamics. *Journal of Fluid Mechanics*, **963**, A22. https://doi.org/10.1017/jfm.2023.347

Miquel, B. (2021). Coral: A parallel spectral solver for fluid dynamics and partial differential equations. *Journal of Open Source Software*, **6**(65), 2978. https://doi.org/10.21105/joss.02978

Nikurashin, M., & Vallis, G. (2011). A theory of deep stratification and overturning circulation in the ocean. *Journal of Physical Oceanography*, **41**(3), 485–502. https://doi.org/10.1175/2010jpo4529.1

Nikurashin, M., & Vallis, G. (2012). A theory of the interhemispheric meridional overturning circulation and associated stratification. *Journal of Physical Oceanography*, **42**(10), 1652–1667. https://doi.org/10.1175/jpo-d-11-0189.1

Phillips, N. (1954). Energy transformations and meridional circulations associated with simple baroclinic waves in a two-level, quasi-geostrophic model. *Tellus*, **6**(3), 274–286. https://doi.org/10.1111/j.2153-3490.1954.tb01123.x

Redi, M. H. (1982). Oceanic isopycnal mixing by coordinate rotation. *Journal of Physical Oceanography*, **12**(10), 1154–1158. https://doi.org/10.1175/1520-0485(1982)012<1154:aimbcr>2.0.co;2

Salmon, R. (1978). Two-layer quasigeostrophic turbulence in a simple special case. *Geophysical & Astrophysical Fluid Dynamics*, **10**(1), 25–52. https://doi.org/10.1080/03091927808242628

Salmon, R. (1980). Baroclinic instability and geostrophic turbulence. *Geophysical & Astrophysical Fluid Dynamics*, **15**(1), 157–211. https://doi.org/10.1080/03091928008241178

Salmon, R. (1998). *Lectures on geophysical fluid dynamics*. Oxford University Press.

Smith, K., & Marshall, J. (2009). Evidence for enhanced eddy mixing at middepth in the southern ocean. *Journal of Physical Oceanography*, **39**(1), 50–69. https://doi.org/10.1175/2008jpo3880.1

Stanley, Z., Bachman, S., & Grooms, I. (2020). Vertical structure of ocean mesoscale eddies with implications for parameterizations of tracer transport. *Journal of Advances in Modeling Earth Systems*, **12**(10), e2020MS002151. https://doi.org/10.1029/2020ms002151

Taylor, G. I. (1915). Eddy motion in the atmosphere. *Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character*, **215**(523–537), 1–26. https://doi.org/10.1175/1520-0493(1915)43<315c:emita>2.0.co;2

Thompson, A., & Young, W. (2006). Scaling baroclinic eddy fluxes: Vortices and energy balance. *Journal of Physical Oceanography*, **36**(4), 720–736. https://doi.org/10.1175/jpo2874.1

Thompson, A., & Young, W. (2007). Two-layer baroclinic eddy heat fluxes: Zonal flows and energy balance. *Journal of the Atmospheric Sciences*, **64**(9), 3214–3231. https://doi.org/10.1175/jas4000.1

Tréguier, A.-M. (1999). Evaluating Eddy mixing coefficients from eddy-resolving ocean models: A case study. *Journal of Marine Research*, **57**(1), 89–108. Retrieved from https://elischolar.library.yale.edu/cgi/viewcontent.cgi?article=3310&context=journal_of_marine_research

Vallis, G. (2017). *Atmospheric and oceanic fluid dynamics*. Cambridge University Press.

Venaille, A., Vallis, G. K., & Smith, K. S. (2011). Baroclinic turbulence in the ocean: Analysis with primitive equation and quasigeostrophic simulations. *Journal of Physical Oceanography*, **41**(9), 1605–1623. https://doi.org/10.1175/jpo-d-10-05021.1

Wolfe, C., & Cessi, P. (2010). What sets the strength of the middepth stratification and overturning circulation in eddying ocean models? *Journal of Physical Oceanography*, **40**(7), 1520–1538. https://doi.org/10.1175/2010jpo4393.1

Yankovsky, E., Zanna, L., & Smith, K. S. (2022). Influences of mesoscale ocean eddies on flow vertical structure in a resolution-based model hierarchy. *Journal of Advances in Modeling Earth Systems*, **14**(11), e2022MS003203. https://doi.org/10.1029/2022ms003203

Young, W. (2012). An exact thickness-weighted average formulation of the Boussinesq equations. *Journal of Physical Oceanography*, **42**(5), 692–707. https://doi.org/10.1175/jpo-d-11-0102.1

Zhang, W., & Wolfe, C. L. (2022). On the vertical structure of oceanic mesoscale tracer diffusivities. *Journal of Advances in Modeling Earth Systems*, **14**(6), e2021MS002891. https://doi.org/10.1029/2021ms002891

## References From the Supporting Information

Bouillaut, V., Miquel, B., Julien, K., Aumaitre, S., & Gallet, B. (2021). Experimental observation of the geostrophic turbulence regime of rapidly rotating convection. *Proceedings of the National Academy of Sciences*, **118**(44), e2105015118. https://doi.org/10.1073/pnas.2105015118

Burns, K., Vasil, G., Oishi, J., Lecoanet, D., & Brown, B. (2020). Dedalus: A flexible framework for numerical simulations with spectral methods. *Physical Review Research*, **2**, 023068. https://doi.org/10.1103/physrevresearch.2.023068

Miquel, B., Bouillaut, V., Aumaitre, S., & Gallet, B. (2020). On the role of the Prandtl number in convection driven by heat sources and sinks. *Journal of Fluid Mechanics*, **900**, R1. https://doi.org/10.1017/jfm.2020.485

Miquel, B., Lepot, S., Bouillaut, V., & Gallet, B. (2019). Convection driven by internal heat sources and sinks: Heat transport beyond the mixing-length or “ultimate” scaling regime. *Physical Review Fluids*, **4**(12), 121501. https://doi.org/10.1103/physrevfluids.4.121501