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:

(2.1)\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

(2.2)\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

(2.3)\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

(2.4)\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

(2.5)\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

(2.6)\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

(3.1)\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

(3.2)\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

(3.3)\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

(4.1a,b)\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

(4.2)\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

(4.3)\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

(4.4)\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:

(4.5a,b)\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

(4.6)\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

(5.1)\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:

(5.2)\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,

(5.3)\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

(5.4)\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.