Condensates in thin-layer turbulence

Abstract

We examine the steady state of turbulent flows in thin layers using direct numerical simulations. It is shown that when the layer thickness is smaller than a critical height, an inverse cascade arises which leads to the formation of a steady state condensate where most of the energy is concentrated in the largest scale of the system. For layers of thickness smaller than a second critical height, the flow at steady state becomes exactly two-dimensional. The amplitude of the condensate is studied as a function of layer thickness and Reynolds number. Bi-stability and intermittent bursts are found close to the two critical points. The results are interpreted based on a mean-field three-scale model that reproduces some of the basic features of the numerical results.

1 Introduction

Turbulent flows in geophysical and astrophysical contexts are often subject to geometrical constraints, such as thinness in a particular direction, that can strongly affect the behaviour of the flow. This occurs, for instance, in planetary atmospheres and oceans (Pedlosky 2013) whose behaviour can strongly deviate from the classical three-dimensional homogeneous and isotropic turbulence. This is related to the well-known fact that the behaviour of flows at large Reynolds numbers $Re$ depends on the dimensionality of the system. In three dimensions, vortex stretching transfers energy to small scales in a direct cascade (Frisch 1995). By contrast, in two dimensions, the conservation of enstrophy in addition to energy gives rise to an inverse energy cascade, a transfer of energy to the large scales (Boffetta & Ecke 2012). Flows in thin layers display properties of both systems, with the large scales behaving like a two-dimensional (2-D) flow and the small scales behaving like a 3-D flow. As a result, such systems are known to cascade energy both to large and to small scales (Smith, Chasnov & Waleffe 1996). In fact, it has been shown in (Celani, Musacchio & Vincenzi 2010; Benavides & Alexakis 2017; Musacchio & Boffetta 2017) that as the height of the layer $H$ is varied, the system transitions from a state where energy cascades only to the small scales for large $H$ , to a state where energy cascades to both large and small scales when $H$ is smaller than approximately half the size of the forcing length scale $\ell$ . In particular Benavides & Alexakis (2017), using a Galerkin truncated model of the full Navier–Stokes equations, were able to provide strong evidence of the criticality of the transition. In addition, they observed a second transition to exact two-dimensionalisation for layers of very small thickness $H\propto \ell Re^{-1/2}$ . This transition had been predicted theoretically using bounding techniques by Gallet & Doering (2015). Similar transitions from a strictly forward cascade to an inverse cascade have been observed in other systems like rotating turbulence (Deusebio et al. 2014), stratified turbulence (Sozza et al. 2015), rotating and stratified flows (Marino, Pouquet & Rosenberg 2015), magneto-hydrodynamic systems (Alexakis 2011; Seshasayanan, Benavides & Alexakis 2014; Seshasayanan & Alexakis 2016) and helically constrained flows (Sahoo & Biferale 2015; Sahoo, Alexakis & Biferale 2017), to mention a few (see Alexakis & Biferale (2018) for a review).

The thin layer, however, remains possibly the simplest model exhibiting such transitions and it thus deserves a detailed study at the different stages of inverse cascade evolution. In the presence of an inverse cascade, for finite systems and in the absence of a large-scale dissipation term, there are two stages in the development of the flow. In the first stage (at early times), energy is transferred to larger and larger scales by the inverse cascade. This process stops, however, when scales comparable to the system size are reached, after which energy starts to pile up at these largest scales. In the long-time limit, the increase of the large-scale energy saturates and a condensate is formed, where nearly all energy is found in the first few Fourier modes. For 2-D Navier–Stokes turbulence, the possibility of such a condensation phenomenon was first conjectured in the seminal paper of Kraichnan (1967), first seen in direct numerical simulations by Hossain, Matthaeus & Montgomery (1983), further explored quantitatively by Smith & Yakhot (1993, 1994) and, more recently, by Chertkov et al. (2007), Bouchet & Simonnet (2009), Chan, Mitra & Brandenburg (2012), Frishman, Laurie & Falkovich (2017), Frishman & Herbert (2018). Spectral condensation has also been studied in other quasi-2-D systems such as quasi-geostrophic flows (see Vallis & Maltrud 1993; Kukharkin, Orszag & Yakhot 1995; Kukharkin & Orszag 1996; Venaille & Bouchet 2011). In terms of the real space flow field, this spectral condensation corresponds to coherent system-size vortices or shear layers. In two dimensions, where the cascade of energy is strictly inverse, a steady state in the condensate regime is realised when the energy of the condensate is so large that the dissipation due to viscosity at large scales balances the energy injection due to the forcing. For split cascading systems, this is not necessarily true due to the presence of non-vanishing 3-D flow variations associated with the direct cascade. Therefore, in this case other processes exist that can redirect the energy back to the small scales where viscous dissipation is more efficient. Such mechanisms have been demonstrated for rotating turbulence, where a flux-loop mechanism has been identified (cf. Bartello, Métais & Lesieur 1994; Alexakis 2015; Seshasayanan & Alexakis 2018). Similar condensates have also been observed in 3-D fast rotating convection (Favier, Silvers & Proctor 2014; Guervilly, Hughes & Jones 2014; Rubio et al. 2014; Favier, Guervilly & Knobloch 2019).

Condensates in thin layers have been observed experimentally: the first study by Sommeria (1986) was followed by the important contributions of Paret & Tabeling (1997, 1998) and more recently by Shats, Xia & Punzmann (2005), Shats et al. (2007), Xia et al. (2008), Xia, Shats & Falkovich (2009), Byrne, Xia & Shats (2011), Xia et al. (2011). An up-to-date review of relevant experiments is presented in Xia & Francois (2017). These experiments operate primarily in the long-time limit in which the condensate is fully developed. This wealth of experimental studies of thin-layer condensates is in striking contrast with the existing numerical results which have focused exclusively on the transient growth of total kinetic energy due to the inverse cascade. In these numerical simulations, the condensate state reached after long time in the thin-layer case has not yet been examined due to the long computation time needed. In this study, we aim to fill this gap and investigate the behaviour of turbulent flow at the condensate stage for a thin layer forced at intermediate scales, using direct numerical simulations (DNS) and low-order modelling. The DNS provide a detailed picture of the behaviour of the full system, while the modelling shines light on the main physical processes involved in the problem.

The remainder of this article is structured as follows. In § 2, we present the set-up and define the quantities we will be measuring. In § 3, we present the results of a large number of direct numerical simulations (DNS) of thin-layer turbulence. Next, in § 4, we discuss the behaviour close to the two critical points and in § 5, we present spectra and spectral fluxes of energy. In § 6, we introduce a low-order model which captures many features of the DNS results. Finally, in § 7, we discuss our results and summarise.

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}$$

where $\boldsymbol{u}$ is the velocity field, $P$ is physical pressure divided by constant density and $\unicode[STIX]{x1D708}$ is (kinematic) viscosity. Energy is injected into the system by $\boldsymbol{f}$ , a stochastic force, that depends only on $x$ and $y$ and has only $x$ and $y$ components, i.e. that is a two-dimensional–two-component (2D2C) field. We make this assumption firstly to specifically force the inversely cascading components of the flow and secondly because it is widely used in previous studies such as Smith et al. (1996), Celani et al. (2010), Deusebio et al. (2014), Gallet & Doering (2015) and thus enables us to compare more easily with the literature. The force is divergence free, hence it can be written as $\boldsymbol{f}=(-\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D713},\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713},0)$ . The spectrum of $\boldsymbol{f}$ is concentrated in a ring of wavenumbers of radius $k_{f}\equiv 2\unicode[STIX]{x03C0}/\ell$ . It is delta correlated in time, which leads to a fixed mean injection of energy $\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{f}\rangle =\unicode[STIX]{x1D716}$ , where $\langle \cdot \rangle$ denotes an ensemble average over infinitely many realisations. We use random initial conditions whose small energy is spread out over a range of wavenumbers. In some cases, in order to compare with previous studies, we used a hyper-viscosity, which amounts to replacing $\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}$ by $-\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. 2011). 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. 2010)), 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. (2011).

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. (2009) 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. 2010) 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 2017). 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. 2010; Benavides & Alexakis 2017).

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 1973; Dallas, Fauve & Alexakis 2015; Cameron, Alexakis & Brachet 2017). 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 2012). 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 2017). 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 (2017). More generally, multistability is observed in many turbulent flows, see Weeks et al. (1997), Ravelet et al. (2004) 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 1985; Platt, Spiegel & Tresser 1993). 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 (2002). 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 (2017) 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. (2011) and figure 3 of Celani et al. (2010).

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 1843) (see Darrigol (2017) for a historical review). Eddy viscosity was quantified for the first time by Boussinesq (1877) and later widely popularised through the works of Taylor (Taylor 1915, 1922), see also Kraichnan (1976). It has been estimated in various limits both in 2-D and 3-D flows (Meshalkin 1962; Sivashinsky & Yakhot 1985; Bayly & Yakhot 1986; Yakhot & Sivashinsky 1987; Hefer & Yakhot 1989; Dubrulle & Frisch 1991; Gama, Vergassola & Frisch 1994; Cameron, Alexakis & Brachet 2016; Alexakis 2018).

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 1991) 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}$$

Note in particular that eddy viscosities do not dissipate energy, but merely redistributes it between different scales. Adding the three model equations leads to

(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.

7 Conclusions

We present the first detailed numerical study of the steady state of thin-layer turbulence as a function of the system parameters using an extensive set of high-resolution simulations.

It is shown that the split cascade observed at early times of the flow evolution (Celani et al. 2010; Benavides & Alexakis 2017; Musacchio & Boffetta 2017) leads to the formation of condensate states in the long-time limit. Three different states were found for large $Re$ . (i) For very thick layers the system saturates in a regular 3-D turbulence state with no inverse cascade and negligible dissipation at large scales. (ii) At intermediate-layer thickness, a flux-loop condensate is formed in which part of the energy transferred to the condensate by the 2-D motions is transferred back to the small scales by the 3-D motions. (iii) For very thin layers, the system becomes two-dimensional and forms a 2-D-turbulence condensate, where the inversely cascading energy is balanced by the dissipation due to viscosity at large scales. The transition from 3-D turbulence to the flux-loop condensate occurs at a critical height $H_{3D}$ $(Q_{3D})$ that is a decreasing (increasing) function of $Re$ , but saturates at a $Re$ independent value for large $Re$ . For values of $H$ slightly smaller than $H_{2D}$ the amplitude of the large-scale velocity $U_{ls}^{2}$ jumps discontinuously to a large value and increases linearly after that. Close to the threshold, a hysteresis diagram was constructed where the system saturates to a different attractor (3-D turbulence or flux-loop condensate) depending on the initial conditions. Whether this hysteresis behaviour persists at larger $Re$ and larger box sizes $1/K$ remains an open question. The flux-loop condensate transitions to a 2-D turbulence condensate at a critical height $H_{2D}$ that scales like $H_{2D}\propto \ell Re^{-3/4}$ unlike the early stages of the development where $H_{2D}\propto \ell Re^{-1/2}$ (Benavides & Alexakis 2017). For the 2-D turbulence condensate, the large-scale energy was found to be inversely proportional to $Re$ and independent from $H$ . The transition from a flux-loop condensate to a 2-D turbulence condensate showed strong spatio-temporal intermittency leading to a scaling of the average 3-D energy as the square of the deviation from onset $U_{3D}^{2}\propto (H-H_{2D})^{2}$ , similarly as in Benavides & Alexakis (2017).

A three-mode model has been proposed which reproduces the DNS scalings of the critical points $H_{2D}$ and $H_{3D}$ as well as the amplitude of the condensate in the 2-D turbulence regime. The model demonstrates the basic mechanisms involved: a 2-D flow that moves energy from the forcing scale to the condensate and a 3-D flow that takes away energy both form the large scales and the forcing scales. The model does not describe bi-stability or discontinuity close to $Q_{3D}$ . Nonetheless, it captures the occurrence of both transitions observed in the DNS and provides several correct quantitative predictions.

We stress once more that the present work is the first numerical study of thin-layer turbulent condensation. Previous studies of the thin-layer problem were restricted to the transient inverse cascade regime due the long computation time needed to reach the condensate state. Therefore, the present study is novel and provides an important first step towards a better understanding of thin-layer turbulent condensates (of which the Earth’s atmosphere and ocean may be viewed as examples, despite the idealised nature of our set-up), many open questions remain. The complexity of the physics involved close to criticality goes beyond the mean-field model and requires further targeted studies. We are convinced that both critical points deserve more detailed investigations by means of numerical simulations, experiments and modelling. Another important remaining open problem is the formation of an inverse cascade from a 3-D forcing. To study this, the amplitude of the 3-D components of the forcing should be varied compared to the 2-D components in a future study. Three-dimensional forcing will make a connection with more natural forcing mechanisms like convection that also display condensates (Favier et al. 2014; Guervilly et al. 2014; Rubio et al. 2014).

Concerning the realisability of the present numerical results in an experiment, it needs to be stressed that this study only considers the triply periodic domain for simplicity. When attempting to transfer the results to no-slip boundary conditions, a word of caution is therefore in order: viscous boundary layers may lead to large-scale drag, which is explicitly left out from the model set-up used here. Also, 3-D turbulence in boundary layers may infect the interior flow, thereby affecting even high wavenumbers and the two-dimensionalisation even in the bulk of the flow. However, the wealth of experimental observations of turbulent condensates in thin layers, as referenced in the introduction and summarised in Xia & Francois (2017), suggests that the condensation phenomenon at finite height is robust between different boundary conditions as well as between the different forcing methods used in experiment and numerical simulations. In particular, it would be very interesting to probe the discontinuity and associated phenomena reported here in an experiment. This has not been done before and experimental studies of thin-layer turbulent condensates have the advantage of allowing higher Reynolds numbers and much better time statistics.