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
\begin{align} & \partial _t \theta + \left ({\boldsymbol{u}} \boldsymbol{\cdot } \boldsymbol{\nabla }\right )\theta = \nu \Delta \theta + F, \quad {\boldsymbol{u}} := -\boldsymbol{\nabla }^\perp \psi , \quad \psi :=(-\varDelta )^{-\frac {1}{2}}\theta , \end{align}
where
$\theta ({\boldsymbol{x}},t)$
represents a scalar field (such as surface temperature or buoyancy) advected in a planar domain
$\varOmega \subset \mathbb{R}^2$
under the influence of forcing
$F(\theta , {\boldsymbol{x}},t)$
and diffusivity
$\nu$
. Here,
$\boldsymbol{\nabla }^\perp = (-\partial _{2},\partial _{1})^\top$
tying the incompressible velocity field
$\boldsymbol{u}$
to the stream function
$\psi$
and
$\varDelta$
respectively denote the perpendicular gradient and the Laplacian operators.
Equation (1.1) describes the transport of an active scalar, where the velocity field
$\boldsymbol{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
$(-\varDelta)^{-1/2}$
associated with the half-Laplacian rather than to the Laplacian. In the full plane, this yields
\begin{equation} \psi ({\boldsymbol{x}}) =\frac {1}{2\pi } \int _{\mathbb R^2} \frac {\theta ({\boldsymbol{x}}')}{|{\boldsymbol{x}}-{\boldsymbol{x}}'|}\,\textrm{d}^2{\boldsymbol{x}}'\quad \text{ and } \quad {\boldsymbol{u}}({\boldsymbol{x}}) =-\frac {1}{2\pi } \int _{\mathbb R^2} \frac { ({\boldsymbol{x}}-{\boldsymbol{x}}')^\perp }{|{\boldsymbol{x}}-{\boldsymbol{x}}'|^3} \theta ({\boldsymbol{x}}')\,\textrm{d}^2{\boldsymbol{x}}', \end{equation}
where
${\boldsymbol{x}}^\perp = (-x_2,x_1)$
. These relations imply, somewhat surprisingly, that
$\theta$
and
$\boldsymbol{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 (Reference Blumen1978), 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. Reference Held, Pierrehumbert, Garner and Swanson1995; Juckes Reference Juckes1995; Hakim, Snyder & Muraki Reference Hakim, Snyder and Muraki2002; Lapeyre & Klein Reference Lapeyre and Klein2006; LaCasce & Mahadevan Reference LaCasce and Mahadevan2006 and the comprehensive review by Lapeyre Reference Lapeyre2017) . 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. Reference Held, Pierrehumbert, Garner and Swanson1995; Bernard et al. Reference Bernard, Boffetta, Celani and Falkovich2007), passive scalar turbulence (Celani et al. Reference Celani, Cencini, Mazzino and Vergassola2004), Burgers turbulence (Valade, Thalabard & Bec Reference Valade, Thalabard and Bec2024) and the 3-D Navier–Stokes equations (Constantin, Majda & Tabak Reference Constantin, Majda and Tabak1994; Sukhatme & Pierrehumbert Reference Sukhatme and Pierrehumbert2002).
In the absence of forcing and dissipation, smooth solutions of the SQG equation conserve two quadratic invariants: the Hamiltonian
$\mathcal{H}= { ({1}/{2})\int _\varOmega \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 _\varOmega \theta ^2=({1}/{2}) \int _\varOmega |{\boldsymbol{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. Reference Smith, Boccaletti, Henning, Marinov, Tam, Held and Vallis2002; Burgess, Scott & Shepherd Reference Burgess, Scott and Shepherd2015).
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 Reference Sukhatme and Pierrehumbert2002; Celani et al. Reference Celani, Cencini, Mazzino and Vergassola2004; Capet et al. Reference Capet, Klein, Hua, Lapeyre and Mcwilliams2008) and shallower ones (Watanabe & Iwayama Reference Watanabe and Iwayama2007). These observations point toward the presence of significant corrections due to intermittency (Kolmogorov Reference Kolmogorov1962; Frisch Reference Frisch1995; Benzi & Toschi Reference Benzi and Toschi2023). 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. (Reference Valadão, Ceccotti, Boffetta and Musacchio2024) provided evidence that large-scale fluctuations can obscure a
$k^{-5/3}$
scaling, with steeper subleading terms appearing in the spectrum.
Figure 1. SQG snapshots of squared vorticity
$\omega =\boldsymbol{\nabla }^\perp \boldsymbol{\cdot } {\boldsymbol{u}}$
, temperature
$\theta$
and velocity
$\boldsymbol{u}$
at a typical time. The zoomed-in regions display
$(1/32)^2$
of the computational domain. Taken from run VI (see table 1).
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^2$
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. Reference Held, Pierrehumbert, Garner and Swanson1995; Smith et al. Reference Smith, Boccaletti, Henning, Marinov, Tam, Held and Vallis2002; Celani et al. Reference Celani, Cencini, Mazzino and Vergassola2004; Watanabe & Iwayama Reference Watanabe and Iwayama2007; Capet et al. Reference Capet, Klein, Hua, Lapeyre and Mcwilliams2008; Foussard et al. Reference Foussard, Berti, Perrot and Lapeyre2017; Jha, Seshasayanan & Dallas Reference Jha, Seshasayanan and Dallas2025). 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 Reference Sukhatme and Pierrehumbert2002; Valade et al. Reference Valade, Thalabard and Bec2024; Valadão et al. Reference Valadão, Ceccotti, Boffetta and Musacchio2024; Valadão et al. Reference Valadão, De Lillo, Musacchio and Boffetta2025).
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:
\begin{equation} \lim _{\nu \to 0} \big\langle {\theta ^2}\big\rangle \lt \infty \quad \text{and } \quad \lim _{\nu \to 0} \big\langle {\nu |\boldsymbol{\nabla } \theta |^2 }\big\rangle \gt 0, \end{equation}
where
$\left \langle {\boldsymbol{\cdot }}\right \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 Reference Watanabe and Iwayama2007; Foussard et al. Reference Foussard, Berti, Perrot and Lapeyre2017), 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 (Reference Kolmogorov1962) refined similarity hypothesis, which links the non-Gaussianity of velocity increments to fluctuations in the dissipation field (Chevillard et al. Reference Chevillard, Garban, Rhodes and Vargas2019; Dubrulle Reference Dubrulle2019). While evidence of intermittency in SQG turbulence has been reported in unforced settings (Sukhatme & Pierrehumbert Reference Sukhatme and Pierrehumbert2002; Valade et al. Reference Valade, Thalabard and Bec2024), 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 Reference Scott and Dritschel2014; Lapeyre Reference Lapeyre2017). As a result, although the scalar field
$\theta$
and the velocity field
$\boldsymbol{u}$
share certain features – such as their variance, dissipation rate and dimensional scaling – they exhibit distinct statistical behaviours. For instance, Capet et al. (Reference Capet, Klein, Hua, Lapeyre and Mcwilliams2008) 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 Reference Yaglom1949). 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,
$\varOmega =\mathbb{T}^2$
, by decomposing the fields into their Fourier components. The Fourier transform of the scalar field is defined as
$\hat {\theta }_{\boldsymbol{k}} = ({1}/{(2\pi )^2}) \int _{\mathbb{T}^2} \theta ({\boldsymbol{x}})\,\textrm{e}^{-\textrm{i}\,{\boldsymbol{k}}\boldsymbol{\cdot }{\boldsymbol{x}}}\,\textrm{d}^2x$
, such that
$\theta ({\boldsymbol{x}}) = \sum _{{\boldsymbol{k}}\in \mathbb{Z}^2} \hat {\theta }_{\boldsymbol{k}} \textrm{e}^{\textrm{i}\,{\boldsymbol{k}}\boldsymbol{\cdot }{\boldsymbol{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
\begin{equation} \hat {{\boldsymbol{u}}}_{\boldsymbol{k}} = -\textrm{i}\frac {{\boldsymbol{k}}^\perp }{|{\boldsymbol{k}}|}\hat {\theta }_{\boldsymbol{k}} \quad \text{and} \quad \hat {\theta }_{\boldsymbol{k}} = \textrm{i}\frac {{\boldsymbol{k}}^\perp }{|{\boldsymbol{k}}|} \boldsymbol{\cdot } \hat {{\boldsymbol{u}}}_{\boldsymbol{k}} \quad \text{for } {\boldsymbol{k}}\in \mathbb{Z}^2\backslash \{\boldsymbol{0}\}, \end{equation}
where the null components are set to zero:
$\hat {{\boldsymbol{u}}}_{\boldsymbol{0}} = \boldsymbol{0}$
and
$\hat \theta _{\boldsymbol{0}} = 0$
. From (2.1), it follows that the Fourier components of
$\theta$
and
$\boldsymbol{u}$
have the same amplitudes:
$|\hat {{\boldsymbol{u}}}_{\boldsymbol{k}}| = |\hat {\theta }_{\boldsymbol{k}}|$
for all
$\boldsymbol{k}$
. As a result, they share the same power spectrum
\begin{align} E(k) = \dfrac {1}{2}\!\!\!\!\!\!\sum _{\quad k\leqslant |{\boldsymbol{k}}|\lt k+1}\!\!\!\!\!\!\!\!\!\!|\hat {\theta }_{\boldsymbol{k}}|^2 = \dfrac {1}{2}\!\!\!\!\!\!\sum _{\quad k\leqslant |{\boldsymbol{k}}|\lt k+1} \!\!\!\!\!\!\!\!\!\! |\hat {{\boldsymbol{u}}}_{\boldsymbol{k}}|^2; \end{align}
as well as their spatially averaged square values and gradient norms
\begin{equation} {\unicode{x2A0F}}\theta ^2 = {\unicode{x2A0F}} |{\boldsymbol{u}}|^2 \quad \text{and} \quad {\unicode{x2A0F}} |\boldsymbol{\nabla } \theta |^2 = {\unicode{x2A0F}} |\boldsymbol{\nabla } {\boldsymbol{u}} |^2, \end{equation}
where we adopt the shorthand
${\unicode{x2A0F}} = \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 ,{\boldsymbol{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
\begin{align} & \hat F_{\boldsymbol{k}} = \sqrt {\dfrac {2 \mathcal I}{c_d\pi k_{\!{f}}^2}}\,\textrm{e}^{-\frac {|{\boldsymbol{k}}|^2}{2k_{\!{f}}^2}}\,\hat \eta _{{\boldsymbol{k}}} - \alpha \,\hat \theta _{{\boldsymbol{k}}}\, {\mathbf 1}_{|{\boldsymbol{k}}| \leqslant k_{\!{f}}}. \end{align}
Here,
$\eta _{{\boldsymbol{k}}}$
are random complex Gaussian coefficients, with statistics prescribed by
$\mathbb{E} [\hat {\eta }_{{\boldsymbol{k}}}(t) ] = 0$
and
$\mathbb{E} [\hat {\eta }_{{\boldsymbol{k}}}(t)\hat {\eta }_{{\boldsymbol{k}}'}(t') ]=\delta _{{\boldsymbol{k}}+{\boldsymbol{k}}'}\,\delta (t-t')$
, where
$\mathbb{E}[\boldsymbol{\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
${\boldsymbol{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 _{{\boldsymbol{k}}} |\hat F_{\boldsymbol{k}}|^2 = \mathcal I - ({\alpha }/{2}) \sum _{|{\boldsymbol{k}}| \leqslant k_f} |\hat \theta _{{\boldsymbol{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 _{|{\boldsymbol{k}}| \leqslant k_{\!{f}}} \hat \theta _{\boldsymbol{k}}\,\textrm{e}^{{i}\,{\boldsymbol{k}} \boldsymbol{\cdot } {\boldsymbol{x}}}$
and
$\psi _{\!{f}} = \sum _{|{\boldsymbol{k}}| \leqslant k_{\!{f}}} |{\boldsymbol{k}}|^{-1} \hat \theta _{\boldsymbol{k}}\,\textrm{e}^{{i}\,{\boldsymbol{k}} \boldsymbol{\cdot } {\boldsymbol{x}}}$
, the following global balances for the energy and the Hamiltonian hold:
\begin{equation} \mathcal I = \varepsilon +\alpha \, {\big \langle {\theta ^2_{\!{f}}}\big \rangle } \quad \text{ and }\quad k^{-1}_{\mathcal I} {\mathcal I}= \gamma + \alpha \, {\langle {\theta _{\!{f}}\, \psi _{\!{f}}}\rangle }, \end{equation}
where
$\varepsilon := \nu {\left \langle {|\boldsymbol{\nabla } \theta |^2}\right \rangle }$
,
$\gamma := \nu {\left \langle { \boldsymbol{\nabla } \theta \boldsymbol{\cdot } \boldsymbol{\nabla } \psi }\right \rangle }$
and the angular brackets
$\left \langle {\boldsymbol{\cdot }}\right \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 (Reference Boffetta and Ecke2012) on the 2-D Navier–Stokes equations, we introduce the characteristic SQG damping and viscous centroid wavenumbers as, respectively,
\begin{equation} k_\alpha = \dfrac {\big \langle {\theta _{\!{f}}^2}\big \rangle }{\langle {\theta _{\!{f}} \,\psi _{\!{f}}} \rangle } \quad \text{and} \quad k_\nu = \dfrac {\varepsilon }{\gamma }. \end{equation}
On phenomenological grounds, we expect
$k_\nu \gt k_{\mathcal I}\gt k_\alpha$
, 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
\begin{equation} k_\nu \gg k_{\mathcal I} \gtrsim k_\alpha , \text{ with } k_I,k_\alpha = O(1). \end{equation}
Under this assumption, the steady-state balances (2.5), combined with the definitions (2.6), yield the following asymptotic estimates:
\begin{equation} \varepsilon \simeq \left (1-\dfrac {k_\alpha }{k_{\mathcal I}}\right ){\mathcal I},\quad {\big \langle {\theta _{\!{f}}^2}\big \rangle } \simeq \dfrac {k_\alpha }{\alpha k_{\mathcal I}} {\mathcal I} \quad \text{and } \gamma \simeq \dfrac {\varepsilon }{k_\nu }. \end{equation}
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 \gt k_{\!{f}}} E(k) \lt \infty$
. Equation (2.8) also implies that the viscous dissipation of the Hamiltonian vanishes in the limit
$\nu \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_{\mathcal 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 _{\mathcal I}$
and dissipation
$\tau _\nu$
as the turnover times at wavenumbers
$k_{\mathcal I}$
and
$k_\nu$
, respectively. The dissipative scales
$k_\nu$
and
$\tau _\nu$
a priori differ from the usual Kolmogorov scales
$k_\eta = (\varepsilon /\nu ^3)^{1/4}$
,
$\tau _\eta = (\nu /\varepsilon )^{1/2}$
found in 3-D turbulence. We postpone further comments on this issue to § 2.3. We denote
$\tau _\alpha =1/\alpha$
as the damping time scale. The SQG turbulence should correspond to a situation where
$\tau _\alpha \gg \tau _{\mathcal I}$
, in order to prevent large-scale fluctuations from dominating the turbulent statistics of the direct cascade (see Homann, Bec & Grauer Reference Homann, Bec and Grauer2013 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 Reference Orszag1971), 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
${\mathcal I}=1$
, and a damping rate
$\alpha =0.1$
, yielding an injection scale of
$\ell _{\mathcal I} = 2\pi /k_{\mathcal I} \simeq 2.8$
. To address finite- diffusivity effects, we define the Péclet number as
\begin{equation} \textit{Pe} := \dfrac {\ell _{\mathcal I} \big \langle {\theta _{\!{f}}^2}\big \rangle ^{1/2}}{\nu }. \end{equation}
The Péclet number in SQG is expected to play a role analogous to the Reynolds number in 3-D turbulence (Sukhatme & Pierrehumbert Reference Sukhatme and Pierrehumbert2002). Table 1 summarises the key 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.
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$
.
The infrared damping can be quantified in terms of a damped Péclet number: here, we define
$\textit{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_{\textit{max}} = N/3$
is chosen to satisfy
$k_{\textit{max}} \simeq \,k_\eta$
, following good practices in numerical studies of 3-D turbulence (Grauer, Homann & Pinton Reference Grauer, Homann and Pinton2012; Iyer, Sreenivasan & Yeung Reference Iyer, Sreenivasan and Yeung2015; Buaria et al. Reference Buaria, Pumir, Bodenschatz and Yeung2019). 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_{\textit{max}} \simeq 0.4 k_\eta$
.
Regarding simulation durations, runs I and II consist of a single realisation over
$T_{\textit{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_{\textit{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
$\textit{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}){\unicode{x2A0F}} \theta ^2$
and the Hamiltonian
$({1}/{2}){\unicode{x2A0F}} \theta \psi$
, as well as their instantaneous dissipation rates
$\nu {\unicode{x2A0F}} |\boldsymbol{\nabla } \theta |^2$
and
$\nu {\unicode{x2A0F}} \boldsymbol{\nabla } \theta \boldsymbol{\cdot } \boldsymbol{\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 _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.
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}}^+}\, {\textrm d}\tau \langle {f(\boldsymbol{\cdot })g(\boldsymbol{\cdot }+\tau )}\rangle / \langle { fg}\rangle$
.
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 the dissipation rate
$\nu {\unicode{x2A0F}} |\boldsymbol{\nabla } \theta |^2$
and the energy
$({1}/{2}){\unicode{x2A0F}} \theta ^2$
(inset of the right panel) reflects the direct cascade mechanism. The cross-correlation peaks at
$\tau \simeq 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
$\textit{Pe}$
. Despite significant statistical fluctuations, our numerical results suggest that the energy
$\mathcal E$
and the Hamiltonian
$\mathcal H$
converge towards asymptotic values approximately given by
\begin{equation} {\mathcal E} \; \simeq 2.7\; \ell _{\mathcal I}^{2/3} {\mathcal I}^{2/3},\quad {\mathcal H} \; \simeq 0.4\; \ell _{\mathcal I}^{ 5/3 } {\mathcal I}^{2/3}. \end{equation}
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
$\textit{Pe}$
. While
$\gamma$
algebraically vanishes as
$\textit{Pe}$
increases,
$\varepsilon$
asymptotes to a finite value
\begin{equation} \varepsilon \simeq 0.34 {\mathcal I},\quad \gamma \propto \textit{Pe}^{-3/4}. \end{equation}
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 \simeq 0.65 k_{\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
$\textit{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.
Figure 3. (a): averaged values of the Hamiltonian
$\mathcal H$
and energy
$\mathcal E$
as a function of
$\textit{Pe}$
. (b): corresponding dissipation rates,
$\gamma$
and
$\varepsilon$
, 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
$\textit{Pe}\gt 5\times 10^4$
and the error bars indicate twice the standard deviations of the statistical ensembles.
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
$\textit{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
\begin{equation} \varepsilon ^{(2q)}:={\nu \,\big \langle {|\boldsymbol{\nabla }(\theta ^{q})|^2}\big \rangle } ={q^2\,\big \langle {\theta ^{2q-2}}\big \rangle {\mathcal I}} -\dfrac {\alpha q^2}{2q-1}{\,\big \langle {\theta ^{2q-1}\theta _{\!{f}}}\big \rangle }. \end{equation}
Similarly to the quadratic case, the dissipation rates exhibit anomalous behaviours: one observes
$\varepsilon ^{(4)}\simeq 1.9 {\mathcal I} \langle {\theta ^2} \rangle$
,
$\varepsilon ^{(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}}\,\varepsilon ^{2/3}k^{-5/3}$
for
$k_{\mathcal I} \ll k \ll k_\nu$
, with an increasing range as
$\textit{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_\nu$
and
$k_\eta$
, share a similar scaling
$\propto \textit{Pe}^{3/4}$
with
$k_\eta /k_\nu \simeq 10$
. This supports the interpretation of
$k_\nu$
as the dominant dissipation scale. To summarise, the SQG K41 picture is mediated by the identifications, holding up to possibly non-universal constant factors
\begin{equation} \varepsilon \propto {\mathcal I},\quad k_\nu \propto k_\eta , \quad \big\langle {\theta _{\!{f}}^2}\big\rangle \propto \ell _{\mathcal I}^{2/3}{\mathcal I}^{2/3}. \end{equation}
These identifications lead to
$\gamma \simeq \varepsilon /k_\nu \propto {\mathcal I}/k_\eta \propto \textit{Pe}^{-3/4}$
, consistent with the algebraic behaviour reported in figure 3.
Figure 5. (a): energy spectra observed from the six runs; the black dashed line indicates the Kolmogorov prediction
$C_{{K}}\varepsilon ^{2/3} k^{-5/3}$
, with a Kolmogorov constant
$C_{{K}}\simeq 0.35$
. (b): Péclet number dependence of the viscous wavenumber
$k_\nu$
and Kolmogorov wavenumber
$k_\eta$
. See the text for definitions.
3. Homogenous isotropic SQG turbulence
From the spectral picture of K41 theory, one might think that the statistics of
$\theta$
and
$\boldsymbol{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
$\boldsymbol{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
$\boldsymbol{u}$
– focusing in particular on the various third-order structure functions.
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 _{\mathcal 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 (Reference Juckes1995) and Appendix A.1.
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.4 u_1}$
. (b): time series of the spatially averaged flatnesses defined as
$\mathcal{F}_\theta (t) = {\unicode{x2A0F}} \theta ^4 / ({\unicode{x2A0F}} \theta ^2)^2$
,
$\mathcal{F}_{u_1}(t) = {\unicode{x2A0F}} u_1^4 / ({\unicode{x2A0F}} u_1^2)^2$
, along with representative fields configurations observed during two extreme fluctuations.
3.2. Yaglom’s law
To investigate two-point spatial statistics, we define the increments of the velocity and scalar fields as
$\delta {\boldsymbol{u}} = {\boldsymbol{u}}({\boldsymbol{x}}+{\boldsymbol \ell },t) - {\boldsymbol{u}}({\boldsymbol{x}},t)$
and
$\delta \theta = \theta ({\boldsymbol{x}}+{\boldsymbol \ell },t) - \theta ({\boldsymbol{x}},t)$
. The longitudinal and transverse increments are denoted as
$\delta u_\parallel =\delta {\boldsymbol{u}} \boldsymbol{\cdot }\hat {\boldsymbol \ell }$
and
$\delta u_\perp = \delta {\boldsymbol{u}} \boldsymbol{\cdot }\hat {\boldsymbol \ell }^\perp$
, with
$\hat {{\boldsymbol \ell }} = {\boldsymbol \ell }/|{\boldsymbol \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({\boldsymbol \ell }) := \langle \delta \boldsymbol{u}\,(\delta \theta ) ^2 \rangle$
, which takes the form of the (vectorial) balance equation
\begin{align} &\boldsymbol{S}_3({\boldsymbol \ell }) = -2 \varepsilon {\boldsymbol \ell } \left [1 -\mathfrak C_{\!{f}}(\ell )-\mathfrak C_\nu (\ell )\right ]\!,\quad \mathfrak C_{{{f}},\nu }(\ell ):= \dfrac {2}{\ell ^2} \int _0^\ell \ell ' C_{{{f}},\nu }(\ell ')\,\textrm{d}\ell '\nonumber \\ &\quad \text{with } C_{\!{f}}(\ell ):=\dfrac {\alpha }{\varepsilon } \left [ \left \langle {\theta (\mathbf{0})\theta _{\!{f}}(\mathbf{ 0})}\right \rangle -\left \langle {\theta (\mathbf{0})\theta _{\!{f}}({\boldsymbol \ell })}\right \rangle \right ] \text{ and } C_\nu (\ell ):=\dfrac {\nu }{\varepsilon } \left \langle {\boldsymbol{\nabla } \theta (\mathbf{0})\boldsymbol{\cdot } \boldsymbol{\nabla } \theta ({\boldsymbol \ell })}\right \rangle\! . \end{align}
Figure 7(a) shows the behaviour of these correlation terms. The inertial range is defined as the range of scales where the correction terms
$\mathfrak C_{\!{f}}$
and
$\mathfrak C_\nu$
are negligible, that is
$\lambda _{\boldsymbol{\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 _\nu$
and
$\ell _{\mathcal I}$
, respectively. Within the inertial range, (3.1) prescribes the scaling of the third-order longitudinal and transverse mixed structure functions
\begin{equation} S^\parallel _3(\ell ) := \big \langle {\delta u_\parallel (\delta \theta )^2}\big \rangle = -2\varepsilon \ell \quad \text{and} \quad S^\perp _3(\ell ) := \big \langle {\delta u_\perp (\delta \theta)^2}\big \rangle = 0. \end{equation}
Figure 7(b) shows the longitudinal component of the structure function
$S^\parallel _3$
, compensated by
$-2\varepsilon \ell$
for all six runs. The results reveal convergence towards an asymptotic profile, with a plateau in the range
$\ell \gg \ell _\nu$
that expands as the Péclet number
$\textit{Pe}$
increases. At small scales
$\ell \lesssim 0.5\ell _\nu$
, the profile is dominated by viscous corrections, captured by
$1-\mathfrak C_\nu (\ell )$
. At scales
$\ell \gtrsim 20 \ell _\nu \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 _\nu \lesssim \ell \lesssim 20 \ell _\nu$
, 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.
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\varepsilon \ell )$
, as a function of the normalised separation
$\ell /\ell _\nu$
for the various runs. The black dotted line represents the viscous correction
$1-\mathfrak C_\nu (\ell )$
in (3.1), measured for run VI. The diamonds indicate the perpendicular contribution to Yaglom’s law.
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. 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
$\boldsymbol{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\varepsilon \ell$
.
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 \varepsilon \ell$
. The shaded area on the left is indicative of the inertial range for the highest-resolution run – see text for details.
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
$\boldsymbol{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
\begin{equation} \zeta _3^{u_\parallel }(\ell ) = \dfrac {{d} \log S_3^{u_\parallel }}{{d} \log \ell }. \end{equation}
Figure 9(a) reveals that in the SQG scaling range, all our runs exhibit a monotonous decay of
$\zeta _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
$\zeta _3^{u_\parallel } \simeq 1.1$
over the range
$5 \ell _\nu \lesssim \ell \lesssim 50 \ell _\nu$
, 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 \overset {\textit{la}w}\sim -\partial _1 u_1$
(see also Appendix A.2.2).
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 no reflect any type of inverse cascade. As this seemingly contradicts the earlier interpretation of Capet et al. (Reference Capet, Klein, Hua, Lapeyre and Mcwilliams2008), 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. Reference Capet, Klein, Hua, Lapeyre and Mcwilliams2008). This contribution arises from the integral identity (1.2) tying
$\theta$
to
$\boldsymbol{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 (Reference Bernard1999) for 2-D Navier–Stokes equations. The formula is most conveniently written in complex notations, promoting in particular the vorticity
$\omega =\boldsymbol{\nabla }^\perp \boldsymbol{\cdot } {\boldsymbol{u}}$
and the scalar
$\theta$
to functions
$\varOmega , \varTheta$
defined on the complex plane
$\mathbb C$
. Writing
$z = x_1 + i x_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
\begin{equation} \dfrac {1}{3}\partial _{\xi \xi }\big\langle (\delta U)^3\big\rangle + \langle A(0) \varOmega (\xi )\rangle = \partial _{\xi \xi ^*}\big\langle \delta U (\delta \varTheta )^2 \big\rangle . \end{equation}
Bernard’s formula involves the complex field
$A = a_1 + ia_1$
, associated with the vector field
$\boldsymbol{a} := (-\varDelta )^{1/2} [{\boldsymbol{u}} \theta ]- {\boldsymbol{u}} \omega$
.
Although not immediately evident, the field
$\boldsymbol{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
$\left \langle {\delta U (\delta \varTheta )^2}\right \rangle$
into a single-field contribution
$ \left \langle {(\delta U)^3}\right \rangle$
and a correction term
$\langle A(0) \varOmega (\xi )\rangle$
. It is a physical-space counterpart of a decomposition proposed by Capet et al. (Reference Capet, Klein, Hua, Lapeyre and Mcwilliams2008) 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 ,{\boldsymbol{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
$\left \langle {\delta U (\delta \varTheta )^2}\right \rangle = -2 \varepsilon \xi$
, implying that (3.4) there reduces to
$\partial _{\xi \xi }\langle (\delta U)^3\rangle = - 3\langle A(0) \varOmega (\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 longitunal structure functions, which takes the explicit differential form
\begin{align} & \mathcal{L}_{u_\parallel } S_3^{u_\parallel } = \mathcal{L}_{\parallel } S_3^\parallel + \left \langle {a_\parallel (\ell \boldsymbol{e}_1)\,\omega (\boldsymbol{0})}\right \rangle\! ,\nonumber \\ \text{with } \mathcal{L}_{u_\parallel } f := -\dfrac {1}{12}\ell ^{-2}&\partial _{\ell } \big [\ell ^{-1} \partial _{\ell }\big [\ell ^5\partial _{\ell }\big ( \ell ^{-1} f \big )\big]\big], \quad \mathcal{L}_{\parallel } f= \dfrac {1}{4}\partial _\ell \big [\ell ^{-1}\partial _\ell (\ell f)\big], \end{align}
naturally involving the vorticity
$\omega$
and the ageostrophic vector field
$\boldsymbol{a}$
. Again, we observe
$\mathcal{L}_{\parallel } [-2\varepsilon \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 \boldsymbol{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
\begin{equation} \langle a_\parallel (\ell \boldsymbol{e}_1) \omega \rangle \sim \dfrac {\varepsilon }{\ell }\log \dfrac {\ell _a}{\ell },\quad \ell \lt \ell _a \simeq 10 \ell _\nu \end{equation}
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}^{-1}_{u_\parallel } [\langle a_\parallel (\ell \boldsymbol{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
$\boldsymbol{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 Reference Pedlosky2013; Vallis Reference Vallis2017). 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
$\boldsymbol{a}$
is a particular example of such ageostrophic corrections; it endows SQG with 3-D-like features.
Figure 10. Snapshots of the scalar field
$\theta$
(a), vorticity
$\omega = \boldsymbol{\nabla }^\perp \boldsymbol{\cdot }{\boldsymbol{u}}$
(b) and the divergence of the ageostrophic term
$\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{a}}$
(c), normalised to unit variance and zoomed over
$(1/16)^2$
of the computational domain. (Data from run V.)
One effect is physical:
$\boldsymbol{a}$
is involved in vorticity stretching. This becomes clear when examining the equation for vorticity
$\omega = \boldsymbol{\nabla }^\perp \boldsymbol{\cdot } {\boldsymbol{u}}=(-\varDelta )^{1/2} \theta$
, directly derived from the SQG dynamics (1.1) as
\begin{equation} \partial _t \omega + {\boldsymbol{u}} \boldsymbol{\cdot } \boldsymbol{\nabla } \omega = -\boldsymbol{\nabla } \boldsymbol{\cdot } {\boldsymbol{a}}+ \nu \Delta \omega + (-\varDelta )^{1/2} F,\quad {\boldsymbol{a}}:= (-\varDelta )^{1/2} [{\boldsymbol{u}}\,\theta ]- {\boldsymbol{u}}\,\omega . \end{equation}
Here, the vector field
$\boldsymbol{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). With the presence of a non-vanishing
$\boldsymbol{\nabla } \boldsymbol{\cdot } {\boldsymbol{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){\unicode{x2A0F}} \omega ^2$
is not conserved in SQG, and its evolution is governed by
\begin{equation} \dfrac {\text{d}}{\text{dt}} \dfrac {1}{2}{\unicode{x2A0F}} \omega ^2 = {\unicode{x2A0F}} {\boldsymbol{a}} \boldsymbol{\cdot }\boldsymbol{\nabla } \omega -\nu {\unicode{x2A0F}} |\boldsymbol{\nabla } \omega |^2 + {\unicode{x2A0F}} \omega (-\varDelta )^{1/2} F, \end{equation}
where the term
${\unicode{x2A0F}} {\boldsymbol{a}} \boldsymbol{\cdot }\boldsymbol{\nabla } \omega = -{\unicode{x2A0F}}\omega \boldsymbol{\nabla } \boldsymbol{\cdot }{\boldsymbol{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
$\boldsymbol{\nabla } \boldsymbol{\cdot } {\boldsymbol{a}}$
and the SQG fields. In particular,
$\boldsymbol{\nabla } \boldsymbol{\cdot } {\boldsymbol{a}}$
exhibits strong values along the sharp fronts of a scalar filament, while remaining nearly zero elsewhere.
Another effect is geophysical, as
$\boldsymbol{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
\begin{equation} \left (\partial _{11}+\partial _{22} +\partial _{33}\right ) \psi = 0, \quad \left . \partial _3\psi \right |_{x_3=0}= \theta ,\quad \psi \underset {x_3\to -\infty }{\to } 0. \end{equation}
This system arises as a first-order perturbation of the primitive equations under the assumption of vanishing interior potential vorticity (Held et al. Reference Held, Pierrehumbert, Garner and Swanson1995). Comparing (3.7) with formula (3.a) in Held et al. (Reference Held, Pierrehumbert, Garner and Swanson1995), we identify
$\boldsymbol{\nabla }\boldsymbol{\cdot } {\boldsymbol{a}} = -\partial _3 u_3|_{x_3=0}$
.
In other words, and up to a rotational component,
$\boldsymbol{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
$\boldsymbol{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 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 _\nu$
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_\parallel$
displaying noticeable skewness.
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.
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_\parallel$
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_\parallel }$
, 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 12. Typical snapshots of the squared velocity gradient (a) and scalar gradient (b), normalised to unit variance, for run VI.
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$
\begin{equation} S_p^{\phi }=\left \langle {|\delta \phi |^p}\right \rangle \propto \ell ^{\zeta _p^\phi },\quad \delta \phi := \phi ({\boldsymbol{x}}+\ell \mathbf{e}_{\boldsymbol{1}}) -\phi ({\boldsymbol{x}}) ,\end{equation}
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.
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
$\textit{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.
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 _\nu \simeq 0.1\ell _{\mathcal I}$
, suggesting that no inertial-range power-law scaling regime holds for any of the three SQG fields. The insets allow us to assess the
$\textit{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
$\textit{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 (Reference Lapeyre2017) 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. Reference Benzi, Biferale, Ciliberto, Struglia and Tripiccione1996), 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 \simeq 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. The mixed structure functions associated with inviscid conservation laws of higher orders are defined by
\begin{equation} S_p^\parallel = \big\langle \delta u_\parallel |\delta \theta |^{p-1} \big\rangle . \end{equation}
For odd
$p=2q+1\gt 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
\begin{equation} S_p^\parallel \propto \ell ^{\zeta ^{\parallel }_p}. \end{equation}
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
$\zeta ^{\parallel }_p$
.
Figure 14. (a): local scaling exponents
$\zeta _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.
This type of scaling laws is similar to those observed in passive scalar turbulence (Gauding, Danaila & Varea Reference Gauding, Danaila and Varea2017). In passive scalar turbulence, however, Watanabe & Gotoh (Reference Watanabe and Gotoh2004) also identified scaling ranges for the family of structure functions constructed as
$S_p^{\textit{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^{\textit{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 Reference Benzi and Toschi2023). 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 observed. This strategy was introduced by Benzi et al. (Reference Benzi, Ciliberto, Baudet, Chavarria and Tripiccione1993a ,Reference Benzi, Ciliberto, Tripiccione, Baudet, Massaioli and Succi 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
\begin{equation} S_p^\phi \propto \big(S_2^\phi\big)^{C^\phi _{p}} \quad (\phi =\theta , u_\parallel , u_\perp ), \end{equation}
with the exponents
$C^\phi _p$
a priori different for each SQG field. The presence of scaling ranges in the usual sense (4.1) implies ESS with
$C_{p}^\phi = \zeta ^\phi _p/\zeta ^\phi _2$
. 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 ) ]^{\zeta ^\phi _p}$
, ESS would still hold, with the exponent
$C_{p}^\phi = \zeta ^\phi _p/\zeta ^\phi _2$
characterising the dominant power-law contribution (see Benzi et al. Reference Benzi, Ciliberto, Baudet, Chavarria and Tripiccione1993a
).
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
$\zeta _p^{u_\parallel }$
normalised by
$C_p$
. The shaded area there indicates the inertial range identified in figure 8 .
The ESS analysis of our SQG runs is summarised in figure 15, showing the scale-dependent ESS exponents
$C_p^\phi (\ell )= \zeta ^\phi _p(\ell )/\zeta ^\phi _2(\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 _\nu \lesssim \ell \lesssim 40 \ell _\nu$
, extending beyond the range over which Yaglom’s law is observed. figure 15(d) points out that when normalised by their ESS exponent
$C^{u_\parallel }_p$
, the scale-dependent exponents
$\zeta ^{u_\parallel }_p(\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^\phi _p})^{1/C^\phi _p}$
are independent of
$p$
.
4.2.3. Refined similarity
The RSS theory of Kolmogorov (Reference Kolmogorov1962) for 3-D turbulence ties deviations from self-similarity to heterogeneities in the scaling behaviours for the coarse-grained dissipation field
\begin{equation} \varepsilon _\ell (\boldsymbol{x},t) := \frac {4}{\pi \ell ^2} \int _{|\boldsymbol{r}|\lt \ell /2}\varepsilon (\boldsymbol{x} + \boldsymbol{r}) d^2\boldsymbol{r},\quad \text{with } \varepsilon ({\boldsymbol{x}},t) = \nu |\boldsymbol{\nabla }\theta |^2. \end{equation}
Refined self-similarity suggests the phenomenology
$\left \langle {\delta \phi ^p}\right \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 Reference Frisch1995), so that
\begin{equation} \left \langle {\delta \phi ^p}\right \rangle \propto \ell ^{\zeta _p^\phi } \ \text{ with } \ \zeta _p^\phi = p/3+\chi _{p/3}. \end{equation}
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 ) = \textrm{d} \log \langle \varepsilon _\ell ^{p/3}\rangle /\textrm{d}\log \ell$
. The exponents remain approximately constant across all scales for
$p\leqslant 5$
, and vary by at most 10 % within the inertial range for
$p=6$
and
$7$
. The shaded bands represent the exponents
$\zeta ^\parallel _p$
defined in § 4.2.1 obtained from mixed structure functions.
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.
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 _\nu \leqslant \ell \leqslant 10\ell _\nu$
. 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
\begin{equation} \langle \log (\varepsilon _\ell /\varepsilon ) \rangle = \mu _{62}\log (\ell /\ell _\nu ), \quad \textrm{Var}\,[\log (\varepsilon _\ell /\varepsilon )] = -2 \mu _{62}\log (\ell /\ell _\nu ), \end{equation}
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 ^{3D}_{62} \simeq 0.13$
commonly measured (Iyer et al. Reference Iyer, Sreenivasan and Yeung2015; Chevillard et al. Reference Chevillard, Garban, Rhodes and Vargas2019). 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
$\zeta ^\parallel _2 C_p^\phi$
, with
$\zeta _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.
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]$
.
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} \simeq 0.16\pm 0.02$
. We observe no sign of saturation for the
$\zeta _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 Reference Iyer, Sreenivasan and Yeung2020; Müller & Krstulovic Reference Müller and Krstulovic2025), and passive scalar studies (Celani et al. Reference Celani, Lanotte, Mazzino and Vergassola2000, Reference Celani, Lanotte, Mazzino and Vergassola2001; Iyer et al. Reference Iyer, Schumacher, Sreenivasan and Yeung2018), 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
\begin{equation} \big \langle {\delta u_\parallel |\delta \theta |^{p-1}}\big \rangle \sim \ell ^{p/3} \big\langle \varepsilon _\ell ^{p/3}\big\rangle, \end{equation}
entailing the identification
$\zeta _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 (Reference Duchon and Robert2000) might suggest the identification
$ \delta u (\delta \theta )^2 \equiv -2\varepsilon _\ell \ell$
entailing rather
$S_p^{\textit{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 {\boldsymbol{u}}| \sim |\delta \theta | \sim \ell ^{1/3}$
. Using systematic numerical simulations with up to
$16\,384^2$
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 (Reference Sukhatme and Pierrehumbert2002) 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
$\boldsymbol{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. Reference Iyer, Schumacher, Sreenivasan and Yeung2018). -
(ii) Another consequence of the duality between
$\theta$
and
$\boldsymbol{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. Reference Valadão, Ceccotti, Boffetta and Musacchio2024).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 Reference Watanabe and Iwayama2007), 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
$\textit{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 Reference Drivas and Eyink2019; Iyer Reference Iyer2023), and a similar behaviour has been reported in the highest-to-date resolved DNS of 3-D NS flows (Iyer et al. Reference Iyer, Drivas, Eyink and Sreenivasan2025). 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.
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 (Reference Juckes1995) for the relevant stability analysis and to figure 18 for a sketch of those specific SQG flows in the inviscid case.
Figure 18. Filaments (a) and vortices (b) as two prototypes of inviscid SQG flows.
A.1.1. Filaments
The SQG filaments of width
$2\delta$
, say aligned along
$x_1$
, are prescribed by the scalar field
\begin{equation} \theta ({\boldsymbol{x}}) = \left \{\!\! \begin{array}{ll} \dfrac {U_0}{2} & \mbox{ } |x_2|\leqslant \delta ,\\[8pt] 0 & \mbox{otherwise}, \end{array} \right . \end{equation}
which in turn generate a horizontal velocity
${\boldsymbol{u}} ({\boldsymbol{x}}) = U_0\log \left | ({x_2 - \delta })/({x_2+\delta })\right | {\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(\boldsymbol{x}) \underset {}{\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
\begin{equation} \theta ({\boldsymbol{x}}) = \left \{\!\! \begin{array}{ll} \theta _0, & \mbox{} |\boldsymbol{x}|\leqslant r_0, \\[6pt] 0 & \mbox{otherwise}, \end{array} \right . \end{equation}
and prescribe the circular flow
\begin{equation} {\boldsymbol{u}}({\boldsymbol{x}}) = u_\perp \hat {\boldsymbol{x}}^\perp ,\quad u_\perp := \theta _0 \int ^{+\infty }_{0} \mathcal{J}_1\left (v \frac {{\boldsymbol{x}}}{r_0}\right )\mathcal{J}_1\left (v\right ) \, \textrm{d}v ,\end{equation}
where
$\mathcal{J}_{0}$
and
$\mathcal{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
\begin{equation} u_\perp ({\boldsymbol{x}}) \underset {r_0^\pm }{\propto } \log | r_0-|{\boldsymbol{x}}|\,|,\quad u_\perp ({\boldsymbol{x}}) \underset {+\infty }{\propto } |\boldsymbol{x}|^{-2}. \end{equation}
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 (Reference Frisch1995), we list below the symmetries holding in particular in the inviscid case:
-
(i) time or space translation:
$t\mapsto t + \tau ,\quad {\boldsymbol{x}} \mapsto {\boldsymbol{x}} + {\boldsymbol{r}}$
with
$\tau \in \mathbb{R}$
and
${\boldsymbol{r}} \in \mathbb{R}^2$
; -
(ii) rotation and reflection:
${\boldsymbol{x}} \mapsto \mathcal O{\boldsymbol{x}},\quad {\boldsymbol{u}} \mapsto \mathcal O {\boldsymbol{u}},\quad \theta \mapsto det(\mathcal{O})\theta ,$
with
$\mathcal O \in O(2)$
. Pure rotations correspond to
$\mathcal O \in SO(2)$
, characterised by
$\det \mathcal O=1$
. Reflections correspond to
$\mathcal O \in O(2) \backslash SO(2)$
, with
$\det \mathcal O=-1$
; -
(iii) Time or space rescaling:
$t\mapsto \alpha t,\quad {\boldsymbol{x}} \mapsto \lambda {\boldsymbol{x}},\quad {\boldsymbol{u}} \mapsto (\lambda /\alpha ) {\boldsymbol{u}},\quad \theta \mapsto (\lambda /\alpha ) \theta$
with
$\alpha , \lambda \gt 0$
; -
(iv) Galilean invariance:
$\theta \mapsto \theta + \theta _0,\quad {\boldsymbol{u}} \mapsto {\boldsymbol{u}} + {\boldsymbol{u}}_0, \quad {\boldsymbol{x}} \mapsto {\boldsymbol{x}} + t{\boldsymbol{u}}_0$
with
$\theta _0 \in \mathbb{R}, {\boldsymbol{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 Reference Jackson2021), ensuring the non-vanishing asymptotics
${\boldsymbol{u}} \sim {\boldsymbol{u}}_0$
,
$\psi \sim - {\boldsymbol{x}}^\perp \boldsymbol{\cdot }{\boldsymbol{u}}_0$
as
$|{\boldsymbol{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
${\boldsymbol{u}}({\boldsymbol{x}},t),\theta ({\boldsymbol{x}},t)$
solves the inviscid SQG system, then so does the pair
\begin{equation} {\boldsymbol{u}}'({\boldsymbol{x}}',t') := (\lambda /\alpha ){\boldsymbol{u}}({\boldsymbol{x}}'/\lambda ,t'/\alpha ),\quad \theta '({\boldsymbol{x}}',t'):=(\lambda /\alpha )\theta '({\boldsymbol{x}}'/\lambda ,t'/\alpha ). \end{equation}
A.2.2. Statistical symmetries
The statistical symmetries are statistical formulation of the classical symmetries, which hold at the level of the turbulent averages
$\left \langle {.}\right \rangle$
. They mean
\begin{equation} {\boldsymbol{u}}'({\boldsymbol{x}}',t') , \theta '({\boldsymbol{x}}',t') \overset {\textit{la}w}{\sim } {\boldsymbol{u}}({\boldsymbol{x}},t) , \theta ({\boldsymbol{x}},t), \end{equation}
where the transforms
${\boldsymbol{x}},t,{\boldsymbol{u}},\theta \mapsto {\boldsymbol{x}}',t',{\boldsymbol{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
(A7)where
\begin{equation} p_\theta (y) = p_\theta (-y),\quad p_{u_i}(y) = p_{u_i}(-y), \end{equation}
$p_\theta (y)= \left \langle {\delta (\theta ({\boldsymbol{x}}_0)-y)}\right \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 ({\boldsymbol{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
(A8)implying that the odd order moments vanish:
\begin{equation} p_{\delta \theta }(y) = p_{\delta \theta }(-y),\quad p_{\delta u_\perp }(y) = p_{\delta u_\perp }(-y), \end{equation}
$ \langle {\delta \theta ^{2p+1}} \rangle = \langle {\delta u_{\perp }^{2p+1}}\rangle =0$
.
-
(iii) Mixed structure functions
(A9)
\begin{equation} \forall p \in \mathbb N\, \left \langle {\delta u_{\parallel }\delta \theta ^{2p+1}}\right \rangle =\left \langle {\delta u_{\parallel }\delta u_\perp ^{2p+1}}\right \rangle =0,\quad \big \langle {{\boldsymbol{u}}({\boldsymbol{x}})\theta ^2}\big \rangle =0,\quad \big \langle {{\boldsymbol{u}}({\boldsymbol{x}}) \omega ^2}\big \rangle =0 .\end{equation}
The second identity stems from the observations
$\boldsymbol{\nabla }\boldsymbol{\cdot } \langle {{\boldsymbol{u}}({\boldsymbol{x}})\theta (0)^2} \rangle =0$
(incompressibility) and
$\boldsymbol{\nabla} ^\perp \boldsymbol{\cdot } \langle {{\boldsymbol{u}}({\boldsymbol{x}})\theta ^2} \rangle = \langle {\omega ({\boldsymbol{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 Reference Tao2016), 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+ i u_2$
, and promoting the vorticity and the streamfunction into the complex fields
$ \varOmega (z,z^*), \varPsi (z,z^*)$
, one checks the relations
\begin{equation} U(z,z*) = -2i \partial _{z^*} \varPsi ,\quad \varOmega = -2 i \partial _z U,\quad \varOmega =-4 \partial _{zz^*} \varPsi , \end{equation}
where the second identity uses incompressibility. We now define the structure functions
\begin{equation} S_3^U(\xi ,\xi ^\star ) = \big \langle \delta U^3 \big \rangle ,\quad S_3^{\varTheta }(\xi ,\xi ^\star )=\big \langle \big (\delta \varTheta \big )^2 \delta U \big \rangle ,\quad S_3^{\varOmega }(\xi ,\xi ^\star )=\big \langle \left (\delta \varOmega \right )^2 \delta U \big \rangle , \end{equation}
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
\begin{equation} S_3^U= 6\,\big\langle {U^2U'}\big \rangle ,\quad S_3^\varTheta = 4\,\left \langle {U \varTheta \varTheta '}\right \rangle\! ,\quad S_3^\varOmega = 4\,\left \langle {U \varOmega \varOmega '}\right \rangle\!. \end{equation}
Let us simply emphasise the crucial role played by the mirror symmetry, which makes the correlators
$ \langle {U'\varOmega ^2} \rangle, \langle {U'\varTheta ^2} \rangle$
vanish! From (B2) and homogeneity, the following calculation string holds:
\begin{equation} \partial _{\xi \xi }\big \langle {(\delta U)^3}\big \rangle = -6\big \langle {(\partial _z'U')\partial _z(U^2)}\big \rangle = 3\left \langle {U\varOmega \varOmega '}\right \rangle . \end{equation}
Using (B3), we then recover the kinematic relation derived in Bernard (Reference Bernard1999), namely
\begin{equation} \partial _{\xi \xi }S_3^U = \dfrac {3}{4} S_3^\varOmega . \end{equation}
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,
\begin{equation} {\boldsymbol{a}} = (-\varDelta )^{1/2}\left ({\boldsymbol{u}}\theta \right )-{\boldsymbol{u}}\omega ,\quad A:= a_1+ia_2= 2\left |\partial _{z}\right |[U\varTheta ]-U\varOmega , \end{equation}
where the shorthand notation
$|\partial _{z}|$
stands for the complex counterpart of
$({1}/{2})(-\varDelta )^{1/2}$
, such that
$\varOmega =2|\partial _{z}| \varTheta$
and
$|\partial _{z}|^2 = -\partial _{zz*}$
.
We then observe
\begin{equation} \partial _{\xi \xi ^*}\left \langle {U\varTheta \varTheta '}\right \rangle = \left \langle {|\partial _{z}|\left [ U \varTheta \right ]|\partial _{z'}|\left [\varTheta '\right ]}\right \rangle = \dfrac {1}{4} \left \langle {\left (A+ U\varOmega \right ) \varOmega '}\right \rangle\! , \end{equation}
pointing out that the lack of a minus sign in the first equality comes from homogeneity. Combining (B7) and (B3) leads to
\begin{equation} \partial _{\xi \xi ^*}S_3^\varTheta = \left \langle {A \varOmega '}\right \rangle +\dfrac {1}{4} S_3^\varOmega . \end{equation}
Combining (B5) and (B8) finally yields
\begin{equation} \partial _{\xi \xi ^*}S_3^\varTheta = \left \langle {A \varOmega '}\right \rangle +\dfrac {1}{3} \partial _{\xi \xi }S_3^U, \end{equation}
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 (Reference Brachet2000), to observe that the the third-order-but-one-separation isotropic tensors
\begin{equation} B_{\textit{ijk}}({{\boldsymbol \ell }}) =\left \langle \delta u_i \delta u_j \delta u_k \right \rangle\! ,\quad \delta u_i:= u_i({{\boldsymbol \ell }})-u_i(0) ,\end{equation}
can be written as
\begin{equation} B_{\textit{ijk}}({{\boldsymbol \ell }}) = 2 \left (b_{ij,k}+b_{ki,j}+b_{jk,i}\right ), \quad b_{ij,k} := \left \langle u_i(0)u_j(0)u_k({{\boldsymbol \ell }})\right \rangle\! , \end{equation}
where homogeneity obviously implies the identity
$b_{ij,k} = - \langle {u_i({{\boldsymbol \ell }})u_j({{\boldsymbol \ell }})u_k(0)} \rangle$
. As an isotropic tensor,
$b_{\textit{ij}, k}$
takes the form
\begin{equation} b_{\textit{ij}, k} = C(\ell )\dfrac {\delta _{\textit{ij}}\ell _k}{\ell }+D(\ell )\dfrac {\delta _{ik} \ell _j+\delta _{jk} \ell _i}{\ell } +F(\ell )\dfrac {\ell _i\ell _j\ell _k}{\ell ^3}, \end{equation}
which with the incompressibility condition
$\partial _k b_{ij,k}=0$
imposes
\begin{equation} -2D = (d-1)C+rC'\quad \&\quad F= -C + \ell \partial _\ell C, \end{equation}
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:
$$\begin{align}\big \langle \delta u_\parallel ^3 \big \rangle & = B_{\textit{ijk}} \hat \ell _i\hat \ell _j\hat \ell _k = 6(1-d)C\end{align}$$
$$\begin{align}\big \langle \delta u_\parallel |{\delta \boldsymbol{u}}|^2 \big \rangle & = \hat \ell _iB_{\textit{ijj}} =2(1-d)\left ((1+d)C+\ell \partial _\ell C\right ),\end{align}$$
$$\begin{align}\big \langle \delta u_\perp |{\delta \boldsymbol{u}}|^2 \big \rangle & = \hat n_iB_{\textit{ijj}}=0.\end{align}$$
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^\varTheta$
defined in (B2) relate to their real counterpart as
\begin{align} & \ell S^U = \big (\big \langle {\|\delta u\|^2\delta u_\parallel }\big \rangle + i\big \langle {\|\delta u\|^2\delta u_\perp }\big \rangle \big )\xi ,\quad \ell S^\varTheta = \big (\big \langle {|\delta \theta |^2\delta u_\parallel }\big \rangle + i\big \langle {|\delta \theta |^2\delta u_\perp }\big \rangle \big )\xi . \end{align}
In line with (B14a ), isotropy and mirror symmetry make the perpendicular contributions vanish, such that the previous equation in fact reduces to
\begin{align} & \ell S^U = \big \langle {|\delta \boldsymbol{u}|^2\delta u_\parallel }\big \rangle \xi ,\quad \ell S^\varTheta = \big \langle {|\delta \theta |^2\delta u_\parallel }\big \rangle \xi . \end{align}
On the one hand, now combining (B13) and (B14a ), we obtain
\begin{equation} C(\ell ) = -\dfrac {1}{6}\big \langle \delta u_\parallel ^3\big \rangle ,\quad S_3^U = 6\dfrac {F}{\ell ^3}\xi ^3, \quad F:= \ell ^{2}\partial _\ell \big (\ell ^{-1}C \big ). \end{equation}
Using
$\partial _\xi f(\ell ) = \xi ^*\partial _\ell f/(2\ell )$
, this in turn yields the first piece of the differential formulation
\begin{align} & \partial _{\xi \xi }S_3^U = \dfrac {3\xi }{2\ell ^3}\partial _\ell \big (\ell ^{-1} \partial _\ell \ell ^3 F \big )=3\dfrac {\xi }{\ell } \mathcal L_{u_\parallel }\big [\big \langle {\delta u_\parallel ^3}\big \rangle \big ],\nonumber\\ &\quad\text{identifying } \mathcal{L}_{u_\parallel } := -\frac {1}{12}\ell ^{-2}\partial _{\ell } \big (\ell ^{-1} \partial _{\ell }\big (\ell ^5\partial _{\ell }\big ( \ell ^{-1} \boldsymbol{\cdot }\big )\big )\big ). \end{align}
On the other hand, we also compute from (B16), using
$\partial _{\xi ^*} f(\ell ) = \xi \partial _\ell f/(2\ell )$
\begin{equation} \partial _{\xi \xi ^*} S_3^\varTheta = \dfrac {\xi }{\ell } \mathcal L_M\big [\big \langle {|\delta \theta |^2\delta u_\parallel }\big \rangle \big ],\quad \mathcal{L}_{M}:= \dfrac {1}{4}\partial _\ell \big (\ell ^{-1}\partial _\ell \left ( \ell \boldsymbol{\cdot }\right )\big ). \end{equation}
Similar to (B15), we also have
\begin{equation} \left \langle {A\varOmega '}\right \rangle = -\left \langle {A'\varOmega }\right \rangle =-\frac {\xi }{\ell }\big ( \big \langle {a_\parallel (\ell \mathbf{e}_{\boldsymbol{1}})\, \omega (0)}\big \rangle +i\big \langle {a_\perp (\ell \mathbf{e}_{\boldsymbol{1}})\, \omega (0)}\big \rangle \big ). \end{equation}
Using (B9) with (B18), (B19) and (B20) yields
$\left \langle {a_\perp (\ell \mathbf{e}_{\boldsymbol{1}})\,\omega (0)}\right \rangle =0$
and the differential formulation (3.5).
Open access
