2 Model description and methodology

We consider an incompressible flow in a thin layer of thickness $H=\unicode[STIX]{x03C0}/q$ in the vertical direction ( $z$ ) and of size $2\unicode[STIX]{x03C0}L$ in the remaining two directions $(x,y)$ with free slip boundary conditions $\unicode[STIX]{x2202}_{z}u_{x}=\unicode[STIX]{x2202}_{z}u_{y}=u_{z}=0$ at $z=\pm H/2$ . The flow velocity $\boldsymbol{u}$ is governed by the Navier–Stokes equations:

where $\unicode[STIX]{x1D708}$ is the kinematic viscosity, $P$ is the pressure and $\boldsymbol{F}$ is a two-dimensional forcing (varying only along $x$ and $y$ ) that acts at a particular horizontal length scale  $\ell _{f}$ .

As discussed in § 1, in this work we do not solve for the full system requiring 3-D numerical simulations, but we rather focus on a model that is obtained by a severe Galerkin truncation in the vertical direction such that only two modes are kept. The first mode $\boldsymbol{u}_{2D}$ has two components that are obtained by vertically averaging $u_{x}$ and $u_{y}$ . This corresponds to a purely 2-D flow which thus satisfies the incompressibility condition $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{2D}=0$ . This allows us to write $\boldsymbol{u}_{2D}$ in terms of a streamfunction $\unicode[STIX]{x1D713}(t,x,y)$ as indicated in (2.3). The second mode $\boldsymbol{u}_{q}$ has all three components and a prescribed vertical dependence as given below:

The vector field $\boldsymbol{u}_{q}$ satisfies free slip boundary conditions at $z=\pm H/2$ and the incompressibility condition $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{q}=0$ that, in terms of $\boldsymbol{v}_{q}(t,x,y)=(v_{x},v_{y},v_{z})$ , is written as $\unicode[STIX]{x2202}_{x}v_{x}+\unicode[STIX]{x2202}_{y}v_{y}=qv_{z}$ . With this notation the truncated Navier–Stokes equations can be written as

where the overbar stands for vertical averaging $\overline{f}\equiv (1/H)\int f\,\text{d}z$ and $p_{q}$ is the partial pressure that guarantees the incompressibility of $\boldsymbol{u}_{q}$ . Note that if one plugs in $\boldsymbol{u}_{q}$ from (2.3) into (2.5), the vertical dependence drops out and we are left with two coupled partial differential equations that depend only on the horizontal directions $(x,y)$ . We note that our model is rigorously valid for very thin layers where $1/q$ is of the order of the viscous length scales, in which case it can be derived based on an asymptotic expansion, using the thinness of the layer as a small parameter. Despite the lack of rigor for thicker layers, as we will show later, this model is able to capture the relevant physical processes and display a rich behaviour.

Furthermore, note that in our system we have included the hypo-dissipation term $\unicode[STIX]{x1D707}\unicode[STIX]{x1D6E5}^{-2}\boldsymbol{u}_{2D}$ that is responsible for saturating the inverse cascade when present. If such a term is absent in the presence of an inverse cascade the energy of the large-scale modes will grow to very large values forming a condensate (Chertkov et al. Reference Chertkov, Connaughton, Kolokolov and Lebedev2007; Boffetta & Ecke Reference Boffetta and Ecke2012; Laurie et al. Reference Laurie, Boffetta, Falkovich, Kolokolov and Lebedev2014) whose amplitude growth is balanced by the small viscous forces. These two cases correspond to two very different situations. In the first, there is a constant flux of energy to the large scales while in the second situation the system reaches an almost thermal equilibrium with a non-constant and very weak (compared to the amplitude of the fluctuations) inverse flux of energy. Although both cases are physically relevant, in this work we only focus on the first situation and in all our runs $\unicode[STIX]{x1D707}$ is tuned so that the largest scale is sufficiently suppressed.

In the absence of forcing and dissipation terms this system conserves the total energy of the flow given by

where angular brackets $\langle \boldsymbol{\cdot }\rangle$ stand for spatial average. For $\boldsymbol{u}_{q}=0$ one recovers the two-dimensional Navier–Stokes equation in which case the enstrophy $\unicode[STIX]{x1D6FA}=\langle |\unicode[STIX]{x1D735}\times \boldsymbol{u}_{2D}|^{2}\rangle /2=\langle |\unicode[STIX]{x0394}\unicode[STIX]{x1D713}|^{2}\rangle /2$ is also conserved. In the presence of forcing and dissipation the system eventually reaches a steady state where the energy injected by the forcing at the averaged rate $I=\langle \boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{u}_{2D}\rangle _{T}$ is balanced by the dissipation rates $\unicode[STIX]{x1D716}_{3D}^{+}$ , $\unicode[STIX]{x1D716}_{2D}^{+}$ , and $\unicode[STIX]{x1D716}_{2D}^{-}$ due to viscous and hypo-viscous forces defined as:

where the brackets $\langle \boldsymbol{\cdot }\rangle _{T}$ indicate space and time average. The dissipation rate $\unicode[STIX]{x1D716}_{2D}^{+}$ measures the rate energy is dissipated at the small scales of the 2-D velocity field $\boldsymbol{u}_{2D}$ , while $\unicode[STIX]{x1D716}_{3D}^{+}$ measures the rate energy is dissipated by the 3-D field $\boldsymbol{u}_{q}$ . Finally $\unicode[STIX]{x1D716}_{2D}^{-}$ measures the rate energy is dissipated at the large scales of the 2-D velocity field. At steady state a balance is reached and we obtain:

Thus, with a large enough scale separation $L\gg \ell _{f}$ , the ratio $\unicode[STIX]{x1D716}_{2D}^{-}/I$ provides us with a measure of what fraction of the energy injected is cascading to the large scales. On the other hand, at high enough Reynolds numbers, the ratios $\unicode[STIX]{x1D716}_{2D}^{+}/I$ and $\unicode[STIX]{x1D716}_{3D}^{+}/I$ provide us with the fraction of the energy that cascades to the small scales.

We need to stress here the importance of considering the large Reynolds number and large-box limit. Here by large-box limit we refer to the horizontal dimensions of the box $L\rightarrow \infty$ and not $H$ , which is considered fixed in this limiting procedure. In 2-D turbulence there is always some forward transfer of energy that is required to sustain the forward cascade of enstrophy. Similarly, in 3-D turbulence in the presence of a hypo-dissipation term there is an inverse transfer of energy that is needed to sustain the energy at large scales against any dissipative mechanism that removes energy from these scales. These forward and inverse transfers, however, depend on the amplitude of the dissipation coefficients and the box size. When the dissipation coefficients are decreased and the box size is increased, their amplitude decreases to zero (see discussions in Eyink (Reference Eyink1996), Boffetta & Musacchio (Reference Boffetta and Musacchio2010), Boffetta & Ecke (Reference Boffetta and Ecke2012)). Thus, in order to distinguish a weak inverse cascade from these finite size effects we need to demonstrate that $\unicode[STIX]{x1D716}_{2D}^{-}$ remains finite as the box size is increased and the hypo-dissipation is decreased. Similarly for $\unicode[STIX]{x1D716}_{2D}^{+}$ and $\unicode[STIX]{x1D716}_{3D}^{+}$ , we need to investigate their large Reynolds number behaviour.

Two different forcing functions were used in this study. The first was a deterministic time-independent forcing explicitly given by

where $F_{0}$ is the forcing amplitude and $k_{f}=\unicode[STIX]{x03C0}/\ell _{f}$ is the forcing wavenumber. It corresponds to forcing a square array of vortices with alternating sign of vorticity. This forcing, although more physical than our second choice, does not inject energy at a constant rate since the injection rate depends on the particular flow realization.

For a better control of the energy injection rate in our system a second type of forcing was used that was designed to inject energy at a given shell of wavenumbers of modulus $k_{f}$ at a constant rate at each instant of time and for every realization. It is written as

where $\tilde{\boldsymbol{u}}_{2D}(\boldsymbol{k})$ stands for the Fourier transform of the field $\boldsymbol{u}_{2D}$ . The sums are over all wavenumbers that satisfy $|\boldsymbol{k}|=k_{f}$ , and $\unicode[STIX]{x1D6FA}_{\boldsymbol{k}}$ is a random frequency that takes values between 1 and $-$ 1 (for $|\boldsymbol{k}|=k_{f}$ and zero otherwise) and leads to a decorrelation of the phases between the different forced modes. The total energy injection rate for this forcing at each instant of time and for each realization is given by $I_{0}$ . This forcing thus allows us to control the energy injection rate without employing a delta correlated random forcing.

The relevant non-dimensional control parameters of our system depend on the domain geometry, dissipation parameters, and on the forcing mechanism and scale. The ratio of the inverse layer thickness to the forcing wavenumber is given by $Q\equiv q/k_{f}$ and is our primary control parameter. The relative scale separation between the forcing scale and the horizontal box size is measured by $k_{f}L$ . The Reynolds number of our system is defined as, $Re=u_{f}/k_{f}\unicode[STIX]{x1D708}$ , where $u_{f}$ is the root mean square (r.m.s.) value of the velocity at the forcing scale

Similarly we define a Reynolds number based on $\unicode[STIX]{x1D707}$ that is related to the inverse cascade as $R_{\unicode[STIX]{x1D707}}=u_{f}k_{f}^{5}/\unicode[STIX]{x1D707}$ . The value of $\unicode[STIX]{x1D707}$ was always tuned so that no large-scale condensate was formed and thus it was always linked to the size of the box $k_{f}L$ . We note however that $u_{f}$ is a quantity that is measured a posteriori and does not truly express a control parameter. For this reason we also define $Re_{f}$ , which is a Reynolds number based on $F_{0}$ or $I_{0}$ , that represents a true control parameter of the given system. For the runs of constant forcing amplitude given in (2.9) we define $Re_{f}=F_{0}^{1/2}/\unicode[STIX]{x1D708}k_{f}^{3/2}$ and $R_{\unicode[STIX]{x1D707}f}=F_{0}^{1/2}k_{f}^{9/2}/\unicode[STIX]{x1D707}$ . For the constant energy injection rate runs of (2.10) we define $Re_{f}=I_{0}^{1/3}/\unicode[STIX]{x1D708}k_{f}^{4/3}$ and $R_{\unicode[STIX]{x1D707}f}=I_{0}^{1/3}k_{f}^{14/3}/\unicode[STIX]{x1D707}$ .

Equations (2.4) and (2.5) were solved using a standard parallel pseudo-spectral code with a fourth-order Runge–Kutta scheme for time integration and a two-thirds dealiasing rule. More details on the parallelization can be found in Gomez, Mininni & Dmitruk (Reference Gomez, Mininni and Dmitruk2005). All runs started from random initial conditions and were carried out long enough so that a statistically steady state was reached. All measurements and averages were made at this state unless otherwise stated. The resolutions used varied from $512^{2}$ to $4096^{2}$ grid points.

Our goal is to investigate how the system transitions to and from a bidirectional cascade as $Q$ is varied in the limit of large box size $L$ (or large-scale separation $L\gg 1/k_{f}$ ) and large Reynolds numbers. Thus, close to the points of criticality, a series of runs were performed varying the parameter $Q$ for fixed values of $k_{f}L$ and $Re_{f}$ . Then the same series of runs were repeated for larger values of $k_{f}L$ or $Re_{f}$ until we observed convergence. The runs were separated into three cases: A, B, and C. See table 1 for a list of our cases and their corresponding parameters. Cases A and B investigated the transition from a purely forward cascade to a bidirectional cascade, with a constant forcing amplitude (2.9) and a constant energy injection rate (2.10), respectively. For these runs we fixed $Re_{f}$ while varying $k_{f}L$ , and measured the inverse energy cascade (which we expected to transition from being zero to being non-zero) via $\unicode[STIX]{x1D716}_{2D}^{-}/I$ . Finally, Case C investigated the transition from a bidirectional cascade to a purely inverse energy cascade. For these runs we fixed $k_{f}L$ , varied $Re_{f}$ , and measured the forward energy cascade (which we expected to transition from being non-zero to zero) via $\unicode[STIX]{x1D716}_{3D}^{+}/I$ .

Figure 2. (a) Runs from Case A. Measure of the inverse cascade for various runs as we increased the box size $k_{f}L$ . Our runs converge to a sharp transition in the large box limit, confirming the existence of $Q_{3D}$ . The dashed vertical line represents the approximate value of $Q_{3D}$ . (b) Runs from Case B. A similar tendency is seen, although convergence is not yet reached.

3 Transition from a forward to a bidirectional cascade

We first investigate the first critical point $Q_{3D}$ that marks the transition from a forward cascade to bidirectional cascade. Figure 2 shows the dependence of the inverse cascade measured by the ratio of $\unicode[STIX]{x1D716}_{2D}^{-}/I$ as a function of $Q$ for flows driven by the deterministic forcing (2.9) in panel (a) and by constant injection of energy forcing (2.10) in panel (b). Since both cases are similar, we focus on panel (a), Case A. The three lines correspond to three different box sizes $k_{f}L$ as marked in the legend. For the smallest box size $k_{f}L=8$ the transition appears smooth with the presence of an inverse cascade at even the smallest values of $Q$ displayed. As the box size is increased the transition appears to become sharper, and at the largest box size the system seems to converge into a critical transition at the value $Q=Q_{3D}\simeq 3.6$ . Close to the critical point the inverse cascade $\unicode[STIX]{x1D716}_{2D}^{-}/I$ appears to scale linearly with the deviation from criticality $\unicode[STIX]{x1D716}_{2D}^{-}\propto (Q-Q_{3D})I$ for $Q>Q_{3D}$ . Looking at Case B in figure 2(b), the transition seems to take longer to converge but the same tendency is seen. For large values of $Q$ the curves saturate to a value slightly smaller than unity, which reflects that there is still a finite value of energy that is dissipated by viscous effects. This is because only moderate values of $Re$ are used. We discuss this issue more in the next section where the second critical point is investigated and $Re$ is varied.

An alternative way to quantify the rate at which energy cascades towards the large scales is by measuring the rate of increase of energy. An inverse cascading system without large-scale dissipation and with a constant energy injection rate is expected to lead to a linear increase of energy with time. This linear increase is expected after a short transient time where small scales reach a quasi-steady state and before the largest scales of the system are exited and a condensate starts to form. The growth rate of the energy (i.e. the linear slope) is due to the inverse cascade and is equal to the inverse energy flux. We tested this method of measurement by performing runs identical to those of Case B (using the forcing with constant energy injection given in (2.10)) but this time without a hypo-dissipation term, $\unicode[STIX]{x1D707}=0$ . The energy evolution for the runs without hypo-dissipation are depicted in figure 3(a) for different values of $Q$ , where a linear increase of energy can be seen. The slope of the linear growth of energy for different values of $Q$ and $k_{f}L$ is calculated. These results are compared with the results obtained from steady-state runs in the presence of hypo-dissipation and are shown in figure 3(b) for three different values of the box size. For the smallest value of $k_{f}L$ the two measurements differ quite a bit, with the steady-state hypo-viscous simulations (dashed) which display a much more smooth behaviour, and the $\unicode[STIX]{x1D707}=0$ runs (solid) having significantly larger error. As the box size is increased the two methods of measuring the amplitude of the inverse cascade converge and the two curves overlap. Although measuring the inverse cascade by the slope of the energy versus time graph converges faster to the large box limit, the calculation using the hypo-dissipation term leads to much smaller error bars due to the long-time averaging that is possible in this case. The smaller error is favourable when attempting to determine very small trends of inverse cascade, as is seen figure 2.

Figure 3. (a) The evolution of the energy for different values of $Q$ in the absence of large-scale dissipation. (b) The results of our constant energy injection runs, comparing the two methods of measuring the inverse cascade. The solid lines represent the slope measurement with no hypo-dissipation, whereas the dashed lines represent the steady-state long-time-averaged measurements of runs with hypo-dissipation (Case B). The plots are shifted to emphasize differences. The horizontal black lines all correspond to $\unicode[STIX]{x1D716}_{2D}^{-}=0$ for the corresponding plot.

We next look how the spectral distribution of energy is changing as $Q$ is varied by looking at the energy spectra $E(k),E_{2D}(k),$ and $E_{q}(k)$ defined as:

where $\tilde{\boldsymbol{u}}_{2D}(\boldsymbol{k})$ and $\tilde{\boldsymbol{v}}_{q}(\boldsymbol{k})$ stand for the Fourier mode amplitude of wavenumber $\boldsymbol{k}$ of the $\boldsymbol{u}_{2D}$ and $\boldsymbol{v}_{q}$ fields respectively. The total energy spectrum $E(k)=E_{2D}(k)+E_{q}(k)$ is shown in figure 4(a) for different values of $Q$ and for $k_{f}L=32$ . These spectra were outputted during the run time of our simulations and then afterwards averaged together during the steady-state period to get the final results. The smaller values of $Q$ are displayed with darker colours while the larger values of $Q$ are displayed with lighter colours. Clearly as $Q$ is increased there is more energy in the large scales and less energy in the small scales. For small values of $Q$ the energy spectrum is close to flat at the large scales $k<k_{f}$ and at the small scales $k>k_{f}$ it is compatible with a $k^{-5/3}$ power-law scaling. As we increase $Q$ the large scales gain more and more energy until we reach the point where a spectrum compatible with $k^{-5/3}$ is observed at the large scales $k<k_{f}$ , and at small scales the spectrum is very steep (steeper than $k^{-3}$ that is predicted for 2-D turbulence). This very steep behaviour is due to the small value of the Reynolds number used. The $k^{-3}$ spectrum is only realized in very high $Re$ (see Boffetta & Musacchio Reference Boffetta and Musacchio2010), and such values are beyond the scope of this work. Figure 4(b) compares the three spectra $E_{q}$ (blue, vertical dashed line) and $E_{2D}$ (red, dashed line) for the intermediate value of $Q=3.44>Q_{3D}$ . At large wavenumbers the energy spectrum is dominated by $E_{q}$ , while at small wavenumbers, which exhibit some inverse cascade, the spectrum is dominated by $E_{2D}$ .

Figure 4. (a) Spectra of the total energy for different values of $Q$ (averaged over many outputs). The colour of the curves become lighter and lighter as we have more and more inverse cascade. (b) The total energy $E(k)$ (black), the 2-D energy $E_{2D}$ (red) and the 3-D energy $E_{q}$ (blue) for $Q=3.44$ .

The direction of the cascade is best described by the direct measurement of the transfer of energy among scales provided by the nonlinear terms in our equations. More precisely, this flux of energy expresses the rate energy is transferred out of a given set of wavenumbers due to the nonlinearities. By performing a filtering of the 2-D velocity field $\boldsymbol{u}_{2D}$ and the 3-D velocity field $\boldsymbol{u}_{q}$ in Fourier space so that only the wavenumbers with modulus smaller than $k$ are kept, denoted by $\boldsymbol{u}_{2D}^{{<}k}$ and $\boldsymbol{u}_{q}^{{<}k}$ , one looks at only the structures of scales larger than $\ell =2\unicode[STIX]{x03C0}/k$ . With these filtered velocity one can calculate the flux of energy, defined to be:

The first flux $\unicode[STIX]{x1D6F1}_{2D}$ expresses the rate energy is transferred from large-scale $\boldsymbol{u}_{2D}$ modes to small-scale $\boldsymbol{u}_{2D}$ modes by self-interaction. The second flux $\unicode[STIX]{x1D6F1}_{q}$ expresses the rate energy is transferred from large-scale $\boldsymbol{u}_{q}$ modes to small-scale $\boldsymbol{u}_{q}$ modes through interactions with the 2-D field $\boldsymbol{u}_{2D}$ . Finally, the last flux $\unicode[STIX]{x1D6F1}_{T}$ expresses the rate energy is transferred to the small scales by a simultaneous exchange of energy from one field to the other. The total energy flux is given by

In the inertial range $\unicode[STIX]{x1D6F1}_{E}(k)$ is constant and positive if the cascade is forward and constant and negative if the cascade is inverse. In the case for which $\boldsymbol{u}_{q}=0$ and enstrophy is conserved we can also define the flux of enstrophy as

We note however that if $\boldsymbol{u}_{q}$ is not exactly zero the enstrophy flux $\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FA}}(k)$ is not constant due to the enstrophy generation by the vorticity stretching term $\overline{\boldsymbol{u}_{q}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{q}}$ .

Figure 5 shows the time-averaged total energy flux $\unicode[STIX]{x1D6F1}_{E}$ normalized by the energy injection rate $I$ for different values of $Q$ . As before, dark lines correspond to small values of $Q$ while light lines correspond to larger values of $Q$ . For the smallest value of $Q$ the flux $\unicode[STIX]{x1D6F1}_{E}$ is almost zero for $k<k_{f}$ and positive for $k>k_{f}$ describing a unidirectional forward cascade. As $Q$ is increased a bidirectional cascade appears with both positive and negative fluxes in either side of the forcing wavenumber indicating that energy cascades to both large and small scales. At the largest value of $Q$ almost all energy cascades inversely demonstrated by the negative values of the flux for $k<k_{f}$ and almost zero values for $k>k_{f}$ . The time-averaged flux shows the predicted constant-in- $k$ behaviour for small $k$ while it is affected by viscosity in the large scales and thus decays with $k$ . Instantaneous fluxes however are strongly fluctuating. This is displayed in figure 5(b), where the time-averaged value of the flux (thick black line) is displayed along with numerous instantaneous fluxes (cyan lines) for the case of $Q=3.75$ . The role of the fluctuations becomes particular important close to the transition where the amplitude of the fluctuations becomes much larger than the mean value.

Figure 5. (a) The total energy flux $\unicode[STIX]{x1D6F1}_{E}$ for different values of $Q$ (averaged over many outputs). The lighter the shade of green, the higher the value of $Q$ . (b) For the case of $Q=3.75$ : The cyan/light blue lines show the instantaneous values of the total energy flux; the solid black line is their average, corresponding to one of the curves on (a).

To demonstrate the role of each field in the cascade of energy we show, in figure 6(a), the decomposition of the total energy flux $\unicode[STIX]{x1D6F1}_{E}$ into the three components $\unicode[STIX]{x1D6F1}_{2D},\unicode[STIX]{x1D6F1}_{q}$ and $\unicode[STIX]{x1D6F1}_{T}$ for the value of $Q=3.75$ . The flux $\unicode[STIX]{x1D6F1}_{2D}$ is negative for all values of $k<k_{f}$ while it is positive but small for $k>k_{f}$ . The remaining fluxes $\unicode[STIX]{x1D6F1}_{q}$ and $\unicode[STIX]{x1D6F1}_{T}$ are positive for all $k$ . Thus the inverse cascade is driven by the $\unicode[STIX]{x1D6F1}_{2D}$ term whereas the forward cascade is driven mostly by $\unicode[STIX]{x1D6F1}_{q}$ and $\unicode[STIX]{x1D6F1}_{T}$ . It is worth noting that neither of these partial fluxes is constant in either inertial range. Furthermore $\unicode[STIX]{x1D6F1}_{q}$ and $\unicode[STIX]{x1D6F1}_{T}$ are positive and finite for $k<k_{f}$ implying that part of the energy transferred to the large scales by the $\boldsymbol{u}_{2D}$ field is brought back to the small scales by interactions with the $\boldsymbol{u}_{q}$ field. As $k$ is decreased the effect of $\unicode[STIX]{x1D6F1}_{q}$ and $\unicode[STIX]{x1D6F1}_{T}$ weakens and at the smallest $k$ only $\unicode[STIX]{x1D6F1}_{2D}$ remains non-zero.

Figure 6(b) shows the flux of enstrophy $\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FA}}$ for different values of $Q$ . As discussed in the previous section enstrophy is only conserved by the nonlinearities when $\boldsymbol{u}_{q}=0$ . For large values of $Q$ , for which $\boldsymbol{u}_{q}$ is small, we expect that enstrophy will be quasi-conserved and its flux will be positive and slowly varying. This is indeed observed in this figure for $Q=8.75$ , where the enstrophy flux is slowly decreasing due to viscous effects. For smaller values of $Q$ , however, the enstrophy flux becomes non-monotonic with a sharp increase close to the dissipation scales. This is due to the generation of enstrophy by the $\boldsymbol{u}_{q}$ field that leads to an excess of enstrophy that is transported to the small scales.

Figure 6. (a) For the case of $Q=3.75$ : the total energy flux is shown in the thick black line and the contributions from each nonlinear term are shown in the other coloured lines. Notice that $\unicode[STIX]{x1D6F1}_{2D}$ is the only one with a negative contribution. One should also note that the green dotted line, corresponding to the term $\unicode[STIX]{x1D6F1}_{q}$ , is always positive for all $Q$ examined but that the blue dash-dot line, corresponding to $\unicode[STIX]{x1D6F1}_{T}$ , can be both negative and positive depending on $Q$ and $k$ . (b) The enstrophy flux $\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FA}}$ for different values of $Q$ (averaged over many outputs) normalized by the enstrophy injection rate $I_{\unicode[STIX]{x1D6FA}}=\langle \unicode[STIX]{x1D735}\times \boldsymbol{u}_{2D}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{F}\rangle _{T}$ . The colour of the curves becomes lighter and lighter as the value of $Q$ increases.

The results so far have demonstrated the presence of a bidirectional cascade in Fourier space. Further insight is obtained by looking at the resulting structures in physical space. Figure 7 shows, in (a–c), the 2-D vorticity $\unicode[STIX]{x1D735}\times \boldsymbol{u}_{2D}$ for three values of $Q$ . Panels (d–f) show the energy density ${\mathcal{E}}_{q}\equiv |\boldsymbol{v}_{q}|^{2}/4$ of the 3-D field at the same instant of time as the vorticity for the same values of $Q$ . The vorticity snapshots show that, as $Q$ is varied from small to large values, the structures change from small scale to structures larger than the forcing scale $k_{f}L=16$ , further demonstrating the growth of the inverse cascade as $Q$ is increased. For the smallest value of $Q$ the structures are small compared to the forcing scale, and one can notice by the sharp changes from blue to red that there are also large gradients of vorticity. For the largest value of $Q$ the flow has very similar structure to that obtained in pure 2-D simulations, where structures grow in size due to vortex clustering and stretching of the iso-vorticity lines (Chen et al. Reference Chen, Ecke, Eyink, Rivera, Wan and Xiao2006; Eyink Reference Eyink2006; Xiao et al. Reference Xiao, Wan, Chen and Eyink2009; Boffetta & Musacchio Reference Boffetta and Musacchio2010; Boffetta & Ecke Reference Boffetta and Ecke2012). Finally, the intermediate value of $Q$ appears to have characteristics of both extreme cases but concentrated in different regions of space. Some regions in space resemble the small-scale eddies and large gradient case of small $Q$ while other regions have the typical 2-D eddy structures. The structures in the 3-D energy density snapshots, shown in figure 7(d–f) are filamentary and are well correlated with the areas of high strain in the 2-D field. Furthermore, the density of these areas of high ${\mathcal{E}}_{q}$ appears to decrease with $Q$ . Combining these observations, one notices that, for small values of $Q$ , the intensity of vorticity is not uniform throughout the domain – some regions appear to be ‘3-D active’ and some regions to be ‘3-D quiet’. As $Q$ is increased the ‘3-D active’ regions appear to become less space filling. A similar behaviour has been observed in the transition from forward to inverse cascade in 2-D magnetohydrodynamic flows, where the role of $\boldsymbol{u}_{q}$ was played by the magnetic field (Seshasayanan et al. Reference Seshasayanan, Benavides and Alexakis2014; Seshasayanan & Alexakis Reference Seshasayanan and Alexakis2016). This behaviour, as we will discuss next, has direct consequences for the behaviour of the flow close to the second critical point.

Figure 7. The 2-D vorticity $\unicode[STIX]{x1D735}\times \boldsymbol{u}_{2D}$ (a–c) and 3-D energy density ${\mathcal{E}}_{q}$ (d–f), for three values of $Q$ . The box size $k_{f}L=16$ is fixed for all cases. The grid pattern which is seen in (a) is due to the forcing.

4 Transition from a bidirectional to an inverse cascade

We now focus on large values of $Q$ for which the flow transitions from a bidirectional cascade to an inverse unidirectional cascade. As shown in the previous section the terms that involve the 3-D field $\boldsymbol{u}_{q}$ are responsible for driving the forward cascade. In their absence we recover the 2-D Navier–Stokes equation that leads to an inverse cascade with all energy dissipated at large scales in the large $Re$ limit. Thus in order to transition to a unidirectional inverse cascade the flow needs to become purely two-dimensional with $\boldsymbol{u}_{q}=0$ . This is indeed possible since a purely 2-D flow is always a solution of the 3-D Navier–Stokes equation, equation (2.1). Therefore if the initial data are exactly two-dimensional the flow will remain two-dimensional for all times cascading energy inversely. However these solutions can be unstable and small perturbations can grow exponentially driving the flow away from the 2-D behaviour and thus alter the direction of cascade. This brings us to the intuitive understanding of the second critical point. The bidirectional cascade will transition to a unidirectional inverse cascade when the layer thickness is so small ( $Q$ is large enough) that all 3-D perturbations are damped and decay exponentially in time, and thus the 2-D solution is an attractor of the system. This is expected to occur when the thickness of the box $H_{2D}$ is such that the viscous dumping rate $\unicode[STIX]{x1D708}/H_{2D}^{2}$ due to the vertical variation alone is similar to the shear rate $u_{f}k_{f}$ that drives the 3-D instability. This argument $\unicode[STIX]{x1D708}/H_{2D}^{2}\sim u_{f}k_{f}$ implies that the second critical point $Q_{2D}$ will satisfy:

A more precise estimate for the location of $Q_{2D}$ can be obtained by looking at the energy evolution of the 3-D component of the flow that reads

Note that the equation above is also valid for the full Navier–Stokes equations, not just our model, with $\boldsymbol{u}_{q}=\boldsymbol{u}-\boldsymbol{u}_{2D}$ . The 2-D solutions will be globally stable if the functional ${\mathcal{H}}$ is negative definite in which case all 3-D perturbations will decay independent of the initial conditions. We can prove that this is the case for a given range of $Q$ using the following rigorous inequalities. Using Hölder’s inequality the vortex stretching term $\langle \boldsymbol{u}_{q}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{u}_{2D})\boldsymbol{\cdot }\boldsymbol{u}_{q}\rangle$ is bounded by

where $\Vert \unicode[STIX]{x1D735}\boldsymbol{u}_{2D}\Vert _{\infty }$ stands for the $L_{\infty }$ norm (both in time and space) of the gradient of $\boldsymbol{u}_{2D}$ . Furthermore Poincaré’s inequality gives, for the dissipation term, that $\langle |\unicode[STIX]{x1D735}\boldsymbol{u}_{q}|^{2}\rangle \geqslant \unicode[STIX]{x1D708}q^{2}\langle |\boldsymbol{u}_{q}|^{2}\rangle$ . These two results lead to the negativity of ${\mathcal{H}}$ if $Q^{2}=q^{2}/k_{f}^{2}\geqslant \Vert \unicode[STIX]{x1D735}\boldsymbol{u}_{2D}\Vert _{\infty }/k_{f}^{2}\unicode[STIX]{x1D708}$ . Provided that $\Vert \unicode[STIX]{x1D735}\boldsymbol{u}_{2D}\Vert _{\infty }$ is bounded, this result guarantees that, beyond some value of $Q$ , the flow becomes exactly two-dimensional. However, deriving an upper bound on $\Vert \unicode[STIX]{x1D735}\boldsymbol{u}_{2D}\Vert _{\infty }$ in terms of the control parameters of the system is not trivial. An upper bound on the long-time-averaged growth rate of $\langle |\boldsymbol{u}_{q}|^{2}\rangle$ was derived in Gallet & Doering (Reference Gallet and Doering2015) where the two-dimensionalization of a flow due to the presence of an external strong magnetic field $B_{0}$ was examined. Their results extend directly to our case by setting $B_{0}=0$ . They were able to show, using properties of the steady-state 2-D Navier–Stokes equations (Alexakis & Doering Reference Alexakis and Doering2006), that the 3-D energy will decay, in the long-time limit if $Q$ satisfies:

where each $c_{i}$ is a positive dimensionless number that depends on forcing shape and other non-universal properties of the flow and $Re$ is based on the r.m.s. velocity of the flow. Note that $Q$ does not scale as we predicted with respect to $Re$ . This is because this bound is very conservative and does not capture all the physics at the transition point. However, despite this loose scaling, this results guarantees that the transition from a bidirectional to a purely inverse cascade will be through a critical point because it guarantees that there is a value of $Q$ above which $\unicode[STIX]{x1D716}_{3D}^{+}$ is exactly zero.

Another less conservative bound which was derived in Gallet & Doering (Reference Gallet and Doering2015) is based on the linear stability analysis of 2-D flows to 3-D perturbations. Linear stability guarantees that the 2-D solutions are locally stable but does not exclude the possibility that a locally attracting 3-D solution also exist for the same parameters. Their analysis shows that 2-D solutions (for harmonic forcing and in the absence of hypo-dissipation term) are linearly stable provided that

Up to logarithmic corrections this result follows the scaling $Q_{2D}\sim Re^{1/2}$ , which is what our scale analysis predicts.

Figure 8. Run set C. $\unicode[STIX]{x1D716}_{3D}^{+}$ normalized by the energy injection rate for various runs as we increased the Reynolds number. Rescaling $Q$ by our scale analysis prediction collapses the curves into one. The dashed vertical line represents the approximate value of $Q_{2D}$ , after which all runs have decaying 3-D energy.

This scaling is clearly demonstrated by our numerical simulations in figure 8 that shows the energy dissipation rate $\unicode[STIX]{x1D716}_{3D}^{+}$ as a function of the rescaled $Q/\sqrt{Re}$ . This rescaling of $Q$ collapses the curves together. The vertical dashed line marks the critical value $Q_{2D}=C_{2D}\sqrt{Re}$ (where $C_{2D}\simeq 1.16$ ) and represents the point after which the amplitude of the 3-D field decays exponentially. In the absence of $\boldsymbol{u}_{q}$ we have a 2-D solution – $\unicode[STIX]{x1D716}_{3D}^{+}$ is zero and most of the energy is dissipated in the large scales. Note that for a given value of $Q$ the exact two-dimensionalization occurs for the values of $Re$ in the range $Re\leqslant Q^{2}/C_{2D}^{2}$ . Since $Re$ is bounded from above there is always some viscous dissipation $\unicode[STIX]{x1D716}_{2D}^{+}$ . Thus, a strict inverse cascade where all the injected energy is dissipated at the large scales $\unicode[STIX]{x1D716}_{2D}^{-}=I_{0}$ and $\unicode[STIX]{x1D716}_{2D}^{+}+\unicode[STIX]{x1D716}_{3D}^{+}=0$ is realized only in the limit $Re\rightarrow \infty$ , $Q\rightarrow \infty$ with $Q\geqslant Q_{2D}=C_{2D}Re^{1/2}$ . On the other hand the small-scale dissipation due to the $\boldsymbol{u}_{q}$ field $\unicode[STIX]{x1D716}_{3D}^{+}$ becomes exactly zero $\unicode[STIX]{x1D716}_{3D}^{+}=0$ for any finite value of $Re$ that satisfies $Re\leqslant Q^{2}/C_{2D}^{2}$ . It is also worth pointing out that the dependence of $\unicode[STIX]{x1D716}_{3D}^{+}$ on the deviation from criticality differs from that of $\unicode[STIX]{x1D716}_{2D}^{-}$ close to the first critical point $Q_{3D}$ that was linear. Close to $Q_{2D}$ the energy dissipation $\unicode[STIX]{x1D716}_{3D}^{+}$ follows the relation

The exponent $\unicode[STIX]{x1D6FD}$ was found to be close to $\unicode[STIX]{x1D6FD}\simeq 2$ . To understand the origin of this exponent we need to look at the temporal and spatial form of the unstable field $\boldsymbol{u}_{q}$ .

Figure 9. A typical time series of the 3-D energy for a $Q$ very close to $Q_{2D}$ . (a) Shows linear scalings and the typical picture of an on–off intermittent signal, whereas (b) shows the same signal but in log–linear scale.

Figure 9 shows a typical signal for the evolution of the energy of the 3-D field $\boldsymbol{u}_{q}$ as a function of time, in linear (a) and log-linear (b) scale. The signal exhibits bursts of energy followed by times with very weak energy. This is a signature of ‘on–off’ intermittency. On–off intermittency is a generic behaviour that appears in the vicinity of an instability in the presence of multiplicative noise (Fujisaka et al. Reference Fujisaka, Ishii, Inoue and Yamada1986; Platt, Spiegel & Tresser Reference Platt, Spiegel and Tresser1993). It has been observed in turbulent dynamo simulations (Alexakis & Ponty Reference Alexakis and Ponty2008), in electronic circuits (Hammer et al. Reference Hammer, Platt, Hammel, Heagy and Lee1994), in electrohydrodynamic convection (John, Stannarius & Behn Reference John, Stannarius and Behn1999) and spin-wave systems (Rödelsperger, Čenys & Benner Reference Rödelsperger, Čenys and Benner1995). In such situations the unstable modes have a growth rate that varies strongly with time, taking both positive and negative values. Our system also falls in this category because for small values of $\boldsymbol{u}_{q}$ it is linearly advected by a turbulent/chaotic flow $\boldsymbol{u}_{2D}$ . Thus, at some instances of time $\boldsymbol{u}_{2D}$ is efficient at exponentially increasing $\boldsymbol{u}_{q}$ while at others $\boldsymbol{u}_{q}$ is exponentially decreased. In its simplest form on–off intermittency models a single mode $X$ that is amplified or decreased by a fluctuating growth rate $\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D709}$ and is restricted from taking very large values by a nonlinearity $-X^{3}$ :

Here $\unicode[STIX]{x1D6FC}$ represents mean growth rate and $\unicode[STIX]{x1D709}$ represents a zero mean random noise. In the thin layer turbulent system the 3-D instabilities of the 2-D flow, whose energy is described by (4.2), have an averaged growth rate $\unicode[STIX]{x1D6FC}$ that is proportional to the deviation from criticality $\unicode[STIX]{x1D6FC}\propto (Q_{2D}-Q)$ . At the same time the role of the turbulent fluctuations of the 2-D field $\boldsymbol{u}_{2D}$ is played by the random multiplicative noise. If the averaged growth rate $\unicode[STIX]{x1D6FC}$ in (4.7) is sufficiently smaller than the fluctuations $\unicode[STIX]{x1D709}$ of the instantaneous growth rate, then the system spends long intervals of time with very small amplitudes (off phase) intervened by short burst where the nonlinearities are effective (on phases). In log scale the amplitude of the unstable mode follows a biased random walk bounded from above by the nonlinearities. The model predicts that, close to the onset, the probability distribution function (PDF) $P(E_{q})$ of the energy of the mode $E_{q}\sim X^{2}$ follows the scaling

where $D$ is the amplitude of the noise. Furthermore, the duration of the off phases $T_{off}$ diverges with the deviation from the onset (here $Q_{2D}-Q$ ) as $T_{off}\propto (Q_{2D}-Q)^{-1}$ while the amplitude $E_{on}$ and the time $T_{on}$ in the on phase becomes independent of ( $Q_{2D}-Q$ ). This implies that the time-averaged energy scales like

where $\langle {\mathcal{E}}_{q}\rangle _{T}$ is also the time average of $E_{q}$ . Note that this is not the behaviour we observe in figure 8, which suggests a scaling closer to $\langle {\mathcal{E}}_{q}\rangle _{T}\propto (Q_{2D}-Q)^{2}$ , since the time averaged $\unicode[STIX]{x1D716}_{3D}^{+}$ is proportional to the time average of $E_{q}$ . Figure 10(a) shows the PDF of the energy $E_{q}$ calculated for different values of $Q$ . In agreement with the model when the critical value $Q_{2D}$ is approached, the PDF becomes singular showing a power-law behaviour $P(E_{q})\propto (E_{q})^{S(Q)}$ . The exponents $S(Q)$ of this power law are shown in figure 10(b). As criticality is approached $Q\rightarrow Q_{2D}$ these exponents tend to $-1$ in agreement again with the model. However this asymptotic value is not approached linearly (i.e. $S(Q)\simeq (Q_{2D}-Q)/D-1$ ) as the model suggests but closer to a quadratic behaviour $S(Q)\simeq (Q_{2D}-Q)^{2}/D^{2}-1$ . Thus, just as with the scaling observed in (4.6) there is a disagreement with this model.

Figure 10. (a) The PDF in time of $E_{q}$ . As before, dark means farther away from $Q_{2D}$ and thus smaller $Q$ and light green is the most two-dimensional (closest to $Q_{2D}$ ). The blue dashed line represents $E_{q}^{-1}$ and the red dashed lines represent where the power-law exponent measurements were taken (if any). (b) Various values of the exponents $S(Q)$ versus the fraction $Q/Q_{2D}$ that measures the deviation from criticality. These values are the slopes of the red dashed lines in (a). Notice the approach to $-1$ , as the model predicts.

Figure 11. The vorticity $\unicode[STIX]{x1D735}\times \boldsymbol{u}_{2D}$ (a) and the 3-D energy density ${\mathcal{E}}_{q}$ (b) for a point $Q/Q_{2D}=0.92$ , right before $Q_{2D}$ . For $k_{f}L=8$ and constant forcing amplitude. Note the very localized structure of the 3-D energy, an extreme of a trend seen in figure 7.

Resolution comes from looking at the spatial behaviour of the unstable modes. Figure 11 shows a snapshot of the vorticity $\unicode[STIX]{x1D735}\times \boldsymbol{u}_{2D}$ in panel (a) and of the 3-D energy density ${\mathcal{E}}_{q}=|\boldsymbol{v}_{q}|^{2}/4$ in panel (b). The vorticity shows the classical behaviour of 2-D turbulence. The 3-D energy, however, shows a very intermittent behaviour in space: most of the 3-D energy is concentrated in a single structure occupying a small fraction of the area of the domain size. Comparing with figure 7(d–f) we see that, as $Q$ is increased and $Q_{2D}$ is approached, there are less and less 3-D structures occupying a smaller and smaller fraction of the domain area. For $Q$ really close to $Q_{2D}$ in figure 11(b) only a single structure or two exist. The unstable solution thus is not only intermittent in time, but it is also intermittent in space!

Figure 12. (a) The PDF in space of ${\mathcal{E}}_{q}$ . As before, dark means farther away from $Q_{2D}$ and thus smaller $Q$ and light green is the most two-dimensional (closest to $Q_{2D}$ ). The blue dashed line represents $E_{q}^{-1}$ and the red dashed lines represent where the measurements of the exponents were taken. (b) Various values of the measured power-law exponents $Z$ versus the fraction $Q/Q_{2D}$ . The measured power laws are demonstrated by the red dashed lines in (a).

Figure 12 shows the PDF in space of ${\mathcal{E}}_{q}$ measured for different values of $Q$ . Just like the PDF of the energy $E_{q}$ the PDF of ${\mathcal{E}}_{q}$ shows a power-law behaviour: as criticality is approached $P({\mathcal{E}}_{q})\propto {\mathcal{E}}_{q}^{Z}$ with $Z\rightarrow -1$ as $Q\rightarrow Q_{2D}$ . Spatial intermittency is not taken into account in the on–off intermittency model (4.7) that takes into account the evolution of only one single mode. In our system the nonlinearity (on phase) is not only visited rarely in time but also rarely in space. The averaged energy in space and time will then satisfy

where $V_{on}$ is the area occupied by the ‘3-D active’ regions, and $V$ is the total area of our system. The scaling of $\langle {\mathcal{E}}_{q}^{n}\rangle _{T}$ on $(Q_{2D}-Q)$ will thus not only depend on the scaling of the time fraction $T_{on}/T_{off}$ but also on the area fraction $V_{on}/V_{off}$ . The model in (4.7) is not sufficient to describe our system and most likely an extended system with random multiplicative noise both in space and time will be required to capture correctly the statistics of our system (Grinstein, Munoz & Tu Reference Grinstein, Munoz and Tu1996; Horsthemke & Lefever Reference Horsthemke and Lefever2006). Such a possibility will be examined in future work.