# Surface quasigeostrophic turbulence: The refined study of an active scalar

Nicolas Valade, Jérémie Bec, Simon Thalabard

► To cite this version:

Nicolas Valade, Jérémie Bec, Simon Thalabard. Surface quasigeostrophic turbulence: The refined study of an active scalar. Journal of Fluid Mechanics, 2025, 1021, pp.A17. 10.1017/jfm.2025.10698. hal-05408199

HAL Id: hal-05408199

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

Submitted on 10 Dec 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
<figure><img src="image_2.png" /></figure>

# Surface quasi-geostrophic turbulence: the refined study of an active scalar

Nicolas Valade<fnref n="1" /> ID, Jérémie Bec<fnref n="1" />,<fnref n="2" /> ID and Simon Thalabard<fnref n="2" /> ID

**Corresponding authors:** Nicolas Valade, nicolas.valade@inria.fr; Jérémie Bec, jeremie.bec@univ-cotedazur.fr; Simon Thalabard, simon.thalabard@univ-cotedazur.fr

(Received 20 March 2025; revised 29 July 2025; accepted 3 September 2025)

Surface quasi-geostrophic (SQG) theory describes the two-dimensional active transport of a scalar field, such as temperature, which – when properly rescaled – shares the same physical dimension of length/time as the advecting velocity field. This duality has motivated analogies with fully developed three-dimensional turbulence. In particular, the Kraichnan – Leith – Batchelor similarity theory predicts a Kolmogorov-type inertial range scaling for both scalar and velocity fields, and the presence of intermittency through multifractal scaling was pointed out by Sukhatme & Pierrehumbert (2002 *Chaos* **12**, 439–450), in unforced settings. In this work, we refine the discussion of these statistical analogies, using numerical simulations with up to 16 384² collocation points in a steady-state regime dominated by the direct cascade of scalar variance. We show that mixed structure functions, coupling velocity increments with scalar differences, develop well-defined scaling ranges, highlighting the role of anomalous fluxes of all the scalar moments. However, the clean multiscaling properties of SQG transport are blurred when considering velocity and scalar fields separately. In particular, the usual (unmixed) structure functions do no follow any power-law scaling in any range of scales, neither for the velocity nor for the scalar increments. This specific form of the intermittency phenomenon reflects the specific kinematic properties of SQG turbulence, involving the interplay between long-range interactions, structures and geometry. Revealing the multiscaling in single-field statistics requires us to resort to generalised notions of scale invariance, such as extended self-similarity and a specific form of refined self-similarity. Our findings emphasise the fundamental entanglement of scalar and velocity fields in SQG turbulence: they evolve hand in hand and any attempt to isolate them destroys scaling in its usual sense.

<footnote n="1">¹Inria, CNRS, Calisto team, Université Côte d'Azur, Sophia Antipolis 06902, France</footnote>

<footnote n="2">²Université Côte d'Azur, CNRS, Institut de Physique de Nice, Nice 06200, France</footnote>
This perspective sheds new lights on the discrepancies in spectra and structure functions that have been repeatedly observed in SQG numerics for the past 20 years.

**Key words:** turbulent flows, intermittency, quasi-geostrophic flows

# 1. Introduction

The surface quasi-geostrophic (SQG) theory describes the dynamics of rapidly rotating and stably stratified flows near flat horizontal boundaries through the two-dimensional (2-D) transport equation

$$ \partial_t \theta + (\mathbf{u} \cdot \nabla) \theta = \nu \Delta \theta + F, \quad \mathbf{u} := -\nabla^{\perp} \psi, \quad \psi := (-\Delta)^{-\frac{1}{2}} \theta, \quad (1.1) $$

where $\theta(\mathbf{x}, t)$ represents a scalar field (such as surface temperature or buoyancy) advected in a planar domain $\Omega \subset \mathbb{R}^2$ under the influence of forcing $F(\theta, \mathbf{x}, t)$ and diffusivity $\nu$. Here, $\nabla^{\perp} = (-\partial_2, \partial_1)^{\top}$ tying the incompressible velocity field $\mathbf{u}$ to the stream function $\psi$ and $\Delta$ respectively denote the perpendicular gradient and the Laplacian operators.

Equation (1.1) describes the transport of an active scalar, where the velocity field $\mathbf{u}$ depends non-locally on the advected scalar $\theta$, analogous to the velocity – vorticity relationship in the 2-D Euler equations. In SQG, however, long-range interactions are mediated by the Green’s function $(-\Delta)^{-1/2}$ associated with the half-Laplacian rather than to the Laplacian. In the full plane, this yields

$$ \psi(\mathbf{x}) = \frac{1}{2\pi} \int_{\mathbb{R}^2} \frac{\theta(\mathbf{x}')}{|\mathbf{x} - \mathbf{x}'|} d^2\mathbf{x}' \quad \text{and} \quad \mathbf{u}(\mathbf{x}) = -\frac{1}{2\pi} \int_{\mathbb{R}^2} \frac{(\mathbf{x} - \mathbf{x}')^{\perp}}{|\mathbf{x} - \mathbf{x}'|^3} \theta(\mathbf{x}') d^2\mathbf{x}', \quad (1.2) $$

where $\mathbf{x}^{\perp} = (-x_2, x_1)$. These relations imply, somewhat surprisingly, that $\theta$ and $\mathbf{u}$ share the same physical dimension. This makes $\theta$ comparable to both a coarse-grained vorticity and a scalar analogue of the velocity field (see figure 1).

Initially proposed by Blumen (1978), the SQG model has emerged as a minimal yet insightful framework for studying the dynamics of upper-oceanic and lower-tropospheric flows, which involve a balance between strong rotation and strong stratification. It is particularly relevant to the problems involving potential vorticity inversion and frontogenesis (see, for example, Held et al. 1995; Juckes 1995; Hakim, Snyder & Muraki 2002; Lapeyre & Klein 2006; LaCasce & Mahadevan 2006 and the comprehensive review by Lapeyre 2017). At a more fundamental level, the SQG system (1.1) has also gathered significant interest in physical and mathematical turbulence research due to its phenomenological and structural analogies to various prototypical fluid systems. These include the 2-D Navier–Stokes (Held et al. 1995; Bernard et al. 2007), passive scalar turbulence (Celani et al. 2004), Burgers turbulence (Valade, Thalabard & Bec 2024) and the 3-D Navier–Stokes equations (Constantin, Majda & Tabak 1994; Sukhatme & Pierrehumbert 2002).

In the absence of forcing and dissipation, smooth solutions of the SQG equation conserve two quadratic invariants: the Hamiltonian $\mathcal{H} = (1/2) \int_{\Omega} \psi \theta$ akin to the energy in 2-D Euler flows, and the scalar variance, which also corresponds to the kinetic energy, $\mathcal{E} = (1/2) \int_{\Omega} \theta^2 = (1/2) \int_{\Omega} |\mathbf{u}|^2$ akin to the enstrophy. In the presence of forcing and dissipation, these invariants lead to a dual cascade scenario, featuring an inverse cascade of the Hamiltonian and a direct cascade of scalar variance (Smith et al. 2002; Burgess, Scott & Shepherd 2015).
<figure><img src="image_3.png" /><figcaption>Figure 1. SQG snapshots of squared vorticity $\omega = \nabla^{\perp} \cdot u$, temperature $\theta$ and velocity $u$ at a typical time. The zoomed-in regions display $(1/32)^2$ of the computational domain. Taken from run VI (see table 1).</figcaption></figure>

The direct transfer of scalar variance – and equivalently of kinetic energy – is known to differ fundamentally from the enstrophy cascade in 2-D turbulence, as it sustains more energetic small scales. In short, SQG turbulence is rough. Similar to 3-D turbulence, dimensional arguments tie the direct SQG cascade to a Kolmogorov energy spectrum proportional to $k^{-5/3}$, where $k$ denotes the wave number. However, previous numerical studies have reported both steeper spectral slopes (Sukhatme & Pierrehumbert 2002; Celani et al. 2004; Capet et al. 2008) and shallower ones (Watanabe & Iwayama 2007). These observations point toward the presence of significant corrections due to intermittency (Kolmogorov 1962; Frisch 1995; Benzi & Toschi 2023). One could also question the extent to which SQG scaling laws (if any) are universal. In particular, the recent work of Valadão et al. (2024) provided evidence that large-scale fluctuations can obscure a $k^{-5/3}$ scaling, with steeper subleading terms appearing in the spectrum.

In this work, we aim to examine in greater detail the Kolmogorovian features of steady-state SQG turbulence through a series of highly resolved direct numerical simulations, employing up to 16 384² collocation points. These simulations incorporate small-scale dissipation and a statistically homogeneous large-scale forcing scheme to sustain turbulence.

Because SQG is based on a large-scale geophysical balance, the physical interpretation of small-scale dissipation is better understood as a form of turbulent eddy diffusivity rather than molecular friction. Consequently, many studies involving SQG pragmatically replace the standard Laplacian dissipation with hyper-dissipation or *ad hoc* small-scale filtering to achieve extended inertial ranges even at moderate numerical resolution (Held et al. 1995; Smith et al. 2002; Celani et al. 2004; Watanabe & Iwayama 2007; Capet et al. 2008; Foussard et al. 2017; Jha, Seshasayanan & Dallas 2025). To draw parallels with 3-D turbulence and to pave the way for future Lagrangian studies, we here choose to rely
on standard Laplacian-mediated dissipation, following the approach of other fundamental SQG studies (Sukhatme & Pierrehumbert 2002; Valade et al. 2024; Valadão et al. 2024; Valadão et al. 2025).

Classical 3-D turbulence theory relies on the phenomenological assumptions that the steady-state velocity fields have finite variance and that the kinetic energy dissipation rate remains non-zero in the limit of vanishing viscosity – a property referred to as anomalous dissipation. For SQG turbulence, these conditions translate into the following requirements:

$$ \lim_{\nu \to 0} \langle \theta^2 \rangle < \infty \quad \text{and} \quad \lim_{\nu \to 0} \langle \nu | \nabla \theta |^2 \rangle > 0, \qquad (1.3) $$

where $\langle \cdot \rangle$ denotes a suitable space – time averaging. We propose a numerical set-up designed to check whether these conditions are satisfied. We point out that their realisation in SQG turbulence is not straightforward: The presence of Hamiltonian inverse transfers requires large-scale damping to establish a statistically steady state. Since interactions between scales in SQG might be non-local (see, e.g. Watanabe & Iwayama 2007; Foussard et al. 2017), this damping could influence the statistics of the direct cascade.

Our numerical simulations enable a refined analysis of the multiscaling properties of SQG. In 3-D turbulence, multiscaling is well captured by the Kolmogorov (1962) refined similarity hypothesis, which links the non-Gaussianity of velocity increments to fluctuations in the dissipation field (Chevillard et al. 2019; Dubrulle 2019). While evidence of intermittency in SQG turbulence has been reported in unforced settings (Sukhatme & Pierrehumbert 2002; Valade et al. 2024), its phenomenology remains poorly understood. Qualitatively, coherent structures play a pivotal role in the SQG dynamics across scales, generating sharp scalar gradients and complex velocity patterns entangling filaments, fronts and vortices, as illustrated in figure 1. These structures are likely coupled through cascades of shear instabilities (Scott & Dritschel 2014; Lapeyre 2017). As a result, although the scalar field $\theta$ and the velocity field $u$ share certain features – such as their variance, dissipation rate and dimensional scaling – they exhibit distinct statistical behaviours. For instance, Capet et al. (2008) observed that the longitudinal velocity increment $\delta u_\parallel(\ell)$ has a positive skewness, while $\delta\theta$ is symmetric. Moreover, the direct cascade of scalar variance imposes a negative correlation between $\delta u_\parallel$ and $(\delta\theta)^2$, stemming from a SQG-specific version of Yaglom's law (Yaglom 1949). As a specificity of SQG turbulence, these structural differences raise questions about its classification relative to two and three dimensions and about the quantification of intermittency. In particular, distinct scaling laws, if they exist, may apply to the scalar and velocity fields individually.

To address these questions, the paper is organised as follows. Section 2 introduces our numerical set-up and discusses the establishment of a statistically steady SQG state in a periodic domain. Section 3 investigates inertial-range scaling, with a particular focus on Yaglom's law and its implications for the statistical properties of both velocity and scalar fields. Section 4 explores intermittency by analysing higher-order structure functions and scalar dissipation statistics within the framework of the refined self-similarity hypothesis. Finally, § 5 summarises our findings and formulates concluding remarks. Some technical but otherwise standard aspects of SQG are compiled in the appendices.

## 2. Steady-state SQG turbulence in a doubly periodic domain

### 2.1. Surface quasi-geostrophic turbulence

#### 2.1.1. Steady-state balances

Our numerical simulations implement the SQG system (1.1) in the doubly $2\pi$-periodic domain, $\Omega = \mathbb{T}^2$, by decomposing the fields into their Fourier components. The Fourier
transform of the scalar field is defined as $\hat{\theta}_k = (1/(2\pi)^2) \int_{\mathbb{T}^2} \theta(\mathbf{x}) e^{-i\mathbf{k}\cdot\mathbf{x}} d^2\mathbf{x}$, such that $\theta(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^2} \hat{\theta}_k e^{i\mathbf{k}\cdot\mathbf{x}}$, and similarly for the components of the velocity field. In Fourier space, the velocity-scalar SQG relationship (Riesz transform) is then mediated through the explicit identities

$$ \hat{\mathbf{u}}_k = -i \frac{\mathbf{k}^{\perp}}{|\mathbf{k}|} \hat{\theta}_k \quad \text{and} \quad \hat{\theta}_k = i \frac{\mathbf{k}^{\perp}}{|\mathbf{k}|} \cdot \hat{\mathbf{u}}_k \quad \text{for } \mathbf{k} \in \mathbb{Z}^2 \setminus \{\mathbf{0}\}, \qquad (2.1) $$

where the null components are set to zero: $\hat{\mathbf{u}}_0 = \mathbf{0}$ and $\hat{\theta}_0 = 0$. From (2.1), it follows that the Fourier components of $\theta$ and $\mathbf{u}$ have the same amplitudes: $|\hat{\mathbf{u}}_k| = |\hat{\theta}_k|$ for all $\mathbf{k}$. As a result, they share the same power spectrum

$$ E(k) = \frac{1}{2} \sum_{k \le |\mathbf{k}| < k+1} |\hat{\theta}_k|^2 = \frac{1}{2} \sum_{k \le |\mathbf{k}| < k+1} |\hat{\mathbf{u}}_k|^2; \qquad (2.2) $$

as well as their spatially averaged square values and gradient norms

$$ f \theta^2 = f |\mathbf{u}|^2 \quad \text{and} \quad f |\nabla \theta|^2 = f |\nabla \mathbf{u}|^2, \qquad (2.3) $$

where we adopt the shorthand $f = \int_{\mathbb{T}^2} / (2\pi)^2$ for the spatial average over $\mathbb{T}^2$. To achieve a statistically steady state, we employ a stochastic forcing scheme $F(\theta, \mathbf{x}, t)$ that combines the superposition of white-in-time Gaussian random Fourier modes with a linear damping operating solely at the largest scales. In Fourier space, the forcing takes the form

$$ \hat{F}_k = \sqrt{\frac{2\mathcal{I}}{c_d \pi k_f^2}} e^{-\frac{|\mathbf{k}|^2}{2k_f^2}} \hat{\eta}_k - \alpha \hat{\theta}_k \mathbf{1}_{|\mathbf{k}| \le k_f}. \qquad (2.4) $$

Here, $\eta_k$ are random complex Gaussian coefficients, with statistics prescribed by $E[\hat{\eta}_k(t)] = 0$ and $E[\hat{\eta}_k(t) \hat{\eta}_k'(t')] = \delta_{k+k'} \delta(t-t')$, where $E[\cdot]$ denotes averaging over random realisations. The numerical factor $c_d \simeq 1.004$ accounts for discretisation effects (i.e. the transformation of integrals over wavenumber space into discrete sums over $\mathbf{k} \in \mathbb{Z}^2$ due to the toroidal geometry). It is defined such that the parameter $\mathcal{I}$ corresponds to the total energy injection rate: $(1/2) \sum_k |\hat{F}_k|^2 = \mathcal{I} - (\alpha/2) \sum_{|\mathbf{k}| \le k_f} |\hat{\theta}_k|^2$.

The damping coefficient $\alpha$ regulates the magnitude of large-scale contributions to the energy and the Hamiltonian. In the stationary state, by letting $\theta_f = \sum_{|\mathbf{k}| \le k_f} \hat{\theta}_k e^{i\mathbf{k}\cdot\mathbf{x}}$ and $\psi_f = \sum_{|\mathbf{k}| \le k_f} |\mathbf{k}|^{-1} \hat{\theta}_k e^{i\mathbf{k}\cdot\mathbf{x}}$, the following global balances for the energy and the Hamiltonian hold:

$$ \mathcal{I} = \varepsilon + \alpha \langle \theta_f^2 \rangle \quad \text{and} \quad k_{\mathcal{I}}^{-1} \mathcal{I} = \gamma + \alpha \langle \theta_f \psi_f \rangle, \qquad (2.5) $$

where $\varepsilon := \nu(|\nabla\theta|^2)$, $\gamma := \nu(\nabla\theta \cdot \nabla\psi)$ and the angular brackets $\langle \cdot \rangle$ denote averages over space, time and forcing realisations. In the Hamiltonian balance, we introduced the centroid wavenumber $k_{\mathcal{I}} = \tilde{c}_d k_f / \sqrt{\pi}$ (with $\tilde{c}_d \simeq 1.301$), which defines an injection scale.

### 2.1.2. Definition of SQG turbulence

Following the review by Boffetta & Ecke (2012) on the 2-D Navier-Stokes equations, we introduce the characteristic SQG damping and viscous centroid wavenumbers as, respectively,

$$ k_{\alpha} = \frac{\langle \theta_f^2 \rangle}{\langle \theta_f \rangle \langle \psi_f \rangle} \quad \text{and} \quad k_{\nu} = \frac{\varepsilon}{\gamma}. \qquad (2.6) $$
On phenomenological grounds, we expect $k_ν > k_I > k_α$, with the magnitudes of the scale separations defining steady-state regimes with different physical properties. In this work, we define SQG turbulence as a steady-state regime satisfying

$$k_ν \gg k_I \gtrsim k_α, \quad \text{with } k_I, k_α = O(1). \qquad (2.7)$$

Under this assumption, the steady-state balances (2.5), combined with the definitions (2.6), yield the following asymptotic estimates:

$$\varepsilon \simeq \left(1 - \frac{k_\alpha}{k_I}\right) \mathcal{I}, \quad \langle\theta_f^2\rangle \simeq \frac{k_\alpha}{\alpha k_I} \mathcal{I} \quad \text{and} \quad \gamma \simeq \frac{\varepsilon}{k_ν}. \qquad (2.8)$$

The estimates (2.8) fulfil the Kolmogorovian requirements (1.3) of anomalous scalar dissipation and ensure the finiteness of the scalar variance. Anticipating on the Kolmogorov 1941 (K41) phenomenology discussed in § 2.3.4, we expect the SQG spectra to have finite ultraviolet capacity, i.e. $\sum_{k>k_f} E(k) < \infty$. Equation (2.8) also implies that the viscous dissipation of the Hamiltonian vanishes in the limit $ν \to 0$. As such, our definition of SQG turbulence through (2.7) formalises a steady-state regime sustaining a direct cascade of scalar variance, with no inverse cascade. At this stage, SQG turbulence already differs from 3-D turbulence, as the energy dissipation rate does not perfectly balance the injected energy. The difference $\simeq \alpha \langle\theta_f^2\rangle = (k_\alpha/k_I)\mathcal{I}$ corresponds to energy leaking into the largest scales. This leakage introduces a degree of non-universality, but it comes as a necessary trade-off for the requirement that both $\alpha$ and $k_\alpha$ be finite, and prevents the scalar variance from blowing up.

### 2.1.3. A comment on SQG scales

From the SQG turbulence asymptotics (2.8), it is natural to define the SQG turnover time at wavenumber $k$ as $\tau := \varepsilon^{-1/3}k^{-2/3}$. We therefore define the time scales of injection $\tau_I$ and dissipation $\tau_ν$ as the turnover times at wavenumbers $k_I$ and $k_ν$, respectively. The dissipative scales $k_ν$ and $\tau_ν$ a priori differ from the usual Kolmogorov scales $k_η = (\varepsilon/ν^3)^{1/4}$, $\tau_η = (ν/\varepsilon)^{1/2}$ found in 3-D turbulence. We postpone further comments on this issue to § 2.3. We denote $\tau_α = 1/\alpha$ as the damping time scale. The SQG turbulence should correspond to a situation where $\tau_α \gg \tau_I$, in order to prevent large-scale fluctuations from dominating the turbulent statistics of the direct cascade (see Homann, Bec & Grauer 2013 for a related discussion on the 3-D Navier-Stokes equations). The conversion between physical length scales and Fourier wavenumbers is mediated through the convention $\ell = 2\pi/k$.

### 2.2. Numerical set-up

Our numerical simulations of the forced-dissipated SQG (1.1) in the doubly $2\pi$-periodic torus $\mathbb{T}^2$ employ a GPU-based pseudo-spectral solver. The solver uses the standard 2/3-rule for dealiasing (Orszag 1971), fourth-order Runge-Kutta time marching and up to $N^2 = 16\,384^2$ collocation points. The forcing scheme (2.4) uses $k_f = 3$, an injection rate $I = 1$, and a damping rate $\alpha = 0.1$, yielding an injection scale of $\ell_I = 2\pi/k_I \simeq 2.8$. To address finite-diffusivity effects, we define the Péclet number as

$$Pe := \frac{\ell_I \langle\theta_f^2\rangle^{1/2}}{\nu}. \qquad (2.9)$$

The Péclet number in SQG is expected to play a role analogous to the Reynolds number in 3-D turbulence (Sukhatme & Pierrehumbert 2002). Table 1 summarises the key
<table><caption>Table 1. Numerical settings. See the text for definitions. In all runs, we set the damping rate to $\alpha = 0.1$ and prescribe $k_f = 3$ and $\mathcal{I} = 1$, corresponding to an injection wavenumber $k_\mathcal{I} \simeq 2.2$.</caption><thead><tr><th>Runs</th><th>N<sup>2</sup></th><th>ν</th><th>Pe</th><th>k<sub>max</sub>/k<sub>η</sub></th><th>k<sub>max</sub>/k<sub>ν</sub></th><th>N<sub>real</sub></th><th>T<sub>max</sub>/τ<sub>I</sub></th><th>⟨θ<sub>f</sub><sup>2</sup>⟩</th></tr></thead><tbody><tr><td>I</td><td>1024<sup>2</sup></td><td>3.0 × 10<sup>-4</sup></td><td>2.3 × 10<sup>4</sup></td><td>1.0</td><td>28</td><td>1</td><td>30</td><td>7.4</td></tr><tr><td>II</td><td>2048<sup>2</sup></td><td>1.2 × 10<sup>-4</sup></td><td>5.8 × 10<sup>4</sup></td><td>1.0</td><td>30</td><td>1</td><td>30</td><td>6.8</td></tr><tr><td>III</td><td>4096<sup>2</sup></td><td>4.7 × 10<sup>-5</sup></td><td>1.5 × 10<sup>5</sup></td><td>1.0</td><td>33</td><td>1</td><td>300</td><td>7.0</td></tr><tr><td>IV</td><td>8192<sup>2</sup></td><td>1.9 × 10<sup>-5</sup></td><td>3.7 × 10<sup>5</sup></td><td>1.0</td><td>35</td><td>4</td><td>7</td><td>7.3</td></tr><tr><td>V</td><td>16 384<sup>2</sup></td><td>7.5 × 10<sup>-6</sup></td><td>9.0 × 10<sup>5</sup></td><td>1.0</td><td>36</td><td>4</td><td>7</td><td>6.8</td></tr><tr><td>VI</td><td>16 384<sup>2</sup></td><td>2.0 × 10<sup>-6</sup></td><td>3.5 × 10<sup>6</sup></td><td>0.4</td><td>7</td><td>4</td><td>7</td><td>6.7</td></tr></tbody></table>

parameters of our simulation series. Across the six runs, the injection and damping parameters remain fixed, while the Péclet number spans over two decades.

The infrared damping can be quantified in terms of a damped Péclet number: here, we define $Pe_\alpha := \langle \theta_f^2 \rangle^{1/2}/(\alpha \ell_\mathcal{I}) \simeq 10$, which remains constant across all simulations. The damping rate $\alpha = 0.1$ to which this corresponds was determined through trial and error, to ensure that large-scale dynamics evolves on time scales significantly slower than the injection time scale, specifically $\tau_\alpha \simeq 10 \simeq 4\tau_\mathcal{I}$, in line with § 2.1.3.

For runs I to V, the maximal wavenumber $k_{max} = N/3$ is chosen to satisfy $k_{max} \simeq k_\eta$, following good practices in numerical studies of 3-D turbulence (Grauer, Homann & Pinton 2012; Iyer, Sreenivasan & Yeung 2015; Buaria et al. 2019). This choice is relatively conservative in the context of SQG turbulence. As noted in § 2.1.3, the relevant dissipative scale for SQG is $k_\nu$, which is amply resolved for these runs. This motivates the parameters of run VI, trading a higher Péclet number for a less conservative scheme with $k_{max} \simeq 0.4k_\eta$.

Regarding simulation durations, runs I and II consist of a single realisation over $T_{max} = 10 \alpha^{-1} \simeq 30 \tau_\mathcal{I}$. Runs IV to VI gather statistics from four realisations, each initialised from different equilibrated states and subjected to distinct forcing realisations, with $T_{max} = 2.5 \alpha^{-1} \simeq 7 \tau_\mathcal{I}$. Run III corresponds to a single, much longer realisation, specifically designed to evaluate statistical convergence, as detailed in § 2.3.

## 2.3. Convergence towards SQG turbulence

We here assess the convergence of our numerical protocol towards SQG turbulence, by verifying that the conditions in (2.7) and (2.8) are satisfied as $Pe$ increases.

### 2.3.1. Steady-state and stationary statistics

To evaluate the relevance of our numerical protocol to reach a statistical steady state, we first focus on run III, which features the longest time series, extending up to 300$\tau_\mathcal{I}$. Figure 2 displays the typical temporal evolution of the energy $(1/2)\int \theta^2$ and the Hamiltonian $(1/2)\int \theta \psi$, as well as their instantaneous dissipation rates $\nu f |\nabla \theta|^2$ and $\nu f \nabla \theta \cdot \nabla \psi$, normalised by their time average values $\mathcal{E}$, $\mathcal{H}$, $\varepsilon$ and $\gamma$, respectively. It shows that all time series are statistically stationary, with well-defined averages and correlation time $\tau_c \sim 2\tau_\mathcal{I}$. With $\tau_\mathcal{I}$ emerging as the characteristic correlation time, it is justified to consider averages over durations of the order of tens of $\tau_\mathcal{I}$ as representative of long-time steady-state statistics.

The statistics also reveal a high degree of cross-correlations. The apparent synchronisation between the energy and the Hamiltonian reflect the global balance between injection and large-scale damping. In turn, the cross-correlation shown between
<figure><img src="image_4.png" /><figcaption>Figure 2. (a): time series of the SQG invariants and their corresponding dissipation rates for run III, normalised by their time-averaged values. (b): temporal self-correlations of scalar variance, Hamiltonian and cross-correlation of the energy and its dissipation rate. The label entries use the notations $\rho_{f,g} = \int_{\mathbb{R}^+} d\tau \langle f(\cdot)g(\cdot+\tau) \rangle / \langle fg \rangle$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. (a): averaged values of the Hamiltonian $\mathcal{H}$ and energy $\mathcal{E}$ as a function of $Pe$. (b): corresponding dissipation rates, $\gamma$ and $\epsilon$, respectively. All these quantities are here normalised by K41 dimensional predictions incorporating the appropriate powers of $\mathcal{I}$ and $\ell_{\mathcal{I}}$. The shaded areas highlight the minimum and maximum over the runs with $Pe > 5 \times 10^4$ and the error bars indicate twice the standard deviations of the statistical ensembles.</figcaption></figure>

the dissipation rate $\nu f |\nabla\theta|^2$ and the energy $(1/2)\int f\theta^2$ (inset of the right panel) reflects the direct cascade mechanism. The cross-correlation peaks at $\tau \approx 0.3\tau_{\mathcal{I}}$ corresponds to the typical cascading time from injection down to dissipation scales.

### 2.3.2. Scalar variance and dissipation

The convergence of SQG dynamics towards a developed turbulent state is examined in figure 3, which displays the steady-state values of the quadratic SQG invariants and their dissipation rates as a function of the Péclet number $Pe$. Despite significant statistical
<figure><img src="image_6.png" /><figcaption>Figure 4. (a): Péclet dependence of higher-order scalar moments. (b): corresponding diffusive dissipation rates normalised by $I^{(2q)} = q^2(\theta^{2q-2})\mathcal{I}$ – see (2.12). The shaded areas and the error bars are computed as in figure 3.</figcaption></figure>

fluctuations, our numerical results suggest that the energy $\mathcal{E}$ and the Hamiltonian $\mathcal{H}$ converge towards asymptotic values approximately given by

$$ \mathcal{E} \approx 2.7 \ell_{\mathcal{I}}^{2/3} \mathcal{I}^{2/3}, \quad \mathcal{H} \approx 0.4 \ell_{\mathcal{I}}^{5/3} \mathcal{I}^{2/3}. \qquad (2.10) $$

The distinct roles played by the energy and the Hamiltonian dissipation rates in establishing this asymptotic turbulent regime are evident from their respective dependencies on $Pe$. While $\gamma$ algebraically vanishes as $Pe$ increases, $\varepsilon$ asymptotes to a finite value

$$ \varepsilon \approx 0.34\mathcal{I}, \quad \gamma \propto Pe^{-3/4}. \qquad (2.11) $$

The $-3/4$ scaling for the Hamiltonian dissipation is consistent with basic phenomenological arguments (discussed below). For the energy dissipation, the observed asymptotic value of $\varepsilon$ implies, together with the steady-state balance (2.8), that $k_{\alpha} \approx 0.65k_{\mathcal{I}}$. This indicates that more than 60% of the injected energy remains at the large scales and is damped by the $\alpha$ term. It is worth noting that a slow (for instance logarithmic) decay of $\varepsilon$ as a function of $Pe$ cannot be entirely ruled out: the high level of statistical fluctuations in the dissipation field at the highest Péclet numbers produce high uncertainties on the average estimates.

### 2.3.3. Higher-order moments

Our simulations achieve a sufficiently strong form of statistical stationarity, where besides the scalar variance, the higher-order moments of $\theta$ also stabilise upon time averaging and asymptote to finite values as $Pe \to \infty$. This behaviour is illustrated in figure 4. These moments are associated with non-quadratic inviscid conservation laws (Casimirs), and are dissipated by both large-scale damping and diffusivity. The steady-state budget for $\langle \theta^{2q} \rangle$ is

$$ \varepsilon^{(2q)} := \nu |\nabla(\theta^q)|^2 = q^2 (\theta^{2q-2}) \mathcal{I} - \frac{\alpha q^2}{2q-1} \langle \theta^{2q-1} \theta_f \rangle. \qquad (2.12) $$
<figure><img src="image_7.png" /><figcaption>Figure 5. (a): energy spectra observed from the six runs; the black dashed line indicates the Kolmogorov prediction $C_K \epsilon^{2/3} k^{-5/3}$, with a Kolmogorov constant $C_K \simeq 0.35$. (b): Péclet number dependence of the viscous wavenumber $k_ν$ and Kolmogorov wavenumber $k_η$. See the text for definitions.</figcaption></figure>

Similarly to the quadratic case, the dissipation rates exhibit anomalous behaviours: one observes $\epsilon^{(4)} \simeq 1.9 \mathcal{I} \langle \theta^2 \rangle$, $\epsilon^{(6)} \simeq 4.8 \mathcal{I} \langle \theta^4 \rangle$, etc.

### 2.3.4. The K41 phenomenology in SQG turbulence

Figure 5 shows the power spectra from the six different runs. Kolmogorov 1941 similarity theory predicts the development of a scaling regime $E(k) \approx C_K \epsilon^{2/3} k^{-5/3}$ for $k_\mathcal{I} \ll k \ll k_\nu$, with an increasing range as Pe grows. The observed spectral slope is slightly steeper than $-5/3$, suggesting the presence of intermittency corrections, as discussed in § 1. As seen in figure 5(a), the viscous and Kolmogorov wavenumbers, $k_ν$ and $k_η$, share a similar scaling $\propto Pe^{3/4}$ with $k_η/k_ν \simeq 10$. This supports the interpretation of $k_ν$ as the dominant dissipation scale. To summarise, the SQG K41 picture is mediated by the identifications, holding up to possibly non-universal constant factors

$$ \epsilon \propto \mathcal{I}, \quad k_\nu \propto k_\eta, \quad \langle \theta_f^2 \rangle \propto \ell_{\mathcal{I}}^{2/3} \mathcal{I}^{2/3}. \qquad (2.13) $$

These identifications lead to $\gamma \simeq \epsilon/k_\nu \propto \mathcal{I}/k_\eta \propto Pe^{-3/4}$, consistent with the algebraic behaviour reported in figure 3.

## 3. Homogenous isotropic SQG turbulence

From the spectral picture of K41 theory, one might think that the statistics of $\theta$ and $u$ are similar. This, however, dismisses the presence of physical-space structures, which are already apparent at the level of one-point statistics. A more explicit distinction between the two SQG fields emerges when considering two-point statistics, which explicitly relate to cascade mechanisms. The SQG version of Yaglom's law ties a mixed third-order structure function to the SQG dissipation, revealing non-trivial correlations between $\theta$ and $u$. This exact law provides a natural framework for evaluating finite-diffusivity effects and identifying inertial-range statistics. After discussing the basic phenomenology of one-point statistics, we examine the validity of Yaglom's law in our numerical simulations and discuss its implications for the statistics of $\theta$ and $u$ – focusing in particular on the various third-order structure functions.
<figure><img src="image_8.png" /><figcaption>Figure 6. (a): probability density functions (PDF) of the scalar field $\theta$ (♦) and of the velocity component $u_1$ (★) standardised to zero mean and unit variance, based on the full space – time dataset. The black dashed line represents a Gaussian distribution, while the grey shaded areas show exponential tails $\propto e^{-2.4u_1}$. (b): time series of the spatially averaged flatnesses defined as $F_{\theta}(t) = f\theta^4/(f\theta^2)^2$, $F_{u_1}(t) = fu_1^4/(fu_1^2)^2$, along with representative fields configurations observed during two extreme fluctuations.</figcaption></figure>

### 3.1. One-point statistics

A basic phenomenology of SQG structures can already be grasped by comparing the one-point statistics of the scalar and velocity fields, as shown in figure 6 for run III. The statistics in figure 6 are computed over a finite time window $\simeq 300\tau_I$: both the scalar and the velocity components exhibit heavy tails, with an exponential decay indicative of intense non-Gaussian events – more pronounced for the velocity field. The extreme events in the two fields reflect different physical structures, which can be identified through the time series of the spatially averaged flatnesses, as shown in figure 6(b). At a qualitative level, large fluctuations in $\theta$ occur when strong vortices form, although these structures do not leave a clear footprint in the velocity fluctuations. By contrast, intense velocity fluctuations arise from the presence of long scalar fronts, along which the velocity becomes very large. In the inviscid case, such scalar filaments lead to logarithmic divergences in the velocity field – see for instance Juckes (1995) and Appendix A.1.

### 3.2. Yaglom's law

To investigate two-point spatial statistics, we define the increments of the velocity and scalar fields as $\delta u = u(x+\ell, t) - u(x, t)$ and $\delta\theta = \theta(x+\ell, t) - \theta(x, t)$. The longitudinal and transverse increments are denoted as $\delta u_\parallel = \delta u \cdot \hat{\ell}$ and $\delta u_\perp = \delta u \cdot \hat{\ell}^\perp$, with $\hat{\ell} = \ell/|\ell|$. Assuming statistical stationarity, homogeneity and isotropy, standard turbulence calculations yield an exact relation – Yaglom's law – for the mixed third-order structure function, $S_3(\ell) := \langle \delta u (\delta\theta)^2 \rangle$, which takes the form of the (vectorial) balance equation
<figure><img src="image_9.png" /><figcaption>Figure 7. (a): correlations involved in the correction terms of Yaglom’s law (3.1) obtained from run VI. (b): compensated structure function $S_3^{\parallel} = \langle \delta u_{\parallel}(\delta\theta)^2 \rangle / (-2\epsilon\ell)$, as a function of the normalised separation $\ell/\ell_v$ for the various runs. The black dotted line represents the viscous correction $1 - \mathcal{E}_v(\ell)$ in (3.1), measured for run VI. The diamonds indicate the perpendicular contribution to Yaglom’s law.</figcaption></figure>

$$ S_3(\ell) = -2\epsilon\ell [1 - \mathcal{E}_f(\ell) - \mathcal{E}_v(\ell)], \quad \mathcal{E}_{f,v}(\ell) := \frac{2}{\ell^2} \int_0^\ell \ell' C_{f,v}(\ell') d\ell' \\ \text{with } C_f(\ell) := \frac{\alpha}{\epsilon} [\langle \theta(\mathbf{0})\theta_f(\mathbf{0}) \rangle - \langle \theta(\mathbf{0})\theta_f(\ell) \rangle] \quad \text{and} \quad C_v(\ell) := \frac{\nu}{\epsilon} \langle \nabla\theta(\mathbf{0}) \cdot \nabla\theta(\ell) \rangle. \quad (3.1) $$

Figure 7(a) shows the behaviour of these correlation terms. The inertial range is defined as the range of scales where the correction terms $\mathcal{E}_f$ and $\mathcal{E}_v$ are negligible, that is $\lambda\nabla\theta \ll \ell \ll \lambda_{\theta_f}$, where the lower bound is set by the correlation length of the scalar gradient and the upper bound by that of the filtered field $\theta_f$. In practice, figure 7 reveals that these correlation lengths are of the order of $\ell_v$ and $\ell_I$, respectively. Within the inertial range, (3.1) prescribes the scaling of the third-order longitudinal and transverse mixed structure functions

$$ S_3^{\parallel}(\ell) := \langle \delta u_{\parallel}(\delta\theta)^2 \rangle = -2\epsilon\ell \quad \text{and} \quad S_3^{\perp}(\ell) := \langle \delta u_{\perp}(\delta\theta)^2 \rangle = 0. \quad (3.2) $$

Figure 7(b) shows the longitudinal component of the structure function $S_3^{\parallel}$, compensated by $-2\epsilon\ell$ for all six runs. The results reveal convergence towards an asymptotic profile, with a plateau in the range $\ell \gg \ell_v$ that expands as the Péclet number $Pe$ increases. At small scales $\ell \lesssim 0.5\ell_v$, the profile is dominated by viscous corrections, captured by $1 - \mathcal{E}_v(\ell)$. At scales $\ell \gtrsim 20\ell_v \simeq 0.1\ell_I$, the contribution from the large-scale damping can no longer be neglected, with deviations from Yaglom’s scaling being at least of the order of 10%. In summary, our numerical results confirm the validity of Yaglom’s law in the range $0.5\ell_v \lesssim \ell \lesssim 20\ell_v$, providing a practical identification of the inertial range in SQG turbulence. For runs V and VI, which have the highest resolution ($N^2 = 16\,384^2$), this range extends over slightly more than one decade.

### 3.3. Third-order structure functions

While the mixed structure function $S_3$ prescribes a negative correlation $\propto \ell$ between the scalar and the velocity field, the behaviour of the third-order structure functions associated with each individual field are different. The first distinction arises from symmetry considerations. In homogeneous isotropic SQG turbulence, both the scalar increment $\delta\theta$ and the transverse velocity increment $\delta u_{\perp}$ have symmetric distribution.
<figure><img src="image_10.png" /><figcaption>Figure 8. (a): skewness of the longitudinal velocity increments for run VI. (b): third-order parallel structure functions for the various runs, compensated by Yaglom's scaling $2\epsilon\ell$. The shaded area on the left is indicative of the inertial range for the highest-resolution run – see text for details.</figcaption></figure>

This implies $S_3^\theta = \langle (\delta\theta)^3 \rangle = 0$ and $S_3^{u\perp} = \langle (\delta u_\perp)^3 \rangle = 0$ (see Appendix A.2.1 for details). A more subtle observation concerns the statistics of longitudinal velocity increments, which exhibit positive skewness over the SQG inertial range (see figure 8(a)). Given that both $u$ and $\theta$ share the same direct cascade spectrum, this positive skewness is somewhat counterintuitive. Another intriguing feature is that Yaglom's scaling does not hold for the amplitude of $S_3^{u\parallel} = \langle (\delta u_\parallel)^3 \rangle$. No clear plateau emerges in figure 8(b), where $S_3^{u\parallel}$ is shown for runs III to VI, normalised by Yaglom's scaling $2\epsilon\ell$.

The fact that the third-order longitudinal velocity structure function $S_3^{u\parallel}$ deviates from the dimensional linear scaling $\propto \ell$ highlights a key specificity of SQG turbulence: the velocity $u$ and scalar field $\theta$ are distinct entities with non-trivial correlations. These differences can be traced back to the specific kinematic structure of SQG, which will be further commented on in § 3.4.

To quantify deviations from linear scaling, we use hereafter the local scaling exponents

$$\xi_3^{u\parallel}(\ell) = \frac{d \log S_3^{u\parallel}}{d \log \ell}. \qquad (3.3)$$

Figure 9(a) reveals that in the SQG scaling range, all our runs exhibit a monotonous decay of $\xi_3^{u\parallel}$ from 1.8 to 1.1, indicating the absence of a well-defined power-law scaling range. For run VI, we cannot entirely rule out a power-law behaviour with exponent $\xi_3^{u\parallel} \simeq 1.1$ over the range $5\ell_v \lesssim \ell \lesssim 50\ell_v$, but this barely overlaps with the previously defined SQG inertial range. One interpretation is that the discrepancy comes from finite-diffusivity effects, which may influence $S_3^{u\parallel}$ over a broader range than for $S_3^\parallel$. Another interpretation is that the apparent scaling range close to the forcing scale is simply a visual artefact. It would then reflect a large-scale inflection point due to the finite damping coefficient that we used to maintain steady-state statistics. While we lean towards the second interpretation, providing a clear-cut answer would require extending our study to tremendously large resolutions.

As a final comment, please observe that the skewness of the longitudinal velocity vanishes at small scales. This is a consequence of incompressibility and isotropy playing out in a 2-D domain, leading to $\partial_1 u_1 = -\partial_2 u_2 \stackrel{\text{law}}{\sim} -\partial_1 u_1$ (see also Appendix A.2.2).
<figure><img src="image_11.png" /><figcaption>Figure 9. (a): scale-dependent scaling exponents for the third-order structure functions, determined through logarithmic derivatives. (b): ageostrophic term. The solid lines represent Bernard's balance (3.5). The shaded areas indicate the inertial range identified in figure 8.</figcaption></figure>

### 3.4. The kinematic nature of the skewness phenomenon

In both 2-D and 3-D Navier-Stokes turbulence, the skewness of the longitudinal velocity increments reflect the cascade directionality – negative for the direct cascades regime of energy (in three dimensions) and enstrophy (two dimensions), positive for the 2-D inverse cascade. However, this direct correspondence breaks down in SQG. While figure 8 shows that the SQG velocity increments are positively skewed, they do not reflect any type of inverse cascade. As this seemingly contradicts the earlier interpretation of Capet et al. (2008), let us examine this effect in greater detail.

In our setting of SQG turbulence, the energy flux is fully prescribed by the mixed structure function $S_3^{\parallel}$, which is negative and follows Yaglom's law (3.1). The velocity structure function $S_3^{u\parallel}$ in turn does not enter by itself any simple scale-by-scale energy budget, and is necessarily compensated by an ageostrophic contribution that counterbalances its sign (see Capet et al. 2008). This contribution arises from the integral identity (1.2) tying $\theta$ to $u$, which provides a purely kinematic relation between the unmixed structure function ($S_3^{u\parallel}$) and the mixed one ($S_3^{\parallel}$).

This relation is mediated through a formula analogous to one pointed out by Bernard (1999) for 2-D Navier-Stokes equations. The formula is most conveniently written in complex notations, promoting in particular the vorticity $\omega = \nabla^{\perp} \cdot u$ and the scalar $\theta$ to functions $\Omega, \Theta$ defined on the complex plane $\mathbb{C}$. Writing $z = x_1 + ix_2$, $U = u_1 + iu_2$ and $\delta U = U(z + \xi) - U(z)$ the (complex) increment over the separation $\xi \in \mathbb{C}$, the SQG version of Bernard's formula reads

$$ \frac{1}{3} \partial_{\xi\xi} \langle (\delta U)^3 \rangle + \langle A(0)\Omega(\xi) \rangle = \partial_{\xi\xi} \ast \langle \delta U (\delta\Theta)^2 \rangle. \quad (3.4) $$

Bernard's formula involves the complex field $A = a_1 + ia_1$, associated with the vector field $a := (-\Delta)^{1/2}[u\theta] - u\omega$.

Although not immediately evident, the field $a$ captures deviations from the leading-order geostrophic balance and thus represents an ageostrophic component. It also admits a physical interpretation in terms of vorticity stretching. These connections are developed further in § 3.5.
Equation (3.4) provides a decomposition of the mixed structure function $\langle \delta U(\delta\Theta)^2 \rangle$ into a single-field contribution $\langle (\delta U)^3 \rangle$ and a correction term $\langle A(0)\Omega(\xi) \rangle$. It is a physical-space counterpart of a decomposition proposed by Capet et al. (2008) for the flux of scalar variance. We stress out that (3.4) is of purely kinematic nature – only relying on the definition of $\theta$, $\omega$, $u$. To couple it with the SQG dynamics, we use Yaglom’s law (3.1).

In the complex setting, Yaglom’s law (3.1) entails the inertial range relation $\langle \delta U(\delta\Theta)^2 \rangle = -2\epsilon\xi$, implying that (3.4) there reduces to $\partial_{\xi\xi}\langle (\delta U)^3 \rangle = -3\langle A(0)\Omega(\xi) \rangle$. In words, the third-order velocity statistics are blind to the direct cascade, in the sense that they do not relate to anomalous dissipation but are only prescribed by the subleading terms featured in (3.1).

The argument can be freed from the complex setting, upon tediously observing that (3.4) implies a balance between longitudinal structure functions, which takes the explicit differential form

$$ \mathcal{L}_{u\parallel} S_3^{u\parallel} = \mathcal{L}_{\parallel} S_3^{\parallel} + \langle a_{\parallel}(\ell e_1) \omega(0) \rangle, $$

$$ \text{with } \mathcal{L}_{u\parallel} f := -\frac{1}{12} \ell^{-2} \partial_{\ell} [\ell^{-1} \partial_{\ell} [\ell^5 \partial_{\ell} (\ell^{-1} f)]], \quad \mathcal{L}_{\parallel} f = \frac{1}{4} \partial_{\ell} [\ell^{-1} \partial_{\ell} (\ell f)], \quad (3.5) $$

naturally involving the vorticity $\omega$ and the ageostrophic vector field $a$. Again, we observe $\mathcal{L}_{\parallel}[-2\epsilon\ell] = 0$, so that the operator $\mathcal{L}_{\parallel}$ annihilates the leading-order term appearing in Yaglom’s law. In the inertial range, the scale-by-scale balance (3.5) therefore relies on the ageostrophic contribution in the right-hand side of (3.5) and the subleading term featured in Yaglom’s law. This is confirmed numerically in figure 9(b), showing $\langle a_{\parallel}(\ell e_1)\omega \rangle$ and $\mathcal{L}_{u\parallel} S_3^{u\parallel}$ matching one another within the inertial range. Besides, the ageostrophic contribution consistently departs from the dimensional scaling $\propto 1/\ell$, with an apparent log-corrected behaviour

$$ \langle a_{\parallel}(\ell e_1)\omega \rangle \sim \frac{\epsilon}{\ell} \log \frac{\ell_a}{\ell}, \quad \ell < \ell_a \approx 10\ell_v \quad (3.6) $$

independent of the Péclet number. This substantiates the idea that the absence of scaling in $S_3^{u\parallel}$ arises from the contribution of $\mathcal{L}_{u\parallel}^{-1}[\langle a_{\parallel}(\ell e_1)\omega \rangle]$, as evidenced by the local slopes in figure 9(a). Essentially, it shows that the velocity skewness in SQG originates from the kinematics imposed by the presence of the two fields $u$ and $\theta$, while the anomalous energy flux appears inconsequential. We refer the reader to Appendix B for the derivations of (3.4) and (3.5).

### 3.5. Vorticity stretching and ageostrophy

The SQG framework is obtained from the Boussinesq equations with a first-order expansion in Rossby number (Pedlosky 2013; Vallis 2017). Although this is implicit in the scalar equation – which involves only geostrophic quantities – first-order corrections appear in the dynamics of the velocity and vorticity, and represent departures from geostrophic balance. The field $a$ is a particular example of such ageostrophic corrections; it endows SQG with 3-D-like features.

One effect is physical: $a$ is involved in vorticity stretching. This becomes clear when examining the equation for vorticity $\omega = \nabla^{\perp} \cdot u = (-\Delta)^{1/2}\theta$, directly derived from the SQG dynamics (1.1) as

$$ \partial_t \omega + u \cdot \nabla \omega = -\nabla \cdot a + \nu \Delta \omega + (-\Delta)^{1/2} F, \quad a := (-\Delta)^{1/2}[u \theta] - u \omega. \quad (3.7) $$

Here, the vector field $a$ appears under a divergence operator, meaning that it is only defined up to a rotational component (i.e. up to the addition of a divergence-free vector field).
<figure><img src="image_12.png" /><figcaption>Figure 10. Snapshots of the scalar field $\theta(a)$, vorticity $\omega = \nabla^{\perp} \cdot u$ (b) and the divergence of the ageostrophic term $\nabla \cdot a$ (c), normalised to unit variance and zoomed over $(1/16)^2$ of the computational domain. (Data from run V.)</figcaption></figure>

With the presence of a non-vanishing $\nabla \cdot a$ term on its right-hand side, unrelated to forcing and dissipation, (3.7) makes apparent the difference between two dimensions and SQG. In particular, and unlike in 2-D Navier-Stokes turbulence, the enstrophy $(1/2)\omega^2$ is not conserved in SQG, and its evolution is governed by

$$ \frac{d}{dt} \frac{1}{2} \int \omega^2 = \int a \cdot \nabla \omega - \nu \int |\nabla \omega|^2 + \int \omega (-\Delta)^{1/2} F, \quad (3.8) $$

where the term $\int a \cdot \nabla \omega = -\int \omega \nabla \cdot a$ acts as a non-viscous source or sink of enstrophy. This is the direct analogue of the vortex stretching term in 3-D Navier-Stokes. Figure 10 illustrates the non-trivial correlations between $\nabla \cdot a$ and the SQG fields. In particular, $\nabla \cdot a$ exhibits strong values along the sharp fronts of a scalar filament, while remaining nearly zero elsewhere.

Another effect is geophysical, as $a$ is directly connected to vertical motions. Recall that the SQG framework governs the evolution of the surface buoyancy $\theta$ and serves as the boundary condition for 3-D potential vorticity inversion problem. Specifically, in the upper ocean, the streamfunction $\psi$ satisfies a Poisson equation with mixed (Robin) boundary conditions, namely

$$ (\partial_{11} + \partial_{22} + \partial_{33}) \psi = 0, \quad \partial_3 \psi|_{x_3=0} = \theta, \quad \psi \xrightarrow{x_3 \to -\infty} 0. \quad (3.9) $$

This system arises as a first-order perturbation of the primitive equations under the assumption of vanishing interior potential vorticity (Held et al. 1995). Comparing (3.7) with formula (3.a) in Held et al. (1995), we identify $\nabla \cdot a = -\partial_3 u_3|_{x_3=0}$.

In other words, and up to a rotational component, $a$ represents an ageostrophic correction to the horizontal velocity field.

Its regions of divergence (convergence) correspond to upward (downward) vertical motions – therefore explicitly related to 3-D motion.

## 4. Intermittency

The distinct behaviours of the third-order structure functions indicate that the SQG fields $\theta$ and $u$, while non-trivially correlated, exhibit different scaling properties. They should be treated as separate statistical entities. Besides, the apparent lack of systematic scaling raises the question of how to quantify intermittency in SQG turbulence. In this section, we report results on intermittency beyond third-order statistics, employing three different
<figure><img src="image_13.png" /><figcaption>Figure 11. (a): probability density function of the scalar (green), longitudinal velocity (blue) and perpendicular velocity (red) increments for various separations, standardised to zero mean and unit variance, from run VI. Curves are shifted vertically for clarity. (b): corresponding flatnesses as a function of the separation $\ell$. The vertical arrows indicates scales at which PDFs of the left panel were evaluated. The shaded area indicates the inertial range identified in figure 8.</figcaption></figure>

but consistent approaches to extract SQG anomalous scaling exponents, based on extended self-similarity, mixed structure functions and refined similarity.

## 4.1. *Observations of intermittency*

### 4.1.1. *The PDFs and flatness*

Figure 11(b) shows the PDFs of standardised increments (zero mean, unit variance) of the scalar field, transverse velocity and longitudinal velocity from run VI. The distributions correspond to increment lengths spanning from the inertial range down to the grid resolution ($\Delta x = 0.05 \ell_v$ for run VI). The scale dependence of the statistics provides clear evidence of intermittency, and highlights differences between the fields. For large separations, the PDF of velocity increments have exponential tails, with the longitudinal component $\delta u_\|$ displaying noticeable skewness.

As $\ell$ decreases, the three PDFs flatten, with their tails becoming closer to stretched exponentials. At the smallest scales $\ell \sim \Delta x$, the three distributions nearly collapse towards a symmetric profile, indicating that the gradient statistics of all three quantities are statistically alike (and unskewed). This identification is captured visually in figure 12, where only a trained eye could distinguish between the squared gradients of the velocity and scalar fields. The widening of the PDF is quantified in particular by the flatnesses $F^\phi(\ell) = \langle \delta\phi^4 \rangle / \langle \delta\phi^2 \rangle^2 (\phi = \theta, u_\| \text{ and } u_\perp)$, which continuously increase as $\ell \to 0$, as seen in figure 11(b).

A scaling behaviour $F^{u_\perp} \propto \ell^{-0.2}$ can be quite clearly identified for the perpendicular velocity increments within the inertial range. The flatnesses of the parallel velocity increments $F^{u_\|}$, and of scalar increments $F^\theta$ are compatible with this trend but show less well-defined scaling ranges.

These observations indicate that SQG intermittency extends beyond the third-order statistics discussed in § 3.

Unlike $S_3$, we are not aware of any kinematic relations linking the fourth-order statistics of the various fields. Consequently, the scaling behaviour $\propto \ell^{-0.2}$ for the flatness is to be considered as an empirical observation.
<figure><img src="image_14.png" /><figcaption>Figure 12. Typical snapshots of the squared velocity gradient (a) and scalar gradient (b), normalised to unit variance, for run VI.</figcaption></figure>

### 4.1.2. The lack of a (usual) scaling range

In turbulence studies, scaling exponents are usually identified through the presence of an inertial range scaling for the structure functions of order $p$

$$S_p^{\phi} = (|\delta\phi|^p) \propto \ell^{\zeta_p^{\phi}}, \quad \delta\phi := \phi(\mathbf{x} + \ell\mathbf{e}_1) - \phi(\mathbf{x}), \quad (4.1)$$

associated with a prescribed scalar field $\phi$. In SQG turbulence, the presence of well-defined scaling exponents is, however, not a generic feature of single-field statistics. This feature is worth emphasising. Figure 13 presents the behaviour with the separation $\ell$ of the local slopes $\zeta_p^{\phi}(\ell) = d \log S_p^{\phi}(\ell)/d \log \ell$ for the second- and fourth-order structure functions associated with the three SQG fields ($\phi = \theta, u_{\parallel}$ and $u_{\perp}$). A (usual) scaling range in the sense of (4.1) should translate into a plateau for the local exponent over the inertial range.

The main panels show the results for run VI. At both orders, the local scaling exponents decrease monotonically with increasing scale within the inertial range – This trend is particularly pronounced for the longitudinal velocity increments. The main observation is the absence of a well-defined plateau for $\ell \lesssim 20\ell_v \simeq 0.1\ell_{\tau}$, suggesting that no inertial-range power-law scaling regime holds for any of the three SQG fields. The insets allow us to assess the Pe-dependency. They superpose the local scaling exponents $\zeta_2^{u\parallel}(\ell)$ and $\zeta_4^{u\parallel}(\ell)$ for runs III to VI. In both cases, the scaling exponents collapse at small scales, revealing a non-trivial scale-dependent master curve. The small plateau observed at intermediate scales can presumably be attributed to finite-size effects: It does not widen with increasing
<figure><img src="image_15.png" /><figcaption>Figure 13. (a): scale-dependent second-order exponents, $\zeta_2^\phi(\ell)$, for the SQG fields $\phi = \theta$, $u_\parallel$, $u_\perp$ from run VI. Inset: convergence of $\zeta_2^{u_\parallel}(\ell)$ with increasing Péclet number Pe. (b): same as the left panel but for the fourth-order structure function. The dashed lines in the insets indicate the K41 values (2/3 and 4/3). The shaded areas indicate the inertial range identified in figure 8.</figcaption></figure>

Pe and its value varies with the Péclet number. Note that those apparent asymptotic values do not coincide with the K41 similarity exponents, nor are bounded in any way by the latter.

On the one hand, the measurements shown in figure 13 suggest that neither the scalar field nor the velocity field exhibits power-law scaling in the usual sense. In particular, this suggests that the energy spectrum slope is not a well-defined quantity: this could explain the difficulties reported for instance by Lapeyre (2017) at consistently determining scaling exponents in SQG turbulence. On the other hand, both the constant skewness in figure 8 and the flatness scaling $\propto \ell^{-0.2}$ of figure 11 are examples of self-scaling properties (Benzi et al. 1996), reflecting a weaker form of scale-invariance through the presence of a relative scaling between different orders. Those self-scaling properties, together with the behaviour of specific mixed structure functions, provide a robust framework to characterise intermittency in SQG turbulence. This is what we discuss next.

## 4.2. Scaling exponents in SQG

We now argue that SQG intermittency is structured by the transport properties of SQG turbulence, which in particular prescribe scaling laws for mixed structure functions. Similarly to the third-order case, single-field statistics are not straightforwardly prescribed by these mixed structure functions. We, however, describe two ways to extract scaling exponents from single-field statistics: extended self-similarity (ESS) and refined similarity (RSS). Both methods provide scaling exponents of order-$p$ consistently deviating from K41 similarity scaling $p/3$. The quadratic corrections are captured by Kolmogorov–Obukhov log-normal dissipation model $\zeta_p \approx p/3 + \mu_{62}(p/3)(1-p/3)$ with $\mu_{62} = 0.16 \pm 0.02$ up to order $p=7$.

### 4.2.1. *Mixed structure functions*

The transport dynamics provides the most natural way to identify power-law scaling in the usual sense, reflecting the behaviour of mixed structure functions connected to conservation laws. The special case $\langle \delta u_\parallel \delta \theta^2 \rangle \propto \ell$ associated with inertial-range constancy of the energy flux was extensively discussed in § 3, in connection to Yaglom’s law.
<figure><img src="image_16.png" /><figcaption>Figure 14. (a): local scaling exponents $\xi_p^\parallel$ associated with the mixed structure functions $\langle |\delta u_\parallel| |\delta \theta|^{p-1} \rangle$ for run VI. (b): same as the left panel but for the Watanabe – Gotoh mixed structure functions $\langle |\delta u_\parallel| \delta \theta^2|^{p/3} \rangle$. In both panels, we show integer $p \in [2, 7]$ (from blue/bottom to red/top). Dash lines indicate their inertial-range mean values and shaded areas their standard deviations.</figcaption></figure>

The mixed structure functions associated with inviscid conservation laws of higher orders are defined by

$$S_p^\parallel = \langle |\delta u_\parallel| |\delta \theta|^{p-1} \rangle. \qquad (4.2)$$

For odd $p = 2q + 1 > 1$, $S_p^\parallel$ are involved in the scale-by-scale budgets of the scalar moments $(1/2q)\langle \theta^{2q} \rangle$, corresponding to conservation laws of inviscid SQG – Yaglom’s law corresponds to $q = 1$. Figure 14(a) shows that the $S_p^\parallel$’s exhibit scaling in the usual sense, namely

$$S_p^\parallel \propto \ell^{\xi_p^\parallel}. \qquad (4.3)$$

For $p \in [2, 7]$, the local scaling exponents measured in run VI are nearly constant within an inertial scaling range increasing with $p$, hereby allowing for a faithful measurement of $\xi_p^\parallel$.

This type of scaling laws is similar to those observed in passive scalar turbulence (Gauding, Danaila & Varea 2017). In passive scalar turbulence, however, Watanabe & Gotoh (2004) also identified scaling ranges for the family of structure functions constructed as $S_p^{WG} = \langle |\delta u_\parallel| \delta \theta^2|^{p/3} \rangle$. This observation does not extend to SQG. Figure 14(b) shows the local scaling exponent of the Watanabe – Gotoh structure functions $S_p^{WG}$, exhibiting monotonous decrease within the inertial range. In particular, for $p = 3$, the mixed structure function $\langle |\delta u_\parallel| \delta \theta^2 \rangle$ does not grow linearly with $\ell$, nor do any of its powers, and this contrasts with $S_3^\parallel = \langle |\delta u_\parallel| \delta \theta^2 \rangle \propto \ell$.

This discrepancy shows that the precise definition of structure functions is crucial to extract scaling behaviour. It also highlights the role of transport: scale-by-scale transport is mediated by the (signed) longitudinal velocity increment $\delta u_\parallel$, which in turn produces both the scaling range and sign cancellations in the structure functions.

### 4.2.2. Extended self-similarity

One way to extract power laws directly from single fields statistics is given by ESS (Benzi & Toschi 2023). It consists in expressing structure functions in terms of a lower-order one in order to extend the range of scales over which power-law behaviour is
<figure><img src="image_17.png" /><figcaption>Figure 15. Extended self-similarity for integer orders $p \in [1, 7]$ (from blue to red symbols) for the three SQG fields (left: $u_\parallel$, centre: $\theta$, and right: $u_\perp$). The coloured areas represent the values obtained in § 4.2.1 using mixed structure functions. The bottom right panel shows the scale-dependent exponent $\xi_p^{u_\parallel}$ normalised by $C_p$. The shaded area there indicates the inertial range identified in figure 8.</figcaption></figure>

observed. This strategy was introduced by Benzi et al. (1993a,b) and proves remarkably efficient to extract scaling exponents in 3-D homogeneous isotropic turbulence, even at relatively low Péclet numbers. Instead of looking for pure power laws, ESS assumes a self-scaling property for the structure functions, namely

$$S_p^{\phi} \propto (S_2^{\phi})^{C_p^{\phi}} \quad (\phi = \theta, u_{\parallel}, u_{\perp}), \qquad (4.4)$$

with the exponents $C_p^{\phi}$ a priori different for each SQG field. The presence of scaling ranges in the usual sense (4.1) implies ESS with $C_p^{\phi} = \xi_p^{\phi}/\xi_2^{\phi}$. The lack of usual scaling ranges does not, however, preclude the self-scaling property (4.4). For example, assuming a generalised scaling of the form $S_p^{\phi}(\ell) \propto [\ell f(\ell)]^{\xi_p^{\phi}}$, ESS would still hold, with the exponent $C_p^{\phi} = \xi_p^{\phi}/\xi_2^{\phi}$ characterising the dominant power-law contribution (see Benzi et al. 1993a).

The ESS analysis of our SQG runs is summarised in figure 15, showing the scale-dependent ESS exponents $C_p^{\phi}(\ell) = \xi_p^{\phi}(\ell)/\xi_2^{\phi}(\ell)$ at integer orders $p \in [1, 7]$, constructed from the scale-dependent exponents of the structure functions. The presence of an ESS scaling range translates into $C_p^{\phi}(\ell)$ being constant over an inertial range of scales. While the scalar increments exhibit no clear ESS scaling range, the velocity structure functions do. The local slopes of both longitudinal and transverse structure functions fluctuate
around constant values $C_p^{u\parallel}$ and $C_p^{u\perp}$ over a broad range of scales $0.1\ell_v \lesssim \ell \lesssim 40\ell_v$, extending beyond the range over which Yaglom's law is observed. figure 15(d) points out that when normalised by their ESS exponent $C_p^{u\parallel}$, the scale-dependent exponents $\xi_p^{u\parallel}(\ell)$ of the velocity field collapse to a profile which is largely independent of the order: this is another way to substantiate the relevance of the ESS relation (4.4), showing that the quantities $(S_p^\phi)^{1/C_p^\phi}$ are independent of $p$.

### 4.2.3. Refined similarity

The RSS theory of Kolmogorov (1962) for 3-D turbulence ties deviations from self-similarity to heterogeneities in the scaling behaviours for the coarse-grained dissipation field

$$ \varepsilon_{\ell}(\mathbf{x}, t) := \frac{4}{\pi \ell^2} \int_{|\mathbf{r}|<\ell/2} \varepsilon(\mathbf{x} + \mathbf{r}) d^2\mathbf{r}, \quad \text{with } \varepsilon(\mathbf{x}, t) = \nu |\nabla\theta|^2. \quad (4.5) $$

Refined self-similarity suggests the phenomenology $\langle\delta\phi^p\rangle \sim \ell^{p/3}\langle\varepsilon_{\ell}^{p/3}\rangle$. If the dissipation field is multifractal, the moments of the coarse-grained field should behave as power laws $\langle\varepsilon_{\ell}^q\rangle \propto \ell^{\chi_q}$ (see, e.g. Frisch 1995), so that

$$ \langle\delta\phi^p\rangle \propto \ell^{\zeta_p^\phi} \quad \text{with} \quad \zeta_p^\phi = p/3 + \chi_{p/3}. \quad (4.6) $$

In SQG, however, such bridging relations do not hold – they would imply scaling in the usual sense, which we recall is not observed (see § 4.1.2). Still, our numerical simulations suggest that the dissipation field is multifractal and that the moments of dissipation behave as power laws: this allows us to use (4.6) to define the RSS scaling exponents as $\zeta_p^\phi = p/3 + \chi_{p/3}$. Figure 16(a) shows the local similarity exponents $p/3 + \chi_{p/3}(\ell)$, with $\chi_{p/3}(\ell) = d \log(\varepsilon_{\ell}^{p/3})/d \log \ell$. The exponents remain approximately constant across all scales for $p \le 5$, and vary by at most 10% within the inertial range for $p = 6$ and 7. The shaded bands represent the exponents $\zeta_p^{\parallel}$ defined in § 4.2.1 obtained from mixed structure functions.

Further analysis indicates that the dissipation field is well approximated by log-normal statistics. Figure 16(c) shows the probability density function of $\log \varepsilon_{\ell}$, standardised to zero mean and unit variance, for scales $0.2\ell_v \le \ell \le 10\ell_v$. While statistics for small dissipation (left branches) barely superpose, those for strong dissipative events largely collapse on a master curve, with only slight deviations from a Gaussian. Figure 16(b) shows that the dependency of the mean and variance of $\log \varepsilon_{\ell}$ with the coarse-graining scale $\ell$, revealing a scaling behaviour reasonably captured by

$$ \langle\log(\varepsilon_{\ell}/\varepsilon)\rangle = \mu_{62} \log(\ell/\ell_v), \quad \text{Var}[\log(\varepsilon_{\ell}/\varepsilon)] = -2\mu_{62} \log(\ell/\ell_v), \quad (4.7) $$

with $\mu_{62} = 0.16 \pm 0.02$. This phenomenology is akin to 3-D homogeneous isotropic turbulence – although the specific value of the constant $\mu_{62}$ is slightly larger than the value $\mu_{62}^{3D} \simeq 0.13$ commonly measured (Iyer et al. 2015; Chevillard et al. 2019). They suggest that the scaling exponents of the dissipation moments behave as $\chi_q \approx \mu_{62} q (1-q)$.

### 4.2.4. Nonlinarity of the scaling exponents

The three types of scaling exponents (mixed, ESS, refined) prove to be remarkably consistent between them. Figure 17 shows their behaviour as function of the order $p$. For ESS, where the value $C_p^\phi$ represents a ratio of exponents, we report the value of $\xi_2^{\parallel} C_p^\phi$, with
<figure><img src="image_18.png" /><figcaption>Figure 16. (a): similarity exponents for integer values of $p \in [1, 7]$. The coloured bands show the values for the mixed exponents. (b): the PDF of the logarithm of the coarse-grained dissipation field standardised to unit variance and zero mean. The shaded area indicates the inertial range identified in figure 8. (c): mean and variance of $\log(\varepsilon_{\ell}/\varepsilon)$ as a function of the coarse-graining scale.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17. Scaling exponents extracted from different methods: ESS ($\delta u_{\parallel}$, $\delta u_{\perp}$), mixed structure functions (●) and coarse-grained dissipation (■). Error bars indicate the root-mean-squared errors. The dashed black line shows K41 scaling and the red shaded area represents the log-normal models $\zeta_p = p/3 + \mu_{62}(p/3)(1 - p/3)$, with $\mu_{62} \in [0.14, 0.18]$.</figcaption></figure>
$\xi_2^\parallel$ obtained the mixed exponent of order 2. For all measurements, error bars represent the standard deviation of the local exponent within the inertial range.

In spite of the lack of usual scaling range, figure 17 suggests that SQG intermittency is characterised by nonlinear scaling exponents deviating from K41 similarity. The log-normal modelling of the dissipation field provides a reasonable fit up to orders $p \sim 6-7$ with $\mu_{62} \sim 0.16 \pm 0.02$. We observe no sign of saturation for the $\xi_p$'s at least up to the order $p=7$ that we computed. We, however, point out that in Navier-Stokes simulations (two- or three-dimensional) (Iyer, Sreenivasan & Yeung 2020; Müller & Krstulovic 2025), and passive scalar studies (Celani et al. 2000, 2001; Iyer et al. 2018), asymptotic behaviours of the exponents become apparent only at higher orders $p \gtrsim 10$.

Finally, the correspondence between dissipation and mixed structure function exponents suggests that, instead of (4.6), a consistent refined K62 phenomenology in SQG is mediated through

$$ \langle \delta u_{\parallel} | \delta \theta |^{p-1} \rangle \sim \ell^{p/3} \langle \varepsilon_{\ell}^{p/3} \rangle, \quad (4.8) $$

entailing the identification $\xi_p^\parallel = \xi_p = \chi_{p/3}(\ell) + p/3$ between the mixed and refined scaling exponents.

Note that (4.8) features the fluxes associated with the Casimirs $\langle \theta^{p-1} \rangle$ on its left-hand side. This, maybe, is counterintuitive. A local version of Yaglom's law in the spirit of Duchon & Robert (2000) might suggest the identification $\delta u(\delta \theta)^2 = -2\varepsilon_\ell \ell$ entailing rather $S_p^{WG} \sim \langle \varepsilon_\ell^{p/3} \rangle \ell^{p/3}$ instead of (4.8). This naive guess is, however, not correct: § 4.2.1 has ruled out scaling behaviours for the Watanabe–Gotoh structure functions.

## 5. Concluding remarks

In this paper, we investigated forced-dissipated SQG turbulence in a regime dominated by the direct cascade of scalar variance. The aim was to refine the analogy between SQG and Navier-Stokes (NS) turbulence, grounded in the Kolmogorov similarity prediction that $|\delta u| \sim |\delta \theta| \sim \ell^{1/3}$. Using systematic numerical simulations with up to 16 384² collocation points, we confirm that the two systems share many key turbulent features, including anomalous dissipation, exact laws for third-order structure functions, a skewness phenomenon and intermittency in the form of anomalous scaling. This substantiates the early observations of Sukhatme & Pierrehumbert (2002) in an unforced setting, which notably showed that SQG flows are strongly intermittent and, in that sense, closer to the 3-D direct cascade than to the 2-D inverse cascade.

However, the idea that $\theta$ is a scalar analogue of a 3-D turbulent velocity needs to be promoted only with great care. The analogy holds only if one considers mixed statistics involving correlations between the velocity and the scalar field. For statistics related to single fields $\theta$, $u_\parallel$ or $u_\perp$, quantifying intermittency requires examining self-scaling exponents that reflect relative scaling between structure functions with different orders. Moreover, the SQG system exhibits a very specific form of the skewness phenomenon. Such specificities do not fit into the usual NS turbulence frameworks, and this may be precisely what makes SQG interesting. As a conclusion, let us highlight three main SQG lessons, along with some speculations on their implications for NS turbulence.

(i) In SQG, there is a clear distinction between the advecting field $u$ and the advected scalar $\theta$, although both are dimensionally similar (and share a common K41 description). The fact that scaling appears only in certain mixed statistics highlights the fundamental role of transport and geometry in shaping turbulent
intermittency – more radically than in NS flows. Such effects may, however, also be present in NS turbulence and could help explain why transverse and longitudinal structure functions have different scaling properties (Iyer *et al.* 2018).

(ii) Another consequence of the duality between $\theta$ and $u$ is the absence of a straightforward relation between the skewness of the fields and cascade properties. This contrasts with NS turbulence, where the skewness of longitudinal increments reflects the cascade direction (negative for the 2-D and 3-D direct cascades, positive for the 2-D inverse cascade). In SQG, this relation breaks down; skewness arises from statistical symmetries and non-trivial kinematic constraints, making the skewness – cascade connection a genuine specificity of NS turbulence.

(iii) The universality of SQG turbulence remains unsettled. On the one hand, the ESS and scaling behaviours emphasised in this study are primarily tied to transport properties, which formally require only scale separation (i.e. the existence of an inertial range). From this perspective, we expect the self-scaling properties of SQG fields to be relatively insensitive to the precise details of forcing. Supporting this view, we note that hallmark transport features of SQG turbulence, such as the approximate $k^{-5/3}$ energy spectrum and the validity of Yaglom's law, are also observed in numerical studies employing comparable (though not identical) forcing protocols (Valadão *et al.* 2024).

On the other hand, due to large-scale damping in our forcing protocol, only a moderate fraction of the injected energy is transferred towards the small scales. It is plausible that different statistical ‘phases’ of SQG turbulence emerge upon varying the degree of energy leakage from the injection to the dissipation scales. In particular, let us keep in mind that scale-to-scale transfers in SQG are non-local (Watanabe & Iwayama 2007), and one may therefore anticipate non-trivial imprints of the large scales on the forward cascade. This may partly explain the lack of classical power-law scaling observed in single-field statistics and the apparent weak decay of dissipation with increasing *Pe*.

While speculative, the latter trend raises broader questions about the universality of turbulent anomalous dissipation across systems. Recent studies have shown that weakly vanishing dissipation remains theoretically consistent with 3-D NS turbulence (Drivas & Eyink 2019; Iyer 2023), and a similar behaviour has been reported in the highest-to-date resolved DNS of 3-D NS flows (Iyer *et al.* 2025). Those studies rely on deterministic forcing protocols that couple large-scale energy injection to dissipation. One can speculate that non-local couplings between scales, which are present in SQG but *a priori* absent in the classical theory of 3-D turbulence, may lead to a weak form of turbulence non-universality.

**Acknowledgements.** The authors acknowledge stimulating discussions with G. Boffetta, C. Campolina, S. Musacchio, É. Simonnet and V. Valadão. We also thank A. Mailybaev and the Instituto de Matemática Pura e Aplicada, Rio de Janeiro, for support and hospitality (through summer visiting grants), where a part of this work was undertaken. The authors are grateful to the Université Côte d'Azur's Center for High-Performance Computing (OPAL infrastructure) for providing resources and support.

**Funding.** This work was supported by the French National Research Agency (ANR project TILT; ANR-20-CE30-0035) and the French government through the France 2030 investment plan managed by the ANR, as part of the Initiative of Excellence Université Côte d'Azur under reference number ANR-15-IDEX-01.

**Declaration of interests.** The authors report no conflict of interest.

**Data availability statement.** The data that support the findings of this study are available upon reasonable request.
<figure><img src="image_20.png" /><figcaption>Figure 18. Filaments (a) and vortices (b) as two prototypes of inviscid SQG flows.</figcaption></figure>

# Appendix A. Standard features of SQG

## A.1. Prototype SQG flows

It is useful to have in mind prototype inviscid SQG flows, in order to depict SQG turbulence as a mixture of filaments and vortices. We refer the reader to Juckes (1995) for the relevant stability analysis and to figure 18 for a sketch of those specific SQG flows in the inviscid case.

### A.1.1. Filaments

The SQG filaments of width $2\delta$, say aligned along $x_1$, are prescribed by the scalar field

$$ \theta(\mathbf{x}) = \begin{cases} \frac{U_0}{2} & |x_2| \le \delta, \\ 0 & \text{otherwise,} \end{cases} \qquad (A1) $$

which in turn generate a horizontal velocity $u(\mathbf{x}) = U_0 \log |(x_2 - \delta)/(x_2 + \delta)| \hat{\mathbf{1}}$. The SQG filaments prescribe two counter-propagating jets, diverging logarithmically at the boundaries of the filament. Far from the filament, the horizontal velocity decreases as $u_1(\mathbf{x}) \sim -2U_0\delta x_2^{-1}$.

### A.1.2. Rankine patches

In analogy to the 2-D Rankine vortex, Rankine SQG patches of radius $r_0$ are prescribed by the scalar field

$$ \theta(\mathbf{x}) = \begin{cases} \theta_0, & |\mathbf{x}| \le r_0, \\ 0 & \text{otherwise,} \end{cases} \qquad (A2) $$

and prescribe the circular flow

$$ \mathbf{u}(\mathbf{x}) = u_{\perp} \hat{\mathbf{x}}^{\perp}, \quad u_{\perp} := \theta_0 \int_{0}^{+\infty} \mathcal{J}_1 \left( v \frac{\mathbf{x}}{r_0} \right) \mathcal{J}_1 (v) \, dv, \qquad (A3) $$
where $J_0$ and $J_1$ are Bessel functions of the first kind, with orders 0 and 1. At the boundary of the patch, the flow has a logarithmic divergence and far from the patch, the velocity decays algebraically, namely

$$u_{\perp}(\mathbf{x}) \propto \log|r_0 - |\mathbf{x}||, \quad u_{\perp}(\mathbf{x}) \stackrel{+\infty}{\propto} |\mathbf{x}|^{-2}. \qquad (A4)$$

## A.2. Symmetries

### A.2.1. Classical symmetries

Similar to NS equations, the SQG system (1.1) considered in $\mathbb{R}^2$ is formally invariant under a certain number of continuous symmetries. Following the formulations of Frisch (1995), we list below the symmetries holding in particular in the inviscid case:

(i) time or space translation: $t \mapsto t + \tau$, $\mathbf{x} \mapsto \mathbf{x} + \mathbf{r}$ with $\tau \in \mathbb{R}$ and $\mathbf{r} \in \mathbb{R}^2$;

(ii) rotation and reflection: $\mathbf{x} \mapsto \mathcal{O}\mathbf{x}$, $\mathbf{u} \mapsto \mathcal{O}\mathbf{u}$, $\theta \mapsto \det(\mathcal{O})\theta$, with $\mathcal{O} \in \mathcal{O}(2)$. Pure rotations correspond to $\mathcal{O} \in \mathcal{O}(2)$, characterised by $\det \mathcal{O} = 1$. Reflections correspond to $\mathcal{O} \in \mathcal{O}(2) \setminus \mathcal{O}(2)$, with $\det \mathcal{O} = -1$;

(iii) Time or space rescaling: $t \mapsto \alpha t$, $\mathbf{x} \mapsto \lambda \mathbf{x}$, $\mathbf{u} \mapsto (\lambda/\alpha)\mathbf{u}$, $\theta \mapsto (\lambda/\alpha)\theta$ with $\alpha, \lambda > 0$;

(iv) Galilean invariance: $\theta \mapsto \theta + \theta_0$, $\mathbf{u} \mapsto \mathbf{u} + \mathbf{u}_0$, $\mathbf{x} \mapsto \mathbf{x} + t\mathbf{u}_0$ with $\theta_0 \in \mathbb{R}$, $\mathbf{u}_0 \in \mathbb{R}^2$. Please note that, to be properly formulated, the Galilean invariance needs to be associated with suitable transformed Green's functions (Jackson 2021), ensuring the non-vanishing asymptotics $\mathbf{u} \sim \mathbf{u}_0$, $\psi \sim -\mathbf{x}^{\perp} \cdot \mathbf{u}_0$ as $|\mathbf{x}| \to \infty$. This proper formulation, however, goes beyond the scope of this work.

Taking for example space-time rescaling, the formulation specifically means that if the pair $\mathbf{u}(\mathbf{x}, t), \theta(\mathbf{x}, t)$ solves the inviscid SQG system, then so does the pair

$$\mathbf{u}'(\mathbf{x}', t') := (\lambda/\alpha)\mathbf{u}(\mathbf{x}'/\lambda, t'/\alpha), \quad \theta'(\mathbf{x}', t') := (\lambda/\alpha)\theta'(\mathbf{x}'/\lambda, t'/\alpha). \qquad (A5)$$

### A.2.2. Statistical symmetries

The statistical symmetries are statistical formulation of the classical symmetries, which hold at the level of the turbulent averages $\langle \cdot \rangle$. They mean

$$\mathbf{u}'(\mathbf{x}', t'), \theta'(\mathbf{x}', t') \stackrel{\text{law}}{\sim} \mathbf{u}(\mathbf{x}, t), \theta(\mathbf{x}, t), \qquad (A6)$$

where the transforms $\mathbf{x}, t, \mathbf{u}, \theta \mapsto \mathbf{x}', t', \mathbf{u}', \theta'$ are specified by the classical symmetries listed in § Appendix A.2.1, and respectively yield stationarity or homogeneity, isotropy or mirror symmetry and scale invariance. We let aside the Galilean invariance, which is not discussed in the present paper.

By construction, the statistical symmetries have consequences for the distributions of the SQG fields. We list below the ones with relevance for the present work.

(i) One-point statistics

$$p_{\theta}(y) = p_{\theta}(-y), \quad p_{u_i}(y) = p_{u_i}(-y), \qquad (A7)$$

where $p_{\theta}(y) = \langle \delta(\theta(\mathbf{x}_0) - y) \rangle$ is the one-point PDF of $\theta$, and $p_{u_1}(y), p_{u_2}(y_0)$ that of each velocity component, with according definitions. For $\theta$, the symmetry stems from the identities $\langle e^{i\theta(\mathbf{x}_0)} \rangle = \langle e^{i\theta(0)} \rangle = \langle e^{-i\theta(0)} \rangle$, employing successively homogeneity and mirror symmetry, and implying that all odd moments vanish. For the velocity, the same computation holds true, from isotropy instead of mirror symmetry.
(ii) Increment statistics

$$p_{\delta\theta}(y) = p_{\delta\theta}(-y), \quad p_{\delta u_{\perp}}(y) = p_{\delta u_{\perp}}(-y), \qquad (A8)$$

implying that the odd order moments vanish: $\langle \delta\theta^{2p+1} \rangle = \langle \delta u_{\perp}^{2p+1} \rangle = 0$.

(iii) Mixed structure functions

$$\forall p \in \mathbb{N} \quad \langle \delta u_{\parallel} \delta \theta^{2p+1} \rangle = \langle \delta u_{\parallel} \delta u_{\perp}^{2p+1} \rangle = 0, \quad \langle \mathbf{u}(\mathbf{x}) \theta^2 \rangle = 0, \quad \langle \mathbf{u}(\mathbf{x}) \omega^2 \rangle = 0. \qquad (A9)$$

The second identity stems from the observations $\nabla \cdot \langle \mathbf{u}(\mathbf{x})\theta(0)^2 \rangle = 0$ (incompressibility) and $\nabla^{\perp} \cdot \langle \mathbf{u}(\mathbf{x})\theta^2 \rangle = \langle \omega(\mathbf{x})\theta^2 \rangle = 0$ (mirror symmetry). Similar arguments apply for the third identity.

## Appendix B. Derivations of Bernard's balances

This section summarises the technical derivation leading to the SQG versions of Bernard's balances of § 3.4, either in their complex or differential formulations.

### B.1. The complex formulation (3.4)

Following standard conventions in complex analysis (see, e.g. Tao 2016), we introduce the (Wirtinger) complex derivatives as $\partial_z = (1/2)(\partial_1 - i\partial_2)$ and $\partial_{z^*} = (1/2)(\partial_1 + i\partial_2)$. Combining the velocity components into a single complex velocity field $U(z, z^*) = u_1 + iu_2$, and promoting the vorticity and the streamfunction into the complex fields $\Omega(z, z^*)$, $\Psi(z, z^*)$, one checks the relations

$$U(z, z^*) = -2i \partial_{z^*} \Psi, \quad \Omega = -2i \partial_z U, \quad \Omega = -4 \partial_{zz^*} \Psi, \qquad (B1)$$

where the second identity uses incompressibility. We now define the structure functions

$$S_3^U(\xi, \xi^*) = \langle \delta U^3 \rangle, \quad S_3^\Theta(\xi, \xi^*) = \langle (\delta\Theta)^2 \delta U \rangle, \quad S_3^\Omega(\xi, \xi^*) = \langle (\delta\Omega)^2 \delta U \rangle, \qquad (B2)$$

where $\delta F = F(z + \xi, z^* + \xi^*) - F(z, z^*)$ denotes complex increments of the field $F$.

With the shorthands $U = U(z, z^*)$, $U' = U(z + \xi, z^* + \xi^*)$, etc., routine calculations employing statistical homogeneity, isotropy and mirror symmetry yield

$$S_3^U = 6 \langle U^2 U' \rangle, \quad S_3^\Theta = 4 \langle U \Theta \Theta' \rangle, \quad S_3^\Omega = 4 \langle U \Omega \Omega' \rangle. \qquad (B3)$$

Let us simply emphasise the crucial role played by the mirror symmetry, which makes the correlators $\langle U'\Omega^2 \rangle$, $\langle U'\Theta^2 \rangle$ vanish! From (B2) and homogeneity, the following calculation string holds:

$$\partial_{\xi\xi} \langle (\delta U)^3 \rangle = -6 \langle (\partial'_z U') \partial_z (U^2) \rangle = 3 \langle U \Omega \Omega' \rangle. \qquad (B4)$$

Using (B3), we then recover the kinematic relation derived in Bernard (1999), namely

$$\partial_{\xi\xi} S_3^U = \frac{3}{4} S_3^\Omega. \qquad (B5)$$

Being of purely kinematic nature, the relation (B5) holds true both in SQG and in 2-D NS. In order to obtain the SQG balance (3.4), we now introduce the ageostrophic field and his complex counterpart as, respectively,

$$\mathbf{a} = (-\Delta)^{1/2} (\mathbf{u}\theta) - \mathbf{u}\omega, \quad A := a_1 + i a_2 = 2 |\partial_z| [U\Theta] - U\Omega, \qquad (B6)$$

where the shorthand notation $|\partial_z|$ stands for the complex counterpart of $(1/2)(-\Delta)^{1/2}$, such that $\Omega = 2|\partial_z|\Theta$ and $|\partial_z|^2 = -\partial_{zz^*}$.
We then observe

$$ \partial_{\xi\xi^*} \langle U\Theta\Theta' \rangle = \langle |\partial_z| [U\Theta] |\partial_{z'}| [\Theta'] \rangle = \frac{1}{4} \langle (A + U\Omega) \Omega' \rangle, \quad (B7) $$

pointing out that the lack of a minus sign in the first equality comes from homogeneity. Combining (B7) and (B3) leads to

$$ \partial_{\xi\xi^*} S_3^{\Theta} = \langle A\Omega' \rangle + \frac{1}{4} S_3^{\Omega}. \quad (B8) $$

Combining (B5) and (B8) finally yields

$$ \partial_{\xi\xi^*} S_3^{\Theta} = \langle A\Omega' \rangle + \frac{1}{3} \partial_{\xi\xi} S_3^U, \quad (B9) $$

which is the desired formula (3.4).

## B.2. The differential formulation (3.5)

Deriving (3.5) from (3.4) relies on algebraic manipulations and the theory of isotropic tensors.

### B.2.1. Isotropic tensors in dimension d

Let us first recall general relations on isotropic tensors. We follow Brachet (2000), to observe that the third-order-but-one-separation isotropic tensors

$$ B_{ijk}(\ell) = \langle \delta u_i \delta u_j \delta u_k \rangle, \quad \delta u_i := u_i(\ell) - u_i(0), \quad (B10) $$

can be written as

$$ B_{ijk}(\ell) = 2 (b_{ij,k} + b_{ki,j} + b_{jk,i}), \quad b_{ij,k} := \langle u_i(0) u_j(0) u_k(\ell) \rangle, \quad (B11) $$

where homogeneity obviously implies the identity $b_{ij,k} = -\langle u_i(\ell)u_j(\ell)u_k(0) \rangle$. As an isotropic tensor, $b_{ij,k}$ takes the form

$$ b_{ij,k} = C(\ell) \frac{\delta_{ij} \ell_k}{\ell} + D(\ell) \frac{\delta_{ik} \ell_j + \delta_{jk} \ell_i}{\ell} + F(\ell) \frac{\ell_i \ell_j \ell_k}{\ell^3}, \quad (B12) $$

which with the incompressibility condition $\partial_k b_{ij,k} = 0$ imposes

$$ -2D = (d-1)C + rC' \quad \& \quad F = -C + \ell \partial_{\ell} C, \quad (B13) $$

with $d$ denoting space dimensionality ($d=2$ for SQG); this leaves the function $C$ as the only unknown. Direct calculation led to the following useful properties:

$$ \langle \delta u_{\parallel}^3 \rangle = B_{ijk} \hat{\ell}_i \hat{\ell}_j \hat{\ell}_k = 6(1-d)C, \quad (B14a) $$

$$ \langle \delta u_{\parallel} | \delta \mathbf{u} |^2 \rangle = \hat{\ell}_i B_{ijj} = 2(1-d) \left( (1+d)C + \ell \partial_{\ell} C \right), \quad (B14b) $$

$$ \langle \delta u_{\perp} | \delta \mathbf{u} |^2 \rangle = \hat{n}_i B_{ijj} = 0. \quad (B14c) $$

### B.2.2. The SQG

We now focus on the case $d=2$, relevant for SQG. We first observe that the complex structure functions $S^U, S^\Theta$ defined in (B2) relate to their real counterpart as

$$ \ell S^U = (\langle \|\delta u\|^2 \delta u_{\parallel} \rangle + i \langle \|\delta u\|^2 \delta u_{\perp} \rangle) \xi, \quad \ell S^{\Theta} = (\langle |\delta \theta|^2 \delta u_{\parallel} \rangle + i \langle |\delta \theta|^2 \delta u_{\perp} \rangle) \xi. \quad (B15) $$
In line with (B14a), isotropy and mirror symmetry make the perpendicular contributions vanish, such that the previous equation in fact reduces to

$$ \ell S^U = (|\delta \mathbf{u}|^2 \delta u_\parallel)\xi, \quad \ell S^\Theta = (|\delta \theta|^2 \delta u_\parallel)\xi. \qquad (B16) $$

On the one hand, now combining (B13) and (B14a), we obtain

$$ C(\ell) = -\frac{1}{6} \langle \delta u_{\parallel}^3 \rangle, \quad S_3^U = 6 \frac{F}{\ell^3} \xi^3, \quad F := \ell^2 \partial_{\ell} (\ell^{-1} C). \qquad (B17) $$

Using $\partial_\xi f(\ell) = \xi^* \partial_\ell f/(2\ell)$, this in turn yields the first piece of the differential formulation

$$ \partial_{\xi\xi} S_3^U = \frac{3\xi}{2\ell^3} \partial_\ell (\ell^{-1} \partial_\ell \ell^3 F) = 3 \frac{\xi}{\ell} \mathcal{L}_{u_\parallel} [\langle \delta u_\parallel^3 \rangle], \\ \text{identifying } \mathcal{L}_{u_\parallel} := -\frac{1}{12} \ell^{-2} \partial_\ell (\ell^{-1} \partial_\ell (\ell^5 \partial_\ell (\ell^{-1} \cdot))). \qquad (B18) $$

On the other hand, we also compute from (B16), using $\partial_{\xi^*} f(\ell) = \xi \partial_\ell f/(2\ell)$

$$ \partial_{\xi\xi^*} S_3^\Theta = \frac{\xi}{\ell} \mathcal{L}_M [ |\delta\theta|^2 \delta u_\parallel ], \quad \mathcal{L}_M := \frac{1}{4} \partial_\ell (\ell^{-1} \partial_\ell (\ell \cdot)). \qquad (B19) $$

Similar to (B15), we also have

$$ \langle A\Omega' \rangle = -\langle A'\Omega \rangle = -\frac{\xi}{\ell} \left( \langle a_\parallel(\ell\mathbf{e}_1)\omega(0) \rangle + i \langle a_\perp(\ell\mathbf{e}_1)\omega(0) \rangle \right). \qquad (B20) $$

Using (B9) with (B18), (B19) and (B20) yields $\langle a_\perp(\ell\mathbf{e}_1)\omega(0) \rangle = 0$ and the differential formulation (3.5).

## REFERENCES

* BENZI, R., BIFERALE, L., CILIBERTO, S., STRUGLIA, M. & TRIPICCIONE, R. 1996 Generalized scaling in fully developed turbulence. *Physica D: Nonlinear Phenom.* **96** (1-4), 162–181.

* BENZI, R., CILIBERTO, S., BAUDET, C., CHAVARRIA, G.R. & TRIPICCIONE, R. 1993a Extended self-similarity in the dissipation range of fully developed turbulence. *Europhys. Lett.* **24** (4), 275.

* BENZI, R., CILIBERTO, S., TRIPICCIONE, R., BAUDET, C., MASSAIOLI, F. & SUCCI, S. 1993b Extended self-similarity in turbulent flows. *Phys. Rev. E* **48** (1), R29.

* BENZI, R. & TOSCHI, F. 2023 Lectures on turbulence. *Phys. Rep.* **1021**, 1–106.

* BERNARD, D. 1999 Three-point velocity correlation functions in two-dimensional forced turbulence. *Phys. Rev. E* **60** (5), 6184.

* BERNARD, D., BOFFETTA, G., CELANI, A. & FALKOVICH, G. 2007 Inverse turbulent cascades and conformally invariant curves. *Phys. Rev. Lett.* **98** (2), 024501.

* BLUMEN, W. 1978 Uniform potential vorticity flow: part I. Theory of wave interactions and two-dimensional turbulence. *J. Atmos. Sci.* **35**, 774–783.

* BOFFETTA, G. & ECKE, R.E. 2012 Two-dimensional turbulence. *Annu. Rev. Fluid Mech.* **44**, 427–451.

* BRACHET, M. 2000 A primer in classical turbulence theory. In *Instabilities and Nonequilibrium Structures VI*, pp. 5–34. Springer.

* BUARIA, D., PUMIR, A., BODENSCHATZ, E. & YEUNG, P.-K. 2019 Extreme velocity gradients in turbulent flows. *New J. Phys.* **21** (4), 043004.

* BURGESS, B., SCOTT, R. & SHEPHERD, T. 2015 Kraichnan–Leith–Batchelor similarity theory and two-dimensional inverse cascades. *J. Fluid Mech.* **767**, 467–496.

* CAPET, X., KLEIN, P., HUA, B.L., LAPEYRE, G. & MCWILLIAMS, J. 2008 Surface kinetic energy transfer in surface quasi-geostrophic flows. *J. Fluid Mech.* **604**, 165–174.

* CELANI, A., CENCINI, M., MAZZINO, A. & VERGASSOLA, M. 2004 Active and passive fields face to face. *New J. Phys.* **6**, 72.

* CELANI, A., LANOTTE, A., MAZZINO, A. & VERGASSOLA, M. 2000 Universality and saturation of intermittency in passive scalar turbulence. *Phys. Rev. Lett.* **84** (11), 2385.
CELANI, A., LANOTTE, A., MAZZINO, A. & VERGASSOLA, M. 2001 Fronts in passive scalar turbulence. *Phys. Fluids* **13** (6), 1768–1783.

CHEVILLARD, L., GARBAN, C., RHODES, R. & VARGAS, V. 2019 On a skewed and multifractal unidimensional random field, as a probabilistic representation of Kolmogorov’s views on turbulence. *Ann. Henri Poincare* **20**, 3693–3741.

CONSTANTIN, P., MAJDA, A.J. & TABAK, E. 1994 Formation of strong fronts in the 2-D quasigeostrophic thermal active scalar. *Nonlinearity* **7**, 1495.

DRIVAS, T. & EYINK, G. 2019 An Onsager singularity theorem for Leray solutions of incompressible Navier–Stokes. *Nonlinearity* **32** (11), 4465.

DUBRULLE, B. 2019 Beyond Kolmogorov cascades. *J. Fluid Mech.* **867**, P1.

DUCHON, J. & ROBERT, R. 2000 Inertial energy dissipation for weak solutions of incompressible Euler and Navier–Stokes equations. *Nonlinearity* **13** (1), 249.

FOUSSARD, A., BERTI, S., PERROT, X. & LAPEYRE, G. 2017 Relative dispersion in generalized two-dimensional turbulence. *J. Fluid Mech.* **821**, 358–383.

FRISCH, U. 1995 *Turbulence: the Legacy of AN Kolmogorov*. Cambridge University Press.

GAUDING, M., DANAILA, L. & VAREA, E. 2017 High-order structure functions for passive scalar fed by a mean gradient. *Intl J. Heat Fluid Flow* **67**, 86–93.

GRAUER, R., HOMANN, H. & PINTON, J.-F. 2012 Longitudinal and transverse structure functions in high-Reynolds-number turbulence. *New J. Phys.* **14** (6), 063016.

HAKIM, G.J., SNYDER, C. & MURAKI, D.J. 2002 A new surface model for cyclone-anticyclone asymmetry. *J. Atmos. Sci.* **59** (16), 2405–2420.

HELD, J., PIERREHUMBERT, R., GARNER, S. & SWANSON, K. 1995 Surface quasi-geostrophic dynamics. *J. Fluid Mech.* **282**, 1–20.

HOMANN, H., BEC, J. & GRAUER, R. 2013 Effect of turbulent fluctuations on the drag and lift forces on a towed sphere and its boundary layer. *J. Fluid Mech.* **721**, 155–179.

IYER, K. 2023 Dependence of the asymptotic energy dissipation on third-order velocity scaling. *Phys. Rev. Fluids* **8** (8), L082601.

IYER, K., DRIVAS, T., EYINK, G. & SREENIVASAN, K. 2025, Whither the zeroth law of turbulence? arXiv preprint arXiv: 2504.13298.

IYER, K., SCHUMACHER, J., SREENIVASAN, K. & YEUNG, P.-K. 2018 Steep cliffs and saturated exponents in three-dimensional scalar turbulence. *Phys. Rev. Lett.* **121** (26), 264501.

IYER, K., SREENIVASAN, K. & YEUNG, P.-K. 2015 Refined similarity hypothesis using three-dimensional local averages. *Phys. Rev. E* **92** (6), 063024.

IYER, K., SREENIVASAN, K. & YEUNG, P.-K. 2020 Scaling exponents saturate in three-dimensional isotropic turbulence. *Phys. Rev. Fluids* **5** (5), 054605.

JACKSON, J. 2021 *Classical Electrodynamics*. Wiley & Sons.

JHA, V.B., SESHASAYANAN, K. & DALLAS, V. 2025 Cascades transition in generalised two-dimensional turbulence. *J. Fluid Mech.* **1008**, A23.

JUCKES, M. 1995 Instability of surface and upper-tropospheric shear lines. *J. Atmos. Sci.* **52**, 3247–3262.

KOLMOGOROV, A.N. 1962 A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. *J. Fluid Mech.* **13** (1), 82–85.

LACASCE, J. & MAHADEVAN, A. 2006 Estimating subsurface horizontal and vertical velocities from sea-surface temperature. *J. Mar. Res.* **64**, 695–721.

LAPEYRE, G. 2017 Surface quasi-geostrophy. *Fluids* **2**, 7.

LAPEYRE, G. & KLEIN, P. 2006 Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory. *J. Phys. Oceanogr.* **36**, 165–176.

MÜLLER, N.P. & KRSTULOVIC, G. 2025 Lack of self-similarity in transverse velocity increments and circulation statistics in two-dimensional turbulence. *Phys. Rev. Fluids* **10** (1), L012601.

ORSZAG, S. 1971 On the elimination of aliasing in finite-difference schemes by filtering high-wavenumber components. *J. Atmos. Sci.* **28** (6), 1074–1074.

PEDLOSKY, J. 2013 *Geophysical Fluid Dynamics*. Springer Science & Business Media.

SCOTT, R.K. & DRITSCHEL, D.G. 2014 Numerical simulation of a self-similar cascade of filament instabilities in the surface quasigeostrophic system. *Phys. Rev. Lett.* **112** (14), 144505.

SMITH, K., BOCCALETTI, G., HENNING, C., MARINOV, I., TAM, C., HELD, I. & VALLIS, G. 2002 Turbulent diffusion in the geostrophic inverse cascade. *J. Fluid Mech.* **469**, 13–48.

SUKHATME, J. & PIERREHUMBERT, R. 2002 Surface quasigeostrophic turbulence: the study of an active scalar. *Chaos* **12**, 439–450.

TAO, T. 2016 246a, notes 1: complex differentiation.
VALADÃO, V.J., CECCOTTI, T., BOFFETTA, G. & MUSACCHIO, S. 2024 Nonequilibrium fluctuations of the direct cascade in surface quasi-geostrophic turbulence. *Phys. Rev. Fluids* **9**, 094601.

VALADÃO, V.J., DE LILLO, F., MUSACCHIO, S. & BOFFETTA, G. 2025 Scaling and predictability in surface quasi-geostrophic turbulence. *J. Fluid Mech.* **1017**, A38.

VALADE, N., THALABARD, S. & BEC, J. 2024 Anomalous dissipation and spontaneous stochasticity in deterministic surface quasi-geostrophic flow. *Ann. Henri Poincaré* **25**, 1–23.

VALLIS, G.K. 2017 *Atmospheric and Oceanic Fluid Dynamics*. Cambridge University Press.

WATANABE, T. & GOTOH, T. 2004 Statistics of a passive scalar in homogeneous turbulence. *New J. Phys.* **6** (1), 40.

WATANABE, T. & IWAYAMA, T. 2007 Interacting scales and triad enstrophy transfers in generalized two-dimensional turbulence. *Phys. Rev. E* **76**, 046303.

YAGLOM, A. 1949 On the local structure of a temperature field in a turbulent flow. *Dokl. Akad. Nauk SSSR* **69**, 743–746.