# Gyre turbulence: Anomalous dissipation in a two-dimensional ocean model

Lennard Miller, Bruno Deremble, Antoine Venaille

► To cite this version:

Lennard Miller, Bruno Deremble, Antoine Venaille. Gyre turbulence: Anomalous dissipation in a two-dimensional ocean model. Physical Review Fluids, 2024, 9 (5), pp.L051801. 10.1103/Phys-RevFluids.9.L051801. hal-04765467

HAL Id: hal-04765467

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

Submitted on 16 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Gyre Turbulence: Anomalous Dissipation in a Two-Dimensional Ocean Model

Lennard Miller¹,²,*, Bruno Deremble²,† and Antoine Venaille¹‡

¹ENS de Lyon, CNRS, Laboratoire de Physique (UMR CNRS 5672),
F-69342 Lyon, France, ² Université Grenoble Alpes, CNRS, INRAE, IRD,
Grenoble-INP, Institut des Géosciences de l'Environnement, Grenoble, France

The exploration of a two-dimensional wind-driven ocean model with no-slip boundaries reveals the existence of a turbulent asymptotic regime where energy dissipation becomes independent of fluid viscosity. This asymptotic flow represents an out-of-equilibrium state, characterized by a vigorous two-dimensional vortex gas superimposed onto a western-intensified gyre. The properties of the vortex gas are elucidated through scaling analysis for detached Prandtl boundary layers, providing a rationalization for the observed anomalous dissipation. The asymptotic regime demonstrates that boundary instabilities alone can be strong enough to evacuate wind-injected energy from the large-scale oceanic circulation.

## Introduction.
In three-dimensional turbulence, the transfer of energy from large to small scales results in an energy dissipation rate that remains independent of viscosity, regardless of its smallness [1, 2]. This dissipative anomaly is a robust empirical observation sometimes referred to as the zeroth law of turbulence [3]. Conversely, two-dimensional flows are subject to an inverse energy cascade [4], which results in the self-organization of the flow at the domain scale [5]. To ensure efficient dissipation in two-dimensional flows, an additional mechanism is required to generate small-scale structures where dissipation can operate, thus disrupting the inverse cascade. The strong shear near lateral boundaries could serve as a means to create such dissipative structures [6]. Numerical studies of the bounded Navier-Stokes equations have examined the relationship between dissipation and viscosity during a dipole-wall collision [7–9]. The possibility for dissipation to remain finite in the inviscid limit has been raised [10, 11], yet remains a topic of debate [12, 13].

The existence of a dissipative anomaly in two-dimensional flows with boundaries would bear significant practical implications, for example in contributing to a deeper understanding of the energy cycle in the ocean [14–17]. In fact, classical linear models for the emergence of western intensified currents [18, 19], such as the Gulf Stream or the Kuroshio, do provide a remarkable example of a dissipative anomaly. Here we show that this dissipative anomaly persists in a nonlinear regime, and unveil a new Gyre Turbulence regime with a western intensified mean flow and finite energy dissipation rate.

## Flow model.
The simplest model describing western intensification of oceanic currents is the rigid-lid barotropic quasigeostrophic model on a closed domain tangent to the Earth [18]:

$$ \partial_t \omega + J(\psi, \omega) + \beta^* \partial_x \psi = \nu^* \Delta \omega - \partial_y \tau, \quad (1) $$

$$ \omega = \Delta \psi, \quad \tau = -\cos(\pi y). \quad (2) $$

The term $J(\psi, \omega) = (u\partial_x + v\partial_y)\omega$ is the advection of vorticity $\omega$ by the streamfunction $\psi$, with $u = -\partial_y\psi$ the zonal (x-direction) and $v = \partial_x\psi$ the meridional (y-direction) velocity components. We solve this two-dimensional model on a square domain with no-slip boundary conditions. Forcing comes from the wind stress curl $-\partial_y\tau$. Time and length have been rescaled such that both the length of the domain $L$ and the maximum value of the wind stress $\tau_0$ are 1. The shape of the forcing corresponds to a single gyre, and is somewhat relevant to the North Atlantic case, with net injection of negative vorticity, as observed in the subtropical gyre (figure 1A). The only difference to the incompressible two-dimensional Navier-Stokes equations is the term $\beta^*\partial_x\psi$. This term comes from the curl of the Coriolis force, assuming linear variations of the Coriolis parameter in the meridional direction $y$. This framework is called the beta-plane approximation, and it captures the effects of differential rotation induced by a rotating planet [18].

Earlier studies of the single gyre model with free-slip boundary conditions exhibited an inertial run-away in the invicid limit: a solution with unrealistically large velocities, and lacking an intensified western boundary current [20–26]. One way to prevent this inertial run-away is to use no-slip conditions [27–29]. It was also noticed that no-slip boundary conditions in two-dimensional turbulence can drastically alter flow organization [30–32]. In the context of a wind-driven shallow water model, [29] observed that no-slip conditions induced intermittent formation of intense vortices through boundary layer detachment. We show below that this mechanism enables the large-scale flow to remain in an oceanic gyre regime when considering the invicid limit of (1–2) (see Supplementary Material [33], movie S1 for a direct comparison of the free-slip against the no-slip regimes).

<footnote>* lennard.miller@univ-grenoble-alpes.fr</footnote>

<footnote>† bruno.deremble@univ-grenoble-alpes.fr</footnote>

<footnote>‡ antoine.venaille@ens-lyon.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: **A:** The β-plane model with an inlet for the streamlines of the time-average sea surface velocity field in the North Atlantic. Colors show the corresponding relative vorticity field. **B:** Temporal average of the streamlines in the Gyre Turbulence regime ($\nu^* = 2.5 \times 10^{-5}$), superimposed on the corresponding relative vorticity field (colors) **C:** Snapshot of relative vorticity in the Gyre Turbulence regime ($\nu^* = 2 \times 10^{-6}$). Blue: anticyclonic vortices; red: cyclonic vortices (see Supplementary Material [33], movie S2 for an animation).</figcaption></figure>

*Linear dynamics.* Linear theories for wind-driven gyres compute steady states of equations (1)-(2), by neglecting the advection term [18, 19]. In the domain bulk, the vorticity equation simplifies into Sverdrup balance, a cornerstone of midlatitude ocean dynamics: $\beta^*v = -\partial_y\tau$, meaning that an injection of negative vorticity is balanced by a southward transport of the fluid. To ensure mass conservation, this interior circulation must be complemented by boundary layers. The majority of this recirculation occurs within the western boundary layer, thereby breaking the East-West symmetry established by the Sverdrup balance. In the viscous solution found by Munk (figure 2A, top left inlet) the boundary layer thickness scales as $\delta_M = (\nu^*/\beta^*)^{1/3}$ [34], which implies that total dissipation is dominated by contributions from the boundary layer while energy injection comes from the domain bulk. The confinement of energy dissipation in a western boundary layer holds when viscosity is replaced by other dissipation mechanisms, such as a linear drag [18]. In this scenario, large-scale gyre patterns and therefore energy injection do not depend on details of those linear boundary layers. This is a strong incentive to look for a turbulent dissipative anomaly in this system, adding back nonlinearities into the problem.

*Parameter regimes.* The linear Munk solution is a limiting case in the $(\nu^*, \beta^*)$ parameter space, corresponding to $\beta^* \gg 1$ and a restricted range of $\nu^*$ to be discussed below. The parameter space can be dissected into four regions (figure 2A): in the limit of weak $\beta^*$, the effects of differential rotation are negligible with respect to other terms and the flow response is equivalent to that of the two-dimensional Navier-

Stokes (NS) equations, with a transition from laminar to turbulent flow when $\nu^*$ decreases [31]. The same transition from laminar to turbulent flow occurs when differential rotation is important ($\beta^* \gg 1$). When $\nu^* = 10^{-2}$, we observe that the flow consists of a domain interior in which a Sverdrup balance holds with and a stationary western boundary layer, similar to the linear Munk solution.

When viscosity is decreased, the boundary layer becomes increasingly inertial and its thickness will be governed by $\delta_I = \beta^{*-1}$ [18, 35]. The transition between the laminar Munk regime and this inertial regime occurs around $\delta_M = \delta_I$, i.e. $\nu^* = 1/\beta^{*2}$. When viscosity is decreased further below this threshold the boundary layer becomes unstable, and those instabilities feed the domain with filaments and vortices. As friction is further reduced the system eventually enters into the Gyre Turbulence regime: while the time mean flow remains close to a Sverdrup interior with western intensified boundary layers (figure 1B), the instantaneous vorticity field is dominated by a vigorous heterogeneous vortex gas which is densest in the north-western corner of the domain (figure 1C, for an animation see Supplementary Material [33], movie S2). The gyre structure is not observed on instantaneous streamfunction fields, which is instead dominated by contributions arising from Rossby basin modes [36] (not shown here), and to a lesser extent by contributions from vortices.

*Energy budget.* The central result of this article is depicted in figure 2B showing that the dissipative anomaly of the linear Munk regime persists in the nonlinear Gyre turbulence regime. Dissipation $\epsilon$ ap-
<figure><img src="image_3.png" /><figcaption>Figure 2: A: Parameter space (no-slip boundary conditions). Insets show typical instantaneous vorticity fields (color) and streamfunction (lines) B: Energy dissipation rate, rescaled by the energy injection through a Sverdrup interior. C: Energy rescaled by the energy of the inertial western boundary layer. Errorbars show the standard deviation.</figcaption></figure>

pears in the energy budget as

$$ \frac{\partial E}{\partial t} = \mathcal{P} - \varepsilon, \quad \mathcal{P} = \int \tau u \, dA, \quad \varepsilon = 2\nu^* Z, \quad (3) $$

where $E = \int (u^2 + v^2) \, dA/2$ is the total energy and $Z = \int \omega^2 \, dA/2$ the total entrophy. We define time average and fluctuations as

$$ \langle \psi \rangle = \lim_{T \to \infty} \frac{1}{T} \int_s^{s+T} \psi \, dt, \quad \psi' = \psi - \langle \psi \rangle, \quad (4) $$

where the integration is started from a moment $s$ in time after which the system is observed to be in statistical equilibrium, for which $\langle \mathcal{P} \rangle = \langle \varepsilon \rangle$.

In Gyre Turbulence, the average dissipation $\langle \varepsilon \rangle$ is observed to be insensitive to a decrease in $\nu^*$ (Fig. 2B). Our rescaling by $\beta^{*-1}$ shows that the dissipation rate remains close to that predicted by energy injection through an interior flow governed by a Sverdrup balance. However, the dissipation mechanism changes drastically as fluctuations become increasingly important at smaller values of $\nu^*$, while the mean flow contributes only a negligible fraction to the total dissipation. The importance of the fluctuations can also be seen in the total energy of the flow (Fig. 2C). The rescaling by $\beta^{*-1}$ reveals that the total energy of the mean flow remains close to the energy in a western inertial boundary layer, while the magnitude of the fluctuations increases when $\nu^*$ decreases. To rationalize these observations, we describe below the mean flow structure (related to to energy injection), and the production of

filaments and vortices (related to energy dissipation).

**Mean flow structure.** In the Gyre Turbulence regime, the time-averaged production term $\langle \mathcal{P} \rangle$ must be independent of $\nu^*$, as it balances the time-averaged dissipation $\langle \varepsilon \rangle$. This constraint, together with the observation in figure 2B that the energy of the mean flow reaches a plateau, suggests that the bulk streamfunction $\langle \psi \rangle$ displayed in figure 1B only has a weak dependence on $\nu^*$ in the Gyre Turbulence regime. While the order of magnitude of the mean flow agrees with Sverdrup balance in the domain interior (figure 3A), the gyre pattern is different than the prediction from linear theory, likely due to nonlinear rectification mechanisms in the presence of Rossby waves [32, 37]. We observe only weak changes in this pattern when lowering viscosity in the Gyre Turbulence regime (see Supplementary material [33], figure 2).

The gyre pattern is connected to inertial boundary layers with a well defined functional relation between streamfunction and potential vorticity $\langle q \rangle = \langle \omega \rangle + \beta y$, such that $J(\langle q \rangle, \langle \psi \rangle) = 0$. We identified two regions where such relations hold (figure 3), which we will call the inflow (blue) and the outflow (red) layer. The inflow layer is consistent with classical theory predicting $\langle q \rangle = -(\beta^*/U_{in})\langle \psi \rangle$, which leads to a western boundary layer thickness $\delta_I = 1/\beta^*$, where $U_{in} < 0$ is the westward inflow scaling as $1/\beta^*$ [18, 35]. The northward velocity in the inertial boundary layer is $U_I \sim 1$, so that mass transport in this layer compensates the southward Sverdrup bulk transport. No-slip
<figure><img src="image_4.png" /><figcaption>Figure 3: A: Contour plot of the mean stream function $\langle\psi\rangle$ at $\nu^* = 2.5 \times 10^{-5}$, with colored areas identified as inertial boundary layers. B: $q-\psi$ relations in the inertial boundary layers.</figcaption></figure>

boundary condition is guaranteed by a Prandtl sublayer with thickness $\delta_P \sim \sqrt{\nu^*}$. The vorticity within this Prandtl layer is

$$ \omega_{max} = \frac{U_I}{\delta_P} \sim \frac{1}{\sqrt{\nu^*}}. \qquad (5) $$

Assuming that dissipation of the mean flow is governed by these viscous sublayers yields

$$ \nu^* \int \langle\omega\rangle^2 dA \sim \nu^* \omega_{max}^2 \delta_P \sim \nu^{*\frac{1}{2}} \qquad (6) $$

which compares well against the observed dissipation due to the mean flow (see figure 2B).

The outflow layer close to the northern boundary corresponds to a meandering jet with velocity $U_{out}$ and a strong cyclonic recirculation. In this area, we observe a negative correlation between $\langle q \rangle$ and $\langle \psi \rangle$. Stationary Rossby wave meanders and a stationary vortex on a beta-plane with mean flow $U_{out}$ both select the size $\delta_{out} \sim \sqrt{U_{out}/\beta^*}$. We assume that this length sets the vortex and jet width, and that jet transport is set by the transport of the western boundary layer, which yields

$$ \delta_{out} \sim \beta^{*-2/3} \quad \text{and} \quad U_{out} \sim \beta^{*-1/3}. \qquad (7) $$

An adaptation of classical Charney theory to this inertial region leads to $\langle q \rangle = -\beta^*/U_{out}\langle\psi\rangle + \beta^*$, which fits well with numerical results (see figure 3B and Supplementary Material [33] for further information on boundary layers).

**Statistics of dissipation.** While the production term $\langle P \rangle$ depends crucially on the mean flow, we showed that, in the turbulent regime, the dissipation term is dominated by contributions from fluctuations of vorticity (figure 2B). Since dissipation is proportional to entrophy, we show the probability distribution function of $\omega^2$ in figure 4A (we checked that

the distribution of $\omega^2$ is close to the distribution of local dissipation $|\nabla u|^2$). The core of the distribution $\Pi(\log(\omega^2))$ is close to a Gaussian distribution, similar to recent observations from more comprehensive ocean models [17].

The peak of $\Pi(\log(\omega^2))$ changes little upon varying the viscosity, revealing that most of the values of vorticity in the bulk are only weakly dependent on $\nu^*$. A much stronger dependence on $\nu^*$ is observed in the tails of the distributions, where we notice important deviations from lognormality. In fact, the average dissipation is dominated by contributions from the tails, which can be seen by plotting the quantity $\omega^2\Pi$. We show these tails in figure 4B after rescaling the vorticity by $\omega_{max}$ defined in (5), collapsing the tails onto a single curve centered around unity. This confirms that the primary mechanism for vorticity injection is the detachment of the Prandtl layer. This scaling for the vorticity can now be used to estimate the energy dissipation rate.

**Scaling Analysis of the Vortex Gas.** The cyclonic recirculation we observe in the mean flow is a signature of the presence of intense cyclonic vortices in this region (see figure 1C). These large vortices result from the coalescence of smaller vortices in the interior. If the vortices grow too large, their drift velocity overcomes the advective velocity from the detached jet, and the vortex collides with the wall, detaching the sublayer over a length $\delta_{out}$ (see Supplementary Material [33], movie S2).

Building upon this observation, the theory for the mean flow and the observed statistics of the dissipation, it is possible to propose a scaling theory that predicts the $\nu^*$-dependency of the energy of the vortex gas. For this, we suppose that the energy is determined by $N$ characteristic vortices of size $\lambda$. If we assume that the vortices are created by the roll-up of the detached sublayer of the western boundary along a length $\delta_{out}$, then the characteristic vortex radius scales as $\lambda \sim \sqrt{\delta_{out}\delta_P}$. Furthermore, in the inviscid limit we observe a balance between the dissipation of these vortices and the injection of energy through the mean Sverdrup flow. The balance between energy dissipation through the characteristic vortices and injection through a Sverdrup interior takes the form

$$ \nu^* N \lambda^2 \omega_{max}^2 \sim \beta^{*-1}. \qquad (8) $$

This determines the area fraction of the vortices as $N\lambda^2 \sim \beta^{*-1}$, which we verify (see Supplementary Material [33], figure 2). The total energy of the vortex gas is then given by

$$ E \sim N \lambda^2 (\omega_{max} \lambda)^2 \sim \nu^{*-1/2} \beta^{*-5/3}, \qquad (9) $$

with additional logarithmic corrections [38]. Although the ideas behind this scaling do not incorporate the fragmentation or coalescence and the resulting inhomogeneity of the vortex gas, a
<figure><img src="image_5.png" /><figcaption>Figure 4: **A:** Probability Distribution functions Π of $log(ω^2)$ for different values of $v^*$. Also shown are the values responsible for dissipation, $ω^2Π$. **B:** Values responsible for dissipation after rescaling by $ω_{max} = 1/√v^*$. The area under the curves is equal to the total dissipation $⟨ε⟩$.</figcaption></figure>

reasonable agreement between simulations and this scaling theory is observed (figure 2C).

*Discussion.* Our study challenges the consensus that a weakly dissipated two-dimensional ocean model would lead to energy accumulation due to the absence of a forward energy cascade [16, 18]. Instead, no-slip boundaries offer a compelling route to dissipation, sustaining low-energy western-intensified gyres with a vigorous vortex gas. This Gyre Turbulence regime is a genuine out-of-equilibrium state defying description using equilibrium statistical mechanics [2, 5], or quasilinear approaches [39, 40].

Originally discovered in freely decaying 2D flows [41], the vortex gas regime emerges in baroclinically unstable stratified quasi-geostrophic flows or any two-dimensional flow with forcing amplifying vorticity extrema [38, 42–45]. Gyre Turbulence, lacking bulk instabilities, generates intense vortex cores solely at boundary layers. The interplay of a western boundary and the $β$ effect also seems to prevent crystallisation previously observed on a polar cap [46], although this feature may yet occur at lower viscosity not attained in the presented simulations. A quantitative description of the large-scale gyre pattern will require an examination of potential vorticity mixing induced by the vortex gas, as in [38, 47, 48].

In the broader context of the real ocean, turbulent

western boundaries have long been recognized as significant energy sinks [15], with the potential for a linear forward cascade mediated by planetary vorticity gradients [49] and intense vortex generation through boundary layer detachment [29]. Our study reveals that incompressible two-dimensional $β$-plane turbulence serves as a minimal model sustaining wind gyres co-existing with mesoscale vortices that contain the majority of the flow energy. However, the two-dimensional model falls short in capturing the fate of these vortices, as evidenced by their energy increase with $v^{*-1/2}$. Three-dimensional and ageostrophic effects [16, 50], interactions with bottom topography [51–53], and air-sea interactions [54] are expected to play vital roles in dissipating these fine-scale structures. We conjecture that a large part of the energy flux toward these fine scales will continue to be governed by larger scales, a hypothesis currently under examination using more comprehensive ocean models.

*Acknowledgements.* This project has received financial support from the CNRS through the 80 Prime program, and was performed using HPC resources from GENCI-TGCC (Grant 2022-A0130112020). We warmly thank G. Roullet for insightful inputs on this topic.

[1] U. Frisch, *Turbulence: the legacy of AN Kolmogorov* (Cambridge university press, 1995).

[2] G. L. Eyink and K. R. Sreenivasan, Onsager and the theory of hydrodynamic turbulence, Reviews of modern physics **78**, 87 (2006).

[3] B. Dubrulle, Beyond kolmogorov cascades, Journal of Fluid Mechanics **867**, P1 (2019).

[4] G. Boffetta and R. E. Ecke, Two-dimensional turbulence, Annual review of fluid mechanics **44**, 427 (2012).

[5] F. Bouchet and A. Venaille, Statistical mechanics of two-dimensional and geophysical flows, Physics reports **515**, 227 (2012).

[6] B. Deremble, W. Dewar, and E. Chassignet, Vorticity dynamics near sharp topographic features, J. Mar. Res **74**, 249 (2016).

[7] H. J. Clercx and G. J. F. van Heijst, Dissipation of kinetic energy in two-dimensional bounded flows, Phys-
ical Review E 65, 066305 (2002).

[8] W. Kramer, H. Clercx, and G. Van Heijst, Vorticity dynamics of a dipole colliding with a no-slip wall, Physics of Fluids 19 (2007).

[9] G. Keetels, W. Kramer, H. Clercx, and G. van Heijst, On the reynolds number scaling of vorticity production at no-slip walls during vortex-wall collisions, Theoretical and computational fluid dynamics 25, 293 (2011).

[10] R. Nguyen van yen, M. Farge, and K. Schneider, Energy dissipating structures produced by walls in two-dimensional flows at vanishing viscosity, Physical Review Letters 106, 184502 (2011).

[11] M. Waidmann, R. Klein, M. Farge, K. Schneider, et al., Energy dissipation caused by boundary layer instability at vanishing viscosity, Journal of Fluid Mechanics 849, 676 (2018).

[12] D. Sutherland, C. Macaskill, and D. Dritschel, The effect of slip length on vortex rebound from a rigid boundary, Physics of Fluids 25, 093104 (2013).

[13] H. Clercx and G. Van Heijst, Dissipation of coherent structures in confined two-dimensional turbulence, Physics of Fluids 29, 111103 (2017).

[14] R. Ferrari and C. Wunsch, Ocean circulation kinetic energy: Reservoirs, sources, and sinks, Annual Review of Fluid Mechanics 41, 253 (2009).

[15] X. Zhai, H. L. Johnson, and D. P. Marshall, Significant sink of ocean-eddy energy near western boundaries, Nature Geoscience 3, 608 (2010).

[16] A. Pouquet and R. Marino, Geophysical turbulence and the duality of the energy flow across scales, Physical review letters 111, 234501 (2013).

[17] B. Pearson and B. Fox-Kemper, Log-normal turbulence dissipation in global ocean models, Physical review letters 120, 094501 (2018).

[18] G. K. Vallis, Atmospheric and oceanic fluid dynamics (Cambridge University Press, 2017).

[19] A.-L. Dalibard and L. Saint-Raymond, Mathematical study of degenerate boundary layers: A large scale ocean circulation problem, American Mathematical Society 253 (2018).

[20] V. Kamenkovich, V. Sheremet, A. Pastushkov, and S. Belotserkovsky, Analysis of the barotropic model of the subtropical gyre in the ocean for finite reynolds numbers. part i, Journal of marine research 53, 959 (1995).

[21] V. Sheremet, V. Kamenkovich, and A. Pastushkov, Analysis of the barotropic model of the subtropical gyre in the ocean for finite reynolds numbers. part ii., Journal of marine research 53, 995 (1995).

[22] G. Ierley and V. Sheremet, Multiple solutions and advection-dominated flows in the wind-driven circulation. part i: Slip, Journal of marine research 53, 703 (1995).

[23] V. Sheremet, G. Ierley, and V. Kamenkovich, Eigenanalysis of the two-dimensional wind-driven ocean circulation problem, Journal of marine research 55, 57 (1997).

[24] F. Dupont and D. N. Straub, Effect of a wavy wall on the single gyre Munk problem, Tellus A: Dynamic Meteorology and Oceanography 56, 387 (2004).

[25] P. S. Berloff, J. C. McWilliams, and A. Bracco, Material transport in oceanic gyres. part i: Phenomenology, Journal of Physical Oceanography 32, 764 (2002).

[26] B. Fox-Kemper and J. Pedlosky, Wind-driven barotropic gyre i: Circulation control by eddy vorticity fluxes to an enhanced removal region, Journal of Marine Research 62, 169 (2004).

[27] P. Cessi, R. V. Condie, and W. Young, Dissipative dynamics of western boundary currents, Journal of Marine Research 48, 677 (1990).

[28] H. Nakano, H. Tsujino, and R. Furue, The kuroshio current system as a jet and twin "relative" recirculation gyres embedded in the sverdrup circulation, Dynamics of atmospheres and oceans 45, 135 (2008).

[29] C. Q. C. Akuetevi and A. Wirth, Dynamics of turbulent western-boundary currents at low latitude in a shallow-water model, Ocean Science 11, 471 (2015).

[30] H. Clercx, S. Maassen, and G. Van Heijst, Decaying two-dimensional turbulence in square containers with no-slip or stress-free boundaries, Physics of Fluids 11, 611 (1999).

[31] H. Clercx, G. Van Heijst, D. Molenaar, and M. Wells, No-slip walls as vorticity sources in two-dimensional bounded turbulence, Dynamics of atmospheres and oceans 40, 3 (2005).

[32] W. Kramer, M. van Buren, H. Clercx, and G. van Heijst, β-plane turbulence in a basin with no-slip boundaries, Physics of Fluids 18, 026603 (2006).

[33] See Supplemental Material at [URL will be inserted by publisher] for details on numerical methods, boundary layer analysis and movies.

[34] W. H. Munk, On the wind-driven ocean circulation, Journal of Atmospheric Sciences 7, 80 (1950).

[35] J. G. Charney, The gulf stream as an inertial boundary layer, Proceedings of the National Academy of Sciences 41, 731 (1955).

[36] J. Pedlosky, Geophysical fluid dynamics (Springer Science & Business Media, 2013).

[37] B. Fox-Kemper, Wind-driven barotropic gyre ii: Effects of eddies and low interior viscosity, Journal of Marine Research 62, 195 (2004).

[38] B. Gallet and R. Ferrari, The vortex gas scaling regime of baroclinic turbulence, Proceedings of the National Academy of Sciences 117, 4491 (2020).

[39] B. Gallet and W. Young, A two-dimensional vortex condensate at high reynolds number, Journal of Fluid Mechanics 715, 359 (2013).

[40] A. Frishman and C. Herbert, Turbulence statistics in a two-dimensional vortex condensate, Physical review letters 120, 204505 (2018).

[41] J. C. McWilliams, The emergence of isolated coherent vortices in turbulent flow, Journal of Fluid Mechanics 146, 21 (1984).

[42] B. K. Arbic and G. R. Flierl, Baroclinically unstable geostrophic turbulence in the limits of strong and weak bottom ekman friction: Application to midocean eddies, Journal of Physical Oceanography 34, 2257 (2004).

[43] A. F. Thompson and W. R. Young, Scaling baroclinic eddy fluxes: Vortices and energy balance, Journal of physical oceanography 36, 720 (2006).

[44] A. Venaille, G. K. Vallis, and K. S. Smith, Baroclinic turbulence in the ocean: Analysis with primitive equation and quasigeostrophic simulations, Journal of Physical Oceanography 41, 1605 (2011).

[45] A. van Kan, B. Favier, K. Julien, and E. Knobloch,
Spontaneous suppression of inverse energy cascade
in instability-driven 2-d turbulence, Journal of Fluid
Mechanics **952**, R4 (2022).

[46] L. Siegelman, W. R. Young, and A. P. Ingersoll, Po-lar vortex crystals: Emergence and structure, Pro-ceedings of the National Academy of Sciences **119**, e2120486119 (2022).

[47] B. Gallet and R. Ferrari, A quantitative scal-ing theory for meridional heat transport in plan-etary atmospheres and oceans, AGU Advances **2**, e2020AV000362 (2021).

[48] L. Siegelman and W. Young, Two-dimensional turbu-lence above topography: Vortices and potential vor-ticity homogenization, Proceedings of the National Academy of Sciences **120**, e2308018120 (2023).

[49] B. Nadiga and D. Straub, Alternating zonal jets and energy fluxes in barotropic wind-driven gyres, Ocean Modelling **33**, 257 (2010).

[50] D. Balwada, J.-H. Xie, R. Marino, and F. Feraco, Di-rect observational evidence of an oceanic dual kinetic energy cascade and its seasonality, Science Advances **8**, eabq2566 (2022).

[51] W. K. Dewar and A. M. Hogg, Topographic inviscid dissipation of balanced flow, Ocean Modelling **32**, 1 (2010).

[52] M. Nikurashin, G. K. Vallis, and A. Adcroft, Routes to energy dissipation for geostrophic flows in the southern ocean, Nature Geoscience **6**, 48 (2013).

[53] L.-F. Zhang and J.-H. Xie, Spectral condensation in quasi-geostrophic turbulence above small-scale to-pography, arXiv preprint arXiv:2311.16612 (2023).

[54] L. Renault, P. Marchesiello, S. Masson, and J. C. McWilliams, Remarkable control of western bound-ary currents by eddy killing, a mechanical air-sea cou-pling process, Geophys. Res. Lett. **46**, 2743 (2019).

[55] B. Cushman-Roisin and J.-M. Beckers, *Introduction to geophysical fluid dynamics: physical and numerical aspects* (Academic press, 2011).
## SUPPLEMENTARY MATERIALS 1: INERTIAL BOUNDARY LAYERS

In the theory of inertial inflow boundary layers the relation between $\psi$ and $q$ is established in an intermediate connection area just outside the boundary layer, where we expect the same functional relation between $\psi$ and $q$ to hold as inside the boundary layer (for details see [18]). In the limit of vanishing boundary layer thickness, the dynamical balance in the boundary layer will be the advection of relative vorticity by the interior flow and the creation of relative vorticity by the $\beta$-term. The thickness then scales like $\delta_I \sim \sqrt{-U_{in}/\beta^*}$, and the outer matching yields $q = -\beta^*\psi/U_{in}$.

In the simulations we identify two areas in which a functional relation between $q$ and $\psi$ holds. To identify these regions, we draw the isolines $u, v = 0$ that cross the centre of the gyre (the maximum of $\langle\psi\rangle$ close to the western boundary). The inflow layer is defined as the area west of the line $v = 0$ and south of $u = 0$, whereas the outflow layer is defined as the area north of $u = 0$ and east of $v = 0$. We also cut out areas very close to the boundaries where we expect viscous effects to be dominant. The resulting areas are depicted in figure 3 in the main text.

For the outflow layer (with positive inflow velocity along the x-direction), a matching with the interior flow is not self-consistent as the boundary layer thickness becomes imaginary and the solution becomes oscillatory. Furthermore, the inertial region does not reach the western part of the domain.

We noticed that this inertial region contains a jet meandering around a closed cyclonic recirculation zone. Calling $\delta_o$ the jet width, calling $U_{out}$ the jet velocity, and assuming continuity of the mass transport between the western boundary layer and the meandering jet results in

$$ \delta_o U_{out} = U_I \delta_I = U_{Sv} L. \qquad (10) $$

We observe that the wavelength of jet meanders and the radius of the cyclonic recirculation zone are of the order of the jet width $\delta_o$. Meanders can be interpreted as a doppler-shifted Rossby wave, which is stationary when

$$ \delta_o = \sqrt{U_{out}/\beta^*}. \qquad (11) $$

This length also corresponds to the typical size of a stationary vortex on a beta plane in the presence of a mean flow $U_{out}$. Combining Eq. (11) and (10) yields

$$ \delta_o = \beta^{*-\frac{2}{3}}, U_{out} = \beta^{*-\frac{1}{3}}. \qquad (12) $$

To determine the functional relation between $\psi$ and $q$ in the outflow layer we consider the (non-inertial) region in the northwestern part of the domain. This region connects the western boundary layer in the inflow layer to the meandering jet of the outflow layer. Close to the outflow layer, we assume that the flow velocity is zonal with velocity $U_{out} > 0$ and we neglect relative vorticity with respect to planetary vorticity (which is only marginally satisfied). In that case we can apply the same reasoning as classical inertial layer theory with $\langle q \rangle = \beta^*y$ and $\langle \psi \rangle = -U_{out}(y - L)$, assuming $\psi = 0$ at the northern wall in the outflow layer. We hence retrieve the relation.

$$ \langle q \rangle = -\beta^*/U_{out}\langle \psi \rangle + \beta^*L. \qquad (13) $$

To test this relation, we diagnose the velocities $U_{in}$ and $U_{out}$ as the average zonal velocities along the isolines $v = 0$ for each part of the gyres connected to the matching regions. This yields values of $U_{out} = 0.18$ and $U_{in} = 0.07$, for which the scaling theory gives the correct order of magnitude ($U_{out} \sim \beta^{*-\frac{1}{3}} = 0.21$, $U_{in} \sim \beta^{*-1} = 0.01$, prefactors of order $\pi$ not taken into account). Using the diagnosed values, the slopes given by inertial layer theory then agree well with the observed functional relations between $\langle \psi \rangle$ and $\langle q \rangle$ in both areas. The slight offset of the relation for the outflow may stem from the fact that the relative vorticity is not negligible in the matching region.

## SUPPLEMENTARY MATERIALS 2: MEAN FLOW VARIATIONS

Although the energy of the mean flow remains constant its structure changes slightly when decreasing $\nu^*$ (figure 5). In the interior, the linear response of the system to the forcing gives the Sverdrup solution of the form
$$ \psi_{Sv} = \frac{1}{\beta^*} \pi(1-x) \sin(\pi y) \qquad (14) $$

If we subtract this from the mean stream function $\langle\psi\rangle$ we observe that the departures from this solution are mainly zonal structures and hence suggest that they results from the non-linear interaction of Rossby basin modes as in [37] (also no inertial relation between $\langle\psi\rangle$ and $\langle q\rangle$ was observed for these zonal structures).

<figure><img src="image_6.png" /><figcaption>Figure 5: Top line: Mean stream function $\langle\psi\rangle$ for different values of $\nu^*$. Bottom line: Deviations from Sverdrup Flow. Contour lines are plotted with intervals $\beta^{*-1}$.</figcaption></figure>

## SUPPLEMENTARY MATERIALS 3: NUMERICAL METHODS

For the numerical implementation of the quasigeostrophic equations we apply standard finite-difference discretization procedure using Basilisk software (http://basilisk.fr). Vorticity and stream function are collocated at cell vertices and we use the Arakawa Jacobian for the advection term [55]. For the inversion of the elliptic equation we use a multigrid method. Time integration is performed with a second order predictor-corrector scheme. Impermeability conditions are achieved by imposing $\psi = 0$ at the edges when inversing $\omega$. The no-slip boundary condition is implemented by specifying the value of the vorticity on the domain boundaries for the viscous operator. By performing a Taylor expansion in the vicinity of the boundary (in this example, the western boundary) we get

$$ \psi_1 = \psi_0 + \Delta \left. \frac{\partial \psi}{\partial x} \right|_0 + \left. \frac{\Delta^2}{2} \frac{\partial^2 \psi}{\partial x^2} \right|_0. \qquad (15) $$

So at this order, the vorticity at the edge is

$$ \omega_0 = \left. \frac{\partial^2 \psi}{\partial x^2} \right|_0 = 2 \frac{\psi_1}{\Delta^2}, \qquad (16) $$

since $\psi_0 = 0$ (no flow) and $\partial\psi/\partial x|_0 = 0$ (no slip).

In practice the methods outlined above work well in the domain interior, we therefore use them to timestep $\omega$ in the interior. We then inverse $\omega$ to obtain $\psi$, from which we can at last calculate the new values of $q$ at the boundaries with equation 16.

Almost all simulations that are discussed in the main text were performed at a numerical resolution of 2048 × 2048 gridpoints. The only exception is the most turbulent run (at $\nu^* = 2 \times 10^{-6}$), which was performed at a resolution of 8192 × 8192 gridpoints. To check for numerical convergence the simulation at $\nu^* = 6.25 \times 10^{-6}$ was relaunched with a doubled numerical resolution. Both its mean energy $\langle E \rangle$ and its mean dissipation $\langle \epsilon \rangle$ changed by less than 5%.
## SUPPLEMENTARY MATERIALS 4: DISSIPATIVE AREA FRACTION

The scaling analysis of the vortex gas presented in the main text is further validated by observing the area fraction of dissipative structures. We define the dissipative area fraction as the fractional area where $\omega > \omega_{max}/10$. A by-product of the scaling analysis is that the total vortex area $A$, which we diagnose with the dissipative area fraction, becomes independent of viscosity. It is given by

$$A \sim N\lambda^2 \sim \beta^{*-1}, \qquad (17)$$

which we verify in figure 6.

<figure><img src="image_7.png" /><figcaption>Figure 6: A: Snapshots of dissipative structures at two values of $\nu^*$. To create the binary images, we compare vorticity values to $\omega_{max}/10$. B Area of dissipative structures against $\nu^*$. Note how, although the structures become smaller (as shown in A), the total area remains constant (the prediction of area fraction from scaling analysis is $\beta^{*-1} = 10^{-2}$).</figcaption></figure>

## MOVIE LEGENDS

**Movie S1.** Comparison of spin-up of flow regimes corresponding to free-slip and no-slip boundary conditions at equal values of $\nu^* = 10^{-5}$ and $\beta^* = 100$. On the left the free-slip simulation develops inertial runaway, with an energetic recirculation at the domain scale. The no-slip simulation exhibits Gyre Turbulence, with energetic small scales.

**Movie S2.** Animation of Gyre Turbulence in the Vortex Gas regime ($\nu^* = 2 \times 10^{-6}$, $\beta^* = 100$, Numerical Resolution: 8192×8192 grid points). Boundary layer detachments occur at all boundaries, but the northwestern corner remains the most active region, mostly due to the persistence of the strong cyclonic (red) vortex. With regards to the scaling theory of the vortex gas, notice the event at the start of the animation, where the cyclonic (red) vortex collides with the western boundary and detaches the viscous sublayer.