1. Introduction
Oceans and planetary atmospheres host currents or jets in thermal-wind balance with meridional buoyancy gradients. This situation is prone to baroclinic instability, however, and the resulting flows are strongly turbulent. In the ocean, this turbulence takes the form of ‘mesoscale’ eddies of size comparable to the Rossby deformation radius, a length scale of the order of 15–20 km in the Southern Ocean. While these vortices are key contributors to heat, salt and carbon transport, they are not resolved in state-of-the-art global climate models, and modellers need to parametrize the turbulent transport instead. It was soon realized that this turbulent transport is ill-described by standard horizontal diffusion (Gent Reference Gent2011). Instead, rapid rotation and strong stratification induce quasi-geostrophic (QG) flows in the ocean interior that transport and mix tracers predominantly along density surfaces. Redi (Reference Redi1982) thus proposed to relate the sub-grid fluxes to the large-scale gradients through a diffusion tensor that mixes tracers along mean density surfaces only. While avoiding spurious cross-isopycnal mixing, such a parametrization scheme is unable to describe the transport of buoyancy. Gent & McWilliams (Reference Gent and McWilliams1990) subsequently proposed to add an extra eddy-induced advective flux, ensuring the transport of buoyancy (Gent et al. Reference Gent, Willebrand, McDougall and McWilliams1995). Following Griffies (Reference Griffies1998), these two contributions are combined conveniently into a ‘Gent–McWilliams/Redi’ (GM/R) diffusion tensor that includes both a symmetric part, Redi's isoneutral diffusion scheme, and an antisymmetric part encoding the GM advective contribution.
A general derivation of the GM/R diffusion tensor was provided by McDougall & McIntosh (Reference McDougall and McIntosh2001). Using a thickness-weighted average (TWA) together with an expansion for small fluctuations, they showed that the (GM/R) diffusion tensor governs the eddy-induced fluxes arising in the evolution equation for the thickness-weighted average of the tracer concentration. The GM flux then corresponds to advection by a ‘quasi-Stokes’ velocity, which is related directly to the GM transport coefficient 
$K_{GM}$ and can be added readily to the Eulerian mean velocity to form the residual velocity (Andrews & McIntyre Reference Andrews and McIntyre1976). More recently, Young (Reference Young2012) relaxed the small-fluctuation assumption and derived an exact TWA formulation of the Boussinesq equations, again framed in terms of a three-dimensional residual velocity that includes the eddy-induced advection. Eddy-forcing of the residual velocity arises through the divergence of Eliassen–Palm vectors, which can be recast in terms of the potential vorticity fluxes using some very general form of the Taylor–Bretherton relation (see also Maddison & Marshall (Reference Maddison and Marshall2013) and the discussion around (5.4) below).
While they should guide the implementation of parametrizations into global ocean models, such very general TWA formulations do not provide prescriptions for the magnitude or depth dependence of the eddy-induced fluxes. Instead, progress can be made regarding the overall magnitude and structure of the eddy-induced transport by leveraging the rapid global rotation and focusing on the simpler QG system. Historically, this approach has led to great insight into the overall magnitude of the transport induced by baroclinic turbulence, based on the study of the two-layer QG model (Phillips Reference Phillips1954; Salmon Reference Salmon1978, Reference Salmon1980; Larichev & Held Reference Larichev and Held1995; Held & Larichev Reference Held and Larichev1996; Arbic & Flierl Reference Arbic and Flierl2004a,Reference Arbic and Flierlb; Thompson & Young Reference Thompson and Young2006, Reference Thompson and Young2007; Arbic & Scott Reference Arbic and Scott2007; Chang & Held Reference Chang and Held2019; Gallet & Ferrari Reference Gallet and Ferrari2020, Reference Gallet and Ferrari2021). In a similar fashion, the goal of the present study is to understand the structure of the diffusion tensor in a simple and strongly idealized situation. We thus consider the eddy-induced meridional and vertical transport arising from the QG dynamics of a three-dimensional horizontally homogeneous vertically sheared zonal flow on the 
$\beta$ plane; see figure 1. Baroclinic instability of the base state leads rapidly to turbulence, and we wish to characterize the transport properties of the resulting equilibrated state. Focusing on a QG system allows one to make progress regarding the structure of the diffusion tensor, specifically as follows.
(1) The TWA coincides with the standard Eulerian average at fixed depth
$z$ in the QG system (see Appendix B), which allows for a particularly compact derivation of the GM/R diffusion tensor.(2) In the QG system, there is an inverse cascade of energy and buoyancy variance in the interior, with no anomalous dissipation or mixing. Hence no diapycnal diffusivity coefficient appears in the Redi tensor in the limit of vanishing viscosity and small-scale diffusivities. Together with the (statistical) zonal invariance, this leads to a GM/R diffusion tensor that involves two vertically dependent transport coefficients only, the GM coefficient
$K_{GM}(z)$ and the Redi diffusivity $K_R(z)$. By contrast, the TWA can be performed for both rotating and non-rotating stratified turbulence, and the theory does not dictate a priori whether diapycnal fluxes arise. (Such diapycnal fluxes certainly arise for non-rotating stratified turbulence; see Linden Reference Linden1979; Peltier & Caulfield Reference Peltier and Caulfield2003; Maffioli, Brethouwer & Lindborg Reference Maffioli, Brethouwer and Lindborg2016; Caulfield Reference Caulfield2021.)(3) The QG expression for potential vorticity leads to the Taylor–Bretherton relation between the profiles of
$K_{GM}$ and $K_R$ in the interior of the domain.(4) Additionally, the QG boundary conditions at top and bottom readily provide constraints on the top and bottom values of
$K_{GM}$ and $K_R$. Namely, these two coefficients are equal at top and bottom for low bottom drag.
Figure 1. (a) 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 profile 
$N(z)$ for the buoyancy frequency. The background zonal shear flow has profile 
$U(z)$. The flow coexists with a background meridional buoyancy gradient as a result of thermal-wind balance. Bottom friction damps the fluctuating kinetic energy. (b) Snapshot of the departure buoyancy field 
$b$ in the equilibrated state of the numerical simulation.
We introduce the theoretical set-up in § 2. We highlight the main conservation relations in § 3, from which we derive the GM/R diffusion tensor in § 4. In § 5, we identify constraints relating the GM coefficient and the Redi diffusivity. We illustrate these constraints using direct numerical simulation in § 6, before concluding in § 7. In Appendix A, we show that the contributions from the small-scale diffusive terms are negligible in the QG regime. In Appendix B, we make the connection between the present QG approach and the more general TWA approach of McDougall & McIntosh (Reference McDougall and McIntosh2001). Finally, the details of the numerical procedure are provided in Appendix C.
2. Quasi-geostrophic dynamics of an idealized three-dimensional patch of ocean
We consider the idealized patch of ocean represented in figure 1. Water occupies a volume 
$(x,y,z)\in [0,L]^2 \times [-H,0]$ with a stress-free boundary at the surface 
$z=0$, and a linear-friction boundary condition at 
$z=-H$. The fluid layer is subject to global rotation around the vertical axis with a local Coriolis parameter 
$f_0+\beta y$, where 
$y$ denotes the meridional (north–south) coordinate. Additionally, the fluid layer is density-stratified with an arbitrary buoyancy frequency profile 
$N(z)$, and we restrict attention to a single stratifying agent. We focus on the rapidly rotating strongly stratified regime for which the fluid motion is governed by quasi-geostrophy (Salmon Reference Salmon1998; Vallis Reference Vallis2006; Venaille, Vallis & Smith Reference Venaille, Vallis and Smith2011). In that limit, the velocity field 
$\boldsymbol {u}=(u,v,w)$ consists of a leading-order horizontal geostrophic flow 
$(u,v)=(-P_y,P_x)$, where the generalized pressure field 
$P$ is defined as the opposite of the streamfunction, together with subdominant vertical velocity 
$w$. The buoyancy field is given by 
$B=f_0 P_z$ as a result of hydrostatic balance.
The base flow consists of an arbitrary zonal velocity profile 
$U(z)=-P_y$. Differentiating with respect to 
$z$ indicates that the zonal flow is in thermal-wind balance with a 
$z$-dependent meridional buoyancy gradient, 
$\partial _y B=-f_0\,U'(z)$. We consider the evolution of arbitrary departures from this base state. 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, and 
$b=f_0 \, p_z$ the departure buoyancy. In the following, we adopt dimensionless variables, with time expressed in units of 
$|f_0|^{-1}$ and lengths in units of 
$H$. For brevity, we use the same symbols for the dimensionless variables. The QG limit is obtained for small isopycnal slope of the base state, or equivalently, in the small Rossby number limit for 
${{O}} (1)$ stratification. Denoting as 
$\epsilon$ the typical magnitude of the isopycnal slope, the QG system can be derived through the following scalings:
\begin{equation} \left.\begin{array}{@{}c@{}} N^2 \sim 1, \quad \beta \sim \epsilon , \quad \partial_x, \partial_y, \partial_z \sim 1 , \quad \partial_t \sim \epsilon , \\ U \sim \epsilon , \quad (u,v) \sim \epsilon , \quad w \sim \epsilon^2 , \quad b \sim \epsilon . \end{array}\right\} \end{equation}
A standard asymptotic expansion of the equations of motion leads to the QG system (Pedlosky Reference Pedlosky1979; Salmon Reference Salmon1998; Vallis Reference Vallis2006), as recalled in Gallet et al. (Reference Gallet, Miquel, Hadjerci, Burns, Flierl and Ferrari2022) for the specific notations and scalings considered here. The evolution then reduces to a conservation equation for the QG potential vorticity (QGPV). The dimensionless QGPV departure 
$q$ is related to the departure pressure 
$p$ through
\begin{equation} q = \varDelta_\perp p + \partial_z \left( \frac{p_z}{N^2(z)} \right) , \end{equation}
where 
$\varDelta _\perp =\partial _{xx}+\partial _{yy}$. The (dimensionless) QGPV conservation equation reads
\begin{equation} \partial_t q + U(z)\,q_x + J(p,q) = [-{\beta}+{\mathcal{S}}'(z)] p_x + {\mathcal{D}}_q , \end{equation}
where the Jacobian is 
$J(g,h)=g_x h_y - g_y h_x$, and we denote the isopycnal slope of the base state as 
${\mathcal {S}}(z)=U'/N^2$. The first term on the right-hand side of (2.3) corresponds to the distortion of the background meridional potential vorticity (PV) gradient by the meridional flow. The second term, 
${\mathcal {D}}_q$, is the contribution from the viscosity and buoyancy diffusivity, which damp the small-scale vorticity and buoyancy fluctuations. This term is detailed in Appendix A.
At the same level of approximation, the evolution equation for the (dimensionless) buoyancy departure 
$b=p_z$ reads
\begin{equation} \partial_t b + U(z)\,b_x + J(p,b) = U'(z)\,p_x - w\,N^2(z) + {\mathcal{D}}_b , \end{equation}
where the diffusive term 
${\mathcal {D}}_b$ is provided in Appendix A. The first two terms on the right-hand side are the sources of buoyancy fluctuations in the system; they correspond to the distortion by the turbulent flow of the background meridional and vertical buoyancy gradients, respectively. At the surface, where 
$w=0$, this equation reduces to
\begin{equation} \partial_t b|_0 + U(0)\,b_{x}|_0 + J(p|_0, b|_0) = U'(0)\,p_x|_0 + {\mathcal{D}}_b|_0 , \end{equation}
where 
$0$ refers to quantities evaluated at 
$z=0$. Quantities evaluated just above the bottom Ekman boundary layer are denoted with the subscript 
$-1^+$. At this depth, the pumping vertical velocity is given by 
$w|_{-1^+}=\kappa \varDelta _{\perp}p|_{-1^+}$, where the friction coefficient 
$\kappa$ can be either related to the vertical viscosity through standard Ekman layer theory over a flat bottom boundary, or specified at the outset as an independent coefficient parametrizing more realistic drag on the ocean floor (see Appendix C and Gallet et al. Reference Gallet, Miquel, Hadjerci, Burns, Flierl and Ferrari2022). The evolution equation for the buoyancy at 
$z=-1^+$ reads
\begin{align} &\partial_t b|_{{-}1^+} + U({-}1)\,b_{x}|_{{-}1^+} + J(p|_{{-}1^+},b|_{{-}1^+}) \nonumber\\ &\quad = U'({-}1)\,p_x|_{{-}1^+} - N^2({-}1)\,\kappa \varDelta_{{\perp}} p|_{{-}1^+} + {\mathcal{D}}_b|_{{-}1^+}. \end{align}
One way to integrate the QG system consists in time-stepping the QGPV conservation equation (2.3) together with the top and bottom buoyancy evolutions (2.5) and (2.6). To infer the pressure field at each time step, one inverts the relation (2.2) using 
$b|_{-1^+}$ and 
$b|_{0}$ as boundary conditions. A desirable feature of this QG approach is that it is fully compatible with periodic boundary conditions in 
$x$ and 
$y$ for the departure fields. We adopt such periodic boundary conditions in the following.
In the bulk of the domain, the buoyancy evolution (2.4) provides a diagnostic relation to infer the subdominant vertical velocity 
$w$, the latter being crucial to parametrize eddy-induced vertical transport.
3. Material invariants: buoyancy, QGPV and the cross-invariant
We denote with an overbar 
$\bar {\cdot }$ a time average together with a horizontal area average. Our goal is to characterize the transport properties of the flow, more specifically the diffusion tensor connecting the eddy-induced fluxes to the large-scale background gradients of some arbitrary tracer, be it active or passive. One can gain insight into the structure of this tensor by focusing on two specific tracers: buoyancy and QGPV. In this section, we thus derive rigorous constraints between the meridional and vertical turbulent fluxes of buoyancy and QGPV: 
$\overline {vb}(z)$, 
$\overline {wb}(z)$, 
$\overline {vq}(z)$ and 
$\overline {wq}(z)$.
The first constraint stems from the conservation of buoyancy variance. Multiplying the buoyancy evolution (2.4) with 
$b$ before averaging over time, 
$x$ and 
$y$ leads to
\begin{equation} N^2\,\overline{wb} = U'\,\overline{vb} , \end{equation}up to diffusive corrections that vanish in the regime of low viscosity and buoyancy diffusivity. A proof that the diffusive contributions indeed vanish is provided in Appendix A, based on the well-known absence of a forward energy cascade in QG turbulence. We recast the equality (3.1) in the form
\begin{equation} \overline{wb} = {\mathcal{S}}(z)\,\overline{vb} , \end{equation}which shows that the mean buoyancy transport is directed along the mean isopycnals.
We derive a second constraint on the fluxes based on the conservation of the cross-invariant 
$bq$. Multiply the QGPV evolution (2.3) with 
$b$ and the buoyancy conservation (2.4) with 
$q$. Summing the resulting equations and averaging over time, 
$x$ and 
$y$ leads to
\begin{equation} \overline{wq} = {\mathcal{S}}(z)\,\overline{vq} - \frac{\beta-{\mathcal{S}}'(z)}{N^2(z)}\,\overline{vb} , \end{equation}with the equality holding in the low-diffusivity limit, where, as shown in Appendix A, the contributions from viscosity and buoyancy diffusivity vanish.
4. Arbitrary tracer: Gent–McWilliams/Redi diffusion tensor
Define the 
$z$-dependent eddy diffusivities 
$K_R(z)$ and 
$K_{GM}(z)$ as
\begin{equation} K_R(z) = \frac{\overline{vq}}{-\beta+{\mathcal{S}}'(z)} \quad \text{and} \quad K_{GM}(z) = \frac{\overline{vb}}{U'(z)} , \end{equation}
namely, 
$K_R$ is the ratio of the meridional PV flux to minus the background meridional PV gradient, while 
$K_{GM}$ is the ratio of the meridional buoyancy flux to minus the background meridional buoyancy gradient. The reasons for the notations 
$K_{GM}$ and 
$K_R$ will become obvious at the end of the derivation below.
Now consider a (passive or active) tracer 
$\tau$ stirred by the three-dimensional flow and subject to horizontally uniform gradients 
$G_y(z)={{O}}(\epsilon )$ and 
$G_z(z)={{O}}(1)$ in the meridional and vertical directions, respectively. That is, the total tracer field reads
\begin{equation} \int_{0}^z G_z(\tilde{z}) \, \mathrm{d}\tilde{z} + y\,G_y(z) + \tau(x,y,z,t) . \end{equation} Under these conditions and with the scalings (2.1), the evolution equation for 
$\tau$ reads
\begin{equation} \partial_t \tau + U(z)\,\tau_x + J(p,\tau) ={-} p_x\,G_y(z) - w\,G_z(z) + {\mathcal{D}}_\tau , \end{equation}
where 
${\mathcal {D}}_\tau$ denotes small-scale diffusion. A few remarks are in order regarding the background meridional and vertical gradients: 
$G_y$ and 
$G_z$ above should be understood as the lowest-order background gradients that enter QG dynamics. Naturally, a subdominant vertical gradient 
$G_z^{(1)}= y\,\partial _z G_y(z)={{O}}(\epsilon )$ exists to ensure the equality of the cross-derivatives 
$\partial _y(G_z+G_z^{(1)})=\partial _z G_y$. However, one can check easily that 
$G_z^{(1)}$ is subdominant and does not arise in the QG evolution (4.3). Similarly, keeping 
$G_z^{(1)}$ on the right-hand side of (4.6) would lead to negligible corrections to the fluxes, of higher order in 
$\epsilon$.
The meridional and vertical fluxes of 
$\tau$ are related to the background meridional and vertical gradients 
$G_y$ and 
$G_z$ through a diffusion tensor
\begin{equation} \left(\begin{matrix} \overline{v\tau} \\ \overline{w\tau} \end{matrix} \right) =\left[\begin{matrix} A_1(z) & A_2(z) \\ A_3(z) & A_4(z) \\ \end{matrix} \right] \left( \begin{matrix} G_y \\ G_z \end{matrix} \right) ,\end{equation}
where the 
$A_i(z)$ are unknown 
$z$-dependent coefficients at this stage. Apply relation (4.4) to the tracers 
$b$ and 
$q$, the associated background gradients being inferred readily from the right-hand-side terms of (2.3) and (2.4): 
$G_y=-U'(z)$ and 
$G_z=N^2$ for 
$\tau =b$, and 
$G_y=\beta - {\mathcal {S}}'(z)$ and 
$G_z=0$ for 
$\tau =q$. We obtain the following fluxes:
\begin{equation} \left( \begin{matrix} \overline{vb} \\ \overline{wb} \end{matrix} \right) = \left( \begin{matrix} -A_1\,U'(z) + A_2 N^2\\ -A_3\,U'(z) + A_4 N^2 \end{matrix} \right),\quad \left( \begin{matrix} \overline{vq} \\ \overline{wq} \end{matrix} \right) = \left( \begin{matrix} A_1 [\beta - {\mathcal{S}}'(z)] \\ A_3 [\beta - {\mathcal{S}}'(z)] \end{matrix} \right) . \end{equation}
There are four constraints on these four fluxes, which allow us to express the four coefficients 
$A_i$ in terms of 
$K_{GM}(z)$ and 
$K_{R}(z)$: the first two constraints are simply the definitions (4.1) of 
$K_{GM}$ and 
$K_{R}$, the third constraint is (3.2), namely the fact that the mean transport of 
$b$ is along the mean isopycnals, and the fourth constraint is the cross-invariant relation (3.3). After a straightforward calculation of the coefficients 
$A_i$, the diffusion tensor connecting the fluxes to the background gradients becomes
\begin{equation} \left( \begin{matrix} \overline{v\tau} \\ \overline{w\tau} \end{matrix} \right) =\left[\begin{matrix} -K_R & (K_{GM}-K_R) {\mathcal{S}}\\ - (K_{GM}+K_R){\mathcal{S}} & - K_R {\mathcal{S}}^2\\ \end{matrix} \right] \left( \begin{matrix} G_y \\ G_z \end{matrix} \right) .\end{equation}
This form for the diffusion tensor corresponds to the GM/R parametrization, where 
$K_{GM}(z)$ denotes the 
$z$-dependent GM coefficient, and 
$K_R(z)$ denotes the Redi diffusivity. The former represents the skew flux associated with adiabatic transport by the eddying flow (Griffies Reference Griffies1998). Using the definition in (4.1) of 
$K_{GM}$, we check in Appendix B that the GM part of the tensor (4.6) corresponds to the QG limit of the advective fluxes associated with the more general quasi-Stokes streamfunction introduced by McDougall & McIntosh (Reference McDougall and McIntosh2001). The Redi part of the tensor represents mixing along the neutral direction in the limit of weak isopycnal slope 
${\mathcal {S}}(z)$. As can be inferred from the QGPV conservation equation, the Redi diffusivity 
$K_R(z)$ also equals the Taylor–Kubo eddy diffusivity coefficient deduced at any height 
$z$ from the Lagrangian correlation function of the horizontal QG velocity field. Several similar estimates for the PV diffusivity have been compared in channel geometry by Abernathey, Ferreira & Klocker (Reference Abernathey, Ferreira and Klocker2013).
5. Constraints on the GM and Redi coefficients
The derivation above allows us to obtain constraints on the GM and Redi coefficients as defined by (4.1). First, at the upper boundary the buoyancy equation (2.5) has exactly the same structure as the QGPV conservation equation (2.3). Both equations correspond to advection by the base zonal flow and the QG flow at the upper boundary, with a source term that corresponds to the distortion of a background meridional gradient. We conclude that the diffusivities relating the meridional flux to the background meridional gradient are equal for 
$q$ and 
$b$ at 
$z=0$ (and given by the Taylor–Kubo eddy diffusivity coefficient associated with the surface horizontal flow). Using the definitions (4.1), this leads to the constraint
\begin{equation} K_{GM}(0) = K_R(0) . \end{equation}The same relation holds near the bottom boundary (just above the Ekman boundary layer) when the drag coefficient is low:
\begin{equation} K_{GM}({-}1^+) = K_R({-}1^+) . \end{equation}
The equality of 
$K_{GM}(z)$ and 
$K_R(z)$ is a common assumption when implementing the GM/R parametrization in global models (Griffies Reference Griffies1998). We put this assumption on a firmer analytical footing by showing that it holds near the upper and lower boundaries, although the two coefficients differ generically in the interior of the fluid column. In the general case, however, we can further relate the vertical dependence of 
$K_{GM}(z)$ and 
$K_R(z)$ through the Taylor–Bretherton relation. Multiplying (2.2) with 
$v=p_x$ before averaging horizontally yields, after a few integrations by parts in the horizontal directions,
\begin{equation} \overline{vq}=\frac{\mathrm{d}}{\mathrm{d} z} \left( \frac{\overline{vb}}{N^2} \right) .\end{equation}This equation corresponds to the horizontal average of a more general relation derived initially by Bretherton (Reference Bretherton1966) and often referred to as the Taylor–Bretherton relation (Taylor Reference Taylor1915; Dritschel & McIntyre Reference Dritschel and McIntyre2008; Young Reference Young2012). Using the definitions (4.1), we recast (5.3) as
\begin{equation} K_R({\mathcal{S}}'-\beta)=\frac{\mathrm{d}}{\mathrm{d} z} (K_{GM} {\mathcal{S}}) . \end{equation}This equality was obtained previously by several authors (e.g. Smith & Marshall Reference Smith and Marshall2009) and is recalled here for the sake of completeness only.
6. A numerical example
With the goal of illustrating the results above, we turn to the direct numerical simulation of an isolated horizontally homogeneous patch of ocean. The set-up is chosen to reproduce conditions in the Antarctic Circumpolar Current (ACC), with surface-intensified stratification, shear and turbulence. The base state consists of a (dimensionless) stratification decreasing linearly with depth, from 
$N^2(0)=400$ at the surface to 
$N^2(-1)=50$ at the bottom, together with an exponential profile for the background shear flow, 
$U(z)=Ro \, \textrm {e}^{2z}$ with 
$Ro=0.03$, and a dimensionless planetary vorticity gradient 
$\beta =4.0 \times 10^{-5}$. These values for 
$Ro$ and 
$\beta$ are ten times smaller than typical values in the ACC; this choice leaves invariant the dissipation-free QG dynamics while ensuring that the numerical simulation is indeed performed in the fully QG regime. In other words, we simulate an idealized horizontally homogeneous and fully QG version of the ACC. We have also kept 
$f_0>0$, which leads conveniently to 
$\overline {vb}>0$ while being equivalent to the ACC situation up to an equatorial symmetry. We solve for the fully nonlinear evolution of the departures from the base state inside a domain of dimensionless size 
$(x,y,z)\in [0,500]^2\times [-1,0]$ using periodic boundary conditions in the horizontal directions and small values for the dissipative coefficients. The numerical procedure is detailed in Appendix C.
After some transient, the system reaches a statistically steady equilibrated state, illustrated in figure 1 through a snapshot of the departure buoyancy field 
$b$. We extract the time and horizontal area averages of the meridional and vertical fluxes of buoyancy and QGPV in this statistically steady state. The corresponding profiles are shown in figure 2, where figure 2(a) provides the diagnosed GM coefficient and Redi diffusivity. In agreement with results from re-entrant channel simulations (Abernathey et al. Reference Abernathey, Ferreira and Klocker2013), 
$K_{GM}(z)$ is monotonic in the interior of the domain, whereas 
$K_R(z)$ exhibits a maximum at mid-depth (Tréguier Reference Tréguier1999; Smith & Marshall Reference Smith and Marshall2009). In line with the constraints (5.1) and (5.2), the interior profiles of 
$K_{GM}(z)$ and 
$K_R(z)$ tend to a common limiting value as we approach the top or bottom boundary. This tendency is disrupted by diffusive boundary layer effects in the immediate vicinity of the boundaries (more strongly so near the surface). These boundary layers shrink as we lower the diffusivities employed in the numerical simulation.
Figure 2. Time and horizontally averaged profiles from the numerical simulation. (a) Transport coefficients as defined in (4.1) (grey region corresponding to the depth where the meridional QGPV gradient vanishes).(b) Vertical buoyancy flux compared to the GM/R prediction, which validates adiabatic transport in the interior. (c) Vertical QGPV flux compared to the GM/R prediction, which validates the cross-invariant relation (3.3).
Having diagnosed 
$K_{GM}(z)$ and 
$K_R(z)$, we turn to the vertical fluxes of buoyancy and QGPV, with the goal of comparing the numerical fluxes to the predictions of the GM/R diffusion tensor. Figure 2(b) shows that the vertical buoyancy flux 
$\overline {wb}$ is captured accurately by the GM/R prediction 
${\mathcal {S}}K_{GM}\,U'(z)$. Because 
${\mathcal {S}}K_{GM}\,U'(z)={\mathcal {S}} \overline {vb}$, this validates the fact that buoyancy is transported adiabatically in the interior, in line with (3.2). Figure 2(c) shows that the vertical QGPV flux 
$\overline {wq}$ is captured accurately by the GM/R prediction 
$-{\mathcal {S}}(K_{GM}+K_R)(\beta -{\mathcal {S}}')$, the latter expression being also equal to the right-hand side of the cross-invariant relation (3.3). The good agreement in figure 2(c) thus validates the simple and exact cross-invariant relation (3.3) in the interior of the domain.
7. Conclusion
We have studied the transport properties of the turbulent QG flow arising from the baroclinic instability of a horizontally homogeneous vertically sheared zonal current. While less general than the TWA formulation of the Boussinesq equations (McDougall & McIntosh Reference McDougall and McIntosh2001; Young Reference Young2012), the QG limit allows one to make progress on the structure of the diffusion tensor relating the eddy-induced fluxes to the background gradients. Based on the conservation of buoyancy variance and of a cross-invariant involving buoyancy and QGPV, we thus derived a particularly simple GM/R form for the diffusion tensor. First, in the interior of the domain, there are no diapycnal fluxes, provided that the viscosity and small-scale diffusivities are small. The diffusion tensor then involves only two vertically dependent coefficients: the GM transport coefficient 
$K_{GM}(z)$, and the Redi diffusivity 
$K_R(z)$. Second, based on the definition of QGPV, one can relate 
$K_{GM}$ and 
$K_R$ through the Taylor–Bretherton relation (5.4). Finally, based on the QGPV and buoyancy evolution equations, one obtains that 
$K_{GM}$ and 
$K_R$ are equal to one another at top and bottom. These results provide some support for the common modelling assumption 
$K_{GM} \simeq K_R$ near the boundaries (Griffies Reference Griffies1998). However, the two coefficients are allowed to depart from one another in the interior of the fluid column, and indeed they do in the present numerical simulation (in line with previous studies; see Abernathey et al. Reference Abernathey, Ferreira and Klocker2013). It would be interesting to investigate whether some equivalent of the boundary relation 
$K_{GM} \simeq K_R$ exists beyond the present idealized QG framework, for a primitive equation or Boussinesq system. TWA would likely play a central role for such an extension.
Acknowledgements
We thank G. Hadjerci, R. Ferrari and W.R. Young for insightful discussions.
Funding
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-A0102A12489, 2022-A0122A12489, 2022-A0122A10803 and 2023-A0142A12489).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Diffusive contributions
We consider the impact of the standard diffusive terms (viscosity and buoyancy diffusivity) within the framework of a QG system. We use different coefficients for the diffusivities in the horizontal and vertical directions: 
$E_{\perp }$ and 
$E_z$, respectively, for the horizontal and vertical dimensionless viscosities (Ekman numbers), and 
$E_{b,\perp }$ and 
$E_{b,z}$, respectively, for the horizontal and vertical dimensionless buoyancy diffusivities. With these notations, the diffusive term in the buoyancy equation (2.4) reads
\begin{equation} {\mathcal{D}}_b = E_{b,\perp} \varDelta_{{\perp}} b + E_{b,z}\,\partial_{zz} b , \end{equation}while the diffusive term in the QGPV conservation equation (2.3) reads
\begin{equation} {\mathcal{D}}_q = E_{{\perp}} \varDelta_{{\perp}}^2 p + E_z \varDelta_{{\perp}} p_{zz} + \partial_z \left[ \frac{1}{N^2} \left( E_{b,\perp} \varDelta_{{\perp}} b + E_{b,z} \partial_{zz} b \right) \right] . \end{equation}
In contrast to standard three-dimensional turbulence, QG dynamics are characterized by an inverse energy cascade (Charney Reference Charney1971), together with a forward cascade of buoyancy variance at the boundaries only: in the limit where the various diffusive coefficients 
$E_i$ are sent to zero simultaneously, there is no ‘anomalous’ energy dissipation and no ‘anomalous’ dissipation of buoyancy variance in the interior (Lapeyre Reference Lapeyre2017). That is, the limit 
$E_i\,\overline {(\varDelta _\perp p)^2} \to 0$ holds for any 
$z$, and the limits 
$E_i\, \overline {|\boldsymbol {\nabla }_\perp b|^2} \to 0$ and 
$E_i\, \overline {(b_{z})^2} \to 0$ hold pointwise for 
$z \neq -1, 0$. Additionally, any 
$z$-derivative of these profiles also vanishes in the vanishing-diffusivity limit for 
$z \neq -1, 0$. The diffusive contribution to the right-hand side of the adiabatic transport relation (3.2) reads
\begin{equation} \overline{b{\mathcal{D}}_b} ={-} E_{b,\perp}\, \overline{|\boldsymbol{\nabla}_\perp b|^2} - E_{b,z}\, \overline{(b_{z})^2} + \frac{E_{b,z}}{2}\,\frac{\mathrm{d}}{\mathrm{d}z} \overline{b^2} . \end{equation}
The three terms vanish in the vanishing-diffusivity limit for 
$z \neq -1, 0$, leading to 
$\overline {b{\mathcal {D}}_b} \to 0$ and (3.2). The diffusive contribution to the right-hand side of the cross-invariant relation (3.3) reads 
$(\overline {b{\mathcal {D}}_q}+\overline {q{\mathcal {D}}_b})/N^2$, where
\begin{align}
\overline{q{\mathcal{D}}_b} &=
\frac{E_{b,\perp}}{2}\,\frac{\mathrm{d}}{\mathrm{d}z}
\overline{(\varDelta_{{\perp}} p)^2} + \frac{E_{b,z}}{2}
\left[- \frac{\mathrm{d}^3}{\mathrm{d}z^3}
\overline{(\boldsymbol{\nabla}_\perp p)^2} +
3\,\frac{\mathrm{d}}{\mathrm{d}z}
\overline{(\boldsymbol{\nabla}_\perp b)^2} \right]
\nonumber\\ &\quad + E_{b,\perp} \left\{ -
\frac{\mathrm{d}}{\mathrm{d}z} \left[
\frac{\overline{(\boldsymbol{\nabla}_\perp b)^2}}{N^2}
\right] + \frac{1}{2N^2}\,\frac{\mathrm{d}}{\mathrm{d}z}
\overline{(\boldsymbol{\nabla}_\perp b)^2}
\right\}\nonumber\\ &\quad +E_{b,z} \left\{
\left[\frac{\mathrm{d}^2}{\mathrm{d}z^2}
\left(\frac{\overline{b^2}}{2}\right) -\overline{(b_z)^2}
\right] \frac{\mathrm{d}}{\mathrm{d}z} \left( \frac{1}{N^2}
\right) +\frac{1}{2 N^2}\, \frac{\mathrm{d}}{\mathrm{d}z}
\overline{b_z^2} \right\},
\end{align}
\begin{align}
\overline{b{\mathcal{D}}_q} & =
\frac{E_\perp}{2}\,\frac{\mathrm{d}}{\mathrm{d}z}
\overline{(\varDelta_\perp p)^2} -
\frac{E_z}{2}\,\frac{\mathrm{d}}{\mathrm{d}z}
\overline{(\boldsymbol{\nabla}_\perp b)^2} \nonumber\\
&\quad + \frac{\mathrm{d}}{\mathrm{d}z}
\left[-\frac{E_{b,\perp}}{N^2}\,
\overline{(\boldsymbol{\nabla}_\perp b)^2} +
\frac{E_{b,z}}{2 N^2}\,\frac{\mathrm{d}^2}{\mathrm{d}z^2}
\overline{b^2} -
\frac{E_{b,z}}{N^2}\,\overline{(b_z)^2}\right] \nonumber\\
&\quad +
\frac{E_{b,\perp}}{2 N^2}\,\frac{\mathrm{d}}{\mathrm{d}z}
\overline{(\boldsymbol{\nabla}_\perp b)^2} -
\frac{E_{b,z}}{2 N^2}\,\frac{\mathrm{d}}{\mathrm{d}z}
\overline{(b_z)^2} .
\end{align}
The terms on the right-hand sides of both expressions vanish in the vanishing-diffusivity limit for 
$z \neq -1, 0$, leading to 
$\overline {q{\mathcal {D}}_b} + \overline {b{\mathcal {D}}_q} \to 0$. Hence the approximate relation (3.3) for low diffusivities.
Appendix B. Connection to the TWA and the residual mean approach
B.1. The evolution equations for the TWA and for the standard fixed-$z$ averaged tracer concentration are identical in the QG limit
For a given tracer 
$\tau$, McDougall & McIntosh (Reference McDougall and McIntosh2001) consider the evolution equations for the TWA 
$\hat {\tau }$ and for the standard fixed-
$z$ average 
$\bar {\tau }$, where the time average is to be understood as an average over a few eddy turnover times. They show that the two evolution equations differ by the divergence of some vector 
$\boldsymbol {E}$, see their equations (53–55). With the QG scalings (2.1), however, this additional term is negligible. The horizontal components of 
$\boldsymbol {E}$ are 
${{O}}(\epsilon ^3)$, smaller than the meridional flux 
$\overline {v \tau }={{O}}(\epsilon ^2)$ discussed in the present study. The vertical component of 
$\boldsymbol {E}$ is the time derivative of some 
${{O}}(\epsilon ^2)$ material invariant, averaged over the slow QG eddy turnover time scale 
$1/\epsilon$. The latter time derivative is thus of order at least 
$\epsilon ^4$, much smaller than the vertical flux 
$\overline {w \tau }={{O}}(\epsilon ^3)$ discussed in the present study. The vector 
$\boldsymbol {E}$ is thus a higher-order term that is negligible at the level of the QG approximation. In other words, the evolution equations for the TWA and the fixed-
$z$ average are identical at the level of the QG approximation.
B.2. The coefficient $K_{GM}$ of the present study describes advection by the quasi-Stokes velocity in the QG limit
In the general TWA formulation of McDougall & McIntosh (Reference McDougall and McIntosh2001), the fluxes arising from the antisymmetric part of the diffusion tensor are expressed in terms of a quasi-Stokesstreamfunction 
$\boldsymbol {\psi }$ as
\begin{equation} \left[\begin{matrix} 0 & -\boldsymbol{\psi} \boldsymbol{\cdot} \boldsymbol{e}_y \\ \boldsymbol{\psi} \boldsymbol{\cdot} \boldsymbol{e}_y & 0 \\ \end{matrix} \right] \left( \begin{matrix} G_y \\ G_z \end{matrix} \right) , \end{equation}
where we restrict attention to the case of zonally invariant statistics. With the QG scalings (2.1), the 
$y$ component of the quasi-Stokes streamfunction provided in McDougall & McIntosh (Reference McDougall and McIntosh2001) reduces to 
$\boldsymbol {\psi }\boldsymbol {\cdot } \boldsymbol {e}_y=-{\overline {vb}}/{N^2} + {{O}}(\epsilon ^3)$. Using our definition in (4.1) for 
$K_{GM}$, we can re-express the right-hand side as 
$- K_{GM} {\mathcal {S}} + {{O}}(\epsilon ^3)$. This shows that (B1) is indeed equivalent to the GM part of the tensor (4.6) of the present study at the QG level of approximation.
Appendix C. Direct numerical simulation
The numerical simulations are performed using an intermediate set of equations between the Boussinesq equations and the QG system. Indeed, on the one hand, the QG system – (2.3) with the boundary conditions (2.5)–(2.6) – is rather impractical for implementation in standard pseudo-spectral solvers. On the other hand, going back to the full primitive equations is also impractical because the latter are incompatible with periodic boundary conditions in the horizontal directions: they involve terms that are proportional to 
$y$, thus breaking the invariance to translations in 
$y$. Fortunately these terms are subdominant in the QG range of parameters, and a convenient way to simulate the QG dynamics of a patch of ocean consists in discarding them from the set of Boussinesq equations. There are two such terms. First, the base state has a 
$z$-dependent meridional buoyancy gradient, and the vertical advection of the associated buoyancy field leads to coefficients that are proportional to 
$y$. We neglect this subdominant vertical advection of the background meridional buoyancy gradient. Second, the planetary vorticity gradient leads to a stretching term of the form 
$(f_0+\beta y)\,\partial _z \boldsymbol {u}$ in the vorticity equation. We neglect the subdominant contribution 
$\beta y\, \partial _z \boldsymbol {u}$ in the following. Specifically, we end up with the following set of dimensionless primitive-like equations for the departure fields:
\begin{align} &\partial_t \boldsymbol{u} + U'(z)\,w\,\boldsymbol{e}_x + U(z)\,\partial_x \boldsymbol{u} + (\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u} + \boldsymbol{e}_z \times \boldsymbol{u} -\beta [\psi \boldsymbol{e}_y + \varDelta_{{\perp}}^{{-}1} \{ \psi_{yz} \}\,\boldsymbol{e}_z] \nonumber\\ &\quad ={-} \boldsymbol{\nabla} p + b \boldsymbol{e}_z + E_\perp \varDelta_\perp \boldsymbol{u} + E_z\,\partial_{zz} \boldsymbol{u} , \end{align}
\begin{align} &\partial_t b + U(z)\,b_x - U'(z)\,v + N^2(z)\,w + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} b = E_{b,\perp} \varDelta_\perp b + E_{b,z}\,\partial_{zz} b , \end{align}
where 
$\boldsymbol {u}=(u,v,w)$ now denotes the velocity departure from the base state. The toroidal streamfunction 
$\psi (x,y,z,t)$ has a vanishing horizontal area average at every depth 
$z$ and is defined as 
$\psi =\varDelta _\perp ^{-1} \{ \partial _y u - \partial _x v \}$. For rapid rotation and strong stratification, the set (C1-C2) reduces to the limiting QG system ((2.3) with the boundary conditions (2.5)–(2.6)). Additionally, the form of the 
$\beta$ term ensures conservation of mechanical energy.
We solve (C1)–(C2) inside a horizontally periodic domain with the pseudo-spectral solver Coral (Miquel Reference Miquel2021), used previously for the Eady model (Gallet et al. Reference Gallet, Miquel, Hadjerci, Burns, Flierl and Ferrari2022) and for turbulent convective flows (Miquel et al. Reference Miquel, Bouillaut, Aumaître and Gallet2020; Bouillaut et al. Reference Bouillaut, Miquel, Julien, Aumaître and Gallet2021), and validated against both analytical results (Miquel et al. Reference Miquel, Lepot, Bouillaut and Gallet2019) and solutions computed with the Dedalus software (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). The background stratification is strong, and the global rotation is fast (low Rossby number), to ensure a strongly QG regime. We use insulated boundary conditions at top and bottom for 
$b$, free-slip boundary conditions at the surface for the velocity departure, and a frictional boundary condition at the bottom, 
$\partial _z u|_{-1}=-(\tilde {\kappa }/E_z) u|_{-1}$, 
$\partial _z v|_{-1}=-(\tilde {\kappa }/E_z) v|_{-1}$ and 
$w|_{-1}=0$. Such parametrized bottom drag is detailed in the study of the Eady model (Gallet et al. Reference Gallet, Miquel, Hadjerci, Burns, Flierl and Ferrari2022) together with the connection between the coefficient 
$\tilde {\kappa }$ and the QG friction coefficient 
$\kappa$ arising in (2.6). The dissipative coefficients in the simulation have values 
$\tilde {\kappa }=4.5\times 10^{-4}$, 
$E_z=3\times 10^{-6}$, 
$E_\perp =0.003$, 
$E_{b,z}=3 \times 10^{-7}$ and 
$E_{b,\perp }=0.001$.
