2.1. Navier–Stokes–Smagorinsky equation

In practice, the closure problem takes as inputs a filtering scale $\delta$ and a filter $G$ . It aims to predict the dynamics of the filtered field as the convolution $\boldsymbol{u}_{\delta } = G_\delta * \boldsymbol{u}$ with the kernel $G_\delta (\boldsymbol{x}) = \delta ^{-3} G(\boldsymbol{x}/\delta )$ , providing implicit or explicit parametrisations for the so-called subgrid-scale tensor $\boldsymbol{\rm R}_{\delta} = (\boldsymbol{u} \otimes \boldsymbol{u} )_{\delta }-\boldsymbol{u}_{\delta }\otimes \boldsymbol{u}_{\delta }$ . We start with the simplest example of the Smagorinsky model, in which ${\boldsymbol{R}}_\delta$ is expressed as the viscous stress with the effective eddy viscosity $\nu _s$ . Dropping $\delta$ from the subscripts, the resulting Navier–Stokes–Smagorinsky (NSS) dynamics for the filtered velocity becomes

(2.1) \begin{equation} \partial _t \boldsymbol{u} +\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla } \boldsymbol{u} = -\boldsymbol{\nabla } p+\boldsymbol{\nabla } \boldsymbol{\cdot} \left (2 \nu _s \boldsymbol{\rm s}\right )+{\boldsymbol{f}}, \quad \boldsymbol{\nabla } \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{equation}

where $p(\boldsymbol{x},t)$ is the pressure and ${\boldsymbol{f}}(\boldsymbol{x},t)$ is a large-scale force with correlation length $L_f\gg \delta$ . The eddy viscosity is defined in terms of the filtered rate-of-strain tensor $\boldsymbol{s}$ as

(2.2) \begin{equation} \nu _s = (c_s \delta )^2 \|\boldsymbol{\rm s}\|,\quad \boldsymbol{\rm s} = \tfrac {1}{2} \left ( \boldsymbol{\nabla } \boldsymbol{u} + (\boldsymbol{\nabla } \boldsymbol{u})^\textrm T\right ), \end{equation}

with the Smagorinsky constant $c_s = 0.1$ and the norm $\|\boldsymbol{\rm s}\|^2 = \sum 2\mathrm{s}_{ij}^2$ . Considering the filtering scale $\delta$ to be much larger than the viscous scale, the kinematic viscosity $\nu \ll \nu _s$ is small and, therefore, omitted in the NSS equations (2.1). The filtering operation leading to (2.2) is formal only. While the NSS dynamics aims to mimic the NS) equations for the resolved scales $\ell \gg \delta$ , the role of eddy viscosity is to replace the energy transfer at the subgrid scales $\ell \lesssim \delta$ by a dissipative scheme (Cruz & Lamballais Reference Cruz and Lamballais2023).

2.2. Kinematic rescaling

Similarly to the NS system previously discussed by Mailybaev (Reference Mailybaev2022a ), Mailybaev & Thalabard (Reference Mailybaev and Thalabard2022), the starting point of our study is the projection of the velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ into the rescaled velocity field $\boldsymbol{U}(\boldsymbol{X},\tau )$ . At a heuristic level, the new field represents a quasi-Lagrangian velocity increment (Belinicher & L’vov Reference Belinicher and L’vov1987) with a dynamically rescaled amplitude and timeframe. The course of proper time $\tau$ is determined by the local level of fluctuations, running fast when the latter are high and slowly otherwise. The rescaling is parametrised by an observational scale $\ell$ and a label $\boldsymbol{x}_0 \in \mathbb{R}^3$ . The label selects a Lagrangian trajectory as

(2.3) \begin{equation} \dfrac {{\rm d}\boldsymbol{x}_*}{{\rm d}t} = {\boldsymbol{u}}(\boldsymbol{x}_*,t),\quad \boldsymbol{x}_*(0)= \boldsymbol{x}_0. \end{equation}

Fluctuations at scale $\ell$ are measured as amplitudes of quasi-Lagrangian increments as

(2.4) \begin{equation} a_\ell (t) = \sqrt {\int _{|\boldsymbol{X}| \leqslant 1} \big |\varDelta _\ell \boldsymbol{u}\big (\boldsymbol{x}_*(t),\boldsymbol{X},t\big ) \big |^2 d^3\boldsymbol{X}}, \quad \varDelta _\ell \boldsymbol{u}(\boldsymbol{x},\boldsymbol{X}, t) = \boldsymbol{u}\big (\boldsymbol{x}+\ell \boldsymbol{X},t\big )-\boldsymbol{u}\big (\boldsymbol{x},t\big ). \end{equation}

The space variable $\boldsymbol{X}$ represents the vectorial separation with respect to the Lagrangian trajectory in units of the observational scale $\ell$ . The velocity field is then rescaled as

(2.5) \begin{equation} \boldsymbol{U}(\boldsymbol{X},\tau ) = \dfrac {\varDelta _\ell \boldsymbol{u}\big (\boldsymbol{x}_*(t),\boldsymbol{X},t\big )}{a_\ell (t)} ,\quad \tau = \frac {1}{\ell }\int _0^t a_\ell (s) \, {\rm d}s. \end{equation}

Here the first formula introduces the rescaled velocity $\boldsymbol{U}$ depending on the rescaled time $\tau$ , which is expressed in terms of the velocity difference $\varDelta _\ell \boldsymbol{u}$ at original time $t$ . The second formula defines the relation between $t$ and $\tau$ , which is one-to-one for $a_\ell (t) \gt 0$ . Physically, $\tau$ represents a Lagrangian clock measuring time in units of the local turnover time $ \ell /a_\ell (t)$ encountered along the trajectory $\boldsymbol{x}_*(t)$ .

Note that the choice of $a_\ell$ in (2.4) is by no means unique. In practice, the essential requirement is that it should be homogeneous of degree 1 with respect to the velocity field and positive along a trajectory. Here, for simplicity, we restrict our discussion to a specific example rather than deriving the hidden symmetry in full generality. We refer to Mailybaev & Thalabard (Reference Mailybaev and Thalabard2022) for more general derivations.

2.3. Switching the observation scale

Changing the scale $\ell \mapsto \ell '$ in (2.4) and (2.5) reduces to a specific transformation of the respective fields $\boldsymbol{U}(\boldsymbol{X},\tau ) \mapsto \boldsymbol{U}'(\boldsymbol{X},\tau ')$ expressed as

(2.6) \begin{equation} \boldsymbol{U}'(\boldsymbol{X},\tau ') = \frac {\boldsymbol{U}(\boldsymbol{X}/\lambda ,\tau )}{A(\tau )}, \ \ \tau ' = \lambda \int _0^\tau A(s) \, {\rm d}s, \ \ A(\tau ) = \sqrt {\int _{|\boldsymbol{X}| \leqslant 1} | \boldsymbol{U}(\boldsymbol{X}/\lambda ,\tau ) |^2 d^3\boldsymbol{X}}, \end{equation}

where $\lambda = \ell /\ell '$ and the second relation introduces the change of time $\tau \mapsto \tau '$ . The kinematic relations (2.5) and (2.6) are demonstrated schematically in figure 2(a) by straight and wavy arrows, respectively. In this figure, we specified the parameters in the subscripts: $\boldsymbol{u}_\delta$ for the original field depending on the filtering scale $\delta$ and $\boldsymbol{U}_{\delta ,\ell ,\boldsymbol{x}_0}$ for the rescaled field depending additionally on the observation scale $\ell$ and the label $\boldsymbol{x}_0$ .

Figure 2.

(a) Kinematic rescaling diagram, where the projections (2.5) are shown by straight arrows for two different scales $\ell$ and $\ell ' = \lambda \ell$ , and the transformation (2.6) is shown by the wavy arrow. (b) Similar diagram for a different filtering scale $\delta ' = \delta /\lambda$ . (c) Manifestation of the hidden self-similarity as the fusion of the previous diagrams in their right-hand sides.

2.4. Navier–Stokes–Smagorinsky under kinematic rescaling

The governing equations for the rescaled velocity field $\boldsymbol{U}(\boldsymbol{X},\tau )$ are obtained by differentiating the definition (2.5) with respect to $\tau$ and using the NSS system (2.1) and (2.2). The lengthy but elementary derivation yields (see Appendix A.1)

(2.7) \begin{equation} \partial _{\tau } \boldsymbol{U} = \Lambda _{\boldsymbol{U}} \big [-\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla } \boldsymbol{U}-\boldsymbol{\nabla } P+2 c_s^2 \,\alpha ^{-2}\, \boldsymbol{\nabla } \boldsymbol{\cdot} (\|\boldsymbol{\rm S}\| \boldsymbol{\rm S})+\boldsymbol{F} \big ], \quad \boldsymbol{\nabla } \boldsymbol{\cdot} {\boldsymbol{U}}=0. \end{equation}

Here P is the rescaled pressure, $\alpha = \ell /\delta$ , the tensor $\boldsymbol{\rm S} = ({1}/{2}) \boldsymbol{\nabla } \boldsymbol{U} + ({1}/{2})(\boldsymbol{\nabla } \boldsymbol{U})^\textrm T$ , and the gradient operator $\boldsymbol{\nabla }$ acts on the new space variable $\boldsymbol{X}$ . The operator $\Lambda _{\boldsymbol{U}}$ is defined for any field ${\boldsymbol{V}}(\boldsymbol{X},\tau )$ as

(2.8) \begin{equation} \Lambda _{\boldsymbol{U}}[{\boldsymbol{V}}] = \tilde {{\boldsymbol{V}}}-\boldsymbol{U} \int _{|\boldsymbol{X}| \leqslant 1} \boldsymbol{U} \boldsymbol{\cdot} \tilde {{\boldsymbol{V}}} d^3\boldsymbol{X}, \quad \tilde {{\boldsymbol{V}}} = {\boldsymbol{V}}(\boldsymbol{X},\tau )-{\boldsymbol{V}}(\boldsymbol{0},\tau ). \end{equation}

The rescaled forcing term is expressed, with the times $t \mapsto \tau$ related by (2.5), as

(2.9) \begin{equation} \boldsymbol{F}(\boldsymbol{X},\tau ) = \frac {\ell \, \Delta _\ell {\boldsymbol{f}}\big (\boldsymbol{x}_*(t),\boldsymbol{X},t\big )}{a_\ell ^2(t)}, \quad \varDelta _\ell {\boldsymbol{f}}(\boldsymbol{x},\boldsymbol{X},t) = {\boldsymbol{f}}(\boldsymbol{x}+\ell \boldsymbol{X},t)-{\boldsymbol{f}}(\boldsymbol{x},t). \end{equation}

3.1. The limit of small observational scales

System (2.7) governs the dynamics of the new velocity and pressure fields $\boldsymbol{U}(\boldsymbol{X},\tau )$ and $P(\boldsymbol{X},\tau )$ . It is, however, not closed with respect to these fields, because the coupling with the original dynamics $\boldsymbol{u}(\boldsymbol{x},t)$ permeates through the dependence of forcing term (2.9) on $a_\ell (t)$ . We now argue that a closed formulation follows asymptotically at small observational scales $\ell \ll L_f$ , in which case the forcing term is negligible. Indeed, we have $\partial _\tau \boldsymbol{U} \sim 1$ for $\boldsymbol{X} \sim 1$ by construction from (2.4) and (2.5). This implies that the left-hand side in (2.7) is ${\sim}1$ . As for the forcing term, we first consider the case $\delta \lesssim \ell \ll L_f$ , describing an observational scale in the resolved (turbulent) range. The dimensional estimate ${\boldsymbol{f}} \sim U_f^2/L_f$ in (2.9) yields $\varDelta _\ell {\boldsymbol{f}} \sim U_f^2 \ell /L_f^2$ , with $U_f$ the integral (forcing-scale) velocity. Using the Kolmogorov 1941 (K41) estimate $a_\ell \sim U_f (\ell /L_f)^{1/3}$ for the amplitude of velocity fluctuations in (2.9), we conclude that the forcing term $\boldsymbol{F} \sim (\ell /L_f)^{4/3} \ll 1$ . This estimate extends to the subgrid range $\ell \lt \delta \ll L_f$ as $\boldsymbol{F} \sim (\delta /L_f)^{4/3} \ll 1$ .

Therefore, under the condition of small scales $\delta , \ell \ll L_f$ , the forcing term can be neglected in the rescaled system (2.7). This gives

(3.1) \begin{equation} \delta , \ell \ll L_f: \quad \partial _\tau \boldsymbol{U} = \Lambda _{\boldsymbol{U}}\big [-\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla } \boldsymbol{U}-\boldsymbol{\nabla } P+ 2 c_s^2 \alpha ^{-2} \boldsymbol{\nabla } \boldsymbol{\cdot} (\|\boldsymbol{\rm S}\| \boldsymbol{\rm S}) \big ], \quad \boldsymbol{\nabla } \boldsymbol{\cdot} {\boldsymbol{U}}=0. \end{equation}

We call the system (3.1) $\alpha$ -NSS, emphasising that the dependence on the parameters $\ell$ and $\delta$ is realised only through their ratio $\alpha = \ell /\delta$ . Based on these properties, we formulate the hypothesis of hidden scale invariance:

3.2. Statistical identities

The diagrammatic interpretation of this type of scale invariance is presented in figure 2. Panels (a) and (b) show two independent diagrams corresponding to different filtering parameters $\delta$ and $\delta ' = \lambda \delta$ . When viewed at small scales and in a statistical sense, these diagrams merge on their right sides as shown in panel (c), where $\boldsymbol{U}_\alpha$ denotes a solution of the $\alpha$ -NSS system (3.1). In other terms, the hidden scale invariance is understood as the statistical equivalence

(3.2) \begin{equation} \boldsymbol{U}_{\delta ,\ell ,\boldsymbol{x}_0} \stackrel {\mathrm{law}}{=} \, \boldsymbol{U}_{ \lambda \delta ,\lambda \ell ,\boldsymbol{x}'_0} \stackrel {\mathrm{law}}{=} \boldsymbol{U}_\alpha \ \ \textrm {for} \ \ \alpha = \ell /\delta . \end{equation}

The diagram in figure 2(c) also reveals that the rescaled statistics for different values of $\alpha$ and $\alpha ' = \alpha /\lambda$ are related by the transformation $h_\lambda : \boldsymbol{U}_\alpha \mapsto \boldsymbol{U}_{\alpha '}$ prescribed by (2.6). This provides an additional statistical identity

(3.3) \begin{equation} \boldsymbol{U}_{\alpha /\lambda } \stackrel {\mathrm{law}}{=} h_\lambda [\boldsymbol{U}_\alpha ]. \end{equation}

3.3. Hidden scale invariance in the inertial interval

A stronger form of HS develops in the inertial interval, i.e. at scales $\ell$ defined by the condition $\delta \ll \ell \ll L_f$ . Here, the Smagorinsky dissipative term becomes negligible and the system (3.1) simplifies into

(3.4) \begin{equation} \delta \ll \ell \ll L: \quad \partial _\tau \boldsymbol{U} + \Lambda _{\boldsymbol{U}} \big [\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla } \boldsymbol{U}+\boldsymbol{\nabla } P \big ]=0, \quad \boldsymbol{\nabla } \boldsymbol{\cdot} {\boldsymbol{U}}=0. \end{equation}

This system is independent of both $\delta$ and $\ell$ , and coincides with the rescaled system derived earlier in the framework of NS turbulence by Mailybaev (Reference Mailybaev2022a ), Mailybaev & Thalabard (Reference Mailybaev and Thalabard2022). This leads to the hypothesis of inertial-interval hidden symmetry:

This suggests that, despite the intermittency, the rescaled NSS and NS statistics coincide at the inertial-interval scales. The statistical identity (3.3) in the inertial interval simplifies into

(3.5) \begin{equation} \boldsymbol{U}_\infty \stackrel {\mathrm{law}}{=} h_\lambda [\boldsymbol{U}_\infty ], \end{equation}

where $\boldsymbol{U}_\infty$ denotes a solution of system (3.4) obtained in the limit $\alpha \to \infty$ .

3.4. Hidden scale invariance in other SGS models

We observe that the concept of HS applies equally to hyperviscous Smagorinsky models, such as $\propto \! \delta ^{2p} |\boldsymbol{\nabla }^2|^p \boldsymbol{\nabla } \boldsymbol{\cdot} (\nu _s \boldsymbol{\rm s})$ . More interestingly, modern strategies in LES use implicit dissipation schemes in order to introduce backscattering effects; see, e.g. Pope (Reference Pope2000), Lamballais, Cruz & Perrin (Reference Lamballais, Cruz and Perrin2021). In practice, these schemes are implemented by under-resolving the NSS through a (genuine) filtering operation, hereby transforming the NSS equations (2.1) into

(3.6) \begin{equation} \partial _t \boldsymbol{u} = \big (-\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla } \boldsymbol{u}-\boldsymbol{\nabla } p+\boldsymbol{\nabla } \boldsymbol{\cdot} \left (2 \nu _s \boldsymbol{\rm s}\right )+{\boldsymbol{f}}\big )_{\delta }, \quad \boldsymbol{\nabla } \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{equation}

where the right-hand side is the filtered field. Consider the sharp spectral filter $G_\delta (\boldsymbol{x}) = (2\pi )^{-3}\int _{|{\boldsymbol{k}}|\leqslant c k_\delta } \! d^3 {\boldsymbol{k}} e^{i{\boldsymbol{k}}\boldsymbol{\cdot} \boldsymbol{x}}$ defined as a Galerkin truncation at a given multiple $c$ of the wavenumber $k_\delta = 2\pi /\delta$ . Then the kinematic rescaling of § 2 yields the filtered $\alpha$ -NSS as

(3.7) \begin{equation} \partial _\tau \boldsymbol{U} = \big (\Lambda _{\boldsymbol{U}} [ -\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla } \boldsymbol{U}-\boldsymbol{\nabla } P+2c_s^2 \alpha ^{-2}\, \boldsymbol{\nabla } \boldsymbol{\cdot} (\|\boldsymbol{\rm S}\| \boldsymbol{\rm S})+\boldsymbol{F} ]\big )_{\alpha }, \quad \boldsymbol{\nabla } \boldsymbol{\cdot} {\boldsymbol{U}} =0 \end{equation}

(see Appendix A.2). Here the filtering operation $(\ldots )_\alpha$ applies with respect to the rescaled variable $\boldsymbol{X}$ at scale $\alpha = \ell /\delta$ . Hence, the HS extends to implicit LES schemes. Following Mailybaev (Reference Mailybaev2022a ) one can relate the HS of SGS models to their invariance with respect to space–time scaling and Galilean symmetry groups. From this point of view, the extension of HS to subgrid (dissipative) scales refers to the time-scale-invariant (quadratic in velocity) SGS model. On the contrary, the HS is limited to the inertial scales of the NS, because physical viscosity breaks the time-scale invariance (Mailybaev & Thalabard Reference Mailybaev and Thalabard2022).

4.1. Numerical set-up

We probe the hidden self-similarity conjectures through direct numerical simulations of the NSS in a triply $2\pi$ -periodic domain. Our simulations use a pseudospectral scheme with 2/3-rule dealiasing (Orszag Reference Orszag1971) and a second-order Adam–Bashforth time integration with step $\varDelta t$ . The forcing ${\boldsymbol{f}}(\boldsymbol{x},t) = \boldsymbol{\nabla } \times \psi$ is specified through the stream vector $\boldsymbol \psi = \sqrt {\epsilon } \sum _{ 1 \leqslant |{\boldsymbol{k}}| \leqslant 2\pi /L_f} e^{i {\boldsymbol{k}} \boldsymbol{\cdot} \boldsymbol{x}} {\boldsymbol \eta ({\boldsymbol{k}},t)}/(\sum _{1 \leqslant |{\boldsymbol{k}}| \leqslant 2\pi /L_f} |{\boldsymbol{k}}|^2)$ . Here $\epsilon = 0.14$ and $L_f = 4\pi /3$ are the energy injection rate and scale, and $\boldsymbol \eta ({\boldsymbol{k}},t)$ are white-in-time independent complex standard Gaussian random vectors satisfying the Hermitian symmetry $\boldsymbol \eta ({\boldsymbol{k}})=\boldsymbol \eta ^*(-{\boldsymbol{k}})$ . We describe two sets of runs with resolution $N$ ranging from $256^3$ to $1024^3$ , referred to as explicit and implicit runs. The explicit runs use $\delta \gg \varDelta x$ ensuring that dissipation in the subgrid range is dominated by the Smagorinsky damping term, while the implicit runs use $\delta \sim \varDelta x$ and allow for energy backscatter; see table 1 for more details. For numerical convenience, we use $|\boldsymbol{X}| = \max _j|X_j|$ in the amplitudes (2.4) corresponding to averages in a cube.

Table 1.

Numerical details for two sets of runs with $k_{max }=N/3$ and $\varDelta x = 2\pi /N$ using the filtered NSS equation (3.6) at resolutions $ N^3 = 256^3, 512^3, 1024^3$ . The total integration time is $T_{max } = 15$ for resolutions up to $512^3$ and $T_{max }=4.6$ for $1024^3$ .

4.2. Temporal averages

The statistics of the rescaled field $\boldsymbol{U}(\boldsymbol{X},\tau )$ can be accessed using different observables. In this work, we limit our analysis to the longitudinal component $U_\parallel (\tau ) = \boldsymbol{U}({\boldsymbol{e}},\tau )\boldsymbol{\cdot} {\boldsymbol{e}}$ with the unit vector ${\boldsymbol{e}} = (1,0,0)$ ; other observables (not presented here) were also tested with similar conclusions. The $\tau$ -average of $|U_\parallel |^p$ is defined as

(4.1) \begin{equation} \left \langle {|U_\parallel (\tau )|^p}\right \rangle _\tau =\lim _{T\to \infty }\dfrac {1}{T}\int _0^T |U_\parallel (\tau )|^p {\rm d}\tau , \end{equation}

from which the usual statistical properties such as variance, skewness, flatness, etc. are estimated. The corresponding probability density function (PDF) is obtained from the average of the Dirac delta function as $p_{\parallel }(y)= \left \langle {\delta (U_\parallel (\tau )-y )}\right \rangle _\tau$ . We emphasise that the average (4.1) is expressed in terms of the rescaled Lagrangian clock $\tau$ . However, it can be efficiently estimated as space–time average over the original fields using a conversion formula relying on spatial homogeneity and incompressibility. For $|U_\parallel |^p$ the conversion formula states

(4.2) \begin{equation} \left \langle {|U_\parallel (\tau )|^p}\right \rangle _\tau = \left \langle {\left |\dfrac {\varDelta _\ell u_\parallel (\boldsymbol{x},t)}{a(\boldsymbol{x},t)}\right |^p\, J(\boldsymbol{x},t)}\right \rangle _{\boldsymbol{x},t}, \end{equation}

where $\left \langle {\boldsymbol{\cdot} }\right \rangle _{\boldsymbol{x},t}$ denotes the space–time average in the original coordinates and $\varDelta _\ell u_\parallel = ( \boldsymbol{u}(\boldsymbol{x} +\ell {\boldsymbol{e}}) - \boldsymbol{u}(\boldsymbol{x}) ) \boldsymbol{\cdot} {\boldsymbol{e}}$ . In other words, (4.2) transforms quasi-Lagrangian averages into usual (Eulerian) ones, that one can easily compute without having to track any fluid trajectory.

Here $J(\boldsymbol{x},t)$ is the Jacobian factor due to the change of time $\tau \mapsto t$ as

(4.3) \begin{equation} J(\boldsymbol{x},t) = \frac {a(\boldsymbol{x},t)}{\left \langle {a(\boldsymbol{x},t)}\right \rangle _{\boldsymbol{x},t}}, \quad a(\boldsymbol{x},t) = \Big (\int _{|\boldsymbol{X}| \leqslant 1} \big | \varDelta _\ell \boldsymbol{u}(\boldsymbol{x},\boldsymbol{X},t) \big |^2 d^3\boldsymbol{X}\Big )^{1/2}. \end{equation}

We refer to Mailybaev & Thalabard (Reference Mailybaev and Thalabard2022) for technical details of this derivation. Expressions similar to (4.2) follow for other single-time observables.

4.3. Numerical results

Figure 3(a) illustrates the HS by verifying the collapse at the level of PDFs of $U_\parallel$ for increasing prescribed values of the parameter $\alpha = \ell /\delta$ . The top and bottom graphs correspond to the explicit and implicit runs, respectively, and we recall that different resolutions feature different $\delta$ ; see table 1. In both cases, the curves collapse at small scales, verifying the hidden self-similarity: the statistics depend only on the ratio $\ell /\delta$ . As $\alpha$ increases, the graphs begin to diverge from one another, with the smallest-resolution run departing from the other two. This signals a sensitivity of the statistics to the large (forcing) scales, and this is compatible with the HS getting broken at scales $ O(L_f/\delta )$ .

Figure 3(b) further characterises this large-scale sensitivity, by showing the variance and flatness of $U_\parallel (\tau )$ as functions of $\alpha = \ell /\delta$ . Looking at the statistics from large to small scales (right to left), the gradual convergence is most visible at the level of the flatness, signalling that the PDF tails require a larger inertial range to converge – see also the top graphs in each of figure 3(a).

Figure 3.

(a) The PDF of $ U_\parallel$ at selected values of $\ell /\delta$ for the various runs. The data is shifted horizontally for clarity with the vertical lines indicating the shifts. The top graphs in each show the corresponding data using logarithmic vertical axis. (b) Variance and flatness of $U_\parallel$ as a function of $\alpha = \ell /\delta$ . The vertical lines indicate the values of $\alpha$ used in (a). Top and bottom graphs use, respectively, the explicit and implicit simulations.

Explicit runs use large $\delta$ and therefore do not develop the inertial range; their HS applies to dissipative scales only. Implicit runs have smaller $\delta$ , and a pronounced inertial interval develops at the resolution $1024^3$ for the scales $6 \lesssim \ell /\delta \lesssim 42$ ; see the plateau in the bottom figure 3(b) for both the variance and the flatness. Note that the flatness curves do not collapse well in this range because of the shorter inertial interval for lower resolutions. The PDFs of $U_\parallel$ for 18 different scales $\ell$ from the inertial interval are shown by solid lines in figure 4(a); these scales $0.01\lt \ell /L\lt 0.1$ correspond to the squares in figure 4(b). They also correspond to the range over which a flatness scaling $\propto \ell ^{-0.12}$ is clearly observable in figure 1(b). In figure 4(a), these PDFs are visually indistinguishable, confirming the inertial-interval HS. The collapse is quantified in figure 4(b), where the black squares show the Euclidean distance ( $L_2$ -norm) between the PDF at scale $\ell$ and reference PDF from the middle of the inertial interval.

Figure 4.

(a) Solid lines represent 18 PDFs of $U_\parallel$ for the $1024^3$ implicit run; see full black square in (b). Red dashed lines present six PDFs of $U_\parallel$ for the Johns Hopkins University (JHU) $8192^3$ NS simulation; see full red circles in the (b). (b) The $L_2$ -distance between the PDFs of (a) and the reference NSS PDF from the middle of the inertial interval.

We further test the universality of inertial-interval HS with respect to the dissipation mechanism. For this purpose, we superpose our results with the data obtained from the NS simulations of the JHU turbulence database at $N = 8192^3$ (Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008; Yeung et al. Reference Yeung, Donzis and Sreenivasan2012, Reference Yeung, Zhai and Sreenivasan2015). The six indistinguishable dashed red lines in figure 4(a) are the PDFs of $U_\parallel$ for different scales from the inertial interval; see the full red dots in panel (b). The close-to-perfect superposition of the PDFs for the NSS and NS simulations confirms the universality of inertial-interval HS. Red circles in figure 4(b) quantify this collapse by measuring the $L_2$ -distance between the NS PDF at scale $\ell$ and the same reference PDF from the NSS data. It is remarkable that the NS inertial interval is shorter despite a much higher resolution.