Transitions of turbulent superstructures in generalized Kolmogorov flow

Self-organized large-scale flow structures occur in a wide range of turbulent flows. Yet, their emergence, dynamics, and interplay with small-scale turbulence are not well understood. Here, we investigate such self-organized turbulent superstructures in three-dimensional turbulent Kolmogorov flow with large-scale drag. Through extensive simulations, we uncover their low-dimensional dynamics featuring transitions between several stable and meta-stable large-scale structures as a function of the damping parameter. The main dissipation mechanism for the turbulent superstructures is the generation of small-scale turbulence, whose local structure depends strongly on the large-scale flow. Our results elucidate the generic emergence and low-dimensional dynamics of large-scale flow structures in fully developed turbulence and reveal a strong coupling of large- and small-scale flow features.

Despite their generic occurrence in a broad range of prototypical flows, the emergence and dynamics of turbulent superstructures, and in particular their interplay with small-scale turbulence, are currently not well understood. However, these aspects are crucial for developing a low-dimensional description of the large-scale dynamics of fully developed turbulent flows. Moreover, Landau pointed out in a famous remark on Kolmogorov's 1941 phenomenology [21] that the large-scale structure of a flow impacts the temporal variation of the rate of energy dissipation, potentially precluding a universal statistical theory of turbulence valid for all flows. Therefore, a characterization of the coupling of large and small scales has far-reaching implications for the assessment of universal features (or lack thereof) of small-scale turbulence [22].
The need to simultaneously resolve the slowly evolving large scales and the rapidly fluctuating small scales renders this problem challenging, both for experiments as well as for simulations. To address this, we here study a prototypical shear flow-a generalized three-dimensional (3D) turbulent Kolmogorov flow, which allows the investigation of turbulent superstructures without the complications imposed by boundaries. Traditional Kolmogorov flow is a simple, two-dimensional (2D) shear flow driven by a single Fourier mode, originally proposed by Kolmogorov to study the onset of turbulence on a periodic domain [23]. Indeed, the onset of turbulence [24] (for 3D see [25]) as well as the chaotic dynamics significantly above the onset [26][27][28] have been studied in detail. In contrast to the simple two-dimensional setting, the direct energy cascade toward smaller scales in three dimensions generates fully developed small-scale turbulence. Computational studies of turbulent 3D Kolmogorov flow revealed pronounced spatio-temporal large-scale intermittency and inhomogeneity [29,30], anisotropy [31], as well as translational symmetry breaking as a function of domain size [32].
Here, we generalize 3D Kolmogorov flow by including a large-scale drag term, which allows the manipulation of the range of scales on which turbulent superstructures occur. We reveal their low-dimensional dynamics in fully developed turbulence encompassing millions to billions of degrees of freedom by means of extensive simulations. In particular, we observe transitions between several largescale structures as a function of the damping parameter. Remarkably, we find that the flow is most effectively driven for strong damping, when the large-scale flow is shaped to resonate with the shear forcing. The generation of small-scale turbulence acts as main dissipation channel for the large-scale flow. A detailed spatiotemporal analysis of the energy transfer between large and small scales reveals a strong coupling of large-and small-scale flow features even in fully turbulent flows. In particular, we find that the mean rate of energy dissipation is well correlated with the energy injection rate by the large-scale forcing with a time lag associated to the energy transfer across scales. This establishes general-ized Kolmogorov flow as a prototypical incarnation of a flow in which the large-scale flow variation impacts the mean rate of kinetic energy dissipation as conjectured by Landau [21].
Our Kolmogorov flow is governed by the incompressible Navier-Stokes equation on a periodic 6π × 2π × 2π domain. Here, u denotes the incompressible velocity field (∇ · u = 0), and p the kinematic pressure. The shear force f = A sin(k f y)e x forces the flow on a single mode, corresponding to the smallest spanwise wave number k f = 1, with A = 1/2 (code units). The linear damping term, controlled by the parameter µ, only affects the sharp-Fourier-filtered field u, which only contains spatial scales larger than the forcing scale. This permits the manipulation of large-scale flow structures through the damping parameter, while allowing for freely evolving small scales. The flow is also subject to viscous dissipation with the kinematic viscosity ν.
To investigate the dynamics of turbulent superstructures, we conduct an extensive series of pseudo-spectral direct numerical simulations (DNS). We confirmed the robustness of our results for Reynolds numbers up to Re ≈ 1.8 × 10 4 [with Re = U L/ν based on the forcing scales L = 2π/k f and U = (LA) 1/2 ; Supplemental Material (SM)]. In the following, we focus on the intermediate Reynolds number of Re ≈ 7180. Important features of the flow can be assessed from the 2D dynamics in the forcing plane, which we obtain by decomposing the velocity field into its z-averaged part and fluctuations, u = v + we z + u . Here, v denotes the z-averaged, 2D velocity field in the forcing plane, and w is the z-averaged vertical (out-of-plane) velocity component, which is expected to be generally small; u denotes the 3D velocity fluctuations. The averaged 2D dynamics follows (see also [33]): where · z denotes z-averaging. Apart from interactions mediated by the 3D fluctuations through the stresses, from which only horizontal (in-plane) velocity fluctuations contribute to this equation, the vertical velocity w does not contribute to the dynamics of the 2D velocity v.
From our simulations, we observe the emergence of large-scale vortex structures in the forcing plane, which we characterize by the vorticity of the 2D averaged velocity, Ω = ∂ x v y − ∂ y v x . Figure 1 shows a range of different large-scale states along with visualizations of the full 3D vorticity field for three different values of µ. For small µ, a single large-scale vortex pair emerges, which can be perceived as a strongly inhomogeneous agglomeration of small-scale vorticity. The emergence of this single vortex pair can be understood as a long-wavelength instability of a flow which is initially proportional to the forcing (k x = 0) and then destabilizes, leading to a spanwise velocity on the smallest possible wave number (k x = k f /3). Interestingly, this wave number is also expected based on a linear stability analysis of laminar 2D Kolmogorov flow [23,34]. By increasing µ, this mode can be stabilized, triggering different flow patterns. For large µ, the large-scale flow predominantly organizes into three vortex pairs. The small-scale turbulence also becomes more homogeneous, resulting in a reduced vorticity flatness when compared to the single-vortex-pair state, see Fig. 1. Since the turbulent Taylor-scale Reynolds number increases from small non-zero values of µ to larger ones, this trend is not due to an increase of small-scale turbulence (SM). Overall, this suggests that different superstructures can induce different levels of small-scale intermittency. For intermediate values of µ, the flow dynamically switches between the three different states (one, two, and three vortex pairs). Therefore, the flow cannot be characterized by a single, universal large-scale flow state for a fixed set of parameters in this range. The spontaneous symmetry breaking in the form of the occurrence of one to three vortex pairs is a feature of the aspect-ratio-three domain. In domains of aspect ratio one [29,30,32] the mean velocity field is proportional to the forcing and therefore invariant with respect to continuous x translations. Larger-aspect-ratio domains, on the other hand, allow for an even greater variety of turbulent superstructures.
To characterize the large-scale dynamics, and in par- ticular the switching between large-scale states, we compute the stream function ψ, given by ∆ψ = −Ω. From its Fourier transformψ(k x , k y ), we compute the relative intensity of the three dominant streamwise Fourier Since I 1 + I 2 + I 3 = 1, the large-scale dynamics takes place in a triangle, whose corners correspond to the pure one, two and three vortex-pair states. By computing a mean streamwise wave number k s = k f (I 1 + 2I 2 + 3I 3 )/3 and averaging the solution over intervals of time where k s (t) is approximately constant, we can identify the well-defined large-scale states shown in Fig. 1. Figure 2 shows the large-scale dynamics in the intensity plane along with the mean streamwise wave number of the large-scale structures as a function of time. For low values of the damping parameter, the large-scale flow organizes into a singlevortex pair, which from time to time destabilizes due to strong, self-induced spanwise flow (see streamlines in . For large values of µT 1.8, we find that only the three-vortex-pair state remains accessible to the dynamics, with bursts destroying the pattern only for it to be reformed with a possible streamwise offset (Supplemental Video 6).
To uncover the role of small-scale turbulence for the large-scale dynamics, we investigate the energy budget of the flow. The total energy E(t) = 1 2 u · u xyz is a sum of the contributions from the z-averaged 2D velocity field E 2D (t) = 1 2 v · v xy , from the z-averaged vertical velocity E w (t) = 1 2 w 2 xy , and from the remaining 3D fluctuations E (t) = 1 2 u · u xyz . The energy of the 2D flow evolves according to Here, ε f (t) = v · f xy denotes the energy input into the flow by the shear forcing, ε µ (t) = µ v · v xy is the dissipation by large-scale damping, and ε ν (t) = ν ∇v : ∇v xy the viscous dissipation. The term Π(t) = − ∇v : u u z xy is the work performed by the 2D flow on the 3D stresses and characterizes the energy transfer between the 2D flow and the 3D turbulent fluctuations [35]. Figure 3(a) shows the time-averaged contributions to the kinetic energy as a function of the damping parameter µ. The total energy splits up into comparable contributions from the 2D energy and the energy in 3D turbulent fluctuations, with a negligible amount of energy contained in the average vertical velocity.
Remarkably, the total energy in the flow increases as we increase µ. At first, this appears counterintuitive given that we increase the strength of a dissipative term. We can explain this effect through the individual terms in the 2D budget which are shown in Fig. 3(b). The first important observation is that the energy input into the flow ε f increases with the damping parameter. Because ε f (t) = v · f xy , the energy input is maximal if v f , i.e. if the 2D flow is a simple shear flow in streamwise direction. Out of the three large-scale states, the threevortex-pair state is closest to a simple shear flow, as seen by comparing the streamlines in Fig. 1 with the forcing. Tuning the flow to this state by increasing the damping therefore maximizes the energy input. Figure 3(b) also shows that the energy input into the 2D flow is almost completely balanced by the energy transfer from the 2D flow to the 3D fluctuations. The generation of small-scale turbulence therefore is the main dissipation channel for the 2D flow. The dissipation by the large-scale damping is small for the entire parameter range, showing that its primary effect is that of shaping the large-scale flow. Because the 2D flow is predominantly large-scale, also viscous dissipation is negligible. We have also tested the robustness of these results across Reynolds numbers (Supplemental Figure 6).
To characterize the interplay of the large-scale flow and small-scale turbulence in more detail, we analyze how the inhomogeneity of the mean flow affects the spatial structure of the local energy transfer term Π(t, x, y) = −∇v : u u z , see Fig. 4. The explicit argument (t, x, y) indicates the coordinates over which no averaging was performed. The time-averaged energy transfer term is predominantly positive throughout the plane, as expected, and the large-scale state is clearly reflected in its spatial structure, with maxima close to saddle points of the mean flow. As this term injects energy into the 3D fluctuations, it leaves a footprint in the distribution of fluctuation energy E (t, x, y) = 1 2 u · u z . Further analysis confirms that, for all three states and across all values of µ, correlation coefficients of at least 60% and up to 90% are reached for the two quantities. We quantify the inhomogeneity of the 3D fluctuations by the relative standard deviation of the energy ρ = E t − E txy  of ρ are close for the two-vortex-pair and three-vortexpair states, but significantly higher for the one-vortexpair state. This is consistent with the vorticity flatness values given above, further highlighting the connection between the large-scale state and small-scale statistics (Supplemental Figures 3-5 and 7). We now turn to the temporal evolution of the energy injection, transfer and dissipation rates. Complementing the energy balance of v, the energy budget of the spatially averaged 3D fluctuations takes the form where ξ ν (t) = ν ∇u : ∇u xyz denotes viscous dissipation of the 3D fluctuations and ξ w (t) = u z u z · ∇w xy denotes the energy transfer to the mean vertical flow. We find ξ w to be two orders of magnitude smaller than Π and ξ ν , thus we only focus on the interplay of ε f , Π and ξ ν in the following. The temporal evolution of the self-organized largescale flow allows us to investigate the delays between large-and small-scale flow dynamics, which are otherwise only accessible through external temporal modulation [36,37]. Figure 5 shows the time traces of the energy injection into the large-scale flow, the energy transfer term, and the viscous dissipation of the 3D fluctuations for a representative µT = 1.1 at Re ≈ 7180 along with the normalized cross-correlations. Strong oscillations are evident, as well as a fairly strong correlation between the different time traces.
Cross-correlations between the energy transfer term and the large-scale energy injection, C(Π, ε f ; τ ), as well as the cross-correlations between the dissipation of 3D fluctuations and the energy transfer term, C(ξ ν , Π; τ ), reach values of about 90% for a temporal offset of approximately 0.6T ; ξ ν (t) is also strongly correlated with ε f (t), with a temporal offset of approximately 1.2T . The inset of Fig. 5 shows the temporal offsets between the different signals as they change with the damping parameter µ. For most values of µ the offset between ε f (t) and ξ ν (t) approximately equals the sum of the other two offsets. This can be explained by the fact that energy is transferred from the forcing scale to the dissipation scale by a direct cascade, in which the transfer between the 2D flow and the 3D fluctuations is an intermediate step. Remarkably, the two offsets are very similar in size. These observations do not hold for very small values of µ, where we also find significantly lower peak correlations between the injection rate and the dissipation rate (ca. 40-50% rather than 80-90%).
Overall, our results draw the following picture: turbulent superstructures in Kolmogorov flow consist of accumulations of small-scale vortices, which emerge through an instability of the largest streamwise scales accessible to the dynamics. These scales can be controlled by large-scale damping. As a function of the damping parameter, we find a rich transition scenario between several large-scale states with single, two and three vortex pairs. For intermediate damping, the system dynamically switches between various superstructures, illustrating that the system parameters do not uniquely determine the large-scale flow state. The generation of 3D turbulent fluctuations is the main dissipation channel for the large-scale flow, dominating over viscous dissipation and dissipation by damping. The 2D mean flow and the 3D fluctuations are strongly coupled through the energy transfer: the topology of the 2D mean flow determines the spatial distribution of the turbulent energy transfer rate, leading to inhomogeneities in the kinetic energy of the 3D fluctuations. Additional analysis shows that these inhomogeneities likely persists even at much higher Reynolds numbers. Generalized Kolmogorov flow is therefore one striking example of a flow in which the rates of energy transfer and ultimately dissipation vary considerably in response to the self-organization of large-scale superstructures, which has implications for the spatio-temporal structure of small-scale turbulence. Since a coupling of large-and small-scale flow features has been widely observed (see e.g. [1,3,38,39]), the scenario outlined by Landau is presumably relevant for a large class of flows. Furthermore, our results suggest the possibility of a low-dimensional description of large-scale superstructures in fully developed turbulent flows. Since the presented mechanisms governing their emergence and dynamics are quite generic, we expect our results to be relevant for a broad range of flows including geophysical and astrophysical flows.
Supplemental Table I. Summary of direct numerical simulation series. Simulations were run on real-space grids with 3N ×N ×N points spanning the 3D periodic box 6π × 2π × 2π between times t0 = 0 and t1. The forcing amplitude A and wave number k f are used to define large-scale units L = 2π/k f , T = (L/A) 1/2 , U = L/T and the Reynolds number Re = U L/ν. To complement the Reynolds number based on the forcing scales, we also compute the Taylor-scale Reynolds number R λ = 2 E t 5/(3ν ε t) using the time-averaged energy of the velocity fluctuations and the time-averaged total viscous dissipation rate. The time step is determined by the Courant-Friedrichs-Lewy condition, and CFL refers to the Courant number. The small-scale resolution of the simulations is characterized through the product of the largest resolved wave number kM and the Kolmogorov length scale η; τη denotes the Kolmogorov time scale.

Supplemental videos 1 and 2 -evolution of the 3D fluctuations
To illustrate the structure of the 3D turbulence in generalized Kolmogorov flow and its dependence on the large-scale flow state (see also Fig. 1 from the main text), we provide Supplemental Videos 1 and 2 for Re ≈ 4180. Dedicated direct numerical simulations were used to generate the required 4D datasets, with N = 256 and k M η values of 2 and 1.9 for µT = 0 and µT = 3.5, respectively.
The videos show intense vortex filaments, visualized through volume renderings of the vorticity field ω = ∇ × u, see also Supplemental Fig. S1. In particular, the Supplemental Figure and the videos show: • the z-averaged z component of the vorticity (xy plane, "back"); • the x-averaged x component of the vorticity (yz plane, "left"); • the y-averaged y component of the vorticity (zx plane, "bottom"); the colormap corresponds to the one in Fig. 1(a)-(c) from the main text; • the volume-rendered vorticity magnitude (colormap is the black-and-white version of the colormap used for Fig. 1(d)-(f) from the main text); • the xyz coordinate system (bottom left corner), with x, y, z components colored red, green and blue, respectively.
The visualizations were generated with VTK [40].
Supplemental Figure S1. Evolution of the 3D fluctuations for generalized Kolmogorov flow: µT = 0 (left panel, Supplemental Video 1) and µT = 3.5 (right panel, Supplemental Video 2). Color maps are comparable to those used for Fig. 1 from the main text.
Supplemental videos 3-6 -evolution of the mean 2D flow The dynamics of the 2D flow v can be illustrated with animations of the z-averaged z component of the vorticity Ω, which we describe here. Supplemental Fig. S2 shows snapshots of the corresponding Supplemental Videos 3-6.
Supplemental Video 3 shows how the vortex pair in the µT = 0 case may become disorganized but always reforms. Supplemental Videos 4 and 5 show the large-scale dynamics for intermediate values of µ with switches between several large-scale states. Supplemental Video 6 shows that for the comparatively large µT = 3.5 the system only visits the three-vortex-pair state.  Fig. 1(a)-(c) from the main text, the visualizations show the z-averaged out-of-plane vorticity along with streamlines of the 2D flow, averaged over a time window in which the large-scale pattern is persistent. As described in the main text, we identify the large-scale state by the mean streamwise wave number. The visualizations show that there is also a dependence of the flow field on µ, as can be seen from differences in the streamlines of the flow field as the damping parameter is changed. Because temporal averages may only be computed over finite intervals, the figures do not always exhibit the expected symmetry (e.g. the different vortices sometimes have different shapes and/or intensities). The length of the temporal averaging interval is determined by the duration for which the individual large-scale structures exist. Because the structures may reform at a shifted streamwise position after a breakdown, averages computed over distinct temporal intervals cannot be merged directly.
The three figures have a tabular structure (with the damping parameter changing over the vertical and the largescale state over the horizontal). Missing panels signify that a particular large-scale state has not been clearly observed for a given value of µ in our data.

Energy budget for various Reynolds numbers
As explained in the main text, we decompose the generalized Kolmogorov flow velocity into a z-averaged, purely 2D component v, the z-averaged z component w and the 3D fluctuations u . The energy distribution among these components, as well as the contributions to the 2D energy budget, Eq. (3) from the main text, do not vary significantly with the Reynolds number, as can be seen in Supplemental Fig. S6. We show here the contributions to the total energy and the 2D energy budget terms (i.e. variations of Fig. 3 from the main text) for Re ≈ 2850 and Re ≈ 18100. The results are qualitatively the same.

Vorticity statistics for various Reynolds numbers
To characterize how the small-scale statistics depend on the damping parameter µ, we present vorticity statistics in Supplemental Fig. S7. As expected, the PDFs become more and more heavy-tailed as the Reynolds number (based on the forcing scales) is increased, in particular the flatness increases with Re. Additionally, we observe a clear trend: At a fixed Reynolds number, the vorticity PDFs exhibit the heaviest tails for small values of µ, and significantly less heavy tails for larger values of µ, despite the fact that the Taylor-scale Reynolds number increases from small non-zero values of µ to larger values of µ (cf. Supplemental Table I). This is also reflected in the corresponding vorticity flatness. We expect this trend to be the result of the pronounced spatio-temporal inhomogeneity of the large-scale flow. As discussed in the main text, for µT = 0 the flow is strongly inhomogeneous in space and time due to the emergence and decay of a single vortex pair. For small values of µ, the vortex pair is stabilized, giving rise to a temporally persistent, spatially inhomogeneous large-scale structure. For larger values of µ, two and three vortex pairs emerge (as seen in Supplemental Figs. S3-S5), which render the large-scale flow less inhomogeneous.