2 Physical set-up

In this section, we describe the set-up to be investigated. We consider the idealised case of forced incompressible three-dimensional flow in a triply periodic box of dimensions $L\times L\times H$ . The thin direction $H$ will be referred to as the vertical ‘ $z$ ’ direction and the remaining two as the horizontal ‘ $x$ ’ and ‘ $y$ ’ directions. The geometry of the domain is illustrated in figure 1. The flow obeys the incompressible Navier–Stokes equation:

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}P+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

$\boldsymbol{u}$

$P$

$\unicode[STIX]{x1D708}$

$\boldsymbol{f}$

$x$

$y$

$x$

$y$

$\boldsymbol{f}=(-\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D713},\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713},0)$

$\boldsymbol{f}$

$k_{f}\equiv 2\unicode[STIX]{x03C0}/\ell$

$\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{f}\rangle =\unicode[STIX]{x1D716}$

$\langle \cdot \rangle$

$\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}$

$-\unicode[STIX]{x1D708}_{n}(-\unicode[STIX]{x1D6FB}^{2})^{n}\boldsymbol{u}$

Figure 1. The domain used is a box of height $H$ and square base of side length $L$ . The forcing is invariant along the thin direction and stochastic with fixed mean rate of energy input, while involving only wavenumbers $\boldsymbol{k}$ with $|\boldsymbol{k}|=k_{f}=2\unicode[STIX]{x03C0}/\ell$ . The thin direction will be referred to as the vertical, the others as the horizontal directions.

The system (2.1b ) is characterised by three non-dimensional parameters: the Reynolds number based on the energy injection rate $Re=(\unicode[STIX]{x1D716}\ell ^{4})^{1/3}/\unicode[STIX]{x1D708}$ , the ratio between forcing scale and domain height $Q=\ell /H$ and the ratio between forcing scale and the horizontal domain size $K=\ell /L$ . The ratio between $K$ and $Q$ gives the aspect ratio $A=K/Q=H/L$ of the domain. The Kolmogorov dissipation length is denoted as $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D708}^{3/4}/\unicode[STIX]{x1D716}^{1/4}=\ell Re^{3/4}$ .

The simulations performed for this work used an adapted version of the Geophysical High-Order Suite for Turbulence (GHOST) which uses pseudo-spectral methods including $2/3$ aliasing to solve for the flow in the triply periodic domain, (see Mininni et al. Reference Mininni, Rosenberg, Reddy and Pouquet2011). The resolution was varied from $256^{2}\times 16$ grid points to $2048^{2}\times 128$ grid points depending on the values of the parameters. To explore the space spanned by these three parameters, we have performed systematic numerical experiments: for a fixed value of $Re$ and $K=1/8$ , different simulations are performed with $Q$ varying from small to large values. The runs are continued until a steady state is reached where all quantities fluctuate around their mean value. This is repeated for eight different values of $Re$ from $Re=203$ (resolution $256^{2}\times 16$ ) to $4062$ (resolution $2048^{2}\times 128$ ) and for one value of hyper-viscosity ( $n=8$ , $\unicode[STIX]{x1D708}_{8}=10^{-38}$ as in (Celani et al. Reference Celani, Musacchio and Vincenzi2010)), as a consistency check, since many of the previous studies of thin-layer turbulence used hyper-viscosity. For $Re=305$ , we also perform a run with $K=1/16$ ( $L\rightarrow 2L$ ). The number of runs performed for each $Re$ are summarised in table 1.

Table 1. Summary of the different runs performed. For each $Re$ and $K$ several runs for different values of $Q$ have been performed. The horizontal resolution is $n_{x},n_{y}$ , while $n_{z}$ stands for the vertical resolution at $Q=2$ . The vertical resolution was changed with $Q$ to maintain an isotropic grid, $Kn_{x}=Kn_{y}=Qn_{z}$ wherever possible.

To quantify the energy distribution among different scales it is convenient to work in Fourier space. The Fourier series expansion of the velocity reads

(2.2a,b ) $$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},t)=\mathop{\sum }_{\boldsymbol{k}}\hat{\boldsymbol{u}}_{\boldsymbol{k}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}},\quad \hat{\boldsymbol{u}}_{\boldsymbol{ k}}=\frac{1}{L^{2}H}\int \boldsymbol{u}(\boldsymbol{x},t)\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\,\text{d}\boldsymbol{x}\end{eqnarray}$$

where $\hat{\boldsymbol{u}}_{\boldsymbol{k}}=(\hat{u} _{\boldsymbol{k}}^{(x)},\hat{u} _{\boldsymbol{k}}^{(y)},\hat{u} _{\boldsymbol{k}}^{(z)})$ and the sum runs over all $\boldsymbol{k}\in ((2\unicode[STIX]{x03C0}/L)\mathbb{Z})^{2}\times (2\unicode[STIX]{x03C0}/H)\mathbb{Z}$ . In the pseudo-spectral calculations, this sum is truncated at a finite $k_{res}$ . Since flow in a thin layer is a highly anisotropic system, it is important to consider quantities in the vertical and horizontal directions separately. For this purpose, we monitor various quantities in our simulations: first of all, the total energy spectrum as a function of horizontal wavenumber

(2.3) $$\begin{eqnarray}E_{tot}(k_{h})=\frac{1}{2}\mathop{\sum }_{{\boldsymbol{k}\atop k_{x}^{2}+k_{y}^{2}=k_{h}^{2}}}|\hat{\boldsymbol{u}}_{\boldsymbol{k}}|^{2}.\end{eqnarray}$$

In addition, we monitor different components of domain-integrated energy, namely the total horizontal kinetic energy

(2.4) $$\begin{eqnarray}\frac{1}{2}U_{h}^{2}=\frac{1}{2}\mathop{\sum }_{{\boldsymbol{k}\atop k_{z}=0}}(|\hat{u} _{\boldsymbol{k}}^{(x)}|^{2}+|\hat{u} _{\boldsymbol{ k}}^{(y)}|^{2})\end{eqnarray}$$

(based on the (vertically averaged) 2D2C field only), the large-scale horizontal kinetic energy

(2.5) $$\begin{eqnarray}\frac{1}{2}U_{ls}^{2}=\frac{1}{2}\mathop{\sum }_{{\boldsymbol{k}\atop {k<k_{max}\atop k_{z}=0}}}(|\hat{u} _{\boldsymbol{k}}^{(x)}|^{2}+|\hat{u} _{\boldsymbol{ k}}^{(y)}|^{2}),\end{eqnarray}$$

where $k_{max}=\sqrt{2}2\unicode[STIX]{x03C0}/L$ , as well as the (vertically averaged) large-scale kinetic energy in the $z$ component

(2.6) $$\begin{eqnarray}\frac{1}{2}U_{z}^{2}=\frac{1}{2}\mathop{\sum }_{{\boldsymbol{k}\atop {k<k_{max}\atop k_{z}=0}}}|\hat{u} _{\boldsymbol{k}}^{(z)}|^{2}\end{eqnarray}$$

and the three-dimensional kinetic energy (3-D energy), defined as

(2.7) $$\begin{eqnarray}\frac{1}{2}U_{3D}^{2}=\frac{1}{2}\mathop{\sum }_{{\boldsymbol{k}\atop {k_{z}\neq 0\atop }}}|\hat{\boldsymbol{u}}_{\boldsymbol{k}}|^{2}.\end{eqnarray}$$

3 Results from the direct numerical simulations

In this section, we present the results obtained from our simulations. For a given set of parameters $Re,Q,K$ , two different behaviours are possible. For thick layers $Q\ll 1$ , 3-D turbulence is observed, i.e. there is no inverse cascade and the energy injected by the forcing is transferred to the small scales where it is dissipated. No system-size structures appear in this case. For thin layers $Q\gg 1$ , a split cascade is present with part of the energy cascading inversely to the large scales and part of the energy cascading forward to the small scales. For these layers, at steady state, coherent system-size vortices appear with very large amplitudes.

Figure 2. Typical flow fields in the steady state of 3-D turbulence (a) and 2-D turbulence (b) regimes, visualised using squared vorticity. The boxes below show the corresponding side views. Note the astonishing similarity between this figure and figure 1(a,b) of the experimental study by Xia et al. (Reference Xia, Byrne, Falkovich and Shats2011).

A visualisation of the flow field in these two different states is shown in figure 2 for the 3-D turbulence and condensate states. Typical time series of $U_{h}^{2}$ for a thick layer (forward cascade) and a thin layer (inverse cascade) are shown in figure 3(a). For the thick layer, the total energy fluctuates around a mean value of order $(\unicode[STIX]{x1D716}\ell )^{2/3}$ , while for the thin layer, the energy saturates to a much larger value. The energy spectra for the two runs of figure 3(a) at the steady state are shown in figure 3(b), showing quantitatively that energy is concentrated in the large scales for the two different cases. In more detail, $U_{h}^{2}$ for the thin layer shows two different stages: first, at early times, there is a linear increase with time and second, there is saturation at late times. Therefore, to fully describe the evolution of the system, we need to quantify the rate of the initial energy increase and the energy at which it saturates. The red-dashed line indicates a fit to the initial linear increase. This slope provides a measurement of the rate $\unicode[STIX]{x1D716}_{inv}$ at which energy cascades inversely. For the steady state stage, the black dashed-dot line indicates the mean value at late times. For all runs, we measure the slope of the $U_{h}^{2}$ curve and the steady state mean values of all corresponding energies defined in the previous section. For the runs of high resolution, to accelerate convergence, the large-scale velocity $\boldsymbol{u}_{k=1}$ (from a run at the early stage) was increased artificially and the run continued. Alternatively, an output of a converged run was used as initial condition. However, all cases were run sufficiently long to demonstrate that they have reached a steady state.

Figure 3. (a) Shows the typical time evolution of $U_{h}^{2}$ for $Q=1.25<Q_{3D}$ and $Q=4>Q_{3D}$ . In the former case, $U_{h}^{2}$ remains small. In the latter, there is an initially linear increase whose slope measures the rate of inverse energy transfer. After long time, $U_{h}^{2}$ reaches its steady state value. Two quantities are measured: the initial slope (red-dashed line) and the condensate value (horizontal black dashed-dot line). A similar evolution observed in an experiment is shown in Xia et al. (Reference Xia, Shats and Falkovich2009) figure 6. (b) Shows the corresponding spectra: in the presence of an inverse cascade there is a maximum at the largest scale, while in its absence the maximum is near the forcing scale.

Figure 4. (a) Initial slopes, measured as indicated in figure 3(a), non-dimensionalised by the energy injection rate $\unicode[STIX]{x1D716}$ , as a function of $Q$ $(\propto 1/H)$ for all $Re$ used. The same symbols are used in all plots in this section. Thick layers are at small $Q$ (left) and thin layers at large $Q$ (right). (b) The same data collapsed by a rescaling of the abscissa by $\sqrt{Re}$ and the coordinate by the maximum value obtained for that Reynolds number.

Figure 5. Estimated value of $Q_{3D}$ as a function of $Re$ . The top line shows the smallest value of $Q$ for which an inverse cascade was observed and the bottom line shows the largest value of $Q$ for which no inverse cascade was observed. The rightmost point indicates the results from the hyper-viscous runs.

Figure 4 shows the slopes of the initial total energy increase $\unicode[STIX]{x1D716}_{inv}$ , measured as illustrated in figure 3(a) for all our numerical simulations. The slopes are non-dimensionalised by the energy input rate $\unicode[STIX]{x1D716}$ and plotted versus $Q$ for all different values of $Re$ including the hyper-viscous runs. The slope at this early stage measures the strength of inverse energy transfer. At small $Q$ (deep layers), the slope vanishes for all runs, showing that no inverse cascade is present. Moving to larger $Q$ , for every $Re$ , there is a critical value $Q_{3D}(Re)$ of $Q$ above which the slope becomes non-zero. This is the birth of the inverse cascade. Figure 5 shows estimates of $Q_{3D}$ as a function of $Re$ : the upper curve shows the smallest $Q$ for which an inverse cascade was observed for that given $Re$ while the lower curve shows the largest $Q$ for which no inverse cascade was observed. The critical value $Q_{3D}$ lies between these two curves. The point $Q_{3D}$ shifts to larger $Q$ as $Re$ is increased but eventually for the two largest $Re$ simulated, namely $Re=2031$ and $Re=4062$ , as well as the hyper-viscous run, $Q_{3D}$ saturates at $Q_{3D}\approx 2.5$ . (Previous findings (Celani et al. Reference Celani, Musacchio and Vincenzi2010) estimated this value to $Q_{3D}\approx 2$ , however in that work too limited a range of values of $Q$ was used to be able to precisely pinpoint $Q_{3D}$ . Another possible reason for the different result is the different value of $1/K$ associated with the different forcing wavenumber $k_{f}=16$ used). The saturation of $Q_{3D}\approx 2.5$ indicates that $Q_{3D}$ converges to this value at large $Re$ . For $Q>Q_{3D}$ , the slope begins increasing linearly $\unicode[STIX]{x1D716}_{inv}\propto Q-Q_{3D}$ . (We note that small slopes are hard to distinguish from zero slope since the difference only becomes apparent after a long simulation time.)

If $Q$ is increased further, a point $Q_{2D}$ is reached beyond which the slope becomes independent of $Q$ . Above this second critical point, the flow becomes exactly two-dimensional (Benavides & Alexakis Reference Benavides and Alexakis2017). The value of $Q_{2D}$ increases with $Re$ as $Q_{2D}\propto Re^{1/2}$ . This scaling is verified in this work as well and shown in figure 4(b). The two critical points $Q_{3D}$ and $Q_{2D}$ at this early stage of development of the inverse cascade have been studied in detail in the past (Celani et al. Reference Celani, Musacchio and Vincenzi2010; Benavides & Alexakis Reference Benavides and Alexakis2017).

Figure 6. (a) $U_{ls}^{2}$ as defined in equation (2.5), non-dimensionalised by $(\unicode[STIX]{x1D716}\ell )^{2/3}/K^{2}$ as a function of $Q$ . (b) The same data (excluding hyper-viscous run), with large-scale energy rescaled by $1/Re$ and plotted versus $Q/Re^{3/4}$ showing a satisfactory data collapse. The value $Q_{3D}/Re^{3/4}$ (where $U_{ls}^{2}$ plateaus) coincides for all $Re$ at $Q/Re^{3/4}\approx 0.09{-}0.1$ .

Figure 7. Zoomed-in version of figure 6(a) showing that there is a discontinuity in $U_{ls}^{2}/(\unicode[STIX]{x1D716}\ell )^{2/3}$ at $Q_{3D}$ for all Reynolds numbers up to the second highest simulated.

Here we mostly focus on the second stage of evolution: the steady state and the properties of the condensate. Figure 6 shows the equilibrium value of $U_{ls}^{2}$ , as defined in equation (2.5), non-dimensonalised by the forcing energy scale $(\unicode[STIX]{x1D716}\ell )^{2/3}$ and multiplied by $K^{2}$ . In the left panel, it is plotted versus $Q$ (figure 6 a) and in the right panel it is rescaled by $1/Re$ and plotted versus $\unicode[STIX]{x1D702}/H=QRe^{-3/4}$ (figure 6 b). First consider figure 6(a). At small $Q$ , there is very little energy in the large scales. This corresponds to the values of $Q$ that displayed no inverse cascade at the initial stage. In the absence of an inverse cascade, the large scales only possess a small non-zero energy and are expected to be in a thermal equilibrium state (Kraichnan Reference Kraichnan1973; Dallas, Fauve & Alexakis Reference Dallas, Fauve and Alexakis2015; Cameron, Alexakis & Brachet Reference Cameron, Alexakis and Brachet2017). For $Q>Q_{3D}$ the energy in the large scale takes larger values. In all cases, the energy increases nearly linearly $U_{ls}^{2}\propto (Q-Q_{3D})$ for $Q_{2D}>Q>Q_{3D}$ . With the chosen coordinate, the close coincidence of the experiments with $K=1/8$ and $K=1/16$ at $Re=305$ indicates a scaling of $U_{ls}^{2}\propto L^{2}$ . If we zoom in on $Q_{3D}$ (see figure 7), we observe clear signs of small but discontinuous jumps of $U_{ls}^{2}$ at $Q_{3D}$ that are not visible in the zoomed out figure 6(a). These cases are examined in more detail in the next section.

The increase of the large-scale energy stops at the second critical point $Q_{2D}$ , where $U_{ls}^{2}$ becomes independent of $Q$ . It is noteworthy that the curves for various values of $Re$ all follow the same straight line between their respective $Q_{3D}$ and $Q_{2D}$ with only some deviations at low $Q$ . Furthermore, both $Q_{2D}$ and the plateau value of $U_{ls}^{2}$ depend on $Re$ . In figure 6(b), the same data are plotted, but with rescaled axes. The rescaling collapses the data well, with some deviations at small $Q$ related to the convergence of $Q_{3D}$ . This indicates that at large values of $Q$ , $U_{ls}^{2}$ scales like $U_{ls}^{2}\propto (\unicode[STIX]{x1D716}\ell )^{2/3}Re$ . This is precisely the scaling of the condensate of 2-D turbulence (Boffetta & Ecke Reference Boffetta and Ecke2012). The critical value where the transition to this maximum value of $U_{ls}^{2}$ occurs is $Q_{2D}Re^{-3/4}=\unicode[STIX]{x1D702}/H_{2D}\approx 0.09-0.1$ .

The scaling allowing us to collapse the data in figure 4 (transient stage) is different from that in figures 6, 8 and 9 (condensate state). This implies that $Q_{2D}\propto Re^{1/2}$ estimated during the early stage of the inverse cascade development is different from $Q_{2D}\propto Re^{3/4}$ estimated at steady state where a condensate is fully developed. The reason for this difference is that the transition to exactly 2-D motion occurs when the maximum shear in the flow (which produces 3-D motion by shear instabilities) is balanced by small-scale viscous dissipation. In the presence of the inverse cascade, an $E(k)\propto \unicode[STIX]{x1D716}^{2/3}k^{-5/3}$ spectrum is formed at $k>k_{f}$ , such that the peak of the enstrophy spectrum $k^{2}E(k)$ is at the forcing scale. Thus the balance between 2-D shear and 3-D damping is

(3.1) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D716}\ell )^{1/3}/\ell \sim \unicode[STIX]{x1D708}/H^{2}, & & \displaystyle\end{eqnarray}$$

implying $H_{2D}\sim \ell Re^{-1/2}$ (Benavides & Alexakis Reference Benavides and Alexakis2017). In the condensate, however, most of the energy and enstrophy are located in the largest scales and are such that energy injection $\unicode[STIX]{x1D716}$ is balanced by large-scale dissipation $\propto \unicode[STIX]{x1D708}U_{ls}^{2}/L^{2}$ . The large-scale shear is thus $U_{ls}/L\propto (\unicode[STIX]{x1D716}/\unicode[STIX]{x1D708})^{1/2}$ which is balanced by the damping rate of 3-D perturbations at onset,

(3.2) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D716}/\unicode[STIX]{x1D708})^{1/2}\sim \unicode[STIX]{x1D708}/H^{2}, & & \displaystyle\end{eqnarray}$$

giving the scaling $H_{2D}\propto \unicode[STIX]{x1D716}^{1/4}\unicode[STIX]{x1D708}^{3/4}\propto \ell Re^{-3/4}$ . We will recover the very same steady state scaling in § 6 from a low-order model. These two scalings imply the interesting possibility that a flow which becomes exactly two-dimensional at the early stages of the inverse cascade for $Q\gtrsim Re^{1/2}$ may develop 3-D instabilities at the condensate state if $Q\lesssim Re^{3/4}$ .

Figure 8. (a) $U_{3D}^{2}$ as defined in equation (2.7), non-dimensionalised by $(\unicode[STIX]{x1D716}\ell )^{2/3}$ and plotted versus $Q$ . (b) The same information as (a), but in terms of the square-root of the 3-D kinetic energy rescaled by $(\unicode[STIX]{x1D716}H)^{2/3}$ , plotted versus $Q/Re^{3/4}$ . This rescaling indicates that $U_{3D}^{2}\propto (Q_{2D}-Q)^{2}$ close to the transition.

Figure 9. (a) $U_{z}^{2}$ as defined in equation (2.6) as a function of $Q$ . (b) The various curves collapse when the abscissa is chosen as $\unicode[STIX]{x1D702}/H=Q/Re^{3/4}$ and the coordinate as $U_{3D}^{2}/(\unicode[STIX]{x1D716}H)^{2/3}$ . Raising the coordinate to the $1/3$ power, the curve becomes linear close to onset. This indicates that close to onset, $U_{v}^{2}$ scales as $U_{z}^{2}\approx (Q_{c}-Q)^{3}$ , where $Q_{c}\approx 0.09{-}0.1\approx Q_{2D}$ . We note that the scaling exponent is different from that found for $U_{3D}^{2}$ .

Figure 8 shows $U_{3D}^{2}$ as defined in equation (2.7). In the left panel it is non-dimensionalised by the forcing energy $(\unicode[STIX]{x1D716}\ell )^{2/3}$ and plotted versus $Q$ (figure 8 a), while in the right panel, it is non-dimensionalised by $(\unicode[STIX]{x1D716}H)^{2/3}$ , raised to the power  $1/2$ and plotted versus $\unicode[STIX]{x1D702}/H=Q/Re^{-3/4}$ . Figure 8(a) shows that beyond some non-monotonic behaviour at small $Q$ , $U_{3D}^{2}$ decreases monotonically with $Q$ until it reaches zero at $Q_{2D}$ and remains zero beyond this point. The 3-D energy increases with $Re$ at a given $Q$ . Under the rescaling in figure 8(b), the various curves collapse nicely. In particular, the point where $U_{3D}^{2}$ vanishes is sharp and identical for all $Re$ , namely $\unicode[STIX]{x1D702}/H\approx 0.1$ . Comparing with figure 6(b), one sees that this point and $Q_{2D}$ coincide within the range of uncertainties. This means that beyond $Q_{2D}$ , not only is $U_{ls}^{2}$ independent of $Q$ , but also $U_{3D}^{2}$ vanishes. This confirms that $Q_{2D}$ corresponds to the point where the motion becomes invariant along $z$ .

Finally, figure 9 shows the vertical kinetic energy, non-dimensionalised by $(\unicode[STIX]{x1D716}\ell )^{2/3}$ , once plotted versus $Q$ and once taken to the $1/3$ power and plotted versus $Q/Re^{3/4}$ . The general features of figure 9(a) are similar to figure 8(a): like 3-D energy, vertical kinetic energy decreases with $Q$ until it reaches zero and it increases with $Re$ . The curves collapse in figure 9(b) and the behaviour close to $Q_{2D}$ becomes linear if the coordinate is raised to the $1/3$ power, indicating an approximate scaling $U_{z}^{2}\approx (Q_{c}-Q)^{3}$ with $Q_{c}\approx Q_{2D}$ . This indicates that the point beyond which the vertical kinetic energy vanishes is close to $Q_{2D}$ , implying that beyond $Q_{3D}$ , the motion is not only invariant along $z$ but also restricted to the $x$ – $y$ plane. Hence, for $Q>Q_{2D}$ , the flow has two-dimensionalised exactly.

4 Behaviour close to the transitions: hysteresis and intermittency

In this section, we discuss the behaviour close to the two transition points $Q_{2D}$ and $Q_{3D}$ . Each transition shows a different non-trivial behaviour. Close to $Q_{3D}$ , we observe discontinuous transitions and hysteresis for some range of parameters, while close to $Q_{2D}$ , we find both spatial and temporal intermittency with localised bursts of 3-D energy.

4.1 Close to $Q_{3D}$ : discontinuity and hysteresis

We begin by discussing the behaviour of the flow for $Q$ close to $Q_{3D}$ where a sharp increase of the large-scale energy was observed. This sharp increase could indicate the presence of a discontinuity that could further imply the presence of hysteresis.

To verify the presence of a discontinuity we need perform many different runs varying $Q$ is small steps as well as verifying sensitivity to initial conditions. To do this, a hysteresis experiment has been performed at $Re=406$ , consisting of two series of runs, that we refer to as the ‘upper branch’ and the ‘lower branch’, see figure 10. On the upper branch, we start with random initial conditions and $Q\approx 2.25$ for which the system reaches a condensate equilibrium with an associated non-zero value of large-scale energy. Once the run has equilibrated, we use that equilibrium state to initialise a run at $Q\rightarrow Q-\unicode[STIX]{x0394}Q$ with $\unicode[STIX]{x0394}Q=0.1$ . By decreasing $Q$ , the physical height of the box is increased. To be able to use the equilibrium state reached at one $Q$ as initial condition for a neighbouring $Q$ , the $z$ -dependence of the velocity field is scaled and the velocity field is projected onto its diverge-free part, formally $v(x,y,z)\rightarrow \mathbb{P}v(x,y,\unicode[STIX]{x1D706}z)$ , where $\mathbb{P}=\mathbb{I}-\unicode[STIX]{x1D6FB}^{-2}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\,)$ . Having changed $Q$ and applied this procedure, we let the system equilibrate to a new condensate state. This is repeated five more times (step size reduced to $\unicode[STIX]{x0394}Q=0.05$ and then $\unicode[STIX]{x0394}Q=0.025$ ) down to $Q\approx 1.9$ . When $Q$ is now lowered $0.025$ further, the condensate decays into 3-D turbulence and the large-scale energy saturates to close to zero. Reducing $Q$ even more, $U_{ls}^{2}$ remains small, indicating a 3-D turbulent state. The lower branch was calculated similarly, with the only difference that the experiment started at low $Q$ and $Q$ was increased in steps of $0.05$ and then of $0.025$ . For small $Q$ , the two branches coincide, while the lower branch remains at low $U_{ls}^{2}$ (3-D turbulence) up to $Q\approx 2.025$ . For $Q$ larger than $Q=2.025$ , the lower branch merges with the upper branch, closing the hysteresis loop and a condensate is spontaneously formed from 3-D turbulence. In other words, for $Re=406$ in the range $1.9\leqslant Q\leqslant 2.025$ , there are multiple steady states and to which state the system will saturate depends on the initial conditions. The flow field for two such states starting from different initial conditions for $Q\approx 1.97$ is visualized in figure 11.

Figure 10. Hysteresis curve of $U_{ls}^{2}$ non-dimensionalised by the forcing energy. Two experiments are shown, the ‘lower branch’ starting from small $Q$ (deep layer) and increasing $Q$ and the ‘upper branch’ starting from large $Q$ (thin layer) and decreasing  $Q$ .

Figure 11. Visualisation of the typical flow field after long simulation time at $Q\approx 1.97$ in the hysteresis experiment on the upper (b) and lower (a) branches. The lower branch flow field shows small-scale structures and no large-scale organisation, reminiscent of 3-D turbulence. By contrast, the upper branch flow field is characterised by two large-scale vortices in addition to smaller-scale structures in between them.

The following remarks are in order: although for each $Q$ we ran the simulations until saturation was achieved, since we are dealing with a noisy system, rare transitions can exist between the two branches of the hysteresis loop. To test this, we picked the point $Q\approx 1.97$ on both lower and upper branches and ran them for a long time (thousands of eddy turnover times $\unicode[STIX]{x1D70F}=(L^{3}/\unicode[STIX]{x1D716})^{1/2}$ ). In neither case did we see a transition between the two branches, indicating that such transitions are rare (if not absent) in the middle of the hysteresis loop. Near the edges of the hysteresis loop at $Q\approx 2.05$ and $Q\approx 1.9$ , the dependence on simulation time is likely to be stronger, but this has not been investigated.

Furthermore, we note that the bifurcation diagram of figure 10 corresponds to a relatively low Reynolds number $Re=406$ . Whether this subcritical behaviour persists at larger $Re$ and/or larger box sizes (smaller $K$ ) is still an open question. Figure 7 suggests that a discontinuity continues to exist at $Q=Q_{3D}$ up to high Reynolds numbers ( $Re=2031$ shown there). In addition, we found more points at higher $Re$ that showed a dependence on initial conditions but without having enough values of $Q$ to create a hysteresis diagram. These findings suggest that subcritical behaviour and hysteresis might survive even at high $Re$ . However, due to the high computational cost at higher resolution and the long duration of the runs required to verify that the system stays in a particular state, we could not investigate this possibility in detail. Further simulations at larger $Re$ and possibly smaller $K$ (larger boxes) are required to resolve this issue. Similar hysteretic behaviour has recently been reported in rotating turbulence, see Yokoyama & Takaoka (Reference Yokoyama and Takaoka2017). More generally, multistability is observed in many turbulent flows, see Weeks et al. (Reference Weeks, Tian, Urbach, Ide, Swinney and Ghil1997), Ravelet et al. (Reference Ravelet, Marié, Chiffaudel and Daviaud2004) as examples.

Figure 12. Plots showing temporal intermittency at $Re=203$ . (a) Shows a typical time series (on linear–log axes) of 3-D energy close to $Q_{2D}$ . Specifically, $Q=5$ , while $Q_{2D}\approx 5.13$ at this value of $Re$ . (b) Shows probability density functions (PDFs) corresponding to this time series as well as for different values of $Q$ (PDFs shifted by a constant factor for better visibility). The different symbols mark different values of $Q$ , while the dotted lines correspond to power laws with exponents $-1$ (bottom), $-0.8$ (middle) and $-0.3$ (top) respectively.

4.2 Close to $Q_{2D}$ : intermittent bursts

Next, we discuss the behaviour of the flow close to the second critical point $Q_{2D}$ . A typical time series of 3-D energy for $Q\lesssim Q_{2D}$ is shown in figure 12(a). One observes bursty behaviour and variations over many orders of magnitude, indicating on–off intermittency (Fujisaka & Yamada Reference Fujisaka and Yamada1985; Platt, Spiegel & Tresser Reference Platt, Spiegel and Tresser1993). On–off intermittency refers to the situation where a marginally stable attractor loses or gains stability due to noise fluctuations. When instability is present, a temporary burst is produced before the system returns to the attractor. On–off intermittency predicts that the unstable mode $X$ follows a power-law distribution $P(X)\propto X^{\unicode[STIX]{x1D6FF}-1}$ for $X\ll 1$ where $\unicode[STIX]{x1D6FF}$ measures the deviation from onset (here $\unicode[STIX]{x1D6FF}\propto (Q_{2D}-Q)/Q_{2D}$ ) and all moments scale linearly with the deviation $\langle X^{n}\rangle \propto \unicode[STIX]{x1D6FF}$ .

In our system, the 2-D flow forms the marginal stable attractor that loses stability to 3-D perturbations depending on the exact realisation of the 2-D turbulent flow. To formulate this, we decompose the velocity field into its 2-D and 3-D parts, $\boldsymbol{u}=\boldsymbol{u}_{2D}+\boldsymbol{u}_{3D}$ , where the 2-D part is defined as the Fourier sum of $\boldsymbol{u}$ restricted to modes with $k_{z}=0$ . Filtering the 3-D component of equation (2.1b ), dotting with $\boldsymbol{u}_{3D}$ and integrating over the domain gives

(4.1) $$\begin{eqnarray}{\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{t}U_{3D}^{2}=-\langle \{\boldsymbol{u}_{3D}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{2D}\}\boldsymbol{\cdot }\boldsymbol{u}_{3D}\rangle -\unicode[STIX]{x1D708}\langle |\unicode[STIX]{x1D735}\boldsymbol{u}_{3D}|^{2}\rangle ,\end{eqnarray}$$

where $\langle \cdot \rangle$ denotes integration over the domain. The chaotic 2-D motions then act as multiplicative noise while the viscous terms provide a mean decay rate. An important physical mechanism for creating 3-D disturbances is 3-D elliptic instability of the 2-D counter-rotating vortex pair forming the condensate, as described in Le Déz & Laporte (Reference Le Déz and Laporte2002). This may explain the presence of the critical value $Q_{2D}$ itself: instability requires small vertical wavenumbers, but the minimum wavenumber increases with decreasing $H$ and $Q_{2D}$ corresponds to the point where 3-D perturbations begin to decay.

Figure 12 shows that temporal intermittency is present in the thin-layer system. Panel (a) shows a typical time series of 3-D energy at $Q\lesssim Q_{2D}$ which fluctuates over six orders of magnitude. In particular, as mentioned before, there are burst-like excursions in 3-D energy. In figure 12(b), PDFs constructed from this time series and similar ones for different values of $Q$ are shown along with dotted lines indicating power laws with exponents $-1$ , $-0.8$ and $-0.3$ . The PDFs are very close to a power law for a significant range of $U_{3D}^{2}$ and the exponent converges to minus one as the transition is approached, in agreement with on–off intermittency predictions. However, the scaling of 3-D energy with deviation from onset shown in figure 8(b) does not follow the linear prediction of on–off intermittency, but rather $\langle U_{3D}^{2}\rangle \propto (Q_{2D}-Q)^{2}$ . For $U_{z}^{2}$ , figure 9 seems to suggest yet a different scaling, namely $\langle U_{z}^{2}\rangle \propto (Q_{2D}-Q)^{3}$ . A similar behaviour was also found in Benavides & Alexakis (Reference Benavides and Alexakis2017) and was attributed to the spatio-temporal character of the intermittency that not only leads to 3-D motions appearing more rarely in time as criticality is approached but also to them occupying a smaller fraction of the available volume. This appears also to be the case in our results, as demonstrated in figures 13 and 14, where $\boldsymbol{u}_{3D}^{2}$ and the vertical velocity $u_{z}$ are plotted for three different values of $Q$ . As $Q$ approaches the critical value $Q_{2D}$ , the structures become smaller for $\boldsymbol{u}_{3D}^{2}$ and $u_{z}$ with the difference that $\boldsymbol{u}_{3D}^{2}$ shows spot-like structures in figure 13(c) which are absent for $u_{z}$ . This difference may be related to the two different scalings observed for $U_{3D}^{2}$ and $U_{z}^{2}$ with $Q_{c}-Q$ small: if the volume fraction of vertical motion depends on $Q_{c}-Q$ to a different power than that of vertical variations, two different behaviours of $U_{z}^{2}$ and $U_{3D}^{2}$ would follow. A more detailed quantitative investigation of the scaling of volume fraction will be needed to clarify this.

In summary, we have found non-trivial behaviour close to both transitions: we have observed hysteresis near $Q_{3D}$ and spatio-temporal intermittency close to $Q_{2D}$ where the temporal behaviour seems to be described by on–off intermittency. Taking into account these effects will be crucial for understanding the exact nature of the observed transitions.

Figure 13. Snapshots of $\boldsymbol{u}_{3D}^{2}$ for $Re=203$ and $Q=2.5$ (a), $Q=5.0$ (b) and $Q=5.125$ (c) (corresponding to figure 12). The colour bar is chosen in each plot such that the maximum value of $\boldsymbol{u}_{3D}^{2}$ is shown in black. As $Q$ increases towards $Q_{2D}\approx 5.13$ , $U_{3D}^{2}$ becomes more and more localised. In (c), $\boldsymbol{u}_{3D}^{2}$ is concentrated in small columnar structures (upper and lower right-hand corner) absent in figure 14(c).

Figure 14. Snapshots of $v_{z}$ for $Re=203$ and $Q=2.0$ (a), $Q=5.0$ (b) and $Q=5.125$ (c). The colour scale on the right is chosen in each plot such that the minimum (negative) is shown in black and the maximum (positive) in red. (same as in figure 12). As in figure 13, $u_{z}$ becomes localised in smaller and smaller areas as $Q$ increases, i.e. there is increasing spatial intermittency as $Q_{2D}\approx 5.13$ is approached.

5 Spectra and fluxes

In this section, we discuss the spectral space properties of the three different regimes described in the previous section. For this purpose, it is necessary to define a few additional quantities. In addition to the total 1-D energy spectrum defined in (2.3), is of interest to consider the two-dimensional energy spectrum in the $(k_{h},k_{z})$ plane.

(5.1) $$\begin{eqnarray}E(k_{h},k_{z})=\frac{1}{2}\mathop{\sum }_{{\boldsymbol{k}^{\prime }\atop {k_{x}^{\prime 2}+k_{y}^{\prime 2}=k_{h}^{2}\atop k_{z}^{\prime }=k_{z}}}}|\hat{\boldsymbol{u}}_{\boldsymbol{k}^{\prime }}|^{2}.\end{eqnarray}$$

Moreover, the total 1-D energy spectrum may advantageously be split up into three components: the energy spectrum of the (vertically averaged) 2D2C field

(5.2) $$\begin{eqnarray}E_{h}(k_{h})=\frac{1}{2}\mathop{\sum }_{{\boldsymbol{k}\atop {k_{x}^{2}+k_{y}^{2}=k_{h}^{2}\atop k_{z}=0}}}(|\hat{u} _{\boldsymbol{k}}^{(x)}|^{2}+|\hat{u} _{\boldsymbol{ k}}^{(y)}|^{2}),\end{eqnarray}$$

the energy spectrum of the (vertically averaged) vertical velocity

(5.3) $$\begin{eqnarray}E_{z}(k_{h})=\frac{1}{2}\mathop{\sum }_{{\boldsymbol{k}\atop {k_{x}^{2}+k_{y}^{2}=k_{h}^{2}\atop k_{z}=0}}}|\hat{u} _{\boldsymbol{k}}^{(z)}|^{2},\end{eqnarray}$$

and the energy spectrum of the 3-D flow defined as

(5.4) $$\begin{eqnarray}E_{3D}(k_{h})=\frac{1}{2}\mathop{\sum }_{{\boldsymbol{k}\atop {k_{x}^{2}+k_{y}^{2}=k_{h}^{2}\atop k_{z}\neq 0}}}|\hat{\boldsymbol{u}}_{\boldsymbol{k}}|^{2},\end{eqnarray}$$

satisfying $E_{tot}(k_{h})=E_{h}(k_{h})+E_{z}(k_{h})+E_{3D}(k_{h})$ . Furthermore, we introduce three different quantities related to spectral energy flux. First, the total energy flux as a function of horizontal wavenumber

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}(k_{h})=\langle \boldsymbol{u}_{k_{h}}^{{<}}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}\rangle ,\end{eqnarray}$$

where the low-pass filtered velocity field is

(5.6) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{k_{h}}^{{<}}=\mathop{\sum }_{{\boldsymbol{k}\atop k_{x}^{2}+k_{y}^{2}<k_{h}^{2}}}\hat{\boldsymbol{u}}_{\boldsymbol{k}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}. & & \displaystyle\end{eqnarray}$$

With this definition, $\unicode[STIX]{x1D6F1}(k_{h})$ expresses the flux of energy through the cylinder $k_{x}^{2}+k_{y}^{2}=k_{h}^{2}$ due to the nonlinear interactions. The 2-D energy flux as a function of $k_{h}$ is defined as

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{2D}(k_{h})=\langle \overline{\boldsymbol{u}}_{k_{h}}^{{<}}\boldsymbol{\cdot }(\overline{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\overline{\boldsymbol{u}}\rangle ,\end{eqnarray}$$

where the over-bar stands for vertical average and expresses the flux through the same cylinder due to only 2D2C interactions. Finally, we define the 3-D energy flux (due to all interactions other than those in (5.7)) as a function of horizontal wavenumber by

(5.8) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{3D}(k_{h})=\unicode[STIX]{x1D6F1}(k_{h})-\unicode[STIX]{x1D6F1}_{2D}(k_{h}).\end{eqnarray}$$

It expresses the flux due to all interactions other than the ones in (5.7).

Figure 15. Logarithmic surface plots of $E(k_{h},k_{z})$ at steady state in the three regimes (a) $Q<Q_{3D}$ , (b) $Q_{3D}<Q<Q_{2D}$ and (c) $Q_{2D}<Q$ .

Figure 16. Three different energy spectra, $E_{h}(k_{h})$ , $E_{z}(k_{h})$ , $E_{3D}(k_{h})$ at $Re=609$ for 3-D turbulence (a), 2-D turbulence (c) and an intermediate case $Q\in (Q_{3D},Q_{2D})$ (b) flux-loop condensate (cf. main text). For 3-D turbulence ( $Q=1.25<Q_{3D}$ ), the 2-D energy spectrum peaks at the forcing scale and is an order of magnitude bigger than the other components. In the flux-loop condensate ( $Q=4$ ), 2-D energy is maximum at $k=1$ and 3-D energy and vertical energy are non-zero. In 2-D turbulence ( $Q=16>Q_{2D}$ ), 2-D energy is maximum at $k=1$ , but 3-D energy and vertical energy vanish.

Figure 17. Three different components of spectral energy flux, $\unicode[STIX]{x1D6F1}(k_{h})$ , $\unicode[STIX]{x1D6F1}_{2D}(k_{h})$ and $\unicode[STIX]{x1D6F1}_{3D}(k_{h})$ , are shown for the same three cases and in the same order as in figure 16.

Figure 15 shows the steady state 2-D energy spectrum in the three different regimes: (a) $Q<Q_{3D}$ , (b) $Q_{3D}<Q<Q_{2D}$ and (c) $Q_{2D}<Q$ . In the 3-D turbulent case (a), the global maximum is at the forcing scale and $k_{z}=0$ , while large $k_{z}$ modes have a relatively larger fraction of total energy than in cases (b,c). In cases (b,c), a condensate is present with a maximum at the largest wavenumber $k_{h}=1,k_{z}=0$ . In case (b), there is still energy in the $k_{z}\neq 0$ modes, while in case (c), the energy is entirely concentrated in the $k_{z}=0$ mode. Figure 16 shows the energy spectra $E_{h}(k_{h}),E_{z}(k_{h})$ and $E_{3D}(k_{h})$ for the same three cases (a–c). In case (a) (3-D turbulence), all three spectra are of the same order, with a small excess of $E_{h}(k_{h})$ in the large scales and an excess of $E_{3D}(k_{h})$ in the small scales. The small scale separation between the forcing and the dissipation scale does not allow us to observe a $k^{-5/3}$ power-law regime. In case (b), $E_{h}(k_{h})$ clearly dominates in the large scales, forming a steep spectrum (close to $E_{h}(k_{h})\propto k_{h}^{-4}$ ). However, at wavenumbers larger than the forcing wavenumber $k_{f}=8$ , $E_{z}(k_{h})$ and $E_{3D}(k_{h})$ become of the same order as $E_{h}(k_{h})$ . In case (c) (2-D turbulence), where $Q>Q_{2D}$ , the spectra $E_{z}(k_{h})$ and $E_{3D}(k_{h})$ have reduced to values close to the round-off error and are not plotted. The 2-D spectrum $E_{h}(k_{h})$ displays again a steep power-law behaviour close to $E_{h}(k_{h})\propto k_{h}^{-4}$ .

Figure 17 shows the energy fluxes as defined in equations (5.5)–(5.8) for the same three cases examined in figure 16. In panel (a), where the case $Q<Q_{3D}$ is examined, there is almost no inverse flux of energy and $\unicode[STIX]{x1D6F1}(k_{h}<k_{f})$ is practically zero. The small inverse flux that is observed for $\unicode[STIX]{x1D6F1}_{2D}(k_{h})$ at $k<k_{f}$ does not reach the largest scale of the system and is nearly completely balanced by $\unicode[STIX]{x1D6F1}_{3D}(k_{h})$ , which is forward. At wavenumbers larger than $k_{f}$ , the total flux is positive and is completely dominated by $\unicode[STIX]{x1D6F1}_{3D}$ . This is to be contrasted with the rightmost panel (c) with $Q>Q_{2D}$ , where at small wavenumbers, the flux is negative and is dominated by the 2-D flow, while at large wavenumbers there is a very small forward flux. For the intermediate case $Q_{3D}<Q<Q_{2D}$ in panel (b), there is an inverse energy flux. This flux can be decomposed into a negative 2-D part  $\unicode[STIX]{x1D6F1}_{2D}(k_{h})$ and a positive 3-D part $\unicode[STIX]{x1D6F1}_{3D}(k_{h})$ . In other words, the 2-D components of the flow bring energy to the largest scales of the system, which is then brought back to the small scales by the 3-D components of the flow associated with a forward energy flux, thus forming a loop for the energy transfer. For this reason, we refer to this case as flux-loop condensate.

Due to finite viscosity, part of the energy that arrives at the largest scale (shown in figure 17 b) is dissipated. Therefore, the two fluxes are not completely in balance. As $Re$ is increased, however, the fraction of the energy that is dissipated in the large scales is decreased and the two opposite fluxes come closer to balancing each other. This is shown in figure 18, where the energy fluxes for the highest $Re$ simulation and for the simulation with hyper viscosity are plotted. The two opposite directed fluxes are closer in amplitude. At $Re\rightarrow \infty$ it is thus expected that the inverse and forward fluxes at large scales will be in perfect balance and all the energy is dissipated in the small scales. It is worth noting, however, that the inverse cascade (negative flux) due to the 2-D components has much stronger fluctuations than the forward cascading flux that has led to the non-monotonic behaviour of the flux observed in figure 18 at small $k$ due to insufficient time averaging.

Figure 18. Flux-loop condensate steady state fluxes for $Re=4062$ in (a) and the hyper-viscous run in (b).

6 A three-mode model

In this section, we formulate and analyse a simple three-scale ordinary differential equation (ODE) model which reproduces certain features of the DNS results described in § 3.

Figure 19. Sketch of three-mode model. Solid curve: energy spectrum $E(k)$ of the condensate state. The energy injected at $k_{f}$ at a rate $\unicode[STIX]{x1D716}$ is distributed between large and small scales. Moreover, energy is transferred from large to small scales. Finally viscous dissipation occurs at all scales since $Re<\infty$ (short arrows on abscissa) and energy is transferred to the dissipation range (arrow at $k_{\unicode[STIX]{x1D708}}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D702}$ ). The spectrum $E(k)$ shares certain features with figure 1(e) of Xia et al. (Reference Xia, Byrne, Falkovich and Shats2011) and figure 3 of Celani et al. (Reference Celani, Musacchio and Vincenzi2010).

As illustrated in figure 19, our model comprises a 2-D mode $U_{2D}$ at the scale $L$ of the domain, a mode $U_{f}$ at the forcing scale $\ell$ and a 3-D mode $U_{3D}$ at the scale of the layer height $H$ , whose interactions are spectrally non-local, thus taking into account a major result from § 5. The model describes the system at steady state where these scales are well separated, but is not expected to capture the transient phase where all intermediate scales between $L$ and $\ell$ participate due to the inverse cascade. As before, let $Q=\ell /H$ , $K=\ell /L$ and $Re=(\unicode[STIX]{x1D716}\ell ^{4})^{1/3}/\unicode[STIX]{x1D708}$ . Interactions between modes are modelled using eddy viscosity, which amounts to modifying the molecular viscosity $\unicode[STIX]{x1D708}$ by terms involving the small-scale velocities, modelling the effect of small-scale motions on large-scale motions as diffusive. The conceptual foundations of eddy viscosity were laid by de Saint Venant in his effective viscosity, (de Saint-Venant Reference de Saint-Venant1843) (see Darrigol (Reference Darrigol2017) for a historical review). Eddy viscosity was quantified for the first time by Boussinesq (Reference Boussinesq1877) and later widely popularised through the works of Taylor (Taylor Reference Taylor1915, Reference Taylor1922), see also Kraichnan (Reference Kraichnan1976). It has been estimated in various limits both in 2-D and 3-D flows (Meshalkin Reference Meshalkin1962; Sivashinsky & Yakhot Reference Sivashinsky and Yakhot1985; Bayly & Yakhot Reference Bayly and Yakhot1986; Yakhot & Sivashinsky Reference Yakhot and Sivashinsky1987; Hefer & Yakhot Reference Hefer and Yakhot1989; Dubrulle & Frisch Reference Dubrulle and Frisch1991; Gama, Vergassola & Frisch Reference Gama, Vergassola and Frisch1994; Cameron, Alexakis & Brachet Reference Cameron, Alexakis and Brachet2016; Alexakis Reference Alexakis2018).

There are two notable cases where the dependence of eddy viscosity $\unicode[STIX]{x1D708}_{E}$ on the flow amplitude $U_{s}$ and length scale $l_{s}$ is known. For $Re\rightarrow \infty$ , one expects that $\unicode[STIX]{x1D708}_{E}$ becomes independent of $\unicode[STIX]{x1D708}$ and the only dimensionally consistent possibility for $\unicode[STIX]{x1D708}_{E}$ is given by

(6.1) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{E}=c_{1}U_{s}l_{s},\end{eqnarray}$$

where $c_{1}$ is a non-dimensional number. In the low- $Re$ limit, on the other hand, an exact asymptotic expansion can be carried out (see Dubrulle & Frisch Reference Dubrulle and Frisch1991) which reveals that

(6.2) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{E}=c_{2}\frac{U_{s}^{2}l_{s}^{2}}{\unicode[STIX]{x1D708}}+O(U_{s}^{4}l_{s}^{4}/\unicode[STIX]{x1D708}^{2}),\end{eqnarray}$$

where the non-dimensional number $c_{2}$ can be evaluated by the expansion. It may seem counter-intuitive that the low- $Re$ limit could have any relevance for the turbulent problem, but since we have established in the DNS that the presence of $Q_{2D}$ is a finite- $Re$ phenomenon ( $Q_{2D}\propto Re^{3/4}$ ), we clearly need to include a finite $Re$ ingredient to describe it and the exact result (6.2) is selected for this purpose. The sign of the prefactors $c_{1},c_{2}$ depends on the exact form of the small-scale flow and in particular its dimensionality. While 2-D flows tend to have negative eddy viscosities and transfer energy upscale, 3-D flows are expected to have positive eddy viscosities and transfer energy downscale. For our model, we are going to consider that interactions among the three different scales $L>\ell >H$ are such that the flow at the smaller scale acts as an eddy viscosity on the flow at the larger scale. These interactions are illustrated in figure 19. In particular, the energy injected at the forcing scale $k_{f}$ at a rate $\unicode[STIX]{x1D716}$ is transferred both to the large scale $L$ (by a negative eddy viscosity $-\unicode[STIX]{x1D707}$ ) and to the small scales (by a positive eddy viscosity $\unicode[STIX]{x1D70E}$ ). The large scales lose energy directly to the small scales (via a positive eddy viscosity term $\unicode[STIX]{x1D702}$ ), while the small scales dissipate energy by transfer to the dissipation range, modelled by a nonlinear energy sink. In addition, viscosity is finite, such that all scales dissipate locally. The set of equations below formalises these ideas:

(6.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}U_{2D}^{2}=-(\unicode[STIX]{x1D708}-\unicode[STIX]{x1D707}+\unicode[STIX]{x1D702})\frac{U_{2D}^{2}}{L^{2}}, & \displaystyle\end{eqnarray}$$

(6.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}U_{f}^{2}=\unicode[STIX]{x1D716}-(\unicode[STIX]{x1D708}+\unicode[STIX]{x1D70E})\frac{U_{f}^{2}}{\ell ^{2}}-\unicode[STIX]{x1D707}\frac{U_{2D}^{2}}{L^{2}}, & \displaystyle\end{eqnarray}$$

(6.3c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}U_{3D}^{2}=\unicode[STIX]{x1D702}\frac{U_{2D}^{2}}{L^{2}}+\unicode[STIX]{x1D70E}\frac{U_{f}^{2}}{\ell ^{2}}-\frac{U_{3D}^{3}}{H}-\unicode[STIX]{x1D708}\frac{U_{3D}^{2}}{H^{2}}. & \displaystyle\end{eqnarray}$$

(6.4) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}(U_{2D}^{2}+U_{f}^{2}+U_{3D}^{2})=\unicode[STIX]{x1D716}-\unicode[STIX]{x1D708}\left(\frac{U_{2D}^{2}}{L^{2}}+\frac{U_{f}^{2}}{\ell ^{2}}+\frac{U_{3D}^{2}}{H^{2}}\right)-\frac{U_{3D}^{3}}{H}, & & \displaystyle\end{eqnarray}$$

showing that the total kinetic energy only changes due to molecular viscosity $\unicode[STIX]{x1D708}$ , energy injection $\unicode[STIX]{x1D716}$ and the sink term representing the 3-D energy cascade to the dissipation range, $U_{3D}^{3}/H$ . Depending on $Re$ , either of the two expressions for eddy viscosity (6.1), (6.2) may be expected to yield an adequate description of the multi-scale interactions in the problem. A model that interpolates smoothly between the large and small $\unicode[STIX]{x1D708}$ limits, thus taking into account the finite- $Re$ information necessary for describing $Q_{2D}$ , is given by

(6.5a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D707}=\unicode[STIX]{x1D6FC}\frac{U_{f}^{2}\ell ^{2}}{\unicode[STIX]{x1D708}+U_{f}\ell },\quad \unicode[STIX]{x1D702}=\unicode[STIX]{x1D6FD}\frac{U_{3D}^{2}H^{2}}{\unicode[STIX]{x1D708}+U_{3D}H},\quad \unicode[STIX]{x1D70E}=\unicode[STIX]{x1D6FE}\frac{U_{3D}^{2}H^{2}}{\unicode[STIX]{x1D708}+U_{3D}H},\end{eqnarray}$$

with $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE}>0$ non-dimensional coupling constants. In the limits $\unicode[STIX]{x1D708}\rightarrow 0$ and $\unicode[STIX]{x1D708}\rightarrow \infty$ , the above expressions converge to the formulae for eddy viscosities described before. The nonlinear dynamical system thus defined possesses a varying number of fixed points depending on parameters. To classify them, first note that $\unicode[STIX]{x1D716}\neq 0\Rightarrow U_{f}\neq 0$ at any fixed point by (6.3b ) and the definition of $\unicode[STIX]{x1D707}$ in (6.5). Hence there are four possibilities:

1. (i) laminar state: $U_{2D}=U_{3D}=0$ (all energy in forcing scale),

2. (ii) three-dimensional turbulence state: $U_{2D}=0$ and $U_{3D}\neq 0$ ,

3. (iii) two-dimensional condensate state: $U_{2D}\neq 0$ and $U_{3D}=0$ and

4. (iv) flux-loop condensate state: $U_{2D}\neq 0$ and $U_{3D}\neq 0$ .

As shown in § A.1, in the zero viscosity limit, there is neither a laminar state nor a 2-D condensate fixed point in the model. This emphasises the importance of including finite- $Re$ information into the model for describing both $Q_{3D}$ and $Q_{3D}$ in a single model. The laminar state appears for values of $Re\equiv (\unicode[STIX]{x1D716}\ell ^{4})^{1/3}/\unicode[STIX]{x1D708}$ below a critical value $Re_{c}$ for which there is no transfer, neither to large nor to small scales. Above this critical value, one of the three other states is stable, depending on the value of $Q=\ell /H$ . For small values of $Q$ (large $H$ ), the system is in the 3-D turbulence state, where energy is only exchanged between the forcing scale $\ell$ and the small scale $H$ . Above the critical value $Q_{3D}$ , the system transitions to the flux-loop condensate state where part of the injected energy is transferred to the large scales and then back to the small scales, thus forming a loop. Finally, at sufficiently large $Q$ above a second critical point $Q_{3D}$ , the system transitions to the 2-D condensate where it follows 2-D dynamics and there is only a transfer of the injected energy to the large scales.

Figure 20. Steady state $U_{3D}^{2}$ from the full model for four different $Re$ as well as for the $Re=\infty$ limit. The $Re=500$ and the $Re=\infty$ cases are almost indistinguishable. The parameters used are $L=\ell =1$ , $\unicode[STIX]{x1D6FC}=1$ , $\unicode[STIX]{x1D6FD}=5$ , $\unicode[STIX]{x1D6FE}=0.5$ .

Figure 21. Model steady state $U_{2D}^{2}$ with and without rescaling. The overall structure of the plots is strongly reminiscent of corresponding DNS results in figure 6(a,b). In (a), $U_{2D}$ vanishes at small $Q$ , increases monotonically between $Q_{3D}$ and $Q_{2D}$ and remains constant for $Q>Q_{2D}$ . (b) Shows the same data as (a) with the abscissa rescaled by $Re^{3/4}$ and the coordinate rescaled by $Re$ as in figure 6(b). The collapse improves with increasing $Re$ . Parameters: $\unicode[STIX]{x1D6FC}=0.51$ , $\unicode[STIX]{x1D6FD}=8$ , $\unicode[STIX]{x1D6FE}=0.1$ , $\ell =L/15=1$ .

From this simple model, three major predictions may be derived:

1. (i) Firstly, the critical point $Q_{3D}$ is predicted to converge to a $Re$ -independent value at large $Re$ as is shown in figure 20. In fact, in the infinite $Re$ limit of the model, there remains only one bifurcation, namely that at $Q_{3D}$ between two-dimensional turbulence and the split cascade state.

2. (ii) Secondly, the critical point $Q_{2D}$ is predicted to obey

   (6.6) $$\begin{eqnarray}Q_{2D}\propto Re^{3/4}.\end{eqnarray}$$

3. (iii) Thirdly, for $Q>Q_{2D}$ , i.e. in the 2-D turbulent state, the steady state energy is predicted to be

   (6.7) $$\begin{eqnarray}U_{2D}^{2}=(\unicode[STIX]{x1D716}\ell )^{2/3}\left(\frac{L}{\ell }\right)^{2}Re.\end{eqnarray}$$

The detailed derivations of these results are given in § A.1. All these three main features are in agreement with the DNS and therefore the diagram that displays the different phases of the model, shown in figure 21(a), resembles the corresponding figure 6(a) from the DNS. Indeed, the same rescaling collapses the curves in both cases, see figures 6(b) and 21(b). We also note that for $0<Q_{2D}-Q\ll 1$ , it is predicted that $U_{3D}^{2}\sim (Q_{2D}-Q)^{2}$ (see § A.2), again in agreement with the DNS.

We understand the present ODE model as a mean-field description which captures the global system behaviour and averaged quantities, but does not take fluctuations into account. Due to the importance of fluctuations near criticality, the ODE model does not reproduce the detailed behaviour there. However, when a fluctuating energy input is taken into account by replacing $\unicode[STIX]{x1D716}\rightarrow \unicode[STIX]{x1D716}+\unicode[STIX]{x1D70E}\unicode[STIX]{x1D701}$ in (6.3b ) ( $\unicode[STIX]{x1D701}$ being white Gaussian noise), on–off intermittency is found close to $Q_{2D}$ where the PDF of $U_{3D}^{2}$ follows a power law with an exponent tending to $-1$ as $Q\rightarrow Q_{2D}$ from below (see § A.3), just as in the DNS. This is a consequence of the structure of the model equations.

To conclude this section, we reiterate that the model presented above successfully captures the location of the critical points (up to a scaling factor) as well as the amplitude of the condensate $U_{2D}^{2}$ , while not producing a hysteresis. The intermittency close to $Q_{2D}$ found in DNS is reproduced by the model when additive noise is included.