Convective mesoscale turbulence at very low Prandtl numbers

Abstract

Horizontally extended turbulent convection, termed mesoscale convection in natural systems, remains a challenge to investigate in both experiments and simulations. This is particularly so for very low molecular Prandtl numbers, such as occur in stellar convection and in the Earth's outer core. The present study reports three-dimensional direct numerical simulations of turbulent Rayleigh–Bénard convection in square boxes of side length $L$ and height $H$ with the aspect ratio $\varGamma =L/H$ of 25, for Prandtl numbers that span almost 4 orders of magnitude, $10^{-3}\le Pr \le 7$, and Rayleigh numbers $10^5 \le Ra \le 10^7$, obtained by massively parallel computations on grids of up to $5.36\times 10^{11}$ points. The low end of this $Pr$-range cannot be accessed in controlled laboratory measurements. We report the essential properties of the flow and their trends with the Rayleigh and Prandtl numbers, in particular, the global transport of momentum and heat – the latter decomposed into convective and diffusive contributions – across the convection layer, mean vertical profiles of the temperature and temperature fluctuations and the kinetic energy and thermal dissipation rates. We also explore the degree to which the turbulence in the bulk of the convection layer resembles classical homogeneous and isotropic turbulence in terms of spectra, increment moments and dissipative anomaly, and find close similarities. Finally, we show that a characteristic scale of the order of the mesoscale seems to saturate to a wavelength of $\lambda \gtrsim 3H$ for $Pr\lesssim 0.005$. We briefly discuss possible implications of these results for the development of subgrid-scale parameterization of turbulent convection.

1. Introduction

Thermal convection in stellar and planetary interiors and their atmospheres is complex because it is driven by several factors in combination, but researchers often study the idealized model of convection – Rayleigh–Bénard convection (RBC) – where a fluid layer bounded between two horizontal plates is heated and cooled uniformly from the bottom and the top, respectively (Tritton 1977; Siggia 1994; Kadanoff 2001; Ahlers, Grossmann & Lohse 2009; Chillà & Schumacher 2012). An important parameter governing RBC is the Prandtl number $Pr$, which is the ratio of the kinematic viscosity $\nu$ and the thermal diffusivity $\kappa$ of the fluid. Stellar and planetary convection is often characterized by very small Prandtl numbers; for example, $Pr \approx 10^{-6}$ in the Sun (Brandenburg & Subramanian 2005; Schumacher & Sreenivasan 2020; Garaud 2021) and $Pr \approx 10^{-2}\unicode{x2013}10^{-1}$ in the Earth's outer core (Calkins et al. 2012; Aurnou et al. 2015; Guervilly, Cardin & Schaeffer 2019). The Rayleigh number $Ra$, which characterizes the strength of the driving buoyancy force relative to viscous and thermal dissipative forces, is another control parameter, and is very high in most natural flows. The ratio of the horizontal and vertical extents of a convective flow is the aspect ratio $\varGamma$; a crucial feature of all natural convective flows is that $\varGamma \gg 1$, which allows the formation of turbulent superstructures – coherent flow patterns with characteristic scale larger than the depth of the convection layer (Cattaneo, Lenz & Weiss 2001; Rincon, Lignières & Rieutord 2005). A possible example is supergranulation on the Sun's surface (Nordlund, Stein & Asplund 2009; Rincon & Rieutord 2018). Even though RBC incorporates a number of simplifications, such as the Boussinesq approximation (Tritton 1977; Schumacher & Sreenivasan 2020), it appears that a study of the flow in extended layers at low Prandtl numbers is highly worthwhile. This is the primary objective of the current work. We are also motivated by the relevance of low-$Pr$ convection for industrial applications that use liquid metals.

Despite this relevance, low-$Pr$ turbulent convection has not been explored extensively in experiments mainly because liquid metals such as mercury, gallium and sodium are difficult to handle and optically opaque (Cioni, Ciliberto & Sommeria 1997; Zürner et al. 2019). Even when these difficulties are circumvented, the lowest Prandtl number that can be explored is $Pr \approx 0.006$ for liquid sodium (Horanyi, Krebs & Müller 1999), which is still some three orders of magnitude higher than in the solar convection zone. Direct numerical simulations (DNS) offer important tools but they, too, are hindered by demanding resolution requirements due to the highly inertial nature of low-$Pr$ convection (Breuer et al. 2004; Schumacher, Götzfried & Scheel 2015; Scheel & Schumacher 2016, 2017; Pandey, Scheel & Schumacher 2018; Zwirner et al. 2020). The fact that the required computational power increases with increasing aspect ratio as $\varGamma ^2$ further limits numerical investigations. Yet, Pandey et al. (2018) were able to perform DNS of RBC in a $\varGamma = 25$ domain by achieving $Pr$ as low as 0.005 at $Ra = 10^5$ to study turbulent superstructures. In the present work, we significantly extend the parameter range from Pandey et al. (2018) and Fonda et al. (2019) by further decreasing the Prandtl number fivefold, while also increasing the Rayleigh number by two orders of magnitude. The highest Reynolds number achieved in the present work is nearly $5.6 \times 10^4$, requiring massively parallel DNS on computational grids of more than $5\times 10^{11}$ points.

This work has three main goals: (i) report trends of heat and momentum transfer with respect to Prandtl number over nearly 4 orders of magnitude; (ii) assess the closeness of small-scale statistical properties in the bulk of the convection layer to the classical Kolmogorov-type behaviour (Kolmogorov 1941b; Frisch 1995). This assessment comprises energy spectra, a test of the 4/5th law and an investigation of the dissipative anomaly (Sreenivasan 1984, 1998) in turbulent convection flow with boundaries; (iii) analyse the large-scale circulation patterns, the turbulent superstructures of convection, particularly their trends with decreasing $Pr$. Note that the convective flows in the Earth's and stellar interiors are also associated with (differential) rotation and magnetic fields that can lead to strong departures from local isotropy at larger and intermediate scales (Aurnou et al. 2015), but we shall here focus on the influence of low Prandtl number. These DNS series will thus provide a unique data base for the parametrization of turbulent transport in mesoscale configurations characterized by a degree of large-scale order, with well-resolved thermal and kinetic energy dissipation rates.

It is becoming increasingly clear that many properties of low-$Pr$ convective flows differ from those at moderate and high Prandtl numbers. For example, the efficacy of low-$Pr$ flows in transporting heat is lower, and that in transporting momentum higher, than in high-$Pr$ flows; the disparity between the two increases as $Pr$ is lowered (Scheel & Schumacher 2017). Due to high (low) thermal (momentum) diffusivity, low-$Pr$ convection exhibits coarser thermal structures but the length scales in the velocity field have a broader distribution. This results in an enhanced separation between the energy injection and energy dissipation scales in low-$Pr$ convective flows (Schumacher et al. 2015). In this regime, the kinetic energy spectrum has been observed to approximate the classical Kolmogorov scaling (Kolmogorov 1941b) with the $-5/3$ power in the inertial range (Lohse & Xia 2010; Mishra & Verma 2010; Bhattacharya, Verma & Samtaney 2021). Here, we analyse the kinetic energy spectra in the bulk region of the flow and show that they indeed exhibit the classical Kolmogorov scaling with an inertial range, showing no tendency towards Bolgiano scaling (Bolgiano 1959), according to which, the conversion of the kinetic to potential energy leads to the steeper $k^{-11/5}$ scaling with the wavenumber k in the inertial range (Lohse & Xia 2010; Verma, Kumar & Pandey 2017; Verma 2018).

Turbulent superstructures of convection can be characterized by a typical spatial scale $\lambda$ and a temporal scale $\tau$; finer scales evolve much faster than $\tau$. Thus, the scales $\lambda$ and $\tau$ help distinguish the coarse and gradually evolving large-scale patterns from the finer (and faster) turbulent fluctuations (Pandey et al. 2018; Krug, Lohse & Stevens 2020). The characteristic scales of superstructures have been observed to depend on $Pr$ and $Ra$ (Hartlep, Tilgner & Busse 2003, 2005; von Hardenberg et al. 2008; Bailon-Cuba, Emran & Schumacher 2010; Emran & Schumacher 2015; Pandey et al. 2018; Schneide et al. 2018; Stevens et al. 2018; Fonda et al. 2019; Green et al. 2020; Krug et al. 2020; Lenzi, von Hardenberg & Provenzale 2021; Pandey, Schumacher & Sreenivasan 2021), as well as on the thermal boundary conditions at the horizontal top and bottom plates (Vieweg, Scheel & Schumacher 2021a). The characteristic length scale of superstructures, which is nearly twice the depth $H$ of the convection layer at the onset of convection, increases with increasing $Ra$ (Pandey et al. 2018; Stevens et al. 2018). This dependence on $Pr$ is complex, with $\lambda (Pr)$ showing a peak near $Pr \approx 7$ and decreasing as $Pr$ departs from this value (Pandey et al. 2018). This decreasing trend of $\lambda (Pr)$ continues to hold up to a $Pr \approx 0.005$ below which the scales seem to level off at a wavelength of $\lambda \gtrsim 3H$.

The remainder of this article is organized as follows. In § 2, we briefly describe the DNS and note the parameter space explored. In § 3, we discuss the flow structures and the scaling of the global transport of heat and momentum, and study in § 4 the vertical profiles of temperature, convective and diffusive transports, as well as dissipation rates. Statistical properties of the flow such as kinetic energy spectra and third-order structure functions are examined in § 5 for their compatibility with Kolmogorov forms. The characterization of turbulent superstructures is presented in § 6. We conclude the main findings in § 7, and present an outlook. Appendices A and B deal with specific tests of sufficient resolution. Appendix C discusses technical detail of the estimates of length and velocity scales of superstructures.

2. Details of DNS

We perform DNS of RBC in a closed rectangular domain with square cross-section of length $L = 25H$, where $H$ is the depth of the convection layer. We solve the following non-dimensionalized equations incorporating the Oberbeck–Boussinesq (OB) approximation:

(2.1)$$\begin{gather} \frac{\partial {\boldsymbol u}}{\partial t} + ({\boldsymbol u} \boldsymbol{\cdot} {\boldsymbol \nabla}) {\boldsymbol u} ={-}{\boldsymbol \nabla}p + T \hat{{\boldsymbol z}} + \sqrt{\frac{Pr}{Ra}} \, \nabla^2 {\boldsymbol u}, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial T}{\partial t} + ({\boldsymbol u} \boldsymbol{\cdot} {\boldsymbol \nabla}) T = \frac{1}{\sqrt{PrRa}} \, \nabla^2 T, \end{gather}$$

(2.3)$$\begin{gather}{\boldsymbol \nabla} \boldsymbol{\cdot} {\boldsymbol u} = 0. \end{gather}$$

Here, ${\boldsymbol u} \equiv (u_x,u_y,u_z)$, $T$ and $p$ are the velocity, temperature and pressure fields, respectively. The Prandtl number $Pr$ and the Rayleigh number $Ra = \alpha g \Delta T H^3/\nu \kappa$, where $\alpha$ is the coefficient of thermal expansion of the fluid, $g$ is the acceleration due to gravity and $\Delta T$ is the imposed temperature difference between the horizontal plates. We have used the layer's depth $H$, the free-fall velocity $u_f = \sqrt {\alpha g \Delta T H}$, the free-fall time $t_f = H/u_f$ and $\Delta T$ as the non-dimensionalizing length, velocity, time and temperature scales, respectively. We use the no-slip condition on all boundaries. We employ the isothermal condition on the horizontal plates and the adiabatic condition on the sidewalls. This allows direct comparison with other simulations at higher $Pr$ (Pandey et al. 2018; Fonda et al. 2019) and controlled laboratory experiments, such as those of Moller et al. (2022), in exactly the same setting to enable a consistent analysis across the $Ra$–$Pr$ parameter plane.

We use two different solvers for the simulations. For moderate and large $Pr$, we use a spectral element solver Nek5000 (Fischer 1997), where the flow domain is divided into a finite number of elements $N_e$. The Lagrangian interpolation polynomials of order $N$ are further used to expand the turbulence fields within each element (Scheel, Emran & Schumacher 2013), thus resulting in a total of $N_e N^3$ mesh cells in the entire flow domain. As low-$Pr$ convection is dominated by inertial forces, the resulting flow acquires increased fine structure whose resolution requires more extensive computational resources (Schumacher et al. 2015). Therefore, we performed those simulations using a second-order finite difference solver that requires significantly less working memory at a given grid size. Here, the flow domain is divided into $N_x \times N_y \times N_z$ non-uniform mesh cells (Krasnov, Zikanov & Boeck 2011; Liu, Krasnov & Schumacher 2018). We have verified that the results obtained from both the solvers agree well with each other by performing two simulations for $Pr = 0.005, Ra = 10^5$ and $Pr = 0.7, Ra = 10^7$ using both solvers. We refer to Appendix A for a direct comparison of globally averaged and horizontally averaged convective heat fluxes and dissipation rates from the two solvers. Important parameters of all the simulations are provided in table 1.

Table 1. Important parameters of the simulations in a rectangular box of $\varGamma = 25$ with square cross-section; the numbers of mesh cells, $N_x \times N_y \times N_z$ and $N_e \times N^3$, are for the finite difference and spectral element solvers, respectively; $N_{tot}$ represents the total number of mesh cells in units of a billion; $Nu$ and $Re$ are the volume and time averaged Nusselt and Reynolds numbers, respectively; $u_{rms}$ is the root-mean-square velocity computed over the entire volume; $\lambda$ and $\tau$ are, respectively, the characteristic length and time scales of turbulent superstructures. Runs with superscript $a$ are taken from Pandey et al. (2018), while those with superscript $b$ are taken from Fonda et al. (2019).

3. Flow morphology and global transport

3.1. Structures of velocity and temperature fields

Because the time scales of heat and momentum diffusion processes are very different in low-$Pr$ convection, the temperature field shows coarser structures than the velocity field. This is illustrated in figure 1, which displays the instantaneous temperature, vertical velocity and local turbulent kinetic energy fields in the mid-horizontal plane, $z = H/2$, for the biggest simulations with $Pr = 0.001, Ra = 10^7$. Panels (a,d,g) show the fields in the entire cross-section, and (b,c,e,f,h,i) depict the marked magnifications to highlight small-scale structures. The flow pattern of the temperature and vertical velocity fields show similarities at large scales, but the velocity field also consists of very fine structures compared with the highly diffusive temperature field. The finest scale of the turbulent velocity field, denoted as the Kolmogorov scale $\eta$, is estimated as $\eta = (\nu ^3/\langle \varepsilon _u\rangle _{V,t})^{1/4}$. Here, $\langle \varepsilon _u\rangle _{V,t}$ is the combined volume–time average of the kinetic energy dissipation rate field per unit mass, computed at each point by

(3.1)\begin{equation} \varepsilon_u({\boldsymbol x},t) = \frac{\nu}{2} \sum_{i,j=1}^3\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)^2, \end{equation}

with $u_i$ representing the velocity component in the direction of coordinate $x_i$. The finest scale of the temperature field is either the Corrsin scale $\eta _C=\eta /Pr^{3/4}$, which marks the end of the inertial–convective range for $Pr<1$, or the Batchelor scale $\eta _B=\eta /Pr^{1/2}$, which marks the end of the viscous–convective range for $Pr>1$; see e.g. Sreenivasan & Schumacher (2010). It is clear that $\eta < \eta _C$ when $Pr < 1$, and $\eta _B < \eta$ when $Pr>1$. Thus, the finest scales in the flow at hand are either $\eta$ for $Pr<1$ or $\eta _B$ for $Pr>1$.

Figure 1. Turbulent superstructures of convection in a low-$Pr$ flow with $Pr = 0.001$ and $Ra = 10^7$. The panels represent instantaneous temperature fields (a–c), vertical velocity (d–f) and turbulent kinetic energy (g–i) in the midplane. In this low-$Pr$ flow, the thermal energy is primarily contained in large-scales, whereas the kinetic energy is distributed over a broad range of scales. Panels (a,d,g), (b,e,h) and (c, f,i) represent fields of view that are $25H \times 25H$, $6.25H \times 6.25H$ and $1.56H \times 1.56H$, respectively.

The Corrsin scale is nearly 178 times larger than the Kolmogorov scale for $Pr = 0.001$. This large difference is clearly visible in figure 1(c, f), where the temperature and vertical velocity fields are shown in a small cross-section of size $1.56H \times 1.56H$. The figures reveal that the smallest length scale of the thermal structures – the length scale over which the temperature variation is significant – is of the order of $H$, whereas that for velocity structures is much finer. We also find that the dominant structures in the kinetic energy field resemble those in the vertical velocity because of its dominance in the midplane (Pandey et al. 2018). This wide range of length scales present in low-$Pr$ convection engenders a broad inertial range in kinetic energy spectrum, which will be discussed in § 5.2.

To see the effects of $Pr$ on flow structures, we show the temperature and the vertical velocity fields for $Pr = 0.001$, $0.7$ and $7$ at $Ra = 10^7$ in figure 2. To accentuate small structures, the fields are shown in a quarter (in linear dimension) of the entire cross-section. With increasing $Pr$, increasingly finer thermal structures are generated due to decreasing thermal diffusivity. On the other hand, the velocity variation becomes progressively regular as the viscosity increases (or the Reynolds number decreases) with increasing $Pr$.

Figure 2. Temperature (a–c) and vertical velocity (d–f) in midplane for $Ra = 10^7$. Panels (a,d) are for $Pr = 0.001$, (b,e) for $Pr = 0.7$ and (c, f) for $Pr = 7$. Fields are displayed in a magnified region of dimensions $6.25 H \times 6.25 H$ around the centre. Finer temperature contours can be observed for higher $Pr$, panel(c), compared with lower $Pr$, panel (a). The velocity field for lower $Pr$, panel (d), exhibits much finer structures than those in the other two panels, because the Reynolds number decreases as $Pr$ increases.

3.2. Heat and momentum transport laws

Convection at low Prandtl numbers differs from its high-$Pr$ counterpart by reduced heat transport and enhanced momentum transport (Schumacher et al. 2015; Scheel & Schumacher 2016, 2017; Pandey et al. 2018; Zürner et al. 2019; Zwirner et al. 2020). Heat transport is quantified by the Nusselt number $Nu$, defined as the ratio of the total heat transport to that by conduction alone. It is computed as

(3.2)\begin{equation} Nu = 1 + \sqrt{RaPr} \, \langle u_z T \rangle_{V,t} , \end{equation}

where $\langle \,\cdot \, \rangle _{V,t}$ denotes again the average over the entire simulation domain and time. We compute $Nu$ for all simulations and plot them, for fixed $Ra$, as a function of $Pr$ in figure 3(a); $Nu$ increases up to $Pr = 0.7$ but does not change significantly thereafter. A similar trend has also been reported in literature for convection in $\varGamma \approx 1$ domains (Verzicco & Camussi 1999; Schmalzl, Breuer & Hansen 2004; van der Poel, Stevens & Lohse 2013) and also for $\varGamma = 0.1$ (Pandey & Sreenivasan 2021). Figure 3(a) indicates that the molecular diffusion becomes an increasingly dominant mode of heat transport as $Pr$ decreases. For $Ra = 10^5$ and $10^6$, we do best fits to the data for $Pr \leq 0.7$. The transport laws $Nu(Pr)$ for both Rayleigh numbers are given, including error bars, in table 2. In summary, we find that the Nusselt number is consistent with the power law $Nu\sim Pr^{0.19}$ for $Ra = 10^5$ and $Nu\sim Pr^{0.18}$ for $Ra=10^6$. These power-law exponents in the extended convection domain are within the range observed in RBC with $\varGamma \lesssim 1$, as shown in Pandey & Sreenivasan (2021), where more discussion of the $Nu-Pr$ scaling exponent can be found.

Figure 3. (a) Nusselt number for $\varGamma = 25$ increases with increasing $Pr$ for low Prandtl numbers but does not change much for $Pr \geq 0.7$, consistent with Pandey & Sreenivasan (2021). (b) The value of $Nu$ as a function of $Ra$ increases approximately as $Ra^{0.29}$. (c) The Reynolds number based on the root-mean-square velocity decreases with increasing $Pr$, with power law exponents consistent with others in the literature. (d) The value of $Re$ as a function of $Ra$ increases as $Ra^{0.5}$ for all the three cases. The dashed lines in panels (b,d) are best fits, summarized in table 2. Solid lines in panels (a,c) are predictions of $Nu$ and $Re$ from the Grossmann–Lohse theory. Legends in panels (a,b) also apply to the corresponding panels (c,d).

Table 2. Summary of scaling relations for global heat and momentum transports as functions of $Ra$ and $Pr$. Note that the scaling laws with respect to $Ra$ have been obtained by fits to three data points only; it is clear that more definitive results require larger number of data points.

Normalized values of global heat transport in RBC increase with increasing thermal driving and the rate of increase depends on the Prandtl number (Scheel & Schumacher 2016, 2017). We plot $Nu$ for $Pr = 0.001, 0.7$ and 7 against $Ra$ in figure 3(b), which shows that $Nu$ values for $Pr = 0.7$ and 7 are similar. There are only three data points, which would not be adequate to establish a new result. Nevertheless, power-law fits to those three points serve to supplement existing results. We obtain approximately $Nu \sim Ra^{0.29}$ (see table 2). The exponents for all the three Prandtl numbers are essentially similar and agree with those observed in convection for $\varGamma \sim 1$ (Bailon-Cuba et al. 2010; Stevens, Lohse & Verzicco 2011; Scheel, Kim & White 2012; Scheel & Schumacher 2014, 2016). It is interesting that the exponent for $Pr = 0.001$ is not lower compared with that for $Pr \geq 0.7$. A slightly lower scaling exponent of $0.27\pm 0.01$ was reported from simulations in closed cylinders for $\varGamma =1$ (Scheel & Schumacher 2017) at $Pr=0.021$. Recent experiments in strongly turbulent liquid metal convection by Schindler et al. (2022), at nearly the same Prandtl number but for Rayleigh numbers up to $Ra=5\times 10^9$, in a cylinder with $\varGamma =1/2$, reported an even smaller scaling exponent of 0.124. It is possible that the constrained large-scale flow in a closed cylinder affects the scaling exponent at low and moderate Rayleigh numbers, see discussion by Pandey & Sreenivasan (2021).

The Reynolds number $Re$ quantifies the momentum transport in RBC. We compute it with $H$ and $u_{rms}$ as the relevant length and velocity scales, as

(3.3)\begin{equation} Re = u_{rms} \, \sqrt{\frac{Ra}{Pr}} \quad \mbox{where}\ u_{rms} = \sqrt{\langle u_i^2 \rangle_{V,t}} \end{equation}

is the root-mean-square (r.m.s.) velocity. The Reynolds number as a function of $Pr$ is plotted in figure 3(c), which reveals that, for fixed $Ra$, the flow loses its effectiveness in transporting momentum as $Pr$ increases (Käpylä 2021). In thermal convection, the power-law exponent of $Re-Pr$ scaling depends on the range of $Pr$; the Reynolds number decreases with increasing $Pr$ even when the flow is dominated by inertia. The detailed power-law fits can be found in table 2. As a summary, we get $Re \sim Pr^{-0.62}$ and $Re \sim Pr^{-0.65}$ for $Ra = 10^5$ and $Ra = 10^6$, respectively. As for the Nusselt number, the exponents of the $Re-Pr$ scaling agree with those observed for $\varGamma \lesssim 1$ (Verzicco & Camussi 1999; Pandey & Sreenivasan 2021).

The Reynolds number variation with $Ra$, plotted in figure 3(d), shows that $Re$ is consistently higher for lower Prandtl numbers, manifesting in the enhanced prefactors of the power-law fits which are summarized in table 2. The best fits yield $Re \sim Ra^{0.53}$, ${Re \sim Ra^{0.49}}$ and $Re \sim Ra^{0.54}$ for $Pr = 0.001, 0.7$ and 7, respectively. Note that the Reynolds number based on the free-fall velocity scales as $Ra^{0.50}$. Thus, these scaling exponents suggest that the free-fall velocity for a fixed $Pr$ does not depend strongly on $Ra$ (see table 1). The scaling exponents are in the same range as in several other studies in the past; the exponent does not, however, decrease with decreasing Prandtl number as found by Scheel & Schumacher (2017). This might be the result of differences in the aspect ratio, but we reiterate that the present fits use only three data points.

Attempts have been made to predict the global transports in RBC as a function of the control parameters (Shraiman & Siggia 1990; Grossmann & Lohse 2000; Pandey & Verma 2016). Grossmann & Lohse (2000) assumed the existence of a large-scale circulation of the order of the size of the convection cell, and proposed a set of coupled equations relating $Nu$ and $Re$ as functions of $Ra$ and $Pr$ (Grossmann & Lohse 2001). The equations also include a set of constant coefficients, whose values depend on the aspect ratio of the domain. Using the coefficients provided in Stevens et al. (2013) for $\varGamma \approx 1$ RBC, we compute $Nu$ as a function of $Pr$ from the Grossmann–Lohse model and show them as solid curves in figure 3(a). The Nusselt numbers thus estimated are somewhat higher than those computed from the DNS for $Ra = 10^5$. On the low-$Pr$ end, this may be attributed to the fact that the temperature fields for these parameter pairs are dominated by diffusion and barely mixed in the bulk. However, the agreement is better for $Ra \geq 10^6$, which indicates that the heat transport in our extended cell is not much different from that in $\varGamma = 1$ cells. We also plot $Re(Pr)$ from the Grossmann–Lohse model in figure 3(c), and find that there is fair agreement; see also Verma (2018).

4. Vertical profiles across the convection layer

4.1. Temperature and heat flux fields

In the conductive equilibrium state, the vertical temperature gradient is a constant; inhomogeneities in the horizontal directions arise in the convective state, leading to a modification of the linear temperature profile. We compute the mean temperature profiles $\langle T \rangle _{A,t}(z)$ and plot them in figure 4. Here, $\langle \,\cdot \, \rangle _{A,t}$ stands for the averaging over the entire horizontal cross-section of $A=25H\times 25H$ at a fixed height $z$ and the full time interval. In a turbulent convective flow, almost the entire temperature drop occurs within the thermal boundary layers (BLs) on the horizontal plates, while the bulk of the flow outside these BLs remains nearly isothermal (and thus well mixed). Figure 4(a) exhibits this feature. However, the slope of the temperature profile in the midplane increases as $Pr$ decreases. We plot the profiles for $Pr = 0.001$ for all the Rayleigh numbers in figure 4(b). The profile for $Ra = 10^5$ departs only weakly from the linear conduction profile despite a high Reynolds number of the flow. The temperature gradient in the central plane decreases with increasing $Ra$, and even a Rayleigh number of $10^7$ is not enough to generate a well-mixed temperature field in the bulk region for this very low $Pr$.

Figure 4. Horizontal and time averages of temperature as a function of the depth for simulations at (a) $Ra = 10^6$ and (b) $Pr = 0.001$. A well-mixed isothermal region away from the walls occurs only for $Pr \geq 0.7$ in (a), whereas a significant temperature gradient in the central region occurs for lower $Pr$. Dashed black line in panel (b) corresponds to the dimensionless conduction temperature profile $T_{cond} = 1 - z$.

In OB convection, the temperature averaged over the entire flow domain is $\Delta T/2$ but fluctuates at each point in the flow. We decompose the temperature field into its mean and fluctuation as

(4.1)\begin{equation} T({\boldsymbol x},t) = \langle T \rangle_{A,t}(z) + \theta({\boldsymbol x},t) . \end{equation}

Even though the temperature field becomes increasingly diffusive as $Pr$ decreases, the fluctuations increase with decreasing $Pr$; see figure 5(a) for $Ra = 10^7$. The depth variation is captured by the planar temperature fluctuation computed as

(4.2)\begin{equation} \theta_{rms}(z) = \sqrt{ \langle [T - \langle T \rangle_{A,t}(z)]^2 \rangle_{A,t} } = \sqrt{ \langle T^2(z) \rangle_{A,t} - \langle T(z) \rangle_{A,t}^2}. \end{equation}

Figure 5. (a) The r.m.s. temperature fluctuation profiles averaged over the top and bottom halves varying with the distance from the plate for $Ra = 10^7$. The peaks in $\theta _{rms}(z)$ occur near the thermal BL edges, which are indicated by dashed vertical lines. (b) Vertical profiles of the convective heat flux for $Pr = 0.001$. Dashed horizontal lines indicate the global heat flux ($Nu$) for each case. The convective flux vanishes at the top and bottom plates and is largest in the central plane.

Figure 5(a) shows that $\theta _{rms}(z)$ vanishes at the plate due to the imposed isothermal boundary condition. With increasing distance from the bottom plate, however, the strength of fluctuations increases within the thermal BL region. The maxima in $\theta _{rms}(z)$ profiles occur near the edge of the thermal BL (computed as $0.5H/Nu$), marked as dashed vertical lines in figure 5(a). This suggests that the thermal plumes retain their temperature, while the temperature of the ambient fluid decreases (increases) with increasing distance from the bottom (top) plate. This leads to an increasing contrast between the two components of the flow and is reflected as an increasing $\theta _{rms}(z)$ within the BL region (Pandey 2021). In the bulk region, however, $\theta _{rms}$ decreases with distance from the plate because the plumes do not retain their identity and begin to mix with the bulk fluid.

The heat transport occurs due to convective as well as diffusive processes, with their ratio varying with depth. To get the total heat flux in a horizontal plane, we average the temperature equation (2.2) in horizontal directions and in time, which leads to

(4.3)\begin{equation} Nu(z) = \sqrt{RaPr} \langle u_z T \rangle_{A,t} - \frac{\partial \langle T \rangle_{A,t}}{\partial z}=\text{const.} \end{equation}

It is clear from the temperature profiles in figure 4 that the diffusive contribution $- \partial \langle T \rangle _{A,t}/\partial z$ should be small in the well-mixed bulk region – increasing towards the plates and becoming largest at the plates. The variation of the convective heat flux $\sqrt {RaPr} \langle u_z T \rangle _{A,t}$ with depth in figure 5(b) for $Pr = 0.001$ confirms this expectation. The magnitudes of the globally averaged heat flux are indicated as dashed horizontal lines in figure 5(b), showing that the diffusive flux (the distance between the solid curves and the corresponding dashed horizontal lines) is not negligible even in the central region for $Ra \leq 10^6$. The diffusive component dominates the total heat flux in the central region for $Ra = 10^5$, which is consistent with the highly inefficient convective heat transport; see table 1. However, for $Ra = 10^7$, the diffusive contribution diminishes in the central plane. Thus, as $Pr$ becomes smaller, one requires increasing $Ra$ before turbulent processes become important.

4.2. Thermal and kinetic energy dissipation rates

While the mean temperature $\langle T\rangle _{A,t}(z)$ varies sharply near the horizontal plates and weakly in the central region, the vertical mean profile of the thermal dissipation rate field, which is the rate of loss of thermal variance that is computed pointwise by

(4.4)\begin{equation} \varepsilon_T({\boldsymbol x},t) = \kappa \left[ \left( \frac{\partial T}{\partial x} \right)^2 + \left( \frac{\partial T}{\partial y} \right)^2 + \left( \frac{\partial T}{\partial z} \right)^2 \right], \end{equation}

is higher in the vicinity of the horizontal plates and decreases towards the centre (Scheel & Schumacher 2016). We also compute the thermal dissipation rate field defined as

(4.5)\begin{equation} \varepsilon_\theta({\boldsymbol x},t) = \kappa \left[ \left( \frac{\partial \theta}{\partial x} \right)^2 + \left( \frac{\partial \theta}{\partial y} \right)^2 + \left( \frac{\partial \theta}{\partial z} \right)^2 \right] \end{equation}

to quantify the spatial variation of the temperature fluctuations. The mean profile of the thermal dissipation rate $\langle \varepsilon _\theta \rangle _{A,t}(z)$ is plotted in figure 6(a) for $Ra = 10^7$. Note that the vertical mean profiles of $\varepsilon _T$ and $\varepsilon _\theta$ are related by

(4.6)\begin{equation} \langle \varepsilon_T \rangle_{A,t}(z) = \langle \varepsilon_\theta \rangle_{A,t}(z) + \varepsilon_{\langle T \rangle}(z), \end{equation}

where $\varepsilon _{\langle T \rangle } = \kappa \,({\rm d} \langle T \rangle _{A,t}/{\rm d} z)^2$ is the dissipation rate corresponding to the mean temperature profile (Emran & Schumacher 2008). In convective flows with well-developed thermal BLs, $\varepsilon _{\langle T \rangle }$ contributes primarily to the BLs and negligibly in the bulk. This rapid decrease of $\varepsilon _{\langle T \rangle }$ outside the thermal BL region shows a shallow kink in the profiles of $\langle \varepsilon _\theta \rangle _{A,t}(z)$. In figure 6(a), we indicate the thermal BL thicknesses for $Pr = 0.7$ and $Pr = 7$ as dashed vertical lines, and note that the kinks are observed near the edge of the thermal BL. The kink does not appear for $Pr = 0.001$ due to the absence of well-developed thermal BLs.

Figure 6. Variation of the horizontally averaged (a) thermal and (b) kinetic energy dissipation rates in the vertical direction for $Ra = 10^7$. The profiles are further averaged over the top and bottom halves of the domain to improve the statistics. The dissipation profiles are largest at the plates and decrease towards the central plane; however, $\varepsilon _u(z)$ is nearly uniform in the bulk region. The (indistinguishable) dashed vertical lines in panel (a) indicate the edges of the thermal BLs for $Pr = 0.7$ and $Pr = 7$.

We find that $\langle \varepsilon _\theta \rangle _{A,t}(z)$ increases with decreasing $Pr$. This is because the volume-averaged thermal dissipation rate is related to the global heat transport (Shraiman & Siggia 1990) as

(4.7)\begin{equation} \langle \varepsilon_T \rangle_{V,t} = \frac{Nu}{\sqrt{RaPr}} . \end{equation}

As we observe $Nu \sim Pr^{0.2}$, this leads to $\langle \varepsilon _T \rangle _{V,t} \sim Pr^{-0.3}$ for a fixed $Ra$. Thus, the decrease of the thermal dissipation rate with increasing $Pr$ is consistent with the $Pr$-dependence of the Nusselt number.

We now plot in figure 6(b) the profiles of the viscous dissipation rate defined in (3.1). Similar to $\langle \varepsilon _\theta \rangle _{A,t}(z)$, the largest values of $\langle \varepsilon _u\rangle _{A,t}(z)$ are found near the horizontal plate owing to the strongly varying velocity field in the vicinity of the plates. Further, the variation of the profiles $\langle \varepsilon _u\rangle _{A,t}(z)$ in the bulk region is almost negligible compared with that in the viscous BL region near the plates. Figure 6(b) shows that $\langle \varepsilon _u\rangle _{A,t}(z)$ increases with decreasing $Pr$ for all $z$. Note that the globally averaged viscous dissipation rate is related to the Nusselt number as

(4.8)\begin{equation} \langle \varepsilon_u \rangle_{V,t} = \frac{Nu-1}{\sqrt{RaPr}} , \end{equation}

and therefore, $\langle \varepsilon _u \rangle _{V,t}$ should decrease with increasing $Pr$.

5. Characterization of the turbulence in the bulk

5.1. Isotropy in the midplane

Vorobev et al. (2005) used the ratios

(5.1)\begin{equation} G_{ij}=\frac{\langle(\partial u_i/\partial z)^2\rangle (1+\delta_{iz})}{\langle(\partial u_i/\partial x_j)^2\rangle (1+\delta_{ij})} \quad\text{with}\ i,j=x,y,z\, \end{equation}

to determine the degree of anisotropy on the level of second-order derivative moments. Flows with no variation in the vertical direction $z$ yield $G_{ij}\to 0$ (and are thus anisotropic), while $G_{ij} = 1$ for perfectly isotropic flows. The coefficient $G_{11}$, relating the in-plane derivative to a transverse derivative with respect to the vertical direction is summarized in three horizontal planes in table 3; $G_{11}$ remains nearly unity in the bulk region for $0.1 \leq z/H \leq 0.9$, but significant departures are found near the horizontal plates. Similar amplitudes follow for other combinations; see also Nath et al. (2016). We thus conclude that a plausible case exists for exploring similarities with Kolmogorov turbulence in the bulk region; see Mishra & Verma (2010) and Verma et al. (2017).

Table 3. The anisotropy coefficient $G_{11}$ in three horizontal planes; for definition see (5.1); $G_{11}$ as well as the other coefficients $G_{ij}$ remain close to unity in the bulk region between $z = 0.1H$ and $z = 0.9H$, but depart significantly as the horizontal plate is approached and the shear effects dominate.

5.2. Kinetic energy spectra

In three-dimensional turbulent flows, the kinetic energy injected at large length scales cascades towards smaller scales and eventually gets dissipated at the smallest scales by viscous action. We shall not consider the connection to the Onsager conjecture that it may be related to singularities in weak solutions of the Euler equations. In the inertial range – the range of length scales far from both the injection as well as dissipation scales – the kinetic energy spectrum $E(k)$, whose integral over all wavenumbers $k$ yields the kinetic energy, follows the standard Kolmogorov scaling

(5.2)\begin{equation} E(k) = K_{K} \varepsilon_u^{2/3} k^{{-}5/3}, \end{equation}

where $K_{K}$ is the Kolmogorov constant and $\varepsilon _u$ denotes the volume- and time-averaged kinetic energy dissipation rate.

For our purposes, it would be useful to study the behaviour of two-dimensional (2-D) energy spectrum in a horizontal plane. The 2-D Fourier transform of a field $f(x,y,z_0)$ in a horizontal plane at $z = z_0$ is defined as

(5.3)\begin{equation} f(x,y,z_0) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \hat{F}(k_x,k_y) \exp({-{\rm i} (k_x x + k_y y)})\,{\rm d}k_x \,{\rm d}k_y, \end{equation}

where $\hat {F}(k_x,k_y)$ is the Fourier mode corresponding to the wavevector ${\boldsymbol k} \equiv (k_x,k_y)$. Thus, the Fourier modes of the velocity field in midplane ${\boldsymbol u}(x,y,z=H/2) \equiv {\boldsymbol U}(x,y)$ are denoted as $\hat {\boldsymbol U}({\boldsymbol k}) \equiv [\hat {U}_x({\boldsymbol k}), \hat {U}_y({\boldsymbol k}), \hat {U}_z({\boldsymbol k})]$. The kinetic energy in a horizontal plane is equal to the sum of the energies of each Fourier mode, i.e.

(5.4)\begin{align} \frac{1}{2} \langle {\boldsymbol U}^2 \rangle_{A,t} &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{1}{2} |\hat{\boldsymbol U}(k_x,k_y)|^2 \,{\rm d}k_x \,{\rm d}k_y \nonumber\\ &= \int_{0}^{\infty} \left[ \int_{0}^{2{\rm \pi}} \frac{1}{2} |\hat{\boldsymbol U}(k,\phi_k)|^2 \,{\rm d}\phi_k \right] k \,{\rm d}k, \end{align}

with

(5.5a,b)\begin{equation} k = \sqrt{k_x^2+k_y^2} \quad\mbox{and}\quad \phi_k = \arctan(k_y,k_x). \end{equation}

Using the horizontal isotropy of the fields, the expression in the square brackets could be readily integrated to yield ${\rm \pi} |\hat {\boldsymbol U}(k)|^2$, where $|\hat {\boldsymbol U}(k)|^2/2$ is the average kinetic energy of all the Fourier modes lying in an annular region between radii $k$ and $k+{\rm d}k$. Thus, the average planar kinetic energy becomes

(5.6)\begin{equation} \frac{1}{2} \langle {\boldsymbol U}^2 \rangle_{A,t} = \int_{0}^{\infty} {\rm \pi}k |\hat{\boldsymbol U}(k)|^2 \,{\rm d}k = \int_{0}^{\infty} E(k) \,{\rm d}k, \end{equation}

where $E(k) = {\rm \pi}k |\hat {\boldsymbol U}(k)|^2$ is the 1-D kinetic energy spectrum in a horizontal plane (Peltier et al. 1996).

We compute the energy spectrum in the midplane of the low-$Pr$ flows for each instantaneous snapshot and then average the instantaneous spectra over all the available snapshots to obtain the mean kinetic energy spectrum. The energy spectra for flows with small-scale universality at different Reynolds numbers should collapse at sufficiently high Reynolds number if they are plotted against $k \eta$. Equation (5.2) in terms of the normalized wavenumber $k \eta$ reads then (Monin & Yaglom 2007) as

(5.7)\begin{equation} E(k\eta) = K_{K} (\varepsilon_u \nu^5)^{1/4} (k\eta)^{{-}5/3}. \end{equation}

We now plot the normalized energy spectra $E(k\eta ) (\varepsilon _u \nu ^5)^{-1/4}$ as a function of $k\eta$ in figure 7. The spectra for $Pr = 0.001$ are shown in figure 7(a) for all three Rayleigh numbers, whereas the spectra for a fixed $Ra = 10^6$ and $Pr = 0.001, 0.005$ and 0.021 are displayed in figure 7(b). The collapse is excellent beyond the wavenumber corresponding to the maximum of $E(k\eta )$. We show the same spectra in the normalized form $E(k) \varepsilon _u^{-2/3} k^{5/3}$ in the insets of figure 7, which confirm that they collapse onto each other for all the low-$Pr$ flows in low to moderately large wavenumber range. The plateau in the dimensionless wavenumber range $k \in [k_0, k_1]$ increases with increasing $Ra$. The inset of figure 7(a) shows that the inertial range $[k_0, k_1]$ can be estimated to be $[2, 60]$, $[4, 150]$ and $[6, 300]$ for $Ra = 10^5, 10^6, 10^7$ and $Pr = 0.001$, respectively. The inertial range increases with increasing $Ra$ and decreasing $Pr$, as expected from the increased Reynolds numbers.

Figure 7. Normalized kinetic energy spectra $E(k\eta )(\varepsilon _u \nu ^5)^{-1/4}$ in the midplane as a function of the normalized wavenumber $k\eta$ for $Pr = 0.001$ (a) and $Ra = 10^6$ (b) collapse well in the inertial as well as viscous ranges. Here, $\varepsilon _u = \langle \varepsilon _u \rangle _{V,t}$. The short dotted vertical lines indicate the normalized wavenumber $k_T \eta$, where $k_T = 2{\rm \pi} /\delta _T$ is the wavenumber corresponding to the mean thermal BL thickness $\delta _T=H/(2 Nu)$. Insets show that the corresponding spectra normalized with $\varepsilon _u^{2/3} k^{-5/3}$ exhibit plateau in the inertial range, which becomes wider with increasing $Ra$ or decreasing $Pr$. The dashed horizontal line in the insets yields the Kolmogorov constant $K_{K} \approx 1.6$ in the planar energy spectra, consistent with the value for isotropic turbulence (Sreenivasan 1995).

The energy spectra normalized in the same way for $Pr = 0.001, 0.005$ and $0.021$ for a fixed $Ra = 10^6$ are shown in the inset of figure 7(b). Again, the normalized spectra collapse quite well. A plateau can be detected for the two lower $Pr$ with the inertial range corresponding to $[4, 150]$ and $[4, 90]$ for $Pr = 0.001$ and $Pr = 0.005$, respectively. The plateau for all cases corresponds to $K_{K} \approx 1.6$, consistent with the experimental and numerical value in isotropic turbulence (Sreenivasan 1995; Yeung & Zhou 1997; Ishihara, Gotoh & Kaneda 2009).

5.3. Third-order structure function

To further explore whether the velocity fluctuations in the bulk of low-$Pr$ convection are close to Kolmogorov turbulence, we compute the third-order longitudinal structure function defined as

(5.8)\begin{equation} S_3(r) = \langle (\delta_L u(r))^3 \rangle, \end{equation}

where $\langle \,\cdot \, \rangle$ denotes an appropriate averaging, and the longitudinal velocity increment is

(5.9)\begin{equation} \delta_L u(r) = [{\boldsymbol u}({\boldsymbol x} + {\boldsymbol r}) - {\boldsymbol u}({\boldsymbol x})] \cdot \hat{\boldsymbol r}, \end{equation}

with $\hat {\boldsymbol r}={\boldsymbol r}/r$. Kolmogorov (1941a) showed that $S_3(r)$ in high-$Re$ homogeneous and isotropic turbulent flow is a universal function of the separation $r$ and varies as

(5.10)\begin{equation} S_3(r) ={-}\frac{4}{5} \varepsilon_u r \end{equation}

in the inertial range; here, and for the remainder of the work, $\varepsilon _u:=\langle \varepsilon _u \rangle _{A,t}(z=1/2)$. The longitudinal structure functions in the $x$- and $y$-directions in midplane for our highest-$Re$ flow show that they are nearly the same for small and moderate increments, so we average in the two directions. We show in figure 8 the averaged third-order structure function in the normalized form $-S_3(r)/(\varepsilon _u r)$ as a function of $r/\eta$. The figure shows that the compensated structure function tends to scale with the analytical form of $r^2$ in the beginning of the viscous range at $r/\eta \sim 1$. It also shows that the normalized structure function exhibits a plateau for an intermediate range of length scales, which implies that $S_3(r)$ indeed approximately varies as $r$ in the inertial range but the numerical value is slightly larger than $4/5$ (Kolmogorov 1941a). One possible reason for this departure is the remnant buoyancy contribution (Yakhot 1992).

Figure 8. Third-order longitudinal structure function (averaged over the $x$- and $y$-directions) in the midplane for $Pr = 0.001, Ra = 10^7$ varies as $r^2$ in the viscous range, whereas it remains nearly a constant in the inertial range. Dashed horizontal line indicates the constant 4/5 appearing in (5.10). Here, $\varepsilon _u = \langle \varepsilon _u \rangle _{A,t}(z=H/2)$.

5.4. Dissipative anomaly

The zeroth law of turbulence (or ‘dissipative anomaly’) states that the mean kinetic energy dissipation rate when scaled by a large-scale velocity such as the r.m.s. and the large length scale becomes a constant for sufficiently high Reynolds numbers (Eyink 1994; Frisch 1995; David & Galtier 2021). Figure 9 displays this rescaled energy dissipation rate

(5.11)\begin{equation} \beta=\frac{\varepsilon_u {\mathcal{L}}}{u_{rms}^3}, \end{equation}

where ${\mathcal {L}}$ is either the height $H$ of the convection layer or the integral scale $\ell$ calculated from the energy spectrum (5.6) in § 5.2 by

(5.12)\begin{equation} \ell=2{\rm \pi}\,\dfrac{\int_0^{\infty} k^{{-}1} E(k) \,{\rm d}k}{\int_0^{\infty} E(k) \,{\rm d}k}. \end{equation}

Note that $u_{rms}$ is also taken with respect to the midplane. The figure summarizes both versions of $\beta$ vs the Taylor microscale Reynolds number

(5.13)\begin{equation} R_{\lambda}= \left( \frac{25 Ra}{9\varepsilon_u^2 Pr} \right)^{1/4} \, u_{rms}^2 \end{equation}

in the bulk region of the flow. The data for different Rayleigh and Prandtl numbers collapse nicely on a curve that saturates at an approximate value of 0.2 (see the inset of the figure). While this is smaller than 0.45 found by Sreenivasan (1998) for isotropic turbulence, the result implies that the strongly disparate viscous and thermal BL widths (and thus the plume stem widths) do not matter for the driving of the turbulence cascade in the bulk of the flow. However, it also highlights the intrinsic difference in dissipation between convection and isotropic turbulence.

Figure 9. Test of the dissipative anomaly in the bulk of the convection layer using the rescaled mean kinetic energy dissipation rate $\beta$, as given by (5.11), vs the Taylor microscale Reynolds number. Asterisks use the integral scale ${\mathcal {L}}=\ell$ in (5.11), and open circles use ${\mathcal {L}}=H$. The inset expands the scale for the case ${\mathcal {L}} = H$. Symbols from left to right in each colour correspond to decreasing $Pr$. The Taylor microscale Reynolds number is given by (5.13).

6. Characteristic lengths and times of turbulent superstructures

The characteristic length scale of the dominant energy-containing structures is of the order of the size of the system when $\varGamma \approx 1$, e.g. a large-scale circulation covering the entire domain (Schumacher et al. 2016; Zürner et al. 2019; Zwirner et al. 2020). However, for $\varGamma \gg 1$, mean circulation rolls with diameters larger than $H$ are observed (Emran & Schumacher 2015); the resulting large-scale patterns of these rolls are termed turbulent superstructures of convection (Pandey et al. 2018; Stevens et al. 2018; Fonda et al. 2019; Green et al. 2020) as we stated already in § 1. Although the superstructures extend all the way from the bottom to the top plate (Pandey et al. 2018), they are conspicuous when the vertical velocity $u_z$, temperature $T$ or the vertical heat flux $u_zT$ field is visualized in the midplane (Fonda et al. 2019; Green et al. 2020). For instance, figure 1 shows turbulent superstructures for low-$Pr$ convection in the form of hot upflows and cold downflows.

The characteristic length $\lambda$ of the superstructures is the typical distance between two consecutive upwelling or downwelling regions and can be estimated using 1-D spectra of the thermal variance, kinetic energy and convective heat flux in a horizontal plane (Hartlep et al. 2003; Pandey et al. 2018; Stevens et al. 2018; Green et al. 2020; Krug et al. 2020). It can also be estimated by computing the two-point auto-correlation function of $u_z$ or $T$ in a horizontal plane and identifying the location of the first minimum, corresponding to $\lambda /2$ (Pandey et al. 2018).

We compute the 1-D power spectra of the vertical velocity, temperature and convective heat flux, all averaged with respect to the azimuthal angle. They are given by

(6.1)$$\begin{gather} S_U(k) = \frac{1}{2{\rm \pi}} \int_0^{2{\rm \pi}} |\hat{U}_z(k,\phi_k)|^2 \,{\rm d} \phi_k, \end{gather}$$

(6.2)$$\begin{gather}S_\varTheta(k) = \frac{1}{2{\rm \pi}} \int_0^{2{\rm \pi}} |\hat{\varTheta}(k,\phi_k)|^2 \,{\rm d} \phi_k, \end{gather}$$

(6.3)$$\begin{gather}S_{U\varTheta}(k) = \frac{1}{2{\rm \pi}} \int_0^{2{\rm \pi}} \mathrm{Re}[\hat{U}_z^*(k,\phi_k) \hat{\varTheta}(k,\phi_k)] \,{\rm d} \phi_k, \end{gather}$$

where $\hat {U}_z(k,\phi _k)$ and $\hat {\varTheta }(k,\phi _k)$ are, respectively, the 2-D Fourier transforms of $u_z(x,y,H/2)$ and $\theta (x,y,H/2)$. Note that the spectrum $S_U(k)$ differs from the energy spectrum $E(k)$ defined in § 5.2. We find that the three spectra exhibit a peak at nearly the same wavenumber $k^\omega _{max}$ corresponding to the maximum of $S_\omega (k)$ (see figure 13), yielding the characteristic spatial scale $\lambda _\omega = 2{\rm \pi} /k^\omega _{max}$ of the superstructures, with $\omega = \{ U, \varTheta, U\varTheta \}$.

Pandey et al. (2018, 2021) found that the characteristic scales $\lambda _U$ and $\lambda _\varTheta$ do not always agree with each other, the contrast being usually larger at moderate and high Prandtl numbers. Here, we extract the spatial scale of superstructures from the power spectrum of the convective heat flux (see also Hartlep et al. 2003; Krug et al. 2020). It was observed in Pandey et al. (2018) that $\lambda$ is a function of $Pr$ and the maximum of $\lambda (Pr)$ was found at $Pr \approx 7$ for $Ra = 10^5$. In the present simulations (in which the $Pr$ range has been extended down to $0.001$ and $Ra$ by two orders of magnitude), the corresponding values of $\lambda$ as a function of $Pr$ are plotted in figure 10(a) for all simulations of table 1. We find that $\lambda$ decreases with decreasing $Pr$ but they seem to level off at $\lambda \simeq 3H$ for the lowest Prandtl numbers of 0.005 and 0.001; see also figure 11 where the magnitude of $\lambda$ is displayed in the streamline plot. It remains to be studied as to how this scale arises starting from the onset of RBC where $\lambda = 2.02 H$ independently of $Pr$ (Chandrasekhar 1981), as indicated in figure 10(a) as a dashed horizontal line. We also refer to Appendix C for further details.

Figure 10. Characteristic spatial (a) and temporal (b) scales of the turbulent superstructures as a function of $Pr$. Both temporal and spatial scales decrease with decreasing Prandtl number, but seem to level off for $Pr \leq 0.005$.

Figure 11. Plots of streamlines for snapshots at (a) $Ra=10^5$, (b) $10^6$ and (c) $10^7$, all at $Pr=0.001$. The view is from the top. Only one half of the cross-section is shown. The corresponding values of $\lambda$ are indicated by a horizontal white bar in each panel.

The characteristic temporal scale ($\tau$) of the superstructures is long compared with the time scale of the turbulent fluctuations. Thus, time averaging the velocity or temperature fields over a time interval $\tau$ yields coarse-grained fields, which are nearly devoid of the small-scale fluctuations. The coarse-grained large-scale structures evolve on the time scales which are of the order of $\tau$ (Pandey et al. 2018; Fonda et al. 2019), related to the time for a fluid parcel to complete a circulation. It is computed as

(6.4)\begin{equation} \tau = 3 \frac{{\rm \pi}(\lambda/4 + H/2)}{u_{rms}}, \end{equation}

where the quantity in the numerator is the circumference of superstructure rolls with elliptical cross-section (Pandey et al. 2018). The factor of three in the above expression arises from the fact that the circulation time in an extended convection flow is not fixed but exhibits a broad distribution with stretched-exponential tails in the Lagrangian frame of reference along massless tracer trajectories (Schneide et al. 2018; Vieweg et al. 2021b). The scale computed using expression (6.4) is plotted in figure 10(b) as a function of $Pr$. We find that $\tau$ increases with $Pr$, consistent with the fact that the Reynolds number, and thus the characteristic velocity $u_{rms}$, decreases with increasing $Pr$, requiring a longer time for fluid parcels to complete a circulation. Note that $\lambda$ in (6.4) also increases with increasing $Pr$, but the increase is not as significant as the decrease in $u_{rms}$.

7. Conclusions and outlook

Our focus here has been the Rayleigh–Bénard convection in a horizontally extended layer for molecular Prandtl numbers as small as $Pr=10^{-3}$, which go beyond those accessible in controlled laboratory experiments and approach astrophysical conditions. We extended the parameter space of previous works by Pandey et al. (2018) and Fonda et al. (2019) by DNS, both towards lower $Pr$ and higher $Ra$, and thus determined more conclusively various parameter dependencies such as global heat and momentum transports, temperature fluctuations, as well as the kinetic and thermal dissipation rates. Among others, these results provide a test for existing predictions by the theory of Grossmann & Lohse (2001) and Stevens et al. (2013). Comparisons show that the predictions for the global heat and momentum transports as a function of the Prandtl number at fixed Rayleigh number are in fair agreement.

We also found that the Nusselt number decreases as $Nu \sim Pr^{0.19}$, whereas the Reynolds number increases approximately as $Re \sim Pr^{-0.6}$ when the Prandtl number decreases from $Pr = 0.7$ to 0.001, as detailed in table 2. The dimensionless mean thermal and kinetic energy dissipation rates also decrease with increasing $Pr$ and their scaling behaviours are consistent with that of the global heat transport. We studied the depth dependence of these quantities, and found that, due to a high diffusive temperature field in convection at very low $Pr$, the bulk fluid is not mixed well and a significant vertical temperature gradient occurs in the bulk region, even for the highest accessible Rayleigh number.

The highly inertial fluid turbulence in the bulk of low-$Pr$ convection layer was studied by examining the kinetic energy spectra, the 4/5th law, local isotropy and a test of the dissipative anomaly. The results suggest that the fluid turbulence in the bulk for the lowest Prandtl number is close to the classical Kolmogorov turbulence. This implies that the temperature field behaves as a passive and highly diffusive scalar stirred by a highly turbulent flow; the impact of thermal plumes, which are the unstable fragments of the thick thermal BLs, can be considered an efficient large-scale forcing for turbulence, dominantly in the lower wavenumber part of the inertial range ($k_T$ in figure 7). Some differences do remain. In particular, the asymptotic value of the rescaled mean kinetic energy dissipation rate falls below that in isotropic turbulence (Sreenivasan 1998), which suggests that BLs do matter, despite the nearness to isotropy in the central region.

Our DNS results are fully resolved and can thus have implications for the modelling of small-scale turbulence in coarser-grid simulation studies of mesoscale convection, particularly for the development of subgrid-scale models that go beyond the mixing length theory of Prandtl (1925). This class of algebraic turbulence models is still a workhorse in astrophysical simulations; see for example discussions of their limitations and extensions in Miesch (2005) and Kupka & Muthsam (2017), or recent extensions by Brandenburg (2016). Finally, we stress that the convection considered here does not incorporate complexities such as rotation, magnetic field, varying molecular transport coefficients or curvature, which are present in geophysical and astrophysical settings. They would have to be included to yield realistic models in these specific instances. The present focus has been the exploration of the effects of low Prandtl numbers, which is an important facet of these flows. A more detailed study of these points in connection with small-scale intermittency in low Prandtl number RBC flows is underway, and will be reported elsewhere.