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 Reference Pedlosky2013) 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 Reference Frisch1995). 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 Reference Boffetta and Ecke2012). 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 Reference Smith, Chasnov and Waleffe1996). In fact, it has been shown in (Celani, Musacchio & Vincenzi Reference Celani, Musacchio and Vincenzi2010; Benavides & Alexakis Reference Benavides and Alexakis2017; Musacchio & Boffetta Reference Musacchio and Boffetta2017) 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 (Reference Benavides and Alexakis2017), 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 (Reference Gallet and Doering2015). Similar transitions from a strictly forward cascade to an inverse cascade have been observed in other systems like rotating turbulence (Deusebio et al. Reference Deusebio, Boffetta, Lindborg and Musacchio2014), stratified turbulence (Sozza et al. Reference Sozza, Boffetta, Muratore-Ginanneschi and Musacchio2015), rotating and stratified flows (Marino, Pouquet & Rosenberg Reference Marino, Pouquet and Rosenberg2015), magneto-hydrodynamic systems (Alexakis Reference Alexakis2011; Seshasayanan, Benavides & Alexakis Reference Seshasayanan, Benavides and Alexakis2014; Seshasayanan & Alexakis Reference Seshasayanan and Alexakis2016) and helically constrained flows (Sahoo & Biferale Reference Sahoo and Biferale2015; Sahoo, Alexakis & Biferale Reference Sahoo, Alexakis and Biferale2017), to mention a few (see Alexakis & Biferale (Reference Alexakis and Biferale2018) 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 (Reference Kraichnan1967), first seen in direct numerical simulations by Hossain, Matthaeus & Montgomery (Reference Hossain, Matthaeus and Montgomery1983), further explored quantitatively by Smith & Yakhot (Reference Smith and Yakhot1993, Reference Smith and Yakhot1994) and, more recently, by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007), Bouchet & Simonnet (Reference Bouchet and Simonnet2009), Chan, Mitra & Brandenburg (Reference Chan, Mitra and Brandenburg2012), Frishman, Laurie & Falkovich (Reference Frishman, Laurie and Falkovich2017), Frishman & Herbert (Reference Frishman and Herbert2018). Spectral condensation has also been studied in other quasi-2-D systems such as quasi-geostrophic flows (see Vallis & Maltrud Reference Vallis and Maltrud1993; Kukharkin, Orszag & Yakhot Reference Kukharkin, Orszag and Yakhot1995; Kukharkin & Orszag Reference Kukharkin and Orszag1996; Venaille & Bouchet Reference Venaille and Bouchet2011). 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 Reference Bartello, Métais and Lesieur1994; Alexakis Reference Alexakis2015; Seshasayanan & Alexakis Reference Seshasayanan and Alexakis2018). Similar condensates have also been observed in 3-D fast rotating convection (Favier, Silvers & Proctor Reference Favier, Silvers and Proctor2014; Guervilly, Hughes & Jones Reference Guervilly, Hughes and Jones2014; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014; Favier, Guervilly & Knobloch Reference Favier, Guervilly and Knobloch2019).
Condensates in thin layers have been observed experimentally: the first study by Sommeria (Reference Sommeria1986) was followed by the important contributions of Paret & Tabeling (Reference Paret and Tabeling1997, Reference Paret and Tabeling1998) and more recently by Shats, Xia & Punzmann (Reference Shats, Xia and Punzmann2005), Shats et al. (Reference Shats, Xia, Punzmann and Falkovich2007), Xia et al. (Reference Xia, Punzmann, Falkovich and Shats2008), Xia, Shats & Falkovich (Reference Xia, Shats and Falkovich2009), Byrne, Xia & Shats (Reference Byrne, Xia and Shats2011), Xia et al. (Reference Xia, Byrne, Falkovich and Shats2011). An up-to-date review of relevant experiments is presented in Xia & Francois (Reference Xia and Francois2017). 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:
The system (2.1b ) is characterised by three non-dimensional parameters: the Reynolds number based on the energy injection rate $Re=(\unicode[STIX]{x1D716}\ell ^{4})^{1/3}/\unicode[STIX]{x1D708}$ , the ratio between forcing scale and domain height $Q=\ell /H$ and the ratio between forcing scale and the horizontal domain size $K=\ell /L$ . The ratio between $K$ and $Q$ gives the aspect ratio $A=K/Q=H/L$ of the domain. The Kolmogorov dissipation length is denoted as $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D708}^{3/4}/\unicode[STIX]{x1D716}^{1/4}=\ell Re^{3/4}$ .
The simulations performed for this work used an adapted version of the Geophysical High-Order Suite for Turbulence (GHOST) which uses pseudo-spectral methods including $2/3$ aliasing to solve for the flow in the triply periodic domain, (see Mininni et al. Reference Mininni, Rosenberg, Reddy and Pouquet2011). The resolution was varied from $256^{2}\times 16$ grid points to $2048^{2}\times 128$ grid points depending on the values of the parameters. To explore the space spanned by these three parameters, we have performed systematic numerical experiments: for a fixed value of $Re$ and $K=1/8$ , different simulations are performed with $Q$ varying from small to large values. The runs are continued until a steady state is reached where all quantities fluctuate around their mean value. This is repeated for eight different values of $Re$ from $Re=203$ (resolution $256^{2}\times 16$ ) to $4062$ (resolution $2048^{2}\times 128$ ) and for one value of hyper-viscosity ( $n=8$ , $\unicode[STIX]{x1D708}_{8}=10^{-38}$ as in (Celani et al. Reference Celani, Musacchio and Vincenzi2010)), as a consistency check, since many of the previous studies of thin-layer turbulence used hyper-viscosity. For $Re=305$ , we also perform a run with $K=1/16$ ( $L\rightarrow 2L$ ). The number of runs performed for each $Re$ are summarised in table 1.
To quantify the energy distribution among different scales it is convenient to work in Fourier space. The Fourier series expansion of the velocity reads
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
In addition, we monitor different components of domain-integrated energy, namely the total horizontal kinetic energy
(based on the (vertically averaged) 2D2C field only), the large-scale horizontal kinetic energy
where $k_{max}=\sqrt{2}2\unicode[STIX]{x03C0}/L$ , as well as the (vertically averaged) large-scale kinetic energy in the $z$ component
and the three-dimensional kinetic energy (3-D energy), defined as
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.
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 4 shows the slopes of the initial total energy increase $\unicode[STIX]{x1D716}_{inv}$ , measured as illustrated in figure 3(a) for all our numerical simulations. The slopes are non-dimensionalised by the energy input rate $\unicode[STIX]{x1D716}$ and plotted versus $Q$ for all different values of $Re$ including the hyper-viscous runs. The slope at this early stage measures the strength of inverse energy transfer. At small $Q$ (deep layers), the slope vanishes for all runs, showing that no inverse cascade is present. Moving to larger $Q$ , for every $Re$ , there is a critical value $Q_{3D}(Re)$ of $Q$ above which the slope becomes non-zero. This is the birth of the inverse cascade. Figure 5 shows estimates of $Q_{3D}$ as a function of $Re$ : the upper curve shows the smallest $Q$ for which an inverse cascade was observed for that given $Re$ while the lower curve shows the largest $Q$ for which no inverse cascade was observed. The critical value $Q_{3D}$ lies between these two curves. The point $Q_{3D}$ shifts to larger $Q$ as $Re$ is increased but eventually for the two largest $Re$ simulated, namely $Re=2031$ and $Re=4062$ , as well as the hyper-viscous run, $Q_{3D}$ saturates at $Q_{3D}\approx 2.5$ . (Previous findings (Celani et al. Reference Celani, Musacchio and Vincenzi2010) estimated this value to $Q_{3D}\approx 2$ , however in that work too limited a range of values of $Q$ was used to be able to precisely pinpoint $Q_{3D}$ . Another possible reason for the different result is the different value of $1/K$ associated with the different forcing wavenumber $k_{f}=16$ used). The saturation of $Q_{3D}\approx 2.5$ indicates that $Q_{3D}$ converges to this value at large $Re$ . For $Q>Q_{3D}$ , the slope begins increasing linearly $\unicode[STIX]{x1D716}_{inv}\propto Q-Q_{3D}$ . (We note that small slopes are hard to distinguish from zero slope since the difference only becomes apparent after a long simulation time.)
If $Q$ is increased further, a point $Q_{2D}$ is reached beyond which the slope becomes independent of $Q$ . Above this second critical point, the flow becomes exactly two-dimensional (Benavides & Alexakis Reference Benavides and Alexakis2017). The value of $Q_{2D}$ increases with $Re$ as $Q_{2D}\propto Re^{1/2}$ . This scaling is verified in this work as well and shown in figure 4(b). The two critical points $Q_{3D}$ and $Q_{2D}$ at this early stage of development of the inverse cascade have been studied in detail in the past (Celani et al. Reference Celani, Musacchio and Vincenzi2010; Benavides & Alexakis Reference Benavides and Alexakis2017).
Here we mostly focus on the second stage of evolution: the steady state and the properties of the condensate. Figure 6 shows the equilibrium value of $U_{ls}^{2}$ , as defined in equation (2.5), non-dimensonalised by the forcing energy scale $(\unicode[STIX]{x1D716}\ell )^{2/3}$ and multiplied by $K^{2}$ . In the left panel, it is plotted versus $Q$ (figure 6 a) and in the right panel it is rescaled by $1/Re$ and plotted versus $\unicode[STIX]{x1D702}/H=QRe^{-3/4}$ (figure 6 b). First consider figure 6(a). At small $Q$ , there is very little energy in the large scales. This corresponds to the values of $Q$ that displayed no inverse cascade at the initial stage. In the absence of an inverse cascade, the large scales only possess a small non-zero energy and are expected to be in a thermal equilibrium state (Kraichnan Reference Kraichnan1973; Dallas, Fauve & Alexakis Reference Dallas, Fauve and Alexakis2015; Cameron, Alexakis & Brachet Reference Cameron, Alexakis and Brachet2017). For $Q>Q_{3D}$ the energy in the large scale takes larger values. In all cases, the energy increases nearly linearly $U_{ls}^{2}\propto (Q-Q_{3D})$ for $Q_{2D}>Q>Q_{3D}$ . With the chosen coordinate, the close coincidence of the experiments with $K=1/8$ and $K=1/16$ at $Re=305$ indicates a scaling of $U_{ls}^{2}\propto L^{2}$ . If we zoom in on $Q_{3D}$ (see figure 7), we observe clear signs of small but discontinuous jumps of $U_{ls}^{2}$ at $Q_{3D}$ that are not visible in the zoomed out figure 6(a). These cases are examined in more detail in the next section.
The increase of the large-scale energy stops at the second critical point $Q_{2D}$ , where $U_{ls}^{2}$ becomes independent of $Q$ . It is noteworthy that the curves for various values of $Re$ all follow the same straight line between their respective $Q_{3D}$ and $Q_{2D}$ with only some deviations at low $Q$ . Furthermore, both $Q_{2D}$ and the plateau value of $U_{ls}^{2}$ depend on $Re$ . In figure 6(b), the same data are plotted, but with rescaled axes. The rescaling collapses the data well, with some deviations at small $Q$ related to the convergence of $Q_{3D}$ . This indicates that at large values of $Q$ , $U_{ls}^{2}$ scales like $U_{ls}^{2}\propto (\unicode[STIX]{x1D716}\ell )^{2/3}Re$ . This is precisely the scaling of the condensate of 2-D turbulence (Boffetta & Ecke Reference Boffetta and Ecke2012). The critical value where the transition to this maximum value of $U_{ls}^{2}$ occurs is $Q_{2D}Re^{-3/4}=\unicode[STIX]{x1D702}/H_{2D}\approx 0.09-0.1$ .
The scaling allowing us to collapse the data in figure 4 (transient stage) is different from that in figures 6, 8 and 9 (condensate state). This implies that $Q_{2D}\propto Re^{1/2}$ estimated during the early stage of the inverse cascade development is different from $Q_{2D}\propto Re^{3/4}$ estimated at steady state where a condensate is fully developed. The reason for this difference is that the transition to exactly 2-D motion occurs when the maximum shear in the flow (which produces 3-D motion by shear instabilities) is balanced by small-scale viscous dissipation. In the presence of the inverse cascade, an $E(k)\propto \unicode[STIX]{x1D716}^{2/3}k^{-5/3}$ spectrum is formed at $k>k_{f}$ , such that the peak of the enstrophy spectrum $k^{2}E(k)$ is at the forcing scale. Thus the balance between 2-D shear and 3-D damping is
implying $H_{2D}\sim \ell Re^{-1/2}$ (Benavides & Alexakis Reference Benavides and Alexakis2017). In the condensate, however, most of the energy and enstrophy are located in the largest scales and are such that energy injection $\unicode[STIX]{x1D716}$ is balanced by large-scale dissipation $\propto \unicode[STIX]{x1D708}U_{ls}^{2}/L^{2}$ . The large-scale shear is thus $U_{ls}/L\propto (\unicode[STIX]{x1D716}/\unicode[STIX]{x1D708})^{1/2}$ which is balanced by the damping rate of 3-D perturbations at onset,
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 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.
The following remarks are in order: although for each $Q$ we ran the simulations until saturation was achieved, since we are dealing with a noisy system, rare transitions can exist between the two branches of the hysteresis loop. To test this, we picked the point $Q\approx 1.97$ on both lower and upper branches and ran them for a long time (thousands of eddy turnover times $\unicode[STIX]{x1D70F}=(L^{3}/\unicode[STIX]{x1D716})^{1/2}$ ). In neither case did we see a transition between the two branches, indicating that such transitions are rare (if not absent) in the middle of the hysteresis loop. Near the edges of the hysteresis loop at $Q\approx 2.05$ and $Q\approx 1.9$ , the dependence on simulation time is likely to be stronger, but this has not been investigated.
Furthermore, we note that the bifurcation diagram of figure 10 corresponds to a relatively low Reynolds number $Re=406$ . Whether this subcritical behaviour persists at larger $Re$ and/or larger box sizes (smaller $K$ ) is still an open question. Figure 7 suggests that a discontinuity continues to exist at $Q=Q_{3D}$ up to high Reynolds numbers ( $Re=2031$ shown there). In addition, we found more points at higher $Re$ that showed a dependence on initial conditions but without having enough values of $Q$ to create a hysteresis diagram. These findings suggest that subcritical behaviour and hysteresis might survive even at high $Re$ . However, due to the high computational cost at higher resolution and the long duration of the runs required to verify that the system stays in a particular state, we could not investigate this possibility in detail. Further simulations at larger $Re$ and possibly smaller $K$ (larger boxes) are required to resolve this issue. Similar hysteretic behaviour has recently been reported in rotating turbulence, see Yokoyama & Takaoka (Reference Yokoyama and Takaoka2017). More generally, multistability is observed in many turbulent flows, see Weeks et al. (Reference Weeks, Tian, Urbach, Ide, Swinney and Ghil1997), Ravelet et al. (Reference Ravelet, Marié, Chiffaudel and Daviaud2004) as examples.
4.2 Close to $Q_{2D}$ : intermittent bursts
Next, we discuss the behaviour of the flow close to the second critical point $Q_{2D}$ . A typical time series of 3-D energy for $Q\lesssim Q_{2D}$ is shown in figure 12(a). One observes bursty behaviour and variations over many orders of magnitude, indicating on–off intermittency (Fujisaka & Yamada Reference Fujisaka and Yamada1985; Platt, Spiegel & Tresser Reference Platt, Spiegel and Tresser1993). On–off intermittency refers to the situation where a marginally stable attractor loses or gains stability due to noise fluctuations. When instability is present, a temporary burst is produced before the system returns to the attractor. On–off intermittency predicts that the unstable mode $X$ follows a power-law distribution $P(X)\propto X^{\unicode[STIX]{x1D6FF}-1}$ for $X\ll 1$ where $\unicode[STIX]{x1D6FF}$ measures the deviation from onset (here $\unicode[STIX]{x1D6FF}\propto (Q_{2D}-Q)/Q_{2D}$ ) and all moments scale linearly with the deviation $\langle X^{n}\rangle \propto \unicode[STIX]{x1D6FF}$ .
In our system, the 2-D flow forms the marginal stable attractor that loses stability to 3-D perturbations depending on the exact realisation of the 2-D turbulent flow. To formulate this, we decompose the velocity field into its 2-D and 3-D parts, $\boldsymbol{u}=\boldsymbol{u}_{2D}+\boldsymbol{u}_{3D}$ , where the 2-D part is defined as the Fourier sum of $\boldsymbol{u}$ restricted to modes with $k_{z}=0$ . Filtering the 3-D component of equation (2.1b ), dotting with $\boldsymbol{u}_{3D}$ and integrating over the domain gives
where $\langle \cdot \rangle$ denotes integration over the domain. The chaotic 2-D motions then act as multiplicative noise while the viscous terms provide a mean decay rate. An important physical mechanism for creating 3-D disturbances is 3-D elliptic instability of the 2-D counter-rotating vortex pair forming the condensate, as described in Le Déz & Laporte (Reference Le Déz and Laporte2002). This may explain the presence of the critical value $Q_{2D}$ itself: instability requires small vertical wavenumbers, but the minimum wavenumber increases with decreasing $H$ and $Q_{2D}$ corresponds to the point where 3-D perturbations begin to decay.
Figure 12 shows that temporal intermittency is present in the thin-layer system. Panel (a) shows a typical time series of 3-D energy at $Q\lesssim Q_{2D}$ which fluctuates over six orders of magnitude. In particular, as mentioned before, there are burst-like excursions in 3-D energy. In figure 12(b), PDFs constructed from this time series and similar ones for different values of $Q$ are shown along with dotted lines indicating power laws with exponents $-1$ , $-0.8$ and $-0.3$ . The PDFs are very close to a power law for a significant range of $U_{3D}^{2}$ and the exponent converges to minus one as the transition is approached, in agreement with on–off intermittency predictions. However, the scaling of 3-D energy with deviation from onset shown in figure 8(b) does not follow the linear prediction of on–off intermittency, but rather $\langle U_{3D}^{2}\rangle \propto (Q_{2D}-Q)^{2}$ . For $U_{z}^{2}$ , figure 9 seems to suggest yet a different scaling, namely $\langle U_{z}^{2}\rangle \propto (Q_{2D}-Q)^{3}$ . A similar behaviour was also found in Benavides & Alexakis (Reference Benavides and Alexakis2017) and was attributed to the spatio-temporal character of the intermittency that not only leads to 3-D motions appearing more rarely in time as criticality is approached but also to them occupying a smaller fraction of the available volume. This appears also to be the case in our results, as demonstrated in figures 13 and 14, where $\boldsymbol{u}_{3D}^{2}$ and the vertical velocity $u_{z}$ are plotted for three different values of $Q$ . As $Q$ approaches the critical value $Q_{2D}$ , the structures become smaller for $\boldsymbol{u}_{3D}^{2}$ and $u_{z}$ with the difference that $\boldsymbol{u}_{3D}^{2}$ shows spot-like structures in figure 13(c) which are absent for $u_{z}$ . This difference may be related to the two different scalings observed for $U_{3D}^{2}$ and $U_{z}^{2}$ with $Q_{c}-Q$ small: if the volume fraction of vertical motion depends on $Q_{c}-Q$ to a different power than that of vertical variations, two different behaviours of $U_{z}^{2}$ and $U_{3D}^{2}$ would follow. A more detailed quantitative investigation of the scaling of volume fraction will be needed to clarify this.
In summary, we have found non-trivial behaviour close to both transitions: we have observed hysteresis near $Q_{3D}$ and spatio-temporal intermittency close to $Q_{2D}$ where the temporal behaviour seems to be described by on–off intermittency. Taking into account these effects will be crucial for understanding the exact nature of the observed transitions.
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.
Moreover, the total 1-D energy spectrum may advantageously be split up into three components: the energy spectrum of the (vertically averaged) 2D2C field
the energy spectrum of the (vertically averaged) vertical velocity
and the energy spectrum of the 3-D flow defined as
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
where the low-pass filtered velocity field is
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
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
It expresses the flux due to all interactions other than the ones in (5.7).
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.
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.
As illustrated in figure 19, our model comprises a 2-D mode $U_{2D}$ at the scale $L$ of the domain, a mode $U_{f}$ at the forcing scale $\ell$ and a 3-D mode $U_{3D}$ at the scale of the layer height $H$ , whose interactions are spectrally non-local, thus taking into account a major result from § 5. The model describes the system at steady state where these scales are well separated, but is not expected to capture the transient phase where all intermediate scales between $L$ and $\ell$ participate due to the inverse cascade. As before, let $Q=\ell /H$ , $K=\ell /L$ and $Re=(\unicode[STIX]{x1D716}\ell ^{4})^{1/3}/\unicode[STIX]{x1D708}$ . Interactions between modes are modelled using eddy viscosity, which amounts to modifying the molecular viscosity $\unicode[STIX]{x1D708}$ by terms involving the small-scale velocities, modelling the effect of small-scale motions on large-scale motions as diffusive. The conceptual foundations of eddy viscosity were laid by de Saint Venant in his effective viscosity, (de Saint-Venant Reference de Saint-Venant1843) (see Darrigol (Reference Darrigol2017) for a historical review). Eddy viscosity was quantified for the first time by Boussinesq (Reference Boussinesq1877) and later widely popularised through the works of Taylor (Taylor Reference Taylor1915, Reference Taylor1922), see also Kraichnan (Reference Kraichnan1976). It has been estimated in various limits both in 2-D and 3-D flows (Meshalkin Reference Meshalkin1962; Sivashinsky & Yakhot Reference Sivashinsky and Yakhot1985; Bayly & Yakhot Reference Bayly and Yakhot1986; Yakhot & Sivashinsky Reference Yakhot and Sivashinsky1987; Hefer & Yakhot Reference Hefer and Yakhot1989; Dubrulle & Frisch Reference Dubrulle and Frisch1991; Gama, Vergassola & Frisch Reference Gama, Vergassola and Frisch1994; Cameron, Alexakis & Brachet Reference Cameron, Alexakis and Brachet2016; Alexakis Reference Alexakis2018).
There are two notable cases where the dependence of eddy viscosity $\unicode[STIX]{x1D708}_{E}$ on the flow amplitude $U_{s}$ and length scale $l_{s}$ is known. For $Re\rightarrow \infty$ , one expects that $\unicode[STIX]{x1D708}_{E}$ becomes independent of $\unicode[STIX]{x1D708}$ and the only dimensionally consistent possibility for $\unicode[STIX]{x1D708}_{E}$ is given by
where $c_{1}$ is a non-dimensional number. In the low- $Re$ limit, on the other hand, an exact asymptotic expansion can be carried out (see Dubrulle & Frisch Reference Dubrulle and Frisch1991) which reveals that
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:
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
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:
(i) laminar state: $U_{2D}=U_{3D}=0$ (all energy in forcing scale),
(ii) three-dimensional turbulence state: $U_{2D}=0$ and $U_{3D}\neq 0$ ,
(iii) two-dimensional condensate state: $U_{2D}\neq 0$ and $U_{3D}=0$ and
(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.
From this simple model, three major predictions may be derived:
(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.
(ii) Secondly, the critical point $Q_{2D}$ is predicted to obey
(6.6) $$\begin{eqnarray}Q_{2D}\propto Re^{3/4}.\end{eqnarray}$$(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. Reference Celani, Musacchio and Vincenzi2010; Benavides & Alexakis Reference Benavides and Alexakis2017; Musacchio & Boffetta Reference Musacchio and Boffetta2017) 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 Reference Benavides and Alexakis2017). 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 (Reference Benavides and Alexakis2017).
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. Reference Favier, Silvers and Proctor2014; Guervilly et al. Reference Guervilly, Hughes and Jones2014; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014).
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 (Reference Xia and Francois2017), 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.
Acknowledgements
We thank three referees for their helpful and valuable suggestions which have helped to improve this paper substantially. This work was granted access to the HPC resources of MesoPSL financed by the Région Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche and the HPC resources of GENCI-TGCC-CURIE & CINES Occigen (Project no. A0010506421, A0030506421 and A0050506421) where the present numerical simulations have been performed. This work has also been supported by the Agence nationale de la recherche (ANR DYSTURB project no. ANR-17-CE30-0004). A.v.K. was supported by Deutscher Akademischer Austauschdienst and Studienstiftung des deutschen Volkes.
Appendix A. Derivation of mean-field model predictions
A.1 Scalings of critical points and condensate amplitude
In the low viscosity limit, the eddy viscosities given in equation (6.5), take the form (6.1) and the resulting system of equations reads
Using this result and considering equation (A 1a ), we can find that the 3-D turbulence fixed point becomes unstable to 2-D perturbations at
and thus we obtain that
Hence, in the low viscosity limit of our three-scale model, there remains only one bifurcation, namely that at $Q_{3D}$ between two-dimensional turbulence and the split cascade state. The second critical point $Q_{2D}$ vanishes to infinity as $Q_{2D}\propto Re^{3/4}$ in this limit. Figure 20 demonstrates close to $Q_{3D}$ that the full model converges to the solution obtained from the asymptotic form of the equations (A 1) as $Re$ increases. This is consistent with the convergence observed in the DNS in figure 6(b).
At finite viscosity, one has to solve the full equations, (6.5) which is difficult analytically for the 2-D condensate state. In order to facilitate analytical progress in deriving predictions from the model, one may formally take the high viscosity limit in which the different eddy viscosities take the form of equation (6.2). The model equations then become
When the latter condition is satisfied and $H$ is sufficiently small ( $Q$ sufficiently large), the system is attracted to the 2-D condensate state, given by
Note that $U_{2D}^{2}$ is inversely proportional to the viscosity and proportional to $L^{2}$ in agreement with the scaling of the data in figure 6(b). The 2-D condensate state ceases to be an attractor of the system when $H$ is sufficiently large such that $U_{3D}$ becomes unstable. This occurs when
Hence, we conclude that
Thus, for moderate values of $Re$ , there is an approximate scaling $Q_{2D}\propto Re^{3/4}\propto \unicode[STIX]{x1D702}$ , the dissipation length (note that $Re^{3}>1/\unicode[STIX]{x1D6FC}$ due to (A 6)), in agreement with the results obtained in § 3, where we showed that the $U_{2D}^{2}$ data points collapse under rescaling such that $QRe^{3/4}=\unicode[STIX]{x1D702}/H$ is on the abscissa and $U_{2D}^{2}K^{2}/[Re(\unicode[STIX]{x1D716}\ell )^{2/3}]$ on the coordinate. Results from the full model equations (6.3) and (6.5) are shown in figure 21 where the same scaling is applied. The corresponding plots for equations (A 5) are very similar. Furthermore, an asymptotic analysis close to $Q_{2D}$ described in § A.2 of this appendix reveals that the scaling for $U_{3D}^{2}\propto (Q_{2D}-Q)^{2}$ which is the same as in the DNS results shown in figure 8(b) although no intermittency is present. The other critical point $Q_{3D}$ , where the 3-D turbulence solution changes stability, can be evaluated numerically and is found to increase with $Re$ indefinitely. This, however, is an artefact of the high viscosity asymptotic form of the eddy viscosities used in this subsection.
A.2 Behaviour of $U_{3D}$ near $Q_{2D}$
Here, we derive the behaviour close to $H_{2D}$ in the three-scale model. First consider $H=H_{2D}(1+\unicode[STIX]{x1D6FF})$ and let $\boldsymbol{x}=(x,y,z)^{T}=(U_{2D}^{2},U_{f}^{2},U_{3D}^{2})^{T}$ , $\tilde{\boldsymbol{x}}=(\tilde{x},{\tilde{y}},\tilde{z})^{T}=\boldsymbol{x}-(x_{2D},y_{2D},0)^{T}$ , where $x_{2D}$ and $y_{2D}$ are the values of $U_{2D}^{2}$ and $U_{f}^{2}$ respectively at $H=H_{2D}$ . Then equations (A 5) can be rewritten exactly in the form
where $C=-\unicode[STIX]{x1D708}/H^{2}+\unicode[STIX]{x1D6FD}H^{2}x_{2}/\unicode[STIX]{x1D708}L^{2}+\unicode[STIX]{x1D6FE}H^{2}y_{2}/\unicode[STIX]{x1D708}\ell ^{2}$ and the specific coefficients of the quadratic term are irrelevant here. By definition of $H_{2D}$ , $C(\unicode[STIX]{x1D6FF}=0)=-\unicode[STIX]{x1D708}/H_{2D}^{2}+\unicode[STIX]{x1D6FD}H_{2D}^{2}x_{2}/\unicode[STIX]{x1D708}L^{2}+\unicode[STIX]{x1D6FE}H_{2D}^{2}y_{2}/\unicode[STIX]{x1D708}\ell ^{2}=0$ . Hence, for small $\unicode[STIX]{x1D6FF}$ , $C\propto \unicode[STIX]{x1D6FF}$ . Specifically,
Hence, considering the $\tilde{z}$ component and balancing the linear term with the $\tilde{z}^{3/2}$ term, we deduce that
This means that $U_{3D}^{2}\propto \unicode[STIX]{x1D6FF}^{2}$ , which is precisely the scaling observed in figure 8(b). It is important to note however that the asymptotic result (A 12) is only valid for very small $\unicode[STIX]{x1D6FF}$ and cannot be extended to $\unicode[STIX]{x1D6FF}\sim O(1)$ where the quadratic terms are dominant.
A.3 On–off intermittency in the three-scale model
When a fluctuating energy injection rate is taken into account in the model by replacing $\unicode[STIX]{x1D716}\rightarrow \unicode[STIX]{x1D716}+\unicode[STIX]{x1D70E}\unicode[STIX]{x1D701}$ , where $\unicode[STIX]{x1D701}\sim {\mathcal{N}}(0,1)$ is Gaussian white noise, on–off intermittency in $U_{3D}^{2}$ can be observed in the three-scale model. This is illustrated in figure 22 in terms of the time series of $U_{3D}^{2}$ and the corresponding PDF, which approaches a power law with exponent $-1$ as $Q\rightarrow Q_{2D}$ .