3.1. Mean fields and the inverse temperature gradient

Figure 1.

Dimensionless time-averaged temperature profiles $\varTheta$ along the vertical centreline for Rayleigh numbers ranging from $\textit{Ra} = 1 \times 10^5$ to $1 \times 10^8$ . The shaded bands A, B and C highlight the three distinct evolutionary regimes of the temperature distribution, corresponding to near-onset, intermediate transition (gradient reversal) and fully developed turbulence, respectively. Representative streamline and temperature contours in the $x$ – $y$ plane demonstrate (b) the double-vortex flow structure at $\textit{Ra} = 3\times 10^5$ and (c) the single LSC at $\textit{Ra} = 3.5\times 10^5$ .

To elucidate the evolution of the mean temperature field and to enable later comparison with classical boundary-layer models (e.g. Shishkina et al. Reference Shishkina, Horn, Emran and Ching2017), we first normalize the dimensional time-averaged temperature $\langle T \rangle _t$ to a dimensionless form $\varTheta = (\langle T \rangle _t - T_c)/\varDelta$ throughout this subsection. Under this normalization, the bottom hot plate corresponds to $\varTheta = 1$ and the top cold plate to $\varTheta = 0$ . With this definition, we anchor our investigation to the theoretical critical Rayleigh number for the onset of convection in a confined cavity, denoted as $\textit{Ra}_{c,\varGamma }$ , which, according to Shishkina (Reference Shishkina2021), is approximated by

(3.1) \begin{equation} \textit{Ra}_{c,\varGamma }\approx (2\pi )^4(1+\varGamma _{x}^{-2})(1+\varGamma _{x}^{-2}/4+\varGamma _{z}^{-2}/4)=23\,378. \end{equation}

Building upon this baseline, we present the evolution of the time-averaged centreline temperature across a comprehensive range of Rayleigh numbers, from $\textit{Ra} = 10^5$ to $10^8$ . The results reveal three primary regimes in the temperature distribution, which are categorized as types A, B and C, corresponding to the highlighted distinct coloured bands in figure 1(a).

In the low- $\textit{Ra}$ regime (type-A), situated near the onset of convection, the temperature profile diverges significantly from those observed at higher parameters. As depicted by the flow structures in figure 1(b), the system is characterized by a double-vortex topology. A significant step towards understanding double-roll convection was taken by Xu et al. (Reference Xu, Xu and Xi2023), who first showed that such topologies could be actively induced by wall shear and that they dramatically enhance heat transfer. In those controlled shear cases, the formation of stable convection channels can dramatically enhance the $\textit{Nu}$ . In the present low- $\textit{Ra}$ confined RB convection, the double-vortex topology emerges naturally without external shear and exhibits an imperfect symmetry. Specifically, the two primary vortices are approximately symmetric with respect to the vertical centreline, whereas the flow is not symmetric with respect to the horizontal midplane. Two small corner vortices appear near the upper corners, while no clear counterpart corner-vortex structures are observed near the lower wall. This top–bottom asymmetry leads to unequal thermal transport in the upper and lower halves of the cell: the upper corner vortices are associated with local recirculation and redistribution of cold fluid near the top plate, whereas the transport of hot fluid near the bottom plate lacks a mirror-symmetric counterpart. Since turbulent mixing remains weak in the low- $\textit{Ra}$ type-A regime, this asymmetric mean-flow topology is retained in the time-averaged centreline temperature profile, resulting in the observed asymmetric, unidirectional development of $\varTheta (y)$ . This double-vortex structure gradually gives way to a single LSC as $\textit{Ra}$ increases (figure 1 c). Beyond this structural change, however, the flow enters an intermediate regime (type-B) once $\textit{Ra}$ exceeds $3.5 \times 10^5$ , in which the temperature profile begins to deviate from the classical boundary-layer form. A hallmark of this regime is the emergence of a persistent temperature inversion in the central core. For instance, at $\textit{Ra} = 9.06\times 10^6$ , the time-averaged centreline temperature along the $y$ -axis exhibits a clear positive gradient over an extended spatial range, even though the system as a whole maintains an overall negative temperature gradient.

Figure 2.

Spatiotemporal evolution and structural characteristics of the temperature field. (a–d) Spatiotemporal sequences of the temperature along the vertical centreline ( $x/H=0.5, z/W = 0.5$ ) for (a) $\textit{Ra}=3\times 10^5$ , (b) $3.5\times 10^5$ , (c) $7.2\times 10^6$ and (d) $4\times 10^7$ , demonstrating the transition from a simple bipartite thermal distribution at lower $\textit{Ra}$ to a multilayered stable blob structure, and eventually to near-boundary confined blobs at higher $\textit{Ra}$ . (e) A snapshot of the instantaneous temperature field on the midspan $x$ – $y$ plane ( $z/W=0.5$ ). The arrows indicate the distinct protrusion of thermal structures, highlighting the penetration of warm and cold fluids into the core region driven by the LSC.

The inverse phenomenon, also called temperature gradient reversal, has been documented in classical RB experiments using water (Chu & Goldstein Reference Chu and Goldstein1973) and air (Gille Reference Gille1967), and its existence has also been reproduced through DNS (Grötzbach Reference Grötzbach1982). While its occurrence in liquid metal systems was first reported by Segawa et al. (Reference Segawa, Naert and Sano1998), those observations were obtained under a limited set of conditions. In contrast, the present numerical investigation systematically characterizes the spatial extent of this reversal across a broad range of Rayleigh numbers. Previous studies have attributed the phenomenon to the accumulation of fluid blobs; our simulations are consistent with this interpretation and, moreover, reveal its dynamical evolution with increasing $\textit{Ra}$ . Spatiotemporal temperature sequences of the centreline confirm that at $\textit{Ra} = 3 \times 10^5$ , the system features a simple flow structure where warm fluid predominantly occupies the lower-half of the cavity and cold fluid occupies the upper-half, as shown in figure 2(a), perfectly corresponding to the thermal signature of the double-vortex structure. However, when $\textit{Ra}$ exceeds $3.5 \times 10^5$ , the flow undergoes a transition into a distinct multilayered configuration.

In this multilayered regime, the LSC advects fluid released from the thermal boundary layers across the bulk. Owing to the strong thermal diffusion of low- $\textit{Pr}$ fluids, these transported fluid parcels would rapidly lose their temperature contrast if they were not carried into the interior sufficiently effectively. In the present type-B regime, however, the LSC can bring warm and cold parcels away from the plates while they still retain a finite thermal signature. This process promotes their accumulation in opposite regions of the bulk rather than directly at the plates: warm blobs preferentially occupy the upper bulk region, whereas cold blobs accumulate in the lower bulk region (as shown in figure 2 b). The finite $\textit{Ra}$ -range of the type-B state indicates that this blob-like accumulation is not governed by the LSC strength alone. With a further increase in $\textit{Ra}$ , turbulent mixing in the bulk is progressively strengthened, which weakens the temperature contrast of the accumulated blobs and reduces the magnitude of the local inverse gradient. Consequently, the blob-like structures are not immediately eliminated, but become increasingly compressed towards the thermal boundary layers. This transition is already visible at intermediate $\textit{Ra}$ and becomes most pronounced near the upper boundary of the type-B regime, where the structures are eventually localized near the thermal boundary layers, as shown in figures 2(c) and 2(d). In the subsequent high- $\textit{Ra}$ regime (type-C), the time-averaged temperature profile exhibits a characteristic overshoot effect, a phenomenon similarly observed in high- $\textit{Ra}$ convection (Ravi et al. Reference Ravi, Henkes and Hoogendoorn1994). This core structure can also be seen from the mean temperature field in the midspan $x$ – $y$ plane, shown in figure 2(e). The superimposed arrows indicate the protrusion of the warm and cold regions into the central core, consistent with the LSC-assisted penetration of thermal fluid from the boundary layers into the bulk. A further constraint comes from the rate at which temperature reversal can be produced and attenuated near the plates. In the universal framework, the present parameter range falls in a boundary-layer-dominated regime for both thermal and kinetic-energy dissipation (Grossmann & Lohse Reference Grossmann and Lohse2000; Xu et al. Reference Xu, Shi and Xi2019). The relevant time scale is therefore taken as the diffusion time across the thermal boundary layer. Following Grossmann & Lohse (Reference Grossmann and Lohse1993), the thermal boundary-layer thickness is estimated as $\delta _T=H/(2\textit{Nu})$ , giving $\tau _{\delta _T}=\delta _T^2/\alpha =H^2/(4\alpha \textit{Nu}^2)$ . This time scale is compared in figure 3(a) with the LSC turnover time ( $\tau _{\textit{LSC}}=H/U_{\textit{LSC}}$ ), where $U_{\textit{LSC}}=\lt u \gt_{t,\textit{max}}$ , and the free-fall time $\tau _{\!f}=H/U_{\!f}$ . To emphasize the relative strength of the three processes, figure 3(b) further shows the corresponding ratios $\varPi _1=\tau _{\delta _T}/\tau _{\textit{LSC}}$ and $\varPi _2=\tau _{\delta _T}/\tau _{\!f}$ .

Figure 3.

Time scale balance for the type-B regime. (a) Characteristic times for LSC turnover, free-fall response and thermal diffusion across the thermal boundary layer as functions of $\textit{Ra}$ . (b) Ratios $\varPi _1=\tau _{\delta _T}/\tau _{\textit{LSC}}$ and $\varPi _2=\tau _{\delta _T}/\tau _{\!f}$ . The horizontal dashed line denotes unity, and the vertical dashed lines mark the two ordering changes at $\textit{Ra}_{c1}\approx 3.33\times 10^5$ and $\textit{Ra}_{c2}\approx 4.21\times 10^7$ .

Figure 3. Long description

Two line graphs depict the relationship between characteristic times and Rayleigh numbers in turbulent Rayleigh-Bénard convection. Panel A: A line graph shows the characteristic times for large-scale circulation (LSC) turnover, free-fall response, and thermal diffusion across the thermal boundary layer as functions of the Rayleigh number (Ra). The x-axis represents the Rayleigh number (Ra) on a logarithmic scale ranging from 10^5 to 10^8. The y-axis represents the characteristic time (τ) on a logarithmic scale ranging from 10^0 to 10^2. Three data series are plotted: τ_f (red circles), τ_LSC (orange triangles), and τ_δT (blue squares). Vertical dashed lines mark the ordering changes at Ra_c1 = 3.33 x 10^5 and Ra_c2 = 4.21 x 10^7. Panel B: A line graph shows the ratios Π_1 = τ_δT/τ_LSC (blue circles) and Π_2 = τ_δT/τ_f (orange triangles) as functions of the Rayleigh number (Ra). The x-axis represents the Rayleigh number (Ra) on a logarithmic scale ranging from 10^5 to 10^8. The y-axis represents the ratio (Π) on a logarithmic scale ranging from 10^-1 to 10^1. A horizontal dashed line denotes unity. Vertical dashed lines mark the ordering changes at Ra_c1 = 3.33 x 10^5 and Ra_c2 = 4.21 x 10^7.

At low $\textit{Ra}$ , ( $\varPi _1\gt 1$ ), indicating that thermal diffusion across the relatively thick thermal boundary layer is slower than one LSC turnover. The flow topology in this range is still characterized by a double-vortex state, and the temperature field is organized primarily as a broad lower-hotter and upper-colder distribution. Under these conditions, warm and cold parcels released from the plates do not separate into persistent accumulations in the bulk; instead, their motion remains tied to the large roll-scale recirculation.

The high- $\textit{Ra}$ limit is constrained in a different way. As $\textit{Ra}$ increases, $\varPi _2$ decreases and eventually falls below unity. The thermal boundary layer then becomes sufficiently thin that diffusion over its thickness, together with the enhanced turbulent mixing in the bulk, rapidly smooths out the temperature difference carried by the transported fluid parcels, making stable bulk accumulation unsustainable. The large-scale motion therefore contributes increasingly to core homogenization, and the blob-like structures become confined to the near-wall regions.

The type-B regime lies between these two limiting behaviours, where the ordering $\varPi _1\lt 1\lt \varPi _2$ , or equivalently $1/\tau _{\textit{LSC}}\lt 1/\tau _{\delta _T}\lt 1/\tau _{\!f}$ , is satisfied. In this interval, warm and cold fluid parcels can be released from the thermal boundary layers and carried into the bulk within one LSC turnover, while diffusion and turbulent mixing are still not strong enough to rapidly homogenize their temperatures. These parcels can therefore remain in the bulk for a sufficiently long time and accumulate in opposite regions, leading to stable blob-like structures and the associated mean-gradient reversal. Ordering changes at $\textit{Ra}_{c1}\approx 3.33\times 10^5$ and $\textit{Ra}_{c2}\approx 4.21\times 10^7$ bracket the observed type-B range, where stable blob-like structures coexist with the mean-gradient reversal.

Despite these observations, classical temperature models (e.g. Shishkina et al. Reference Shishkina, Horn, Emran and Ching2017), which are primarily optimized for isothermal bulk states, typically do not account for the temperature gradient reversal. To address this gap, we propose a semiempirical model for the mean temperature profile. While classical boundary layer theories (e.g. Prandtl–Blasius, Falkner–Skan solutions) and recent models provide excellent predictions for fully developed turbulent convection ( $\textit{Pr} \gtrsim 1$ and high $\textit{Ra}$ , characterized by an isothermal bulk), they encounter new challenges when applied to transitional or low- $\textit{Pr}$ flows where the bulk is non-isothermal. Under these specific conditions, the dimensionless mean temperature profile $\varTheta (y^*)$ does not adhere strictly to the traditional ‘boundary layer plus isothermal core’ paradigm. The DNS data frequently reveals a non-negligible temperature gradient in the bulk, and even a ‘bulk reversal’ in the central core, which traditional models are constrained by the assumption $\partial \varTheta /\partial y^*|_{\textit{bulk}} \approx 0$ , which precludes the description of a temperature reversal. Following established theoretical frameworks (Shishkina & Thess Reference Shishkina and Thess2009; Wang et al. Reference Wang, He and Tong2016; Shishkina et al. Reference Shishkina, Horn, Emran and Ching2017; Xu et al. Reference Xu, Wang, He, Wang, Schumacher, Huang and Tong2021), our analysis here shifts from local centreline probes to the horizontally area-averaged and time-averaged dimensionless temperature profiles $\langle \varTheta \rangle _{A}(y^*)$ in the $x$ – $z$ plane.

To precisely characterize this global vertical transition, we introduce a phenomenological hybrid model based on the principle of superposition, which integrates the near-wall boundary layer decay with a linear gradient correction in the bulk. To align the present dimensional coordinate $y$ with the dimensionless coordinate used in the cited references, we define $y^* = y/H$ . We first define a base function $\varPhi (y^*)$ to characterize the thermal boundary layer decay, adopting the generalized exponential form suggested by Shishkina & Thess (Reference Shishkina and Thess2009),

(3.2) \begin{align} \varPhi (y^* ;a,b,c)=\exp \big (\!-ay^* +by^{*2}+cy^{*3}\big ),\quad y^*\in (0,1), \end{align}

where $a$ , $b$ and $c$ are shape parameters controlling the boundary layer thickness and curvature. The exponent includes a cubic term ( $y^{*3}$ ) to capture the near-wall scaling of the eddy thermal diffusivity as derived by Shishkina & Thess (Reference Shishkina and Thess2009), Shishkina et al. (Reference Shishkina, Horn, Wagner and Ching2015), while efficiently avoiding the computational cost of numerical integration. To guarantee antisymmetry about the central plane ( $y^{*} = 0.5$ ) and satisfy the boundary conditions $\varTheta (0)=1$ and $\varTheta (1)=0$ , an antisymmetric base profile $\varTheta _{\textit{base}}$ is constructed:

(3.3) \begin{align} \varTheta _{\textit{base}}(y^*)=\varPhi (y^{*})-\varPhi (1-y^{*}). \end{align}

However, because $\varTheta _{\textit{base}}$ inherently decays to zero in the bulk (implying $\varTheta \rightarrow 0.5$ ), it cannot describe a non-isothermal core. To accurately capture the bulk gradients observed in the DNS data, we introduce a linear correction term $L(y^{*}) = 1 - 2y^{*}$ , weighted by a fitted bulk-gradient correction factor, $S$ . The final formulation of the hybrid model is thus given by

(3.4) \begin{align} \varTheta (y^*) = \frac {1}{2} + \left (\frac {1}{2} - S\right )[\varPhi (y^{*}) - \varPhi (1-y^{*})] + S(1-2y^{*}). \end{align}

Figure 4.

Validation of the proposed hybrid model against DNS data and classical theory. (a–c) Quantitative comparisons of the area-averaged dimensionless mean temperature profiles ( $\langle \varTheta \rangle _A$ ) among the present DNS data (blue circles), the strict integral solution by Shishkina et al. (Reference Shishkina, Horn, Emran and Ching2017) (red solid line) and our proposed hybrid S-Fit model (blue solid line) at (a) $\textit{Ra}=3.5\times 10^5$ , (b) $2.57\times 10^6$ and (c) $4\times 10^7$ . The inset in (c) provides a magnified view of the central core, illustrating the S-Fit model’s capability to capture the temperature gradient reversal. (d) The evolution of the absolute bulk gradient factor, $|S|$ , as a function of the $\textit{Ra}$ . The inset in (d) further shows the variation of the fitted parameters $a$ , $b$ and $c$ in (3.2) with $\textit{Ra}$ .

In this framework, $S$ is used as a fitted phenomenological indicator, rather than as a first-principles order parameter, to characterize the non-isothermal structure of the bulk. Specifically, $S$ is introduced to represent the approximately linear correction associated with the inverse-gradient core region and is determined, together with $a$ , $b$ and $c$ , through a least-squares fit to the horizontally and temporally averaged temperature profile over $0 \leqslant y^* \leqslant 1$ . When $S \to 0$ , the model correctly degenerates into the classical isothermal bulk profile ( $\varTheta \approx 0.5$ at $y^* = 0.5$ ), which is applicable to fully developed turbulence at high $\textit{Ra}$ . Conversely, when $S \ne 0$ , the term $S(1-2y^*)$ introduces a linear slope in the core.

We validated the efficacy of this hybrid model utilizing DNS data spanning $\textit{Ra} = 3.5 \times 10^5$ to $10^8$ . A detailed comparative analysis, presented in figure 4(a–c), juxtaposes the DNS data (blue circles) with predictions from both models. It demonstrates that while the strict integral solution of Shishkina et al. (Reference Shishkina, Horn, Emran and Ching2017) (red solid line) accurately captures near-wall gradient characteristics based on local similarity, its theoretical assumption of an isothermal bulk ( $\varTheta (\infty ) = 0.5$ ) naturally leads to a different behaviour in the transition region between the boundary layer and the bulk. Consequently, this framework leaves the characterization of the temperature gradient reversal as an open area for further extension. In contrast, our proposed hybrid model (blue solid line), via the introduction of the $S$ factor, precisely delineates the complete profile transition from the wall to the core (as shown in figure 4 c).

To systematically evaluate the overall strength of this non-isothermal core regardless of its direction, we examine the absolute magnitude of the bulk gradient factor, $|S|$ . As illustrated in figure 4(d), $|S|$ progressively decreases with increasing $\textit{Ra}$ , indicating a transition towards fully developed turbulence in which energetic mixing homogenizes the core. For completeness, the corresponding fitted parameters $a$ , $b$ and $c$ in (3.2), determined from the S-Fit procedure at each $\textit{Ra}$ , are also presented in the inset of figure 4(d). When $\textit{Ra}$ approaches $10^8$ , $|S|$ asymptotes to zero, physically reflecting the regression of the core gradient back to an isothermal state. This confirms that the hybrid model accurately degenerates to the classical near-isothermal bulk profile under such highly turbulent conditions. Although $|S|$ is a fitted parameter and does not directly represent the local heat flux, it nevertheless serves as a phenomenological indicator of the strength of the non-isothermal core within the present hybrid model.

3.2. Spatial characterization of the bimodal temperature fluctuations

To accurately quantify the absolute intensity of turbulent thermal mixing, the temperature standard deviation $\sigma _T = \sqrt {\langle (T - \langle T \rangle _t)^2 \rangle _t}$ is evaluated directly using the dimensional temperature field throughout the remainder of this paper. This rigorous representation reflects the true scale of thermal fluctuations driven by the flow dynamics. The vertical profile of $\sigma _T$ reveals a distinct, anomalous bimodal structure across the entire $\textit{Ra}$ range investigated. As illustrated in figure 5, this dual-peak feature becomes progressively more pronounced as $\textit{Ra}$ increases from $2.57 \times 10^6$ to $1.14 \times 10^7$ . This bimodal profile stands in stark contrast to the classical unimodal distribution typically observed in standard RB convection. Existing theoretical frameworks – such as the ‘multilayered boundary layer structure’ proposed for horizontal convection (Yan & He Reference Yan and He2023), the inference that strong shear effects in the viscous boundary layer cause the thermal boundary layer to split into sublayers (Xu et al. Reference Xu, Wang, He, Wang, Schumacher, Huang and Tong2021), or the twin-peak horizontal velocity fluctuations induced by wall jets in slim-box convection (Zhou & Chen Reference Zhou and Chen2021) – cannot adequately explain the secondary peak observed near the core region in the present study. Therefore, new physical perspectives are required to elucidate the underlying turbulent heat transport mechanisms.

Figure 5.

Vertical profiles of the temperature standard deviation, $\sigma _T$ , along the centreline for various Rayleigh numbers ranging from $\textit{Ra}=2.57\times 10^6$ to $1.14\times 10^7$ . The profiles highlight the progressive emergence and development of the anomalous bimodal structure as $\textit{Ra}$ increases within this regime.

Figure 6.

The physical mechanism associated with the bimodal temperature fluctuations. (a) Three-dimensional visualization of the flow topology, illustrating the spiral precession of corner vortices forced by spanwise confinement. (b) Conditional response of the off-wall peak regions to the regional vortical activity. The samples are sorted according to $C_x(t)$ and divided into ten equal-population bins, indexed by $k=1,\ldots ,10$ . The plotted quantity is the normalized conditional response $R_q(k)=\langle q\rangle _k/\langle q\rangle _1$ , where $q=W_z$ or $S_T$ . Here $W_z$ and $S_T$ denote the regional r.m.s. amplitudes of the spanwise velocity fluctuation and temperature fluctuation, respectively. (c–f) Vertical centreline profiles comparing the spanwise velocity fluctuation, $\langle w^{\prime 2} \rangle _t$ (blue solid spheres), and the temperature standard deviation, $\sigma _T$ (yellow solid spheres), at (c) $\textit{Ra} = 4.54 \times 10^6$ , (d) $7.2 \times 10^6$ , (e) $9.06 \times 10^6$ and (f) $1.14 \times 10^7$ , respectively.

The dynamic behaviour of corner vortices, known to be crucial for flow reversals and mode transitions in RB convection (Xi & Xia Reference Xi and Xia2007; Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010; Xu et al. Reference Xu, Chen, Wang and Xi2020), provides a plausible flow–structural mechanism for the bimodal temperature fluctuation profile. Spanwise geometric confinement imposed by the front and back walls suppresses the free spanwise development of these corner vortices (figure 6 a). Upon detaching, they are forced to intrude horizontally into the central region in a coherent, spiral precession within the $x$ – $y$ plane, generating intense spanwise ( $z$ -direction) velocity fluctuations within the plume-shearing and impinging zones. To examine whether this vortical activity is associated with the enhanced spanwise velocity disturbance and temperature fluctuation, we performed a conditional-average analysis within the off-wall $\sigma _T$ peak regions. These regions are used as peak-containing sampling bands because they cover the elevated- $\sigma _T$ zones where the anomalous fluctuation peaks occur; their precise bounds are defined in Appendix A.

At each instant, we evaluate the regional vortical activity $C_x(t)$ , the regional root-mean-squared (r.m.s.) spanwise velocity-fluctuation amplitude $W_z(t)$ and the regional r.m.s. temperature-fluctuation amplitude $S_T(t)$ within the selected bands. The samples are sorted by $C_x(t)$ and divided into 10 equal-population bins, indexed by $k=1,\ldots ,10$ . Figure 6(b) shows the normalized conditional response $R_q(k)=\langle q\rangle _k/\langle q\rangle _1$ at $\textit{Ra} = 7.2 \times 10^6$ , where $q=W_z$ or $S_T$ , with the weakest- $C_x$ bin used as the reference. Both responses increase with $k$ , indicating that stronger regional vortical activity is accompanied by simultaneous increases in the regional spanwise-velocity and temperature-fluctuation amplitudes. The corresponding enhancement amplitudes, confidence intervals and robustness checks against the band selection and the conditioning variable are also reported there.

Figure 6(c–f) further compare the centreline profiles of the spanwise velocity variance $\langle w^{\prime 2}\rangle _t$ and the temperature standard deviation $\sigma _T$ at different $\textit{Ra}$ . The two quantities display the same bimodal organization, with near-core extrema occurring at nearly identical vertical locations. This spatial agreement is consistent with the conditional-average result in figure 6(b) and supports an association between the off-wall temperature-fluctuation peaks and localized spanwise velocity disturbances. Although $\langle w^{\prime 2}\rangle _t$ and $\sigma _T$ have similar large-scale profile shapes, their near-wall peak locations are not exactly identical. The peak of $\langle w^{\prime 2}\rangle _t$ lies closer to the wall than that of $\sigma _T$ , which is consistent with the thinner velocity boundary layer in low- $\textit{Pr}$ liquid metals. The velocity profiles also contain more small-scale variations, reflecting the stronger sensitivity of $\langle w^{\prime 2}\rangle _t$ to local turbulent disturbances.

Figure 7.

Statistical characterization of the dimensional temperature fluctuations at the two peak locations for $\textit{Ra} = 4.54\times 10^6$ . (a) Temperature time series at the near-wall first peak ( $y = \delta _{\sigma -1}$ ). (b) The PDF of the normalized temperature fluctuations at $\delta _{\sigma -1} = 0.04H$ , showing a near-Gaussian distribution. (c) Corresponding time series at the near-core second peak ( $y = \delta _{\sigma -2}$ ). (d) The PDF evaluated at $\delta _{\sigma -2} = 0.15H$ , displaying a strongly non-Gaussian distribution with pronounced positive skewness indicative of intermittent plume arrival; the inset shows the time-averaged temperature profile along the centreline at this parameter, with the red dot indicating the position and mean temperature of the shoulder.

Figure 7. Long description

Panel A: A line graph shows the temperature time series at the near-wall first peak. The x-axis is labeled t/tf and the y-axis is labeled T. The data appears to fluctuate around a mean value. Panel B: A probability density function (PDF) graph displays the normalized temperature fluctuations at the near-wall first peak. The x-axis is labeled (T - <T>ρ)/σT and the y-axis is labeled p.d.f. The data points are shown as blue circles, and a Gaussian fit is represented by an orange line. Panel C: Another line graph shows the temperature time series at the near-core second peak. The x-axis is labeled t/tf and the y-axis is labeled T. The data shows more pronounced fluctuations compared to Panel A. Panel D: A PDF graph displays the normalized temperature fluctuations at the near-core second peak. The x-axis is labeled (T - <T>ρ)/σT and the y-axis is labeled p.d.f. The data points are shown as blue circles, and a Gaussian fit is represented by an orange line. The inset shows the time-averaged temperature profile along the centerline, with the red dot indicating the position and mean temperature of the shoulder.

To further decode the distinct physical nature of the two $\sigma _T$ peaks, we examine the temperature time series and their corresponding PDFs at these specific spatial locations. Figures 7(a) and 7(c) display the temporal evolution of temperature at the near-wall peak ( $\delta _{\sigma -1} = 0.04H$ ) and the near-core peak ( $\delta _{\sigma -2} = 0.15H$ ), respectively. The signal at $\delta _{\sigma -1}$ is characterized by high-frequency fluctuations, which is highly consistent with the chaotic nature of boundary-layer instabilities (He & Xia Reference He and Xia2019; Kim et al. Reference Kim, Schindler, Vogt and Eckert2024). In contrast, the signal at $\delta _{\sigma -2}$ is dominated by lower-frequency, large-amplitude intermittent events.

The PDFs further clarify this distinction. At $\delta _{\sigma -1}$ (figure 7 b), the PDF is approximately Gaussian with only slight positive skewness, a typical signature of the intermittent detachment of thermal plumes from the boundary layer, as documented by He & Xia (Reference He and Xia2019) and Kim et al. (Reference Kim, Schindler, Vogt and Eckert2024) for near-wall regions. At $\delta _{\sigma -2}$ (figure 7 d), by contrast, the PDF exhibits a strongly non-Gaussian profile with pronounced positive skewness, consistent with the observations of He & Xia (Reference He and Xia2019) in plume-abundant regions: the intermittent arrival of hot fluid parcels from the boundary layer skews the local temperature distribution towards higher values. A simple temperature reconstruction helps identify the likely origin of these parcels. The pronounced ‘shoulder’ in the positive tail of the second-peak PDF (figure 7 d) corresponds to a dimensionless temperature fluctuation value $(T-\langle T \rangle _t)/\sigma _T$ of roughly 2.0. Multiplying this value by the standard deviation $\sigma _{T}$ and subsequently adding the local mean temperature $\langle T\rangle _{t}$ at the second peak yields an absolute temperature of approximately 3.81. This value closely matches the time-averaged temperature at the first peak, $\langle T \rangle _t \approx 3.87$ . Thus, the hot fluid parcels responsible for the large variance at the second peak do not originate locally; instead, they are energetic plumes extracted from the boundary layer and transported deep into the core region by the spanwise motion of the corner vortices.

To fully characterize the spatial extent of this bimodal structure, we extend our analysis to the full-field DNS data. Figure 8(a–c) present the contour maps of $\sigma _T$ on $x$ – $y$ cross-sections at fixed spanwise locations ( $z/W = 0.5, 0.7$ and $0.9$ ). It is important to clarify that the apparent tilt of the domains in these visualizations is strictly a graphical choice for optimal viewing of the internal planes; there is no physical alteration to the set-up, and gravity remains oriented exactly along the vertical $y$ -axis. Near the spanwise walls (e.g. $z/W \lt 0.2$ and $z/W \gt 0.8$ ), the strong suppression of turbulent fluctuations by the viscous sublayer leads to severe temperature fluctuation attenuation, resulting in a more regular, constrained contour structure. However, moving towards the spanwise midplane ( $z/W = 0.5$ ), the bimodal structure persists over a broad spatial region.

Figure 8.

Spatial distribution and 3-D evolution of the bimodal $\sigma _T$ structure at $\textit{Ra}=7.2\times 10^6$ . (a–c) Contour plots of $\sigma _T$ on vertical $x$ – $y$ planes at different spanwise depths: $z = 0.5W$ , $0.7W$ and $0.9W$ , respectively. Solid and dashed wireframes in panel (a) indicate the positions of the cross-sectional planes within the 3-D domain. In panels (a)–(c), the view is rotated only to show the internal cross-sections; the physical coordinate system is unchanged, and gravity acts along the vertical $y$ -direction. (d) Vertical profiles of $\sigma _T$ evaluated at various streamwise locations ( $x/L_x$ ) along the mid-depth plane, illustrating the horizontal migration of the bimodal feature. (e) Vertical profiles of $\sigma _T$ at various spanwise locations ( $z/W$ ) along the vertical centreline, highlighting the influence of the spanwise end-walls.

The spatial evolution of this structure along the streamwise ( $x$ ) and spanwise ( $z$ ) directions is intimately linked to its physical origin. Figure 8(d) plots the vertical profiles of $\sigma _T$ at various $x$ -positions on the spanwise midplane. As the interrogation line moves horizontally from the plume-impinging sidewall ( $x \sim 0$ ) to the plume-ejecting sidewall ( $x \sim L_x$ ), the $\sigma _T$ profile undergoes a continuous topological transformation: it begins as an asymmetric bimodal shape where the peak is strictly confined to the lower-half-axis ( $y/H \lt 0.5$ ), transitions to a symmetric bimodal profile in the absolute centre, and eventually mirrors its initial asymmetric shape extending into the upper-half-axis ( $y/H \gt 0.5$ ) before vanishing. Crucially, the spatial presence of the bimodal structure overlaps almost perfectly with the plume-impact and shear-dominated regions, but disappears entirely in the pure plume-ejection region. This spatial localization physically manifests the life cycle of the corner vortices: they form and develop in the impinging zone, migrate horizontally while dissipating kinetic energy due to viscosity, and eventually decay before reaching the ejection region, causing the bimodal structure to disintegrate. Furthermore, the evolution along the spanwise direction $z$ (figure 8 e) confirms the multiscale regulatory effect of the spanwise sidewalls. When the sampling line moves from the spanwise endwall at $z=0$ towards the spanwise midplane ( $z/W=0.5$ ), the $\sigma _T$ profile first appears as a unimodal curve because the no-slip sidewall strongly attenuates local temperature fluctuations. It then develops a weak off-wall shoulder and eventually becomes a clear bimodal profile near $z=0.5W$ . The same transition occurs in reverse when the sampling line approaches the opposite spanwise endwall at $z=W$ . This spanwise evolution indicates that the major plume-induced temperature fluctuations are concentrated around the spanwise midplane and decay gradually towards the sidewalls. Near $z=0.5W$ , the plume-induced temperature fluctuations associated with the corner-vortex-driven transport remain sufficiently strong to form two distinct off-wall peaks in $\sigma _T$ . Closer to the sidewalls, the no-slip boundary condition damps the spanwise velocity fluctuations and weakens the redistribution of thermal anomalies. The secondary near-core peak therefore decreases gradually and finally merges into a conventional unimodal profile dominated by near-wall thermal fluctuations.

The numerical investigation is extended to higher Rayleigh numbers to examine the ultimate fate of the bimodal phenomenon. As shown in figure 9(a), increasing $\textit{Ra}$ elevates the absolute magnitude of $\sigma _T$ substantially, with fluctuations at $\textit{Ra} = 1.0 \times 10^8$ drastically exceeding those at lower values. Visually, however, the bimodal nature disappears at this extreme parameter, where the near-core profile degenerates into a mere inflection point.

To systematically quantify how $\textit{Ra}$ affects the structural topology of the bimodal profile, a min–max normalized temperature standard deviation is introduced,

(3.5) \begin{align} \sigma _{T}^{*}=\frac {\sigma _T-\sigma _{T,\textit{min}}}{\sigma _{T,\textit{max}}-\sigma _{T,\textit{min}}}, \end{align}

where $\sigma _{T,\textit{max}}$ and $\sigma _{T,\textit{min}}$ are the local maximum and minimum of the standard deviation along the profile. As figure 9(b) clearly shows, increasing $\textit{Ra}$ (downward arrow) causes the secondary near-core peak to collapse progressively. At $\textit{Ra} \approx 4 \times 10^7$ , the secondary peak region flattens completely, and by $\textit{Ra} = 1.0 \times 10^8$ , the bimodal characteristic is fully suppressed, reverting to a quasiunimodal state.

Figure 9.

Effect of high Rayleigh numbers on the bimodal temperature fluctuations. (a) Vertical profiles of the absolute $\sigma _T$ , for $\textit{Ra}$ ranging from $7.2 \times 10^6$ up to $1.0 \times 10^8$ . (b) Corresponding profiles of the min–max normalized temperature standard deviation, $\sigma _T^*$ .

Figure 9. Long description

Two line graphs depict vertical profiles of temperature fluctuations at varying Rayleigh numbers. Panel A: The line graph shows vertical profiles of the absolute temperature standard deviation (σT) for Rayleigh numbers (Ra) ranging from 7.2 x 10^6 to 1.0 x 10^8. The x-axis represents the normalized height (y/H) from 0 to 1, and the y-axis represents the temperature standard deviation (σT) from 0 to 14. The graph includes multiple lines, each representing a different Rayleigh number, with an arrow indicating increasing Ra. Panel B: The line graph shows vertical profiles of the min-max normalized temperature standard deviation (σT*) for the same range of Rayleigh numbers. The x-axis represents the normalized height (y/H) from 0 to 1, and the y-axis represents the normalized temperature standard deviation (σT*) from 0 to 1.2. Similar to Panel A, multiple lines represent different Rayleigh numbers, with an arrow indicating increasing Ra.

This topological collapse corroborates the proposed mechanism based on corner-vortex interference. In a Q2D confined cavity, increasing $\textit{Ra}$ causes the primary LSC to strengthen and expand, penetrating deeper into the near-wall regions while the corner rolls become spatially compressed and kinematically weakened. A relatively weak LSC at lower $\textit{Ra}$ grants the corner vortices sufficient spatial freedom to thrive and execute the precessing motion that drives the secondary $\sigma _T$ peak. At sufficiently high $\textit{Ra}$ , by contrast, the overwhelmingly powerful primary roll compresses the corner vortices into very small spaces and restricts their dynamical influence, cutting off the spanwise velocity fluctuations that transport plumes to the core and causing the bimodal temperature fluctuation profile to vanish. This dynamic competition between the primary and corner rolls with increasing $\textit{Ra}$ is well documented in both experimental and numerical RB studies (Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010; Xu et al. Reference Xu, Wang, He, Wang, Schumacher, Huang and Tong2021).

3.3. Role of spanwise confinement in regulating the bimodal temperature fluctuations

Geometric confinement from the spanwise walls underpins the proposed mechanism; relaxing this confinement should accordingly attenuate or eliminate the bimodal structure. To test this geometric sensitivity, we perform simulations at varying spanwise aspect ratios while fixing $\textit{Ra} = 7.2 \times 10^6$ . Figure 10(a–d) present the 3-D mean velocity streamlines, coloured by the $\sigma _T$ , to illustrate the structural evolution of the flow. At the smallest aspect ratio, $\varGamma _z = 0.2$ (figure 10 a), the corner rolls experience severe spanwise constraint. Such strong confinement forces their motion into a coherent, spiral-like precession confined primarily to the streamwise-vertical $x$ – $y$ plane, which in turn generates strong spanwise ( $z$ -direction) velocity fluctuations. As $\varGamma _z$ is progressively increased through 0.4, 0.6 and 0.8 (figure 10 b–d), this geometric confinement gradually relaxes. The corner rolls gain more spatial freedom for spanwise development and their previously ordered, tightly constrained spiral motion gives way to more 3-D turbulent mixing with enhanced directional selectivity: the flow no longer stays locked into a specific precession mode but instead can develop more freely along multiple directions. The transition described above is schematically summarized in figure 10(e). Under strong confinement, the corner vortex remains restricted by the spanwise sidewalls, precessing mainly along the streamwise direction within a narrow vertical plane. As $\varGamma _z$ increases, however, the confinement is gradually lifted: the same corner-vortex motion is no longer confined to a single plane; it can now extend in the spanwise direction, and the originally ordered spiral motion thus gives way to a more 3-D trajectory.

Figure 10.

Effect of the spanwise aspect ratio $\varGamma _z$ on the flow structure and temperature standard deviation $\sigma _T$ at $\textit{Ra} = 7.2 \times 10^6$ . (a–d) Mean velocity streamlines for varying aspect ratios $\varGamma _z = 0.2$ , $0.4$ , $0.6$ and $0.8$ , respectively, coloured by the local magnitude of $\sigma _T$ . The transition from constrained spiral precession to 3-D turbulent mixing is visually evident as $\varGamma _z$ increases. (e) Schematic illustration of the corner-vortex motion under increasing spanwise aspect ratio. Strong confinement keeps the motion nearly planar and streamwise-precessing, whereas relaxed confinement permits spanwise development and leads to a more 3-D vortex path. (f) Vertical profiles of $\sigma _T$ for the $\varGamma _z = 0.8$ case. Dashed lines ( $z_{ct}$ ) represent data extracted along the central vertical plane ( $x/H = 0.5$ ), while solid lines ( $z_{co}$ ) denote data extracted along vertical lines near the corner ( $x/H = 0.08$ ). The values in the legend indicate the specific spanwise depth ( $z$ ) for each extracted line. The profiles are predominantly unimodal across the domain, with the bimodal feature preserved only at specific highly localized corner positions (e.g. $z_{co} = 0.9W$ ).

This change in flow topology is directly reflected in the vertical profiles of $\sigma _T$ (figure 10 f for $\varGamma _z = 0.8$ ). To quantify this spatial dependence rigorously, profiles are extracted along specific vertical lines, distinguished by their horizontal and spanwise coordinates. Here, $z_{ct}$ denotes data extracted along the vertical centreline plane ( $x/H = 0.5$ ), while $z_{co}$ corresponds to data extracted near the corner ( $x/H = 0.08$ ). For example, $z_{co} = 0.9W$ represents the vertical profile at $x/H = 0.08$ and $z/W = 0.9$ , spanning $y$ from $0$ to $H$ . At $\varGamma _z = 0.8$ , the profiles show that the bimodal structure vanishes over most of the domain, collapsing into the classical unimodal shape characteristic of convection. A residual, asymmetric bimodal signature persists only in highly localized regions near the corner ( $e.g.$ , $z_{co} = 0.9W$ ), where the weakened corner roll still exerts influence. At this relatively large aspect ratio ( $\varGamma _z = 0.8$ ), the LSC exhibits some finite-amplitude wandering, which is typical for such wide domains (Sun et al. Reference Sun, Xi and Xia2005; Xi et al. Reference Xi, Zhou and Xia2006; Chen et al. Reference Chen, Xie, Yang and Ni2023). This wandering is already incorporated into the time-averaged $\sigma _T$ profiles and does not alter their extraction, as shown in Appendix B. Thus, the geometric constraint imposed by the spanwise walls is necessary for sustaining the bimodal temperature fluctuations; relaxing this confinement inevitably disrupts the coherent 3-D topology of the corner rolls, thereby eliminating the structured spanwise disturbances required to maintain the secondary thermal peak.

To isolate and quantify the influence of the near-wall regions, we introduce an exclusion parameter $\phi \in [0, 0.5]$ to define a truncated analysis domain

(3.6) \begin{align} \varOmega (\phi ) = [\phi L_x, (1-\phi )L_x] \times [0, H] \times [\phi L_z, (1-\phi )L_z] , \end{align}

with truncated cross-sectional area $A(\phi ) = (1-2\phi )^2 L_x L_z$ . The corresponding area averaged temperature standard deviation is then

(3.7) \begin{align} \langle \sigma _T\rangle _A(y;\phi )=\sqrt {\frac {1}{A(\phi )}\iint _{\varOmega (\phi )\cap \{y=\textit{const}\}}{\langle T^{\prime }\rangle ^2\mathrm{d}x\mathrm{d}z}}. \end{align}

The area-averaged temperature standard deviation profiles $\langle \sigma _T \rangle _A$ exhibit a strong dependence on $\phi$ , as illustrated in figure 11(a). When the extraction width covers the central region and excludes less than 50 % of the near-wall span (i.e. $\phi \leqslant 25\,\%$ ), the profile remains fundamentally unimodal, dominated by the core flow characteristics. As indicated by the downward arrow in figure 11(a), increasing $\phi$ causes the analysis region to retract inward, progressively excluding near-wall regions to isolate the influence of the central core. The arrow indicates the direction of increasing $\phi$ , demonstrating that restricting statistical analysis to the interior away from spanwise boundaries causes a bimodal structure to emerge gradually. Once the exclusion exceeds 50 % ( $\phi \gt 25\,\%$ ), probing deeper into corner-roll-dominated domains, the bimodal structure becomes evident. As the extraction region contracts farther inward, moving farther from the spanwise boundaries, the two distinct off-wall peaks grow increasingly pronounced. This progressive emergence suggests that the bimodal phenomenon is a robust statistical consequence of 3-D flow interactions under wall-confined conditions, as the area-averaged procedure eliminates random contributions from individual coherent structures.

Figure 11.

The essential role of spanwise walls and near-wall regions in the formation of bimodal $\sigma _T$ profiles. (a) Area-averaged temperature standard deviation profiles $\langle \sigma _T \rangle _A$ for various exclusion parameters $\phi$ . The arrow indicates the direction of increasing $\phi$ , demonstrating the gradual emergence of the bimodal structure as the statistical analysis is restricted to the interior away from the spanwise boundaries. (b) Comparison of vertical $\sigma _T$ profiles among the original confined 3-D simulation, a 3-D simulation with periodic boundary conditions in the spanwise direction (cyclic) and a pure two-dimensional (2-D) simulation. The distinct bimodal peaks are exclusive to the confined 3-D case.

Figure 11. Long description

Panel A: A line graph shows area-averaged temperature standard deviation profiles for various exclusion parameters. The x-axis is labeled y/H and ranges from 0 to 1. The y-axis is labeled <σT>_A and ranges from 0 to 0.8. Different lines represent different exclusion parameters, with colors and line styles indicating specific values. The arrow indicates the direction of increasing exclusion parameter, showing the gradual emergence of a bimodal structure as the analysis moves away from the spanwise boundaries. Panel B: Another line graph compares vertical profiles among different simulation types. The x-axis is labeled y/H and ranges from 0 to 1. The y-axis is labeled σT and ranges from 0 to 1.2. Three lines represent different simulations: a dashed line for the original confined 3-D simulation, a solid orange line for a 3-D simulation with periodic boundary conditions, and a solid blue line for a pure 2-D simulation. The confined 3-D case shows distinct bimodal peaks.

A more direct demonstration that the bimodal structure inherently relies on the physical presence of confining spanwise walls is provided in figure 11(b). Here, we compare the vertical $\sigma _T$ profile from the original confined 3-D simulation against two modified configurations: a purely 2-D simulation ( $\sigma _{T { (2\hbox{-}\textrm{D})}}$ ) and a 3-D simulation utilizing periodic boundary conditions in the $z$ -direction ( $\sigma _{T(\textit{cyclic})}$ ). Except for the spanwise boundary conditions, all numerical settings in these modified cases, including the vertical and streamwise grid resolution, solver configuration and bounding-wall thermal conditions, are identical to those in the original 3-D confined case. The resulting profiles clearly show that the pronounced off-wall bimodal peaks, which are a defining feature of the original 3-D confined case, are completely absent in both modified configurations, which instead exhibit standard unimodal or flattened profiles. In the 2-D and cyclic configurations, while corner vortices natively form, the absence of rigid spanwise confinement liberates them from the forced spiral precession. Consequently, they extend freely and uniformly in the spanwise direction. Without the geometric forcing imposed by the spanwise walls, the corner rolls fail to execute the coherent spiral motion, and the resulting strong spanwise kinetic disturbances vanish. As a result, the active transport mechanism that sweeps boundary-layer plumes into the core is severed, rendering the flow unable to sustain the bimodal temperature fluctuations. Consequently, neither modified condition reproduces the bimodal characteristics observed in the original flow. Together with the aspect-ratio study and the wall-exclusion analysis, this comparison further supports the proposed mechanism: spanwise confinement organizes the corner vortices into a coherent 3-D motion, which generates spanwise velocity fluctuations that transport boundary-layer plumes into the core and sustain the bimodal $\sigma _T$ profiles.

Figure 12.

Scaling of the global $\textit{Nu}$ with the $\textit{Ra}$ for low-Prandtl-number fluids under varying spanwise confinements ( $\varGamma _z = 0.2, 0.6$ , and $1.0$ ). The solid orange line indicates the power-law fit for the strongly confined Q2D case ( $\varGamma _z = 0.2$ ), yielding an effective scaling of $\textit{Nu} \sim \textit{Ra}^{0.332}$ . The inset displays the compensated plot, $\textit{Nu}\textit{Ra}^{-0.332}$ versus $\textit{Ra}$ , which highlights the robustness of the $0.332$ scaling for $\varGamma _z = 0.2$ . For geometries with relaxed confinement ( $\varGamma _z = 0.6$ and $1.0$ ), the scaling exponents decrease to 0.278 and 0.269, respectively, evident from the downward trend in the compensated plot. The grey dashed line represents the classical $\textit{Ra}^{1/4}$ scaling predicted by the GL theory. The black open square and grey open diamond denote the low- $\textit{Pr}$ cylindrical-cell data of Schumacher et al. (Reference Schumacher, Bandaru, Pandey and Scheel2016) and Scheel & Schumacher (Reference Scheel and Schumacher2017), respectively, at aspect ratio unity, and are included for comparison with the relaxed-confinement cases.

3.4. Heat transfer scaling of low-Prandtl-number fluids under Q2D confinement

Figure 12 presents the dependence of global heat transfer efficiency on $\textit{Ra}$ under Q2D confinement. For the severely confined case ( $\varGamma _z = 0.2$ ), the $\textit{Nu}$ exhibits an effective scaling of $\textit{Nu} \sim \textit{Ra}^{0.332}$ during the initial turbulent regime. This scaling is further corroborated by the distinct horizontal plateau in the compensated plot ( $\textit{Nu}\textit{Ra}^{-0.332}$ versus $\textit{Ra}$ ) shown in the inset. This exponent noticeably deviates from the classical Grossmann–Lohse (GL) prediction of $\textit{Nu} \sim \textit{Ra}^{1/4}$ for low- $\textit{Pr}$ fluids. To isolate the effect of confinement, supplementary simulations with relaxed spanwise boundaries ( $\varGamma _z = 0.6$ and $1.0$ ) were conducted. As the geometric constraint weakens and the flow regains full three-dimensionality, the scaling exponent systematically decreases to 0.278 and 0.269, respectively, thereby gradually recovering the classical GL scaling. This transition is distinctly reflected in the inset, where the compensated data for $\varGamma _z = 0.6$ and $1.0$ display a pronounced negative slope. The confinement-induced steepening of the scaling law aligns closely with the experimental observations by Glazier et al. (Reference Glazier, Segawa, Naert and Sano1999), who reported an elevated exponent of 0.285 in a $\varGamma _z = 0.5$ cell. In contrast, high-resolution numerical studies of low- $\textit{Pr}$ convection in standard aspect-ratio-unity cylinders by Schumacher et al. (Reference Schumacher, Bandaru, Pandey and Scheel2016) and Scheel & Schumacher (Reference Scheel and Schumacher2017) exhibit agreement with both our relaxed-confinement cases and the classical GL scaling.

While the strongly confined configuration yields the highest effective scaling exponent in the present simulations, it does not produce the maximum absolute heat transport. To clarify this distinction and to compare with the convention of Chong & Xia (Reference Chong and Xia2016), a supplementary analysis based on $\textit{Nu}-1$ is provided in Appendix C. Although the currently observed exponent of 0.332 nominally approaches the asymptotic limit ( $\textit{Nu} \sim \textit{Ra}^{1/3}$ ) proposed by Malkus (Reference Malkus1954) based on the ‘maximum heat flux’ principle, the underlying physical mechanisms are fundamentally distinct. Rather than indicating a transition towards an idealized unconfined ultimate state, this elevated exponent reflects an altered sensitivity of $\textit{Nu}$ to $\textit{Ra}$ driven entirely by spanwise restrictions.

Figure 13.

Instantaneous temperature fields illustrating the effect of spanwise confinement on thermal plume morphology. Panels (a)–(c) and (g)–(i) display the temperature contours at the spanwise midplane ( $z/W = 0.5$ ) for $\textit{Ra} = 2.57 \times 10^6$ and $1.5 \times 10^7$ , respectively, across aspect ratios $\varGamma _z = 0.2, 0.6$ , and $1.0$ . Panels (d)–(f) and (j)–(l) show the corresponding horizontal cross-sections taken near the edge of the thermal boundary layer ( $y = \delta _{th} = H/(2\textit{Nu})$ ). Under severe geometric constraint ( $\varGamma _z = 0.2$ , panels (d) and (j), plumes are more concentrated, whereas higher $\textit{Ra}$ (panels (g)–(l)) and larger $\varGamma _z$ promote spanwise dispersion.

Figure 13. Long description

Panel A: A heat map showing temperature distribution at the spanwise midplane for an aspect ratio of 1. Panel B: A heat map showing temperature distribution at the spanwise midplane for an aspect ratio of 2. Panel C: A heat map showing temperature distribution at the spanwise midplane for an aspect ratio of 4. Panel D: A heat map showing horizontal cross-sections of temperature distribution near the edge of the thermal boundary layer for an aspect ratio of 1. Panel E: A heat map showing horizontal cross-sections of temperature distribution near the edge of the thermal boundary layer for an aspect ratio of 2. Panel F: A heat map showing horizontal cross-sections of temperature distribution near the edge of the thermal boundary layer for an aspect ratio of 4. Panel G: A heat map showing temperature distribution at the spanwise midplane for an aspect ratio of 1 under higher conditions. Panel H: A heat map showing temperature distribution at the spanwise midplane for an aspect ratio of 2 under higher conditions. Panel I: A heat map showing temperature distribution at the spanwise midplane for an aspect ratio of 4 under higher conditions. Panel J: A heat map showing horizontal cross-sections of temperature distribution near the edge of the thermal boundary layer for an aspect ratio of 1 under higher conditions. Panel K: A heat map showing horizontal cross-sections of temperature distribution near the edge of the thermal boundary layer for an aspect ratio of 2 under higher conditions. Panel L: A heat map showing horizontal cross-sections of temperature distribution near the edge of the thermal boundary layer for an aspect ratio of 4 under higher conditions. Each heat map uses a color scale ranging from blue to red, indicating lower to higher temperature values. The x and y axes represent normalized coordinates, with the color intensity showing the temperature distribution across the planes.

The macroscopic alterations in heat transport are intimately linked to the morphological evolution of the thermal plumes, as visualized in figure 13. At a moderate Rayleigh number ( $\textit{Ra} = 2.57 \times 10^6$ ), severe confinement ( $\varGamma _z = 0.2$ ) forces the thermal structures to coalesce into concentrated plumes, a phenomenon readily observed at the spanwise midplane (figure 13 a) and near the edge of the thermal boundary layer (figure 13 d). In stark contrast, as the spanwise constraint is relaxed to $\varGamma _z = 0.6$ and $1.0$ , the plumes regain their spanwise degrees of freedom, exhibiting a more fragmented distribution in both the midplane (figures 13 b and 13 c) and the boundary-layer regions (figures 13 e and 13 f). Although this structural condensation under Q2D confinement is consistent with the enhanced boundary friction effects reported previously (Huang et al. Reference Huang, Kaczorowski, Ni and Xia2013), the confinement in the present cell is insufficient to trigger the formation of highly efficient ‘giant plumes’ (Chong et al. Reference Chong, Huang, Kaczorowski and Xia2015; Chong & Xia Reference Chong and Xia2016). Consequently, this restricted spanwise mobility, coupled with the absence of compensatory giant plumes, suppresses the vigour of convective transport and leads to the observed reduction in absolute $\textit{Nu}$ .

As the Rayleigh number increases further to $1.5 \times 10^7$ , the intensified buoyancy driving forces enable the plumes to partially overcome the geometric suppression. The thermal structures in the Q2D cell become noticeably more dispersed (figures 13 g and 13 j), weakening the condensation effect. This transition towards a more disordered state is further corroborated by the results in the less-confined geometries ( $\varGamma _z = 0.6, 1.0$ ), where the plumes appear increasingly chaotic and spatial correlations are further diminished (figures 13 h, 13 i and 13 k, 13 l). From a structural perspective, this dispersed morphology at high $\textit{Ra}$ functionally mimics the plume states observed at lower $\textit{Ra}$ but larger $\varGamma _z$ . This morphological recovery towards a more 3-D state provides a physical rationale for why the scaling exponent eventually reverts towards the classical 3-D predictions at higher $\textit{Ra}$ regimes.