1. Introduction
In developing statistical theories of turbulence, a certain degree of universality is often assumed. However, even the precise values of supposed universal constants such as the Kolmogorov constant for the energy spectrum are still up for debate (Sreenivasan Reference Sreenivasan1995; Yeung & Zhou Reference Yeung and Zhou1997; Donzis & Sreenivasan Reference Donzis and Sreenivasan2010; Ishihara et al. Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2016). In part, this may be because very different types of flows (including those from simulations) are attempted to be described under one umbrella. While it is agreed that the large scales of turbulence differ between flow types, there is still the expectation that the small scales of turbulence approach some universal statistical state at large Reynolds numbers, as first hypothesized by Kolmogorov in his classical theory of turbulence (hereafter K41, see Kolmogorov (Reference Kolmogorov1941)). In his theory, the central role is played by the mean dissipation rate 
$\langle \varepsilon \rangle$ that determines the average energy flux from large scales to small scales. If we assume that the value of this average energy flux is the only information that propagates through the cascade dynamics from large to small scales, then the small-scale statistics have to take the universal form dictated by the K41 theory. However, the idealized assumptions of the K41 theory are not quite correct. In particular, it turned out that there is a build-up of fluctuations around the average flux of energy towards smaller scales (Cerutti & Meneveau Reference Cerutti and Meneveau1998; Buzzicotti et al. Reference Buzzicotti, Linkmann, Aluie, Biferale, Brasseur and Meneveau2018; Yasuda & Vassilicos Reference Yasuda and Vassilicos2018), leading to extreme fluctuations at the smallest scales. This is known as small-scale intermittency (Sreenivasan & Antonia Reference Sreenivasan and Antonia1997). While this observation contradicts the K41 theory, it does still allow for universality of small-scale statistics for a given Reynolds number.
Complementarily, a limitation of applicability of the K41 theory was brought up shortly after publication by Landau, who remarked that sufficiently slow (or sufficiently spatially extended) flow variations cannot be expected to be filtered out by the cascade dynamics and will lead to alterations of small-scale statistics (Landau & Lifshitz Reference Landau and Lifshitz1987; Kraichnan Reference Kraichnan1974; Frisch Reference Frisch1995). This is now known as large-scale intermittency (see Chien, Blum & Voth (Reference Chien, Blum and Voth2013) for a discussion of the terminology). A way to rationalize his criticism is to assert that turbulent flows have finite spatio-temporal correlations. Therefore, the small-scale statistics at a certain point in space and time can depend on the flow only within a certain spatio-temporal region. Any variations on larger scales will generate dissimilar, statistically independent domains of the flow. Aggregated statistics over these domains will necessarily depend on the large-scale variations. Contrary to small-scale intermittency, the large-scale variations are expected to differ from flow to flow.
Technically, such large-scale variations could be considered in conflict with Kolmogorov's assumption of local isotropy (Kolmogorov Reference Kolmogorov1941), implying that K41 theory does not apply here. While one could therefore exclude such very large-scale variations from the discussion of universality (Kraichnan Reference Kraichnan1974), it is important to note that they are not uncommon in real-world turbulent flows: For example, large-scale variations in atmospheric turbulence were observed to enhance the effect of small-scale intermittency (Böttcher et al. Reference Böttcher, Renner, Waldl and Peinke2003; Muschinski, Frehlich & Balsley Reference Muschinski, Frehlich and Balsley2004; Böttcher, Barth & Peinke Reference Böttcher, Barth and Peinke2007; Feraco et al. Reference Feraco2021). Similarly, many experimental set-ups for turbulence produce large-scale fluctuations (Mouri et al. Reference Mouri, Takaoka, Hori and Kawashima2006; Mouri, Hori & Takaoka Reference Mouri, Hori and Takaoka2009; Blum et al. Reference Blum, Bewley, Bodenschatz, Gibert, Gylfason, Mydlarski, Voth, Xu and Yeung2011), such as von Kármán flow (Voth et al. Reference Voth, La Porta, Crawford, Alexander and Bodenschatz2002; Mordant, Lévêque & Pinton Reference Mordant, Lévêque and Pinton2004b; Lawson & Dawson Reference Lawson and Dawson2014). Certain types of simulated flows can also have emerging large-scale variations, such as Kolmogorov flow (Borue & Orszag Reference Borue and Orszag1996; Goto & Vassilicos Reference Goto and Vassilicos2016; Lalescu & Wilczek Reference Lalescu and Wilczek2021), or they can feature shifting turbulent/non-turbulent interfaces such as turbulent jet flows (Gauding et al. Reference Gauding, Bode, Brahami, Varea and Danaila2021a,Reference Gauding, Bode, Denker, Brahami, Danaila and Vareab). By contrast, many numerical turbulence simulations are designed to display a minimum of large-scale fluctuations, e.g. by enforcing a fixed total energy (Ishihara et al. Reference Ishihara, Kaneda, Yokokawa, Itakura and Uno2007). In practice, spatial or temporal averages are often performed to improve sample sizes, but this may hide the presence (or absence) of the large-scale fluctuations.
Kolmogorov (Reference Kolmogorov1962) and Oboukhov (Reference Oboukhov1962) partly addressed Landau's remark in their refined similarity hypotheses (see also Mouri et al. (Reference Mouri, Takaoka, Hori and Kawashima2006) and Chien, Blum & Voth (Reference Chien, Blum and Voth2013)). Instead of basing their theory on the mean energy dissipation rate 
$\langle \varepsilon \rangle$, they postulated universal statistics on a length scale 
$r$ as a function of the energy dissipation spatially averaged over scale 
$r$. This theory constituted an improvement over the previous one by incorporating small-scale intermittency, acknowledging the fact that energy dissipation fluctuates strongly in space and time (Meneveau & Sreenivasan Reference Meneveau and Sreenivasan1991). However, the theory assumed universal log-normal large-scale fluctuations and did not consider external variations in the sense of Landau, which may depend on the flow type.
In the literature, the impact of large-scale intermittency on statistics is not often discussed explicitly. There are some parts of the literature speaking about external intermittency as the random switching between turbulent and non-turbulent flow, which can affect inertial-range properties and specifically higher-order statistics (Kuznetsov, Praskovsky & Sabelnikov Reference Kuznetsov, Praskovsky and Sabelnikov1992; Mi & Antonia Reference Mi and Antonia2001; Gauding et al. Reference Gauding, Bode, Brahami, Varea and Danaila2021a,Reference Gauding, Bode, Denker, Brahami, Danaila and Vareab). Here, we consider fully developed turbulence, and we understand large-scale intermittency as fluctuations on scales comparable to or larger than the integral scales. In order to assess the impact of such large-scale fluctuations, Monin and Yaglom proposed a two-state model explaining how the mixing of flow regions with different dissipation rates changes the coefficients of structure functions while retaining their scaling (Monin & Yaglom Reference Monin and Yaglom1975; Davidson Reference Davidson2004; Chien et al. Reference Chien, Blum and Voth2013).
One approach to quantify the impact of large-scale fluctuations is to compute conditional statistics. It was found that in many cases, second-order structure functions conditional on the large-scale velocity indeed depend on these large scales (Praskovsky et al. Reference Praskovsky, Gledzer, Karyakin and Zhou1993; Sreenivasan & Stolovitzky Reference Sreenivasan and Stolovitzky1996; Sreenivasan & Dhruva Reference Sreenivasan and Dhruva1998; Blum et al. Reference Blum, Kunwar, Johnson and Voth2010, Reference Blum, Bewley, Bodenschatz, Gibert, Gylfason, Mydlarski, Voth, Xu and Yeung2011; Chien et al. Reference Chien, Blum and Voth2013). Some studies also considered conditional higher-order statistics in this context (Praskovsky et al. Reference Praskovsky, Gledzer, Karyakin and Zhou1993; Blum et al. Reference Blum, Kunwar, Johnson and Voth2010; Gauding et al. Reference Gauding, Bode, Brahami, Varea and Danaila2021a,Reference Gauding, Bode, Denker, Brahami, Danaila and Vareab). Böttcher et al. (Reference Böttcher, Renner, Waldl and Peinke2003, Reference Böttcher, Barth and Peinke2007) found that the increment statistics of atmospheric turbulence matches better with laboratory data when conditioning on the mean wind speed and proposed an ensemble model based on this idea. Shnapp (Reference Shnapp2021) found that also the Lagrangian statistics in canopy flows are affected simultaneously by small-scale and large-scale intermittency.
Complementarily, introducing artificial variations of the large scales allows one to study their impact more systematically. Knebel, Kittel & Peinke (Reference Knebel, Kittel and Peinke2011) used an active grid to generate statistics resembling intermittent atmospheric data. Chien et al. (Reference Chien, Blum and Voth2013) studied how the small-scale statistics of a flow between oscillating grids changes with different degrees of large-scale variations. We here adopt a similar approach by introducing artificial variations into direct numerical simulations. Given that the non-trivial behaviour of higher-order statistics is a hallmark feature of turbulence, we focus on the effect of large-scale intermittency specifically on higher-order statistics at the example of the flatness and propose modelling strategies.
In a parallel line of research, temporal variations of the energy input (akin to the ones that we analyse in the present work) have been used to study the response and resonance behaviour of turbulence. Given periodic driving, sometimes a resonant response of the flow can be observed (Cadot, Titon & Bonn Reference Cadot, Titon and Bonn2003; von der Heydt, Grossmann & Lohse Reference von der Heydt, Grossmann and Lohse2003a,Reference von der Heydt, Grossmann and Lohseb; Kuczaj, Geurts & Lohse Reference Kuczaj, Geurts and Lohse2006; Kuczaj et al. Reference Kuczaj, Geurts, Lohse and van de Water2008; Cekli, Tipton & van de Water Reference Cekli, Tipton and van de Water2010). Some studies considered periodically kicked turbulence (Lohse Reference Lohse2000; Jin & Xia Reference Jin and Xia2008). It was also found that periodic driving can affect transport properties (Jin & Xia Reference Jin and Xia2008; Bos & Rubinstein Reference Bos and Rubinstein2017; Yang Y. et al. Reference Yang, Chahine, Rubinstein and Bos2019; Yang R. et al. Reference Yang, Chong, Wang, Verzicco, Shishkina and Lohse2020). When driving turbulence with periodically varying shear, there is a critical frequency above which the turbulent flow can no longer be sustained (Yu & Girimaji Reference Yu and Girimaji2006; Hamlington & Dahm Reference Hamlington and Dahm2009). Furthermore, studying variations of the energy input can provide insight into the workings of the energy cascade (Bos, Clark & Rubinstein Reference Bos, Clark and Rubinstein2007; Fang & Bos Reference Fang and Bos2023). While we here study a similar set-up with periodic driving, we differ from this line of research by taking the modulated driving as a model for large-scale intermittency and investigating its impact on time-aggregated, higher-order statistics compared to simulations without the modulation.
In this study, we ask how much and by which mechanisms the higher-order flow statistics are affected by slow fluctuations at the large scales. Then we ask whether these effects can be explained by models based on ensembles of statistically stationary flows. To address these questions, we introduce sinusoidal variations of the energy injection rate into direct numerical simulations of turbulence. While these variations are to be understood as a model for the large-scale fluctuations that can occur in various forms in any flow, we focus on sinusoidal injection-rate signals, which allow for a systematic investigation using periodic averaging. For different frequencies of the oscillations, we analyse the statistics of velocity, vorticity and Lagrangian acceleration. We observe an amplification of flatness, which we attribute to three different effects. The first and largest contribution comes from mixing distributions with varying width and can simply be estimated from the time series of mean energy and dissipation rate. The second contribution results from variations of higher-order flow statistics and can be described by an ensemble of stationary flows. A third contribution occurs at certain frequencies of the input, where we observe that time-averaged statistics are not simply a mix of statistics corresponding to the various injection rates. Instead, the second moment of the dissipation rate displays stronger excursions than expected. By considering the time series of both mean and variance of the dissipation rate, we construct an ensemble of stationary flows that can capture all three effects and thus accurately predict flatness values. Finally, we also explore how our results are affected by the Reynolds number.
2. Direct numerical simulations
We analyse data from direct numerical simulations (DNS) of the forced, incompressible Navier–Stokes equation in the vorticity formulation
\begin{equation} \partial_t \boldsymbol{\omega} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\omega} = \boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} + \nu\,\nabla^2 \boldsymbol{\omega} + \boldsymbol{f} \end{equation}
solved by a pseudo-spectral code (Lalescu et al. Reference Lalescu, Bramas, Rampp and Wilczek2022). Here, 
$\boldsymbol {\omega }(\boldsymbol {x}, t)$ is the vorticity field on a periodic domain 
$\boldsymbol {x} \in [0, 2{\rm \pi} ]^3$ (discretized on 
$512^3$ grid points), 
$\boldsymbol {u}(\boldsymbol {x}, t)$ is the velocity field with 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u} = 0$, 
$\nu$ is the kinematic viscosity and 
$\boldsymbol {f}(\boldsymbol {x}, t)$ is the forcing. The flow is forced on the large scales by amplifying a discrete band of Fourier modes 
$k \in [1.4, 2.3]$ (DNS units), imposing a prescribed energy injection rate 
$\xi (t)$ (specified below). Due to the forcing, the simulations are approximately statistically homogeneous and isotropic. The spatial resolution varies between 
$2.4 \leqslant k_M\eta \leqslant 6.4$ across the simulations, where 
$k_M$ is the maximum wavenumber, and 
$\eta$ is the Kolmogorov length scale. Each simulation lasts between 
$131\,T_{int}$ and 
$262\,T_{int}$ (where 
$T_{int}$ is the integral time scale). In order for the statistics to be independent of the initial condition (which is a statistically stationary turbulent field), we use only data after a transient period of at least 
$8\, T_{int}$.
The main focus of this study is on unsteady simulations with periodically varying injection rate
\begin{equation} \xi(t) = \xi_0 + A_\xi \sin(2{\rm \pi} t/ P), \end{equation}
where 
$\xi _0$ is the base injection rate, 
$A_\xi = 0.95\xi _0$ is the amplitude of the oscillations, and the period time 
$P$ is varied from 
$0.5\, T_{int}$ up to 
$32.8\, T_{int}$. Throughout this paper, all reference values such as the integral scales and the Kolmogorov scales are taken from a simulation with constant injection rate 
$\xi (t) = \xi _0$ at Taylor-scale Reynolds number 
$R_{\lambda }^{ref} \approx 103$ (with the exception of the resolution criterion 
$k_M\eta$, for which 
$\eta$ is computed separately in each simulation).
In the course of this work, we will associate the statistics of the oscillating simulations with the statistics of simulations at various constant injection rates 
$\xi (t) = \alpha \xi _0$. After the transient, these simulations become approximately statistically stationary. An ensemble of 23 simulations on a 
$512^3$ grid with 
$\alpha$ linearly spaced between 
$0.05$ and 
$2.14$ (
$R_\lambda$ between 60 and 120) allows us to cover all of the injection rates reached by the oscillating simulations. An extension of three simulations on a 
$1024^3$ grid with 
$\alpha \in \{4, 9, 18\}$ (
$R_\lambda$ between 134 and 173, each lasting 
$41\, T_{int}$) will help us in §§ 4.4 and 4.5 to model features of high-Reynolds-number turbulence that we observe in the oscillating simulations. A separate set of higher-Reynolds-number simulations is described and analysed in § 5 in order to evaluate the Reynolds-number dependence of our results.
3. Results
To illustrate the flow response in the oscillating simulations, figure 1 shows time series of various flow characteristics. The energy and the dissipation rate resemble the injection rate signal with a delay. The amplitude of their oscillations reduces at smaller period times, indicating that the energy cascade acts as a low-pass filter (Cadot et al. Reference Cadot, Titon and Bonn2003; Kuczaj et al. Reference Kuczaj, Geurts and Lohse2006; Bos et al. Reference Bos, Clark and Rubinstein2007). This filtering effect is analysed in more detail in figure 2, where we show the relative amplitude and the time delay of the variations of the dissipation rate compared to the injection rate. Only for the slowest variations of the injection rate (largest 
$P$) does the dissipation rate oscillate at the same amplitude. As soon as the variation time scales become comparable to the integral time of the flow, the amplitude drops. The decrease is consistent with the 
$P^1$ and 
$P^3$ scalings reported in the literature (Kuczaj et al. Reference Kuczaj, Geurts and Lohse2006; Bos et al. Reference Bos, Clark and Rubinstein2007), where the transition threshold between these two regimes depends on the Reynolds number. Similarly, the time delay saturates at a well-defined cascade time of approximately 
$2.4\, T_{int}$ only for the slowest oscillations. Due to the method that we use to measure the time delay 
$\tau _c$, it cannot be larger than the period time 
$P$. Instead, for small period times, we observe that it follows 
$P/2$ approximately, which corresponds to a phase shift of 
${\rm \pi}$ (compare Kuczaj et al. Reference Kuczaj, Geurts and Lohse2006; Bos et al. Reference Bos, Clark and Rubinstein2007).
Time series of different flow characteristics in simulations with oscillating energy injection rate at period times 
$P$. All quantities are normalized by their average value from the reference simulation with injection rate 
$\xi (t) = \xi _0$, average energy 
$E^{ref}$, and Reynolds number 
$R_{\lambda }^{{ref}} \approx 103$. The time-dependent Taylor-scale Reynolds number is computed as 
$R_\lambda (t) = E(t)\,(20/(3\nu \varepsilon (t)))^{1/2}$ and thus peaks when energy is large compared to dissipation rate. The fact that the amplitude of the flow response decreases with the oscillation frequency shows that the flow acts as a low-pass filter. The code and post-processed data used to generate this figure can be explored at https://www.cambridge.org/S0022112024007006/JFM-Notebooks/files/Figure-1.ipynb.
Amplitude 
$A_\varepsilon$ and time delay 
$\tau _c$ of the variations of dissipation rate, determined by finding the maximum point of the cross-correlation between injection rate and dissipation rate, 
$A_\xi A_\varepsilon /2 = \max _\tau \overline {\langle \xi (t-\tau ) \varepsilon (t) \rangle } - \xi _0^2$ and 
$\tau _c = \mathrm {argmax}_\tau \overline {\langle \xi (t-\tau ) \varepsilon (t) \rangle }$ with 
$0 \leqslant \tau \leqslant P$. Note that these formulas give the exact amplitude and phase shift for sinusoidal signals. The code and post-processed data used to generate this figure can be explored at https://www.cambridge.org/S0022112024007006/JFM-Notebooks/files/Figure-2.ipynb.
3.1. Periodically averaged statistics
We investigate three exemplary single-point flow statistics: components of velocity, vorticity and Lagrangian acceleration, denoted by 
$u$, 
$\omega$ and 
$a$, respectively. Since we are specifically interested in the effect of the periodically oscillating energy injection rate, we introduce a periodic average 
$\langle \cdot \rangle$ that averages over all degrees of freedom apart from the phase of the oscillation. All other degrees of freedom, namely the spatial dependence, the three vector components, and full periods of the oscillations, are expected to be statistically equivalent and thus averaged over. Hence a periodically averaged quantity 
$\langle \cdot \rangle$ depends only on the phase 
$0\leqslant t\ {\rm mod}\ P < P$. We further increase the number of samples by aggregating samples of similar phase in linearly spaced temporal bins 
$0 = t_0 < \dots < t_n = P$ such that the 
$i$th bin is given by all samples with
\begin{equation}t_{i-1} \leqslant t\ {\rm mod}\, P < t_i.\end{equation}
We choose a suitable number of bins 
$n$ (as indicated throughout the paper) depending on the need. For short period times, we reduce 
$n$ such that temporal bins are not shorter than 
$\tau _\eta /2$ (where 
$\tau _\eta$ is the Kolmogorov time scale). The corresponding phase-dependent probability density functions (p.d.f.s) based on the periodic average are shown for 
$P \approx 8.2\, T_{int}$ in figures 3(a–c). While velocity distributions are generally close to Gaussian, vorticity and acceleration distributions display heavier tails. Over the period of the oscillations, the width of the p.d.f.s varies strongly. With the artificial oscillations, we imitate the effect of the large-scale fluctuations that may occur in various types of turbulent flows. If such temporal variations are simply averaged over, then this may change flow statistics compared to flows without such variations. For this study, we compare time-averaged statistics of oscillating flows, such as the distributions obtained as the mean of the oscillating p.d.f.s in figures 3(a–c) (grey dashed lines), with statistics from the reference flow at stationary injection rate 
$\xi (t) = \xi _0$ (grey dotted lines). After characterizing the impact of the temporal variations, we pose the question of whether it can be modelled statistically as a mixture of an ensemble of stationary flows.
Periodically averaged statistics of the oscillating flow with period time 
$P \approx 8.2\,T_{int}$. (a–c) The p.d.f.s of velocity, vorticity and acceleration. Here, 
$u_{rms}$, 
$\omega _{rms}$ and 
$a_{rms}$ denote the root mean square of the velocity, vorticity and acceleration, respectively, all computed from the reference simulation with 
$\xi (t) = \xi _0$. A directed circle indicates the arrow of time. (d) Visualization of the colour code for the temporal bins. For the computation of the p.d.f.s, the oscillation period is split into 
$n=6$ equally spaced bins (see top of image). For the spectra, 
$n=8$ temporal bins are used (see bottom of image). The injection rate (grey) is shown for reference. (e) Compensated instantaneous energy spectra. In the course of each period, the energy is first injected at small wavenumbers, then transported to small scales, and then dissipated. The code and post-processed data used to generate this figure can be explored at https://www.cambridge.org/S0022112024007006/JFM-Notebooks/files/Figure-3.ipynb. Animated versions of this figure can be found in the supplementary material (available at https://doi.org/10.1017/jfm.2024.700) as movies 1, 2 and 3 for the simulations with 
$P \approx 2.1\,T_{int}$, 
$8.2\,T_{int}$ and 
$32.8\,T_{int}$, respectively. Movie 4 shows a three-dimensional visualization of the energy density in the simulation with 
$P \approx 8.2\,T_{int}$.
Another quantity that serves to illustrate the different states of the oscillating flows is the compensated energy spectrum shown in figure 3(e). Over the first half of an oscillation period, energy is injected at small wavenumbers 
$k$. During this time, the compensated energy spectrum peaks at those small wavenumbers. In the second half of the oscillation period, the energy cascades towards larger wavenumbers while the energy injection is low. During this time, the compensated energy spectrum peaks at large wavenumbers and decays.
3.2. Impact of oscillations on time-averaged statistics
In order to quantify the probability of extreme events, we consider the total (i.e. time-aggregated) (
$k$th-order hyper-)flatness, defined as
\begin{equation} F_{X,k}^{tot} = \frac{\overline{\langle X(t)^k \rangle}}{\overline{\langle X(t)^2 \rangle}^{k/2}} \end{equation}
for the three quantities 
$X = u, \omega, a$ (components of velocity, vorticity and acceleration) and even 
$k\geqslant 4$. Note that all three quantities have zero mean. The overbar denotes temporal averaging:
\begin{equation} \overline{f(t)} = \frac{1}{P} \int_0^P {\rm d}t\,f(t). \end{equation}
Since the periodic average 
$\langle \cdot \rangle$ depends only on the phase of the oscillation 
$t\ {\rm mod}\ P$, the temporal average ranges only from 
$0$ to 
$P$. For simplicity, we do not write the dependence on space and on vector component indices, which are averaged over by the periodic ensemble average 
$\langle \cdot \rangle$ (described above). We will consider the cases of regular flatness (or kurtosis) with 
$k=4$ and the hyper-flatness with 
$k=6$. Since the high-order moment puts emphasis on the tails of distributions, the flatness is large whenever the tails have comparably high probability. Since the flatness is non-dimensional, it does not change when linearly rescaling the p.d.f., e.g. Gaussian distributions always have regular flatness 3, and 6th-order hyper-flatness 15.
The flatness and hyper-flatness are shown in figure 4 for varying period times (solid black line with dots). When changing the period time, the system transitions between two asymptotic behaviours (black dotted and dashed lines). At short period times, the oscillations are too fast for the flow to respond. Only the average energy input is relevant, which is the same as in the statistically stationary reference flow. Flatness values are the same as in this reference flow (black dashed line). For long period times, the system approaches the quasi-stationary behaviour. Since the injection rate varies over time scales much longer than the dynamics of the flow, the flow has time to equilibrate to the instantaneous injection rate. Therefore, the instantaneous statistics are the same as given by a statistically stationary flow at the same injection rate. Time-aggregated statistics are then given by a superposition of the instantaneous flow statistics, weighted by the distribution of injection rates. This construction, which we call the ‘injection ensemble’ (black dotted line), is detailed in § 4.3. Finally, the transient regime between the two limits depends more crucially on the dynamical response of the flow to the changes in the forcing.
Flatness and hyper-flatness of velocity, vorticity and acceleration at different period times of the oscillations (black dots and solid line). For small period times, the flatness converges to the one of the stationary reference simulation (black dashed line). For large period times, it converges to the injection ensemble (black dotted line), which superposes stationary turbulence statistics weighted according to the distribution of injection rates. The coloured lines show different approximations of the flatness increase (based on 
$12\leqslant n \leqslant 256$ temporal bins per period). The error bars are standard errors obtained by assuming that moments computed from large temporal subintervals are i.i.d. Gaussian. These subintervals were chosen to be multiples of the period time of length at least 
$8\, T_{int}$. The code and post-processed data used to generate this figure can be explored at https://www.cambridge.org/S0022112024007006/JFM-Notebooks/files/Figure-4.ipynb.
In the next section, we show that temporal mixing always increases flatness, with the main effect coming from mixing distributions with varying variance. In order to fully capture the transition of flatness, we then employ an ensemble approach. In a first step, we include the low-pass filter effect (see discussion of figure 2 above). This leads to the ‘mean dissipation ensemble’ (discussed in § 4.3), which correctly captures both limiting behaviours. As we detail in §§ 4.4 and 4.5, by additionally modelling the instantaneous flow states more precisely (‘fluctuating dissipation ensemble’), we achieve good quantitative agreement in capturing the transition.
4. Analysis and modelling
4.1. How temporal mixing increases flatness
In order to understand the impact of temporal mixing, consider the instantaneous (
$k$th-order hyper-)flatness
\begin{equation} F^{inst}_{X,k}(t) = \frac{\langle X(t)^k \rangle}{\langle X(t) ^2 \rangle^{k/2}}. \end{equation}
We will compare this instantaneous flatness to the total flatness 
$F_{X,k}^{tot}$. The difference is that the instantaneous flatness is computed for each phase of the oscillation, while for the total flatness, the averages are computed over all samples, aggregated over the phase of the oscillations. It is common practice to consider such time-averaged statistics in order to obtain large sample sizes.
Figure 5 compares the instantaneous flatness (black solid line) with the total flatness (red solid line) for the simulation at 
$P \approx 8.2\,T_{int}$. For the velocity, the instantaneous flatness does not vary notably within a period of the oscillations, staying slightly below the Gaussian value 3 (figure 5a). Such sub-Gaussian velocity statistics are in line with more recent observations (Jiménez Reference Jiménez1998; Gotoh, Fukayama & Nakano Reference Gotoh, Fukayama and Nakano2002; Mouri et al. Reference Mouri, Takaoka, Hori and Kawashima2002; Wilczek, Daitche & Friedrich Reference Wilczek, Daitche and Friedrich2011). The total flatness, however, takes a value above 3, purely as a result of the temporal mixing. Our results suggest that the difference between slightly sub-Gaussian and super-Gaussian velocity statistics can be attributed to large-scale flow variations. Mouri et al. (Reference Mouri, Takaoka, Hori and Kawashima2002) observed a transition from sub-Gaussian to super-Gaussian velocity statistics along decaying grid turbulence (accompanied by an increase of velocity gradient flatness), while the Reynolds number decreased. Our present observations corroborate their hypothesis that such behaviour may be explained by an increasing degree of large-scale fluctuations further away from the grid.
Instantaneous flatness with 
$k=4$ (a–c, black solid) of velocity, vorticity and acceleration for the simulation at 
$P \approx 8.2\,T_{int}$. A grey solid line indicates the flatness of the reference simulation. The temporal weight (d–f, blue solid) is used to compute the weighted average (a–c, blue dashed), already resulting in a flatness value larger than the reference value of the stationary flow (grey solid). Additionally, the multiplier 
$M_{X,4}$ raises the flatness significantly (black arrows and red solid line). Note that the values depend slightly on the number of temporal bins, which here is 
$n=64$. The code and post-processed data used to generate this figure can be explored at https://www.cambridge.org/S0022112024007006/JFM-Notebooks/files/Figure-5.ipynb.
For our vorticity and acceleration data, we observe notable variations of flatness in the course of a period (figures 5b,c). As in the case of velocity, the total flatness surpasses any of the instantaneous values. This shows that the presence of large-scale fluctuations can significantly impact the flatness of these quantities, too. In the literature, acceleration flatness values scatter considerably, even for the same Reynolds numbers (Vedula & Yeung Reference Vedula and Yeung1999; La Porta et al. Reference La Porta, Voth, Crawford, Alexander and Bodenschatz2001; Voth et al. Reference Voth, La Porta, Crawford, Alexander and Bodenschatz2002; Mordant, Crawford & Bodenschatz Reference Mordant, Crawford and Bodenschatz2004a; Bec et al. Reference Bec, Biferale, Boffetta, Celani, Cencini, Lanotte, Musacchio and Toschi2006; Yeung et al. Reference Yeung, Pope, Lamorgese and Donzis2006; Ishihara et al. Reference Ishihara, Kaneda, Yokokawa, Itakura and Uno2007; Bentkamp, Lalescu & Wilczek Reference Bentkamp, Lalescu and Wilczek2019). (Vedula & Yeung (Reference Vedula and Yeung1999) report pressure gradient flatness.) A possible explanation for this may be the differing degree of large-scale variations across simulations and experiments. In our simulations, we see that large-scale variations can indeed raise the acceleration flatness by roughly a factor of 2.
How does the additional averaging operation lead to this increase in flatness? We can rewrite the total flatness (3.2) as a function of the instantaneous flatness (4.1) as
\begin{align} F_{X,k}^{tot} &= \frac{\overline{F^{inst}_{X,k}(t)\,\langle X(t)^2 \rangle^{k/2}}}{\overline{\langle X(t)^2 \rangle}^{k/2}} \end{align}
\begin{align} &= M_{X,k}\,\overline{w_{X,k}(t)\,F^{inst}_{X,k}(t)}, \end{align}with the temporal weight
\begin{equation} w_{X,k}(t) = \frac{\langle X(t)^2 \rangle^{k/2}}{\overline{\langle X(t)^2 \rangle^{k/2}}} \geqslant 0, \quad \overline{w_{X,k}(t)} = 1, \end{equation}and the multiplier
\begin{equation} M_{X,k} = \frac{\overline{\langle X(t)^2 \rangle^{k/2}}}{\overline{\langle X(t)^2 \rangle}^{k/2}} \geqslant 1. \end{equation}
This means that the total flatness is a weighted time average of the instantaneous flatness, multiplied by a number 
$M_{X,k} \geqslant 1$. The inequality is a result of Jensen's inequality for powers, 
$\overline {f(t)^p} \geqslant \overline {f(t)}^p$, being valid for 
$p\geqslant 1$ and any non-negative time-dependent variable 
$f(t)$. Note that instead of time-dependent quantities, the statements about the flatness can be phrased analogously for any kind of conditional statistics such as considered in Yeung et al. (Reference Yeung, Pope, Lamorgese and Donzis2006). Furthermore, similar expressions appear in superstatistical modelling (e.g. Beck, Cohen & Swinney Reference Beck, Cohen and Swinney2005).
Remarkably, this shows that temporal mixing always increases flatness. The only way to arrive at low flatness would be to have 
$w_{X,k}(t)$ emphasize periods of time where the instantaneous flatness is low, or to have no variations at all. In any case, an absolute lower bound is given by 
$F_{X,k}^{tot} \geqslant \min _t F^{inst}_{X,k}(t)$. In our examples, however, the flatness lies significantly above this lower bound, even far above the value for the stationary flow. Based on (4.3), the increase can be decomposed into contributions from the multiplier and from the weighted average, as illustrated in figure 5.
In the remainder of this study, we investigate in more detail the mechanisms that lead to the total change in flatness as seen in figures 4 and 5. We discern three different effects, each leading to an amplification of flatness, which we capture using increasingly complex models:
(i) Variations of variance. Over time, the instantaneous variance of all considered quantities fluctuates. This automatically leads to an increase in flatness measured by the multiplier
$M_{X,k}$. In the case of self-similar distributions, the instantaneous flatness becomes independent of time. As a result, (4.3) reduces to $F_{X,k}^{tot} = M_{X,k}F^{inst}_{X,k}$. So if the instantaneous distributions are self-similar (which is approximately the case for the velocity), then the temporal variations of the variance are the only effect leading to changes of flatness. Similarly, for vorticity and acceleration, variations of variance are a large contribution to flatness increase (see black arrows in figure 5).(ii) Mixing of non-self-similar statistics. If the instantaneous flatness
$F^{inst}_{X,k}(t)$ changes over time, then the weighted average in (4.3) will impact the flatness, too. To first approximation for slow variations, we can understand this as mixing of flow states that essentially look like stationary turbulence at different injection rates. So in order to understand the impact of the weighted average, we have to measure flatness and variance of stationary flows and superpose them in appropriate ways (we will call these ‘injection ensemble’ and ‘mean dissipation ensemble’).(iii) Dynamical effects. Finally, for injection rate variations on time scales comparable to flow scales, there may be dynamical effects altering the instantaneous flatness values. For period times of approximately
$6$–$12$ integral times, we indeed observe such deviations from simple mixing, leading to even stronger instantaneous fluctuations and thus higher flatness. This will be modelled using the ‘fluctuating dissipation ensemble’.
In the following, we will go through these effects in more detail.
4.2. Variations of variance
The conceptually simplest contribution to the flatness increase is given by the multiplier 
$M_{X,k}$, generated by variations of variance. It is a dominant contribution for all quantities that we consider. Since the velocity p.d.f. is in good approximation self-similar over time, its total flatness can be computed as
\begin{equation} F_{u,k}^{tot} \approx M_{u,k} F_{u,k}^{ref} = \frac{\overline{\langle u^2(t) \rangle^{k/2}}}{\overline{\langle u^2(t) \rangle}^{k/2}}\,F_{u,k}^{ref}= \frac{\overline{E(t)^{k/2}}}{\overline{E(t)}^{k/2}}\,F_{u,k}^{ref} \end{equation}
from the periodically averaged time series of mean energy 
$E(t) = 3\langle u^2(t) \rangle /2$ and the flatness 
$F_{u,k}^{ref}$ of the stationary reference simulation. Here, we have assumed that the total velocity flatness of the stationary reference flow is the same as the (almost constant) instantaneous flatness in the oscillating flows. Figure 4 shows that this multiplier indeed fully explains the increase of velocity flatness and hyper-flatness (blue line).
For vorticity and acceleration, the multiplier can serve as a first estimate of the flatness increase. Given the time series of the mean dissipation rate 
$\langle \varepsilon (t) \rangle$, the instantaneous vorticity variance is given exactly by 
$\langle \omega (t)^2 \rangle = \langle \varepsilon (t) \rangle/(3\nu)$, and the acceleration variance can be estimated by the Heisenberg–Yaglom prediction, 
$\langle a(t)^2 \rangle \approx a_0 \langle \varepsilon (t) \rangle ^{3/2}\nu ^{-1/2}$ (Monin & Yaglom Reference Monin and Yaglom1975). Approximating the weighted average in (4.3) with the value from the stationary flow, we then have
\begin{equation} F_{\omega,k}^{{tot}} \approx M_{\omega,k} F_{\omega,k}^{ref} = \frac{\overline{\langle \varepsilon(t) \rangle^{k/2}}}{\overline{\langle \varepsilon(t) \rangle}^{k/2}}\, F_{\omega,k}^{ref} \end{equation}and
\begin{equation} F_{a,k}^{{tot}} \approx M_{a,k} F_{a,k}^{ref} \approx M_{a,k}^{HY} F_{a,k}^{ref} = \frac{\overline{\langle \varepsilon(t) \rangle^{3k/4}}}{\overline{\langle \varepsilon(t) \rangle^{3/2}}^{k/2}}\,F_{a,k}^{ref}. \end{equation}
Here, 
$F_{\omega,k}^{ref}$ and 
$F_{a,k}^{ref}$ denote the flatness values of the stationary reference simulation. In the case of the acceleration, the multiplier is further approximated as 
$M_{a,k}^{HY}$ using the Heisenberg–Yaglom prediction. These estimates of the flatness capture a large part of the amplification but underestimate it (figure 4, blue lines). Nevertheless, they are a simple means to evaluate whether differences in flatness between two flows of comparable Reynolds number could be explained by large-scale intermittency.
For the sinusoidal modulations used in this study, the long-period limit of the vorticity and acceleration multipliers can even be computed analytically. In this limit, the periodically averaged dissipation rate 
$\langle \varepsilon (t) \rangle$ converges to the sinusoidal form of the injection rate 
$\xi (t)$. Therefore, we can write
\begin{align} M_{\omega, k} \xrightarrow{P\to\infty} M_{\omega, k}^\infty &= \frac{\overline{\xi(t)^{k/2}}}{\overline{\xi(t)}^{k/2}} = \frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}}{\rm d}t\,(1 + \alpha \sin t)^{k/2} \end{align}
\begin{align} &= 1 + 2\sum_{m=1}^{k/4} \binom{k/2}{2m} \binom{2m-1}{m} \left(\frac{\alpha}{2}\right)^{2m}, \end{align}
where 
$\alpha =0.95$ is the relative amplitude of the modulations of the injection rate and 
$\left(\begin{smallmatrix}\cdot \\ \cdot\end{smallmatrix}\right)$ denotes binomial coefficients. The sum is taken over all integers 
$1\leqslant m \leqslant k/4$. For the flatness, this amounts to 
$M_{\omega, 4}^\infty = 1+{\alpha ^2}/{2} \approx 1.45$, and for the 6th-order hyper-flatness, we get 
$M_{\omega, 6}^\infty = 1 + {3\alpha ^2}/{2} \approx 2.35$. Analogously, for the acceleration, evaluation of the integrals for the Heisenberg–Yaglom estimate, (4.8), yields 
$M_{a,4}^{HY} \xrightarrow {} M_{a,4}^{{HY},\infty } \approx 1.69$ and 
$M_{a,6}^{HY} \xrightarrow {} M_{a,6}^{{HY},\infty } \approx 3.23$. These factors, multiplied with the value for the reference flow, are shown as blue dotted lines in figure 4 and agree well with the data.
4.3. Mixing of non-self-similar statistics
In the previous subsection, we found that the main contribution to the flatness amplification comes from the multiplier. In order to get a more accurate estimate, however, we have to include effects of the weighted average in (4.3). In order to understand how this average leads to an increase in flatness, we here want to approximate the temporal mixing by the mixing of an ensemble of stationary turbulent flows. This is motivated by the limit of infinitely slow variations: if the variations of the injection rate are slower than any dynamics in the flow, then the flow will be statistically quasi-stationary, continually relaxing to a state of stationary turbulence at the current injection rate.
To assess this idea, we compare the statistics to an ensemble of 23 simulations with stationary energy injection rates (described in § 2). Instead of mixing statistics in the temporal domain, we will mix statistics of the various statistically stationary flows. Figure 6 shows how the ensemble can be weighted so that its statistics match those of the various oscillating simulations.
Visualization of how to construct the injection and mean dissipation ensembles. Left-hand plot: injection rate time series (black dashed line) and periodically averaged dissipation rate time series (coloured solid lines) for the oscillating flows from 
$P\approx 33\,T_{int}$ (violet) to 
$P\approx 3\,T_{int}$ (yellow), using 
$69 \leqslant n \leqslant 256$ temporal bins. Due to the low-pass filtering of the energy cascade (see discussion of figure 2), the amplitude of the dissipation rate variations decreases with decreasing period times. Centre plots: histograms of those periodically averaged curves. These histograms are used to construct the injection ensemble (leftmost histogram) and the mean dissipation ensembles (other histograms). Dotted lines indicate the distribution (4.11) as a comparison, with rescaled 
$A_\xi$ to match the variance of the dissipation rate. Right-hand plot: injection rates (black dashed lines) and dissipation rates (coloured solid lines) of the statistically stationary ensemble members as functions of time. The code and post-processed data used to generate this figure can be explored at https://www.cambridge.org/S0022112024007006/JFM-Notebooks/files/Figure-6.ipynb.
The simplest way to do this is to match the distribution of injection rates. The sinusoidal injection rate (2.2) takes values distributed according to the p.d.f.
\begin{equation} f(\xi) = \begin{cases}
\dfrac{1}{\rm \pi} [A_\xi^2 - \left(\xi -
\xi_0\right)^2]^{{-}1/2} & \text{if } |\xi - \xi_0| <
A_\xi, \\[8pt] 0 & \text{otherwise}.
\end{cases} \end{equation}
Weighting the ensemble members by this distribution gives rise to what we call the ‘injection ensemble’. As can be seen in figure 4 (dotted black line), it accurately predicts the flatness in the limit of long period times. Since the distribution of injection rates (4.11) does not depend on the period time, this ensemble cannot be used to capture the transition.
In order to capture the transition of vorticity and acceleration flatness, we introduce the ‘mean dissipation ensemble’. Instead of matching the distribution of injection rates, we here match the distribution of dissipation rates. To this end, we compute the mean dissipation rate for each temporal bin, then the distribution of these dissipation rates over the oscillation period. As shown in figure 6, this leads to narrower ensemble weightings at short period times and to wider ensemble weightings at long period times due to the low-pass filter effect of the energy cascade. In the limit 
$P\to \infty$, the periodically averaged dissipation rate equals the injection rate, and the mean dissipation ensemble becomes identical to the injection ensemble described above (leftmost histogram in figure 6). Accordingly, the flatness of the mean dissipation ensemble transitions between the value of stationary turbulence at short period times and the value of the injection ensemble at long period times (see figure 4, orange line). However, while giving the correct limits and better results than the pure multipliers, the mean dissipation ensemble still misses some of the flatness increase in the transition regime.
The mean dissipation ensemble is very much in the spirit of the K41 theory. It relies on the assumption that the instantaneous small-scale statistics in the form of vorticity and acceleration flatness are tied to the instantaneous value of the mean dissipation rate. The fact that it fails to predict flatness values correctly indicates that the small scales depend on more than just the average dissipation rate. This will be addressed in the next subsection, where we show that the combined values of dissipation mean and variance are needed to match statistics of the oscillating flows with ensembles of stationary turbulence.
4.4. Dynamical effects
For period times of the order of 6–12 integral times, figure 4 still shows notable deviations between vorticity and acceleration flatness of the mean dissipation ensemble compared to the oscillating flows. According to (4.2), this must be due to a mismatch in instantaneous variance or flatness. Based on how we constructed the ensemble, it means that the instantaneous statistics of vorticity and acceleration are not uniquely related to the instantaneous mean dissipation rate when comparing between the oscillating flows and stationary turbulence.
This is investigated more closely in figure 7, showing how the various quantities are related by parametric plots, both in the oscillating case (yellow to violet) and in the stationary case (red with dots, each dot representing one simulation). In addition to the vorticity and acceleration flatness, we show the normalized second moment of dissipation, 
$\langle \varepsilon (t)^2 \rangle /\langle \varepsilon (t) \rangle ^2$, as a measure of the fluctuations around the mean dissipation rate (see also Gauding et al. Reference Gauding, Bode, Denker, Brahami, Danaila and Varea2021b). None of these three quantities appears to be determined uniquely by the instantaneous mean dissipation rate 
$\langle \varepsilon (t)\rangle$ (figures 7a,c,e). For very slow, quasi-static changes of the injection rate (infinite period time), we would expect the oscillating curves (yellow to violet) to collapse onto the one for stationary turbulence (red with dots). However, indications of this become visible only at the longest period times. Instead, for all curves we observe a larger flatness on their intensifying branch (increasing dissipation rate) than on their decaying branch (decreasing dissipation rate), most pronouncedly at intermediate period times such as 
$P \approx 11.6\,T_{int}$. This indicates that intensifying turbulence is characterized by heavier-tailed dissipation, vorticity and acceleration statistics than its stationary equivalent.
Parametric plots of instantaneous statistics of the oscillating flows (moments periodically filtered with Gaussian filters of width 
$P/64$), compared with statistics of the statistically stationary flows (red line with dots). (a,c,e) For the oscillating simulations, flatness values of vorticity and acceleration as well as the normalized second moment of dissipation behave similarly as a function of the mean dissipation. (b,d) All flatness curves approximately collapse when plotting them as a function of the normalized second moment of dissipation. Note that we have added simulations of stationary turbulence at 
$\xi (t) = \{4, 9, 18\}\xi _0$ to extend the red lines (larger dots). The code and post-processed data used to generate this figure can be explored at https://www.cambridge.org/S0022112024007006/JFM-Notebooks/files/Figure-7.ipynb.
However, plotting vorticity and acceleration flatness against 
$\langle \varepsilon (t)^2 \rangle /\langle \varepsilon (t) \rangle ^2$ (figures 7b,d) shows that the excursions seen in figures 7(a,c,e) share some similarities. Curves for different period times approximately collapse, which means that vorticity and acceleration flatness and normalized second moment of dissipation scale in the same way, independent of the oscillation frequency. In fact, they even behave similarly to the scaling of stationary turbulence (red line with dots). Note that in order for the red lines to reach this far, the base ensemble was extended by three stationary turbulence simulations on a 
$1024^3$ grid with significantly higher injection rates (larger red dots; for more details see § 2).
4.5. Fluctuating dissipation ensemble
We observed that the instantaneous flows display features of higher-Reynolds-number turbulence, manifesting in large values of the instantaneous flatness and the normalized second moment of dissipation. Based on this, we can construct a refined ensemble model, which we call ‘fluctuating dissipation ensemble’. For each instantaneous value of 
$\theta (t) = \langle \varepsilon (t)^2 \rangle /\langle \varepsilon (t) \rangle ^2$, we select a simulation of stationary turbulence that matches this value most closely (see figure 8). Then we know, based on the approximate collapse in figures 7(b,d), that the selected stationary flow will approximate the instantaneous values of vorticity and acceleration flatness, too. However, we know that total flatness is determined not only by instantaneous flatness but also by instantaneous variance (see (4.3)). Hence we do not expect this ensemble model to predict the total flatness well, yet.
Visualization of the procedure to construct the fluctuating dissipation ensemble for the oscillating flow at 
$P \approx 8.2\,T_{int}$. The green and red curves are the same as in figure 7(e). For each instantaneous value 
$\theta (t)$ of the oscillating flow, the stationary flow (red dots) with the closest value of 
$\theta$ is selected. Then it is rescaled at constant Reynolds number and viscosity, such that also the mean dissipation rate 
$\mu$ is matched (grey arrows). The superposition of the statistics of all the rescaled stationary flows selected in this way then constitutes the fluctuating dissipation ensemble. The code and post-processed data used to generate this figure can be explored at https://www.cambridge.org/S0022112024007006/JFM-Notebooks/files/Figure-8.ipynb.
In order to extend the set of simulations of stationary turbulence that we have, we perform Reynolds-similar rescalings of the data. This is justified by the fact that every solution of the Navier–Stokes equation such as our simulation data can be rescaled in space and time and will still be a solution of the Navier–Stokes equation (Frisch Reference Frisch1995). The rescaling factors in space, 
$\alpha _S$, and in time, 
$\alpha _T$, can be chosen freely. In the rescaled system, each quantity is given by its original value multiplied with the rescaling factors according to its units, e.g. the rescaled dissipation rate is given by
\begin{equation} \varepsilon' = \frac{\alpha_S^2}{\alpha_T^3}\varepsilon. \end{equation}
Dimensionless quantities such as the flatness and the Reynolds number remain constant. This method allows us to generate a family of different flows from a single simulation run, connected by Reynolds similarity. Since we have two degrees of freedom for the rescaling, we can decide to keep the viscosity (units 
$\mathrm {Length}^2\ \mathrm {Time}^{-1}$) constant by choosing 
$\alpha _S^2 = \alpha _T \equiv \alpha ^2$. The data that we generate by this method are the same as what would be obtained by running a simulation at the same Reynolds number and viscosity but with different forcing scale and injection rate.
Let us now use this method to construct our advanced ensemble model. For the mean dissipation ensemble, we matched the distribution of mean dissipation rate 
$\mu (t)=\langle \varepsilon (t)\rangle$ between the oscillating flow and the ensemble. Now we are matching the distribution of 
$\theta (t) = \langle \varepsilon (t)^2 \rangle /\langle \varepsilon (t) \rangle ^2$. However, in order to get a good quantitative match between ensemble model and oscillating flow, it would be better to match both quantities at the same time, i.e. the full p.d.f. of 
$\mu$ and 
$\theta$,
\begin{equation} f(\mu, \theta) = \overline{\delta\left(\mu - \langle \varepsilon(t) \rangle\right) \delta\left(\theta - \frac{\langle \varepsilon(t)^2 \rangle}{\langle \varepsilon(t) \rangle^2} \right)}. \end{equation}
Using the rescaling approach, this is done in the following way: For a given combination of 
$\mu$ and 
$\theta$, we select the stationary flow that matches the value of 
$\theta$ most closely (see figure 8). In general, the mean dissipation rate of the 
$\theta$-matched simulation 
$\langle \varepsilon \rangle _\theta$ will not be equal to 
$\mu$. We can, however, choose a rescaling factor 
$\alpha$ such that
\begin{equation} \mu = \frac{\alpha_S^2}{\alpha_T^3} \langle\,\varepsilon \rangle_\theta = \alpha^{{-}4}\,\langle \varepsilon \rangle_\theta. \end{equation}
After the rescaling, both 
$\mu$ and 
$\theta$ and thus both mean and variance of the dissipation rate are matched (figure 8, grey arrows). Based on the observations in this study, we can hope that this will allow us to approximate both lower- and higher-order instantaneous statistics of vorticity and acceleration, too.
Finally, we superpose the statistics of the rescaled statistically stationary flows to imitate the time-averaged data of the oscillating flow. For example, for the acceleration p.d.f., which rescales with 
$\alpha _S/\alpha _T^2= \alpha ^{-3}$, we have
\begin{equation} f^{ens}(a) = \int {\rm d}\mu \int {\rm d}\theta \, f(\mu, \theta)\,\alpha(\mu, \theta)^{{-}3}\,f^{stat}_\theta(a \alpha(\mu, \theta)^3). \end{equation}
Here, 
$f^{ens}$ denotes the acceleration p.d.f. of the ensemble, and 
$f^{stat}_\theta$ denotes the acceleration p.d.f. (before rescaling) of the stationary flow with the given value of 
$\theta$. The rescaling factor 
$\alpha (\mu, \theta )$ can be computed from 
$\mu$ and 
$\theta$ by (4.14). Any other flow statistics can be computed likewise. In particular, figure 4 shows vorticity and acceleration flatness and hyper-flatness values of this ensemble (red lines).
The fluctuating dissipation ensemble excellently captures vorticity flatness and 6th-order hyper-flatness (see figure 4). In particular, it successfully reproduces the increase in flatness due to dynamical effects that was missing from the mean dissipation ensemble. Although we incorporated only the first two moments of dissipation in the formulation of the model (which are related to second and fourth moments of vorticity), the 6th-order hyper-flatness of vorticity is matched perfectly, too. This indicates that a simple measure of dissipation fluctuations such as 
$\langle \varepsilon (t)^2 \rangle /\langle \varepsilon (t) \rangle ^2$ may be sufficient for accurate modelling of even higher-order small-scale statistics. For acceleration statistics, flatness is captured quite well, while there are deviations for the hyper-flatness. Some of the inaccuracies of the model can already be surmised from the slight mismatch of flatness between oscillating and stationary flows in figures 7(b,d). Overall, the ensemble models are well suited to capture amplification of higher-order moments due to large-scale intermittency. When the large-scale fluctuations are happening on time scales comparable to a few integral scales of the flow, it becomes crucial to include also the second moment of the dissipation rate into the modelling.
5. Comparison to higher Reynolds number
Achieving good convergence for the high-order moments of the oscillating flows requires very long simulation times, making the full analysis including the ensemble modelling feasible only for simulations on 
$512^3$ grid points. In order to evaluate the Reynolds number dependence of our results, we additionally simulated six oscillating flows on a 
$1024^3$ grid, accompanied by a reference flow with constant injection rate at 
$R_\lambda^{\textit{ref}} \approx 227$ and spatial resolution 
$k_M \eta \approx 2.0$. The duration of the simulations is between 
$44\,T_{int}$ and 
$88\,T_{int}$ (with 
$T_{int}$ measured in the reference flow). The simulation parameters are the same as described for the main data set in § 2, only at smaller viscosity. The period times 
$P$ are varied between 
$3.1\,T_{int}$ and 
$17.5\,T_{int}$. The spatial resolution of the oscillating flows varies between 
$1.7 \leqslant k_M\eta \leqslant 3.1$. Due to the computational constraints, we did not consider any ensemble of stationary flows in this case.
Figure 9(a) shows the periodically averaged energy spectrum of the flow at 
$P\approx 8.8\,T_{int}$, which can be compared directly to that shown in figure 3. While the overall dynamics looks very similar, at the higher Reynolds number, we can observe a more pronounced inertial range. Figures 9(b–g) are analogous to figure 4. They show how the flatness and hyper-flatness increase as a function of the period time of the oscillations. As a direct comparison, we show the main data from figure 4 as a grey line in the background. In order to make them comparable, we divided those data by the values from the reference flow of the main data set and multiplied them by the values of the reference flow from the higher-Reynolds-number data set.
(a) Periodically averaged compensated energy spectra for the oscillating flow with period time 
$P\approx 8.8\,T_{int}$ in the higher Reynolds number data set. The bins are the same as used in figure 3. (b–g) Flatness and hyper-flatness of velocity, vorticity and acceleration as functions of period time for the oscillating flows in the higher Reynolds number data set (black dots and solid line). Error bars are computed as in figure 4. The reference flow with stationary injection rate at 
$R_\lambda^{\textit{ref}} \approx 227$ (black dashed line) corresponds to the limit of short period times. Multipliers are shown as blue lines as in figure 4, computed using 
$147 \leqslant n \leqslant 256$ temporal bins per period. Ensemble models are not available for this data set. The flatness of the oscillating flows from the main data set (as in figure 4, black dots and solid line) is shown in grey in the background, multiplied by the ratio of the flatness values of the respective reference flows. The code and post-processed data used to generate this figure can be explored at https://www.cambridge.org/S0022112024007006/JFM-Notebooks/files/Figure-9.ipynb.
Given that we compensate the data only by the respective integral time scales 
$T_{int}$ and reference flatness values 
$F_{X,k}^{ref}$, the increase in flatness is remarkably similar between the two Reynolds numbers, within the error bars of our data. In particular, the factor by which the flatness increases is very close. This can be explained by the multiplier mechanism detailed in § 4.2. There, we showed that we can compute exactly the infinite-period-time limit of the vorticity multiplier and of the Heisenberg–Yaglom estimate of the acceleration multiplier, and we found that they are independent of Reynolds number. As in figure 4, these are shown as blue dotted lines. Since the multipliers account for a large part of the (hyper-)flatness increase, it is expected that we observe a comparable increase at different Reynolds numbers. While the remaining effect of the weighted average in (4.3) may have non-trivial Reynolds number dependence, these are good reasons to believe that the effect of large-scale intermittency will be of similar magnitude at even higher Reynolds numbers.
Furthermore, it is very interesting that the critical period time of the transition (in units of 
$T_{int}$) is almost identical at this higher Reynolds number. As described in the discussion of figure 2, we found a low-pass filter effect of turbulence, which dampens the effect of large-scale intermittency at shorter period times. The fact that we observe a collapse at two different Reynolds numbers shows that it is mainly the integral flow scales determining the time scale of this filter.
We conclude this section by discussing some practical implications. Overall, the preceding discussion indicates that both the time scales of the transition (in integral units) and the magnitudes of the multipliers will remain comparable even at higher Reynolds number. This means that it is possible to estimate the magnitude of the effect of large-scale intermittency on higher-order moments by measuring the amplitude and time scales of the large-scale variations in a turbulent flow. If the frequency of the variations is comparable to or lower than 
${\sim }1/(5T_{int})$, then it is likely that the flatness measurements are affected by them. Moreover, the factor by which the flatness is likely increased due to large-scale intermittency can be estimated by a computation such as the one detailed in § 4.2.
6. Conclusion
The idea of universality in turbulence relies on the assertion that the smallest scales of the flow are the result of its internal nonlinear dynamics and therefore become independent of the external forcing mechanism that generates the turbulence. Next to the Reynolds number, it was put forward by Landau (Landau & Lifshitz Reference Landau and Lifshitz1987) that sufficiently slow, large-scale variations will constrain universality. In order to ensure comparability across different turbulent flows, it is therefore important to identify and characterize the properties of the flow that affect small-scale universality. Here, we evaluated the relevance of Landau's remark for the statistics of velocity, vorticity and acceleration.
We computed the flatness and 6th-order hyper-flatness of these quantities in direct numerical simulations of unsteady, homogeneous, isotropic turbulence subject to sinusoidal variations of the energy injection rate with different period times. Analytically, we found that temporal mixing always increases flatness. For the velocity, large-scale fluctuations can make the difference between sub-Gaussian and super-Gaussian statistics. For the acceleration, we observed a flatness increase by approximately a factor of 2.
We identified three different mechanisms that each lead to an increase in flatness. First, any variations of variance over time lead directly to an increase in flatness that can be measured by the multiplier 
$M_{X,k}$. This fully explains the increase of velocity flatness. The multipliers can be either computed directly or approximated from the time series of mean energy and dissipation rate. Second, through weighted averaging, flatness increases when large instantaneous flatness values coincide with large variance. Together, these two effects can be modelled by an appropriate superposition of statistics of stationary turbulence, identifying each instantaneous state of the oscillating flows with a flow of stationary turbulence based on the mean dissipation rate. The third effect occurs only for period times of the order of 6–12 integral times. In that regime, we observed that the fluctuating input generates stronger excursions of the instantaneous flatness than expected from the identification with stationary turbulence. We found that this is also reflected in the normalized second moment of dissipation, which can act as a good predictor for flatness of both vorticity and acceleration. Within the scope of the quantities considered here, we found that dissipation mean and variance together are a good representation of the instantaneous statistical state of the small scales, and that they can be used to construct an ensemble model that accurately captures both regular flatness and 6th-order hyper-flatness.
While we focused here on the flatness of velocity, vorticity and acceleration, we put forward ensemble descriptions that capture the full flow statistics at once. Since these models combine the data from simpler turbulent flows to model more complex turbulent flows, they require only assumptions about how to superpose these statistics. Meanwhile, they can make predictions about many different flow statistics. An interesting direction of future work will be to investigate whether similar ensembles of lower-Reynolds-number flows can be used to model higher-Reynolds-number turbulence.
We introduced here artificial variations of the energy injection rate into statistically homogeneous and isotropic flows and found that these variations had a strong effect on higher-order statistics. With this work, we hope to spark interest in further studies on large-scale intermittency, including investigations of the strength and impact of ‘naturally occurring’ large-scale fluctuations in experimental and other simulated flows. Quantifying their effect as a correction to the Reynolds number effects may be crucial to explain the observed scatter of higher-order statistics in the literature. This may be achieved by accompanying reports of higher-order statistics by quantities such as the multipliers used in our study or coefficients of inertial-range scaling laws as proposed by Chien et al. (Reference Chien, Blum and Voth2013). If such an analysis allows us to match data across different flow types more precisely, then this would make a strong point for a formulation of universality explicitly conditional on large-scale variations. Overall, we hope that this work can contribute to improve comparability in turbulence by better characterizing the effect of large-scale intermittency.
Supplementary material
Supplementary movies and computational notebook files are available at https://doi.org/10.1017/jfm.2024.700. The simulations were conducted using our code TurTLE (Lalescu et al. Reference Lalescu, Bramas, Rampp and Wilczek2022), which is available at https://gitlab.mpcdf.mpg.de/TurTLE/turtle. The post-processing code and the post-processed data can be accessed at https://www.cambridge.org/S0022112024007006/JFM-Notebooks along with interactive notebooks generating each of the figures in this study.
Acknowledgements
We would like to acknowledge interesting and useful discussions with M. Carbone and G.A. Voth, as well as technical support from C.C. Lalescu. We thank the anonymous reviewers for their constructive comments, which helped to improve the manuscript. GitHub Copilot was used during code development.
Funding
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 101001081). The authors gratefully acknowledge the scientific support and HPC resources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU) under the NHR project ‘EnSimTurb’. NHR funding is provided by federal and Bavarian state authorities. NHR@FAU hardware is partially funded by the German Research Foundation (DFG) – 440719683.
Declaration of interests
The authors report no conflict of interest.
Author contributions
M.W. designed the study. L.B. carried out the numerical simulations and analysis. Both authors analysed the data and wrote the manuscript.
Open access
