2.1. Model and parameters

We consider a horizontal fluid layer of depth $H$, kinematic viscosity $\nu$, thermal diffusivity $\kappa$, specific heat capacity ${C_p}$ and thermal expansion coefficient $\alpha$, all of which are considered constant. The layer is subjected to the vertical gravitational acceleration of (uniform) magnitude $g$. We consider the Boussinesq approximation, where the density $\rho$ is constant and uniform, except in the expression of the buoyancy force. Volumic heat sources and sinks are present in the bulk of the fluid. In the experiment by Lepot etal. (Reference Lepot, Aumaître and Gallet2018) and Bouillaut etal. (Reference Bouillaut, Lepot, Aumaître and Gallet2019), a powerful spotlight shines visible light through the transparent bottom of the tank, into dyed water where absorption occurs over a typical length $\ell$. Following Beer–Lambert's law, the intensity of light decreases exponentially with the height $Z$, measured from the bottom surface. Thus, we model the experimental radiation-induced heating with an exponentially decreasing volumic heat source in the present numerical work

(2.1)\begin{equation} Q_L(Z) = \frac{P}{\ell}\frac{\exp\left(-Z/\ell\right)}{1-\exp\left(-H/\ell\right)}, \end{equation}

where $P$ is the incident light flux at the bottom of the tank (power per unit area), which gets deposited into the fluid through absorption. In the laboratory experiment, this internal heat source yields ‘secular heating’ of the fluid, which amounts to a quasi-stationary state with homogeneous internal cooling $Q_c(Z) = P/H$ balancing the radiative heat input, so that the net power deposited in the fluid layer vanishes: $\int _0^{H} ( Q_L(Z) - Q_c(Z) ) \,\mathrm {d}Z = 0$. We stress the fact that, motivated by the experimental realisations, we are concerned with a distribution of heat sources that takes on finite values at the boundaries. In that respect, the problem of interest here crucially differs from previous studies by Barker, Dempsey & Lithwick (Reference Barker, Dempsey and Lithwick2014), who consider a distribution of heat sources and sinks that vanishes at the boundaries.

The characteristics of the set-up can be cast into three dimensionless control parameters, namely the Prandtl number ${Pr}$, the flux-based Rayleigh number ${Ra_Q}$ and the dimensionless light absorption length $\tilde{\ell}$

(2.2a–c)\begin{equation} {Pr} = \frac{\nu}{\kappa},\quad {Ra_Q}=\frac{\alpha g PH^{4}}{\rho {C_p}\kappa^{2} \nu}, \quad \tilde{\ell} = \frac{\ell}{H}. \end{equation}

The response of the fluid to the imposed flux may be characterised by the magnitude of the temperature difference ${\rm \Delta} T$ that appears between the bottom plate $(Z=0)$ and mid-depth $(Z=H/2)$, and the characteristic flow velocity $U$. To make contact with the numerous RBC literature, one can then define the Reynolds number ${Re}$, the temperature-based Rayleigh number ${Ra}$ and the Nusselt number ${Nu}$

(2.3a–c)\begin{equation} {Re} = \frac{UH}{\nu}, \quad {Ra} = \frac{\alpha g {\rm \Delta} T H^{3}}{\nu\kappa}, \quad {Nu} = \frac{PH}{\rho {C_p} \kappa {\rm \Delta} T}, \quad \mathrm{with}\ {Ra}_Q = {Ra}{Nu}. \end{equation}

2.2. Dimensionless equations and parameters

The Boussinesq equations are written in a dimensionless form by counting lengths in units of the depth $H$, time in units of diffusive time $H^{2}/\kappa$ and temperature in units of $\nu \kappa / \alpha g H^{3}$. Thus, the governing equations for the solenoidal velocity field ($\boldsymbol {\nabla \cdot u}=0$) and temperature field $\theta$ are

(2.4a)\begin{gather} {Pr}^{-1} \left(\partial_t \boldsymbol{u} + \boldsymbol{u\cdot \nabla u}\right) = - \nabla p + \theta\,\hat{\boldsymbol{e}}_z + \nabla^{2} \boldsymbol{u}, \end{gather}

(2.4b)\begin{gather} \partial_t \theta + \boldsymbol{u\cdot \nabla}\theta = \nabla^{2} \theta + Ra_Q \left(\frac{\exp\left(- z/\tilde{\ell}\right)}{\tilde{\ell}\left[1-\exp\left(-1/\tilde{\ell}\right)\right]} - 1\right), \end{gather}

where $z=Z/H$. These governing equations are supplemented with adiabatic boundary conditions for temperature: $\partial _z \theta =0$ at $z=0,\, 1$. Kinematic boundary conditions are always impenetrable ($w=0$) and either no slip ($u=v=0$) or stress free ($\partial _z u = \partial _z v = 0$) at the top and bottom. Finally, we consider periodic boundary conditions in the horizontal.

2.3. Numerical procedure

Direct numerical simulations of equations (2.4) are performed with the High-Performance-Computing code Coral, which employs a Chebyshev–Fourier–Fourier decomposition in space and a family of implicit–explicit time-stepping schemes. The boundary conditions in the vertical direction $z$ are enforced using Galerkin basis recombination. Computations are initialised with either small-amplitude noise or a previous solution and are carried out until a statistically stationary regime is clearly identified. Coral was already used for studying CISS in Miquel etal. (Reference Miquel, Lepot, Bouillaut and Gallet2019), where it is benchmarked against analytical solutions. The non-periodic version of Coral has also been benchmarked against the numerical solutions obtained by Lepot etal. (Reference Lepot, Aumaître and Gallet2018) using the solver Dedalus (Burns etal. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), for ${Pr}=1$ and $7$.

3.1. Bulk transport

As a fluid element circles around a turbulent convection cell, it rapidly heats up in the region $z \lesssim \tilde{\ell}$, while it slowly cools down in the region $z \gtrsim \tilde{\ell}$. In other words, the heat input inside the region $z \lesssim \tilde{\ell}$ is carried outside the heating region by the flow, which amounts to a balance between the source term and the horizontal advection term in the heat equation

(3.1)\begin{equation} \rho {C_p}U \partial_x T \sim \frac{P}{\ell}. \end{equation}

Upon inserting the free-fall velocity estimate, the balance above yields an estimate ${\rm \Delta} T_{bulk}$ for the temperature difference between the heating region $z \lesssim \tilde{\ell}$ and the bulk region $z \gtrsim \tilde{\ell}$. In dimensional and dimensionless form:

(3.2a,b)\begin{equation} \left( {\rm \Delta} T_{bulk} \right) ^{3} \sim \frac {P^{2}H}{\alpha g \rho^{2} C_{p}^{2} \ell^{2}},\quad {Nu} \sim \tilde{\ell} \sqrt{{Ra}{Pr}}. \end{equation}

The reasoning above and the associated ${\rm \Delta} T_{bulk}$ do not take into account the velocity boundary conditions at the bottom plate. In the next section, we show that the latter can induce a possibly stronger temperature difference ${\rm \Delta} T_{BL}$ between the immediate vicinity of the bottom boundary and the heating region $z \lesssim \tilde{\ell}$. The total temperature difference between the bottom plate and the bulk region $z \gtrsim \tilde{\ell}$ is then estimated as ${\rm \Delta} T={\rm \Delta} T_{BL}+{\rm \Delta} T_{bulk}$.

3.2. Stagnant layers near no-slip boundaries

The scaling law (3.2) may be modified in the case of a no-slip bottom boundary condition. Indeed, the stagnant layer in the immediate vicinity of the bottom plate accumulates heat, and without efficient advection outside this region, the heat can only be diffused away. This leads to a thermal boundary layer of thickness $\delta _{\Theta }$ . Balancing the heat source term with the diffusive one estimated inside the boundary layer leads to

(3.3a)\begin{equation} \rho {C_p} \kappa \frac{{\rm \Delta} T_{BL}}{\delta^{2}_\Theta} \sim \frac{P}{\ell}, \end{equation}

yielding

(3.3b)\begin{equation} {\rm \Delta} T_{BL} = \frac{PH^{2}}{\rho {C_p} \kappa \ell }\frac{\delta_\Theta^{2}}{H^{2}}. \end{equation}

A final relation between the Nusselt number and the Rayleigh and Prandtl numbers is obtained by supplementing (3.3) with a scaling estimate of the boundary-layer thickness $\delta _\Theta$. The latter depends crucially on the Prandtl number: in the case of low Prandtl fluids, a naive laminar estimate would lead to a thermal boundary layer thicker than the kinematic one. However, one expects vigorous turbulent mixing outside the kinematic boundary layer, so that both boundary layers end up having the same thickness. We obtain $\delta _\Theta \sim H / \sqrt {{Re}}$. Upon estimating the Reynolds number using the free-fall velocity scaling (based on ${\rm \Delta} T_{BL}$ in a worst-case scenario), equation (3.3) becomes

(3.4)\begin{equation} \left( {\rm \Delta} T_{BL} \right) ^{3} \sim \frac {P^{2}H}{\alpha g \rho^{2} C_{p}^{2} \ell^{2}} \frac{\nu^{2}}{\kappa^{2}}, \end{equation}

which remains negligible compared with the bulk temperature jump from equation (3.2) when ${Pr} \ll 1$. As a consequence, the scaling law ${Nu} \sim \tilde{\ell} \sqrt {{Ra} {Pr}}$ remains valid and is not modified by the stagnant bottom layer for ${Pr} \ll 1$.

The situation is crucially different for high Prandtl fluids. The thermal boundary layer is much thinner than, and nested into, the kinetic one. The flow in the vicinity of the no-slip bottom surface is perceived by the thermal boundary layer as a parallel shear flow, the bottom shear $S(0)$ being estimated as $S(0) \sim U / \delta _\nu \sim U \sqrt {{Re}}/ H \sim \sqrt {{U^{3}}/{H\nu }}$. The boundary-layer thickness is set by a balance between slow advection by the shear, at a speed $S(0)\delta _\Theta$ , and vertical diffusion

(3.5a)\begin{equation} \rho {C_p}S(0) \delta_\Theta \frac{{\rm \Delta} T_{BL}}{H} \sim \frac{\rho {C_p} \kappa {\rm \Delta} T_{BL}}{\delta_\Theta^{2}}, \end{equation}

so that

(3.5b)\begin{equation} \frac{\delta_\Theta}{H} \sim\left(\frac{\kappa}{S(0)H^{2}}\right)^{1/3} \sim \frac{1}{{Re}^{1/2}{Pr}^{1/3}}. \end{equation}

Inserting this expression for the thermal boundary-layer thickness into equation (3.3), one gets the temperature jump associated with the sheared thermal layer and the associated dimensionless heat flux

(3.6 a,b)\begin{equation} \left( {\rm \Delta} T_{BL} \right) ^{3} \sim \frac {P^{2}H}{\alpha g \rho^{2} C_{p}^{2} \ell^{2}}\left(\frac{\nu}{\kappa}\right)^{2/3},\quad {Nu} \sim \tilde{\ell} {Ra}^{1/2}{Pr}^{1/6}. \end{equation}

Noticeably, a comparison with equation (3.2) informs us that the temperature jump in the stagnant boundary layer dominates the temperature jump in the bulk: ${\rm \Delta} T \simeq {\rm \Delta} T_{BL}$ when $Pr \gg 1$. This observation results in a modification of the Nusselt number scaling law due to boundary-layer corrections, see equation (3.6).

The last case that requires our attention is the case of high Prandtl number flows with a stress-free bottom surface. Within the thermal boundary layer, the flow is well approximated by a uniform parallel flow with velocity $U\sim \sqrt {\alpha g {\rm \Delta} T_{BL} H}$. Here, in view of finding hypothetical boundary-layer corrections, we assume again that the overall temperature difference is dominated by the contribution from the boundary layer ${\rm \Delta} T_{BL}$, before testing this assumption a posteriori. A triple balance between horizontal advection, diffusion through the layer and heating yields

(3.7)\begin{equation} \rho {C_p} \sqrt{\alpha g {\rm \Delta} T_{BL} H} \frac{{\rm \Delta} T_{BL}}{H} \sim \rho {C_p} \kappa \frac{{\rm \Delta} T_{BL}}{\delta_\Theta^{2}} \sim \frac{P}{\ell}. \end{equation}

This balance is similar to the bulk balance (3.1). Some straightforward algebra confirms that a transport law identical to (3.2) is found. We conclude that no boundary-layer-induced corrections to the heat-transport scaling law are necessary in this case.

To summarise, the scaling theory results in the following predictions for the dimensionless heat flux:

(3.8)\begin{equation} {Nu} = \left\{\begin{array}{ccc} {c_1} \, \tilde{\ell} \sqrt{{Ra} {Pr}} & \mbox{with stress-free b. c.,} & \forall {Pr},\\ {c_2} \,\tilde{\ell} \sqrt{{Ra} {Pr}} & \mbox{with no-slip b. c.,} & {Pr} \ll {Pr}_*, \\ {c_3} \,\tilde{\ell} {Ra}^{1/2} {Pr}^{1/6} & \mbox{with no-slip b. c.,} & {Pr}\gg {Pr}_*, \end{array}\right. \end{equation}

where the dimensionless prefactors $c_i$ remain to be determined. The transition Prandtl number ${Pr}_*$ is estimated by equating the two no-slip predictions at $Pr=Pr_*$, which yields $Pr_*=(c_3/c_2)^{3}$.

4.1. Heat flux scalings

Guided by the prediction (3.8) and the observations by Lepot etal. (Reference Lepot, Aumaître and Gallet2018) for ${Pr}=1$ and $7$, we represent in figure 1(a) the Nusselt number ${Nu}$ compensated by $\sqrt {Ra}$ for data obtained in a suite of DNS runs with almost logarithmically equispaced Prandtl number values: ${Pr} \in \{0.003,\, 0.01,\, 0.03,\, 0.1,\, 0.3,\, 1,\,3,\,10\}$. We observe that, for a fixed Prandtl number (a given colour on the figure), the data form a plateau as the thermal forcing ${Ra_Q}$ varies. This confirms that the scaling exponent $\gamma =0.5$ is universal and does not depend on the Prandtl number, in line with prediction (3.8).

Figure 1. (a) Compensated heat flux ${Nu} / \sqrt {{Ra}}$ as a function of the Rayleigh number ${Ra}$ for no-slip boundary conditions. The colour code is for the Prandtl number ${Pr}$, while the heating intensity as measured by ${Ra_Q}/{Pr}^{2}$ is encoded in the shape of the symbols. For convenience, we provide a correspondence between the parameter ${Ra_Q}/{Pr}^{2}$ and the Reynolds number ${Re}$. This colour/shape code is employed consistently thereafter throughout the paper. (b) Compensated heat flux ${Nu}/ \sqrt {{Ra} {Pr}}$ as a function of ${Pr}$ for no-slip (filled symbols) and stress-free (open symbols) bottom boundaries.

We now turn to the central question of this paper – the Prandtl number dependence – which amounts to analysing the height of the plateaus in figure 1(a). In this logarithmic representation, the plateaus are separated by a distance that seems to vary with the Prandtl number.

To better characterise this effect, we turn to figure 1(b) where we plot the heat flux compensated by $\sqrt {{Ra}{Pr}}$ – the scaling prediction (3.8) in the absence of boundary-layer correction – as a function of the Prandtl number. Upon this rescaling, the data points for no-slip boundary conditions (filled symbols) fall on a master curve that exhibits: (i)a plateau for ${Pr}\lesssim 0.04$, providing evidence for the ${Nu}\sim \sqrt {{Ra}{Pr}}$ scaling at low Prandtl numbers; (ii) a scaling law with slope $-1/3$ for ${Pr} \gtrsim 0.04$ which, when combined with the rescaling of the figure, suggests the ${Nu} \sim {Ra}^{1/2}{Pr}^{1/6}$ scaling predicted at high Prandtl numbers.

By contrast, the data points obtained for stress-free boundary conditions (open symbols) reasonably form a plateau, as they exhibit variations of less than $30\,\%$ when ${Pr}$ spans three and a half decades. This observation validates the mixing-length scaling ${Nu} \sim \sqrt {{Ra} {Pr}}$ at all Prandtl numbers for stress-free boundary conditions.

To summarise, the numerical data are fully compatible with the prediction (3.8), with the following values for the dimensionless prefactors and transition Prandtl number: $c_1=c_2 \simeq 0.38$, $c_3\simeq 0.13$, which leads to $Pr_*=(c_3/c_2)^{3} \simeq 0.04$.

4.2. Boundary-layer structure: no-slip case

We have established that, in the high Prandtl number regime ${Pr}\gtrsim 0.04$, thermal transfer laws differ for a no-slip and a stress-free bottom boundary and are accurately captured by the scaling prediction (3.8). To validate a posteriori the underlying assumptions of the model, we now examine the structure of the boundary layers in the DNS flows.

The structure of the thermal boundary layer is illustrated in figure 2 for a value of the Prandtl number ${Pr}=3$ that lies comfortably within the high Prandtl number regime, where a clear departure is observed between no-slip and stress-free bottom boundaries. In this figure, we display the vertical temperature profiles $\bar{\theta}$, where the overline denotes an average along the horizontal directions and over time. Upon examination, the temperature profiles are strikingly similar for no-slip and stress-free cases in the bulk of the flow, while the no-slip case develops strong boundary layers associated with a large temperature jump at the bottom of the layer.

Figure 2. (a) Temperature profiles averaged in time and along horizontal directions $x$ and $y$, for ${Pr}=3$, ${Ra_Q}=3\times 10^{10}$, with either a no-slip (blue) or a stress-free (orange) bottom boundary. (b) Time and horizontally averaged contributions to the total heat flux. Solid lines represent the diffusive flux $-\overline {\partial _z \theta }$, whereas dashed lines represent the convective flux $\overline {w\theta }$. Coloured circles materialise $z_{\textit{max}}$, the altitude where the diffusive flux is maximum. For a no-slip bottom boundary (top panel), the thin horizontal black line represents the altitude $z_*$ where the diffusive and the convective fluxes are equal (see text).

Beyond the mere temperature profiles, we analyse the thermal energy fluxes in thevicinity of the bottom bounding surface in figure 2(b). For no-slip boundary conditions,the diffusive flux $\overline {\partial _z \theta }(z)$ grows linearly with the distance from the wall $z$. As a consequence, the diffusive flux strongly dominates the purely convective flux $\overline {w\theta }(z)$, which exhibits a slower quadratic growth with $z$. As $z$ increases, the diffusive flux reaches a maximum at distance $z_{max}$ and decreases again. Advection and diffusion contribute in equal proportions to the heat transfer at a comparable height ${z_*}$, where $\overline {\partial _z \theta } ({z_*}) = \overline {w\theta }({z_*})$. Above ${z_*}$ is the bulk, where pure convection proves the most efficient mechanism.

Being able to extract a thermal boundary-layer thickness from the numerical solutions, we now turn our attention to the validity of the sheared boundary-layer assumption (3.5) near a no-slip boundary. The dimensionless bottom shear $S(0)H^{2}/\kappa$ is estimated as

(4.1)\begin{equation} {S(0)H^{2}}/{\kappa} = \overline{\partial_z\sqrt{u^{2}(x,y,z,t) + v^{2}(x,y,z,t)}}|_{z=0}. \end{equation}

The boundary-layer thickness measured by $z_{max}$ is displayed in figure 3 as a function of this dimensionless bottom shear, for all the data gathered with no-slip boundary conditions. Agood collapse is observed, and the data for large Prandtl number validate the shear layer structure assumed in (3.5). The approach to this asymptote proves rather slow, because the kinetic boundary layer is only moderately thicker than the thermal one for the highest $Pr$ achieved in this study (approximately $3$ times thicker for $Pr=10$).

Figure 3. Structure of the boundary layers as ${Pr}$ and ${Ra_Q}$ vary. We represent the dimensionless thickness of the thermal boundary layer $\delta _\Theta /H$ defined as the depth where the diffusive flux is maximum (see figure 2b), represented as a function of the dimensionless shear at the bottom plate $S(0)H^{2}/\kappa =\overline {\partial _z\sqrt {u^{2}+v^{2}}}$. Dashed line: prediction from (3.5).

4.3. Boundary-layer structure: stress-free case

Returning to figure 2(b), the examination of the energy flux profiles in the stress-free case yields a much different scenario. In contrast with the no-slip case, both the convective flux $\overline {w\theta }$ and the diffusive flux $\overline {\partial _z \theta }$ grow linearly with the distance $z$ to the bottom boundary. As a first consequence, the ratio between these two fluxes remains asymptotically of $O(1)$ as the distance to the wall decreases. Moreover, energy transfers are dominated by pure convection as $\overline {w\theta }>\overline {\partial _z\theta }$ for any altitude $z$. All these observations are in strong support of the working hypotheses used for deriving prediction (3.8), and thus shed light on the (non-purely diffusive) nature of the boundary layers that result in the ${Nu} \sim \sqrt {{Ra} {Pr}}$ scaling observed for all ${Pr}$ in figure 1.