3.1. Mean flow statistics

The heat transfer characteristics of natural convection flows are commonly quantified using the Nusselt number, which, based on the momentum integral thickness, is defined as

(3.1) \begin{equation} {{Nu}}_\delta = \frac {q_{w}\delta }{\rho C_p \kappa \theta _w}, \end{equation}

where $C_p$ represents the isobaric specific heat and $q_w=-\rho C_p\kappa (\partial \overline {\theta }/\partial y )_w$ is the wall heat flux. Figure 2 compares the Nusselt number evolution obtained from our temporally developing DNS with that of steady, spatially developing flows (Vliet & Liu Reference Vliet and Liu1969; Tsuji & Nagano Reference Tsuji and Nagano1988; Tsuji & Kajitani Reference Tsuji and Kajitani2006) and other temporally evolving cases (Abedin et al. Reference Abedin, Tsuji and Hattori2009; Ke et al. Reference Ke, Williamson, Armfield and Komiya2023) for $ \textit{Pr} = 0.71, 4.16$ and $6$ . In the laminar regime, since the mean temperature and velocity profiles of the NCBL follow the analytical solution (2.3), the momentum integral thickness $\delta$ and the corresponding Nusselt number can be obtained as

(3.2a,b) \begin{equation} \delta = 2\mathcal{C}\sqrt {\kappa t} , \quad {{Nu}}_\delta =\frac {2\mathcal{C}}{\theta _w\sqrt {\pi }}, \end{equation}

where $\mathcal{C}$ is an integral constant resulting from integrating the streamwise velocity (2.3b ),

(3.3) \begin{equation} \mathcal{C} = \frac {2\left (1-\sqrt { \textit{Pr} }\right )}{3\sqrt {\pi } \left [{\text{erfc}}(\eta _m)-{\text{erfc}}(\eta _m/\sqrt { \textit{Pr} })\right ]}, \end{equation}

and $\eta _m$ is the maximum velocity location in the similarity coordinate $\eta$ that satisfies

(3.4) \begin{equation} \sqrt { \textit{Pr} }\,{\text{ierfc}}({\eta _m}) = {\text{ierfc}}\left (\frac {\eta _m}{\sqrt { \textit{Pr} }}\right ). \end{equation}

Both $\mathcal{C}$ and $\eta _m$ are implicit functions of $ \textit{Pr}$ , for which no closed-form analytical expressions are available. As a result, ${{Nu}}_\delta$ also inherits an implicit dependence on $ \textit{Pr}$ in the laminar regime. Despite the absence of analytical solutions for $\eta _m$ and $\mathcal{C}$ , we can derive the asymptotic limits of these similarity constants, along with the corresponding behaviour of ${{Nu}}_\delta$ in low- and high-Prantdl number regimes. These results are summarised in equation (3.5) later, with detailed derivations provided in Appendix A:

(3.5a,b,c) \begin{equation} \eta _m=\sqrt { \textit{Pr} \mathcal{W}_0\left (\frac {1}{2\sqrt { \textit{Pr} }}\right )},\quad \mathcal{C}\rightarrow \frac {2}{3\sqrt {\pi }}, \quad {{Nu}}_\delta \rightarrow \frac {4}{3\pi }, \quad \text{as } \textit{Pr} \rightarrow 0 \end{equation}

(3.5d,e,f) \begin{align} \eta _m=\sqrt {\mathcal{W}_0\left (\frac {\sqrt { \textit{Pr} }}{2}\right )}, \quad & \mathcal{C}\rightarrow \frac {2}{3\sqrt {\pi }}\left (\sqrt { \textit{Pr} }-1\right )\propto \sqrt { \textit{Pr} },\nonumber\\ & {{Nu}}_\delta \rightarrow \frac {4}{3\pi }\left (\sqrt { \textit{Pr} }-1\right )\propto \sqrt { \textit{Pr} }, \text{as } \textit{Pr} \rightarrow \infty . \end{align}

Here, $\mathcal{W}_0(x)$ is the principle branch of the Lambert W function. In the intermediate $ \textit{Pr}$ range, the similarity constant $\eta _m$ can be accurately approximated by

(3.6) \begin{equation} \eta _m \approx \sqrt {A\mathcal{W}_0\left (B \textit{Pr} ^{M}\right )}, \end{equation}

with which one can obtain $\mathcal{C}$ and ${{Nu}}_\delta$ using (3.3) and (3.2). Here, $A=0.525$ , $B=0.512$ and $M=0.688$ are all empirical constants obtained by calibrating the approximation in (3.6) against the numerical solution of (3.4) in the range $10^{-3}\leqslant \textit{Pr} \leqslant 10$ . The detailed derivation and a full comparison between the approximation and numerical results are provided in Appendix A. Table 3 lists the numerical solution to (3.4) for selected values of $ \textit{Pr}$ along with the corresponding $\mathcal{C}$ and ${{Nu}}_\delta$ predictions for the laminar NCBL flow. From figure 2, our DNS results are seen to follow the laminar values of ${{Nu}}_\delta$ listed in table 3, as indicated by the coloured dotted lines. This agreement confirms that the early-time flow dynamics are well described by the laminar similarity solution, which captures the $ \textit{Pr}$ dependence of heat transfer prior to transition. The onset of laminar–turbulent transition, however, is marked by a clear deviation from these laminar predictions, occurring at progressively higher ${{Ra}}_\delta$ for increasing $ \textit{Pr}$ in figure 2. For the range of Prandtl numbers considered ( $0.71\leqslant \textit{Pr} \leqslant 6$ ), the transition begins between ${{Ra}}_\delta \approx 10^4$ and ${{Ra}}_\delta \approx 10^5$ , which is broadly consistent with the data reported by Abedin et al. (Reference Abedin, Tsuji and Hattori2009). As the flow undergoes transition, all $ \textit{Pr}$ cases in figure 2 exhibit qualitatively similar trends, with a $ \textit{Pr}$ -dependent offset in the magnitude of ${{Nu}}_\delta$ . Beyond the transitional regime, the flow enters a classical turbulent regime (Ke et al. Reference Ke, Williamson, Armfield and Komiya2023; Wells Reference Wells2023), in which all $ \textit{Pr}$ cases converge towards an empirical 1/3-power-law scaling, indicated by the grey dotted line. The precise ${{Ra}}_\delta$ marking the end of transition (and therefore the beginning of the classical turbulent regime) is difficult to determine accurately, owing to the gradual changes in the ${{Nu}}_\delta$ scaling behaviour and statistical fluctuations inherent to turbulence. Nevertheless, it is clear from figure 2 that the completion of transition exhibits a strong $ \textit{Pr}$ dependence, occurring between ${{Ra}}_\delta \approx 10^6$ and ${{Ra}}_\delta \approx 10^7$ , with the transition ending at progressively higher ${{Ra}}_\delta$ as $ \textit{Pr}$ increases.

Table 3. Selected values of the Prandtl number $ \textit{Pr}$ , along with the corresponding $\eta _m$ and $\mathcal{C}$ and laminar Nusselt number ${{Nu}}_\delta$ computed from numerical solutions of (3.4) and (3.3).

To gain further insight into the underlying flow behaviour across Prandtl numbers, we now examine the mean velocity and temperature profiles in the wall-normal direction. While the global heat transfer characteristics were discussed in terms of the Rayleigh number ${{Ra}}_\delta$ , it is more appropriate to compare local flow structures at matched ${{Gr}}_\delta$ . This choice ensures that the boundary layers have comparable instantaneous (local) momentum boundary layer thicknesses $\delta$ , hereby enabling a more meaningful comparison of mean profile structures and isolating the effect of $ \textit{Pr}$ on the shape of the profiles.

Figure 3. Comparison of mean profiles at matched Grashof number ${{Gr}}_\delta \approx 1.58\times 10^6$ for different Prandtl numbers plotted against wall-normal coordinate scaled by the momentum boundary layer thickness $y/\delta$ . (a) Normalised streamwise velocity profiles $\overline {u}/\overline {u}_m$ ; (b) mean temperature profiles $\theta$ . Results from the present DNS for $ \textit{Pr} =4.16$ and $ \textit{Pr} =6$ are shown alongside the data from Tsuji & Kajitani (Reference Tsuji and Kajitani2006), Abedin et al. (Reference Abedin, Tsuji and Hattori2009), Ke et al. (Reference Ke, Williamson, Armfield and Komiya2023) for $ \textit{Pr} =0.71$ and $ \textit{Pr} =6$ .

Figure 3 compares the mean flow statistics at ${{Gr}}_\delta \approx 1.58\times 10^6$ , with all profiles scaled by the momentum integral thickness $\delta$ . We note that while figure 3 shows normalised profiles, the streamwise velocity magnitudes are seen to increase with increasing $ \textit{Pr}$ at the same ${{Gr}}_\delta$ . Once this $ \textit{Pr}$ dependence is accounted for by normalising the profiles using the instantaneous maximum $\overline {u}_m$ , the mean velocity profiles $\overline {u}/\overline {u}_m$ , as shown in figure 3(a), collapse remarkably well over the entire wall-normal extent for all $ \textit{Pr}$ cases considered. In contrast, the temperature profiles, depicted in figure 3(b), show a strong dependence on the $ \textit{Pr}$ . With increasing $ \textit{Pr}$ , the thermal field becomes progressively more confined near the wall and the profiles develop sharper gradients near the wall as reflected by a higher ${{Nu}}_\delta$ in figure 2. Such a $ \textit{Pr}$ dependence on the thermal field is not surprising, as the flows are compared at the same momentum boundary layer thickness and, thus, experience similar viscous transport across the $ \textit{Pr}$ ; whereas increasing $ \textit{Pr}$ effectively reduces the thermal diffusivity, resulting in a thinner thermal boundary layer with less efficient molecular heat transport across the boundary layer.

Figure 4. Comparison of mean profiles at matched Grashof number ${{Gr}}_\delta \approx 1.58\times 10^6$ for different Prandtl numbers plotted against wall-normal coordinate scaled by the momentum boundary layer thickness $y/\delta$ . (a) Reynolds shear stress $\overline {u^\prime v^\prime }$ normalised by its instantaneous maximum $(\overline {u^\prime v^\prime })_m$ ; and (b) temperature turbulence intensity $\overline {\theta ^\prime \theta ^\prime }$ normalised by its instantaneous maximum $(\overline {\theta ^\prime \theta ^\prime })_m$ .

A similar trend is found in the second-order statistics, as shown in figure 4, where both the Reynolds shear stress $\overline {u^\prime v^\prime }$ and the temperature variance $\overline {\theta ^\prime \theta ^\prime }$ are normalised by their instantaneous maximum to remove the $ \textit{Pr}$ dependence on the magnitude. Notably, at ${{Gr}}_\delta =1.58\times 10^6$ , the flow is in the classical (weakly) turbulent regime where the near-wall boundary layer remains laminar-like despite the presence of turbulence in the outer bulk (Ke et al. Reference Ke, Williamson, Armfield and Komiya2023). In this regime, the turbulent momentum transport near the wall is still emerging and remains relatively weak, and is particularly sensitive to the local dynamic balances: at different $ \textit{Pr}$ , the relative strengths of viscous and buoyancy forces vary, altering the onset and extent of turbulence generation in the near-wall region. Nevertheless, figure 4(a) shows excellent qualitative agreement in the Reynolds shear stress profiles across all $ \textit{Pr}$ cases throughout the boundary layer, with only minor deviations in the near-wall region ( $y/\delta \lt 0.1$ ). Figure 4(b), however, shows clear $ \textit{Pr}$ dependence on the temperature variance profiles $\overline {\theta ^\prime \theta ^\prime }$ . With increasing $ \textit{Pr}$ , the peak of the temperature variance shifts progressively closer to the wall, indicating a thinner thermal boundary layer due to reduced thermal diffusivity. Although there exist some discrepancies between our DNS data and the measurements of Tsuji & Kajitani (Reference Tsuji and Kajitani2006) and Abedin et al. (Reference Abedin, Tsuji and Hattori2009), the peak location of the temperature variance in their data still shows quantitatively consistent agreement with our $ \textit{Pr} =6$ case. The broader profiles they report may be attributed to differences in how the statistics were obtained: both studies used ensemble averaging across multiple simulations with varying initial conditions, whereas our results are based on instantaneous fields averaged over the homogeneous plane. Both Tsuji & Kajitani (Reference Tsuji and Kajitani2006) and Abedin et al. (Reference Abedin, Tsuji and Hattori2009) also observed that initial conditions can have a significant effect on the resulting flow statistics in this classical turbulent regime. This sensitivity is particularly evident for the thermal field, which exhibits steeper gradients and a thinner near-wall layer at higher $ \textit{Pr}$ , making $\overline {\theta ^\prime \theta ^\prime }$ more sensitive to variations in sampling and averaging than the momentum statistics such as $\overline {u^\prime v^\prime }$ .

3.2. Thermal boundary layer thickness

The normalised statistics in figures 3 and 4 suggest that while the momentum integral thickness $\delta$ characterises the bulk momentum structure in the turbulent NCBL, it does not capture the growth of the thermal field. This limitation is evident in figure 2, where the transition locations vary with $ \textit{Pr}$ despite matching ${{Ra}}_\delta$ . The lack of universality in ${{Nu}}_\delta$ scaling highlights that the thermal boundary layer evolves differently from the momentum integral thickness and requires a separate, thermally relevant length scale to characterise the heat transfer. For vertical NCBL flows, definitions of the thermal boundary layer thickness have been less consistent and often motivated by practical considerations rather than rigorous scaling arguments. Here, we follow a similar approach to the momentum integral thickness and define a temperature integral thickness to characterise the thermal boundary layer, given by

(3.7) \begin{equation} \delta _\theta \equiv \int _0^\infty \frac {\overline {\theta }}{\theta _w }\, \mathrm{d}y. \end{equation}

This also allows the corresponding Rayleigh number and Nusselt number to be expressed as

(3.8a,b) \begin{equation} {{Ra}}_\theta = {{Gr}}_\theta Pr = \frac {g\beta \theta _w\delta _\theta ^3}{\nu \kappa } , \quad \,{{Nu}}_\theta =\frac {q_w\delta _\theta }{\rho C_p \kappa \theta _w}={{Nu}}_\delta \frac {\delta _\theta }{\delta }. \end{equation}

Analogous to the momentum integral thickness, (3.7) integrates the mean temperature profile to capture the effective extent of thermal transport from the wall. A similar definition was adopted by Warner (Reference Warner1966), Warner & Arpaci (Reference Warner and Arpaci1968) and Talluru et al. (Reference Talluru, Pan, Patterson and Chauhan2020) to approximate the boundary layer thickness in the absence of direct momentum measurements. Figure 5(a) compares the ratio of the momentum to thermal integral thickness, $\delta /\delta _\theta$ , as a function of the Rayleigh number based on the thermal boundary layer thickness ${{Ra}}_\theta$ . In the laminar regime, $\delta _\theta$ admits an analytical solution obtained by integrating the similarity solution (2.3a ) for the temperature field,

(3.9) \begin{equation} \delta _\theta = \int _0^\infty {\text{erfc}} \left ( \eta \right )\,\mathrm{d} y = \frac {2\sqrt {\kappa t}}{\sqrt {\pi }}, \end{equation}

yielding a $ \textit{Pr}$ -dependent constant for the thickness ratio, $\delta /\delta _\theta =\mathcal{C}\sqrt {\pi }$ , and $\mathcal{C}$ is the integration constant in (3.3), as seen in figure 5(a) for ${{Ra}}_\theta \lt 3\times 10^3$ . At larger ${{Ra}}_\theta$ , the thickness ratio departs from this constant, marking the onset of laminar–turbulent transition. Unlike in figure 2 where the transition location exhibits a clear $ \textit{Pr}$ dependence, the use of ${{Ra}}_\theta$ removes this sensitivity and collapses the transition point across all $ \textit{Pr}$ cases at ${{Ra}}_\theta \approx 3.5\times 10^3$ (see the Nusselt number in figure 6). In the turbulent regime ( ${{Ra}}_\theta \gt O(10^4)$ ), all cases show a similar scaling behaviour despite differences in $ \textit{Pr}$ . Notably, the experimental measurements of Miyamoto et al. (Reference Miyamoto, Kajino, Kurima and Takanami1982) and Tsuji & Nagano (Reference Tsuji and Nagano1989), both obtained at $ \textit{Pr} =0.71$ for spatially developing flows, deviate from results at the same $ \textit{Pr}$ and instead align more closely with the DNS data for $ \textit{Pr} =4.16$ . Such a discrepancy may be due to flow development effects inherent in spatial configurations (e.g. entrainment), as well as a slight increase in ambient temperature (stratification) along the streamwise direction in experiments (Tsuji & Nagano Reference Tsuji and Nagano1989), both of which can suppress the thermal boundary layer growth and result in larger $\delta /\delta _\theta$ values in the turbulent regime. This is supported by the observation that their experimental data show qualitatively better agreement with our DNS data and the analytical solution at $ \textit{Pr} =0.71$ for ${{Ra}}_\theta \lt 10^3$ , where entrainment and ambient stratification effects are significantly weaker in the laminar regime.

Figure 5. Development of the ratio of momentum to thermal integral boundary layer thickness with (a) ${{Ra}}_\theta$ and (b) time $t$ . The (coloured) $ \textit{Pr}$ -dependent constants $\mathcal{C}$ are taken from table 3 for the respective $ \textit{Pr}$ values. In both panels (a) and (b), the scaling exponent is taken as $\xi =0.6156$ from Ke et al. (Reference Ke, Williamson, Armfield, Komiya and Norris2021) for all $ \textit{Pr}$ cases. Experiment measurements of Miyamoto et al. (Reference Miyamoto, Kajino, Kurima and Takanami1982) and Tsuji & Nagano (Reference Tsuji and Nagano1989) are taken from spatially developing flows with $ \textit{Pr} =0.71$ . The inset in panel (b) shows the compensated thickness ratio $\delta /\delta _\theta /( Pr ^{1/3}t^{1-\xi })$ development with time and the dotted lines indicate constant values of 0.52 (pastel blue) and 0.346 (grey).

Figure 6. Trends of Nusselt number ${{Nu}}_\theta$ versus Rayleigh number ${{Ra}}_\theta$ based on the thermal integral thickness. Dotted lines indicate the laminar analytical solution ${{Nu}}_\theta =2/\pi$ (black), the classical turbulent 1/3-power-law scaling (grey) and the log-law corrected 0.381-power-law scaling (violet) in the ultimate regime.

The integral definition of $\delta _\theta$ can also be modelled by the scaling arguments derived from the plume-like integral model of turbulent NCBL (Ke et al. Reference Ke, Williamson, Armfield, Komiya and Norris2021), since (3.7) effectively integrates the buoyancy across the boundary layer and therefore is expected to scale with the integral buoyancy $\hat {b}\hat {l}/\theta _w$ . In that framework, the turbulent bulk flow is dominated by the outer plume whose integral width $\hat {l}$ sets the characteristic wall-normal extent of the momentum-dominated region so that $\delta \propto \hat {l}\propto t^{\xi +1}$ ; while its characteristic buoyancy $\hat {b}$ and velocity $\hat {u}$ scales evolve with time ( $t$ ) as $\hat {b}\propto t^{\xi -1}$ and $\hat {u}\propto t^{\xi }$ , respectively. The scaling exponent $\xi$ results from the power law solution to the integral plume modelling (Ke et al. Reference Ke, Williamson, Armfield, Komiya and Norris2021, their Appendix A), and it varies only very weakly with the Rayleigh number and asymptotically approaches unity at extremely high Rayleigh numbers. As a result, the boundary layer thickness ratio, $\delta /\delta _\theta$ , follows

(3.10) \begin{equation} \frac {\delta }{\delta _\theta } \propto \frac {t^{\xi +1}}{t^{2\xi }} = t^{1-\xi } \propto {{Ra}}_\theta ^{{(1-\xi )}/{6\xi }}. \end{equation}

Evidently, given $0\lt \xi \leqslant 1$ , (3.10) suggests that the momentum boundary layer $\delta$ is always growing faster than the thermal boundary layer $\delta _\theta$ in the turbulent regime (i.e. $t\gg 0$ ). This scaling is consistent with the experimental measurements of Tsuji & Nagano (Reference Tsuji and Nagano1989) where they showed that for a spatially developing turbulent NCBL in air ( $ \textit{Pr} =0.71$ ), the thermal boundary layer defined by (3.7) grows at a slower rate than its momentum counterpart. For the current Rayleigh number range ( ${{Ra}}_\delta \lt O(10^8)$ ) considered, Ke et al. (Reference Ke, Williamson, Armfield, Komiya and Norris2021) demonstrated that $\xi$ can be treated as a constant $\xi =0.6156$ for $ \textit{Pr} =0.71$ . Here, we adopt the same value of $\xi =0.6156$ for the present cases at $ \textit{Pr} =4.16$ and $ \textit{Pr} =6$ since $\xi$ is inherently linked to the large-scale plume dynamics – as we will show later in § 3.4, the bulk plume-like region is largely $ \textit{Pr}$ -independent.

Figure 5(b) compares the development of $\delta /\delta _\theta$ across $ \textit{Pr}$ as a function of time $t$ . In the turbulent regime, all $ \textit{Pr}$ cases follow the power-law trend predicted by (3.10), albeit with a $ \textit{Pr}$ -dependent prefactor $c( Pr )$ , confirming that the growth of both $\delta$ and $\delta _\theta$ is governed by the plume-like outer dynamics. For the Prandtl numbers considered, this prefactor is empirically found to vary from $c=0.462$ to $c=0.562$ and $c=0.625$ for $ \textit{Pr} =0.71$ , $ \textit{Pr} =4.16$ and $ \textit{Pr} =6$ , respectively. In classical Rayleigh–Bénard convection, a similar thickness ratio has been analytically investigated using the Prandtl–Blasius formulation, albeit for wall-gradient based boundary thicknesses (Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010). In particular, Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010) showed that in classical Rayleigh–Bénard convection, the prefactors of the thickness ratio based on wall gradients would follow $c( Pr )=E Pr ^{1/3}$ for $ \textit{Pr} \gt 3$ , where the empirical constant $E$ is given by $E\approx 0.982^{-1}$ . Although our present vertical NCBL differs in gravity orientation and in the use of integral rather than slope thicknesses, the large- $ \textit{Pr}$ behaviour is remarkably similar. As shown in the inset of figure 5(b), the compensated thickness ratio showed that the empirical constant $E$ for our vertical NCBL may be approximated by

(3.11) \begin{equation} E = \frac {c( Pr )}{ Pr ^{1/3}}=\frac {\delta /\delta _\theta }{ Pr ^{1/3}t^{1-\xi }}\approx 0.346\pm 0.003 \quad \quad \text{for }\quad Pr = 4.16 \,\text{ and }\, Pr = 6, \end{equation}

such that the prefactor $c( Pr )$ appears consistent with the $ \textit{Pr} ^{1/3}$ dependence predicted by Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). This agreement implies that, in both Rayleigh–Bénard convection and vertical NCBL, the leading $ \textit{Pr}$ dependence (for large $ \textit{Pr}$ ) of $\delta /\delta _\theta$ is controlled by the conductive thermal sublayer, whose balance is insensitive to the precise gravity orientation, while the outer plume dynamics set the common temporal evolution exponent $\xi$ . However, for $ \textit{Pr} \lt 1$ , the two configurations differ more fundamentally in the conductive sublayer. In Rayleigh–Bénard convection, buoyancy acts normal to the plates and drives a large-scale circulation, so for $ \textit{Pr} \lt 1$ , the higher thermal diffusivity leads to a thermal layer that is thicker than the viscous one and, hence, $\delta /\delta _\theta \lt 1$ ; whereas in vertical NCBL, gravity is aligned with the wall and the buoyancy force extends throughout the region where the temperature departs from the ambient, driving the flow motion wherever the thermal field penetrates. Consequently, the momentum boundary layer is always locked to the thermal field, such that $\delta /\delta _\theta \geqslant 1$ even for $ \textit{Pr} \ll 1$ (Gebhart Reference Gebhart1988, p. 53).

Figure 6 presents the Nusselt number based on the thermal integral thickness, ${{Nu}}_\theta$ , plotted against the corresponding Rayleigh number ${{Ra}}_\theta$ . In the laminar regime, the wall heat flux follows directly from the analytical solution (2.3a ) so that $q_w = 1/\sqrt {\pi \kappa t}$ . With (3.9), it results in a $ \textit{Pr}$ -independent Nusselt number ${{Nu}}_\theta =2/\pi$ for all cases. Although the flows are initialised at different Rayleigh numbers (table 2), they all undergo transition at nearly the same ${{Ra}}_\theta$ , with the onset locations collapsing closely around ${{Ra}}_\theta \approx 3.5\times 10^3$ . It is important to note that the exact location of transition can depend on the perturbation method and amplitude, with larger amplitude disturbances typically lead to earlier transition (Ke et al. Reference Ke, Williamson, Armfield, McBain and Norris2019). However, in the present study, all cases are initialised with a fixed-amplitude broadband white noise ( $A_0=1$ ) in the temperature field, as listed in table 2, such that the velocity field will respond directly through buoyancy coupling, allowing the thermal and velocity fields to evolve self-consistently. This ensures that any observed differences in transition behaviour are attributed to intrinsic Prandtl number effects rather than external excitation. The consistent transition onset across $ \textit{Pr}$ indicates that this process is controlled by the thermal field rather than by the momentum structure – in the sense that it arises from a buoyancy-driven instability, rather than a shear-driven one. This is further supported by the recent linear stability results, which reveal that the buoyancy-driven modes dominate the instability for the NCBL flows, with shear-driven modes appearing only for $ \textit{Pr} \lt 1$ under large-scale perturbations and sufficiently high Rayleigh number (Ke et al. Reference Ke, Armfield and Williamson2024, Reference Ke, Armfield and Williamson2025). As the flow undergoes laminar–turbulent transition, ${{Nu}}_\theta$ increases quickly until it reaches a plateau around ${{Ra}}_\theta \approx 10^4$ and ${{Nu}}_\theta \approx 3$ . We note that such a plateau follows naturally from combining the behaviour of ${{Nu}}_\delta$ in figure 2 and $\delta /\delta _\theta$ in figure 5; and reflects an adjustment from the transition to the turbulent regime rather than a separate scaling regime. In this ${{Ra}}_\theta$ range, the outer plume has become turbulent and sets the extent of $\delta$ , so that $\delta _\theta /\delta$ decreases with ${{Ra}}_\delta$ as per (3.10); while the heat transfer is still transitional (see ${{Ra}}_\delta \approx 3\times 10^6$ for $ \textit{Pr} =4.16$ and $ \textit{Pr} =6$ in figure 2) and the ${{Nu}}_\delta$ – ${{Ra}}_\delta$ relation has a noticeably smaller effective exponent before it converges to the 1/3-power-law in the turbulent regime. The product of these two factors in ${{Nu}}_\theta$ leads to a partial cancellation of their ${{Ra}}_\theta$ dependences, yielding only a very weak net variation of ${{Nu}}_\theta$ , which appears as a short plateau in figure 6. In the turbulent regime, the collapsed data recover a power-law scaling in the form ${{Nu}}_\theta \propto {{Ra}}_\theta ^m$ . In figure 6, we show both $m=1/3$ , as commonly seen in the classical turbulent regime (Versteegh & Nieuwstadt Reference Versteegh and Nieuwstadt1999; Abedin et al. Reference Abedin, Tsuji and Hattori2009; Nakao et al. Reference Nakao, Hattori and Suto2017), and $m=0.38$ , which is representative of a steeper log-law corrected scaling in the ultimate turbulent regime (where $m$ increases slowly with ${Ra}$ , Ng et al. Reference Ng, Ooi, Lohse and Chung2017; Ke et al. Reference Ke, Williamson, Armfield, Komiya and Norris2021) for comparison. Over the limited ${{Ra}}_\theta$ range accessible to the present DNS, the distinction between these two power laws is subtle and an accurate determination of the exponent is hindered by turbulent fluctuations due to the instantaneous nature of the flow. The identification of the ultimate regime therefore rests primarily on the near-wall mechanisms (which will be discussed in §§ 3.3 and 3.4). Nonetheless, figure 6 shows an excellent collapse across $ \textit{Pr}$ in all flow regimes. This indicates that the thermal integral thickness $\delta _\theta$ provides a unified and physically consistent scaling to characterise heat transfer in NCBL. In vertical NCBL, the wall heat flux is controlled by a thin near-wall region in which the mean temperature falls from the wall value to the bulk, and $\delta _\theta$ is constructed to capture the thickness of this region and, hence, to provide the characteristic length scale for heat transfer and buoyancy, irrespective of whether the local dynamics are laminar or turbulent. This interpretation admits a convenient consistency check for $ \textit{Pr} =4.16$ and $ \textit{Pr} =6$ using (3.11), analogous to the large $ \textit{Pr}$ relation in Rayleigh–Bénard convection (Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010), such that ${{Ra}}_\theta ={{Ra}}_\delta ({\delta }/{\delta _\theta } )^{-3} \approx {{Ra}}_\delta E^{-3} Pr ^{-1} = {{Gr}}_\delta E^{-3}$ removes the explicit Prandtl number dependence. In addition, comparison of figures 2 and 6 shows that the Nusselt–Rayleigh relation has the same qualitative behaviour in all regimes for both definitions, indicating that using $\delta _\theta$ does not alter the underlying ${Nu}$ – ${Ra}$ scaling.

3.3. Skin friction scaling

In this subsection, we examine the near-wall momentum transport, characterised by the skin friction coefficient $C_{\kern-1.5pt f}$ given by

(3.12) \begin{equation} C_{\kern-1.5pt f} = \frac {2\tau _w}{\rho \overline {u}_m^2}, \end{equation}

where $\tau _w= \rho \nu (\partial \overline {u}/\partial y)_w$ is the wall shear stress. Here, in (3.12), we adopted a definition analogous to those used in the forced flows, in which the characteristic velocity is typically taken to be the free stream or centreline velocity that drives the flow. In our NCBL, however, the wall shear balances both the shear by which the outer bulk plume imparts momentum into the near-wall boundary layer (characterised by the local plume scale; Ke et al. Reference Ke, Williamson, Armfield, Komiya and Norris2021, Reference Ke, Williamson, Armfield and Komiya2023), and the accumulated buoyancy across the near-wall and the plume (see (4.62) of Ke et al. Reference Ke, Williamson, Armfield, Komiya and Norris2021). We use the maximum streamwise velocity, $\overline {u}_m$ , as the characteristic velocity scale representative of the plume-dominated region in the NCBL. The bulk Reynolds number is then defined as

(3.13) \begin{equation} Re_m \equiv \frac {\delta _m \overline {u}_m}{\nu }, \end{equation}

where $\delta _m$ is the location at which the maximum mean streamwise velocity occurs. Although the momentum integral thickness $\delta$ effectively collapses the mean momentum statistics across $ \textit{Pr}$ in § 3.1, here we instead employ $\delta _m$ to define the bulk Reynolds number, as it provides a clear marker separating the near-wall boundary layer structure from the outer plume-like region (Ke et al. Reference Ke, Williamson, Armfield and Komiya2023), and offers a more direct measure of the near-wall momentum layer in both numerical and experimental datasets. Similar outer scale based Reynolds numbers have also been employed in differentially heated channels (Ng et al. Reference Ng, Ooi, Lohse and Chung2015) and Rayleigh–Bénard convections (Grossmann & Lohse Reference Grossmann and Lohse2001, Reference Grossmann and Lohse2011). Figure 7 presents the skin friction coefficient $C_{\kern-1.5pt f}$ as a function of the bulk Reynolds number $Re_m$ . From figure 7, in the laminar regime, all cases show excellent collapse across different $ \textit{Pr}$ values and follow a linear scaling, $C_{\kern-1.5pt f}\propto Re_m^{-1}$ , as seen in the Prandtl–Blasius–Pohlhausen scaling (Landau & Lifshitz Reference Landau and Lifshitz1987). This laminar scaling persists up to a $ \textit{Pr}$ -dependent transition point, consistent with observations in figure 2. Such a $ \textit{Pr}$ dependence highlights that the onset of laminar–turbulent transition is driven by the thermal field and thus preventing a universal scaling based solely on momentum-based controlling parameters such as ${{Ra}}_\delta$ and $Re_m$ . Although figure 7 shows the flow completes transition at distinct $Re_m$ across $ \textit{Pr}$ , due to differences in the thermal field, it subsequently converges to a $ \textit{Pr}$ -independent linear scaling similar to that observed in the laminar regime, albeit with a slightly different prefactor. This linear scaling is suggestive of a weakly turbulent regime (Ng et al. Reference Ng, Ooi, Lohse and Chung2015; Ke et al. Reference Ke, Williamson, Armfield and Komiya2023), in which the near-wall region remains laminar-like.

Figure 7. Development of the skin friction coefficient $C_{\kern-1.5pt f}$ for the NCBL. The dashed lines indicate the $ \textit{Pr}$ -independent linear scaling $C_{\kern-1.5pt f}\propto Re_m^{-1}$ , characteristic of laminar boundary layers of Prandtl–Blasius–Pohlhausen type, in the laminar (violet) and classical turbulent (turquoise) regimes. The dotted lines indicate the log-law-type scaling (3.14) with $ \textit{Pr}$ -dependent constant $D$ . Inset shows the compensated skin friction coefficient for the log-law scaling (3.14), $C_{\kern-1.5pt f}^* = C_{\kern-1.5pt f}/[2D^2/\ln ^2(0.06Re_m)]$ , with $D$ empirically obtained to be $D=0.31$ , 0.34 and 0.35 for $ \textit{Pr} =0.71$ , 4.16 and 6, respectively.

At higher $Re_m$ , however, the flow gradually transitions into the ultimate turbulent regime, where the near-wall shear becomes strong enough to sustain turbulence (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2011; Ng et al. Reference Ng, Ooi, Lohse and Chung2015, Reference Ng, Ooi, Lohse and Chung2017; Ke et al. Reference Ke, Williamson, Armfield, Komiya and Norris2021, Reference Ke, Williamson, Armfield and Komiya2023). In this regime, the skin friction $C_{\kern-1.5pt f}$ departs from the linear Prandtl–Blasius–Pohlhausen scaling, and starts to follow a log-law type scaling that is commonly seen in fully turbulent shear layers (Kestin & Persen Reference Kestin and Persen1962; White Reference White1991),

(3.14) \begin{equation} C_{\kern-1.5pt f} = \frac {2D^2}{\ln ^2\left (0.06Re_m\right )}, \end{equation}

where $D$ is an empirical constant. For the data presented in figure 7, this coefficient $D$ shows a weak $ \textit{Pr}$ dependence, with $D=0.31$ for $ \textit{Pr} =0.71$ , $D=0.34$ for $ \textit{Pr} = 4.16$ and $D=0.35$ for $ \textit{Pr} =6$ . While the transition from the linear scaling in the classical regime to (3.14) is gradual, Ke et al. (Reference Ke, Williamson, Armfield and Komiya2023) identified the onset of this ultimate regime for $ \textit{Pr} =0.71$ by setting thresholds for Reynolds stresses and near-wall streak developments, placing it around $Re_m\approx 500{-} 550$ . Here, using additional data spanning a wider range of $ \textit{Pr}$ , we observe that the ultimate regime consistently emerges at this range (see inset of figure 7 where the log-law compensated skin friction coefficient $C_{\kern-1.5pt f}^*=1$ ) and is largely independent of $ \textit{Pr}$ considered here. This apparent $ \textit{Pr}$ insensitivity over $0.71\leqslant Pr \leqslant 6$ is consistent with the notion that the ultimate regime is driven by a shear instability, marking the point at which turbulence in the near-wall region becomes self-sustaining and, thus, primarily governed by the momentum field rather than the thermal field. This critical bulk Reynolds number is close to the critical Reynolds number $Re_c=420$ in the classical forced boundary layer theories (Schlichting Reference Schlichting1968; Landau & Lifshitz Reference Landau and Lifshitz1987), and similar concepts have been previously adapted to buoyancy-driven flows by Grossmann & Lohse (Reference Grossmann and Lohse2000) and Wells & Worster (Reference Wells and Worster2008) in the limit of asymptotically high Rayleigh number.

Figure 8. Mean temperature profiles for the turbulent NCBL at different ${{Gr}}_\delta$ and $ \textit{Pr}$ . Profiles are coloured light to dark to represent increasing ${{Gr}}_\delta$ (arrow direction). Data for $ \textit{Pr} =0.71$ are obtained from Ke et al. (Reference Ke, Williamson, Armfield and Komiya2023). The grey dotted lines indicates the near-wall profiles $\theta ^+ = Pr y^+$ ; and the black dotted straight lines indicates the temperature log-law for the flows with their respective $ \textit{Pr} _t$ values.

The emergence of this log-law-type scaling in the skin friction is intrinsically linked to the development of an established logarithmic region in the thermal field, with which the momentum field must respond accordingly to follow a buoyancy-modified log-law (Ke et al. Reference Ke, Williamson, Armfield, Norris and Komiya2020). The corresponding development of the mean velocity field in the ultimate regime has been examined in detail by Ke et al. (Reference Ke, Williamson, Armfield, Norris and Komiya2020), who showed that the friction-scaled profiles exhibit only a very short buoyancy-modified logarithmic region extending up to the local velocity maximum. The present higher- $ \textit{Pr}$ cases display similar velocity development (not shown) behaviour to that reported by Ke et al. (Reference Ke, Williamson, Armfield, Komiya and Norris2021, their figure 7) and Ke et al. (Reference Ke, Williamson, Armfield and Komiya2023, their figure 3), so here we focus on the thermal field, where the self-similar structure is much clearer. Figure 8 shows the mean temperature profiles in wall units. In figure 8, all $ \textit{Pr}$ cases display a clear conductive sublayer close to the wall ( $y^+\leqslant 5$ ), in which the mean temperature follows a Prandtl number dependent linear relation $\theta ^+= Pr y^+$ (Tsuji & Nagano Reference Tsuji and Nagano1988; Ng et al. Reference Ng, Ooi, Lohse and Chung2017). At moderate Grashof number ( ${{Gr}}_\delta \approx 1.6\times 10^6$ , as reported in figure 3), a log-law-like region is already discernible for all $ \textit{Pr}$ cases, but the vertical offset of the profiles exhibits strong sensitivity to ${{Gr}}_\delta$ (or equivalently ${{Ra}}_\theta$ ). As ${{Gr}}_\delta$ increases further, however, the bulk flow evolves towards a statistically self-similar state and the temperature log-law becomes increasingly insensitive to Grashof number; a well-defined logarithmic region is observed for each case, following the canonical log-law form,

(3.15) \begin{equation} \theta ^+ = \frac { Pr _t}{K_v} \ln (y^+)+B_\theta , \end{equation}

where $ \textit{Pr} _t$ is the turbulent Prandtl number and $K_v=0.41$ is the von Kàrmàn constant. Here, $\theta ^+ = \overline {\theta }/\theta _\tau$ and $y^+= y u_\tau /\nu$ are the mean temperature profile and wall-normal coordinate in wall units; $u_\tau =\sqrt {\tau _w/\rho }$ and $\theta _\tau = q_w/u_\tau$ are the friction velocity and temperature. The slope of the logarithmic region, $ \textit{Pr} _t/K_v$ , shows a subtle Prandtl number dependence, with fitted values of $ \textit{Pr} _t=0.85$ (see Yaglom Reference Yaglom1979; Ng et al. Reference Ng, Ooi, Lohse and Chung2017), 1.1 and 1.15 for $ \textit{Pr} =0.71$ , 4.16 and 6, respectively. The intercept of the log-law, $B_\theta$ , varies more strongly with $ \textit{Pr}$ , increasing from 7.5 to 37 and 54 for the same cases. Such an offset in $B_\theta$ may be attributed to the near-wall conductive sublayer, in which the temperature follows a $ \textit{Pr}$ -dependent linear relation. The near-wall $\theta ^+$ is therefore higher for larger $ \textit{Pr}$ at the edge of the conductive sublayer, which shifts the connection into the logarithmic region (buffer region) and produces higher intercept $B_\theta$ .

Figure 9. Development of the bulk Reynolds number $Re_m$ with the thermal Rayleigh number ${{Ra}}_\theta$ for $ \textit{Pr} =0.71$ , 4.16 and 6. The grey shaded area indicates the onset threshold for the ultimate turbulent regime at $Re_m=500\pm 50$ .

Figure 9 shows the development of the bulk Reynolds number, $Re_m$ , plotted against the thermal Rayleigh number ${{Ra}}_\theta$ . In the laminar regime, the boundary layer is self-similar, such that all relevant boundary layer thicknesses (e.g. $\delta$ , $\delta _\theta$ and $\delta _m$ ) scale with $\sqrt {\kappa t}$ , while the velocity scales linearly with $t$ (Schetz & Eichhorn Reference Schetz and Eichhorn1962; Patterson & Imberger Reference Patterson and Imberger1980). The bulk Reynolds number therefore scales with $Re_m\propto t^{3/2}$ and follows a linear scaling with ${{Ra}}_\theta$ , where the $ \textit{Pr}$ effects are reflected in the prefactors due to differences in the thermal field and resulting velocity responses. As the flow begins to transition at ${{Ra}}_\theta \approx 3.5\times 10^3$ , the DNS data deviate from the linear scaling with clear $ \textit{Pr}$ dependence inherited from the laminar regime. Such $ \textit{Pr}$ effects continue to influence the flow throughout the classical weakly turbulent regime ( ${{Ra}}_\theta \gt O(10^4)$ ), during which the flow remains sensitive to initial conditions and retains the imprint of its transitional history. As the flow approaches the critical Reynolds number $Re_m\approx 500$ , the diverged DNS data begin to collapse and the $ \textit{Pr}$ dependence tends to diminish, suggesting that the thermal forcing required to generate sufficient near-wall shear for self-sustained turbulence is nearly $ \textit{Pr}$ independent.

3.4. Ultimate regime and near-wall streaks

A key signature of the ultimate regime is the appearance of a near-wall peak in the velocity spectra, in addition to the already developed outer peak associated with large-scale plume-like structures (Ke et al. Reference Ke, Williamson, Armfield and Komiya2023). This near-wall peak signifies the emergence of coherent streaks forming in the wall-parallel plane – analogous to those observed in canonical forced-wall bounded flows, which are most prominent at $10\lt y^+\lt 30$ (Tomkins & Adrian Reference Tomkins and Adrian2003) with a typical streamwise length of $\lambda _x^+ \approx 1000$ and a spanwise spacing of $\lambda _z^+\approx 100$ – $120$ (Kline et al. Reference Kline, Reynolds, Schraub and Runstadler1967; Hutchins & Marusic Reference Hutchins and Marusic2007a , Reference Hutchins and Marusicb ). Figure 10(a–c) compares the premultiplied spectra $k_z E_{uu}$ ( $k_z = 2\pi /\lambda _z$ is the spanwise wave number and $\lambda _z$ is the spanwise wavelength) in wall units for $ \textit{Pr} =0.71$ , $4.16$ and $6$ at $Re_m\approx 600$ , which is representative of the ultimate regime identified from the skin friction scaling. Note that the coloured contours of $k_zE_{uu}$ are normalised by their instantaneous maximum to emphasise the emergence and relative prominence of the near-wall peak and therefore allows a consistent comparison across $ \textit{Pr}$ , independent of differences in absolute turbulence intensity. From figure 10, although all three $ \textit{Pr}$ cases considered show a clear near-wall spectral peak located at $y^+\approx 18$ , the most energetic spanwise wavelengths demonstrate evident $ \textit{Pr}$ dependence – increasing from $\lambda _z^+\approx 130$ for $ \textit{Pr} =0.71$ (Ke et al. Reference Ke, Williamson, Armfield and Komiya2023) to $\lambda _z^+\approx 170$ and $260$ for $ \textit{Pr} =4.16$ and $6$ , respectively. At the same wall-normal location, weaker secondary local maxima also appear in the premultiplied spanwise spectra: this is most evident in figure 10(b) for $ \textit{Pr} =4.16$ (at $\lambda _z^+\approx 110, 340$ ), while in figure 10(a, c), the corresponding local maxima are less obvious due to smaller amplitude (at $\lambda _z^+\approx 90, 100$ for $ \textit{Pr} =0.71$ and at $\lambda _z^+\approx 65, 95$ for $ \textit{Pr} =6$ ). These features are interpreted as secondary peaks within the same broadband near-wall streak band, associated with less energetic streaks of shorter spacings and with the buoyancy-modified larger-scale motions (see later in figure 11). Although the detailed mechanism by which the spectra reorganises with $ \textit{Pr}$ is beyond the scope of this study, these secondary local maxima remain subdominant for all three cases and the wavelength of the dominant near-wall spectral peak at $y^+\approx 18$ increases monotonically with $ \textit{Pr}$ over $0.71\leqslant Pr \leqslant 6$ . The corresponding wall-parallel visualisations of the instantaneous streamwise velocity field at $y^+\approx 18$ , shown in figure 10(d–f), display elongated regions of alternating high- and low-speed fluid packets, confirming the emergence of energetic spanwise motions identified in panels (a)–(c) across all $ \textit{Pr}$ cases. Notably, the outer spectral peak located at $\lambda _z^+ \approx 1000$ , associated with the large-scale plume-like motions, remains consistently positioned across all $ \textit{Pr}$ . This invariance in the outer spectral structure suggests that the energy-containing scales in the plume-like region appear insensitive to Prandtl number, consistent with the outer plumes reported by Ke et al. (Reference Ke, Williamson, Armfield and Komiya2023). These outer plume-like structures are analogous to the very-large-scale motions (VLSMs) seen in canonical wall-bounded flows, both in their spectral footprint (Lee & Moser Reference Lee and Moser2015; Wang, Wang & He Reference Wang, Wang and He2018, with spanwise wavelengths $\lambda _z^+ \sim \mathcal{O}(10^3)$ ) and in their potential to influence near-wall dynamics. In canonical turbulent boundary layers, VLSMs are known to exert a modulatory effect on the near-wall cycle, altering the amplitude and spatial arrangement of near-wall streaks through amplitude and frequency modulation (Hutchins & Marusic Reference Hutchins and Marusic2007a ; Mathis, Hutchins & Marusic Reference Mathis, Hutchins and Marusic2009). As seen from figure 10, the universality of the plume-like spectral structure across $ \textit{Pr}$ implies that any modulation exerted by these large energetic structures onto the near-wall region should also be largely independent of $ \textit{Pr}$ . This suggests that the observed broadening of near-wall streak spacings with increasing $ \textit{Pr}$ points to a fundamentally different control mechanism – while in isothermal turbulent boundary layers, the velocity streak spacing remains approximately constant at $\lambda _z^+ \approx 100$ , independent of $ \textit{Pr}$ (Li et al. Reference Li, Schlatter, Brandt and Henningson2009), it has been demonstrated that external forcing can systematically alter streak width and spacing (Bai et al. Reference Bai, Zhou, Zhang, Xu, Wang and Antonia2014, Reference Bai, Zhou, Zhang and Antonia2018). In our NCBL, the progressive broadening of near-wall streak spacing is not imposed externally but instead arises intrinsically from the Prandtl-dependent buoyancy remnant in the near-wall, which modifies the local turbulence regeneration cycle while the outer plume modulatory forcing remains statistically similar across $ \textit{Pr}$ .

Figure 10. (a–c) Premultiplied energy spectra of the velocity flucutation $k_zE_{uu}$ in the spanwise direction, normalised by its instantaneous maximum: (a) $ \textit{Pr} =0.71$ at ${{Ra}}_\theta =5.9\times 10^4$ ( $Re_m=605$ ); (b) $ \textit{Pr} =4.16$ at ${{Ra}}_\theta =6.8\times 10^4$ ( $Re_m=587$ ); and (c) $ \textit{Pr} =6$ at ${{Ra}}_\theta =6.4\times 10^4$ ( $Re_m=615$ ). The red dotted lines indicate the most energetic structure, which is representative of the streak spacing, in the near-wall region at $y^+\approx 18$ , with (a) $\lambda _z^+=130$ ; (b) $\lambda _z^+=170$ ; and (c) $\lambda _z^+=260$ . (d–f) Corresponding (truncated) visualisations of the streamwise velocity field in a wall-parallel plane at $y^+\approx 18$ , showing the emergence of near-wall streaks. The green lines indicate the most energetic spanwise wavelength at $y^+\approx 18$ in panels (a)–(c).

Figure 11. (a–c) Reynolds shear stress $-\overline {u^\prime v^\prime }$ budget for the same cases as in figure 10: (a) $ \textit{Pr} =0.71$ at ${{Ra}}_\theta =5.9\times 10^4$ ( $Re_m=605$ ); (b) $ \textit{Pr} =4.16$ at ${{Ra}}_\theta =6.8\times 10^4$ ( $Re_m=587$ ); and (c) $ \textit{Pr} =6$ at ${{Ra}}_\theta =6.4\times 10^4$ ( $Re_m=615$ ). All terms are normalised by the respective local maximum of the instantaneous shear production $\mathcal{P}_m$ . (d) Buoyancy–production ratio $\mathcal{R}$ in the wall-normal direction; the vertical black dotted line indicates $y^+=18$ and the shaded area represents the region $10\lt y^+\lt 30$ where the near-wall streaks are most prominent.

This interpretation is further supported by the Reynolds shear stress budget, which, for our unsteady parallel flow, takes the form

(3.16) \begin{equation} \frac {\partial {\big (-\overline {u^\prime v^\prime }\big )}}{\partial t} = \mathcal{P} + \mathcal{B} + \mathcal{T}_t + \mathcal{D} - \varepsilon + \mathcal{T}_p, \end{equation}

where the shear production $\mathcal{P}$ , buoyancy production $\mathcal{B}$ , turbulent transport $\mathcal{T}_t$ , viscous diffusion $\mathcal{D}$ , pseudo dissipation $\varepsilon$ and the pressure transport $\mathcal{T}_p$ are given by

(3.17a,b,c,d) \begin{equation} \mathcal{P} = \overline {v^\prime v^\prime } \frac {\partial \overline {u}}{\partial y}, \quad \mathcal{B} = -\overline {v^\prime \theta ^\prime }, \quad \mathcal{T}_t = \frac {\partial \overline {v^\prime u^\prime v^\prime }}{\partial y}, \quad \mathcal{D} = -\nu \frac {\partial ^2 \overline {u^\prime v^\prime }}{\partial y^2},\end{equation}

(3.17e,f) \begin{equation} \varepsilon =2\nu \overline {\frac {\partial v^\prime }{\partial x_i}\frac {\partial u^\prime }{\partial x_i}}, \quad \mathcal{T}_p = \left (\frac {\partial \overline {p^\prime u^\prime }}{\partial y}+\frac {\partial \overline {p^\prime v^\prime }}{\partial x}\right )-\overline {p^\prime \left (\frac {\partial u^\prime }{\partial y}+\frac {\partial v^\prime }{\partial x}\right )}. \end{equation}

Note, in (3.16), we formulate the budget for $-\overline {u^\prime v^\prime }$ , since the near-wall streaks arise and are sustained by coherent sweep ( $u^\prime \gt 0$ , $v^\prime \lt 0$ ) and ejection ( $u^\prime \lt 0$ , $v^\prime \gt 0$ ) events, both of which contribute to a net negative Reynolds shear stress $\overline {u^\prime v^\prime }\lt 0$ . Figure 11 shows the individual contributors in (3.16) to the near-wall turbulent momentum transfer $-\overline {u^\prime v^\prime }$ , where all terms are normalised by the local peak value of the shear production in the near-wall region, such that the maximum of $\mathcal{P}$ becomes unity. This normalisation enables consistent comparison across $ \textit{Pr}$ by highlighting the relative contributions from the other budget terms with respect to the dominant turbulence transport generating mechanism. As seen from figure 11, all $ \textit{Pr}$ cases exhibit qualitatively similar distributions, in which shear production ( $\mathcal{P}$ ) dominates the Reynolds stress generation in the near-wall region, with a peak located near $y^+ \approx 20$ . Meanwhile, the buoyancy term ( $\mathcal{B}$ ) is negative throughout the near-wall region for all $ \textit{Pr}$ , indicating that buoyancy consistently acts as a sink opposing the turbulent momentum transfer sustained by the sweep–ejection cycle. As $ \textit{Pr}$ increases from panel (a) to (c), the amplitude of the buoyancy sink remains appreciable relative to the other terms (e.g. $\mathcal{D}$ and $\varepsilon$ ), suggesting that buoyancy-induced suppression of shear-driven turbulence is a persistent and dynamically significant feature across all cases.

To quantify the relative importance of buoyancy suppression against shear production, we define the buoyancy–production ratio $\mathcal{R}$ based on the magnitudes of their respective contributions to the Reynolds shear stress budget,

(3.18) \begin{equation} \mathcal{R} = \left |\frac {\mathcal{B}}{\mathcal{P}}\right |. \end{equation}

Although our flow is not stratified, this ratio is analogous to a flux Richardson number, measuring the local efficiency with which buoyancy extracts energy from shear generated turbulence and suppresses the near-wall momentum transfer associated with the sweep–ejection cycle. Similar flux-Richardson-number-like measurements have also recently been used in mixed turbulent convection in a Poiseuille–Rayleigh–Bénard channel to characterise the local balance between buoyancy and shear (Xu, Li & Xi Reference Xu, Li and Xi2025), in which they revealed a shear-dominated near-wall region and a buoyancy-dominated core. Figure 11(d) plots $\mathcal{R}$ as a function of wall-normal coordinate $y^+$ for all $ \textit{Pr}$ considered. For the higher $ \textit{Pr}$ flows ( $ \textit{Pr} =4.16$ and $ \textit{Pr} =6$ ), the buoyancy-shear ratio is highest immediately adjacent to the wall ( $\mathcal{R}\sim O(1)$ ) as the reduced thermal diffusivity confines temperature gradients to the first few viscous units. This indicates that buoyancy effects in these cases are sharply localised and act most strongly within the streak-forming layer. In contrast, for $ \textit{Pr} =0.71$ , $\mathcal{R}$ is smaller ( $\mathcal{R}\approx 0.058$ ) near the wall and begins to increase beyond $y^+\approx 5$ , suggesting a further spread of temperature fluctuations that allows buoyancy to remain active over a wider region. Notably, since the momentum and thermal layers remain dynamically locked for all $ \textit{Pr}$ in vertical NCBL (§ 3.2), here, we do not expect a singular change of behaviour at $ \textit{Pr} =1$ , as would be seen in Rayleigh–Bénard convection, but rather a smooth variation from the $ \textit{Pr} =0.71$ case with slightly elevated near-wall buoyancy effects.

In addition, within the range $10 \lt y^+ \lt 30$ , where near-wall streaks are most prominent (shaded region), $\mathcal{R}$ remains finite and non-negligible for all $ \textit{Pr}$ , with systematically larger values at higher $ \textit{Pr}$ . At $y^+ = 18$ , $\mathcal{R}$ increases from approximately 0.0519 for $ \textit{Pr} = 0.71$ to 0.0736 and 0.09 for $ \textit{Pr} = 4.16$ and $6$ , respectively. This trend is consistent with the increasingly confined temperature gradients and intensified near-wall buoyancy forcing at higher $ \textit{Pr}$ , indicating that buoyancy contributes increasingly to the suppression of the small-scale turbulence and local momentum transport in the near-wall region. Comparable behaviour has also been observed in other buoyancy-influenced flows. In Rayleigh–Bénard convection with background Couette shear, Blass et al. (Reference Blass, Tabak, Verzicco, Stevens and Lohse2021) demonstrated that increasing $ \textit{Pr}$ leads to stronger near-wall gradients but reduced small-scale turbulence, promoting broader coherent structures in the near-wall region. At a fixed pair of ${Ra}$ and $Re$ , the meandering streaks systematically broaden in spanwise with increasing $ \textit{Pr}$ (see, figure 2 of Blass et al. Reference Blass, Tabak, Verzicco, Stevens and Lohse2021) – though in that configuration, the wall shear is an externally imposed control parameter rather than a flow response. Likewise, Howland et al. (Reference Howland, Ng, Verzicco and Lohse2022) examined vertical convection in a differentially heated channel and showed that larger $ \textit{Pr}$ smooth the wall-shear field and give rise to wider near-wall coherent structures (see, figure 2 of Howland et al. Reference Howland, Ng, Verzicco and Lohse2022). While the three $\mathcal{R}$ profiles in figure 11(d) converge towards similar values near $y^+\approx 30$ , this convergence does not imply that the streak spacing becomes identical at that height, for the streak structure is a result of the accumulated imprint of buoyancy–shear interactions across the entire inner layer. For high- $ \textit{Pr}$ flows, the stronger local suppression of shear-driven transport near the wall leads to a fundamentally altered turbulence structure that persists throughout the inner layer. Beyond this streak-prominent region, $ \textit{Pr}$ diverges across cases with the lower $ \textit{Pr}$ cases eventually overtaking the $ \textit{Pr} =6$ profile at $y^+\approx 30$ and $y^+\approx 45$ , respectively. However, such an increase is driven not only by buoyancy persistence (since the buoyancy extends farther out for the lower $ \textit{Pr}$ cases), but also by the weakening of shear production, which approaches zero and begins to reverse sign for the bulk flow.