Zonal jets and eddy tilts in barotropic geostrophic turbulence

Abstract

The spontaneous formation of zonal jets is a distinctive feature of geostrophic turbulence with the phenomenon witnessed in numerous numerical studies. In such systems, strong rotation anisotropises the spectral evolution of the energy density such that zonal modes are favoured. In physical space, this manifests as eddies zonally elongating and forming into zonal jets. In the presence of large scale dissipation, the flow may reach statistical stationarity such that the zonal structure persists in the zonal and time mean, and is supported by a flux of eddy momentum. What is unclear is how the excitation of Rossby waves arranges the underlying eddy momentum stresses to support the mean flow structures. To study this, we examine a steady-state flow in the so-called ‘zonostrophic’ regime, in which characteristic scales of geostrophic turbulence are well separated and there are several alternating zonal jets that have formed spontaneously. We apply a geometric eddy ellipse formulation, in which momentum fluxes are cast as ellipses that encode information about the magnitude and direction of flux; the latter is described using the tilt angle. With the aid of a zonal filter, it is revealed that the scales responsible for providing the momentum fluxes associated with the jet structure are much smaller than the characteristic scales identified, and occupy a region of the energy spectrum that has been typically associated with isotropic dynamics.

1. Introduction

Spontaneous zonal jet formation is observed all over the world's oceans (Treguier et al. 2003; Galperin et al. 2004; Maximenko, Bang & Sasaki 2005) and is often attributed to geostrophic turbulence. This is a fluid regime that has a high Reynolds number and is in approximate geostrophic balance, which is a fitting approximation for the mesoscale ocean. Such a fluid may exhibit the inverse cascade of energy (Kraichnan 1967; Batchelor 1969) in which energy is transferred successively upscale by eddies of larger size. Under strong rotation, Rossby waves become excited, and the transfer upscale is anisotropised favouring zonal modes, resulting in the formation of zonal jets (Rhines 1975; Holloway & Hendershott 1977). Spontaneous zonal jet formation is also observed in outputs from eddy-resolving ocean models (Treguier & Panetta 1994; Nakano & Hasumi 2005). The signature banded patterns of the Solar System's gas giants may also be attributed to geostrophic turbulence (Williams 1978). Indeed, satellite altimetry data revealed the inverse cascading process in the South Pacific Ocean (Scott & Wang 2005), and analysis of data from Cassini's Jupiter flyby demonstrates the existence of geostrophic turbulence in the Jovian atmosphere (Choi & Showman 2011; Galperin et al. 2014), as well as the forward and inverse cascade (Young & Read 2017).

Since mesoscale ocean eddies are not currently resolved in global climate and Earth system models used for long-range climate projections, accurate parametrisations are crucial for capturing their impact on larger-scale processes. The Gent–McWilliams scheme (Gent & McWilliams 1990), developed to represent vertical eddy momentum fluxes, has become a staple in current global ocean and climate models, effectively capturing the net transport across the Antarctic Circumpolar Current. However, finding a suitable parametrisation of horizontal momentum fluxes in such models, which incorporates the crucial nonlinear effects of the inverse cascade, remains a challenge (Sukoriansky, Galperin & Chekhlov 1999; Eden 2010; Marshall, Maddison & Berloff 2012).

Eddies of larger scale in the inverse cascade are also the longest lived. Rhines (1975) theorised that the anisotropisation of the cascade under strong rotation leads to these eddies zonally elongating into coherent zonal jets with mode number given by the Rhines scale:

(1.1)\begin{equation} k_{jet} = k_R \equiv \sqrt{\frac{\beta}{2U}}, \end{equation}

where $k$ is the inverse scale, $k_{R}$ and $k_{jet}$ are the Rhines and jet scales, respectively, $\beta$ is the northward gradient of the Coriolis frequency, and $U$ is some characteristic velocity that we take to be the root mean square flow velocity.

Jets that form in $\beta$-plane turbulence have been found to possess an asymmetric zonal structure, with sharp eastward jets interleaved by broad westward jets. This is to maintain barotropic stability (Vallis & Maltrud 1993; Danilov & Gurarie 2004; Scott & Dritschel 2019). Dritschel & McIntyre (2008) discussed how this jet structure leads to regions of homogenised potential vorticity separated by sharp gradients that in some instances appear in a ‘staircase’ pattern (Scott & Dritschel 2012). In practice, however, jet meandering reduces the sharpness of the potential vorticity gradient (Srinivasan & Young 2012). Eddy–eddy interactions are necessary for stirring potential vorticity in this way (Scott & Dritschel 2019). Zonostrophic turbulent jets and the eddy–mean-flow interactions that give rise to them will be the focus of this paper.

Motivated by Rhines’ seminal paper, several subsequent numerical studies have identified other characteristic scales and spectra that are important in geostrophic turbulent jet formation. In the anisotropisation process, energy is funnelled into zonal modes, causing a steep spectrum ($k^{-5}$) to develop near the meridional wavenumber axis, referred to here as the zonal spectrum (Chekhlov et al. 1996). This spectrum has been witnessed in numerous numerical studies (Chekhlov et al. 1996; Galperin, Sukoriansky & Huang 2001; Sukoriansky, Galperin & Dikovskaya 2002; Galperin et al. 2006, 2019; Galperin, Sukoriansky & Dikovskaya 2008). In addition to the Rhines scale, another characteristic scale was identified, namely $k_\beta$, the ‘cross-over’ wavenumber (Pelinovsky 1978; Maltrud & Vallis 1991; Vallis & Maltrud 1993), which was later found to be the wavenumber at which the flow anisotropises (Sukoriansky, Dikovskaya & Galperin 2007). Galperin et al. (2006, 2008) and Sukoriansky et al. (2007) argued that between these two characteristic scales lies a new flow regime called ‘zonostrophic turbulence’. It is in this flow regime that jets are found. Sukoriansky et al. (2002) found that spectra of gas giants agree with the steep zonal spectrum, and data from the Cassini mission suggest that the circulation of Jupiter is in agreement with this regime (Galperin et al. 2014). Recently, the zonostrophic regime has also been observed in laboratory experiments (Lemasquerier, Favier & Le Bars 2021, 2023). In this paper, we will examine a flow in which these characteristic scales hold, and verify these scaling relationships.

Though there is no obvious demarcation between the underlying eddy fields and the coherent structures into which they develop, it is useful to examine how eddies interact with the mean flow. A simple method of examining this interaction is to perform Reynolds decomposition (Reynolds 1895), separating the dynamical equations into their mean components and departures from the mean. If flows reach a steady state, then these Reynolds-averaged equations demonstrate a balance between the regions of converging eddy momentum fluxes and the mean flow structures that they support (Starr 1968). Efforts to parametrise horizontal momentum fluxes have followed the residual-mean approach by describing eddy fluxes of potential vorticity as a divergence of the eddy stress tensor (Marshall et al. 2012). One way to visualise momentum fluxes is to cast eddy velocity correlations as a variance ellipse (Preisendorfer 1988). Visualising momentum fluxes through the variance ellipses of eddy velocity correlations provides valuable insights into the strength and direction of the underlying eddy field's momentum fluxes. This approach has been applied successfully to examine eddy–mean-flow interactions in western boundary jets (Waterman & Hoskins 2013) and barotropically unstable jets (Tamarin et al. 2016). In this paper, we demonstrate how this formulation may be used to study jets in a fully turbulent system. Huang & Robinson (1998) followed a similar approach for studying eddy–mean-flow interactions in turbulent zonal jets in a global barotropic model, namely a transient–stationary decomposition (McWilliams 1984). Using calculations of spectral energy fluxes, Huang & Robinson (1998) made the intriguing observation that eddies at scales close to the forcing scale appear to support the jet structures through local shearing of eddies, which results in energy being transferred directly to the jet scale in a non-local transfer mechanism that bypasses the inverse cascade (Srinivasan & Young 2012; Galperin et al. 2019). Our present work offers a unique perspective by employing the geometric eddy ellipse framework. This framework not only allows us to identify the scales responsible for supporting the zonal jets but also provides further insights into the spatial arrangement of underlying eddy momentum fluxes that contribute to the formation and maintenance of these jets.

In this paper, we will investigate a flow within the zonostrophic regime, examining the characteristic scales exhibited by the system, and discuss how these relate to other perspectives of spontaneous jet formation in the literature. Furthermore, the study will explore the eddy–mean-flow interactions that support the formation of zonal jets using the geometric eddy ellipse formulation. By examining the interplay between eddy momentum fluxes and the mean flow, we gain a deeper understanding of the scales and spatial organisation of eddy tilt patterns in zonostrophic jets.

2. Characteristic scales in jets

In addition to the inverse cascade of energy, Rhines theorised the existence of a steep $k^{-5}$ power law. This was found in Chekhlov et al. (1996) to manifest in the vicinity of the meridional wavenumber axes ($k_y$) (Huang, Galperin & Sukoriansky 2001; Sukoriansky et al. 2002). As the energy along $k_y$ follows different scaling laws to the rest of the spectrum, Sukoriansky et al. (2007) divided the total energy spectrum into zonal $E_Z(k)$ and residual $E_R(k)$ (the inverse cascade) components:

(2.1)

where $E(k)$ is the energy density as function of wavenumber $k$, $\epsilon$ is the energy injection rate, and $C_{\beta }$ and $C_K$ are the dimensionless Rhines and Kolmogorov constants, empirically determined to be $C_{\beta } \sim 0.5$ and $C_K \sim 5\unicode{x2013}6$ (Sukoriansky et al. 2002). Here, $\phi = \arctan (-k_x/k_y)$ is the angle measured anticlockwise from the $k_y$-axis, and $k_x$ is the zonal wavenumber. We chose a threshold $\phi = {\rm \pi}/8$ as the size of the wedge around the $k_y$-axis that contains zonal energy. By equating (2.1a) and (2.1b), Sukoriansky et al. (2007) found the transition point in the total spectrum, i.e. the wavenumber at which the energy contained within zonal modes exceeds that of the residual modes:

(2.2)\begin{equation} k_\beta = \left(\frac{C_\beta}{C_K}\right)^{{3}/{10}} \left(\frac{\beta^3}{\epsilon}\right)^{1/5} \approx 0.5\left(\frac{\beta^3}{\epsilon}\right)^{1/5}. \end{equation}

This cross-over wavenumber (Pelinovsky 1978; Maltrud & Vallis 1991; Vallis & Maltrud 1993) marks the point at which the spectrum anisotropises in $\beta$-plane turbulence. Note that this is not to say that for length scales where $k>k_\beta$ the system is isotropic. Instead, we expect energy near $k_x=0$ (the meridional axis) to dominate for length scales exceeding $k_\beta$.

After $k_{\beta }$ is crossed, most of the energy in the system resides in the zonal modes (2.1). By integrating this over wavenumber space, Sukoriansky et al. (2007) obtained the total energy of the system, and found that most of the energy resides in the largest occupied scale on the spectrum, namely the Rhines scale given by (1.1). The size of the jets scales according to the largest, most energetic scale, confirming the Rhines (1975) theory. When there is a sufficiently strong large-scale dissipation such that the energy in the system reaches equilibrium before the flow approaches the domain scale, the Rhines scale settles at the equilibrium wavenumber given by

(2.3)\begin{equation} k_{jet} = k_R = \left(\frac{\beta^2}{8E_{eq}}\right)^{{1}/{4}}, \end{equation}

where $E_{eq}$ is the energy at equilibrium (Danilov & Gurarie 2002; Sukoriansky et al. 2007). For a system to exhibit jet formation, Sukoriansky et al. (2007) argued using a range of numerical experiments that characteristic scales must be well separated as follows:

(2.4)\begin{equation} k_0 \ll k_R \ll k_\beta \ll k_f \ll k_D, \end{equation}

where $k_0$ is the largest allowable scale, $k_f$ is the forcing scale, and energy dissipation occurs at small scales with wavenumber greater than $k_D$. Thus the Rhines scale must be much larger than the cross-over wavenumber to allow significant anisotropy to develop, and must also be much smaller than the domain scale so that domain effects do not interfere. Furthermore, the anisotropic wavenumber must not be close to the forcing scale, so that the inverse cascade may develop. The last inequality arises from classic turbulence arguments (e.g. Kraichnan 1967) that the small-scale dissipation must not be significant in the inertial range. Flows obeying (2.4), particularly where $k_\beta /k_R \gtrsim 2$, are described by Galperin et al. (2008) to be in the zonostrophic regime. Independently, Scott & Dritschel (2012) found that flows where $k_\beta /k_R \gtrsim 4$ demonstrated sharp eastward and broad westward jet structures. The greater this ratio, the more the jet structures approached the idealised potential vorticity staircase structure (Dritschel & McIntyre 2008).

3. Geometric eddy ellipse formulation

The Reynolds decomposed shallow-water momentum equation is given by

(3.1)\begin{equation} \frac{\partial \bar{\boldsymbol{u}}}{\partial t} + \overline{f\boldsymbol{k}\times\boldsymbol{u}} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}} - g\,\boldsymbol{\nabla}\bar{\eta} + \bar{\boldsymbol{F}} + \bar{\boldsymbol{D}}, \end{equation}

where $\boldsymbol {u}\equiv (u,v)$ is the horizontal velocity, with zonal and meridional velocity components denoted $u$ and $v$, respectively, $f = f_0+\beta y$ is the Coriolis frequency under the $\beta$-plane approximation, $\boldsymbol {k}$ is the vertical unit vector, $g$ is the vertical component of the acceleration due to gravity, $\eta$ is the free surface of the fluid, $\boldsymbol {\nabla }$ is the horizontal gradient operator, and $\boldsymbol {F}$ and $\boldsymbol {D}$ are the horizontal forcing and dissipation mechanisms, respectively. We choose here to define the ‘eddying’ portion of the flow as fluctuations with respect to the zonal mean since jets formed in geostrophic turbulence are statistically zonal. The zonal mean is denoted with an overline, and primes indicate fluctuations. This is a crude assumption because in this fully developed turbulent system, there are many interacting scales of motion, so it is difficult to separate exactly what constitutes transient eddies and the long-lived zonal structure. This is something that we will explore. The advective terms vanish upon zonal averaging. We find that eddies influence the mean flow through the divergence of the eddy velocity correlation tensor:

(3.2)\begin{equation} \boldsymbol{\mathsf{T}} = \boldsymbol{\mathsf{T}}(y) = \begin{pmatrix} \overline{u'u'} & \overline{u'v'}\\ \overline{u'v'} & \overline{v'v'} \end{pmatrix}. \end{equation}

Equation (3.2) may be separated into isotropic (trace-only) and anisotropic (trace-free) parts:

(3.3)\begin{equation} \boldsymbol{\mathsf{T}} = K\boldsymbol{\mathsf{I}} + \boldsymbol{\mathsf{E}}, \end{equation}

where $\boldsymbol{\mathsf{I}}$ is the identity matrix, and the zonally averaged kinetic energy density associated with the eddying portion of the flow field is given by

(3.4)\begin{equation} K = K(y) = \tfrac{1}{2}(\overline{u'^2} + \overline{v'^2}), \end{equation}

and we have introduced the eddy momentum stress tensor

(3.5)\begin{equation} \boldsymbol{\mathsf{E}} = \begin{pmatrix} M & N\\ N & -M \end{pmatrix}. \end{equation}

The normal stress difference is given by

(3.6)\begin{equation} M = M(y) = \tfrac{1}{2}(\overline{u'^2} - \overline{v'^2}), \end{equation}

and the shear stress is given by

(3.7)\begin{equation} N = N(y) = \overline{u'v'}. \end{equation}

In a quasi-geostrophic, incompressible fluid system, $\bar {v} = 0$. Thus the dynamical evolution of (3.1) is governed by the divergence of the sheer stress:

(3.8)\begin{equation} \frac{\partial \bar{u}}{\partial t} ={-}\frac{\partial N}{\partial y} \end{equation}

in the absence of forcing and dissipation. The Reynolds stresses (3.6) and (3.7) describe the single-point correlations between eddy velocities $u'$ and $v'$. In the geometric eddy ellipse formulation, if there are strong correlations, then the mean pattern traced in $u'\unicode{x2013}v'$ space will be elliptical, and if no correlations exist, then the trace will be circular. We can recover these ellipse properties from the covariance tensor given by (3.2) by performing a principal axes decomposition (Preisendorfer 1988), diagonalising (3.5) under a rotation through angle $\theta$, following Waterman & Lilly (2015). This is equivalent to rotating the ellipse such that its semi-major axis aligns with the zonal direction. This allows us to rewrite (3.5) as

(3.9)\begin{equation} \boldsymbol{\mathsf{E}} = L\begin{pmatrix} \cos{2\theta} & \sin{2\theta}\\ \sin{2\theta} & -\cos{2\theta} \end{pmatrix}, \end{equation}

where we have defined

(3.10)\begin{equation} L = L(y) = \sqrt{M^2+N^2}. \end{equation}

Here, $L$ is the excess energy in the direction of the major axis of the ellipse compared to the minor axis, and $\theta$ is given by

(3.11)\begin{equation} \theta = \theta(y) = \frac{1}{2}\arctan{\left(\frac{N}{M}\right)}, \end{equation}

which is the tilt of the ellipse with respect to the zonal direction. The geometric eddy ellipse is not a description of the shape of individual eddies themselves, but a description of their net fluxes. Through examining the distribution of eddy ellipses and observing their patterns, we can determine the direction and magnitude of these fluxes (Waterman & Hoskins 2013). Flows may demonstrate significant anisotropy between the zonal and meridional components of their flow velocities, such that $M\neq 0$, but this may not necessarily give rise to a zonal momentum tendency. This kind of behaviour corresponds to neutral ellipse tilts, where $N=0$. Conversely, when $M=0$, flow velocities do not demonstrate any significant anisotropy between the zonal and meridional components, but there may be a zonal momentum tendency where $N\neq 0$. An eddy ellipse tilted at $\theta = \pm ({{\rm \pi} }/{4})$ is achieved in this instance.

By examining the eddy ellipse shapes, we can tell which eddy distributions give rise to momentum fluxes into the mean flow, and which do not. An example of an eddy shape that demonstrates a zonal momentum flux is the banana-shaped eddy depicted by Wardle & Marshall (2000) (see their figure 5), where the eddy shape is such that momentum is fluxed eastwards. Examining the distribution of the geometric ellipses visualises these flux directions more intuitively than examining the velocities. Circular eddy ellipses are achieved only when $M=0$ and $N=0$, and there is no significant anisotropy, either through correlations between $u'$ and $v'$ or where $\overline {u'u'}\neq \overline {v'v'}$. In this scenario, only the isotropic component of the tensor given by (3.3) remains. Examples of other eddy ellipse distributions that give rise to zonal momentum fluxes can be found in figures 3 and 4 of Tamarin et al. (2016). In this paper, we will evolve our flow using the barotropic vorticity equation given by

(3.12)\begin{equation} \frac{\partial\zeta}{\partial t} ={-} u\,\frac{\partial \zeta}{\partial x} - v\, \frac{\partial \zeta}{\partial y} - \beta v + F - r\zeta - \mathcal{D}\,\nabla^4\zeta, \end{equation}

where $r$ and $\mathcal {D}$ are the Rayleigh friction and small-scale biharmonic dissipation parameters, respectively. Here, the zonal and meridional velocity components are given by

(3.13)\begin{equation} (u,v) \equiv \left(-\frac{\partial\psi}{\partial y},\frac{\partial\psi}{\partial x}\right), \end{equation}

written here in terms of the stream function $\psi$. The relative vorticity is given by

(3.14)\begin{equation} \zeta \equiv \left(\frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}\right) \equiv \nabla^2\psi, \end{equation}

and we introduce rotation into our system through the $\beta$-plane. Rayleigh friction is chosen so that the flow reaches statistical equilibrium before the energy approaches the largest scales, such that characteristic scales are well separated according to (2.4).

The Reynolds-averaged inviscid unforced form of the barotropic vorticity equation is

(3.15)\begin{equation} \frac{\partial \bar{\zeta}}{\partial t} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}(\overline{\boldsymbol{u}'\zeta'}), \end{equation}

where eddies influence the mean flow in the barotropic vorticity equation through the divergence of an eddy vorticity flux $\overline {\boldsymbol {u}'\zeta '}$. From this, we can relate (3.2) to the zonal mean eddy vorticity flux in (3.15) through the Taylor–Bretherton identity (Taylor 1915; Bretherton 1966; Plumb 1986):

(3.16)\begin{equation} \overline{\boldsymbol{u}'\zeta '} = \left(-\frac{\partial M}{\partial y}, -\frac{\partial N}{\partial y}\right). \end{equation}

Thus only the components of the eddy momentum stress tensor given by (3.5) appear explicitly in the eddy forcing of the mean barotropic vorticity equation.

4. Model

We will solve (3.12) numerically on a zonally periodic, meridionally bounded channel model with zonal and meridional dimensions $\pm L_x/2$ and $\pm L_y/2$, respectively, where $L_x = 6000\,\mathrm {km}$ and $L_y = 3000\,\mathrm {km}$. A laterally bounded model is chosen so that jets do not migrate meridionally. We must ensure that the appropriate boundary conditions are applied to maintain consistency between (3.12) and the shallow-water equations from which it is derived. We determine $\psi$ by inverting (3.14), and obtain $u,v$ using (3.13). There must be no flow through the meridional walls so $v=0$ such that $\psi '=0$. This imposes no conditions on the zonal mean flow on the boundary, so we are free to choose a mixture of boundary conditions: no-slip on the zonal mean modes, and free-slip boundary conditions on the non-mean modes (e.g. Esler 2008; Williams et al. 2009) given by

(4.1)\begin{equation} \left.\begin{gathered} \frac{\partial\bar{\zeta}}{\partial y} = \frac{\partial^3\bar{\zeta}}{\partial y^3} = 0,\\ \zeta' = \frac{\partial^2 \zeta'}{\partial y^2} = 0. \end{gathered}\right\} \end{equation}

Solving for $\psi$ will also produce $\bar {\psi } = \mathrm {constant}$ on each of the lateral boundaries. These constants are, in general, different. The first is determined by setting $\bar {u}=0$ on each of the lateral boundaries, and the second is found through mass conservation. Setting $\bar {u}=0$, together with the boundary conditions (4.1), trivially satisfies consistency conditions between the barotropic vorticity equation and the shallow-water equations through conserving circulation (McWilliams 1977; Graef & Müller 1996).

We discretise (3.12) on a channel model with $513\times 257$ grid points such that the zonal and meridional grid spacing is given by $\Delta x = \Delta y = 11.7\,\mathrm {km}$, and the time step is $\Delta t = 112.5\,\mathrm {s}$. The flow is spun up from rest and forced using a stochastic isotropic forcing function (e.g. Chekhlov et al. 1996) with a constant enstrophy injection rate $\xi = 8\times 10^{-18}\,\mathrm {s}^{-3}$ at the forcing scale given by $k_f = 64$, such that the constant energy injection rate is $\epsilon = 1.78\times 10^{-8}\,\mathrm {m}^2\,\mathrm {s}^{-3}$. The scale-selective biharmonic diffusion term is used to remove energy at small scales with a coefficient $\mathcal {D} = 5\times 10^{10}\,\mathrm {m}^4\,\mathrm {s}^{-1}$, and the flow is damped through Rayleigh friction with a coefficient $r = 1\times 10^{-8}\,\mathrm {s}^{-1}$. We choose $\beta = 4\beta _0$, where $\beta _0 = 2.29\times 10^{-11}\,\mathrm {m}^{-1}\,\mathrm {s}^{-1}$ is the terrestrial equatorial value. Though the $\beta$-value that we are using is much greater than the terrestrial value, the parameters chosen possess dynamic similarity to terrestrial problems. Note that all wavenumbers quoted have been scaled by $2{\rm \pi} /L_y$.

5. A zonostrophic jet

The parameters for this barotropic channel model have been chosen in such a way that the flow is in the zonostrophic regime and exhibits multiple zonal jets, obeying the scaling relationships presented in Sukoriansky et al. (2007). In this section, we will discuss these characteristic scales in detail, and show how they relate to the structure and stability of the jets.

We observe the evolution of the flow as it is forced isotropically and spun up from rest, eventually approaching its equilibrium jet configuration (figure 1). Early in the flow evolution, zonal structures emerge that merge together as time evolves. What we are essentially observing is the evolution of the Rhines scale as the system increases in energy and $U$ increases. Using (2.3), we calculate the final destination of the Rhines scale, the mode number $k_{eq} \approx 4$. This agrees well with the observed number of peak to peak jets, in which there are $5$ sharp eastward jets interleaved by $4$ broad westward jets. We also note a subtle interplay between the energy and jet stability. As noted in Vallis & Maltrud (1993), sharp eastward and broad westward jet profiles correspond to barotropically stable structures. Rearranging a jet configuration then requires breaking the stability, and we see a fast change in energy as this occurs. So as the flow progresses, there are a few metastable jet configurations that the flow assumes, before approaching its final configuration.

Figure 1. (a) Hovmöller diagram of $\bar {u}$ normalised by the instantaneous maximum zonal velocity. Superimposed is the evolution of the domain average energy as a function of time (red). (b) The zonal time average velocity profile, $[\bar {u}]$ (thick black) and the time average divergence of momentum flux, $-r^{-1}\,\partial _y N$ (thin orange). The vertical black line in (a) marks the beginning of the $20\,000\,\mathrm {h}$ period over which the time averages in (b) are taken.

The tendency in the zonal mean flow is governed by a divergence of the shear stress (3.8). Upon introducing Rayleigh friction, the flow will equilibrate as energy is removed from the largest scales. We introduce this in (5.2), by setting

(5.1)\begin{equation} \bar{\boldsymbol{D}} ={-}r\bar{\boldsymbol{u}}. \end{equation}

Since $\bar {v}=0$, the dynamics is governed solely by the zonal component of the shallow-water equation,

(5.2)\begin{equation} \frac{\partial\bar{u}}{\partial t} ={-} \frac{\partial N}{\partial y} -r\bar{u}. \end{equation}

When the flow reaches a steady state,

(5.3)\begin{equation} \bar{u} ={-} \frac{1}{r}\,\frac{\partial N}{\partial y}. \end{equation}

So the jets will be maintained by the divergence of the shear stresses. Of course, in our case we also have a biharmonic diffusion term, so this relationship is approximate. We examine this relationship in figure 1(b). The time average (denoted by square brackets) is taken over the last $20\,000\,\textrm {h}$ of the flow evolution. Over this period, the divergence in the shear stress broadly matches the zonal mean structure with some small variation. Shorter time intervals produce less agreement, with the flow demonstrating short-time-scale variations that produce a tendency in the mean jet profile such that the relationship given by (5.3) no longer holds. The fact that the divergence of the shear stress produces steady jet patterns after a long-time mean suggests that a scale separation may exist between scales producing the jet patterns and the scales producing temporal variations.

We now demonstrate that these jets are characteristic of those in the literature and lie within the zonostrophic regime. The energy spectrum $E(k_x,k_y)$ demonstrates a cascade of energy towards the largest scales (figure 2a), and shows the characteristic ‘dumbbell’ shape that was first described in Maltrud & Vallis (1991). This is associated with the anisotropisation of the flow due to the presence of Rossby waves. A large amount of energy is seen to have built up along the $k_y$-axis, which can be understood by examining the angularly averaged spectra (figure 2b).

Figure 2. (a) The time average energy spectrum $\log E(k_x,k_y)$. (b) The time and angular averaged spectra: total $E(k)$ (solid black), zonal $E_Z(k)$ (solid orange), and residual $E_R(k)$ (turquoise). The vertical dashed lines from right to left are the forcing scale $k_f$ (plum red), cross-over wavenumber $k_\beta$ (blue), and Rhines scale $k_R$ (green). The jet profile spectrum is marked with a thin yellow line. The slopes of the theoretical residual spectrum ($\propto k^{5/3}$, dashed red) and zonal spectrum ($\propto k^{-5}$, dot-dashed red) are also given. The time average is taken over a $20\,000\,\mathrm {h}$ equilibrium period.

The total spectrum (figure 2b), initially traces the Kolmogorov slope given by $k^{-5/3}$, associated with isotropic dynamics. We calculate the cross-over wavenumber given by (2.2): $k_\beta \approx 16$. We find that as the cascade approaches $k_\beta$, $E(k)$ steepens and possesses several sharp spikes of energy that trace the $k^{-5}$ slope, confirming this as the approximate point at which the spectrum anisotropises. This corresponds to the energy build up along the $k_y$-axis in $E(k_x,k_y)$ (figure 2a).

In splitting the spectrum into the zonal and residual spectra, we find that at smaller scales, most of the energy is concentrated in the residual spectrum, i.e. the region of $E(k_x,k_y)$ that lies outside of the $k_y$-axis, and traces the $k^{-5/3}$ slope. The zonal spectrum also initially traces the $k^{-5/3}$ slope, which is thought to be indicative of the absence of Rossby waves at these scales. The zonal spectrum steepens and surpasses the residual spectrum at $k \approx k_\beta$. At the largest scales, the energy contained in the zonal modes dwarfs that of the residual such that the amplitude of the zonal spectrum is approximately equal to that of the total spectrum (figure 2a). This possesses sharp, discrete spikes of energy that are harmonics of the largest energy-containing scales, i.e. mode numbers $4$ and $5$, corresponding to the number of westward and eastward jets, respectively (Danilov & Gurarie 2004).

Examining figure 2(b), we see that all characteristic scales of interest (i.e. $k_R\approx 4$, $k_\beta \approx 16$ and $k_f=64$) are well separated according to (2.4). However, we note that the biharmonic diffusion causes a strong suppression of the enstrophy cascade throughout the inertial range, so the last inequality in (2.4) is not obeyed by this system.

We also examine the jet profile spectrum, which is the spectrum corresponding to the zonal mean jet structure (figure 1b). Here, we see that the spectrum associated with the jet profile possesses the discrete harmonics observed in the zonal and total energy spectrum (figure 2b). The asymmetric sharp eastward and broad westward jets are required for the jet profile to be barotropically stable (Vallis & Maltrud 1993); we argue that these discrete harmonics are then a necessary consequence of barotropic stability.

6. Eddy tilts

In the previous section, we have seen how the parameters chosen have led to jets forming in the zonostrophic regime, with characteristic scales adequately well separated according to (2.4) and an equilibrium jet profile consisting of sharp eastward and broad westward jets representing a barotropically stable structure. We will now turn our attention to examining eddy–mean-flow interactions using the geometric eddy ellipse formulation, attempting to identify the scales responsible for fluxing momentum into the jet structure, balancing the mean flow.

Tamarin et al. (2016) examined a barotropically unstable zonal jet profile that successively strengthened or weakened as instabilities on the jet flanks decayed and grew, respectively. Here, the jet structures that form as a result of spectral anisotropisation by the $\beta$ effect have barotropically stable jet profiles that satisfy the Rayleigh–Kuo stability criterion. We expect the eddy ellipse patterns corresponding to the jets to persist in time in order to maintain the jet structure. It is peculiar that we observe the relationship given by (5.3) only after a relatively long time averaging, when the jet structures themselves are persistent and stable. To understand this, we examine the Hovmöller plots of the momentum fluxes $N, M$ and the ellipse quantity $\theta$ in figure 3 over a $20\,000\,\textrm {h}$ period in which the energy is in a statistical equilibrium. Though the energy has reached equilibrium, we see that there are significant time-dependent processes occurring in the shear stress $N$, where momentum is fluxed north and south in an alternating pattern (figure 3a). The tendency in $\bar {u}$ is governed by (5.2), so these time-dependent processes result in short-term fluctuations in jet intensity. We see, however, that over a time average, the shear stress possesses a pattern that is consistent with the jet structure (figure 3b). In particular, we find that eastward jets are flanked by a negative flux of eastward momentum to their north and a positive flux of eastward momentum to their south. The converse is true for westward jets. This leads to a characteristic pattern in the stress that is half a jet width offset from the jet structure.

Figure 3. Hovmöller diagrams of: (a) the shear stress $N(y)$ and (b) its normalised time average $[N(y)]/10^{-3}$ (blue solid); (c) the normal stress difference $M(y)$ and (d) its normalised time average $[M(y)]/10^{-2}$ (blue solid); and (e) the corresponding ellipse tilt angle $\theta (y)$ and (f) its time average $[\theta ]$ normalised by ${\rm \pi} /2$ (blue solid). The normalised zonal-time average zonal velocity $[\bar {u}]/U_{max}$ is also plotted for comparison (black dashed), where $U_{max}$ is the maximum value of $[\bar {u}]$. The time averages are taken over a $20\,000\,\mathrm {h}$ equilibrium period.

The difference in normal stresses $M$ is in general an order of magnitude larger than $N$ (figure 3c). The temporal evolution of $M$ shows a zonal structure that is time-dependent and pulsing regularly. The pulses are always positive and coincide with regions where the $N$ changes sign in time. These pulses are located at some of the jet cores and at domain walls, but do not appear to correlate with any strengthening pattern in the zonal mean structure itself. We speculate that the intermittent signal of pulses in $M$ can be attributed to the action of the Rayleigh friction, which becomes important as the jet strengthens. In both the time evolution and the time average, $M$ is always non-negative, such that $\overline {u'^2}> \overline {v'^2}$. Also, although there is a zonal jet pattern in the normal stress field, this structure does not match the zonal mean jet pattern in $\bar {u}$ (figure 3d).

The resulting Hovmöller plot of the eddy tilt angle $\theta$, calculated from $M$ and $N$ using (3.11), does not look particularly instructive (figure 3e). It demonstrates a time alternating pattern of net northward and net southward momentum fluxing which possesses a meridional structure that is strongest at the jet cores, though it is not obvious how these fluxes interact with the jet. However, the time average shows a very clear eddy–mean-flow interaction, consistent with the jet pattern, where eddy ellipses are tilted towards the jet structure in an alternating pattern on the flanks of the jets (figure 3f). We see that there are time-dependent momentum fluxes that dominate the flow as discovered in Huang & Robinson (1998), which mask a more subtle momentum flux that supports the jet structure. The flow converges to this structure only approximately over a 20 000 h period. Our task is therefore to find exactly what scales are responsible for supporting the jet.

6.1. Zonal filter

In the previous section, we calculated a number of characteristic scales. We found that scales where $k > k_\beta$ are dominated by energy in the residual spectrum (2.1b), and scales where $k < k_\beta$ are dominated by energy accumulation along the $k_y$-axis corresponding to the development of the zonal spectrum (2.1a). The latter we associate with the formation of jets. We also know that in the steady state, the number of jets is equal to the Rhines scale $k_R$. Consequently, we might expect that eddies supporting the jet structures will reside between the Rhines scale $k_R$ and the anisotropy scale $k_\beta$. Since the flow favours zonal modes in the zonostrophic regime, we will develop a filter that separates the velocity correlations that contribute to the eddy stresses $M$ and $N$ by zonal scale.

We can exploit the homogeneity in $x$ to find the contributions of eddy velocities of different zonal scales to the eddy stress quantities $M$ and $N$. Since our flow is periodic in $x$, we take a Fourier transform of the eddy velocities in the $x$-direction, and substitute these into the eddy velocity correlation tensor (3.2). Thus

(6.1) \begin{equation} \boldsymbol{\mathsf{T}}(y) = \frac{2}{L_x}\begin{pmatrix} \displaystyle\int_1^{k_D} |\widetilde{u'}|^2 \,{\rm d} k_x & \displaystyle\int_1^{k_D} \widetilde{u'}\widetilde{v'}^\ast \,{\rm d} k_x\\ \displaystyle\int_1^{k_D} \widetilde{v'}\widetilde{u'}^\ast \,{\rm d} k_x & \displaystyle\int_1^{k_D} |\widetilde{v'}|^2 \,{\rm d} k_x \end{pmatrix}, \end{equation}

where the tilde represents the Fourier transform in the zonal direction, the asterisk is the complex conjugate, and we take the limit of integration from $k_x=1$ as the eddying velocities denoted by a prime have the zonal mean subtracted. From this, we can rewrite $M$ as

(6.2a)$$\begin{gather} M = M_l + M_h, \end{gather}$$

(6.2b)$$\begin{gather}\text{where}\quad M_l =\frac{1}{L_x}\int_1^{k_c-1} \left(\widetilde{|u'|}^2 - \widetilde{|v'|}^2\right){\rm d} k_x \end{gather}$$

(6.2c)$$\begin{gather}\text{and}\quad M_h = \frac{1}{L_x}\int_{k_c}^{k_D} \left(\widetilde{|u'|}^2 - \widetilde{|v'|}^2\right){\rm d} k_x, \end{gather}$$

and $N$ as

(6.3a)$$\begin{gather} N = N_l + N_h, \end{gather}$$

(6.3b)$$\begin{gather}\text{where}\quad N_l = \frac{1}{L_x}\int_0^{k_c-1} \widetilde{u'}\widetilde{v'}^\ast \,{\rm d} k_x \end{gather}$$

(6.3c)$$\begin{gather}\text{and}\quad N_h = \frac{1}{L_x}\int_{k_c}^{k_D} \widetilde{u'}\widetilde{v'}^\ast \,{\rm d} k_x, \end{gather}$$

where $k_c$ is some threshold wavenumber, and the subscripts $l$ and $h$ represent low-pass and high-pass filters, respectively.

We examine the scales associated with the eddy–mean-flow relationship given by (5.3), separating $N$ into low-pass ($N_l$) and high-pass ($N_h$) filtered signals as in (6.3a), for different cut-off wavenumbers $k_c$ (figure 4). The low-pass filter on $N$ is only able to capture the jet pattern in its divergence for large values of $k_c$, with the structure beginning to emerge at relatively high wavenumbers $k_c \geq 32$ (figure 4a). Below this, the signal has an amplitude comparable to the jet scale, but the pattern does not match the jet profile. In contrast, the jet pattern is evident for nearly all $k_c$, in the high-pass filter, with a faint, low-amplitude jet pattern still present when wavenumbers $k\geq k_f$ are retained (figure 4b). As $k_c$ varies, the transition between the noisy low wavenumbers and the appearance of the jet pattern in the high wavenumbers is relatively smooth. Though we cannot pinpoint a sharp transition scale, what is clear is that the divergence of $N$ does not capture the jet pattern for scales where $k_c \leq k_\beta$. This is counterintuitive because this is the scale associated with spectral anisotropy; when $k>k_\beta$, it is thought that isotropic dynamics dominate. Choosing $k_c=k_\beta$ does appear to be somewhat representative of the transition between the two regimes of eddy–mean-flow interactions, so we will proceed using this choice.

Figure 4. The time-average zonal velocity profile $[\bar {u}]$ over the $20\,000\,\mathrm {h}$ equilibrium period (black dashed) compared to (a) low-pass filtered divergence of the shear stress, $-r^{-1}\partial _yN_l$, and (b) high-pass filtered divergence of the shear stress, $-r^{-1}\partial _yN_h$. Here, $N_l$ is the low-pass filter on $N$, retaining zonal wavenumbers from $k_x=1$ up to cut-off wavenumber $k_c = 4,8,16,32,64,128,256$ (thin light grey to thick dark grey); $N_h$ is the corresponding high-pass filter from the same cut-off wavenumber $k^c_x$ to $k_x=256$ (thick dark grey to thin light grey) such that $N = N_l+N_h$ for each $k_c$. The special case $k_c = k_\beta = 16$ is shown as the thick blue line.

In figure 5 we examine the Hovmöller diagrams for the low-pass filtered stresses $M_l$ and $N_l$, and the associated tilt angle; these appear almost identical to the diagrams for the unfiltered quantities in figure 3, where they are dominated by noisy, time-dependent processes.

Figure 5. Hovmöller diagrams of the low-pass filtered, large zonal scale: (a) shear stress $N_l(y)$ and (b) its normalised time average $[N_l(y)]/10^{-3}$ (blue solid); (c) normal stress difference $M_l(y)$ and (d) its normalised time average $[M_l(y)]/10^{-2}$ (blue solid); and (e) the corresponding ellipse tilt angle $\theta _l(y)$ and (f) its time average $[\theta _l]$ normalised by ${\rm \pi} /2$ (blue solid). The normalised zonal-time average zonal velocity $[\bar {u}]/U_{max}$ is also plotted for comparison (black dashed), where $U_{max}$ is the maximum value of $[\bar {u}]$. The time averages are taken over a $20\,000\,\mathrm {h}$ equilibrium period.

The Hovmöller diagram of $N_l$ demonstrates the same alternating pattern of momentum flux observed in the shear stress. Unlike the unfiltered case, the time mean over $N_l$ (figure 5b) does not average to reveal a correlation with the jet structure. Thus we can confirm that the scales retained in the low-pass filter are not responsible for the time-mean structure in $N$ (figure 3a). The zonal pulses observed in the Hovmöller diagrams of $M$ (figure 3c) are present in $M_l$ (figure 5c). We see in the time mean of $M_l$ (figure 5d) that these scales provide the zonal structure observed in the time mean of $M$ (figure 3d). As we have noted before, these pulses are not entirely consistent with the jet pattern, though they do broadly resemble some jet structure. They do, however, seem to correlate with the alternating streaking pattern in $N_l$ (figure 5a). The Hovmöller diagram of $\theta _l$ (figure 5e) demonstrates the same alternating north and south flux direction found in $\theta$ (figure 3e), with no discernible tilting pattern in the time mean (figure 5f).

Momentum flux quantities filtered at high pass, where wavenumbers $k \geq k_c = k_\beta$ are retained (figure 6a), demonstrate the momentum flux structure that we have been searching for.

Figure 6. Hovmöller diagrams of the high-pass filtered, small and intermediate zonal scale: (a) shear stress $N_h(y)$ and (b) its normalised time average $[N_h(y)]/10^{-3}$ (blue solid); (c) normal stress difference $M_h(y)$ and (d) its normalised time average $[M_h(y)]/10^{-2}$ (blue solid); and (e) the corresponding ellipse tilt angle $\theta _h(y)$ and (f) its time average $[\theta _h]$ normalised by ${\rm \pi} /2$ (blue solid). The normalised zonal-time average zonal velocity, $[\bar {u}]/U_{max}$ is also plotted for comparison (black dashed), where $U_{max}$ is the maximum value of $[\bar {u}]$. The time averages are taken over a $20\,000\,\mathrm {h}$ equilibrium period.

From the Hovmöller diagram of $N_h$, we see that high-pass filtering reveals an alternating band structure that does not vary in time, with the exception of some noise (figure 6a). This zonal structure is consistent with the zonal jets. The strengths of these bands are much smaller than the magnitudes of the streaks in $N_l$ (figure 5a), which is why they are not apparent in the unfiltered signal $N$. We can see that these bands are a half-jet-width offset from the jet profile, providing the zonal structure observed in the time mean of $N$ with a magnitude that supports the jets (figure 6b).

The Hovmöller diagram of $M_h$ demonstrates small, negative, persistent band structures that align with the jet structure (figure 6c). The sign of $M_h$ is negative everywhere inside the domain, and is positive only at the domain walls. Also, $M_h$ has its highest amplitude, and is most negative within the jet cores. The pattern in the time mean correlates with the jet structure, though the jet structure is purely zonal. This suggests that the eddies at the scales supporting the jet have strong meridional velocity fluxes (figure 6d).

The high-pass tilt angle $\theta _h$ demonstrates a persistent banded pattern consistent with the jet (figure 6e). The pattern in the time mean is consistent with the jet structure (figure 6f) and also reveals some more subtle features. Jets must demonstrate a sharp eastward and broad westward jet structure for barotropic stability (Vallis & Maltrud 1993). A sharp eastward jet requires a sharp change in tilt angle, changing from strongly negative on its northern flank to strongly positive to the south of the jet maximum. A broad westward jet requires a less sharp change in the tilt angle, gradually changing from positive to negative from the north to the south. This results in the sawtooth shape observed in the time-mean structure of the tilt angle.

7. Discussion

We have examined a statistically stationary flow with a jet configuration of $5$ sharp eastward jets interleaved by $4$ broad westward jets. This is in agreement with the Rhines scale, which we calculate to be $k_R \sim 4$. We have also extracted the zonal and residual spectra given by (2.1a) and (2.1b), respectively, and found that they cross at wavenumber $k_\beta = 16$. These jets meet the empirically defined criteria for a zonostrophic jet as described in Galperin et al. (2006, 2008) and Sukoriansky et al. (2007), obeying the scale separation relationship in (2.4). In these studies, it was observed that strong jets formed not only when $k_R$ and $k_\beta$ are well separated, but also when $k_\beta \gtrsim 2k_R$. Independently, Scott & Dritschel (2012) found that as the ratio between $k_\beta$ and $k_R$ increases, the more the jets approach an idealised barotropically stable jet profile in which potential vorticity resembles a staircase pattern. This manifests as sharp harmonics in the total, zonal and jet spectra as discussed in Danilov & Gurarie (2004), which we have also observed here. Thus a significant scale separation between $k_R$ and $k_\beta$ is required for several harmonic peaks of the jet to develop, which in turn ensures the stability of the jet profile. Hence these viewpoints studying jets from the perspectives of stability and spectra are complementary. Eddy–eddy interactions are known to drive potential vorticity mixing (Scott & Dritschel 2019), so the small-scale eddy momentum fluxes may act to sharpen the jets and to maintain barotropic stability. However, what is unclear in these numerical experiments is why these harmonics are found to broadly trace the zonal spectrum, with the spectral peaks increasing in intensity $\propto k^{-5}$. Understanding how the spectral harmonics develop along the zonal spectrum, and the link between the spectral harmonics, the underlying small-scale eddy shearing and jet stability would be an interesting avenue of future exploration.

We have investigated eddy–mean-flow interactions of quasi-steady turbulent jets. We have confirmed the balance between the shear stress divergence and Rayleigh friction in (5.3), and we have explored this relationship in more detail using the geometric eddy ellipse formulation. On closer examination and with the aid of a zonal filter, we have found rather counterintuitively that eddies in the vicinity of the zonal spectrum, where $k_x< k_\beta$, do not support the jet structure. Rather, eddies in which $k_x>k_\beta$ are responsible for maintaining the jet. Huang & Robinson (1998) came to a similar conclusion by examining the spectral energy fluxes; here, we have been able to visualise these interactions using spatial filtering, and have been able to reveal the structure of these fluxes using the geometric eddy ellipse framework.

Whilst studies have consistently confirmed the characteristic scales of geostrophic turbulence, separated by the chain inequality (2.4), it is not clear how the anisotropisation process causes such an arrangement of eddy momentum fluxes at scales where $k \geq k_\beta$, or how this relates to these characteristic scales. Here we have found that the zonal scales that support and maintain the jet structure appear to be much smaller than the Rhines scale, and though we have arbitrarily selected $k_\beta$ as the cut-off wavenumber, it is clear that there are momentum fluxes contributing to supporting the jet that occupy higher wavenumbers than this. These scales are derived on the assumption that once turbulent frequencies are low enough to excite linear Rossby waves, the flow follows linear dynamics at scales larger than this. We can plainly see, from the highly turbulent nature of our flow fields, that this is incorrect. Strong nonlinearity exists throughout the flow dynamics, in turn presenting a slew of complications, e.g. that all linear Rossby waves can be excited at once, that their nonlinear interactions may excite additional waves called zonons (Sukoriansky, Dikovskaya & Galperin 2008; Galperin et al. 2019), and that turbulence and Rossby waves coexist over a wide range of scales. Investigating these non-local dynamics further may hold the key to understanding how these eddy–mean-flow interactions are able to flux momentum directly to the jet, seemingly bypassing the cascade. The largest and most energetic zonal scales are not responsible for maintaining the jet structure. However, $M$ still demonstrates a zonal structure with regular zonal pulses in which $\overline {u'^2}\gg \overline {v'^2}$ (figure 3). We also know through (5.2) that the tendency in the mean jet profile is governed by the divergence of the momentum flux. Other works have found that eddies close to the jet scale are responsible for jet meanders (Srinivasan & Young 2012). To better understand these large-scale processes, it would be useful to examine other eddy ellipse quantities such as $K$ and $L$, as well as studying how these evolve as the jets form.

There have been very few studies that have explored the full two-dimensional energy spectrum (e.g. Huang & Robinson 1998), focusing instead on the angularly averaged spectrum or the zonal spectrum given by (2.1a). It is clear that the scales responsible for maintaining the jet lie within the less explored residual spectrum given by (2.1b), which is often associated with isotropic dynamics. In this study, we have been limited to filtering zonal modes. Filtering meridional modes are not as trivial since the flow is not meridionally homogeneous such that the meridional mean is not equivalent to an ensemble average. Further studies will focus on developing a filter that faithfully captures the role of meridional scales. Furthermore, it would be useful to examine how anisotropy develops as a function of radial wavenumber $k$ and the angle $\phi$ of the wedge around the $k_y$-axis that defines the zonal spectrum.

8. Conclusion

This study has investigated a jet structure that formed in $\beta$-plane turbulence using a barotropic channel model. We observed that an isotropically forced flow developed several alternating jets as its energy approached scales in which $k< k_\beta$ and the flow attained statistical stationarity through the influence of Rayleigh friction. The final jet number $k_{jet} \sim 5$ was close to $k_{R} \sim 4$, the Rhines scale. We calculated the zonal and residual spectra, and found that they crossed at $k_\beta = 16$, with the zonal spectrum steepening and broadly tracing the $k^{-5}$ slope, possessing several harmonics. The jet spectrum $k_{jet}$ also consisted of these harmonics, which corresponded to the asymmetric pattern of sharp eastward and broad westward jets necessary for a barotropically stable jet profile. These characteristic scales were well separated and obeyed the chain inequality in (2.4), thus this system was found to be in the zonostrophic regime.

We then investigated the eddy momentum fluxes over a period in which the flow had reached statistical stationarity, and observed that the eddy momentum flux into the jet structure balanced the Rayleigh friction. However, the eddy–mean-flow relationship described by (5.3) held true only after a long-time averaging with significant time-dependent processes producing a tendency in the mean at shorter time scales. We employed the geometric eddy ellipse formulation to study the eddy–mean-flow relationship, and analysed Hovmöller diagrams of the shear stress $N$, the normal stress difference $M$, and the ellipse tilt angles $\theta$. Our findings revealed that large amplitude, short time scale momentum fluxes dominated the momentum flux field. These existed on the most energetic zonal scales in the vicinity of the Rhines scale, where $k< k_\beta$. In contrast, we found that the small to intermediate scales of eddies, where $k>k_\beta$, exhibited a clear, persistent momentum flux pattern. These momentum fluxes had a low amplitude compared to the short time scale momentum fluxes that dominated the larger scales, hence they were masked. Despite their low amplitude, it was these small to intermediate scale momentum fluxes, which occupied scales typically associated with isotropic dynamics, that played the crucial role of supporting the jet structure. These findings demonstrate the complexity of eddy–mean-flow interactions in geostrophic turbulent jets, and confirm the significance of non-local interactions, in which small scales eddies directly funnel momentum to the largest scales.