2. Numerical methods

We study the canonical model system of rotating Rayleigh–Bénard convection, where the fluid is contained between a hot bottom and cold top boundary with a temperature difference $\Delta T$ , driving the fluid motion. The flow is simultaneously subject to a background rotation $\varOmega$ , aligned with the gravitational axis, leading to a Coriolis force. The system is described by three dimensionless numbers: the Rayleigh number ${Ra}$ , quantifying the strength of the buoyant forcing, the Ekman number ${Ek}$ , capturing the (inverse) strength of rotation, and finally the Prandtl number ${Pr}$ , describing the diffusive properties of the fluid, given by

(2.1) \begin{equation} {Ra}\equiv \frac {g\gamma \Delta \textit{TH}^3}{\nu \kappa }, \qquad {Ek}\equiv \frac {\nu }{2\varOmega H^2}, \qquad {Pr}\equiv \frac {\nu }{\kappa }. \end{equation}

Here, $g$ denotes the gravitational acceleration, $H$ is the domain height, while the fluid properties $\gamma$ , $\nu$ and $\kappa$ represent the thermal expansion coefficient, kinematic viscosity and thermal diffusivity of the fluid, respectively.

We employ the finite difference code developed by Verzicco & Orlandi (Reference Verzicco and Orlandi1996) (see also Ostilla-Monico et al. Reference Ostilla-Monico, Yang, van der Poel, Lohse and Verzicco2015) to perform DNS of the full Navier–Stokes and heat equations in the Oberbeck–Boussinesq approximation. The domain is a $L \times L \times H$ cuboid with Cartesian basis $(\boldsymbol{e}_x, \boldsymbol{e}_y, \boldsymbol{e}_z)$ that is horizontally periodic and has stress-free, constant temperature top and bottom boundaries. The full governing equations evolving the velocity field $\boldsymbol{u}$ and temperature field $T$ in time $t$ are given in convective units by

(2.2a) \begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + \left (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right )\boldsymbol{u} + \frac {\left (\textit{Pr}/\textit{Ra}\right )^{1/2}}{\textit{Ek}}\boldsymbol{e}_z\times \boldsymbol{u} &= -\boldsymbol{\nabla }p + T\boldsymbol{e}_z + \left (\frac {\textit{Pr}}{\textit{Ra}}\right )^{1/2}\Delta \boldsymbol{u} - \underbrace { \nu _\alpha \Delta _h^{-\alpha } \boldsymbol{u}_h}_{\substack {\text{adapted}\\\text{hypoviscosity}}}, \end{align}

(2.2b) \begin{align} \frac {\partial T}{\partial t} + \left (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right )T &= \frac {1}{\left (\textit{Ra}\textit{Pr}\right )^{1/2}}{\Delta} T, \end{align}

(2.2c) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0. \end{align}

Notice that in the last term of (2.2a ) we have included the adapted hypoviscosity. This is the large-scale friction term that takes care of dissipating the inverse energy flux. In a doubly periodic domain as in the (rotating) convection system at hand, we need to be somewhat careful how we can define an effective hypoviscosity, as we cannot simply use the isotropic inverse Laplacian that is commonly applied in triply periodic domains. Instead, motivated by the aim of dissipating the inverse cascade that resides in the 2D horizontal manifold, we apply a 2D variant of the hypoviscosity that acts on the horizontal velocity and the horizontal inverse Laplacian, applied across the full height of the domain. This is conveniently computed from a 2D Fourier transform in the horizontal periodic directions $\mathcal{F}_h\{\ldots \}$ , taking the horizontal velocity $\boldsymbol{u}_h\equiv u_x \boldsymbol{e}_x+u_y\boldsymbol{e}_y$ from the real space $\boldsymbol{u}_h(x,y,z,t)$ to the Fourier space $\hat {\boldsymbol{u}}_h(k_x,k_y,z,t)$ . Then, we can compute the hypoviscous term as

(2.3) \begin{equation} \nu _\alpha \Delta _h^{-\alpha } \boldsymbol{u}_h \equiv \nu _\alpha \mathcal{F}_h^{-1} \left\{\frac {1}{\left (k_x^2+k_y^2\right )^\alpha } \mathcal{F}_h \left \{ \boldsymbol{u}_h \right \} \right\}. \end{equation}

In practice, we need to avoid the singularity at $k_x=k_y=0$ , so we approximate

(2.4) \begin{equation} \nu _\alpha \Delta _h^{-\alpha } \boldsymbol{u}_h \approx \nu _\alpha \mathcal{F}_h^{-1} \left \{\frac {1}{\textrm {max}\left ((2\pi /L)^2,k_x^2+k_y^2\right )^\alpha } \mathcal{F}_h \left \{ \boldsymbol{u}_h \right \} \right \}, \end{equation}

clipping the horizontal wavenumber from below by the smallest non-zero wavenumber $2\pi /L$ .

A different approach for suppression of the LSV in rotating convection is put forward by Oliver et al. (Reference Oliver, Jacobi, Julien and Calkins2023), based on a modification of the nonlinear term. In Appendix A, we compare the different approaches for suppression of the LSV for an exemplary case.

We consider a fluid with ${Pr}=1$ , rotating at ${Ek}=10^{-4}$ and vary ${Ra}$ over the range in which the flow is known to transition into and out of the inverse cascade regime, in agreement with the parameter range considered in de Wit et al. (Reference de Wit, Aguirre Guzmán, Clercx and Kunnen2022a ) for the case with LSV (without hypoviscosity). See table 1 for the full set of input parameters considered in this work. The power $\alpha$ and strength $\nu _\alpha$ of the hypoviscosity are tuned in such a way that the full inverse energy flux is dissipated before it reaches the system scale, while simultaneously ensuring that its action remains mostly restricted to the largest scales (see also next section). We consider three series of runs with different horizontal domain sizes. The aspect ratio of the first series is chosen in agreement with de Wit et al. (Reference de Wit, Aguirre Guzmán, Clercx and Kunnen2022a ), employing $L=10\mathcal{L}_c$ , with $\mathcal{L}_c \equiv 4.8{Ek}^{1/3}H$ , the most unstable wavelength at onset of convection (Chandrasekhar Reference Chandrasekhar1961). The second and third series of runs progressively double the horizontal extent in order to increase the scale separation between the system scale and the scales that dominate the energy injection. The hypoviscosity parameters are kept constant throughout each series of runs. The grid is uniform in the horizontal directions, but non-uniform in the vertical direction, narrowing the grid close to the top and bottom to accurately resolve the boundary layers. The resolution $N_x \times N_y \times N_z$ is chosen with respect to the Kolmogorov scale $\eta$ , ensuring that the grid spacing never exceeds $2\eta$ , which is well below the threshold of $4\eta$ that is obtained for this code in Verzicco & Camussi (Reference Verzicco and Camussi2003).

Table 1. Input parameters that are used for the simulations in this work.

$^*$ Covered in de Wit et al. (Reference de Wit, Aguirre Guzmán, Clercx and Kunnen2022a ).

Simulations are initialised with zero velocity and a linear temperature profile connecting the top and bottom boundaries, superposed by a tiny statistical perturbation to initiate convection. Statistics are then collected after 1000 convective time units when a statistically steady state is reached at all scales.

3. Rotating convection with large-scale friction

Analogous to many other quasi-2D flow systems, rotating convection can give rise to a split cascade scenario (Alexakis & Biferale Reference Alexakis and Biferale2018). In this scenario, from the intermediate scales at which energy is injected, part of the energy is cascaded to small scales as in regular 3D turbulence, while another fraction of the energy is cascaded from the intermediate scales to larger scales, owing to the quasi-two-dimensionalisation. This gives rise to a net forward flux of energy at smaller scales, but a net inverse flux of energy at larger scales. While the forward flux is balanced by regular viscous dissipation, we shall here balance the inverse flux of energy with hypoviscosity at the large scales, to suppress the formation of the LSV.

To show the influence of the hypoviscosity on the flow we consider the kinetic energy spectrum

(3.1) \begin{equation} E(k) \equiv \frac {L}{2\pi } \int _0^H \sum _{k\leq \sqrt {k_x^2+k_y^2} \lt k+2\pi /L} \tfrac {1}{2} \left | \hat {\boldsymbol{u}}\left (k_x,k_y,z\right ) \right |^2 \mathrm{d}z. \end{equation}

This is shown for an exemplary case in figure 1. Indeed, we notice that the largest scales (smallest wavenumbers), in which the LSV resides, are strongly damped by the presence of the hypoviscosity, whereas the intermediate and small scales remain largely unaffected.

Figure 1. (a) Instantaneous kinetic energy spectra $E(k)$ for the cases with ${Ra}=8\times 10^6$ with and without hypoviscosity and for different horizontal domain sizes (see table 1). (b, c) The corresponding snapshots of kinetic energy $|\boldsymbol{u}|^2/2$ for the cases with $L/H=2.235$ at the mid-plane $z=H/2$ .

To tune the hypoviscosity exponent $\alpha$ and strength $\nu _\alpha$ we need to ensure that two conditions are met for all runs in each series.

1. (i) The hypoviscosity needs to be able to dissipate the full inverse energy flux before it condenses into an LSV. In practice, this means that the hypoviscosity must be strong enough such that the energy contained in the smallest (system scale) wavenumber does not dominate the spectrum. This can be translated into a condition on the spectral slope at the smallest wavenumber, demanding that $({\mathrm{d}\log E(k)})/({\mathrm{d}\log k}) \big |_{k=2\pi /L} \gt 0$ .

2. (ii) The action of the hypoviscosity needs to remain restricted to the smallest wavenumbers. We ensure this by checking that the spectral slope at small wavenumbers does not exceed $({\mathrm{d}\log E(k)})/({\mathrm{d}\log k}) \big |_{k=2\pi /L} \lt 2\alpha$ , such that the action of the hypoviscosity is indeed dominated by the smallest wavenumbers.

Then, if there is sufficient scale separation between the large scales where the hypoviscosity acts and the scales that dominate the energy injection, conservation of energy ensures that the rate of energy dissipation by the hypoviscosity will match exactly with the inverse energy flux, independent of the precise values of the hypoviscosity parameters, as long as the aforementioned conditions are met. This is verified numerically for an exemplary case in Appendix B.

Under these conditions, fluctuation fields that are dominated by intermediate- and small-scale structures remain unaffected by the presence of hypoviscosity, see for example figure 2.

Figure 2. Instantaneous snapshots of the temperature fluctuation field $T-\tilde {T}(z)$ , where $\tilde {T}(z)$ is the horizontally and temporally averaged temperature profile, for the cases with ${Ra}=8\times 10^6$ and $L/H=2.235$ without hypoviscosity (a) and with hypoviscosity (b).

4. Transition towards inverse cascade

With the hypoviscosity in place and tuned, we are well equipped to study the transition towards the inverse energy cascade under the influence of large-scale friction, in absence of the LSV. To that extent, we consider the strength of the inverse energy flux $\varepsilon _{ {inv}}$ as we vary the control parameter ${Ra}$ throughout the regime where the inverse cascade is known to occur.

As argued above, when the action of the hypoviscosity is correctly tuned, the energy dissipated by the hypoviscosity coincides with the inverse energy flux. We can thus measure the mean inverse energy flux in this case as

(4.1) \begin{equation} \varepsilon _{ {inv}}^{{(hypo)}} \equiv \left \langle \nu _\alpha \boldsymbol{u}_h \boldsymbol{\cdot }\Delta _h^{-\alpha } \boldsymbol{u}_h \right \rangle , \end{equation}

where we use $\langle \ldots \rangle$ to denote the average over the full spatial domain and time.

For comparison, we compute the inverse flux for the cases without hypoviscosity directly from the triadic nonlinear interactions that transfer kinetic energy from all scales into the system-scale 2D mode in which the LSV resides (Verma Reference Verma2019). The inverse flux is then computed as

(4.2) \begin{equation} \varepsilon _{ {inv}}^{ {(no\hbox{-}hypo)}} \equiv \sum _Q T_{\textrm {3D}} \left(K=\tfrac {2\pi }{L},Q\right) + T_{\textrm {2D}} \left(K=\tfrac {2\pi }{L},Q\right)\!, \end{equation}

where, following Rubio et al. (Reference Rubio, Julien, Knobloch and Weiss2014), Aguirre Guzmán et al. (Reference Aguirre Guzmán, Madonia, Cheng, Ostilla-Mónico, Clercx and Kunnen2020) and de Wit et al. (Reference de Wit, Aguirre Guzmán, Clercx and Kunnen2022a ), we distinguish baroclinic contributions

(4.3a) \begin{equation} T_{\textrm {3D}}(K,Q)\equiv -\left \langle \boldsymbol{u}^{\textrm {2D}}_K\boldsymbol{\cdot }\left(\overline {\boldsymbol{u}^{\textrm {3D}} \boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u}^{\textrm {3D}}_Q}\right)\right \rangle , \end{equation}

and barotropic contributions

(4.3b) \begin{equation} T_{\textrm {2D}}(K,Q)\equiv -\left \langle \boldsymbol{u}^{\textrm {2D}}_K\boldsymbol{\cdot }\left(\boldsymbol{u}^{\textrm {2D}}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^{\textrm {2D}}_Q\right)\right \rangle , \end{equation}

where the overbar $\overline {\vphantom {a}\ldots }$ denotes vertical averaging and with $\boldsymbol{u}^{\textrm {2D}} \equiv \overline {\boldsymbol{u}_h}$ and $\boldsymbol{u}^{\textrm {3D}}\equiv \boldsymbol{u}-\boldsymbol{u}^{\textrm {2D}}$ , and with subscript $K$ denoting Fourier filtering at mode $K$ .

Figure 3. The total inverse flux of kinetic energy $\varepsilon _{ {inv}}$ normalised by the total energy flux $\varepsilon = ({Nu}-1)/\sqrt {\textit{Pr}\textit{Ra}}$ as a function of ${Ra}$ for ${Ek}=10^{-4}$ and ${Pr}=1$ for the different series of runs with and without hypoviscosity. Panel (b) is a zoom of panel (a) that focuses on the cases with hypoviscosity, indicated by the dashed box in (a). For the runs without hypoviscosity, black diamonds indicate cases that show LSV formation. The figure reveals that, while the transition to the inverse cascade without hypoviscosity is discontinuous and hysteretic (lower hysteretic branch is depicted by the dotted line), the transition for the case with hypoviscosity is continuous.

The results are shown in figure 3. While the transition to the inverse cascade regime without large-scale friction, in presence of the LSV, is known to be discontinuous and hysteretic (Favier et al. Reference Favier, Guervilly and Knobloch2019; de Wit et al. Reference de Wit, Aguirre Guzmán, Clercx and Kunnen2022a ), we find that this behaviour is entirely different for the case with large-scale friction, in the absence of the LSV. Here, we observe that the transition in this case is in fact smooth and the discontinuity disappears (and thereby also the hysteretic behaviour, see Appendix C). This is in agreement with what is found in other hypoviscous quasi-2D systems (Seshasayanan et al. Reference Seshasayanan, Benavides and Alexakis2014; Benavides & Alexakis Reference Benavides and Alexakis2017; Pestana & Hickel Reference Pestana and Hickel2019; van Kan & Alexakis Reference van Kan and Alexakis2020). We furthermore find that the strength of the inverse cascade in absence of the LSV is much smaller than for the cases that have an LSV.

Figure 4. Time-averaged kinetic energy transport maps from 3D (a, c) and 2D (b, d) modes $Q$ to 2D modes $K$ , i.e. $T_{\mbox{3D}}(K, Q)$ and $T_{\mbox{2D}}(K, Q)$ , respectively, for the cases with ${Ra}=8\times 10^6$ and $L/H=2.235$ without hypoviscosity (a, b) and with hypoviscosity (c, d). Panels (a–d) also show the respective sums over the donating scales $Q$ , obtained as $\mathcal{T}_{\mbox{3D}}(K)\equiv \sum _Q T_{\mbox{3D}}(K,Q)$ and $\mathcal{T}_{\mbox{2D}}(K) \equiv \sum _Q T_{\mbox{2D}}(K,Q)$ (blue lines). Panel (e) depicts the total 2D energy flux $\varPi _{\mbox{2D}}(K) \equiv -\sum _{K'\lt K} \mathcal{T}_{\mbox{2D}}(K)$ , normalised by the total energy flux $\varepsilon$ .

A point of attention is raised by the observation that the inverse cascade does not seem to decay to zero for large Ra for the case with hypoviscosity. However, we argue that this is an effect of insufficient scale separation between the large scales at which the hypoviscosity acts and the scales that dominate the energy injection, such that the hypodissipation not only measures the inverse energy flux, but also the fraction of energy that is directly injected at these large scales. The energy injection is known to be dominated by features of typical size $\mathcal{L}_c$ , but the self-organising nature of the forcing injects energy at a broad range of scales around this dominant scale. The runs without hypoviscosity (red line) were performed at $L\approx 10\mathcal{L}_c$ and the runs with hypoviscosity I (blue line) employ the same aspect ratio in order to have a one-to-one comparison. However, to increase the scale separation, we have added the series of runs in the wider domains $L \approx 20 \mathcal{L}_c$ (II, green line) and $L \approx 40 \mathcal{L}_c$ (III, yellow line), which also allows for a less aggressive exponent of the hypoviscosity (see table 1). Indeed, we notice that the measured energy flux that is dissipated by the hypoviscosity remains mostly unchanged, but reaches a smaller value at large ${Ra}$ for the wider domains. We thus expect that the inverse flux will decay to values close to zero in the limit of increasingly large domain size as the scales become more cleanly separated.

The findings regarding the inverse cascade transition under the influence of hypoviscosity not only teach us something about what happens to convective systems that are subject to a large-scale friction mechanism, but through comparison, they also teach us about the influence of the LSV in the case without large-scale friction. Indeed, since the hypoviscosity only suppresses the LSV at large scales, leaving the other scales of the flow unaffected, we understand that the discontinuity observed in the inverse cascade transition in the case without large-scale friction is in fact a consequence of interactions between the LSV and the background turbulence.

To further understand the influence of the presence/absence of the LSV on the nature of the inverse cascade, we also consider the scale-to-scale kinetic energy transfer maps. By looking at $T_{\textrm {3D}}(K,Q)$ and $T_{\textrm {2D}}(K,Q)$ , we can study how baroclinic and barotropic kinetic energy is transferred from donating wavenumber $Q$ to receiving wavenumber $K$ . The results are provided in figure 4 for an exemplary case with and without hypoviscosity. It shows that in presence of the LSV (upper panels), the upscale energy transfer is highly non-local, with a direct transfer from all wavenumbers $Q$ into the largest wavenumber $K=2\pi /L$ , as is also known from Rubio et al. (Reference Rubio, Julien, Knobloch and Weiss2014), Aguirre Guzmán et al. (Reference Aguirre Guzmán, Madonia, Cheng, Ostilla-Mónico, Clercx and Kunnen2020) and de Wit et al. (Reference de Wit, Aguirre Guzmán, Clercx and Kunnen2022a ). In absence of the LSV (lower panels), on the contrary, we find that the upscale energy transfer has a much more local signature, with the energy transfer concentrating more around the diagonal of the map, indicating transfer between more like-sized scales. Note furthermore in figure 4(e), that the total 2D energy flux $\varPi _{\textrm {2D}}(K)$ towards low $K$ for the case with hypoviscosity is found to be $\approx 0.003 \varepsilon$ , comparable to the total inverse energy flux measured by the hypodissipation (see figure 3 b). This indicates that the inverse cascade that develops in the 2D manifold is fully dominating the measured total inverse flux in the case with hypoviscosity, while for the case without hypoviscosity (with LSV), the direct (non-local) baroclinic interactions between the large-scale 2D mode and the 3D modes are much stronger contributors (see figure 4 a).

Finally, we can also study what the influence of the LSV is on the total thermal throughput of the convective system, quantified through the Nusselt number $\textit{Nu}$ , given by

(4.4) \begin{equation} {Nu} \equiv \sqrt {\textit{Pr}\textit{Ra}}\langle T\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{e}_z \rangle + 1. \end{equation}

These results are depicted in figure 5. As first observed by Guervilly et al. (Reference Guervilly, Hughes and Jones2014), when the LSV emerges, the heat transfer slightly decreases. Without hypoviscosity, however, the $\textit{Nu}$ for the case without LSV can only be measured in the transient phase before the LSV emerges, at contrast with the statistically steady state obtained here. Comparing the runs with and without hypoviscosity, we find that the suppression of $\textit{Nu}$ by the LSV remains of the order of 3 %–8 % throughout the full range of ${Ra}$ in which the LSV is active for the ${Ek}=10^{-4}$ considered here. In the cases where there is no LSV, there is a good agreement between the runs with and without hypoviscosity.

Figure 5. (a) The average Nusselt number $\textit{Nu}$ as a function of ${Ra}$ for ${Ek}=10^{-4}$ and ${Pr}=1$ for the series of runs with and without hypoviscosity and $L/H=2.235$ . (b) The same data compensated by the $\textit{Nu}$ of the runs with hypoviscosity. For the runs without hypoviscosity, black diamonds denote cases that show LSV formation. This indicates that the LSV lowers the $\textit{Nu}$ by around 3 %–8 %.

It should be noted that the current simulations are conducted with stress-free boundary conditions. When considering no slip boundaries, the effect of Ekman pumping is known to increase $\textit{Nu}$ (Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014), which could counteract the suppressive effect of the LSV.

Oliver et al. (Reference Oliver, Jacobi, Julien and Calkins2023) have studied the scaling of $\textit{Nu}$ in absence of the LSV in the limit of ${Ek}\to 0$ using asymptotically reduced equations. Remarkably, they find in that case that $\textit{Nu}$ approaches the diffusivity-free scaling ${Nu}\sim {Pr}^{-1/2} {Ra}^{3/2} {Ek}^2$ . In the current work, the ${Ek}$ is not sufficiently low to retrieve such scaling.