1 Introduction
Turbulence in geophysical and astrophysical systems is a problem of major importance. Three-dimensional (3D) flows favour a forward cascade, i.e., energy flows from large to small scales, as described by the well-known Kolmogorov theory. In contrast, two-dimensional (2D) flows exhibit an inverse energy transfer, from small to large scales, that manifests itself in the appearance of large-scale structures in the flow (Kraichnan Reference Kraichnan1967; Boffetta & Ecke Reference Boffetta and Ecke2012). However, many systems arising in geophysical and astrophysical fluid dynamics fail to be fully 3D because of the presence of strong restraints, although they are far from being 2D either. These restraints may arise from geometrical confinement (Smith, Chasnov & Waleffe Reference Smith, Chasnov and Waleffe1996; Celani, Musacchio & Vincenzi Reference Celani, Musacchio and Vincenzi2010; Benavides & Alexakis Reference Benavides and Alexakis2017; Xia & Francois Reference Xia and Francois2017), rapid rotation (Smith & Waleffe Reference Smith and Waleffe1999; Pouquet & Marino Reference Pouquet and Marino2013; Campagne et al. Reference Campagne, Gallet, Moisy and Cortet2014), strong stratification (Bartello Reference Bartello1995; Smith & Waleffe Reference Smith and Waleffe2002; Oks et al. Reference Oks, Mininni, Marino and Pouquet2017) or the presence of strong magnetic fields (Alexakis Reference Alexakis2011; Favier et al. Reference Favier, Godeferd, Cambon, Delache and Bos2011). The detailed nature of the energy cascade in constrained 3D flows and its applications to geophysical fluid dynamics remain an open problem (Alexakis & Biferale Reference Alexakis and Biferale2018).
Recent numerical simulations of rapidly rotating Rayleigh–Bénard convection, hereafter RRRBC, have shown that the turbulent state is susceptible to the evolution of large-scale vortex (LSV) structures despite the presence of 3D fluctuations on all scales (Favier, Silvers & Proctor Reference Favier, Silvers and Proctor2014; Guervilly, Hughes & Jones Reference Guervilly, Hughes and Jones2014), in agreement with the predictions of an asymptotic description of the system valid in the limit of vanishingly small Rossby numbers (Julien et al. Reference Julien, Rubio, Grooms and Knobloch2012; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014). Together these studies reveal the presence of an efficient energy transfer mechanism that extracts energy from small-scale 3D fluctuations and deposits it in a box-scale barotropic (i.e. 2D) mode, bypassing the inverse energy cascade familiar from 2D dynamics. This process operates in parallel with energy transfer to large-scale 2D modes by barotropic–barotropic interactions (Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014) but dominates all aspects of the problem. Julien et al. (Reference Julien, Rubio, Grooms and Knobloch2012) and Julien, Knobloch & Plumley (Reference Julien, Knobloch and Plumley2018) have suggested that the presence of the LSV introduces essential correlations among the phases of the small-scale 3D fluctuations that facilitate direct energy extraction from these scales by the large-scale mode, leading to a runaway that is only arrested by additional processes omitted from the simplest problem formulation. Such a runaway is characteristic of subcritical dynamics, where it is triggered by finite-amplitude perturbations. The present paper is therefore devoted to a search for such subcritical behaviour in RRRBC. We do not address the question of whether geostrophic turbulence is itself linearly unstable to the generation of an LSV structure, although our simulations suggest that at large enough rotation rates and large enough Rayleigh numbers this is in fact the case. We emphasize that phase correlations are missed in studies that focus on energy spectra alone, and that such correlations are inevitably absent from flows driven by a prescribed small-scale force such as those studied by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007) and Bouchet & Simonnet (Reference Bouchet and Simonnet2009). In these systems the only possible correlations are between the applied force and the resulting velocity field. In contrast, in the present 3D system the forcing of the large-scale 2D flow can itself be dynamically affected through the action of the LSV on the small convective scales. The LSV observed in RRRBC may thus be a consequence of the proximity of the flow to 2D turbulence or due to the ability of the LSV to shape the correlations among the small-scale fluctuations that appear to drive its formation.
In the present work we provide the first evidence for subcritical dynamics in turbulent RRRBC by demonstrating the coexistence of two numerically stable turbulent states at identical parameter values, one with an LSV structure and one without. Such bistability in turbulent flows is rare, although it has also been found in rapidly rotating turbulence (Alexakis Reference Alexakis2015; Yokoyama & Takaoka Reference Yokoyama and Takaoka2017), thin-layer turbulence (van Kan & Alexakis Reference van Kan and Alexakis2019), Couette flows (Mujica & Lathrop Reference Mujica and Lathrop2006; Zimmerman, Triana & Lathrop Reference Zimmerman, Triana and Lathrop2011; Huisman et al. Reference Huisman, van der Veen, Sun and Lohse2014; Xia et al. Reference Xia, Shi, Cai, Wan and Chen2018) and von Kármán flows (Ravelet et al. Reference Ravelet, Marié, Chiffaudel and Daviaud2004).
2 Mathematical formulation
2.1 Model and governing equations
We consider the evolution of a layer of incompressible fluid, bounded above and below by two impenetrable, fixed-temperature, stress-free horizontal walls, a distance
$h$
apart. The layer rotates about the
$z$
-axis, pointing vertically upwards, with a constant angular velocity
$\unicode[STIX]{x1D734}=\unicode[STIX]{x1D6FA}\boldsymbol{e}_{z}$
, while gravity points downwards:
$\boldsymbol{g}=-g\boldsymbol{e}_{z}$
. The kinematic viscosity
$\unicode[STIX]{x1D708}$
and thermal diffusivity
$\unicode[STIX]{x1D705}$
are assumed to be constant.
In the Boussinesq approximation, using the thermal diffusion time
$h^{2}/\unicode[STIX]{x1D705}$
as a unit of time and the depth
$h$
of the layer as a unit of length, the dimensionless equations are
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{Pr}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\right)=-\unicode[STIX]{x1D735}p-\sqrt{Ta}\,\boldsymbol{e}_{z}\times \boldsymbol{u}+Ra\unicode[STIX]{x1D703}\,\boldsymbol{e}_{z}+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}=w+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{u}\equiv (u,v,w)$
is the velocity,
$p$
is the pressure and
$\unicode[STIX]{x1D703}$
is the temperature fluctuation with respect to a linearly decreasing background temperature. The parameters are the Rayleigh number
$Ra=\unicode[STIX]{x1D6FC}g\unicode[STIX]{x0394}Th^{3}/(\unicode[STIX]{x1D708}\unicode[STIX]{x1D705})$
, the Taylor number
$Ta=4\unicode[STIX]{x1D6FA}^{2}h^{4}/\unicode[STIX]{x1D708}^{2}$
and the Prandtl number
$Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$
. These dimensionless quantities involve
$\unicode[STIX]{x1D6FC}$
, the coefficient of thermal expansion, and
$\unicode[STIX]{x0394}T$
, the imposed temperature difference between the two horizontal plates. For simplicity, and following earlier studies of the formation of large-scales structures in this system, we take
$Pr=1$
. In the two horizontal directions, all variables are assumed to be periodic with the same spatial period
$\unicode[STIX]{x1D706}$
in both
$x$
and
$y$
directions. The boundary conditions at the upper (
$z=1$
) and lower (
$z=0$
) plates are
$\unicode[STIX]{x2202}_{z}u=\unicode[STIX]{x2202}_{z}v=w=\unicode[STIX]{x1D703}=0$
.
We solve (2.1)–(2.3) using the same mixed Fourier fourth-order finite-difference scheme as used in Favier et al. (Reference Favier, Silvers and Proctor2014). We confirm the robustness of the present results by additionally running equivalent simulations using the fully spectral approach of Guervilly et al. (Reference Guervilly, Hughes and Jones2014) and the open-source spectral-element code Nek5000 developed by Fischer, Lottes & Kerkemeier (Reference Fischer, Lottes and Kerkemeier2008) at the Argonne National Laboratory.
2.2 Finite-amplitude initial conditions
In order to explore possible subcritical behaviour in turbulent RRRBC, we consider a particular set of initial conditions. Anticipating that the non-local upscale energy transfer eventually saturates by creating a barotropic vortex dipole at the box scale (Julien et al. Reference Julien, Rubio, Grooms and Knobloch2012; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014), we consider a depth-invariant initial condition, given by
$$\begin{eqnarray}\boldsymbol{u}(t=0)=(A\sin (2\unicode[STIX]{x03C0}y/\unicode[STIX]{x1D706}),-A\sin (2\unicode[STIX]{x03C0}x/\unicode[STIX]{x1D706}),0)\quad \text{and}\quad \unicode[STIX]{x1D703}(t=0)=0,\end{eqnarray}$$
and parametrized by the amplitude
$A$
. This initial condition corresponds to a symmetric vortex dipole at the box scale with a cyclone located at the box centre
$x=\unicode[STIX]{x1D706}/2$
and
$y=\unicode[STIX]{x1D706}/2$
and an anticyclone located at the corners of the periodic domain. In addition, we add random perturbations of small amplitude
$\pm 2.5\times 10^{-2}$
to the temperature field in order to initiate the Rayleigh–Bénard instability. The total kinetic energy density of this initial flow is
$$\begin{eqnarray}K_{0}=\frac{1}{V}\int _{V}\frac{1}{2}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}\,\text{d}V=\frac{1}{2}A^{2},\end{eqnarray}$$
where
$V$
is the total volume. Note that
$K_{0}$
is independent of the aspect ratio of the box. The purely viscous decay of this initial condition is given by
$$\begin{eqnarray}K(t)=K_{0}\exp (-8\unicode[STIX]{x03C0}^{2}Pr~t/\unicode[STIX]{x1D706}^{2}).\end{eqnarray}$$
2.3 Flow decomposition
We define the depth-averaged 2D horizontal flow, subsequently called the 2D mode, as
$$\begin{eqnarray}\langle u\rangle _{z}(x,y)=\int _{0}^{1}u(x,y,z)\,\text{d}z\quad \text{and}\quad \langle v\rangle _{z}(x,y)=\int _{0}^{1}v(x,y,z)\,\text{d}z,\end{eqnarray}$$
where
$u$
and
$v$
are the velocity components in the
$x$
and
$y$
directions, respectively, and fast 3D fluctuations, subsequently called the 3D mode, as
$$\begin{eqnarray}\displaystyle & \displaystyle u^{\prime }(x,y,z)=u(x,y,z)-\langle u\rangle _{z}(x,y), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle v^{\prime }(x,y,z)=v(x,y,z)-\langle v\rangle _{z}(x,y). & \displaystyle\end{eqnarray}$$
We can then define the (purely horizontal) kinetic energy density associated with the slow 2D mode as
$$\begin{eqnarray}K_{2D}=\frac{1}{2\unicode[STIX]{x1D706}^{2}}\iint \left(\langle u\rangle _{z}^{2}+\langle v\rangle _{z}^{2}\right)\,\text{d}x\,\text{d}y\end{eqnarray}$$
and the kinetic energy density associated with the fast 3D mode as
$$\begin{eqnarray}K_{3D}=\frac{1}{2\unicode[STIX]{x1D706}^{2}}\iiint \left(u^{\prime 2}+v^{\prime 2}+w^{2}\right)\,\text{d}x\,\text{d}y\,\text{d}z.\end{eqnarray}$$
3 Results
In this paper, we fix
$Ta=10^{8}$
which is sufficient to sustain a spontaneous non-local inverse cascade for
$5\times 10^{6}\lesssim Ra\lesssim 2\times 10^{7}$
(Favier et al.
Reference Favier, Silvers and Proctor2014; Guervilly et al.
Reference Guervilly, Hughes and Jones2014). For
$Ra<5\times 10^{6}$
the flow is not turbulent enough, while for
$Ra>2\times 10^{7}$
the flow is insufficiently constrained by rotation and so not anisotropic enough, and only the traditional forward energy cascade is observed. In between these two limits, the flow is both turbulent and dynamically constrained by rotation, leading to the spontaneous emergence of a LSV from purely 3D perturbations driven by the Rayleigh–Bénard instability. In order to better understand the nature of the transition from a state with LSV to a state without, we consider here the particular case
$Ra=3\times 10^{7}$
, for which no systematic LSV were observed starting from random infinitesimal temperature perturbations (Favier et al.
Reference Favier, Silvers and Proctor2014; Guervilly et al.
Reference Guervilly, Hughes and Jones2014), and focus on the behaviour with aspect ratio
$\unicode[STIX]{x1D706}=4$
. This aspect ratio is chosen to ensure a clear scale separation between the convective eddies and the box size (see § 4 for a discussion of the effects of varying
$\unicode[STIX]{x1D706}$
). For this set of control parameters, the flow is dominated by 3D small-scale turbulent fluctuations (see figure 2
a,c), characterized by the Reynolds number
$Re=\sqrt{\langle w^{2}\rangle }\approx 535$
and the micro-Rossby number
$Ro_{\unicode[STIX]{x1D714}}=\sqrt{\langle \unicode[STIX]{x1D714}_{z}^{2}\rangle }/\sqrt{Ta}\approx 1.4$
. Here
$\unicode[STIX]{x1D714}_{z}$
is the vertical component of the vorticity and the brackets denote spatial and temporal averaging. Using this simulation as a reference, we ran a number of additional simulations, continuously varying the amplitude
$A$
of the initial vortex dipole defined in (2.4).
Figure 1. (a) Time evolution of the kinetic energy density
$K_{2D}$
for different initial vortex amplitudes
$A$
. The grey area corresponds to the transition where
$K_{2D}\approx K_{3D}$
. (b) Time evolution of the ratio between
$K_{2D}$
and the total kinetic energy
$K_{2D}+K_{3D}$
for different amplitudes
$A$
.
Figure 2. Visualizations of the quasi-steady states at
$t=0.3$
. (a,b) Vertical vorticity component. (c,d) Velocity field streamlines coloured with the velocity amplitude. (a,c) No initial vortex dipole
$A=0$
. (b,d) Initial vortex dipole amplitude
$A=1600$
.
Figure 1(a) shows the temporal evolution of the kinetic energy density
$K_{2D}$
of the 2D flow as defined by (2.10). We observe a clear transition as the initial amplitude
$A$
of the vortex dipole increases. For small amplitudes, typically
$A\lesssim 800$
, the kinetic energy of the 2D flow decreases until it reaches the equilibrium value corresponding to the reference case
$A=0$
(i.e., no initial vortex). Note that this decay closely follows the purely viscous decay of the initial condition as given by (2.6), shown as dashed lines in figure 1(a), indicating that there is no significant energy transfer from the 3D fluctuations to the 2D flow. For larger amplitudes however, typically
$A\gtrsim 800$
, we observe an initial decay of
$K_{2D}$
followed by an approximately linear increase until the energy eventually saturates at very long times. Note that close to the transition threshold, see case
$A=875$
for example, it is not yet clear whether the vortex will grow or decay. In view of figure 1(a), which shows that the large-scale vortex has not yet reached saturation, additional and longer simulations are required to determine any residual dependence of the saturated state on
$A$
, although we expect that all solutions involving a growing LSV will eventually saturate at the same amplitude.
These results clearly point towards subcritical behaviour in the transition between a weak subdominant 2D flow and a strong non-local inverse energy transfer efficiently feeding energy into the largest available spatial scale of the domain. Note that while the two states correspond to the same control parameters, the LSV state has a total kinetic energy density approximately eight times that of the 3D state. The ratio between
$K_{2D}$
and the total kinetic energy
$K_{2D}+K_{3D}$
is shown in figure 1(b). It appears that the finite-amplitude transition occurs once that
$K_{2D}(t=0)\approx K_{3D}$
, although it depends very subtly on the initial condition and can occur even after a long transient, as in the case
$A=875$
. A more accurate estimate of the critical amplitude
$A$
is beyond the scope of this paper, however, since it would require running many realizations only changing the initial temperature perturbation (see § 4 below for a brief discussion concerning the non-deterministic nature of this transition). Visualizations of the saturated states for the same control parameters, but two different initial conditions, are shown in figure 2. Without the initial LSV, or when the amplitude
$A$
is too small, the equilibrium state is dominated by a 3D flow at small scales while the 2D flow remains marginal. These small-scale fluctuations are fully 3D, as expected from the moderate value of the Rossby number,
$Ro=0.55$
(see table 1). In contrast, above the critical value of
$A$
, the LSV is continuously amplified while remaining at the box scale.
Table 1. Summary of the parameters considered:
$Ta$
is the Taylor number,
$Ra$
is the Rayleigh number,
$Ro=\sqrt{Ra/(PrTa)}$
is the input Rossby number,
$\widetilde{Ra}=RaTa^{-2/3}$
is a scaled Rayleigh number and
$\unicode[STIX]{x1D706}$
is the horizontal aspect ratio. The Prandtl number
$Pr=1$
in all simulations.
We now focus on the spectral statistics of the two different states. The kinetic energy spectra for each component of the flow, as defined by equations (2.7)–(2.9) (see also Favier et al. (Reference Favier, Silvers and Proctor2014) and Guervilly et al. (Reference Guervilly, Hughes and Jones2014)), are shown in figure 3. These spectra are averaged over depth and time. With no initial vortex, the 3D flow is dominant at virtually all scales, except for the smallest available wavenumber where most of the energy is contained in the 2D mode. This subdominant 2D flow is in equilibrium, meaning that there is a balance between viscous dissipation and baroclinic forcing (see below). There is no systematic growth of the 2D mode and no large-scale condensate is reached. Above the critical initial amplitude, however, the LSV is able to extract energy efficiently from the small-scale 3D flow, leading to rapid growth of the 2D mode with horizontal wavenumber
$k_{h}\leqslant 3$
, while the 3D flow remains largely unchanged at all scales.
Figure 3. Kinetic energy spectra averaged over depth (
$0<z<1$
) and time (
$0.3<t<0.4$
) as a function of the horizontal wavenumber
$k_{h}\equiv \sqrt{k_{x}^{2}+k_{y}^{2}}$
. (a)
$A=0$
. (b)
$A=1600$
.
We examine next the energetics of the 2D depth-invariant flow. Starting from (2.1)–(2.2), the governing equations for the purely horizontal 2D flow
$\langle \boldsymbol{u}\rangle _{z}$
are (see also Benavides & Alexakis (Reference Benavides and Alexakis2017))
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\langle \boldsymbol{u}\rangle _{z}}{\unicode[STIX]{x2202}t}+\langle \boldsymbol{u}\rangle _{z}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{h}\langle \boldsymbol{u}\rangle _{z}=-Pr\unicode[STIX]{x1D735}_{h}\langle p\rangle _{z}+Pr\unicode[STIX]{x1D735}_{h}^{2}\langle \boldsymbol{u}\rangle _{z}-\langle \boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }\rangle _{z},\end{eqnarray}$$
where
$\unicode[STIX]{x1D735}_{h}$
is the horizontal gradient operator. The last term in (3.1) corresponds to the forcing of the 2D mode by the 3D fluctuations. Taking the scalar product of (3.1) with
$\langle \boldsymbol{u}\rangle _{z}$
and volume-averaging leads to the energy density equation of the 2D mode
$$\begin{eqnarray}\frac{\text{d}K_{2D}}{\text{d}t}=\underbrace{\frac{Pr}{\unicode[STIX]{x1D706}^{2}}\iint \langle \boldsymbol{u}\rangle _{z}\boldsymbol{\cdot }\unicode[STIX]{x1D6FB}_{h}^{2}\langle \boldsymbol{u}\rangle _{z}\,\text{d}x\,\text{d}y}_{\displaystyle {\mathcal{D}}}+\underbrace{\left(-\frac{1}{\unicode[STIX]{x1D706}^{2}}\iint \langle \boldsymbol{u}\rangle _{z}\boldsymbol{\cdot }\langle \boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }\rangle _{z}\,\text{d}x\,\text{d}y\right)}_{\displaystyle {\mathcal{F}}}.\end{eqnarray}$$
The first term on the right-hand side corresponds to the viscous dissipation rate of the 2D flow while the last term represents the 2D energy production from the 3D fluctuations. Looking at (3.2), it is clear that the 2D flow is in equilibrium only when the viscous dissipation
${\mathcal{D}}$
is balanced by the 3D forcing
${\mathcal{F}}$
. It follows that growth of
$K_{2D}$
from a given equilibrium state can only be achieved by reducing the dissipation or increasing the forcing. Figure 4(a) shows that for
$A=0$
, the forcing is not zero but is exactly balanced by viscous dissipation. The sum
${\mathcal{D}}+{\mathcal{F}}$
oscillates rapidly around zero, a fact consistent with the quasi-constant value of
$K_{2D}$
observed for this case in figure 1(a). For
$A=1600$
, however, both dissipation and forcing increase in magnitude, and the sum is on average positive at least initially, corresponding to the growth of the 2D kinetic energy observed in figure 1(a). Note that the dissipation increases slowly with time while the baroclinic forcing term, while strongly fluctuating, remains quasi-constant (although it does grow very slightly). This increase in the dissipation is only observed for the 2D barotropic component; the dissipation associated with the 3D fluctuations remains largely unchanged (not shown). This is consistent with the condensation mechanism observed in purely 2D turbulence (Smith & Yakhot Reference Smith and Yakhot1994; Chertkov et al.
Reference Chertkov, Connaughton, Kolokolov and Lebedev2007; Gallet & Young Reference Gallet and Young2014) whenever no large-scale damping term is introduced to balance the inverse energy flux, and is a consequence of the slow growth of the dominant energy-containing scale.
Figure 4. Time evolution of the forcing and dissipation terms for the 2D mode, as defined by (3.2). The instantaneous value is plotted as a thin transparent line while the thick opaque line corresponds to the time average over a window spanning
$0.01$
vertical viscous time. (a)
$A=0$
. (b)
$A=1600$
. The horizontal dotted lines correspond to the time-averaged values in the reference case
$A=0$
.
Figure 5. (a) Amplitude of the 2D flow, vertical average of the 3D fluctuation amplitude, amplitude of the fluctuating forcing and its rate of working for
$A=0$
and
$A=1600$
at
$t=0.3$
. (b) Power spectrum of the forcing
$\langle \boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }\rangle _{z}$
at the same time.
At this stage, it is clear that the presence of an initially imposed vortex dipole somehow modifies the small-scale convective flow, enhancing energy transfer into the 2D mode. One possible explanation for this behaviour is that the 2D flow becomes significantly correlated with the 3D forcing, a correlation that can be quantified by looking at the angle
$\unicode[STIX]{x1D6FE}$
defined by
$\cos \unicode[STIX]{x1D6FE}=\langle \boldsymbol{u}\rangle _{z}\boldsymbol{\cdot }\langle \boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }\rangle _{z}/(|\langle \boldsymbol{u}\rangle _{z}|\,|\langle \boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }\rangle _{z}|)$
. However, our computations failed to reveal any significant differences in the probability density function of
$\cos \unicode[STIX]{x1D6FE}$
between cases with and without LSV (not shown), implying that the transfer enhancement cannot be explained by an increase in the correlation between the 3D forcing and the 2D flow. Figure 5 provides a clue. The figure reveals a clear imprint of the large-scale 2D structure
$\langle \boldsymbol{u}\rangle _{z}$
in the vertically-averaged 3D fluctuations
$\langle |\boldsymbol{u}^{\prime }|\rangle _{z}$
and in the forcing term
$\langle \boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }\rangle _{z}$
: the fluctuations are locally suppressed, likely by the strong, relatively ordered large-scale shear or vorticity associated with the vortex structure, implying enhanced correlations in the phases of the small-scale field. This suppression of the small-scale 3D fluctuations is consistent with the observed reduction in the Reynolds number (from
$Re\approx 535$
for
$A=0$
to
$Re\approx 512$
for
$A=1600$
) and of the Nusselt number (from
$Nu\approx 28.8$
for
$A=0$
to
$Nu\approx 27.6$
), as already noted by Guervilly et al. (Reference Guervilly, Hughes and Jones2014) for control parameters for which the LSV emerges spontaneously. It was also observed in thin-layer turbulence experiments (Xia et al.
Reference Xia, Byrne, Falkovich and Shats2011) although no subcritical behaviour was reported in this study. This suppression in turn enhances interactions between two large-
$k_{h}$
3D modes that transfer energy into a small-
$k_{h}$
2D mode, bypassing the standard inverse cascade and enhancing the power spectrum of the forcing at these large scales, as shown in figure 5(b). Note that, irrespective of the amplitude
$A$
, the forcing always peaks at
$k_{h}\approx 30$
, which is approximately the Taylor microscale characteristic of the small-scale vorticity field. The 2D energy is, however, very small at these scales (figure 3) and when
$A=0$
this small-scale forcing is in equilibrium with the dissipation (figure 4
a). For
$A=1600$
, however, the suppression of the small-scale fluctuations leads directly to enhanced forcing of the 2D flow at the box scale, allowing for a runaway growth. The non-local nature of this upscale energy transfer has already been discussed in previous studies (Favier et al.
Reference Favier, Silvers and Proctor2014; Rubio et al.
Reference Rubio, Julien, Knobloch and Weiss2014), and we expect similar non-local energy fluxes in our subcritical state, although this remains to be fully explored in future studies.
The resulting positive feedback mechanism differs in an important respect from that present in 2D turbulence driven by a small-scale prescribed stochastic force: in the present case the fluctuating force is itself dynamically modified by the LSV, and not just its rate of working as in Gallet & Young (Reference Gallet and Young2014) for example. In addition, the positive feedback revealed in figure 5 provides an indication that the LSV state will not in general spontaneously jump back to the lower LSV-free state: as soon as a transition starts to take place and the LSV is observed to grow, its backreaction on the 3D fluctuations favours energy transfer into the 2D component, leading to runaway dynamics which can only be arrested by viscous condensation at the box scale or by other large-scale effects not included in our simple model.
4 Discussion
We now discuss how the subcritical behaviour identified here depends on the various control parameters. Thus far we have focused on aspect ratio
$\unicode[STIX]{x1D706}=4$
, but similar simulations were performed with
$\unicode[STIX]{x1D706}=2$
(table 1). We observe very similar results: the initial LSV is amplified only above a finite critical amplitude. For this value of the aspect ratio, however, the existence of a non-local inverse energy transfer is not systematic. When some of the simulations were repeated with different random initial temperature perturbations, some cases exhibited an inverse cascade while others did not (figure 6
a) in a manner reminiscent of pipe flow (Darbyshire & Mullin Reference Darbyshire and Mullin1995). This is a consequence of the stochastic nature of the forcing term arising from the 3D fluctuations (figure 4), which can drastically affect the properties of the transition (Fauve et al.
Reference Fauve, Herault, Michel and Pétrélis2017). Indeed the transition seems much less robust for
$\unicode[STIX]{x1D706}=2$
than for
$\unicode[STIX]{x1D706}=4$
, an effect we ascribe to an increase in the amplitude of the fluctuations that arises from the smaller domain size, indicating that a reasonable scale separation between the initial vortex and the small-scale flow is required for a robust and reproducible transition. We confirmed this observation with simulations at
$\unicode[STIX]{x1D706}=6$
, for which a case with
$A=800$
decays while a case with
$A=1200$
eventually grows (although with the available computational resources neither of these cases can be run until saturation owing to the presence of much longer transients, see figure 6
b). This confirms that the observed transition is robust with respect to change in the aspect ratio of the vortex, although additional simulations (for example, decoupling the aspect ratio of the numerical domain and the initial size of the LSV) are required for a definitive conclusion about its effect on the threshold.
The results discussed above only apply for
$Ta=10^{8}$
and
$Ra=3\times 10^{7}$
. We chose
$Ta=10^{8}$
for numerical convenience, as increasing the Taylor number would be numerically prohibitive for the large-aspect-ratio domains and long-time dynamics considered here. Interestingly, the range of Rayleigh numbers for which a LSV is known to spontaneously emerge increases with the Taylor number (Favier et al.
Reference Favier, Silvers and Proctor2014; Guervilly et al.
Reference Guervilly, Hughes and Jones2014). It is reasonable, therefore, to assume that subcritical transitions will also be observed at higher Taylor numbers, and might even be more prominent. However, this hypothesis remains to be confirmed.
Figure 6. Time evolution of the kinetic energy density
$K_{2D}$
of the 2D flow for different initial vortex amplitudes
$A$
and aspect ratio
$\unicode[STIX]{x1D706}=2$
(a) and
$\unicode[STIX]{x1D706}=6$
(b). The shaded area corresponds to the values explored by
$K_{3D}$
.
We have also explored the effect of varying the Rayleigh number. First, it is known that a LSV cannot be sustained when the convection is insufficiently turbulent. For
$Ta=10^{8}$
,
$Ra\geqslant 5\times 10^{6}$
is required: when we repeated our simulations for
$Ra=4\times 10^{6}$
, varying the amplitude
$A$
of the vortex, we found no subcritical behaviour and all simulations eventually converged towards the same equilibrium dominated by 3D fluctuations. It appears that the subcritical transition is only present for sufficiently turbulent flows, well above the threshold of the linear instability. Finally, the same experiments were repeated for higher Rayleigh numbers (namely
$Ra=4\times 10^{7}$
and
$Ra=5\times 10^{7}$
) with results that are more subtle: for
$\unicode[STIX]{x1D706}=2$
, we did not observe any subcriticality, obtaining a 3D state irrespective of
$A$
. Surprisingly, for
$\unicode[STIX]{x1D706}=4$
and
$Ra=4\times 10^{7}$
, we recovered a LSV state that was not present for
$\unicode[STIX]{x1D706}=2$
. This fact provides further evidence that a clear scale separation between the vortex and the small-scale 3D flow favours a subcritical transition. A detailed study of the scaling of the critical amplitude
$A_{c}$
with
$Ta$
,
$Ra$
and
$\unicode[STIX]{x1D706}$
is, however, beyond the scope of this paper, although it represents an interesting line of investigation not only from a fundamental point of view, but also for possible applications to LSV structures in geophysical and astrophysical flows.
Our choice of initial conditions (2.4) is arbitrary and other options could be explored. In our moderate Rossby number simulations, symmetry breaking favouring cyclonic motions has been observed, in contrast with the symmetry between cyclonic and anticyclonic vortices found in the limit
$Ro\rightarrow 0$
(Julien et al.
Reference Julien, Rubio, Grooms and Knobloch2012). Our initial condition is, however, perfectly symmetric, requiring a significant transient phase to break up the large but unstable anticyclone, something that could be avoided if the calculations were to be initialized with a dominant cyclonic structure. In this respect, the reduced description of Julien et al. (Reference Julien, Rubio, Grooms and Knobloch2012) may be useful for exploring the detailed mechanism behind the subcritical transition discovered in this paper and, in particular, the presence of hysteresis in
$Ra$
, provided, of course, that the reduced equations exhibit similar behaviour.
We believe that the mechanism responsible for the observed subcritical transition is not specific to rotating convection and is likely to occur in other systems with multiple cascade scenarios, including thin-layer and magnetohydrodynamic turbulence. In particular, we emphasize that rotation is not required to observe LSV in 3D flows. Small-scale anisotropy, and its possible enhancement by large-scale flows, is the key and is present in all turbulent systems with multiple cascade scenarios, from thin-layer to magnetohydrodynamic turbulence. We also emphasize that the transition identified here separates two fully turbulent states, in contrast to the classical subcritical transition from laminar to turbulent shear flow in, for example, pipe flow (Darbyshire & Mullin Reference Darbyshire and Mullin1995; Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007). Some concepts developed in this field may nevertheless prove useful in the present context, such as the search for optimal perturbations (Kerswell Reference Kerswell2018).
Acknowledgements
This work was initiated during the workshop ‘Rotating Convection: from the Lab to the Stars’ organized and supported by the Lorentz Center at the University of Leiden (http://www.lorentzcenter.nl/). This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the programme ‘Investissements d’Avenir’ supervised by the Agence Nationale de la Recherche. Computations were also conducted with the support of the HPC resources of GENCI-IDRIS (grants no. A0020407543 and A0040407543), the Rocket High Performance Computing service at Newcastle University and the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). C.G. was supported by the UK Natural Environment Research Council under grant NE/M017893/1. B.F. acknowledges funding by the European Research Council under the European Union’s Horizon 2020 research and innovation programme through grant no. 681835-FLUDYCO-ERC-2015-CoG.
