Discontinuous transitions towards vortex condensates in buoyancy-driven rotating turbulence

Abstract

Using direct numerical simulations of rotating Rayleigh–Bénard convection, we explore the transitions between turbulent states from a three-dimensional (3-D) flow state towards a quasi-2-D condensate known as the large-scale vortex (LSV). We vary the Rayleigh number $Ra$ as control parameter and study the system response (strength of the LSV) in terms of order parameters, assessing the energetic content in the flow and the upscale energy flux. By sensitively probing the boundaries of the region of existence of the LSV in parameter space, we find discontinuous transitions and we identify the presence of a hysteresis loop as well as a memoryless abrupt growth dynamics. We show furthermore that the creation of the condensate state coincides with a discontinuous transition of the energy transport into the largest mode of the system.

1. Introduction

A hallmark feature of three-dimensional (3-D) turbulence is the forward energy cascade, transporting kinetic energy from large scales to ever smaller scales, as described by the celebrated theory of Kolmogorov (1941). In many geophysical and astrophysical flows, however, velocity fluctuations are largely suppressed in one direction as a consequence of, for example, confinement (Benavides & Alexakis 2017; Musacchio & Boffetta 2017, 2019), strong magnetic fields (Alexakis 2011; Seshasayanan, Benavides & Alexakis 2014) or fast rotation (Smith, Chasnov & Waleffe 1996; Seshasayanan & Alexakis 2018; Pestana & Hickel 2019; van Kan & Alexakis 2020), rendering the flow quasi-two-dimensional. This leads to the development of an inverse energy flux, akin to fully 2-D turbulence (Kraichnan 1967; Batchelor 1969), transporting energy from smaller to larger scales. Ultimately, this can lead to accumulation of kinetic energy at the largest available scales, followed by condensation into a vertically coherent large-scale vortex (LSV) structure at the domain size (see Alexakis & Biferale (2018) for a recent review). These LSVs are believed to play a crucial role in, for example, the formation of the Earth's magnetic field (Roberts & King 2013; Aurnou et al. 2015; Guervilly, Hughes & Jones 2015).

Following the framework that is brought forward in Alexakis & Biferale (2018), we aim to classify the transition from a 3-D forward cascading state to the condensate state by considering the behaviour of an order parameter that measures the strength of the LSV as a function of a control parameter of the flow throughout this transition. Then one can observe either a smooth transition, a continuous transition with discontinuous derivative or a discontinuous transition. This categorisation of the transition into the condensate state has shown to be an insightful approach in various other quasi-2-D flow systems (Alexakis 2015; Yokoyama & Takaoka 2017; Seshasayanan & Alexakis 2018; van Kan & Alexakis 2019).

These earlier works, however, have focused on more artificial, idealised flow models, where the turbulent forcing occurs at a single well-defined length scale. Here, we characterise the LSV transition in a natural, broadband-forced system of rotating convection, which is ubiquitous in geophysical and astrophysical flows. In this system, Favier, Guervilly & Knobloch (2019) have shown a bistability of an LSV with a non-LSV state at the same parameter values, depending on the initial conditions. The natural buoyant forcing over a broad range of scales obfuscates the strict separation of the injection, dissipation and condensation scales. Although one may expect that, in such natural and vigorously fluctuating turbulent systems, any transitions between different states are washed out and become gradual, we find that the transition towards the condensate state is in fact sharp and discontinuous.

Such abrupt transitions between turbulent states in a more general sense are a remarkable feature of fluid turbulence and have received much recent interest, being observed in various different flow settings, e.g. in torque measurements of Taylor–Couette and Von Kármán flows (Ravelet et al. 2004; Saint-Michel et al. 2013; Huisman et al. 2014), in states of stochastically forced 2-D and 3-D turbulence (Bouchet & Simonnet 2009; Iyer et al. 2017; Bouchet, Rolland & Simonnet 2019) and in reversals of the large-scale dynamics in thin layers (Sommeria 1986; Michel et al. 2016; Dallas, Seshasayanan & Fauve 2020). These types of abrupt transition are surmised to play an important role in, for example, climate research (Weeks et al. 1997; Jackson & Wood 2018; Herbert, Caballero & Bouchet 2020) and in understanding the geomagnetic reversal (Berhanu et al. 2007; Pétrélis et al. 2009).

2. Numerical approach

We consider the canonical system of rotating Rayleigh–Bénard convection, in which the flow is driven by buoyancy through a temperature difference $\Delta T$ between the bottom and top of the domain, whilst being simultaneously affected by strong background rotation $\varOmega$ along the vertical axis. The input space to this problem is described by three dimensionless numbers: the Rayleigh number $Ra=g\alpha \Delta T H^3/(\nu \kappa )$, quantifying the strength of the thermal forcing, the Ekman number ${Ek}=\nu /(2\varOmega H^2)$, representing the (inverse) strength of rotation, and the Prandtl number $Pr=\nu /\kappa$, containing the diffusive properties of the fluid. Here, $g$ denotes gravitational acceleration, $H$ is the domain height and $\alpha$, $\nu$ and $\kappa$ represent the thermal expansion coefficient, kinematic viscosity and thermal diffusivity of the fluid, respectively. The system is non-dimensionalised into convective units using $H$, $\Delta T$ and the free-fall velocity $U=\sqrt {g\alpha \Delta T H}$.

We solve the full governing set of Boussinesq Navier–Stokes equations through direct numerical simulation, employing the finite-difference code described in Verzicco & Orlandi (1996) and Ostilla-Monico et al. (2015) on a Cartesian grid with periodic sidewalls and stress-free boundary conditions at the top and bottom. For the width $D$ of the domain, we choose $D/H=10\mathcal {L}_c$ with $\mathcal {L}_c=4.8{Ek}^{1/3}$ the most unstable wavelength at onset of convection (Chandrasekhar 1961). The complete set of input parameters as well as the employed resolutions are provided in table 1. Note that we include the convective Rossby number $Ro=U/(2\varOmega H)={Ek}\sqrt {Ra/Pr}$ for reference. A validation of the grid resolution is provided in Appendix A.

Table 1. The different series of input parameters used in this work.

For numerical convenience, we use $Pr=1$ and ${Ek}=10^{-4}$, for which stable LSVs have been observed in earlier direct numerical simulations over a limited range of $Ra$ (Favier, Silvers & Proctor 2014; Guervilly, Hughes & Jones 2014; Favier et al. 2019). Upon increasing $Ra$ from the onset of convection, the two boundaries of the region of existence of the LSV are crossed. At the low-$Ra$ transition, the LSV develops as sufficient turbulent forcing is obtained to set up the upscale transport into the condensate, whereas at the high-$Ra$ transition, the LSV breaks down as too strong thermal forcing renders the flow insufficiently rotationally constrained, breaking the quasi-2-D conditions for upscale energy flux (Favier et al. 2014; Guervilly et al. 2014). We carry out a total of 46 runs at varying $Ra$, concentrated around both LSV transitions.

In order to analyse the LSV, we decompose the flow field $\boldsymbol {u}=u\boldsymbol {e}_x+v\boldsymbol {e}_y+w\boldsymbol {e}_z$ into a 2-D (vertically averaged) barotropic flow and a 3-D (depth-dependent) baroclinic flow (following Julien et al. 2012; Favier et al. 2014, 2019; Rubio et al. 2014; Aguirre Guzmán et al. 2020; Maffei et al. 2021), i.e. $\boldsymbol {u}=\boldsymbol {u}^{2D}+\boldsymbol {u}^{3D}$, where

(2.1a,b)\begin{equation} \boldsymbol{u}^{2D}=\overline{u}\boldsymbol{e}_x+\overline{v}\boldsymbol{e}_y,\quad \boldsymbol{u}^{3D}=\underbrace{(u-\overline{u})}_{u'}\boldsymbol{e}_x +\underbrace{(v-\overline{v})}_{v'}\boldsymbol{e}_y+w\boldsymbol{e}_z, \end{equation}

where the overbar $\overline {{\cdots }}$ denotes vertical averaging. Since the LSV is a largely vertically coherent structure (see figure 1), it resides primarily in the 2-D field, whereas the turbulent baroclinic fluctuations are captured by the 3-D field. Accordingly, we decompose the total kinetic energy $E_{tot}=\tfrac {1}{2}\langle |\boldsymbol {u}|^2\rangle$ into 2-D and (horizontal and vertical) 3-D contributions $E_{tot}=E^{2D}+E_H^{3D}+E_V^{3D}$ as

(2.2a–c)\begin{equation} E^{2D}=\tfrac{1}{2}\langle\overline{u}^2+\overline{v}^2\rangle,\quad E_H^{3D}=\tfrac{1}{2}\langle u'^2+v'^2\rangle,\quad E_V^{3D}=\tfrac{1}{2}\langle w^2\rangle, \end{equation}

where angular brackets $\langle \cdots \rangle$ represent an average over the full spatial domain. We also consider the energy spectrum of the 2-D flow from its Fourier transform $\hat {\boldsymbol {u}}^{2D}_{k_x k_y}$ as

(2.3)\begin{equation} \hat{E}^{2D}(K)=\sum_{K\leq\sqrt{k_x^2+k_y^2}< K+1} \frac{1}{2}|\hat{\boldsymbol{u}}^{2D}_{k_x k_y}|^2, \end{equation}

where we normalise the wavenumber $K$ by the box-size mode $2{\rm \pi} /D$, such that the LSV occupies the $K=1$ mode of the spectrum.

Figure 1. Snapshot of horizontal kinetic energy (in units of $U^2$) of the LSV-forming case $Ra=1.7\times 10^7$, truncated at three-quarter height to reveal a cross-section of the LSV.

3. Results and discussion

In figure 2, the different components of kinetic energy are provided over the range of considered $Ra$, crossing both LSV transitions. From the (largest mode of) 2-D energy that captures the LSV, we find a substantial discontinuity at both boundaries of the LSV state. At the high-$Ra$ transition, we find that this transition is hysteretic: the morphology of the flow depends on its history i.e. its initial conditions. These findings are in line with Favier et al. (2019), showing this bistability of an LSV and a non-LSV state for one case in this parameter range. To study this hysteresis loop, we initialise simulations using flow snapshots from a preceding $Ra$. For decreasing $Ra$ from a non-LSV state, the lower branch of the hysteresis loop is followed (open diamonds), whilst for increasing $Ra$ from an LSV state, the flow adheres to the upper hysteretic branch (filled squares), see figure 2. Note the remarkably large discontinuity in the lower hysteretic branch, where the flow transitions directly from a non-LSV state into nearly the strongest LSV forming state. Hysteresis in the LSV transition has also been observed in a rotating flow system with a sharp bandwidth Taylor–Green forcing (Yokoyama & Takaoka 2017).

Figure 2. Averaged kinetic energy components as a function of $Ra$. Filled squares: upper hysteretic branch of the high-$Ra$ transition; open diamonds: lower hysteretic branch. The cyan dashed-dotted line denotes the low-$Ra$ LSV transition, whilst the magenta and red lines denote the LSV transition of the lower and upper branch of the high-$Ra$ transition, respectively.

At the low-$Ra$ transition, on the other hand, no hysteresis is observed (using increments in $Ra$ of $\sim$2 %). Considering the cases directly above the LSV transition, however, we find that the growth of the LSV from a non-LSV state is non-monotonic: the flow shows an evident plateau during which the LSV has not yet developed, before finally growing relatively suddenly into the stable LSV state, see figure 3(a). The flow in this plateau state shows morphological similarities to the jet state observed in rectangular domains (Guervilly & Hughes 2017; Julien, Knobloch & Plumley 2018), but alternates between being predominantly in the $x$- and $y$-direction. We hypothesise that the peculiar growth behaviour found here signifies a memoryless abrupt growth process, much akin to the nucleation and growth type of dynamics that is observed in plentiful different systems throughout physics (Matsumoto, Saito & Ohmine 2002; Watanabe, Suzuki & Ito 2010; Garmann, Goldfain & Manoharan 2019; Metaxas et al. 2019). To substantiate this conjecture, we simulate an ensemble of 100 additional runs at $Ra=6\times 10^6$ with statistically perturbed initial conditions, using a reduced resolution of $128\times 128\times 72$ for computational feasibility. The hypothesised abrupt memoryless growth would then predict an exponential distribution of the waiting time spent in the metastable plateau state (stage B in figure 3a).

Figure 3. (a) Time evolution of horizontal kinetic energy of seven of the runs in the ensemble of the LSV-forming case $Ra=6\times 10^6$ close to the transition. We distinguish four stages: an initialisation phase (A), a metastable underdeveloped plateau state showing alternating jets (B), (quick) growth of the LSV (C) and the stably developed LSV state (D). The horizontal dashed line denotes the threshold to be crossed and sustained for 2000 convective time units ($U^{-1}H$) before the LSV is said to be completely developed, which defines the time point $t_{LSV}$, depicted by the blue cross. (b) Cumulative distribution of $t_{LSV}$ of the ensemble (blue solid line) with 95 % confidence bounds (blue dotted lines). It is fitted by an exponential distribution (red dashed line) as provided by (3.1). Inset shows the histogram corresponding to the same distribution.

To investigate the distribution of these waiting times, we define a time point $t_{LSV}$ at which the LSV is said to have stably developed once a threshold of horizontal kinetic energy is surpassed and sustained for 2000 convective time units, see figure 3(a). The obtained empirical cumulative distribution is then fitted with an exponential distribution

(3.1)\begin{equation} CDF(t_{LSV})=1-\exp\left(-\frac{t_{LSV}-t_0}{\underbrace{\mu(t_{LSV})-t_0}_{\tau_W}}\right), \end{equation}

where the fit parameter $t_0$ can be interpreted as the (fixed) contributions of the initialisation, growth and stable LSV stages (A, C and D in figure 3a). Here, $\mu (t_{LSV})$ denotes the mean of $t_{LSV}$, providing the maximum-likelihood estimate for the typical waiting time $\tau _W$ in this distribution. Figure 3(b) shows that there is excellent agreement between the hypothesised and the obtained distributions: the exponential distribution remains everywhere in between the 95 % confidence bounds of the empirical distribution. Signs of the exponentiality of waiting times have been observed in the LSV-forming system of sharp bandwidth forced thin-layer turbulence (van Kan, Nemoto & Alexakis 2019). Our findings indicate that, indeed, turbulent fluctuations randomly trigger the growth of the LSV, giving rise to this memoryless abrupt growth dynamics. Similar sudden growth behaviour is also found near the ends of the hysteresis loop in the high-$Ra$ transition, although an extended analysis of how the mean transition time evolves for changing $Ra$, repeating the ensemble average approach for all points in the considered parameter space, is currently out of computational reach for rotating convection.

Considering the observed hysteretic behaviour as well as the exponentially distributed waiting times, an analogy of this transition with first-order phase transitions in equilibrium statistical mechanics seems appropriate (Binder 1987). However, although it is conjectured that in particular the relatively more weakly dissipative large scales of the flow may show resemblance to thermal equilibrium states (Bouchet & Venaille 2012; Alexakis & Biferale 2018), ultimately, the chaotic and dissipative nature of turbulence makes the analogy with equilibrium statistical mechanics indirect. A more immediate interpretation of the transition in this fluctuating dynamical system would be in terms of nonlinear bifurcations. Then, the condensate transition as observed here can be interpreted as a subcritical nonlinear bifurcation, giving rise to two distinct attractors (indeed, the LSV state and the non-LSV state) which remain separated in phase space. Such noise-induced transitions between attractors are also known to be memoryless, yielding exponentially distributed waiting times (Kraut, Feudel & Grebogi 1999).

To understand how the LSV is energetically sustained, we compute the mode-to-mode kinetic energy transfer (see Dar, Verma & Eswaran 2001; Alexakis, Mininni & Pouquet 2005; Mininni, Alexakis & Pouquet 2005, 2009; Verma, Kumar & Pandey 2017; Verma 2019), distinguishing the 3-D to 2-D (baroclinic to barotropic) transport (following Rubio et al. 2014; Aguirre Guzmán et al. 2020)

(3.2a)\begin{equation} T_{3D}(K,Q)={-}\langle\boldsymbol{u}^{2D}_K\boldsymbol{\cdot}(\overline{\boldsymbol{u}^{3D} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^{3D}_Q})\rangle, \end{equation}

and the 2-D to 2-D (barotropic to barotropic) transport

(3.2b)\begin{equation} T_{2D}(K,Q)={-}\langle\boldsymbol{u}^{2D}_K\boldsymbol{\cdot} (\boldsymbol{u}^{2D}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}^{2D}_Q)\rangle, \end{equation}

describing the energetic transport into the Fourier-filtered 2-D flow field $\boldsymbol {u}^{2D}_K$ of mode $K$ from 3-D and 2-D modes $Q$ through triadic interactions arising from the advective term of Navier–Stokes. If $T_{3D},T_{2D} > 0$, there is a net transfer of kinetic energy from mode $Q$ to mode $K$ and vice versa. We also consider the transport into 3-D mode $K$ from the full (unfiltered) flow components $\mathcal {T}_{3D}(K)=\sum _Q T_{3D}(K,Q)$ and $\mathcal {T}_{2D}(K)=\sum _Q T_{2D}(K,Q)$, by summing over the donating modes $Q$.

The results for the shell-to-shell energy transfer throughout the LSV transitions are provided in figure 4. The two main energy fluxes are apparent in both $T_{3D}$ and $T_{2D}$. Near the diagonal, one can observe the direct forward cascade, transporting energy from $Q$ to slightly higher modes $K$. Note here that, while the $T_{2D}$ self-interaction must be symmetric by definition $T_{2D}(K,Q)=-T_{2D}(Q,K)$, this does not apply to $T_{3D}$ as it describes the energetic interactions between scales of the 2-D component and 3-D component of the flow. In the bottom row $K=1$, on the other hand, the upscale energy flux into the LSV can be appreciated. This energy flux is non-local: energy is transported directly from virtually all scales in the system into the box scale of the LSV, without participation of intermediate scales.

Figure 4. Time-averaged kinetic energy transport from 3-D (a,c,e,g) or 2-D (b,d, f,h) modes $Q$ to 2-D modes $K$, i.e. $T_{3D}(K,Q)$ and $T_{2D}(K,Q)$, respectively, in units of $U^3H^{-1}$. Blue lines denote $\mathcal {T}_{3D}(K)$ (a,c,e,g) and $\mathcal {T}_{2D}(K)$ (b,d, f,h). The low-$Ra$ transition is crossed from (a,b) ($Ra=5.6\times 10^6$) to (c,d) ($Ra=5.7\times 10^6$), whilst the high-$Ra$ transition of the upper hysteretic branch is crossed from (e, f) ($Ra=3.10\times 10^7$) to (g,h) ($Ra=3.13\times 10^7$), as also depicted in figure 5.

Figure 5 shows the energetic transport into the box-size mode as a function of $Ra$. Note that this considers the transport from the full, unfiltered 3-D and 2-D flow components into the LSV, that is, a sum over the donating scales in the bottom row $K=1$ of the transfer maps in figure 4. It makes clear that also the upscale transport into the LSV exhibits an evident discontinuous transition, both in $\mathcal {T}_{3D}(K=1)$ as well as, albeit to a lesser degree, in $\mathcal {T}_{2D}(K=1)$. Importantly, the figure indicates that it is the 3-D transport that is the dominant component in the forcing of the LSV.

Figure 5. Total 3-D $\mathcal {T}_{3D}(K=1)$ (blue) and 2-D $\mathcal {T}_{2D}(K=1)$ (red) transport of kinetic energy into the 2-D $K=1$ mode averaged over time as a function of $Ra$. Symbols and vertical lines are as in figure 2. Labels (a–h) depict the corresponding transfer maps in figure 4.

We argue that this upscale transport provides a clue to understanding the physical mechanism behind the observed sudden growth and hysteretic behaviour. As also detailed in Rubio et al. (2014) and Favier et al. (2019), the upscale transport contains a positive feedback loop, where the presence of the LSV itself enhances the upscale transport into the box-size mode. This agrees with our observation that the energetic transport into the LSV increases discontinuously as the LSV is created. The exact nature of how the LSV interacts energetically with its 3-D turbulent background is an interesting non-trivial question to explore in future work, which goes beyond considerations purely in Fourier space as well, by looking at individual vortex interactions, for example. The existence of such a positive feedback loop, however, seems intuitive: the predominantly cyclonic LSV locally increases the total vorticity (background rotation + flow vorticity), thereby strengthening the quasi-2-D conditions and hence the upscale transport. This mechanism allows the LSV to develop once its growth is triggered by rare turbulent fluctuations and allows the LSV to remain stably self-sustained over the hysteresis loop.

Since we consider specifically a cross-section of the full parameter space for varying $Ra$, the influence of other parameters on the transitional dynamics, such as ${Ek}$, $Pr$ and also the aspect ratio $D/H$ (since the domain width is the principal dynamical scale of the LSV), remains an open question. The morphology of the LSV in the asymptotically reduced model for ${Ek}\to 0$ is studied in its total parameter space in more detail in Maffei et al. (2021). For the transition specifically, however, one can argue that, as both attractors are expected to shift in continuous fashion through the phase space, only quantitative changes to the observed transitional dynamics are expected as the other control parameters are varied in vicinity to the computationally tractable values considered here. Nonetheless, the possibility that the discontinuity in the transition vanishes in a certain limit if both attractors would shift to coincide cannot be ruled out from the current simulations; the asymptotically reduced model seems appropriate to investigate this premise for the limit ${Ek}\to 0$.

4. Conclusions

We have described the fluid turbulence transition into a quasi-2-D condensate state in a natural broadband-forced system of rotating Rayleigh–Bénard convection, where the transition is sharply discontinuous, in spite of the lack of a clear separation of scales. We provide evidence of memoryless abrupt growth dynamics and hysteresis in these transitions, raising the picture of a double attractor phase space with a subcritical noise-induced transition between them. Furthermore, the correspondence of our findings with certain aspects of the LSV transition in other, artificially forced, flow systems, as remarked in the text, ultimately shows that this peculiar type of transition is a relevant and robust phenomenon that is expected to survive even in the geo- and astrophysically relevant flow of rotating convection, being one of the principal sources of fluid motion in nature.