1. Introduction
Buoyancy, i.e. the interplay of gravity with mass density inhomogeneities that are typically caused by thermal heterogeneities, is, howsoever introduced, the essential mechanism that drives heat transport in many natural flows. Examples for such natural convection processes can be found on Earth throughout its layers from mantle convection (Christensen Reference Christensen1995) over deep ocean convection (Maxworthy & Narimousa Reference Maxworthy and Narimousa1994) up to convection in its atmosphere (Atkinson & Wu Zhang Reference Atkinson and Wu Zhang1996), eventually determining local and global aspects of weather and climate.
Natural thermal convection flows reveal often a hierarchy of different flow structures such as clusters of clouds over the warm ocean in the tropics of Earth (Mapes & Houze Reference Mapes and Houze1993). The probably most prominent and thoroughly studied example of a hierarchy formation might be given by the solar convection zone in the outer 
$30\,\%$ of the Sun (Schumacher & Sreenivasan Reference Schumacher and Sreenivasan2020). In this case, so-called granules are superposed to larger flow structures termed supergranules: although both of them are driven by the heat flux at the solar surface (Rincon & Rieutord Reference Rincon and Rieutord2018; Schumacher & Sreenivasan Reference Schumacher and Sreenivasan2020), they offer very different lifetimes and horizontal extensions. Unfortunately, our understanding of such hierarchies’ origins is still far from complete (Hanson et al. Reference Hanson, Duvall, Birch, Gizon and Sreenivasan2020) and simpler set-ups of convection become necessary to improve it systematically.
Rayleigh–Bénard convection represents the simplest conceivable set-up and, thus, the paradigm of naturally forced, thermally driven turbulence. Here, a fluid layer of thickness 
$H$ is confined between a heated horizontal plane at the bottom and a cooled one at the top: because of the variation of density with temperature, the layer becomes unstable once subjected to gravity. As a result of intense research over the past decades, it is well-known that such convection systems organise themselves even in the fully turbulent regime into prominent long-living large-scale flow structures. Although clearly distinguishable from the universal smaller-scale turbulence or fluctuations on significantly shorter time scales, the nature of these large-scale flow structures is complex and depends instead on various external factors such as the strength of the thermal driving, the working fluid or the presence of additional physical mechanisms (Vieweg Reference Vieweg2023).
Interestingly, only very recent research identified thermal boundary conditions as the key factor in determining the nature of these long-living large-scale flow structures given a horizontally extended domain. In a nutshell, either so-called turbulent superstructures with characteristic horizontal extensions of 
$\varLambda \sim O(H)$ form (Pandey, Scheel & Schumacher Reference Pandey, Scheel and Schumacher2018; Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018; Krug, Lohse & Stevens Reference Krug, Lohse and Stevens2020; Käufer et al. Reference Käufer, Vieweg, Schumacher and Cierpka2023), or a so-called gradual supergranule aggregation takes place that might result in a domain-sized flow structure with 
$\varLambda \gg O (H)$ if not being interrupted by additional mechanisms such as rotation around the vertical axis (Vieweg, Scheel & Schumacher Reference Vieweg, Scheel and Schumacher2021a; Vieweg et al. Reference Vieweg, Scheel, Stepanov and Schumacher2022). Although the former establish whether the horizontal planes offer uniform temperatures (so-called Dirichlet conditions), the latter correspond to planes that prescribe a uniform vertical temperature gradient or, in other words, a spatially constant heat flux (Neumann conditions). Furthermore, the supergranules are superposed to significantly smaller (yet large-scale) granule-like flow structures, so a hierarchy of different horizontally extended flow structures may establish even in a simple turbulence configuration. This effect of thermal boundary conditions extends also to the Lagrangian material transport and the present coherent features in the flow (Vieweg et al. Reference Vieweg, Schneide, Padberg-Gehle and Schumacher2021b; Schneide et al. Reference Schneide, Vieweg, Schumacher and Padberg-Gehle2022; Vieweg et al. Reference Vieweg, Klünker, Schumacher and Padberg-Gehle2024). Remarkably, these different self-organisations of the flows persist across the entire numerically accessible range of Rayleigh numbers 
${Ra} \lesssim 10^{8}$ (which quantify the strength of the thermal driving) (Vieweg et al. Reference Vieweg, Scheel and Schumacher2021a, Reference Vieweg, Scheel, Stepanov and Schumacher2022). Hence, the way how buoyancy effects are prescribed at the planes or boundaries seems to eventually determine the large-scale nature of the flows in between.
Exceeding a critical value of thermal driving, the buoyancy-induced destabilisation leads to an onset of convection. Although this critical value depends on the thermal boundary conditions, the latter modify in particular the horizontal extension of the emerging flow structures. In more detail, this primary instability leads to the emergence of convection rolls that exhibit predominantly one particular horizontal extension: depending on the mechanical boundary conditions, the corresponding critical wave numbers are 
$k_{{h,crit}} = [2.22, 3.13]$ (Rayleigh Reference Rayleigh1916; Pellew & Southwell Reference Pellew and Southwell1940) or 
$k_{{h,crit}} = 0$ (Hurle et al. Reference Hurle, Jakeman and Pike1967) for applied uniform temperatures or vertical temperature gradients, respectively. This latter value is further supported by the secondary instability slightly above the onset of convection, revealing that ‘each mode is unstable to one of longer wavelength than itself, so that any long box will eventually contain a single roll’ (Chapman & Proctor Reference Chapman and Proctor1980). In other words, any convection roll (of arbitrary size) is, at least slightly above the onset of convection, unstable to a more extended convection roll if buoyancy is introduced via a constant heat flux. Given that this result is obtained from a nonlinear evolution equation for the two-dimensional leading-order temperature perturbation, it is remarkable that a three-dimensional leading Lyapunov vector stability analysis discovered for a Prandtl number 
${Pr} = 1$ (which defines the working fluid) that the gradual supergranule aggregation is, even far beyond the onset of convection, driven by such secondary instabilities (Chapman & Proctor Reference Chapman and Proctor1980; Vieweg et al. Reference Vieweg, Scheel and Schumacher2021a), see also figure 1. Once the numerically finite horizontal extent of the domain is reached, the final statistically stationary state resembles essentially a finite-size relic of critical mode and thus shares similarities with the primary instability. Crucially, the latter is independent of the working fluid, whereas secondary and subsequent instabilities depend at least in the classical case of prescribed temperatures strongly on the working fluid (Busse Reference Busse1978, Reference Busse2003). In the case of a prescribed heat flux, the authors studying secondary instabilities stated that their ‘results hold quite generally for all Prandtl numbers’ (Chapman, Childress & Proctor Reference Chapman, Childress and Proctor1980) but simultaneously ‘do not expect the theory to remain accurate for very small 
${Pr}$’ (Chapman & Proctor Reference Chapman and Proctor1980). As the final supergranule results from the preceding transient supergranule aggregation, clarifying this uncertainty becomes crucial especially due to the strongly varying Prandtl numbers in geophysical and astrophysical flows.
Figure 1. Gradual supergranule aggregation. Although secondary instabilities are essential for the transient growth of the supergranules, the final flow resembles a state described by the primary instability (Hurle, Jakeman & Pike Reference Hurle, Jakeman and Pike1967). This time series visualises a flow at 
${Pr} = 10^{-2}$ (table 1) across the entire horizontal cross-section of aspect ratio 
$\varGamma = 60$ at 
$z_{0} = 1 - \delta _{T} /2$ with the thermal boundary-layer thickness 
$\delta _{T} = 1 / (2 {Nu})$.
In the present work, we conduct direct numerical simulations across an extended range of fluids applicable to geophysical and astrophysical convection systems while prescribing constant vertical temperature gradients at the horizontal top and bottom planes. Providing extraordinarily long evolution times of up to the order of 
$O(10^{4})$ convective time units, we confirm that supergranule aggregation is an omnipresent feature independently of the working fluid. Despite its involved hierarchy of different large-scale flow structures, the global heat and momentum transport of the flows shares clear analogies with the complementary turbulent superstructures that manifest in the case of applied constant temperatures. Interestingly, the bulk stratification might manifest qualitatively differently depending on the working fluid.
2. Numerical method
We consider the simplest conceivable scenario of convection based on the Oberbeck–Boussinesq approximation (Oberbeck Reference Oberbeck1879; Boussinesq Reference Boussinesq1903) where the key idea is that the dependence of material parameters on ‘pressure is unimportant and that even the variation with temperature may be disregarded except in so far as it modifies the operation of gravity’ (Rayleigh Reference Rayleigh1916). As a consequence, the mass density 
$\rho$ becomes a linear function of only the temperature when it acts together with gravity but is constant or incompressible otherwise.
The three-dimensional equations of motion are solved by the spectral-element method Nek5000 (Fischer Reference Fischer1997; Scheel, Emran & Schumacher Reference Scheel, Emran and Schumacher2013). The equations are made dimensionless based on the layer height 
$H$ and the applied constant vertical temperature gradient 
$\beta$ at the plates, resulting in 
$\beta H$ as the characteristic temperature scale. Together with the free-fall inertial balance, the free-fall velocity 
$U_{{f}} = \sqrt {g \alpha \beta H^{2}}$ and free-fall time scale 
$\tau _{{f}} = 1 / \sqrt {\alpha g \beta }$ establish as characteristic units. Here, 
$\alpha$ is the volumetric thermal expansion coefficient at constant pressure and 
$g$ the acceleration due to gravity. This translates the equations eventually into
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u} ={-} \boldsymbol{\nabla} p + \sqrt{\frac{{Pr}}{{Ra}}} \nabla^{2} \boldsymbol{u} + T \boldsymbol{e}_{z}, \end{gather}
\begin{gather}\frac{\partial T}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}) T = \frac{1}{\sqrt{{Ra} {Pr}}} \nabla^{2} T \end{gather}
with 
$\boldsymbol {u}$, 
$T$ and 
$p$ representing the velocity, temperature and pressure field, respectively. The relative strength of the individual terms in these equations is controlled by the Rayleigh and Prandtl number,
\begin{equation} {Ra} := \frac{\alpha g \beta H^{4}}{\nu \kappa} \quad \textrm{and}\quad {Pr} := \frac{\nu}{\kappa}, \end{equation}
only. The quantities 
$\nu$ and 
$\kappa$ denote the viscosity and thermal diffusivity, respectively, and thus define the strength of molecular diffusion processes.
Independently of 
${Ra}$ and 
${Pr}$, (2.1)–(2.3) are complemented by a three-dimensional domain with a square horizontal cross-section 
$A = \varGamma \times \varGamma$ and an aspect ratio 
$\varGamma := L / H = 60$ where 
$L$ is the horizontal periodic length of the domain. We apply at the top and bottom planes mechanical free-slip boundary conditions
\begin{equation} u_{z} (z \in \{0, 1\}) = 0,\quad \frac{\partial u_{x, y}}{\partial z} (z \in \{0, 1\}) = 0, \end{equation}as well as thermal constant heat flux boundary conditions
\begin{equation} \frac{\partial T}{\partial z} (z \in \{0, 1\}) ={-} 1. \end{equation}In spite of our interest in large-scale flow structures, our direct numerical simulations resolve all dynamically relevant scales of the flows ranging from the domain size down to the dissipation scales based on a (refined) Grötzbach criterion (Scheel et al. Reference Scheel, Emran and Schumacher2013). These dissipation scales are given by the so-called Kolmogorov and Batchelor scale (Batchelor Reference Batchelor1959; Kolmogorov Reference Kolmogorov1991; Sreenivasan Reference Sreenivasan2004),
\begin{equation} \eta_{{K}} := \frac{{Pr}^{3/8}}{{Ra}^{3/8} \varepsilon^{1/4}}\quad \textrm{and} \quad \eta_{{B}} := \frac{\eta_{{K}}}{\sqrt{{Pr}}}, \end{equation}
for the velocity and scalar temperature field, respectively, where 
$\varepsilon := ( 1 / 2 ) \sqrt {{Pr} / {Ra}} [(\boldsymbol {\nabla } \boldsymbol {u} ) + (\boldsymbol {\nabla } \boldsymbol {u} )^{{\rm T}}]^{2}$ represents the kinetic energy dissipation rate. Note that whereas the Batchelor scale 
$\eta _{{B}} \leq \eta _{{K}}$ applies for 
${Pr} \geq 1$, the Corrsin scale 
$\eta _{{C}} := \eta _{{K}} / {Pr}^{3/4}$ (Corrsin Reference Corrsin1951) is here not of particular interest as it applies only at 
${Pr} \leq 1$ where 
$\eta _{{C}} \geq \eta _{{K}}$.
3. Results
In contrast to our previous work (Vieweg et al. Reference Vieweg, Scheel and Schumacher2021a), we fix here the Rayleigh number 
${Ra} \approx 2.0 \times 10^{5}$ but vary instead the Prandtl number 
${Pr} \in [ 10^{-2}, 10^{2} ]$ across 
$4$ orders of magnitude centred around 
${Pr} = 1$. The precise parameters are summarised for all our simulation runs in table 1.
Table 1. Simulation parameters of the direct numerical simulations at different Prandtl numbers 
${Pr}$: the Rayleigh number 
${Ra} = 203, 576$, aspect ratio 
$\varGamma = 60$ and free-slip as well as constant heat flux boundary conditions are applied for all runs. The table also contains the total number of spectral elements 
$N_{{e}}$ in the simulation domain, the polynomial order 
$N$ on each spectral element, the total runtime of the simulation 
$t_{{r}}$ in units of the corresponding free-fall times 
$\tau _{{f}}$, a subsequent translation of these runtimes into vertical diffusion times 
$\tau _{\nu,\kappa }$, as well as the resulting integral length scale 
$\varLambda _{T}$ of the temperature field at midplane, Nusselt number 
${Nu}$, Reynolds number 
${Re}$ and the mean Kolmogorov and Batchelor scale, 
$\langle \eta _{{K}} \rangle _{V, t}$ and 
$\langle \eta _{{B}} \rangle _{V, t}$, respectively. 
$\varLambda _{T}$, 
${Nu}$, 
${Re}$, 
$\langle \eta _{{K}} \rangle _{V, t}$ and 
$\langle \eta _{{B}} \rangle _{V, t}$ are obtained from the last 
$500 \tau _{{f}}$ (
$5 \tau _{{f}}$ for 
${Pr} = 10^{-2}$, 
$1000 \tau _{{f}}$ for 
${Pr} = 10^{2}$) of each simulation run.
3.1. Ubiquitous gradual supergranule aggregation
Initialised with its fluid at rest possessing a randomly perturbed linear diffusive equilibrium profile, i.e. 
$\boldsymbol {u} ( t = 0 ) = 0$ and 
$T ( t = 0 ) = T_{{lin}} + \varPsi$ together with 
$T_{{lin}} := 1 - z$ and 
$0 \leq \varPsi (\boldsymbol {x}) \leq 10^{-3}$ (Vieweg Reference Vieweg2023; Vieweg et al. Reference Vieweg, Klünker, Schumacher and Padberg-Gehle2024), every simulation is run as long as necessary to indicate a stationarity of the large-scale flow structure formation. This can be captured, for instance, by (i) the thermal variance 
$\varTheta _{{rms}}$ with the temperature deviation 
$\varTheta = T - T_{{lin}}$ or (ii) the integral length scale (Parodi et al. Reference Parodi, von Hardenberg, Passoni, Provenzale and Spiegel2004) of the temperature field
\begin{equation} \varLambda_{T} (z_{0}, t) := 2 {\rm \pi}\frac{\displaystyle\int_{k_{{h}}} [E_{TT} (k_{{h}}, z_{0}, t) / k_{{h}}]\,{\rm d} k_{{h}}}{\displaystyle\int_{k_{{h}}} E_{TT} (k_{{h}}, z_{0}, t)\,{\rm d} k_{{h}}} \end{equation}
based on the azimuthally averaged Fourier energy spectrum at midplane, 
$E_{TT} (k_{{h}}, z_{0} = 0.5, t)$, as shown in Vieweg et al. (Reference Vieweg, Scheel, Stepanov and Schumacher2022). Note here that neither the Reynolds nor the Nusselt number (see (3.2) and (3.3), respectively) reflect the transient large-scale structure formation properly (Vieweg et al. Reference Vieweg, Scheel and Schumacher2021a, Reference Vieweg, Scheel, Stepanov and Schumacher2022). Running our simulations reveals two particularly interesting results.
First, the gradual supergranule aggregation, first reported by Vieweg et al. (Reference Vieweg, Scheel and Schumacher2021a), sets in even beyond 
${Pr} = 1$ at all accessible Prandtl numbers as both 
$\varTheta _{{rms}}$ and 
$\varLambda _{T}$ increase over time, see also figure 2. Yet, the varying diffusivities affect the pace of the dynamics and thus the necessary simulation runtime 
$t_{{r}}$ to reach a statistically stationary large-scale pattern size, see table 1. Although 
$t_{{r}}$ is by far largest for the upper investigated limit of 
${Pr}$, we find a similar trend in the opposing lower limit. This observation confirms that the efficiency of the aggregation process depends on the interplay of the velocity and temperature field, being in line with our previous results (Vieweg et al. Reference Vieweg, Scheel, Stepanov and Schumacher2022) which trace the (thermal) supergranule aggregation basically back to an advective transfer of thermal variance. Consequently, these runtimes do not support any relation to the diffusive time scales 
$\tau _{\nu } = H^{2} / \nu \equiv \sqrt {{Ra} / {Pr}} \tau _{{f}}$ and 
$\tau _{\kappa } = H^{2} / \kappa \equiv \sqrt {{Ra} {Pr}} \tau _{{f}}$ as contrasted in table 1. Interestingly, while more simulations are required to draw firm conclusions on the interplay of diffusion processes concerning the pace of the aggregation process, the increase of necessary runtime is clearly larger in the direction 
${Pr} \rightarrow \infty$.
Figure 2. Signs of the transient supergranule aggregation. (a) While neither 
${Nu}$ nor 
${Re}$ are affected significantly, (b) 
$\varLambda _{T}$ and 
$\varTheta _{{rms}}$ do indicate the transient supergranule aggregation. The data correspond to 
${Pr} = 10^{-2}$, see also figure 1.
Second, this process ceases, independently of 
${Pr}$, only once the horizontal domain size is reached, implying that thermal variance has significantly aggregated on the scale of the horizontal domain size. Consequently, the integral length scale 
$\varLambda _{T}$ converges in any simulation run towards 
$\varGamma$ as indicated by table 1 and figure 2(b). Figure 3 visualises the temperature and vertical velocity field in horizontal planes within the upper thermal boundary layer for these final states of the flows. In particular, figure 3(a,i,k,o,q) depict the temperature fields across the entire horizontal cross-sections of the domains, whereas figure 3(c,f,m,p,r) exemplary contrast them to the velocity field with respect to its vertical component. The circumstance that the supergranules grow in every run without any upper physical limit confirms that the secondary instability mechanism (Chapman et al. Reference Chapman, Childress and Proctor1980; Chapman & Proctor Reference Chapman and Proctor1980) rules the formation of long-living large-scale flow structures even far beyond the onset of convection independently of 
${Pr}$.
Figure 3. Supergranulation across four orders of Prandtl numbers. Although the velocity field exhibits successively smaller features for decreasing Prandtl numbers 
${Pr}$, the supergranule aggregation can still easily be observed in the temperature field. Panels (a,c,f,i,k,m,o,p,q,r) visualise the entire cross-section at 
$z_{0} = 1 - \delta _{T} /2$. To highlight the vast scale-separation between the temperature and (vertical) velocity field for small 
${Pr}$, panels (b,d,g,j,l,m) show enlarged regions of interest of size 
$15 \times 15$. Panels (e,h) underline this fact by additional magnifications of regions of size 
$4 \times 4$.
Albeit the gradual supergranule aggregation seems to be a ubiquitous feature across all covered fluids, the variation of the Prandtl number still modifies other aspects of the flow. While they display well-ordered stems of localised up- and down-flow regions for large 
${Pr}$, they become increasingly disordered for increasingly smaller 
${Pr}$ due to the reduced importance of molecular friction. Consequently, the ranges of observable scales or details diverge when comparing the temperature and vertical velocity field: this is highlighted in figure 3 by magnifications of fractions of the flows. In the case of 
${Pr} = 1$, both fields offer an equivalent richness of details which is shown in figure 3(l,n). This changes once the Prandtl number moves off unity and the diffusivities of momentum and the scalar temperature or the mean Kolmogorov (Kolmogorov Reference Kolmogorov1991) and Batchelor (Batchelor Reference Batchelor1959) scales differ. On the one hand, the temperature field becomes successively diffuse or imprecise for increasingly smaller 
${Pr}$, cf. figure 3(b,j,l). On the other hand, the velocity field becomes simultaneously successively more chaotic as directly contrasted in figure 3(d,g,n). The tremendous scale separation between the two fields is ultimately highlighted by further magnifications of even smaller regions in figure 3(e,h), underlining the vast complexity of low-
${Pr}$ thermal convection flows. Table 1 quantifies this visual scale separation by including the mean Kolmogorov and Batchelor scale for each simulation.
The increasing local disparity of the temperature and velocity field due to the different time scales of the underlying diffusion processes indicates that the impact of a variation of 
${Pr}$ on the global transport of momentum and heat should be investigated in the following.
3.2. Global transport properties and the role of stratification
An alternative perspective on the response of the dynamical system on its buoyancy-induced forcing is provided by its global momentum and heat transport as can be measured by the Reynolds and Nusselt number, respectively. While the former is given by
\begin{equation} {Re} (t) := \sqrt{\frac{{Ra}}{{Pr}}} u_{{rms}} \quad \textrm{with } u_{{rms}} := \sqrt{\langle \boldsymbol{u}^{2} \rangle_{V}}, \end{equation}the latter quantifies the strength of convective heat transport, by comparing the total heat transport across the fluid layer to a state of pure heat conduction, and results (in the present case of an applied constant heat flux) in (Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002)
\begin{equation} {Nu} (t) = \frac{1}{{\rm \Delta} T_{{N}}} \quad \textrm{with } {\rm \Delta} T_{{N}} := \langle T(z = 0) - T(z = 1) \rangle_{A} \leq 1, \end{equation}
where 
${\rm \Delta} T_{{N}}$ is the dynamically manifesting mean temperature drop across the fluid layer.
Figure 4 visualises via dark markers the dependence of these global transport measures on the Prandtl number for the final flow states, see again figure 3. On the one hand, the Reynolds number can be found to increase steadily when the Prandtl number is decreased. This is in accordance with the vanishing role of viscous diffusion, allowing for higher velocities and leading to successively more inertial flows. As this holds for the entire range of covered Prandtl numbers, it implies that the flow laminarises for 
${Pr} \gg 1$. On the other hand, the Nusselt number shows a more complex behaviour. For decreasing Prandtl numbers in the range 
${Pr} \lesssim 1$, thermal diffusion gains relevance as the disorder in the flow intensifies (see 
${Re}$). In contrast, 
${Nu}$ stagnates for 
${Pr} \gtrsim 1$: this effect might be induced by the full nesting of the thermal boundary layer into the viscous one (Chillà & Schumacher Reference Chillà and Schumacher2012) (the latter of which might be estimated to be 
$\delta _{u} \sim {Pr} \delta _{T}$ based on diffusion arguments), so buoyancy effects are suppressed or protracted by viscous diffusion and thermal plumes detach less frequently.
Figure 4. Global momentum and heat transport for different fluids. (a) While the global momentum transport increases with decreasing 
${Pr}$, (b) the convective heat transport loses importance (relative to purely diffusive heat transport) only for 
${Pr} \lesssim 1$. The dark markers correspond to supergranule data from the late state of this study's flows as described by table 1. In contrast, the bright markers represent turbulent superstructure data (i.e. different thermal boundary conditions) as outlined in the discussion. In a nutshell, the present study differs from 
$\ddagger$ (Pandey et al. Reference Pandey, Krasnov, Sreenivasan and Schumacher2022) as follows: thermal Neumann boundary conditions vs Dirichlet conditions, horizontally periodic domain of 
$\varGamma = 60$ vs closed box of 
$\varGamma = 25$, mechanical free-slip boundary conditions at the top and bottom planes vs no-slip conditions. Note that the series at 
${Ra} \approx 2.0 \times 10^{5}$ and 
${Ra}_{{D}} = 10^{5}$ can be related (Vieweg Reference Vieweg2023).
Thermal plume detachments are fundamentally caused by the applied (inverse or) unstable density stratification introduced at the heated bottom and cooled top plane. These ascending and descending plumes leave consequently the boundary layers and travel, driven by buoyancy, into or even through the bulk, leading to turbulent mixing once the flow is sufficiently inertial. Remarkably, our previous study (Vieweg et al. Reference Vieweg, Scheel and Schumacher2021a) observed a slightly stable density stratification for any constant heat flux-driven convection flow in the bulk region independently of 
${Ra}$ given 
${Pr} = 1$. In other words, the flow structures established a density stratification that was counter-directed to the applied one. Although the strength of this stratification decreased with increasing 
${Ra}$, it remained stable for all accessible 
${Ra} \lesssim 10^{8}$. In the following, we address the question of whether such a stable stratification is a unique feature of every flow that exhibits the effect of supergranule aggregation.
Therefore, we contrast the temperature profiles of all present runs in figure 5(a). Note that the temperature fields are re-scaled here via 
$T_{{rs}} = (T - \langle T \rangle _{V})/{\rm \Delta} T_{{N}} + \langle T \rangle _{V}$ (which does not affect the stratification properties) to allow for a direct comparison. Unlike in our previous study, we find here stably as well as unstably stratified bulks despite the presence of supergranules for any 
${Pr}$. While it is stable for 
${Pr} \geq 1$ and converges for 
${Pr} \gtrsim 7$, it is increasingly unstable for successively smaller Prandtl numbers 
${Pr} < 1$. Interestingly, these trends coincide with the above findings regarding the scaling of 
${Nu} ({Pr})$, suggesting a relation of the bulk stratification with plume detachments. Hence, a stable stratification in the bulk is no omnipresent result of the emergence of supergranular flow structures, while the potentially forming local peaks in the temperature profile can be seen as the consequence of a competition between the protracted overshooting thermal plumes and the opposite boundary layers close to the top and bottom plane.
Figure 5. Stratification and the dominance of supergranular flow structures in different fluids. (a) A stable stratification in the bulk is no omnipresent result of (or necessity for) the emergence of the supergranule but rather related to plume detachments. (b) Simultaneously, the supergranules at 
$k_{{h, min}} = 2 {\rm \pi}/ \varGamma \approx 0.1$ become, evaluating here the midplane, weaker relative to smaller-scale structures (such as granules) for decreasing 
${Pr}$. The temperature fields are re-scaled as described in the text, and both data correspond to the late state of the flows as described in the caption of table 1.
4. Discussion and perspective
Introducing buoyancy in a simple Rayleigh–Bénard convection configuration via a constant heat flux at the top and bottom planes leads without any additional physics to the emergence of a hierarchy of different long-living large-scale flow structures (Vieweg et al. Reference Vieweg, Scheel and Schumacher2021a; Vieweg Reference Vieweg2023). While this hierarchy consists of so-called granules and supergranules as separate stages, the latter are driven by secondary instabilities and might grow until the horizontal domain size is reached (see again § 1 and figure 1). The present study raises the question if mechanisms similar to this secondary instability in constant heat flux-driven Rayleigh–Bénard convection (Chapman et al. Reference Chapman, Childress and Proctor1980; Chapman & Proctor Reference Chapman and Proctor1980) act beyond 
${Pr} = 1$ (Vieweg et al. Reference Vieweg, Scheel and Schumacher2021a) independently of the working fluid and, in particular, even down to very small Prandtl numbers such as found in the solar convection zone (Rincon & Rieutord Reference Rincon and Rieutord2018; Schumacher & Sreenivasan Reference Schumacher and Sreenivasan2020). We therefore conducted a series of simulations across four orders of 
${Pr}$ given a fixed thermal driving in a horizontally extremely extended periodic domain of 
$\varGamma = 60$. We confirmed the presence of the gradual supergranule aggregation as a particular mechanism of self-organisation of long-living large-scale flow structures in naturally forced convection flows independently of the working fluid. Our observations thus suggest that these secondary instabilities dominate any basic heat flux-driven convection flow, leading to a robust hierarchy of different large-scale flow structures.
This omnipresent appearance is in accordance with the accessibility of large-scale 
$k_{z} = 0$ spectral modes in the temperature field for this particular thermal boundary condition. Note that these modes are not accessible in the classical case of applied constant temperatures (Vieweg et al. Reference Vieweg, Scheel, Stepanov and Schumacher2022). As the present configuration corresponds to a ratio of thermal diffusivities 
$\kappa _{{f}} / \kappa _{{s}} \rightarrow \infty$ between the fluid and the above or below solid, this mechanism can be seen as the result of a relaxation of thermal perturbations that happens much quicker in the fluid compared to in the solid plates (Hurle et al. Reference Hurle, Jakeman and Pike1967; Käufer et al. Reference Käufer, Vieweg, Schumacher and Cierpka2023). Moreover, the strength of buoyancy effects is in the heat flux-driven scenario not limited by prescribed temperatures at the boundaries but instead via only the less-restrictive mixing of fluid in between. Hence, these arguments allow and demand eventually the formation of horizontally extended flow structures that might even span across the entire domain to advectively transfer the thermal variance.
Given the fact that the variation of the working fluid affects the relative strength of thermal diffusion as described by (2.4a,b), one might expect a decrease of thermal variance for decreasing 
${Pr}$ due to an increase of 
$\kappa$. However, it turns out that the thermal variance increases: this is also indicated by the colour scales in figure 3. This observation can be explained as follows: smaller 
${Pr}$ result in larger 
${Re}$ and thus in an increased local mixing (with 
$\kappa$ acting also locally). As the flow is increasingly disordered, the large-scale supergranule becomes less dominant compared with smaller-scale velocity structures which is confirmed by the spectral analysis captured in figure 5(b). Consequently, the horizontal mixing on large scales (see also the previous paragraph) becomes successively less effective for smaller 
${Pr}$, leading eventually to an increased thermal variance in the horizontally extended domain. As becomes clear when contrasting the present results with Vieweg et al. (Reference Vieweg, Scheel and Schumacher2021a), (i) the vanishing stable stratification is not an effect of the increased Reynolds number and (ii) the relative heat transport of supergranules compared with smaller-scale structures (such as granules) loses similarly importance when increasing 
${Re}$ via 
${Ra}$ given a fixed 
${Pr}$.
Interestingly, despite the fundamentally different long-living large-scale flow structures between the cases of applied constant temperatures and vertical temperature gradients (see again § 1), their response on a variation of the working fluid shares clear analogies: compare therefore the bright and dark markers in figure 4, respectively. Note here that while the Rayleigh number 
${Ra}_{{D}} := \alpha g {\rm \Delta} T H^{3} / ( \nu \kappa )$ in case of an applied constant temperature difference 
${\rm \Delta} T$, this is related via 
${Ra}_{{D}} \equiv {Ra} / {Nu}$ (Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002; Foroozani, Krasnov & Schumacher Reference Foroozani, Krasnov and Schumacher2021; Vieweg et al. Reference Vieweg, Scheel and Schumacher2021a) to (2.4a,b). In particular, this allows to relate the present 
${Pr} = 1$ run to the corresponding no-slip and 
${Ra}_{{D}} = 10^{5}$ one from Pandey et al. (Reference Pandey, Krasnov, Sreenivasan and Schumacher2022) as described in Vieweg (Reference Vieweg2023). Thus, the particular kind of thermal boundary condition seems not to be of great significance when it comes to qualitative trends of the classical global measures of heat and momentum transport with respect to 
${Pr}$. In other words, different large-scale flow structures respond qualitatively similarly on a variation of the working fluid if judged via global measures of heat or momentum transport. Moreover, this underlines that diffusion processes are primarily locally important and do not rule the large-scale pattern formation.
The omnipresence of supergranule aggregation across all accessible Rayleigh and Prandtl numbers highlights the importance of an understanding of secondary (and subsequent) instabilities (Chapman et al. Reference Chapman, Childress and Proctor1980; Chapman & Proctor Reference Chapman and Proctor1980) slightly above the onset of convection. It is intriguing that such mechanisms survive even into the fully turbulent states of the flows (Vieweg et al. Reference Vieweg, Scheel and Schumacher2021a) where patterns are typically highly susceptible to the influence of instabilities and defects on each other (Busse Reference Busse1978, Reference Busse2003). Moreover, additional physical mechanisms are required to stop the gradual supergranule aggregation before reaching the numerically finite domain size. Weak rotation around the vertical axis has turned out to effectively interrupt this process in the turbulent regime (Vieweg et al. Reference Vieweg, Scheel, Stepanov and Schumacher2022) while also the primary instability changes qualitatively with 
$k_{{h,crit}} > 0$ once rotation is sufficiently strong (Dowling Reference Dowling1988; Takehiro et al. Reference Takehiro, Ishiwatari, Nakajima and Hayashi2002). Interestingly, the ratio of thermal diffusivities 
$\kappa _{{f}} / \kappa _{{s}}$ seems to promise similar effects (Hurle et al. Reference Hurle, Jakeman and Pike1967). This is of particular importance to better resemble the motivating geophysical and astrophysical flows and will be addressed in future studies.
Acknowledgements
The author thanks Professor J. Schumacher for valuable comments on the early manuscript.
Funding
The author is supported by the Deutsche Forschungsgemeinschaft within the Priority Programme DFG-SPP 1881 ‘Turbulent Superstructures’ as well as grant no. 1410/31-1. He gratefully acknowledges the Gauss Centre for Supercomuting e.V. (www.gauss-centre.eu) for funding this work by providing computing resources on the GCS supercomputer SuperMUC-NG at Leibnitz Supercomputing Centre within project pn68ni and through the John von Neumann Institute for Computing (NIC) on the GCS supercomputer JUWELS at Jülich Supercomputing Center (JSC) within projects mesoc and nonbou. In addition, he acknowledges the computing centre of the Technische Universität Ilmenau for providing access to, as well as computing and storage resources on its compute cluster MaPaCC4.
Declaration of interests
The author reports no conflict of interest.
