1. Introduction
Turbulent flows are the fundamental basis of many applications in engineering (Pope Reference Pope2001), geophysics (Thorpe Reference Thorpe2007; Gill Reference Gill2016) and astrophysics (Biskamp Reference Biskamp2003) among others. Despite the considerable improvement in computational power, performing direct numerical simulations (DNS) that resolve all active scales in realistic turbulent fluids remains out of reach for most industrial and natural flows (Moin & Mahesh Reference Moin and Mahesh1998; Pope Reference Pope2001). This is related to the fact that the number of degrees of freedom grows fast with the size of the system. To circumvent this prohibitive cost, large-eddy simulations (LES) are widely adopted in engineering and planetary atmospheric models. The key idea is to take the evolution only down to a specified filtering scale, with the influence of smaller, unresolved subfilter motions taken into account through explicit modelling (Smagorinsky Reference Smagorinsky1963; Deardorff Reference Deardorff1970; Germano Reference Germano1992; Meneveau & Katz Reference Meneveau and Katz2000; Sagaut Reference Sagaut2005). Because the equations remain unclosed, due to the nonlinear Navier–Stokes operator coupling all scales, empirical models grounded in turbulence theory are employed. Many of these models account for the influence of the small-scale motions through the subgrid-scale stress tensor
${\tau }_{\textit{ij}}$
, which captures the action of unresolved scales on the resolved ones. Typically,
$\tau _{\textit{ij}}$
is expressed as a function of the large-scale velocity gradients, consistent with the symmetries of the flow, in order to provide a deterministic model closure for turbulence (Meneveau & Katz Reference Meneveau and Katz2000). The models by Smagorinsky (Reference Smagorinsky1963) and Clark, Ferziger & Reynolds (Reference Clark, Ferziger and Reynolds1979) are some of the most well-known ones. Being quasi-empirical, these models require extensive testing and validation through experiments and simulations (Piomelli et al. Reference Piomelli, Cabot, Moin and Lee1991; Domaradzki, Liu & Brachet Reference Domaradzki, Liu and Brachet1993; Härtel et al. Reference Härtel, Kleiser, Unger and Friedrich1994; Liu, Meneveau & Katz Reference Liu, Meneveau and Katz1994). These models have been found to reasonably reproduce even some intermittency measures Linkmann, Buzzicotti & Biferale (Reference Linkmann, Buzzicotti and Biferale2024) and statistical symmetries Magacho et al. (Reference Magacho, Thalabard, Buzzicotti, Bonaccorso, Biferale and Mailybaev2024), however, more severe tests reveal often significant discrepancies between modelled and measured stress components (Meneveau & Katz Reference Meneveau and Katz2000; Moser, Haering & Yalla Reference Moser, Haering and Yalla2021). In particular, looking at the full turbulent cascade process, the eddy-viscosity approximation may lead to oversimplifying the picture with large errors in fluctuations. In previous works (Borue & Orszag Reference Borue and Orszag1998; Tsinober Reference Tsinober2009; Alexakis & Biferale Reference Alexakis and Biferale2018; Alexakis & Chibbaro Reference Alexakis and Chibbaro2020), the local energy flux transferred from unfiltered to filtered scales was investigated. The mean was found to be consistent with the predictions of the simple Smagorinsky model, however, significant deviations were observed. It is possible that, since these deviations originate from the small, unresolved scales, their impact on the resolved scales could be captured using stochastic modelling. While stochastic subgrid models have been proposed in previous studies (Mason & Thomson Reference Mason and Thomson1992; Gicquel et al. Reference Gicquel, Givi, Jaberi and Pope2002; Carati, Winckelmans & Jeanmart Reference Carati, Winckelmans and Jeanmart2007; Marstorp, Brethouwer & Johansson Reference Marstorp, Brethouwer and Johansson2007), the statistical properties of the associated noise terms were typically prescribed in an ad hoc fashion, and notably are not meant to capture extreme events. The need for an accurate stochastic description of the small scale effects is emphasised by recent results proposing that the Navier–Stokes equations can develop randomness at large scales in a finite time (Eyink & Bandak Reference Eyink and Bandak2020; Thalabard, Bec & Mailybaev Reference Thalabard, Bec and Mailybaev2020; Bandak et al. Reference Bandak, Mailybaev, Eyink and Goldenfeld2024).
In the present work, we aim to quantify the fluctuations around the mean subgrid-scale energy flux and compare these with the predictions of classical models. To this end, we employ data from DNS of the Navier–Stokes equations and follow two complementary approaches: (i) the analysis of several conditional means, as suggested in previous works (Buaria, Pumir & Bodenschatz Reference Buaria, Pumir and Bodenschatz2020; Johnson & Wilczek Reference Johnson and Wilczek2024); notably we shall analyse probabilities conditioned on
$\textit{QR}$
, the invariants of the velocity gradient tensor, which have not been investigated so far. (ii) The joint probability density function (PDF) of model predictions and ground truth, as in Alexakis & Chibbaro (Reference Alexakis and Chibbaro2020), to examine the distribution of subgrid-scale fluxes conditioned on a given model. The ultimate goal is to guide and constrain the development of improved stochastic subgrid-scale models for turbulent flow simulations. To this end, we restrict our analysis of fluctuations to the simplest, and most used deterministic closures, namely the Smagorinsky and Clark models. Although it should be possible in principle to refine the accuracy of such closures by including higher-order velocity gradients in the expansion of the subgrid-scale stress tensor (Carati, Winckelmans & Jeanmart Reference Carati, Winckelmans and Jeanmart2001; Eyink Reference Eyink2006; Chang, Yuan & Wang Reference Chang, Yuan and Wang2022), we deliberately adopt simple forms in order to emphasise the stochastic component, and to provide insights for further practical models. Moreover, with regard to possible future machine learning approaches, the key point is to provide a well-conditioned framework for learning. Whether more accurate deterministic closures are required to achieve this is beyond the scope of the present work.
2. Theoretical background
The stage is set by the incompressible Navier–Stokes equations, governing the evolution of
$u$
, divergence-free velocity field,
\begin{equation} \big [\partial _t-\nu \partial _{\textit{jj}}\big ]u_i+ \partial _{\!j}\big (u_{\!j} u_i \big )+\partial _i\! p=f_i. \end{equation}
Here,
$\nu$
is the kinematic viscosity and
$p$
is the pressure computed through incompressibility. An external forcing term is added on the right-hand side, acting at the largest scales of the system. This forcing is statistically homogeneous and injects energy into the flow at a mean rate
$\epsilon$
. The domain is a triply periodic cubic box of size
$2\pi\! L$
, which ensures the absence of boundary-layer effects. We focus on the high Reynolds number regime, with
$\textit{Re} \equiv \epsilon ^{1/3} L^{4/3}/\nu$
, so that a wide separation of scales emerges between the forcing scale
$L$
and the viscous dissipation scale
$\eta \equiv \nu ^{3/4} \epsilon ^{-1/4}$
.
2.1. Filtering approach
A typical starting point for turbulence modelling is to separate the components describing variations at large and small scales, although their dynamics remain irrevocably coupled because of the nonlinear term. Following Germano (Reference Germano1992) we work with the coarse graining at scale
$l$
of
$\boldsymbol{u}(\boldsymbol{x},t)$
, the solution of the three-dimensional incompressible Navier–Stokes problem (2.1), obtained by suppressing high-frequency fluctuations through local averaging
\begin{equation} \overline {\boldsymbol{u}}^l(\boldsymbol{x},t)=\int\! G^l(\boldsymbol{r})\boldsymbol{u}(\boldsymbol{x}-\boldsymbol{r},t){\rm d}\boldsymbol{r}, \end{equation}
where
$G^l(\boldsymbol{r})=l^{-3}G(\boldsymbol{r}/l)$
is a localised smooth filtering function such that
$\int G(\boldsymbol{r}){\rm d}\boldsymbol{r}=1$
. In particular, we use the Gaussian filter, given in Fourier space by
\begin{equation} \widehat {G}(\boldsymbol{k})=\text{exp}\!\left [-\frac {k^2}{2}\right ]\!. \end{equation}
The dynamics of the resulting large-scale field is not closed because of the nonlinearity
\begin{equation} \big [\partial _t-\nu \partial _{\textit{jj}}\big ]\overline {u_i}^{\,l}+\big [\delta _\textit{ki}-\partial _\textit{ki}{\nabla} ^{-2}\big ]\partial _{\!j}\big (\overline {u_{\!j}}^{\,l} \overline {u_k}^{\,l} \big )=\overline {f_i}^{\,l}-\partial _{\!j}\tau ^l(u_i,u_{\!j}). \end{equation}
On the left-hand side, the standard Navier–Stokes operator acts on the coarse-grained field, while the right-hand side depends on the velocity field at all scales through the subgrid-scale tensor
\begin{equation} \tau ^l(u_i,u_{\!j})=\overline {u_i u_{\!j}}^{\,l}-\overline {u_i}^{\,l}\overline {u_{\!j}}^{\,l}=:\tau ^l_{\textit{ij}}, \end{equation}
the stress exerted on scales larger than
$l$
by the chaotic fluctuations at smaller scales. This tensor depends on both the filtered velocity field and the underlying small-scale fluctuations, and it determines the Galilean-invariant local energy transfer from the large-scale to the small-scale component
\begin{equation} \varPi ^l(\boldsymbol{x},t)=-\tau ^l_{\textit{ij}}\partial _{\!j} \overline {u_i}^{\,l}. \end{equation}
We note that the value of local energy flux here depends on the choice of filter and technical differences between filters have been noted (Alexakis & Chibbaro Reference Alexakis and Chibbaro2020). An alternative (filter-independent) definition of the energy flux within the structure function framework has been explored (Yao et al. Reference Yao, Schnaubelt, Szalay, Zaki and Meneveau2024). Here, we will refrain from examining further filters or alternative definitions and we will limit ourselves to the Gaussian filter. A model of turbulence can be obtained via the closure of (2.4) by expressing the unknown right-hand side in terms of the resolved large-scale field. Traditionally, the subgrid-scale stress tensor
$\tau$
is modelled as a function of the first-order gradients of the filtered velocity. For example, the model by Clark et al. (Reference Clark, Ferziger and Reynolds1979)
\begin{equation} \tau ^l_{\textit{ij}}\sim l^2\partial _k\overline {u}_i^{\,l} \partial _k\overline {u}_j^{\,l}, \end{equation}
which may be obtained as a first term of an expansion of (2.5) in
$l$
(Borue & Orszag Reference Borue and Orszag1998; Carati et al. Reference Carati, Winckelmans and Jeanmart2001; Eyink Reference Eyink2006; Johnson Reference Johnson2021). Another example is the traditional eddy-viscosity model by Smagorinsky (Reference Smagorinsky1963)
\begin{equation} \tau ^l_{\textit{ij}}-\frac {\delta _{\textit{ij}}}{3}\textrm {Tr}[\tau ^l]\sim -l^2\sqrt {\overline {S}^{\,l}_\textit{kh}\overline {S}^{\,l}_\textit{hk}}\overline {S}^{\,l}_{\textit{ij}}, \end{equation}
where
$S$
is the symmetric part of the velocity gradient tensor, the strain rate. Such models lead to a deterministic evolution of the filtered field.
2.2. Velocity gradient invariants
The tensor
$\partial _i u_{\!j}=S_{\textit{ij}}+\varOmega _{\textit{ij}}$
(where
$S_{\textit{ij}}$
is the symmetric and
$\varOmega _{\textit{ij}}$
the antisymmetric part of
$\boldsymbol{\nabla \!}\boldsymbol{u}$
) encapsulates the spatial variations of the velocity field around a given point and offers key insights into the local flow topology, intermittency and the mechanisms driving the energy cascade. For isotropic turbulence,
$\tau$
should depend only on the properties of
$\boldsymbol{\nabla \!}\boldsymbol{u}$
and higher-order derivatives that are invariant under rotations (Carati et al. Reference Carati, Winckelmans and Jeanmart2001; Eyink Reference Eyink2006). In models, however, only information from the gradient tensor
$\boldsymbol{\nabla \!}\boldsymbol{u}$
is usually kept. In total, the gradient tensor is fully characterised by
$5$
scalar-invariant quantities, that are independent of the orientation of the coordinate system. This reflects the fact that the tensor
$\boldsymbol{\nabla \!}\boldsymbol{u}$
possesses five effective degrees of freedom, obtained by starting from its nine components and accounting for the incompressibility constraint and the assumption of isotropy (Meneveau Reference Meneveau2011; Johnson & Wilczek Reference Johnson and Wilczek2024). To provide a synthetic description of the local fluid configurations without referring directly to the tensor components, all
$5$
invariants are in principle needed. However, here, we limit ourselves to
$Q$
and
$R$
, the invariants under reference frame rotations and reflections that characterise the spectrum of
$\boldsymbol{\nabla \!}\boldsymbol{u}$
. In general, for a traceless real matrix
$M$
we will denote
\begin{equation} Q_M=-\frac {1}{2}\textrm {Tr}(M^2),\quad R_M=-\frac {1}{3}\textrm {Tr}(M^3). \end{equation}
From now on, we simply refer to
$Q$
and
$R$
as invariants of the velocity gradient
$\boldsymbol{\nabla \!}\boldsymbol{u}$
(Chong, Perry & Cantwell Reference Chong, Perry and Cantwell1990)
\begin{equation} Q=Q_{{\boldsymbol{\nabla \!}\boldsymbol{u}}}=-\frac {1}{2} \partial _i u_{\!j} \partial _{\!j} u_i, \quad R=R_{{\boldsymbol{\nabla \!}\boldsymbol{u}}}=-\frac {1}{3} \partial _i u_{\!j} \partial _{\!j} u_k \partial _k u_i. \end{equation}
Combining the multiscale description, the invariants of velocity gradient and strain rate are sufficient to express the energy fluxes predicted by the introduced models (2.7) and (2.8)
\begin{equation} \varPi ^l_{\textit{Clark}}=l^2\big [4R_{\overline {S}^{\,l}}-R_{\overline {\boldsymbol{\nabla \!}\boldsymbol{u}}^{\,l}}\big ],\quad \varPi ^l_{\textit{Smag.}}=l^2\big [-Q_{\overline {S}^{\,l}}\big ]^{3/2}. \end{equation}
A striking difference between the models is that Smagorinsky’s prediction, unlike Clark’s, is non-negative. It is important to note that knowledge of the two invariants of the sole velocity gradient tensor does not uniquely determine those of its symmetric part
$S$
, and therefore model predictions (2.11) are not uniquely defined in the
$(Q,R)$
space. Yet, to link the energy flux to the invariants
$Q,R$
has proven useful to get insights into the dynamical process (Johnson & Wilczek Reference Johnson and Wilczek2024), and also for the construction of LES models (Borue & Orszag Reference Borue and Orszag1998). Indeed,
$Q$
and
$R$
provide information on the local geometry of the flow: depending on their sign, the eigenvalue problem of
$\boldsymbol{\nabla \!}\boldsymbol{u}$
admits either two complex conjugate roots and one real root, positive or negative, or three real roots. Taking into account the incompressibility constraint, this leads to four distinct regions in the
$(Q,R)$
plane, where the upper-right quadrant is associated with vortex compression. For more details we refer readers to the literature (Chong et al. Reference Chong, Perry and Cantwell1990; Meneveau Reference Meneveau2011; Johnson & Wilczek Reference Johnson and Wilczek2024).
3. Results
The present study relies on a numerical solution of the Navier–Stokes equations (2.1). The pseudo-spectral code GHOST (Mininni et al. Reference Mininni, Rosenberg, Reddy and Pouquet2011) was employed, using a fourth-order Runge–Kutta scheme and adapting the time step according to the Courant–Friedrichs–Lewy condition. The code considers a periodic cubic box of size
$2\pi$
with resolution
$1024$
in each direction. Viscosity was set to
$\nu =5\boldsymbol{\times }10^{-4}$
. In order to fix the mean energy injection rate
$\epsilon$
, we imposed the forcing (as, for example, in Alexakis & Chibbaro Reference Alexakis and Chibbaro2020)
\begin{align} \widehat {\kern -2pt \boldsymbol{f}}(\boldsymbol{k},t)=\left [\frac {\epsilon }{\sum _{p\leqslant 2}\vert \widehat {\boldsymbol{u}}(\boldsymbol{p},t)\vert ^2}+i\varphi _{\boldsymbol{k}}\right ]\widehat {\boldsymbol{u}}(\boldsymbol{k},t), \end{align}
where the phases
$\varphi _{\boldsymbol{k}}$
are randomly chosen at initialisation and kept fixed in time. In particular, we put
$\epsilon =1$
. After a transient, a steady state was reached, statistically consistent with the classical isotropic turbulence. Figure 1 shows, on the left, the energy spectrum, which provides a Kolmogorov–Obukhov
$5/3$
spectrum over almost two decades. On the right, it displays the mean energy transfer
$\varPi$
across scale
$k$
, normalised by the mean energy injection rate. The actual inertial range of scales is identified by the region where the normalised transfer is close to
$1$
. The blue line represents the spatial and temporal average of (2.6) using the sharp spectral filter (i.e. the Fourier formulation), while the red triangles correspond to the estimation based on Gaussian filtering. We select the wavenumber
$16$
as a representative point for the inertial range.
Kolmogorov energy spectrum
$E(k)$
in (a) and fraction of the mean energy transfer across scales relative to the mean energy injection rate (b).
3.1. Clark vs DNS
3.1.1. The QR diagrams
To compare a priori predictions given by models of turbulence we propose to look at the conditional mean
$\langle \varPi \vert Q,R\rangle$
. This observable has been first used in a recent work (Vela-Martín & Jiménez Reference Vela-Martín and Jiménez2021), yet in a different context. In figure 2 we consider
$Q$
and
$R$
normalised by the mean dissipation at the corresponding scale, and showing the mean flux for a given
$(Q,R)$
with the colour scale: positive values in red indicate a mean energy transfer to smaller scales, while negative values in blue correspond to a net transfer toward larger scales. White contours indicate the isolines of the probability density in the
$(Q,R)$
space, corresponding to levels 10
$^{-4}$
, 10
$^{-5}$
, 10
$^{-6}$
and 10
$^{-7}$
, starting from the curve closest to the origin, which corresponds to the maximum of the distribution. We have chosen the observation window such that the probability isoline corresponding to 10
$^{-7}$
is largely included, and for events with probability below 10
$^{-8}$
we have not shown the mean energy flux (black regions) because of the lack of examples to compute the averages. The violet curve marks the isoline where the mean flux vanishes. The left panel presents the results obtained from (2.6), whereas the right panel shows the prediction of the Clark model (2.7), calculated using the complete DNS fields. First, it is noteworthy that, in the DNS, the backward cascade (backscatter) is confined to the upper-right quadrant, associated with vortex compression configurations. Then, the Clark model successfully captures the qualitative distribution of the mean energy flux across the
$(Q,R)$
diagram. However, quantitatively, it underestimates the right value. It is possible to see that a region with small backscatter flow is observable in the lower-left quadrant, as highlighted by the violet curve. We can discuss a little more the error made by the Clark model from a physical point of view. As discussed in previous works (Eyink Reference Eyink2006), the turbulent subgrid stress may be written as an infinite series of terms in which the Clark model represents the first term. In particular, for a Gaussian filter as the one used here, it may be written (Johnson Reference Johnson2021)
\begin{equation} \tau ^l_{\textit{ij}}\sim l^2\partial _k\overline {u}_i^{\,l} \partial _k\overline {u}_j^{\,l} +\int _0^{l^2}\! {\rm d}\theta \left (\overline {\overline {\boldsymbol{\nabla} _i \boldsymbol{u}_k}^{\sqrt {\theta }} \overline {\boldsymbol{\nabla} _{\!j} \boldsymbol{u}_k}^{\sqrt {\theta }}}^{\phi } -\overline { \overline {\boldsymbol{\nabla} _i \boldsymbol{u}_k}^{\sqrt {\theta }}}^{\phi } \, \overline { \overline {\boldsymbol{\nabla} _{\!j} \boldsymbol{u}_k}^{\sqrt {\theta }}}^{\phi } \right )\!, \end{equation}
where
$\phi =\sqrt {l^2-\theta }$
. The first term is scale local and corresponds to the Clark model, which therefore neglects the higher-order non-local terms involving scales smaller than the filter one. The difference between the Clark model and DNS in the energy flux is related to this term. We emphasise that the discussion is made at the filter scale
$q=16$
, that is in the inertial range. We can see that the local-in-scale model globally underestimates the mean. In figure 3, we move beyond the analysis of the mean and examine the distribution of the Clark prediction versus the ground-truth energy flux conditioned on five representative points in the
$(Q,R)$
-space: the origin and one point in each of the four quadrants, highlighted in figure 2. This assessment is severe, but it permits us to fully recognise the capability to reproduce the cascade process by a model. This type of conditioning is particularly interesting since it corresponds to imposing a constraint on the large-scale flow topology. The solid lines show the conditioned PDF for the measured energy flux directly from the DNS, while the dashed lines give the conditioned PDFs for the predicted energy flux from the Clark model. We recall that
$Q,R$
do not uniquely determine the energy flux in the Clark model and that, for the same
$Q,R$
, a distribution of energy fluxes is observed. Globally speaking, the Clark model makes a decent job overall. The essential features of the flux are captured, and in particular we observe a good agreement in the left/backscatter tail. Nonetheless, the skewness is not correct and the right/positive tail is systematically underestimated. The shape of the PDF appears similar, but a more careful look at the right tails shows that the scaling is also different. Overall, stronger fluctuations are observed in the DNS than the Clark model predictions. Considering the decomposition in scales (3.2), the non-local terms appear important in the positive tails of the flux distributions.
Colour scale represents the mean local energy flux conditioned on
$(Q,R)$
configurations at filtering scale
$q=16$
(corresponding to
$l=1/q$
). Positive fluxes (forward cascade) are shown in red, and negative fluxes (backscatter) in blue. Black regions correspond to events with probability less than 10
$^{-8}$
. The (a) displays DNS results, while the (b) shows the Clark model prediction. White contours indicate isolines of the probability in the
$(Q,R)$
plane at levels 10
$^{-4}$
, 10
$^{-5}$
, 10
$^{-6}$
and 10
$^{-7}$
, moving outward from the origin. The purple curve denotes the isoline of zero flux. Coloured dots in the (b) indicate the
$(Q,R)$
locations used for the conditional analysis in figure 3.
Comparison of local energy flux PDFs conditioned on
$(Q,R)$
configurations, between ground-truth DNS (solid lines) and Clark model predictions (dashed lines). The corresponding
$(Q,R)$
locations are indicated in figure 2.
Colour scale represents the joint PDFs of (2.6) and (2.7). White contours indicate probability isolines at levels 10
$^{-4}$
, 10
$^{-5}$
, 10
$^{-6}$
and 10
$^{-7}$
, progressing outward from the origin. Green upward tripod markers denote the peaks of the conditional PDFs obtained by fixing the model prediction (horizontal cuts), while blue circles represent the corresponding mean values.
First and third rows: statistical cumulants (mean, variance, skewness and flatness) of the DNS density, computed conditionally on a given value of the Clark flux. Horizontal solid lines indicate the unconditional values. Middle row: the conditional variance shown in logarithmic scale (b, to highlight variability) and the relative statistical error
$\sigma /|\mu |$
(a). In the top-left panel (mean), two linear fits with different slopes, respectively for positive and negative values, are indicated by green triangles and red circles.
3.1.2. Joint PDFs
To provide a more precise characterisation of fluctuations, we introduce the joint probability density of the local energy flux (2.6) and the Clark model prediction (2.7). This joint PDF serves as the appropriate observable to assess the fluctuations around the model: for a fixed value of the Clark prediction, the complete flux exhibits a distribution of values, reflecting the influence of the chaotic small-scale dynamics. Our primary objective is to investigate the significance of these fluctuations. In figure 4 we present this joint PDF. The white contours have the same meaning as in the previous plot. Green markers indicate the maxima of the conditional PDFs for fixed model values (horizontal cuts), while blue circles represent the corresponding mean values. The discrepancies between the mode and the mean provide a first indication of non-Gaussian behaviour. Furthermore, to a good approximation the means lie on a linear curve, although the slope differs slightly between the negative and positive ranges, as quantified in figure 6. These evidences corroborate the previous conclusion that the Clark model performs well: it is possible to adjust a multiplicative constant (or, more precisely, two constants depending on whether the cascade is forward or backward) in front of (2.7) so that the model matches the average of the DNS configurations. In the inertial range, the linear behaviour with slope of order one shows that the contributions of the non-local terms are negligible with respect to the mean flux (Johnson Reference Johnson2021). Examples of selected conditional distributions, given the model prediction (i.e. the resolved field), are presented in figure 5. These correspond to horizontal cuts in figure 4. Large fluctuations with a variance of the order of the mean are observed. Then, the PDFs are characterised by pronounced non-Gaussian fat tails for each value of the flux. These heavy-tailed statistics indicate the occurrence of rare but intense deviations from the average behaviour, highlighting the intermittent and strongly nonlinear nature of the system. The observed behaviour suggests that the tails could follow a distribution intermediate between exponential and power-law decay. Yet, it appears that there is not a universal shape for the PDF, and we have not been able to rescale toward a sole master curve. To quantify these observations, we compute the first four cumulants of the conditional PDFs, varying the value of the model. These are mean and variance
\begin{equation} \mu =\int\! \varPi \,{\rm d}\mathbb{P}[\varPi \vert \varPi _\textit{model}],\quad \sigma ^2=\int (\varPi -\mu )^2\,{\rm d}\mathbb{P}[\varPi \vert \varPi _\textit{model}], \end{equation}
skewness and flatness
\begin{equation} \gamma _1=\frac {1}{\sigma ^3}\int (\varPi -\mu )^3\,{\rm d}\mathbb{P}[\varPi \vert \varPi _\textit{model}],\quad \gamma _2=\frac {1}{\sigma ^4}\int (\varPi -\mu )^4\,{\rm d}\mathbb{P}[\varPi \vert \varPi _\textit{model}], \end{equation}
where
$\varPi$
denotes the DNS value (2.6). As a first non-trivial result, we note that, on average, the DNS flux has the same sign as the Clark model prediction. This allows us to fit two linear slopes, depending on the sign, both constrained to pass through the
$(0,0)$
flux point, yielding
$1.4 \pm 0.2$
for the forward cascade and
$0.5 \pm 0.2$
for the backward case. Continuing the analysis, and interpreting the variance as a measure of statistical uncertainty, we observe that extreme events (i.e. those with large absolute values) are associated with greater variability. Interestingly, the ratio
$\sigma /\vert \mu \vert$
appears to decrease as one moves away from zero. The skewness, however, clearly indicates that the distribution remains not symmetric. Gaussian flatness is recovered for extreme Clark model predictions, whereas fat-tailed behaviour dominates near the most probable events.
3.2. Smagorinsky vs DNS
For completeness we perform a parallel analysis for the classical eddy-viscosity model of Smagorinsky, described by expression (2.8). Unlike the Clark model, the Smagorinsky formulation predicts only non-negative energy flux, thereby completely failing to account for the backscatter. In figure 7, we present the joint PDF of the ground truth and model prediction: a multiplicative constant can be fitted to align either the maximum probability or the mean value, as shown in figure 9. As in figure 5, figure 8 displays the conditional PDFs of the DNS energy flux for fixed values of the Smagorinsky output. These distributions exhibit prominent heavy tails, characteristic of rare but intense events. A fit of the conditional distributions reveals that the tails decay more slowly than exponentially, supporting the idea that rare events dominate the flux statistics and cannot be neglected in closure strategies. To further quantify these features, we look at the first four statistical cumulants of the conditional distributions, shown in figure 9. As with the Clark model, the conditional mean has a linear trend versus the model prediction: a linear fit gives
$1.1 \pm 0.4$
, where the uncertainty corresponds to the root-mean-square deviation. Yet the variance increases with increasing model flux much more for the Smagorinsky case, indicating greater uncertainty for more intense events. However, in that regime, the relative uncertainty
$\sigma /\vert \mu \vert$
decreases. The skewness exhibits a pronounced positive trend, highlighting the strong asymmetry of the distributions. As for the flatness, values close to the Gaussian reference (i.e.
$3$
) are approached at large predicted values, while extreme flatness, indicative of intermittent fat-tailed fluctuations, dominates in a region of width
$\epsilon$
around
$\varPi _{\textit{Smag.}}=0$
, where
$\epsilon$
denotes the mean energy injection rate. These findings illustrate that the behaviour near the Smagorinsky prediction is not Gaussian.
Colour scale represents the joint PDFs of (2.6) and (2.8). White contours indicate probability isolines at levels 10
$^{-4}$
, 10
$^{-5}$
, 10
$^{-6}$
and 10
$^{-7}$
, progressing outward from the origin. Green upward tripod markers denote the peaks of the conditional PDFs obtained by fixing the model prediction (horizontal cuts), while blue circles represent the corresponding mean values.
First and third rows: statistical cumulants (mean, variance, skewness and flatness) of the DNS density, computed conditionally on a given value of the Smagorinsky prediction. Horizontal solid lines indicate the unconditional values. Middle row: the conditional variance shown in logarithmic scale (b, to highlight variability) and the relative statistical error
$\sigma /|\mu |$
(a). In the top-left panel (mean), a linear fit constrained to the
$(0,0)$
flux point is indicated by green triangles.
4. Discussion
We have investigated different conditional statistics of the filtered turbulent velocity field in three dimensions. The analysis of the mean energy flux in the QR diagram shows unambiguously that the backscatter inverse cascade is correlated solely with the vortex compression, while positive flux is due to all quadrants. Results highlight the role of the non-local terms, which are found to be responsible for an increase of the positive extreme events of the energy flux. With respect to the LES framework, we have found that the Clark model reproduces satisfactorily the mean energy flux and provides an approximate reconstruction of the fluctuations around the mean, but it remains insufficient to fully capture the complex statistical structure of the energy flux in turbulent flows. While the Clark model can be considered a good first-order approximation, we have highlighted the presence of heavy-tailed distributions and marked asymmetries, hallmarks of the intermittent nature of the modelled unresolved scales.
These findings underscore the limitations of classical deterministic closures based solely on filtered quantities and point toward the need for stochastic components to better account for the influence of unresolved scales. The strong skewness and enhanced flatness observed at moderate model outputs demonstrate that rare events and intermittent bursts remain statistically significant and cannot be neglected. In this context, fluctuations in the subgrid stress tensor cannot be modelled by simple Gaussian noise. Starting from available stochastic modelling for LES (Gicquel et al. Reference Gicquel, Givi, Jaberi and Pope2002), our results provide a quantitative foundation for the development of stochastic subgrid-scale models that can incorporate both the mean behaviour and the statistical variability of energy transfer.
We need to remark that, in this work, we have focused on stationary statistics to construct time-independent PDFs. An important direction for future research is to investigate temporal correlations, with the aim of constructing fully dynamical stochastic processes for the subgrid-scale energy flux. Temporal correlations indeed constitute a critical point. For instance, in the limiting case of a delta-correlated-in-time closure, a dependence on the integration time step arises: if the random variable is sampled much more frequently than the eddy turnover time of the LES field, the deterministic dynamics effectively perceives only its mean. As discussed in previous works on stochastic modelling (Dreeben & Pope Reference Dreeben and Pope1998; Minier, Chibbaro & Pope Reference Minier, Chibbaro and Pope2014; Innocenti, Marchioli & Chibbaro Reference Innocenti, Marchioli and Chibbaro2016), physical and numerical constraints on the time step may be needed. We regard this as an important question to be addressed in future research.
We conclude by outlining the perspective of a data-driven approach to capture fluctuations around simple deterministic models. The main present challenge is the curse of dimensionality in the three-dimensional setting: collecting sufficient data to train a generative model for full three-dimensional Eulerian fields is practically impossible, and convergence of the training is also difficult to achieve. A more tractable strategy, in our view, is to reproduce the fluctuations around a first-order approximation expressed in terms of the velocity gradient, rather than attempting to model the entire subgrid-scale stress tensor. The neglected contributions can be written in terms of higher-order gradients, and locality allows the training to be well conditioned when restricted to a neighbourhood of data, thereby improving stability. Unfortunately, the powerful neural networks underlying state-of-the-art generative models, such as diffusion models, remain essentially uninterpretable from a physical point of view, leaving them as promising but opaque tools.
Acknowledgements
This work received financial support from the CNRS through the MITI interdisciplinary initiatives, under its exploratory research program. S.C. acknowledges the funding by ANR SPEED project ANR-20-CE23-0025-01. It was granted access to the HPC resources of GENCI-TGCC & GENCI- CINES (Project No. A0170506421). Finally, we would like to thank L. Biferale for useful discussions.
Declaration of interests
The authors report no conflict of interest.
Open access
