1. Introduction
Turbulence is a complex and chaotic phenomenon characterised by a large number of interacting degrees of freedom organised hierarchically across multiple scales of motion (Frisch Reference Frisch1995). Determining whether a turbulent flow remains predictable or, at each scale, retains any degree of predictability has been a longstanding challenge, tracing back to the pioneering works by Lorenz, Ruelle, Leith and Kraichnan (Lorenz Reference Lorenz1969; Leith Reference Leith1971; Leith & Kraichnan Reference Leith and Kraichnan1972; Ruelle Reference Ruelle1979; Deissler Reference Deissler1986). The nonlinear amplification of small-scale perturbations led to the formulation of the famous ‘butterfly effect’ (Lorenz Reference Lorenz1963). When extended to multiscale systems such as turbulence, these ideas give rise to the concept of the cascade of errors, in which small-scale perturbations progressively amplify and propagate to larger scales, gradually spoiling predictability at larger and larger scales (Boffetta et al. Reference Boffetta, Celani, Crisanti and Vulpiani1997; Rotunno & Snyder Reference Rotunno and Snyder2008; Boffetta & Musacchio Reference Boffetta and Musacchio2017; Palmer Reference Palmer2024).
A key mathematical tool for studying predictability is the Lyapunov exponent and its finite-time version (FTLE). Ruelle predicted that the Lyapunov exponent in turbulence is proportional to the inverse of the smallest time, the Kolmogorov time, and therefore it increases with the flow’s Reynolds number (
$\textit{Re}$
) (Ruelle Reference Ruelle1979). Nonetheless, a turbulent flow at large
$\textit{Re}$
remains predictable at large scales since the error grows with the characteristic turnover time of that scale, which is independent of the Reynolds number. This property is familiar in oceanography and atmospheric flows where the smallest time scale can be very small (seconds) (Garratt Reference Garratt1994) while the weather remains predictable for days.
In this work, we consider the surface quasi-geostrophic (SQG) equation, a model that describes the flow governed by the conservation of buoyancy at the surface of a rotating, stratified fluid (Blumen Reference Blumen1978; Pierrehumbert, Held & Swanson Reference Pierrehumbert, Held and Swanson1994). Beyond its geophysical relevance, the SQG model has gained attention in the fluid dynamics community due to its striking similarities to three-dimensional (3-D) Navier–Stokes (NS) turbulence while keeping some properties of two-dimensional (2-D) flows. In particular, the SQG model has two inviscid quadratic invariants, similar to 2-D turbulence (Celani et al. Reference Celani, Cencini, Mazzino and Vergassola2004; Lapeyre Reference Lapeyre2017; Valade, Thalabard & Bec Reference Valade, Thalabard and Bec2024) with one of the two, the surface potential energy, displaying, in the presence of forcing and dissipation, a direct cascade à la Kolmogorov towards the small scales similar to 3-D turbulence (Valadão et al. Reference Valadão, Ceccotti, Boffetta and Musacchio2024). Despite this similarity, previous numerical studies (Pierrehumbert et al. Reference Pierrehumbert, Held and Swanson1994; Ohkitani & Yamada Reference Ohkitani and Yamada1997; Sukhatme & Pierrehumbert Reference Sukhatme and Pierrehumbert2002) reported that the scaling exponent of the spectrum of the surface potential energy deviates from the Kolmogorov value
$-5/3$
predicted by dimensional arguments.
Based on very high-resolution direct numerical simulations of the SQG model, we study the predictability of the direct cascade in relation to the scaling properties of the surface buoyancy field. At large Reynolds numbers, we observe the recovery of the Kolmogorov scaling in the surface potential energy spectrum. Further, we measure the finite-time distribution of the Lyapunov exponent as a function of
$\textit{Re}$
and we find an anomalous scaling law in which the Lyapunov exponent grows faster than what is predicted on dimensional grounds, similar to what was observed in 3-D NS turbulence (Boffetta & Musacchio Reference Boffetta and Musacchio2017; Mohan, Fitzsimmons & Moser Reference Mohan, Fitzsimmons and Moser2017; Berera & Ho Reference Berera and Ho2018). Despite this anomaly, we find that the distribution of the FTLE follows an almost universal function, independent of the Reynolds number of the flow.
The remainder of this paper is organised as follows. In § 2 we review the basic definitions of the SQG equation and discuss its statistical properties in the turbulent regime. Section 3 discusses the main results obtained in this work through the use of extensive numerical simulations covering a large range of Reynolds numbers. We split the results into two subsections: § 3.1 explores the Reynolds dependence on the dimensional scaling properties of SQG, principally the scaling exponent of the surface potential energy spectrum; § 3.2 addresses the Eulerian predictability and the statistics of the FTLE as functions of Reynolds number. In § 4, we summarise our results, pointing out directions for future research.
2. The surface quasi-geostrophic model
The SQG model describes the large-scale dynamics of a rapidly rotating, stably stratified flow using a 2-D equation for the surface buoyancy field
$\theta ({\boldsymbol x}, t)$
(Juckes Reference Juckes1994; Pierrehumbert et al. Reference Pierrehumbert, Held and Swanson1994; Lapeyre & Klein Reference Lapeyre and Klein2006; Lapeyre Reference Lapeyre2017; Vallis Reference Vallis2017; Siegelman et al. Reference Siegelman2022):
\begin{equation} \partial _t \theta + \boldsymbol {v} \boldsymbol{\cdot }{\boldsymbol{\nabla }} \theta = \nu {\nabla} ^2 \theta - \mu {\nabla} ^{-2}\theta + f\kern-0.5pt. \end{equation}
The incompressibility condition
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {v}=0$
is enforced through the stream function
$\psi (\boldsymbol x,t)$
which defines the velocity field as
$\boldsymbol {v}({\boldsymbol x}, t) = (-\boldsymbol{\nabla} _y \psi , \boldsymbol{\nabla} _x \psi )$
. The relation between the buoyancy field and the stream function is given by
$\psi ={\nabla} ^{-1} \theta$
or, in Fourier space,
$\hat \psi =\hat \theta /k$
(where
$k\equiv |\boldsymbol k|$
) and consequently the velocity field can be written in terms of buoyancy field as
\begin{equation} \hat {\boldsymbol {v}}({\boldsymbol k}) = \left (-\frac {i k_y}{k}, \frac {i k_x}{k}\right ) \hat {\theta }({\boldsymbol k}) \end{equation}
from which one observes that
$\theta$
has the dimension of a velocity.
In the absence of forcing and dissipation (
$f = 0$
,
$\nu = 0$
,
$\mu = 0$
), the SQG equation (2.1) conserves two quadratic quantities, the vertically integrated energy (VIE),
\begin{equation} V = \tfrac {1}{2} \langle \psi \theta \rangle, \end{equation}
and the surface potential energy (SPE),
\begin{equation} E = \tfrac {1}{2} \langle \theta ^2 \rangle, \end{equation}
where angle brackets stand for spatial average.
The dissipative and forcing terms in the SQG equation represent the effects of scales not resolved by the model, their specific form being somewhat arbitrary (Smith et al. Reference Smith, Boccaletti, Henning, Marinov, Tam, Held and Vallis2002; Lapeyre & Klein Reference Lapeyre and Klein2006; Burgess, Scott & Shepherd Reference Burgess, Scott and Shepherd2015). In (2.1) we introduce a diffusivity
$\nu$
and a large-scale damping
$\mu$
chosen to be active at, respectively, small and large scales only, while the forcing is active on a narrow range of scales around
$\ell _f$
.
Under these conditions, one expects that a turbulent flow develops with a double cascade phenomenology (Blumen Reference Blumen1978; Pierrehumbert et al. Reference Pierrehumbert, Held and Swanson1994). Within this scenario, SPE is primarily transferred from the forcing scale to smaller ones (
$\ell \lt \ell _f$
), producing the direct cascade which is eventually dissipated by viscosity at the diffusive scale
$\ell _\nu$
. Meanwhile, VIE undergoes an inverse cascade, transferring energy to scales larger than the forcing (
$\ell \gt \ell _f$
) until it is dissipated at the friction scale
$\ell _\mu$
. In the statistically stationary state, the SPE and VIE balances are given by
\begin{align} \varepsilon _I &= \varepsilon _\nu + \varepsilon _\mu , \end{align}
\begin{align} \eta _I &= \eta _\nu + \eta _\mu , \end{align}
where
$\varepsilon _I = \langle \theta f\rangle$
and
$\eta _I=\langle \psi f\rangle$
are the input rates,
$\varepsilon _\nu = \nu \langle |{\boldsymbol{\nabla }} \theta |^2\rangle$
and
$\eta _\nu =\nu \langle \boldsymbol{\nabla }\psi \boldsymbol{\cdot }\boldsymbol{\nabla }\theta \rangle$
are the small-scale dissipation rates while
$\varepsilon _\mu = \mu \langle \theta {\nabla} ^{-2} \theta \rangle$
and
$\eta _\mu =\mu \langle \psi {\nabla} ^{-2} \theta \rangle$
are the large-scale dissipation rates of SPE and VIE, respectively.
Under the assumptions of statistical homogeneity and isotropy, Blumen (Reference Blumen1978) first predicted the power-law behaviour of the energy spectrum
$E(k)\equiv \langle |\hat \theta ({\boldsymbol k})|^2\rangle /2$
for sufficiently large scale separations
$\ell _\nu \ll \ell _f\ll \ell _\mu$
. In this case, the spectral energy densities follow
\begin{align}&\, E(k) \simeq \eta _\mu ^{2/3} k^{-1}, \quad 1/\ell _\mu \ll k \ll 1/\ell _f , \end{align}
\begin{align}& E(k) \simeq \varepsilon _\nu ^{2/3} k^{-5/3}, \quad 1/\ell _f \ll k \ll 1/\ell _\nu \end{align}
and the diffusive and friction scales are determined on dimensional grounds as
\begin{equation} \ell _\nu \equiv \left (\frac {\nu ^3}{\varepsilon _\nu } \right )^{1/4}, \quad \ell _\mu \equiv \left (\frac {\eta _\mu }{\mu ^3} \right )^{1/9}\kern-1.3pt. \end{equation}
The ratio between the dissipative and forcing scales, i.e. the extension of the inertial range, defines the Reynolds numbers associated with the flow. For the direct cascade of SPE, in analogy to 3-D NS turbulence, we define the Reynolds number as
\begin{equation} \textit{Re} \equiv \frac {\varepsilon _I^{1/3} \ell _f^{4/3}}{\nu } \simeq \left (\frac {\ell _f}{\ell _\nu } \right )^{4/3}\kern-1.3pt. \end{equation}
Note that in (2.10)
$\textit{Re}$
is based on
$\varepsilon _I$
and it is therefore defined a priori. Alternatively, we could use the VIE dissipation
$\varepsilon _{\nu }$
resulting in a slightly smaller value of
$\textit{Re}$
, as discussed below.
The scaling laws in (2.7)–(2.8) can also be obtained from the analogous for SQG of the exact four-fifths law of turbulence (Valadão et al. Reference Valadão, Ceccotti, Boffetta and Musacchio2024; Valade, Bec & Thalabard Reference Valade, Bec and Thalabard2025) and by assuming self-similarity of the statistics (Frisch Reference Frisch1995). Following this approach, one expects velocity and surface buoyancy to have the same scaling exponent
$1/3$
, which corresponds to that of the velocity field in 3-D NS turbulence. This observation allows us to adapt to the SQG model the dimensional arguments developed by Ruelle for the predictability of 3-D NS turbulence (Ruelle Reference Ruelle1979). According to the latter, the Lyapunov exponent
$\lambda$
characterising the exponential growth of an infinitesimal perturbation of a solution of (2.1) should be proportional to the inverse of the smallest dynamical time, i.e. the Kolmogorov time
$\tau _{\nu }\equiv (\nu /\varepsilon _{\nu })^{1/2}$
:
\begin{equation} \lambda \simeq \frac {1}{\tau _\nu }\simeq \frac {1}{\tau _f} \textit{Re}^{1/2}. \end{equation}
In the following section, we numerically investigate the scaling properties of the surface buoyancy field and the prediction (2.11).
3. Numerical simulations and results
We explore the statistical properties and the Eulerian predictability of the direct cascade in SQG at different Reynolds numbers by numerically integrating (2.1) at high resolution with a pseudo-spectral, GPU-accelerated code. Simulations are performed in a square domain of size
$L_x = L_y = 2\pi$
with periodic boundary conditions, using a regular grid with resolution
$N\times N$
ranging from
$N = 1024$
to
$N=16\,384$
. Simulations cover more than two decades in the diffusion coefficient, corresponding to a Reynolds number (2.10) which varies from
$\textit{Re}=600$
to
$\textit{Re}=1\,59\,000$
, while the large-scale dissipation coefficient
$\mu$
is fixed. For all runs, the system is driven by a constant-amplitude forcing with random phases, active within a narrow circular shell in wavevector space centred on
$k_f=3.5$
and with a small width
$\Delta k=0.5$
. This forcing provides constant SPE and VIE injection rates
$\varepsilon _I$
and
$\eta _I$
, respectively, with
$\varepsilon _I \approx \eta _I k_f$
since
$\Delta k \ll k_f$
. Specific details of the GPU code performances can be found in Valadão et al. (Reference Valadão, Boffetta, De Lillo, Musacchio and Crialesi-Esposito2025).
The most relevant parameters for the simulations are listed in table 1. All the simulations are performed in statistically stationary states, including the subset of simulations for computing the Lyapunov exponent. We also performed a careful study on the sensitivity of the following results to the maximum resolved wavenumber
$k_{\textit{max }}\ell _\nu$
by increasing resolution at fixed
$\textit{Re}$
. We found independence of the results on the resolution for
$k_{\textit{max }}\ell _\nu \gtrsim 1.5$
.
Table 1. Relevant parameters of the simulation: Reynolds number, diffusive scale
$\ell _\nu =\nu ^{3/4}\varepsilon _\nu ^{-1/4}$
, diffusive time
$\tau _\nu =\sqrt {\nu /\varepsilon _\nu }$
, total length of the Lyapunov simulations
$T_{\textit{tot}}$
(
$9$
independent realisations of length
$287 \tau _{f}$
for runs
$D$
). Common parameters for all simulations: forcing wavenumber
$k_f=3.5$
and width
$\varDelta k_f=0.5$
, surface potential energy input
$\varepsilon _I=24$
, friction coefficient
$\mu =1.0$
, characteristic time at the forcing scale
$\tau _f=\varepsilon _I^{-1/3}\ell _f^{2/3}=0.51$
, maximum resolved wavenumber
$k_{\textit{max }}=N/3$
(2/3 dealiasing rule).
3.1. Statistics of the direct cascade
For large separations between the forcing and the dissipative scales, one expects that almost all the SPE is transferred and dissipated at small scales, while VIE is dissipated at large scales. This is the essence of the argument developed by Fjørtoft for 2-D turbulence (Fjørtoft Reference Fjørtoft1953) and verified by numerical simulations of 2-D NS double cascade (Boffetta & Musacchio Reference Boffetta and Musacchio2010). In the present case, using (2.5) and (2.6) and the scaling relation
$\eta _{\ell } \simeq \ell \varepsilon _{\ell }$
this argument gives
\begin{equation} \frac {\eta _\nu }{\eta _\mu }= \left (\frac {\ell _\mu -\ell _f}{\ell _f-\ell _\nu }\right )\frac {\ell _\nu }{\ell _f} \simeq \textit{Re}^{-3/4}, \end{equation}
where we have used (2.10) to express
$\ell _{\nu }/\ell _f$
as a function of
$\textit{Re}$
.
Figure 1. The SPE relative dissipation
$\varepsilon _\nu / \varepsilon _I$
(filled circles) and VIE relative dissipation
$\eta _\nu / \eta _I$
(open circles), both as functions of
$\textit{Re}$
. The inset shows
$\eta _\nu /\eta _I$
as a function of
$\textit{Re}$
on a log–log scale.
Figure 1 shows the fraction of small-scale dissipation of the inviscid invariants as a function of the Reynolds number of the flow in statistically stationary conditions. Indeed, we find that the small-scale relative dissipation of VIE vanishes in the limit of large
$\textit{Re}$
following the prediction (3.1). On the contrary, since we do not resolve the inverse cascade and
$\ell _{\mu } \sim \ell _f$
, there remains a constant large-scale SPE dissipation
$\varepsilon _\mu$
for large
$\textit{Re}$
. Therefore, the direct cascade transfers only a fraction of the total injected energy equivalent to
$\varepsilon _\nu \approx 0.45\varepsilon _I$
for
$\textit{Re}\gtrsim 10^4$
.
Figure 2. Normalised SPE direct cascade fluxes
$\varPi (k)$
of all runs of table 1. The dashed line represents
$\varepsilon _\nu / \varepsilon _I=0.45$
.
Stationary fluxes of SPE in Fourier space are presented in figure 2 for different Reynolds numbers. At moderate
$\textit{Re} \lesssim 3000$
the fluxes for
$k\gt k_f$
decay quickly as a consequence of the viscous dissipation. For large
$\textit{Re}$
, however, a plateau of constant flux emerges at a level corresponding to the viscous dissipation rate
$\varepsilon _\nu$
. We emphasise that SQG turbulence exhibits large fluctuations in the flux of the direct cascade, as studied in detail in Valadão et al. (Reference Valadão, Ceccotti, Boffetta and Musacchio2024). These fluctuations arise from the interplay between the accumulation of energy in large-scale structures and intense dissipative events triggered by the formation of filamentary shocks that transfer energy from large to small scales over short time intervals. Thus, very long integrations are necessary to observe the convergence to the constant flux plateau of figure 2.
Figure 3. Time-averaged spectra
$E(k)$
for all simulations. Colour coding is the same as that of figure 2. Inset:
$C_K$
as a function of the Reynolds number.
The time-averaged spectra
$E(k)$
of SPE are shown in figure 3 for all the runs in table 1. All spectra exhibit power-law behaviour,
$E(k) \propto k^{-\beta }$
, in an intermediate range of scales, which becomes wider as
$\textit{Re}$
increases. We observe that at moderate
$\textit{Re}$
, when a scaling range is already clearly observable, the scaling exponent
$\beta$
deviates significantly from the dimensional prediction
$5/3$
, a feature already reported by previous investigations at comparable Reynolds numbers (Pierrehumbert et al. Reference Pierrehumbert, Held and Swanson1994; Ohkitani & Yamada Reference Ohkitani and Yamada1997; Sukhatme & Pierrehumbert Reference Sukhatme and Pierrehumbert2002). Nonetheless, we find that when
$\textit{Re}$
is sufficiently large,
$\textit{Re} \gtrsim 10^4$
, the scaling of the dimensional prediction (2.8) is closely recovered.
To quantify this important result, we measured the correction
$\xi$
to the dimensional scaling exponent by fitting the intermediate range of the spectra with
\begin{equation} E_K(k) = C_K \varepsilon _\nu ^{2/3} k^{-5/3} \left (\frac {k}{k_f}\right )^{-\xi}\kern-3pt, \end{equation}
where
$\xi$
and
$C_K$
are the fitting parameters. In order to estimate the robustness of the fit, we adopted the following procedure. For each run, we fit the data with (3.2) in a range of wavenumbers
$k \in [k_0, k_1]$
with varying
$k_0\in [3k_f, 5k_f]$
and
$k_1\in [8k_f, 10k_f]$
. This produces a set of parameters
$\xi$
and
$C_K$
for each run, from which we compute the mean using twice the standard deviation as an estimation of the error.
In figure 4, we plot the dependence of
$\xi$
on
$\textit{Re}$
, together with the spectra of figure 3 compensated with the expression (3.2). It is evident that, while for
$\textit{Re} \lesssim 10^4$
, the exponent correction
$\xi$
depends on
$\textit{Re}$
(approximately as
$\textit{Re}^{-3/4}$
), for larger values of
$\textit{Re}$
, the correction decreases much faster and becomes smaller than
$5\, \%$
for
$\textit{Re}\gt 2 \times 10^4$
. In this limit, we also observe the convergence of the dimensionless constant to
$C_K = 5.05 \pm 0.11$
(see figure 3).
We remark that this behaviour, which suggests the existence of a minimum Reynolds number for the recovery of dimensional scaling, is very different from what was observed in the direct cascade of 3-D NS turbulence, where Kolmogorov scaling is observed as soon as the spectrum displays a power-law behaviour.
3.2. Predictability of the direct cascade
In this section, we investigate the predictability of the cascade by computing how two solutions
$\theta (\boldsymbol{x},t)$
and
$\theta '(\boldsymbol{x},t)$
separate in time on average. We consider infinitesimally close solutions so that the average separation rate is given by the maximal Lyapunov exponent of the flow.
Starting from a solution
$\theta (\boldsymbol{x},t)$
of (2.1) in a statistically stationary state, we generate a perturbed solution as
$\theta '(\boldsymbol{x},t)=\theta (\boldsymbol{x},t)+2\sqrt {\varDelta } W({\boldsymbol x})$
, where
$W({\boldsymbol x})$
is a Gaussian random white noise with zero mean and unit variance while
$\varDelta$
is a small parameter. The SPE error
$E_{\varDelta }$
is defined, for any time, as
\begin{equation} E_{\varDelta }(t)=\tfrac {1}{2} \langle \delta \theta ({\boldsymbol x},t) ^2 \rangle ,\end{equation}
where the difference field is
$\delta \theta =(\theta '-\theta )/\sqrt {2}$
and the normalisation coefficient
$1/\sqrt {2}$
ensures that
$E_{\varDelta }=E$
for two completely uncorrelated fields. At initial time, by definition, we have
$E_{\varDelta }(t)=\varDelta$
. We measure the FTLE by computing the growth rate of the error
\begin{equation} \gamma _{\tau }(t)=\frac {1}{2\tau } \ln \left (\frac {E_\varDelta (t+\tau )}{\varDelta } \right ) \end{equation}
and then by rescaling the perturbed field to the initial SPE error
\begin{equation} \theta ' \leftarrow \theta - \sqrt {\frac {\varDelta }{E_\varDelta }} (\theta - \theta '). \end{equation}
By repeating the steps (3.4) and (3.5) over many time intervals of the same length
$\tau$
, we obtain a distribution of the FTLE along the trajectory. The rescaling procedure (3.5) ensures the permanence of the perturbation in the exponential growth regime when
$\varDelta$
and
$\tau$
are sufficiently small (Vulpiani, Cecconi & Cencini Reference Vulpiani, Cecconi and Cencini2009).
From the definition (3.4) one can compute the FTLE for any time multiple of
$\tau$
,
$T=n \tau$
, simply by averaging
\begin{equation} \gamma _{T}(t) = \frac {1}{n} \sum _{k=1}^{n} \gamma _{\tau }(t+k \tau ) \end{equation}
and the Lyapunov exponent is given by the average of the FTLE over a very long trajectory (and becomes independent of the initial condition):
\begin{equation} \lambda = \lim _{T \to \infty } \gamma _{T}(t). \end{equation}
In general, the distribution of the FTLE around the Lyapunov exponent, for sufficiently large
$T$
, follows the large deviation principle (Vulpiani et al. Reference Vulpiani, Cecconi, Cencini, Puglisi and Vergni2014) which states that
\begin{equation} \rho (\gamma _T) = \frac {1}{N_T} {\rm e}^{-T C(\gamma _T) }, \end{equation}
where
$N_T$
is a normalising factor and
$C(\gamma _T)$
is the Cramér function, independent of
$T$
which, in general, vanishes at
$\gamma _T=\lambda$
and is positive for
$\gamma _T \neq \lambda$
(Boffetta et al. Reference Boffetta, Cencini, Falcioni and Vulpiani2002). For not too large fluctuations, the Cramér function can be approximated by a quadratic form:
\begin{equation} C(\gamma _T)\approx \frac {(\gamma _T-\lambda )^2}{2\varOmega } , \end{equation}
where
$\varOmega$
, proportional to the variance of the distribution
$\rho (\gamma _T)$
, is obtained from
\begin{equation} \varOmega = \lim _{T \to \infty }T\langle (\gamma _T-\lambda )^2\rangle _{\gamma }, \end{equation}
where the angle brackets represent the average over the distribution (3.8) while
$\varOmega$
is expected to be independent of
$T$
in the limit of large
$T$
.
We computed the FTLE for simulations at different Reynolds numbers corresponding to runs
$B_2$
,
$B_4$
,
$C_2$
,
$C_3$
,
$C_4$
,
$C_5$
,
$D_1$
,
$D_2$
and
$D_3$
of table 1. For all runs, we excluded from the statistics the initial transient during which the perturbation aligns with the most unstable direction of the system. For the three cases at the highest Reynolds numbers, we compensated for the increased computational cost by averaging over nine independent, shorter simulations run in parallel, each of which with different realisations of the forcing
$f(\boldsymbol{x})$
and initial perturbation noise
$W(\boldsymbol{x})$
. All rescaling times
$\tau$
were kept around the Kolmogorov time scale of the simulation
$\tau _\nu$
, more precisely, between
$\tau /\tau _\nu \in [0.5,0.8]$
depending on the run.
Figure 5 shows a representative realisation of the field
$\theta (\boldsymbol{x})$
along with its corresponding perturbation field
$\delta \theta (\boldsymbol{x})$
. The errors accumulate predominantly in filamentary zones between coherent structures. These regions are dominated by small-scale structures formed by energy transfer from larger scales (Pierrehumbert et al. Reference Pierrehumbert, Held and Swanson1994). Since such structures appear intermittently in time, the convergence of the FTLE statistics requires very long simulations, as discussed below.
Figure 5. Fields (a)
$\theta (\boldsymbol x)$
and (b)
$\delta \theta (\boldsymbol x)$
of run
$D_3$
. The colour code of (a) is based on
$\theta (\boldsymbol x)/\sup _{\boldsymbol x}|\theta (\boldsymbol x)|$
. The colour code of (b) is based on a log scale
$(=\log _{10}(|\delta \theta (\boldsymbol x)|/\sup _{\boldsymbol x}|\delta \theta (\boldsymbol x)|))$
to facilitate visualisation.
Figure 6. Energy spectra
$E(k)$
(full lines) and error spectra
$E_\varDelta (k)$
(dotted lines) for the different runs as functions of
$k\ell _\nu$
. The error spectra are computed from the difference field
$\delta \theta$
before the rescaling (3.5) to the initial error.
The distribution of the error field among scales, right before the rescaling, is shown in figure 6 where we plot, together with the energy spectra
$E(k)$
, the error spectra defined in a similar way as
$E_{\varDelta }(k)=\langle |\delta \hat {\theta }(\boldsymbol {k})|^2 \rangle /2$
. We observe that for all simulations, the error is concentrated at small scales
$k \ell _{\nu } \simeq 0.1$
close to the dissipative range. We remark that the relative magnitude
$E_{\varDelta }(k)/E(k)$
is very small, ensuring that the perturbation remains in the linear regime.
Figure 7 presents the FTLE for the different runs as a function of the average time
$T$
. We see that in all the cases, the average FTLE converges, after a long transient and for
$T \gtrsim 200\tau _f$
, to the asymptotic value, which represents the Lyapunov exponent of the flow.
Figure 7. Convergence of the FTLE as a function of the average time. Light to dark colours represent runs at increasing Reynolds numbers:
$B_2$
,
$B_4$
,
$C_2$
,
$C_3$
,
$C_4$
,
$C_5$
,
$D_1$
,
$D_2$
and
$D_3$
.
From figure 7, it is evident that the Lyapunov exponent increases with the Reynolds number. As discussed in § 2, we expect that
$\lambda$
follows the Ruelle scaling (2.11). The latter relies on the spectrum having a Kolmogorov-like power-law behaviour. As shown in figure 4, the scaling exponent
$5/3$
is recovered only for
$\textit{Re} \gt 2 \times 10^4$
. Ruelle’s prediction is therefore expected to hold only for large
$\textit{Re}$
.
In figure 8, we plot the Lyapunov exponents of our simulations as functions of
$\textit{Re}$
. We find that
$\lambda$
grows with
$\textit{Re}$
faster than that predicted by (2.11) and the best fit gives
$\lambda \tau _f \simeq \textit{Re}^{0.7}$
or, equivalently,
$\lambda \tau _{\nu } \simeq \textit{Re}^{0.2}$
. Remarkably, the scaling persists also at
$\textit{Re}\lt 10^4$
, where the spectrum displays a significant correction to the Kolmogorov exponent
$5/3$
(see figure 4).
In figure 8 we also plot the scaled variance
$\varOmega$
as a function of
$\textit{Re}$
, which displays a scaling law compatible with that of the Lyapunov exponent
$\varOmega \tau _f \simeq \textit{Re}^{0.7}$
. This indicates that, in the range of
$\textit{Re}$
investigated here, the ratio
$\varOmega /\lambda$
is approximately constant (
$0.24\pm 0.02$
) and that the central part of the Cramér function has a self-similar evolution with
$\textit{Re}$
.
Figure 8. Reynolds scaling of the mean FTLE
$\langle \gamma _T\rangle =\lambda$
. Non-dimensionalisation is made with (a)
$\tau _f$
and (b)
$\tau _\nu$
.
A qualitatively similar behaviour has been observed for the Lyapunov exponent of 3-D turbulence (Boffetta & Musacchio Reference Boffetta and Musacchio2017; Mohan et al. Reference Mohan, Fitzsimmons and Moser2017; Berera & Ho Reference Berera and Ho2018; Ge, Rolland & Vassilicos Reference Ge, Rolland and Vassilicos2023) with a correction to the dimensional scaling (2.11) slightly smaller than in the present case,
$\lambda \propto \textit{Re}^{0.64}$
. We remark that the origin of this correction in 3-D turbulence is still unclear since it cannot be simply attributed to intermittency. Indeed, the multifractal extension of the dimensional prediction (2.11) predicts, for 3-D turbulence, an exponent smaller than
$1/2$
(Aurell et al. Reference Aurell, Boffetta, Crisanti, Paladin and Vulpiani1996).
Figure 9. (a) Cramér function for run
$C_4$
and different times. (b) Cramér function for different
$\textit{Re}$
at fixed
$T\lambda =7$
. The dashed line represents the quadratic form
$(x-1)^2/(2 \varOmega /\lambda )$
with
$\varOmega /\lambda =0.24$
.
In figure 9(a), we show the Cramér function
$C(\gamma _{T})$
obtained from the distribution of the FTLE (3.8) as
$C(\gamma _{T}) = - ({1}/{T}) \log (\rho (\gamma _T)/\rho (\lambda ))$
at different times
$T$
at
$\textit{Re} = 21\,200$
(run
$C_4$
). We notice that the convergence towards the asymptotic behaviour of the left tail (
$\gamma _T \lt \lambda$
) is much faster than that of the right tail (
$\gamma _T \gt \lambda$
). At
$T \gtrsim 2/\lambda$
, the left tail is already converged, while the right tail achieves convergence only for
$T \gtrsim 7/\lambda$
. Figure 9(b) shows the comparison of the asymptotic behaviour of the Cramér functions
$C(\gamma _{T})$
computed at
$T\lambda = 7$
for different
$\textit{Re}$
in the range
$[10\,600, 63\,500]$
.
Since
$\lambda$
and
$\varOmega$
scale with
$\textit{Re}$
with the same exponent, we rescale both
$\gamma _T$
and
$C(\gamma _T)$
by
$\lambda$
which, according to the quadratic approximation (3.9), predicts the collapse to the function
$(x-1)^2/(2 \varOmega /\lambda )$
.
Figure 9 shows that the quadratic approximation (3.9) with
$\varOmega /\lambda = 0.24$
describes well the core of
$C(\gamma _{T})$
(for
$\gamma \simeq \lambda$
), while we observe significant deviations for
$\gamma _{T}\gt \lambda$
. This means that the situations in which the dynamics of the system is strongly unpredictable (i.e. with
$\gamma _{T} \gg \lambda$
) are more frequent than what is expected by Gaussian statistics. Nonetheless, it is remarkable that the probability of these extreme deviations seems to be independent of the Reynolds numbers, as shown by the collapse of the right tail of
$C(\gamma _{T})$
for different
$\textit{Re}$
.
4. Conclusions
In this study, we investigated the statistical properties and predictability of turbulence in SQG model using high-resolution direct numerical simulations across a wide range of Reynolds numbers. Our analysis focused on two central aspects: the scaling behaviour of the energy spectrum in the direct cascade of SPE and the chaotic dynamics characterised by FTLEs.
Our results indicate that, for
$\textit{Re}\gtrsim 2\times 10^4$
, the energy spectrum approaches the Kolmogorov-like scaling
$E(k) \propto k^{-5/3}$
. This observation, together with the convergence in Reynolds of the prefactor
$C_K$
, suggests that SQG turbulence exhibits a well-defined inertial range – similar to that of 3-D NS turbulence – but it does so only at very large Reynolds numbers. This regime was not observed by earlier studies at moderate Reynolds numbers.
In terms of predictability, we showed that the Lyapunov exponent scales anomalously with the Reynolds number as
$\lambda \propto \textit{Re}^{0.7}$
, exceeding the dimensional prediction
$\lambda \propto \textit{Re}^{1/2}$
. The presence of an anomalous scaling is reminiscent of similar behaviour observed in 3-D NS turbulence (Boffetta & Musacchio Reference Boffetta and Musacchio2017; Mohan et al. Reference Mohan, Fitzsimmons and Moser2017; Berera & Ho Reference Berera and Ho2018; Ge et al. Reference Ge, Rolland and Vassilicos2023). In particular, it has been observed that
$\lambda \propto \textit{Re}^{0.64}$
(Boffetta & Musacchio Reference Boffetta and Musacchio2017). This is in contrast with the scaling of the Lyapunov exponent of Lagrangian trajectories
$\lambda _L$
, which has been found to be in agreement with the Ruelle dimensional prediction
$\lambda _L \propto \tau _\nu$
(Bec et al. Reference Bec, Biferale, Boffetta, Cencini, Musacchio and Toschi2006) up to minor corrections which can be estimated from multifractal arguments.
The difference between Eulerian and Lagrangian predictability in NS turbulence can be ascribed to the fact that the Lagrangian Lyapunov exponent is determined by the viscous time scale
$\tau _\nu$
, while the Eulerian predictability can be influenced also by the sweeping time scale
$\tau _E = \ell _\nu /U$
, where
$U$
is the root-mean-square velocity (Ge et al. Reference Ge, Rolland and Vassilicos2023). Considering that
$\tau _E \ll \tau _\nu$
and that
$\tau _\nu /\tau _E \simeq \textit{Re}^{1/4}$
, the dependence of
$\lambda$
on
$\textit{Re}$
can be rewritten as a combination of two time scales:
\begin{align} \lambda = \tau _\nu ^{-1} \textit{Re}^\xi = \tau _\nu ^{-1+4\xi }\tau _E^{-4\xi }.\\[-28pt]\nonumber \end{align}
The Ruelle scaling is recovered for
$\xi =0$
, while for
$\xi =1/4$
the predictability is determined by the sweeping time. Our results in SQG give
$\xi ^{\textit{SQG}}=0.2$
which is larger than the value observed in 3-D NS,
$\xi ^{NS} = 0.14$
, suggesting that the sweeping effect has a stronger influence on the Eulerian predictability in SQG than in NS.
It is worth noting that the dimensional arguments discussed above do not consider the effects of intermittency corrections. It would be interesting to develop a multifractal-like approach for the direct cascade of SQG to check how the observed intermittency (Valade et al. Reference Valade, Bec and Thalabard2025) affects the dimensional scaling of the Lyapunov exponent.
Beyond the average growth rate of infinitesimal perturbations, we investigated the statistical properties of FTLEs. Notably, both the variance and the shape of the associated Cramér function exhibit self-similar behaviour across Reynolds numbers when appropriately rescaled. The ratio
$\varOmega / \lambda$
remains approximately constant across
$\textit{Re}$
, indicating a form of universality in the core of the FTLE distribution.
Although we have studied the predictability problem from the point of view of the exponential growth of infinitesimal perturbations, it would be very interesting to investigate in detail the complementary regime of large errors and the statistics of finite-size Lyapunov exponents (Boffetta & Musacchio Reference Boffetta and Musacchio2017). Recent results have been obtained in the case of a decaying SPE cascade in SQG (Valade et al. Reference Valade, Thalabard and Bec2024), where the authors were able to connect the hyperdiffusive behaviour of Lagrangian fluid parcels with the anomalous diffusion of the system.
Acknowledgements.
We acknowledge HPC CINECA for computing resources within the INFN-CINECA Grants INFN24-FieldTurb and INFN25-FieldTurb.
Funding.
This work was supported by Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union - NextGenerationEU - and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2 within the activities of Spoke 3 (Astrophysics and Cosmos Observations).
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 request.
Author contributions.
V.J.V. performed the simulations and analysed the data. All authors contributed equally in the discussions to reach conclusions and in writing the paper.
Open access
