Supergranule aggregation: a Prandtl number-independent feature of constant heat flux-driven convection flows

Abstract

Supergranule aggregation, i.e. the gradual aggregation of convection cells to horizontally extended networks of flow structures, is a unique feature of constant heat flux-driven turbulent convection. In the present study, we address the question if this mechanism of self-organisation of the flow is present for any fluid. Therefore, we analyse three-dimensional Rayleigh–Bénard convection at a fixed Rayleigh number ${Ra} \approx 2.0 \times 10^{5}$ across $4$ orders of Prandtl numbers ${Pr} \in [10^{-2}, 10^{2}]$ by means of direct numerical simulations in horizontally extended periodic domains with aspect ratio $\varGamma = 60$. Our study confirms the omnipresence of the mechanism of supergranule aggregation for the entire range of investigated fluids. Moreover, we analyse the effect of ${Pr}$ on the global heat and momentum transport, and clarify the role of a potential stable stratification in the bulk of the fluid layer. The ubiquity of the investigated mechanism of flow self-organisation underlines its relevance for pattern formation in geophysical and astrophysical convection flows, the latter of which are often driven by prescribed heat fluxes.

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 1995) over deep ocean convection (Maxworthy & Narimousa 1994) up to convection in its atmosphere (Atkinson & Wu Zhang 1996), 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 1993). 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 2020). 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 2018; Schumacher & Sreenivasan 2020), they offer very different lifetimes and horizontal extensions. Unfortunately, our understanding of such hierarchies’ origins is still far from complete (Hanson et al. 2020) 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 2023).

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 2018; Stevens et al. 2018; Krug, Lohse & Stevens 2020; Käufer et al. 2023), 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 2021a; Vieweg et al. 2022). 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. 2021b; Schneide et al. 2022; Vieweg et al. 2024). 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. 2021a, 2022). 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 1916; Pellew & Southwell 1940) or $k_{{h,crit}} = 0$ (Hurle et al. 1967) 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 1980). 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 1980; Vieweg et al. 2021a), 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 1978, 2003). 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 1980) but simultaneously ‘do not expect the theory to remain accurate for very small ${Pr}$’ (Chapman & Proctor 1980). 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 1967). 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 1879; Boussinesq 1903) 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 1916). 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 1997; Scheel, Emran & Schumacher 2013). 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

(2.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

(2.2)\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}

(2.3)\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,

(2.4a,b)\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

(2.5a,b)\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

(2.6)\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. 2013). These dissipation scales are given by the so-called Kolmogorov and Batchelor scale (Batchelor 1959; Kolmogorov 1991; Sreenivasan 2004),

(2.7a,b)\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 1951) 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. 2021a), 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 2023; Vieweg et al. 2024), 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. 2004) of the temperature field

(3.1)\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. (2022). 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. 2021a, 2022). Running our simulations reveals two particularly interesting results.

First, the gradual supergranule aggregation, first reported by Vieweg et al. (2021a), 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. 2022) 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. 1980; Chapman & Proctor 1980) 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 1991) and Batchelor (Batchelor 1959) 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

(3.2)\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. 2002)

(3.3)\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 2012) (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. 2022) 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 2023).

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. 2021a) 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. 2021a; Vieweg 2023). 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. 1980; Chapman & Proctor 1980) act beyond ${Pr} = 1$ (Vieweg et al. 2021a) 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 2018; Schumacher & Sreenivasan 2020). 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. 2022). 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. 1967; Käufer et al. 2023). 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. (2021a), (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. 2002; Foroozani, Krasnov & Schumacher 2021; Vieweg et al. 2021a) 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. (2022) as described in Vieweg (2023). 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. 1980; Chapman & Proctor 1980) 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. 2021a) where patterns are typically highly susceptible to the influence of instabilities and defects on each other (Busse 1978, 2003). 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. 2022) while also the primary instability changes qualitatively with $k_{{h,crit}} > 0$ once rotation is sufficiently strong (Dowling 1988; Takehiro et al. 2002). Interestingly, the ratio of thermal diffusivities $\kappa _{{f}} / \kappa _{{s}}$ seems to promise similar effects (Hurle et al. 1967). This is of particular importance to better resemble the motivating geophysical and astrophysical flows and will be addressed in future studies.