5.1. $ \textit{Ra}$ dependence of $ \textit{Nu}$ with varying BCs

Figure 3.

(a) $ \textit{Nu}$ versus $ \textit{Ra}$ in logarithmic coordinates for different boundary conditions from 2-D and 3-D DNS results. Open symbols denote the 3-D DNS data for $ \textit{Pr} = 10$ and $ \textit{Pr} = 0.1$ , from Shishkina (Reference Shishkina2017). (b) Compensated heat transfer coefficient $ \textit{Nu}\textit{Ra}^{-1/4}$ as a function of $ \textit{Ra}$ . The dashed line, the solid line and the dash-dotted line indicate the $ \textit{Nu} \sim Ra^{1/5}$ , $ \textit{Nu} \sim Ra^{1/4}$ and $ \textit{Nu} \sim Ra^{1/5}$ scalings, respectively.

Figure 3. Long description

Panel A: A scatter plot shows the relationship between Nusselt number (Nu) and Rayleigh number (Ra) in logarithmic coordinates. The x-axis represents Ra, ranging from 10^8 to 10^12, and the y-axis represents Nu, ranging approximately from 10 to 200. Different symbols and colors denote various boundary conditions from 2-D and 3-D DNS results. Open symbols represent 3-D DNS data for specific Prandtl numbers. Panel B: A scatter plot shows the compensated heat transfer coefficient as a function of Ra. The x-axis represents Ra, ranging from 10^8 to 10^12, and the y-axis represents NuRa^-1/4, ranging from 0.2 to 0.35. Dashed, solid, and dash-dotted lines indicate different scalings: Nu ~ Ra^1/5, Nu ~ Ra^1/4, and Nu ~ Ra^1/5, respectively.

Figure 3 shows the relationship between $ \textit{Nu}$ and $ \textit{Ra}$ . Over the present parameter range, the data exhibit an effective scaling exponent lying between $1/5$ and $1/4$ , consistent with the Shishkina–Grossmann–Lohse (SGL) framework for HC (Shishkina et al. Reference Shishkina, Grossmann and Lohse2016), which predicts a transition from the classical laminar regime ( $ \textit{Nu} \sim Ra^{1/5}$ ) to a higher- $ \textit{Ra}$ regime characterised by $ \textit{Nu} \sim Ra^{1/4}$ at moderate $ \textit{Pr}$ . Under the free-slip condition, heat transfer follows $ \textit{Nu} \sim Ra^{1/5}$ , whereas the no-slip data over the same range approach the $ \textit{Ra}^{1/4}$ regime. This difference may be attributed to $ \textit{Pr}$ dependence: the present value $ \textit{Pr}=4.3$ lies between the two Prandtl numbers considered by Shishkina (Reference Shishkina2017), namely $ \textit{Pr}=0.1$ and $ \textit{Pr}=10$ . According to the existing phase diagram (Shishkina & Wagner Reference Shishkina and Wagner2016; Passaggia & Scotti Reference Passaggia and Scotti2024), $ \textit{Pr}=4.3$ may fall in the same laminar heat-transfer regime $I_l$ as the $ \textit{Pr}=0.1$ no-slip data, which also display a $1/5$ scaling over a comparable $ \textit{Ra}$ interval. However, this interpretation is not supported by the no-slip comparison at the same aspect ratio $\varGamma$ . As shown in figure 3, the $ \textit{Pr}=4.3$ data collapse with the $ \textit{Pr}=10$ branch of Shishkina (Reference Shishkina2017), indicating that $ \textit{Pr}=4.3$ is already large enough to place the flow in the high- $ \textit{Pr}$ regime $I_{l}^{*}$ .

We therefore attribute the $ \textit{Nu} \sim Ra^{1/5}$ scaling under free-slip primarily to the change in velocity boundary condition. The free-slip condition does not simply modify the prefactor of the heat-transfer relation, but rather alters the regime organisation of HC by shifting the transition between regimes. The present results indicate that while the scaling exponents themselves remain robust, the locations of regime transitions are sensitive to the imposed velocity boundary conditions. In particular, free-slip boundaries shift the onset of the transitional regime without altering its fundamental scaling behaviour.

5.2. Velocity profiles for varying BCs

The upper boundary of the HC system comprises the heating plate, the cooling plate and the intermediate adiabatic region between them. Near the thermally driven plates, the boundary condition is isothermal, while in the interior region between the plates, both experiments and simulations employ adiabatic boundary conditions. The mean velocity profiles exhibit distinct patterns depending on both the thermal regime and the imposed slip condition, and the velocity field evolves progressively from the heating plate through the intermediate region towards the turbulent zone beneath the cooling plate.

Figure 4.

Velocity profiles for varying slip and thermal boundary conditions at different horizontal locations $x/l_0$ and $ \textit{Ra}=10^9$ : (a) no-slip, isothermal (Configuration I); (b) no-slip, adiabatic (Configuration II); (c) free-slip, isothermal (Configuration III); (d) free-slip, adiabatic (Configuration IV).

Figure 4 presents four representative examples of the velocity distribution in the isothermal and adiabatic regions. Figure 4(a,b) shows the mean velocity profiles at various streamwise locations under no-slip boundary conditions, while figure 4(c,d) shows the corresponding free-slip cases. Notably, within a region adjacent to the upper boundary, the velocity profiles at different streamwise positions exhibit self-similar characteristics. The horizontal pressure gradient, however, is confined to the thermal boundary-layer region. Once the depth exceeds a threshold value, the horizontally averaged temperature in the bulk becomes nearly uniform. Consequently, the pressure-driven horizontal flow decelerates and eventually reverses direction. Owing to the no-slip condition on the remaining boundaries, the flow gradually decreases from the bulk towards the bottom and undergoes oscillatory adjustment to zero. This oscillatory behaviour is related to the aspect ratio $\varGamma$ of the domain.

To ensure continuity of the velocity boundary condition while minimising transitional effects between adjacent regions, we first consider configurations in which both the isothermal and adiabatic sections are assigned the same velocity boundary condition – either uniformly no-slip or uniformly free-slip.

5.3. Velocity and temperature profiles in Configuration I: no-slip and isothermal boundary

Figure 5.

Vertical profiles of the mean streamwise velocity for (a,b) Configuration I with no-slip, isothermal boundary conditions and (c,d) Configuration II with no-slip, adiabatic boundary conditions. (a,c) Normalised velocity, $U^+ \equiv U/u_{\tau }$ , as a function of $z^+$ at several representative streamwise locations $x/l_0$ in each configuration. The dashed lines denote the linear viscous relation, $U^+ = z^+$ . (b,d) Normalised velocity, $U/U_0$ , plotted against $\xi$ using the same data as in panels (a,c), respectively. The solid lines denote the analytical profile $\operatorname {sech}^2[a_1(\xi -\xi _p)]$ for $\xi \geqslant \xi _p$ , where the fitting parameter $a_1$ and the peak location $\xi _p$ are indicated in the figure.

Figure 5. Long description

Panel A: A line graph shows normalized velocity (U+) as a function of z+ at several streamwise locations (x/l0 = 0.2, 0.5, 0.8) for Configuration I with no-slip, isothermal boundary conditions. The dashed line represents the linear viscous relation U+ = z+. Panel B: A line graph plots normalized velocity (U/U0) against ξ using the same data as Panel A. The solid line denotes the analytical profile for U/U0, with fitting parameter 1.72 and peak location 0.49 indicated. Panel C: A line graph shows normalized velocity (U+) as a function of z+ at several streamwise locations (x/l0 = 2, 5, 8) for Configuration II with no-slip, adiabatic boundary conditions. The dashed line represents the linear viscous relation U+ = z+. Panel D: A line graph plots normalized velocity (U/U0) against ξ using the same data as Panel C. The solid line denotes the analytical profile for U/U0, with fitting parameter 1.54 and peak location 0.44 indicated.

Before introducing the analytical model, we first examine the near-wall velocity structure to assess the validity of the laminar boundary-layer assumption. Figure 5(a,c) presents the normalised streamwise velocity, $ U^+ \equiv U/u_\tau$ (with $u_{\tau }=\sqrt {\nu \,\partial U/\partial z|_0}$ ) as a function of the wall-normal distance $z^+ \equiv z u_\tau / \nu$ at varying $x/l_0$ for no-slip boundary condition. Despite of the convection flow in the bulk, the velocity profiles largely follow the linear viscous scaling $ U^+ = z^+$ , indicating that the boundary layer retains the characteristic structure of a laminar wall layer and, thereby, supporting the laminar formulation in the subsequent theoretical analysis.

It is worth noting that the linear viscous relation $ U^+ = z^+$ persists only up to $z^+ \approx 1$ , substantially smaller than the viscous sublayer thickness in canonical wall-bounded shear flows, where the linear scaling typically extends to $z^+ \approx 5$ . This difference can be understood from the near-wall momentum balance. Expanding the streamwise velocity for the no-slip wall gives

(5.1) \begin{align} U(z) = \left . \frac {\partial U}{\partial z} \right |_{0} z + \frac {1}{2} \left . \frac {\partial ^2 U}{\partial z^2} \right |_{0} z^2 + O \big(z^3 \big). \end{align}

Since both velocity components vanish at the wall (shown in figure 17), the convective terms in the boundary-layer momentum equation are identically zero at $ z = 0$ , yielding

(5.2) \begin{align} \left . \nu \frac {\partial ^2 U}{\partial z^2} \right |_{z=0} = \left . \frac {1}{\rho } \frac {\partial P}{\partial x} \right |_{z=0} , \end{align}

with $\partial P/\partial x \lt 0$ near the wall. Expressed in wall units,

(5.3) \begin{align} U^+ = z^+ + \frac {1}{2} \left ( \frac {\nu }{\rho u_\tau ^3} \frac {\partial P}{\partial x} \right ) (z^+)^2 + O \big((z^+)^3 \big), \end{align}

showing that the departure from the linear viscous law is governed by the dimensionless pressure-gradient term $ \nu (\partial P / \partial x) / (\rho u_\tau ^3)$ . In canonical wall-bounded shear flows, this quantity is generally small, so the quadratic correction remains negligible over the viscous sublayer. In HC, however, the flow is driven primarily by the buoyancy-induced horizontal pressure gradient $\partial P/\partial x$ , while the wall shear velocity $u_{\tau }$ is comparatively weak. As a result, the coefficient of the quadratic term is significantly larger, thereby causing the velocity profile to bend more rapidly towards the near-wall velocity maximum.

Figure 5(b,d) shows the corresponding velocity profiles normalised by the local maximum velocity $ U_0$ . After normalisation, the profiles obtained at different streamwise locations collapse onto a nearly self-similar distribution, indicating that the streamwise evolution is well described by a single similarity variable. The solid line represents the analytical approximation,

(5.4) \begin{align} \frac {U}{U_0} = \operatorname {sech}^2[a_1(\xi - \xi _p)] \end{align}

for $\xi \geqslant \xi _p$ . The parameter values are shown in the figure. The excellent agreement with the DNS data confirms that the proposed streamfunction accurately captures the self-similar wall-jet structure and therefore provides the kinematic foundation for deriving the analytical temperature solutions presented in the following sections.

Under the no-slip and isothermal boundary conditions, the normalised velocity profiles exhibit the characteristic structure of a plane wall jet. Accordingly, the streamfunction introduced in § 3 is adopted to describe the self-similar velocity field. A plane wall jet comprises two fundamental components: the jet mode and the wall mode. In the jet mode, the characteristic velocity scale is defined by the maximum centreline velocity $U_0$ and the jet half-width $\sigma$ . For a wall jet, two distinct locations exist at which the velocity reaches half its maximum value; in this analysis, the position farther from the boundary is used to define the half-width. As demonstrated in figure 5(b), the jet mode resembles a free jet, but is shifted away from the wall to the position $z_p$ at which $U_0$ is attained. As shown in figure 2, the $U/U_0$ profiles at different $x/l_0$ beneath the no-slip heated plate are also well captured by (4.2).

Figure 6.

Normalised mean temperature profiles $\theta (\xi )$ in the vicinity of the heated plate as a function of $\xi = z/\lambda$ in Configuration I with no-slip and isothermal boundary conditions. Data are horizontally averaged over $0.1 \leqslant x/l_0 \leqslant 0.9$ for 2-D DNS at $ \textit{Ra} = 5\times 10^9$ (circles), $5\times 10^{10}$ (squares) and $5\times 10^{11}$ (triangles). Crosses denote the horizontal-area-averaged profile over the same $x$ -range and $0 \leqslant y/w_0 \leqslant 1$ from the 3-D DNS at $ \textit{Ra} = 5\times 10^{11}$ . The solid line in panel (a) represents the theoretical fit given by (4.5) with $a = 0.63$ , $k_1 = 5.48$ and $k_2 = 2.46$ . (b) Residuals $\theta _{\textit{data}} - \theta _{\textit{fit}}$ corresponding to each dataset in panel (a).

Figure 6(a) shows the normalised mean temperature profiles $\theta$ beneath the heated plate, at which the boundary conditions are no-slip and isothermal. The profiles obtained from both the 2-D and 3-D DNS collapse closely onto a single self-similar distribution over more than two decades in Rayleigh number, indicating that the thermal boundary-layer structure is largely independent of Ra after appropriate normalisation. The solid line represents the analytical solution given by (4.5), which agrees well with the DNS data throughout the entire thermal boundary layer. To assess the fitting quality more quantitatively, figure 6(b) shows the corresponding residuals, $\theta _{\textit{data}} - \theta _{\textit{fit}}$ for each dataset. The residuals remain within approximately 1.5 % over the whole profile, demonstrating the high accuracy of the analytical solution. The residual distributions exhibit smooth, systematic variations rather than random scatter, indicating that the reduced analytical formulation successfully captures the dominant boundary-layer physics. The remaining discrepancies most likely arise from higher-order effects neglected in the present reduced model, such as weak departures from exact self-similarity and finite streamwise-development effects.

Figure 7.

Ratio of the velocity characteristic length $\sigma$ to the thermal characteristic length $\lambda$ along the streamwise direction in Configuration I with no-slip and isothermal boundary conditions.

As shown in figure 7, the ratio of the velocity characteristic length $\sigma$ to the thermal characteristic length $\lambda$ approaches a constant value of approximately 0.9 in the isothermal region. The near-constant ratio indicates that the velocity and thermal boundary layers evolve proportionally in the streamwise direction, consistent with the local self-similarity observed in figure 5(b). This coupled spatial development supports the reduced boundary-layer description adopted throughout the analytical framework.

5.4. Velocity and temperature profiles in Configuration II: no-slip and adiabatic boundary

Figure 8.

Normalised mean velocity profiles in Configuration II with no-slip and adiabatic boundary conditions at (a) $ \textit{Ra}=10^9$ and (b) $ \textit{Ra}=10^{11}$ . The solid line shows a fit of (4.2) with $k_1=3.57$ and $k_2=1.6$ , and the dashed line shows (4.2) with $k_1=3.79$ and $k_2=1.7$ .

In the intermediate region between the heated and cooled plates, both experiments and simulations employ adiabatic boundary conditions. The presence of the adiabatic boundary renders the conventional definition of the local thermal boundary-layer thickness – defined by the intersection of the temperature-gradient tangent with the bulk temperature level – inappropriate for this configuration. It is found that when the characteristic length scale is chosen as $0.7$ times the depth at which the bulk temperature is reached, the normalised temperature profiles at different streamwise positions collapse onto a single curve, exhibiting self-similar behaviour. Hence, $\lambda = z_{0.7}$ is adopted as the reference thermal length scale in this section.

Figure 9.

Mean temperature profiles $\theta$ (panels ai–ci) and the corresponding differences $\theta _{\textit{data}} - \theta _{\textit{fit}}$ from the theoretical profiles (panels aii–cii) obtained in Configuration II with no-slip and adiabatic boundary conditions. (ai,aii) Profiles at various streamwise locations for $ \textit{Ra}=10^{11}$ . The solid line denotes equation (4.7) with parameters $H(\epsilon , \gamma , \delta , \alpha , q)=H(-4.54,2.13,-0.13,-1.23,-0.61)$ and coefficients $c_1=-0.91, c_2=-0.002$ . (bi,bii) Experimental profiles at various streamwise locations for $ \textit{Ra}=2.55\times 10^{12}$ . The solid line represents equation (4.7) with parameters $H(-4.13,2.03,-0.03,0.54,0.27)$ and coefficients $c_1=-1.17, c_2=0.002$ . (ci,cii) Horizontally averaged temperature profiles over the interval $3 \leqslant x/l_0 \leqslant 7$ from the 2-D DNS at $ \textit{Ra}=5\times 10^{9}$ (circles), $5\times 10^{10}$ (squares), and $5\times 10^{11}$ (triangles). Crosses denote the horizontal-area-averaged profile over the same $x$ -range and $0 \leqslant y/w_0 \leqslant 1$ from the 3-D DNS at $ \textit{Ra}=5\times 10^{11}$ . The solid line represents (4.7) with parameters $H(-2.81,1.70,0.30,-1.29,-0.64)$ and coefficients $c_1=-1.06, c_2=0.06$ . The dashed line represents the corresponding fit to the 3-D data, with parameters $H(-3.89,1.97,0.03,-0.06,-0.03)$ and coefficients $c_1=-1.11, c_2=0.02$ .

Figure 9. Long description

Panel A: Two line graphs. Panel Ai: The x-axis is labeled with ξ and the y-axis is labeled with θ(ξ). The graph shows multiple data series for different streamwise locations (x/l0 = 3, 4, 5, 6, 7) represented by different symbols (circles, squares, diamonds, triangles, stars). A solid line represents the theoretical prediction using equation (4.7) with specific parameters. Panel Aii: The x-axis is labeled with ξ and the y-axis is labeled with θ_data - θ_fit. The graph shows the differences between the data and the theoretical fit for the same streamwise locations. Panel B: Two line graphs. Panel Bi: The x-axis is labeled with ξ and the y-axis is labeled with θ(ξ). The graph shows experimental data for different streamwise locations (x/l0 = 3, 5, 7) represented by different symbols (circles, squares, triangles). A solid line represents the theoretical prediction using equation (4.7) with specific parameters. Panel Bii: The x-axis is labeled with ξ and the y-axis is labeled with θ_data - θ_fit. The graph shows the differences between the experimental data and the theoretical fit for the same streamwise locations. Panel C: Two line graphs. Panel Ci: The x-axis is labeled with ξ and the y-axis is labeled with θ(ξ). The graph shows horizontally averaged temperature profiles from 2-D DNS at different Rayleigh numbers (Ra = 5 x 10^9, 5 x 10^10, 5 x 10^11) represented by different symbols (circles, squares, triangles). Crosses denote the horizontal-area-averaged profile from 3-D DNS. A solid line represents the theoretical prediction using equation (4.7) with specific parameters. Panel Cii: The x-axis is labeled with ξ and the y-axis is labeled with θ_data - θ_fit. The graph shows the differences between the data and the theoretical fit for the same Rayleigh numbers.

In the adiabatic region with no-slip boundary conditions, the mean velocity profiles also exhibit wall-jet characteristics and are well described by the streamfunction (4.2), as shown in figure 8. Based on the established velocity structure model, the temperature structure can be further validated. Figure 9(a) presents the normalised temperature profiles at various streamwise locations under no-slip adiabatic boundary conditions for $ \textit{Ra}=10^{11}$ , showing good agreement with the analytical solution (4.7). The experimental temperature profiles shown in figure 9(b) likewise match the theoretical predictions well. Figure 9(c) demonstrates that the theoretical solution also agrees well with surface-averaged temperature profiles over a range of $ \textit{Ra}$ . As in the isothermal region, the normalised temperature profiles in the adiabatic region collapse when scaled with the reference thermal length scale. The model captures the distinct thermal characteristics near the boundary when comparing adiabatic and isothermal conditions: specifically, it correctly predicts the non-monotonic temperature gradient characterised by a slope that initially increases before decreasing, which is a key feature absent in the isothermal case. In addition, after spatial and temporal averaging, the 3-D mean temperature profiles remain close to the 2-D profiles under the adopted normalisation. This indicates that the dominant wall-normal thermal structure is governed primarily by the same coupling between the bulk and the boundary layer captured by the 2-D description.

Figure 10.

Ratio of the velocity characteristic length $\sigma$ to the thermal characteristic length $\lambda =z_{0.7}$ along the streamwise direction in Configuration II with no-slip and adiabatic boundary conditions.

The velocity profiles in the adiabatic region do not collapse as tightly as those in the isothermal region when normalised by the half-width $\sigma$ . Beneath the heated plate, the sustained heat flux maintains a stable near-wall shear layer and the mean flow retains a jet-like structure to a good approximation. In the adiabatic region, however, the flow emerging from the heated section is shaped by both the pressure gradient and the large-scale circulation. Consequently, the velocity profiles are jet-like in shape rather than true wall jets. The ratio $\sigma /\lambda$ is shown in figure 10; at moderately low $ \textit{Ra}$ , it remains approximately constant.

5.5. Velocity and temperature profiles in Configuration III: free-slip and isothermal boundary

In oceanic systems, the air–liquid interface is commonly approximated as a free-slip boundary and the maximum horizontal velocity $U_0$ is attained at the surface. Within the self-similar region, the normalised velocity $U/U_0$ collapses approximately onto a single curve within one characteristic length-scale region, suggesting that the near-surface velocity profile can be described by a universal function of $\xi$ .

Figure 11.

Mean velocity profiles normalised by $U_0$ and the corresponding length scale $\sigma$ in Configuration III with free-slip and isothermal boundary conditions. Various markers correspond to different streamwise positions at (a) $ \textit{Ra}=10^9$ and (b) $ \textit{Ra}=10^{10}$ . The solid and dashed lines show (4.11) with $k=0.82$ and $0.8$ , respectively.

Figure 11 presents the normalised vertical profiles of mean horizontal velocity $U/U_0$ at $ \textit{Ra} = 10^9$ and $ \textit{Ra} = 10^{10}$ . Beneath the thermally driven boundary, all velocity profiles obtained across the range $0.2 \le x/l_0 \le 0.8$ are well collapsed up to $z/\sigma \approx 1$ and they closely resemble a plane free jet propagating along a boundary. The solid and dashed lines denote the empirical jet-flow equation (4.11) with $k = 0.82$ at $ \textit{Ra} = 10^9$ and $k = 0.80$ at $ \textit{Ra} = 10^{10}$ . Unlike no-slip boundaries where viscous retardation forces $U(0)=0$ and shifts the velocity peak into the interior, the free-slip condition allows maximum momentum transport along the top surface, causing the shear layer to broaden unimpeded downward into the bulk fluid.

The corresponding vertical temperature distribution beneath the isothermal plate is evaluated using the integral solution derived in § 4.3. Figure 12(a) compares the normalised temperature profiles $\theta (\xi )$ against 2-D DNS data across three orders of magnitude in $ \textit{Ra}$ ( $10^8 \le Ra \le 10^{10}$ ), as well as 3-D DNS data at $ \textit{Ra} = 10^{10}$ . All profiles are horizontally averaged over $0.1 \le x/l_0 \le 0.9$ . The theoretical curves, calculated via (4.12) with parameter sets listed in the legend, show robust agreement with the numerical data across the region $\xi \le 4$ . For comparison, the dash-dotted line illustrates the corresponding theoretical profile in Configuration I with no-slip and isothermal conditions. The free-slip profile exhibits a visibly sharper approach to the outer ambient value $\theta \to 1$ , reflecting the absence of shear-driven thermal mixing near $z=0$ .

Figure 12.

(a) Normalised mean temperature profiles $\theta (\xi )$ in the vicinity of the heated plate as a function of $\xi = z/\lambda$ in Configuration III with free-slip and isothermal boundary conditions. Data are horizontally averaged over $0.1 \leqslant x/l_0 \leqslant 0.9$ for 2-D DNS at $ \textit{Ra} = 10^{8}$ (circles), $10^{9}$ (squares) and $10^{10}$ (triangles). Crosses denote the horizontal-area-averaged profile over the same $x$ -range and $0 \leqslant y/w_0 \leqslant 1$ from the 3-D DNS at $ \textit{Ra} = 10^{10}$ . Curves of (4.12) are shown for $k=2.6,\ c=-0.5$ (solid line), $k=3.61,\ c=-0.34$ (dashed line), and $k=6.36,\ c=-0.19$ (dotted line). For comparison, the dash-dotted line represents the theoretical profile for no-slip and isothermal boundaries shown in figure 6. (b) Residuals $\theta _{\textit{data}} - \theta _{\textit{fit}}$ corresponding to each dataset in panel (a).

Figure 12. Long description

Panel A: Two line graphs showing normalized mean temperature profiles as a function of xi in the vicinity of the heated plate. The horizontal axis is labeled xi, and the vertical axis is labeled theta(xi). The graph includes data points for different Rayleigh numbers: Ra = 10^8 (red circles), Ra = 10^9 (green squares), Ra = 10^10 (blue triangles), and Ra = 10^10 (3-D) (orange crosses). Several curves represent different equations: solid line for Eq. (4.12) with k = 2.6, c = -0.5; dashed line for Eq. (4.12) with k = 3.61, c = -0.34; dotted line for Eq. (4.12) with k = 6.36, c = -0.19; and dash-dotted line for Eq. (4.5) with a = 0.63, k1 = 5.48, k2 = 2.46. Panel B: A line graph showing residuals corresponding to each dataset in Panel A. The horizontal axis is labeled xi, and the vertical axis is labeled theta_data - theta_fit. The graph includes the same data points for different Rayleigh numbers: Ra = 10^8 (red circles), Ra = 10^9 (green squares), Ra = 10^10 (blue triangles), and Ra = 10^10 (3-D) (orange crosses).

To rigorously quantify the fidelity of this theoretical framework, figure 12(b) plots the residuals $\theta _{\textit{data}} - \theta _{\textit{fit}}$ . Across all 2-D DNS cases, the local deviations remain well within $\pm 0.01$ throughout the near-wall region $\xi \le 4$ . For the 3-D DNS dataset (crosses), a mild deviation of $\approx 0.027$ emerges around $\xi = 1.7$ . This slight mismatch stems from secondary three-dimensional spanwise plumes that perturb the thermal transition layer, a feature inherently absent in 2-D formulations. Nevertheless, the maximum error across all datasets remains below $4\,\%$ , confirming that (4.12) captures the essential coupling between the surface jet velocity and conductive-convective heat transport under free-slip and isothermal conditions.

As shown in figure 13, the characteristic length-scale ratio $\sigma /\lambda$ approaches a near-constant plateau in the central region of the isothermal plate. Its $ \textit{Ra}$ dependence, however, differs markedly from that in Configuration II, as shown in figure 10. Under the no-slip boundary condition, $\sigma /\lambda$ increases gradually with increasing $ \textit{Ra}$ , whereas under the free-slip boundary condition (Configuration III), it decreases. This opposite trend reflects the fundamentally different mechanisms governing momentum transport. Under no-slip conditions, wall shear promotes viscous diffusion of momentum away from the boundary, causing the momentum layer to broaden relative to the thermal boundary layer as $ \textit{Ra}$ increases. In contrast, under free-slip conditions, the absence of tangential shear confines the pressure-driven surface jet closer to the interface, resulting in a reduction of the momentum thickness relative to the thermal boundary-layer thickness. The contrasting $ Ra$ -dependence of $ \sigma / \lambda$ therefore highlights the distinct responses of the momentum boundary layer to increasing thermal forcing under no-slip and free-slip boundary conditions.

Figure 13.

Ratio of the velocity characteristic length $\sigma$ to the thermal characteristic length $\lambda$ along the streamwise direction in Configuration III with free-slip and isothermal boundary conditions.

5.6. Velocity and temperature profiles in Configuration IV: free-slip and adiabatic boundary

Figure 14.

Normalised mean velocity profiles in Configuration IV with free-slip and adiabatic boundary conditions at different $x/l_0$ at (a) $ \textit{Ra}=10^9$ and (b) $ \textit{Ra}=10^{10}$ . The solid and dashed lines show (4.11) with $k=0.84$ and $k=0.85$ , respectively.

Figure 14. Long description

Panel A and Panel B: Two line graphs compare normalized mean velocity profiles in Configuration 4 with free-slip and adiabatic boundary conditions at different positions. The x-axis represents the normalized height (z/σ) ranging from 0 to 1. The y-axis represents the normalized mean velocity (U/U0) ranging from 0.5 to 1.0. Different symbols represent various positions (x/l0 = 3, 4, 5, 6, 7). Solid and dashed lines show the fit of Equation (4.11) with k = 0.84 and k = 0.85, respectively. The data points and lines indicate how the velocity profiles change with height under the specified conditions.

As the fluid beneath the heated plate develops free jet-like characteristics, the pressure-driven flow continues to evolve through the adiabatic interior region. Analogously to jet development downstream of a nozzle, the velocity profiles in the intermediate adiabatic region can be described by the free-jet model, where the half-width provides an appropriate scaling parameter for self-similarity, as demonstrated in figure 14.

Figure 15.

(ai–ci) Mean temperature profiles $\theta$ and the corresponding deviations, $\theta _{\textit{data}}-\theta _{\textit{ fit}}$ , from the (aii–cii) theoretical profiles for Configuration IV with free-slip and adiabatic boundary conditions. In panels (ai–ci), the theoretical curves are given by (4.14), with the corresponding fitting parameters indicated in each panel. (ai,aii) Temperature profiles at various streamwise locations $x/l_0$ from the 2-D DNS at $ \textit{Ra}=10^{10}$ . (bi,bii) Experimental temperature profiles at various $x/l_0$ for $ \textit{Ra}=2.55\times 10^{12}$ . (ci,cii) Horizontally averaged temperature profiles over $3\le x/l_0\le 7$ from the 2-D DNS at $ \textit{Ra}=10^{8}$ (circles), $10^{9}$ (squares) and $10^{10}$ (triangles). Crosses denote the horizontal-area-averaged profile over the same streamwise interval and $0\le y/w_0\le 1$ from the 3-D DNS at $ \textit{Ra}=10^{10}$ . The solid, dashed, dotted and dash-dotted curves represent the theoretical fits to the circle, square, triangle and cross datasets, respectively.

Figure 15(a) presents the normalised temperature profiles at various streamwise locations under free-slip adiabatic boundary conditions for $ \textit{Ra}=10^{10}$ , showing good agreement with the analytical solution (4.14). The experimental results shown in figure 15(b) are likewise well reproduced. However, as shown in figure 15(c), the agreement for surface-averaged temperature profiles at different $ \textit{Ra}$ is less satisfactory under the free-slip condition than under no-slip, particularly in the region where the profiles transition into the bulk. In the adiabatic region, the mean temperature equation contains an additional streamwise evolution term $\text{d}\psi '(1-\theta )$ , which is absent in the isothermal case. Under free-slip conditions, the near-wall flow is more susceptible to modulation by outer recirculation and entrainment, as no tangential stress is imposed at the boundary. Consequently, increasing $ \textit{Ra}$ shifts the position and curvature of the transitional layer of $\theta (\xi )$ , consistent with the observation that free-slip boundaries shift the onset of the transitional regime.

It is noted that in figure 15(a,c), the fitting coefficient $c_2=0$ . To satisfy the adiabatic boundary condition, the solution (4.14) requires

(5.5) \begin{align} \theta '(0)=c_2kF_2(\alpha -\gamma +1,\beta -\gamma +1,2-\gamma ,0)=0, \end{align}

which implies $c_2=0$ . In the experimental measurements, however, the extremely thin near-surface region of negligible temperature variation is difficult to resolve and setting $c_2=0$ does not yield a satisfactory fit to the experimental data, as shown in figure 15(b). Nevertheless, for ideal free-slip adiabatic conditions, the mean temperature profile reduces to the compact expression:

(5.6) \begin{align} \theta (\xi )=1 - F_1(\alpha ,\beta ,\gamma ;h(\xi )) . \end{align}

Figure 16.

Ratio of the velocity characteristic length $\sigma$ to the thermal characteristic length $\lambda =z_{0.7}$ along the streamwise direction in Configuration IV with free-slip and adiabatic boundary conditions.

The ratio $\sigma /\lambda$ is shown in figure 16. As in Configuration II, it remains approximately constant at low $ \textit{Ra}$ , but increases progressively with $x$ as $ \textit{Ra}$ is raised. This reflects the sustained entrainment and enhanced shear instability of the emerging flow, which promotes vertical momentum exchange. By contrast, the zero-flux constraint of the adiabatic boundary limits the growth of the thermal layer, resulting in a divergence between the momentum and thermal broadening rates.

5.7. Validation of the reduced analytical framework

The analytical solutions derived in § 4 contain a small number of coefficients that characterise the local velocity and temperature fields. These coefficients are determined independently for each configuration by fitting the analytical profiles to the corresponding collapsed DNS or experimental profiles presented in the preceding figures. Specifically, the normalised velocity profiles are first determined to identify the appropriate functional form of the streamfunction and the relation between its shape parameters. The velocity-related coefficients ( $k$ or $k_2$ ) and the temperature-related coefficients ( $a$ and, where applicable, $d$ ) are subsequently obtained from the corresponding normalised temperature profiles. Once these coefficients are specified, the integration constants in the analytical solutions are uniquely determined by the prescribed thermal boundary conditions. No additional adjustable parameters are introduced.

To quantify the fitting accuracy, each analytical profile is compared with its corresponding fitted dataset, from which we calculate the fitting residuals $\theta _{\textit{data}} - \theta _{\textit{fit}}$ and plot together with the corresponding analytical profiles in figures 6, 9, 12 and 15. In addition, we calculate the relative $L_2$ errors, defined by

(5.7) \begin{align} L_2 = \sqrt {\frac {\sum _{i=1}^N \left ( \phi _i^{\textit{DNS}} - \phi _i^{\textit{fit}} \right )^2}{\sum _{i=1}^N \left ( \phi _i^{\textit{DNS}} \right )^2}}, \end{align}

the root-mean-square (r.m.s.) residual $R_2$ , defined by

(5.8) \begin{align} R_2 = \sqrt {\frac {1}{N} \sum _i \left(\phi _i^{\textit{DNS}} - \phi _i^{\textit{fit}} \right)^2}, \end{align}

and maximum pointwise deviations, defined by

(5.9) \begin{align} E_{max } = \max \left |\phi _i^{\textit{DNS}} - \phi _i^{\textit{fit}}\right |. \end{align}

obtained for each fitted profile are summarised in table 6. Here, $N$ denotes the number of data points within the prescribed evaluation interval for the selected profile. The reported values demonstrate consistently small errors across all four boundary-condition configurations. The relative $L_2$ errors remain below approximately $3\,\%$ , the $R_2$ values are all below $2.2\,\%$ and the $E_{max }$ are generally below $3.7\,\%$ . Considering that the analytical profiles are derived from reduced boundary-layer equations rather than empirical interpolations, these results demonstrate that the proposed framework captures the essential structure of the DNS and experimental profiles with high accuracy.

Table 6.

Relative $L_2$ errors, r.m.s. residuals $R_2$ and maximum pointwise deviations $E_{max }$ for the DNS or experimental data shown in previous figures. The metrics are calculated for each profile over the corresponding range of $\xi$ given in the table. The values reported in each panel represent the worst-case profile.

Table 6. Long description

A table with six columns and eight rows. The columns are labeled Configuration, Figures, Range of ξ, L2, R2, and Emax. The rows are labeled with Roman numerals I to IV and include specific figures and ranges of ξ. Configuration I includes Figure 6 with a range of ξ from 0 to 4, L2 value of 0.0087, R2 value of 0.0074, and Emax value of 0.0148. Configuration II includes Figures 9(ai), 9(bi), and 9(ci) with ranges of ξ from 0 to 4, 0 to 4, and 0 to 2.5 respectively, L2 values of 0.0201, 0.0245, and 0.0296, R2 values of 0.0145, 0.0177, and 0.0205, and Emax values of 0.0303, 0.0332, and 0.0317. Configuration III includes Figure 12 with a range of ξ from 0 to 4, L2 value of 0.0285, R2 value of 0.0215, and Emax value of 0.0370. Configuration IV includes Figures 15(ai), 15(bi), and 15(ci) with ranges of ξ from 0 to 3, 0 to 3, and 0 to 2 respectively, L2 values of 0.0219, 0.0251, and 0.0203, R2 values of 0.0126, 0.0169, and 0.0104, and Emax values of 0.0236, 0.0400, and 0.0231.

To evaluate the robustness of the fitted coefficients, each coefficient is independently perturbed by $\pm 5\,\%$ , while all remaining coefficients are kept fixed. Since $k_1 \approx 2.23\,k_2$ , only $k_2$ is independently perturbed in the sensitivity analysis. The resulting error ranges are summarised in table 7. Although the degree of sensitivity varies among the four configurations, the overall agreement remains good under all perturbations. These results indicate that the predictive capability of the proposed analytical framework is robust.

Table 7.

Sensitivity of the analytical temperature profiles to independent $\pm 5\,\%$ perturbations in the listed model coefficients. For Configurations I–IV, evaluations are performed using the datasets corresponding to figures 6, 9(ci), 12 and 15(ci), respectively. Each entry reports the range of $L_2$ , $R_2$ and $E_{max }$ resulting from the $\pm 5\,\%$ perturbations.

Table 7. Long description

The table presents the sensitivity analysis of analytical temperature profiles to perturbations in model coefficients for four configurations. It has ten rows and five columns. The columns are labeled Configuration, Coefficient, L2, R2, and Emax. The rows are labeled with Roman numerals I, II, III, and IV, each containing different coefficients and their corresponding error ranges. For Configuration I, the coefficients k2 and a have specific ranges for L2, R2, and Emax. Configuration II includes coefficients k2, a, and d with their respective error ranges. Configuration III lists coefficients k and a with their error ranges. Configuration IV includes coefficients k, a, and d with their error ranges. Each entry reports the range of L2, R2, and Emax resulting from the perturbations.

Finally, it should be emphasised that the fitting procedure determines only the numerical values of a small number of characteristic coefficients. The mathematical forms of the analytical solutions are uniquely prescribed by the governing equations together with the corresponding hydrodynamic and thermal boundary conditions. Consequently, the fitted coefficients merely adjust the characteristic scales of the boundary layer, whereas the analytical profile itself is determined by the underlying boundary-value problem.

5.8. Validation of self-similarity assumptions and streamwise scalings

Figure 17.

Budget of terms in the governing momentum equation (3.2): (a) horizontal momentum balance comparing the pressure gradient $\partial _x P$ , convective terms ( $U \partial _x U, V \partial _z U$ ) and viscous diffusion $-\nu \partial _z^2 U$ ; (b) vertical momentum balance between the vertical pressure gradient $-\partial _z P$ and buoyancy profile $T$ . Profiles are horizontally averaged over the region $0.1 \le x/l_0 \le 0.9$ at $ \textit{Ra} = 10^{11}$ with no-slip boundary condition.

Figure 17. Long description

Panel A: A line graph shows the horizontal momentum balance comparing the pressure gradient, convective terms, and viscous diffusion. The x-axis represents the values of the horizontal momentum-balance terms in units of 10^-3, and the y-axis represents the normalized height (H-z)/H. The red line represents the pressure gradient, the yellow dashed line represents the convective term U∂xU, the green dotted line represents the convective term V∂zU, and the blue dashed line represents the viscous diffusion term -ν∂^2zU. Panel B: A line graph shows the vertical momentum balance between the vertical pressure gradient and buoyancy profile. The x-axis represents the values of the vertical momentum-balance terms, and the y-axis represents the normalized height (H-z)/H. The red line represents the vertical pressure gradient -∂zP, and the blue dashed line represents the buoyancy profile T.

Figure 17 demonstrates that the pressure-gradient and viscous-diffusion terms constitute the dominant balance throughout the analysed boundary layer, thereby providing quantitative justification for the reduced analytical formulation adopted in the present study. The quantified pressure gradient also provides the physical basis for the early departure from the linear viscous scaling discussed in (5.1)–(5.3). Unlike the asymptotic infinite- $ \textit{Pr}$ limit (Chiu-Webster et al. Reference Chiu-Webster, Hinch and Lister2008), the inertial terms remain finite rather than identically zero. Therefore, they are retained in the governing equations, allowing the unified analytical framework to capture the systematic differences in the boundary-layer structure and temperature profiles arising from different hydrodynamic and thermal boundary conditions.

For the boundary-layer ordinary differential equation (3.7) to admit self-similar solutions, all coefficients $a = \lambda (U_{0x}\lambda +U_0\lambda _x)/\nu$ , $b = -(U_{0x}\lambda ^2)/\nu$ and $c = -(\lambda ^2\rho ^{-1}P_x)/(\nu U_0)$ must be constant, a condition consistent with a laminar boundary layer. This ensures that the dependence of $\psi$ on the streamwise coordinate $x$ is expressed exclusively through the length scale $\lambda (x)$ , while the functional form of $\psi (\xi )$ remains invariant along the streamwise direction.

Figure 18.

Power-law dependence of physical parameters on streamwise position in log–log coordinates. (a) Scaling of $U_0=m_1(x/l_0)^{n_1}$ . In the interval $x/l_0 \in (0.1,\,1)$ , the exponent is $n_1=1.2$ for all $ \textit{Ra}$ , with $m_1=0.02$ , $0.0135$ , $0.028$ for no-slip $ \textit{Ra}=10^9$ , no-slip $ \textit{Ra}=10^{11}$ and free-slip $ \textit{Ra}=10^9$ , respectively. In the interval $x/l_0 \in (1,\,9)$ , the corresponding parameters are $n_1=0.25$ , $0.2$ , $0.18$ and $m_1=0.0247$ , $0.0165$ , $0.04$ . (b) Scaling of $P_x=m_2(x/l_0)^{n_2}$ ; the absolute value of $P_x$ is shown for representation on a logarithmic scale. In the interval $x/l_0 \in (0.1,\,1)$ , $n_2=1.15$ , $1.3$ , $1.95$ and $m_2=0.009$ , $0.0045$ , $0.01$ . In the interval $x/l_0 \in (1,\,9)$ , $n_2=-1$ , $-1.1$ , $-1.35$ and $m_2=0.035$ , $0.0165$ , $0.025$ . (c) Scaling of the local Reynolds number $ \textit{Re}=m_3(x/l_0)^{n_3}$ . In the interval $x/l_0 \in (0.1,\,1)$ , $n_3=1.1$ for all $ \textit{Ra}$ , with $m_3=2$ , $5$ , $1.2$ . In the interval $x/l_0 \in (1,\,9)$ , $n_3=0.5$ , $1.1$ , $0.8$ and $m_3=2.9$ , $3.8$ , $1.5$ .

A dimensional analysis of the coefficients reveals that $\lambda$ , $U_0$ and $P_x$ satisfy the scaling relations $U_{0x}\lambda ^2 \propto x^0$ , $U_0\lambda \lambda _x \propto x^0$ and $\lambda ^2 P_x/U_0 \propto x^0$ . We further hypothesise that these quantities each follow a power-law distribution in the streamwise direction, as demonstrated in figure 18(a,b). Previous studies (Mullarney et al. Reference Mullarney, Griffiths and Hughes2004; Shishkina Reference Shishkina2017) show that the velocity profiles attain a maximum in the wall-normal direction regardless of boundary conditions. In RBC, this maximum velocity is horizontally uniform, whereas in HC, the pressure-driven flow causes the peak to vary with streamwise position. Assuming $U_0(x) \propto x^n$ , the above-mentioned scaling relations imply $\lambda (x) \propto x^{(1-n)/2}$ and $P_x(x) \propto x^{2n-1}$ , indicating that the pressure $P$ scales quadratically with $U_0$ , consistent with the conclusion that $P \sim \rho U^2$ in RBC (Ching et al. Reference Ching, Leung, Zwirner and Shishkina2019).

Figure 19.

(a) Log–log plot of $\Delta T_b(x)$ against $x/l_0$ for various $ \textit{Ra}$ and boundary conditions. Lines denote reference power-law fits. (b) Streamwise distribution of the parameters $a$ and $d$ from the reduced temperature equation (3.9) at $ \textit{Ra}=10^{10}$ for no-slip and free-slip conditions. In both panels (a) and (b), the thermal boundary is adiabatic in the range $1 \leqslant x/l_0 \leqslant 9$ , and isothermal in the ranges $0 \leqslant x/l_0 \leqslant 1$ and $9 \leqslant x/l_0 \leqslant 10$ .

Figure 19. Long description

Panel A: A log-log plot shows the relationship between ΔTb and x/l0 for various Rayleigh numbers (Ra) and boundary conditions. The data series include no-slip conditions at Ra = 10^9 and Ra = 10^11, and free-slip conditions at Ra = 10^9. Reference power-law fits are denoted by lines. Panel B: A line graph illustrates the streamwise distribution of parameters a and d from the reduced temperature equation at Ra = 10^10 under no-slip and free-slip conditions. The thermal boundary is adiabatic in the range 1<x/l_0<9 and isothermal in the ranges 0<x/l_0<1 and 9<x/l_0<10.

Under these local power-law assumptions, the coefficients $a$ and $d$ appearing in the reduced temperature equation (3.9) become independent of the streamwise coordinate to leading order. Specifically, the term $\Delta T_b(x)$ in the coefficient $d$ is found to follow an approximate power law over the adiabatic self-similar interval, $\Delta T_b(x) \propto x^{p}$ . It is shown in figure 19(a). Therefore, $({\mathrm d(\Delta T_b)})/({\Delta T_b\,\mathrm dx}) \propto x^{-1}$ . Using $\lambda ^2 U_0 \propto x$ , we obtain $d \propto x^0$ . Accordingly, the streamwise dependence of the effective coefficients $a$ and $d$ is examined in figure 19(b).

Since the streamfunction solution varies with the velocity boundary condition and direct solution of the governing equations is analytically intractable, we propose a flow model comprising vertically oriented plumes and horizontally oriented jets. Differential heating induces plume structures analogous to those in RBC (Sanmiguel Vila et al. Reference Sanmiguel Vila, Discetti, Carlomagno, Astarita and Ianiro2016), which produce density variations and thereby a source of buoyancy. The resulting pressure gradient drives horizontal flow, acting as the momentum source for the jet. A jet is typically defined as fluid ejected from a nozzle and allowed to develop freely. Although no physical nozzle exists in the HC system, the horizontal pressure gradient can be regarded as an equivalent distributed momentum source that drives jet-like flow. This interpretation is supported by the observation that the pressure gradient is confined to an extremely thin region and primarily governs the velocity boundary-layer development, as shown in figure 17. The bulk-region flow, however, deviates from jet dynamics owing to thermal homogenisation.

It is worth noting that transitions between different boundary conditions induce localised deviations in the scaling of physical parameters. Before the jet develops into its self-similar mode, there is a flow-establishment zone. As shown in figures 18 and 19, the assumed power-law relationship remains valid only within a limited streamwise domain. As the flow approaches the cold plate, the intensifying turbulence triggers a transition to turbulent plume dynamics, causing the pressure distribution to depart from its initial behaviour. This accounts for the observed velocity enhancement and reduction in thermal boundary-layer thickness near the cold plate. To exclude flow-establishment effects and mitigate edge effects from the heated and cooled plates, all data are drawn from the isothermal interval $0.1 \leqslant x/l_0 \leqslant 0.9$ and the adiabatic interval $2 \leqslant x/l_0 \leqslant 8$ .

For an ideal jet, self-similarity requires linear spreading and momentum conservation implies that $U_0$ decays with streamwise distance. In the HC system studied here, however, the velocity increases with streamwise distance. By defining the local Reynolds number as $ \textit{Re}(x)=U_0(x)\,\sigma (x)/\nu$ , as shown in figure 18(c), it is found that $ \textit{Re}(x) \sim x^{1.1}$ in the isothermal region, whereas $ \textit{Re}(x) \sim x^{\gamma }$ in the adiabatic region, where $\gamma \geqslant 0.5$ is consistent with the scaling of a turbulent plane jet. These distinct scalings indicate differences in the lateral momentum transport. Consequently, the HC flow is best described as jet-like in its self-similar structure, rather than as a true jet.

5.9. Influence of surface slip on the vertical temperature structure

Figure 20.

Mean temperature profiles under different velocity boundary conditions in the adiabatic region at $ \textit{Ra}=10^{11}$ . (a) Vertical temperature distributions at several streamwise positions. (b) Normalised near-surface temperature profiles: the solid curve is given by (4.7) with the same parameters as in figure 9(a); the dotted curve is given by (5.6) with parameters $F(\alpha ,\beta ,\gamma )=F(0.5,\,-112.2,\,0.5)$ . The inset shows fitting functions $\theta (\xi )=1.05\,\xi ^2$ (dash-dotted line) and $\theta (\xi )=1.9\,\xi ^3$ (dashed line).

Despite the distinct velocity boundary conditions, the normalised temperature profiles discussed in §§ 5.4 and 5.6 exhibit similar shapes. To further examine how the velocity boundary condition modulates the temperature distribution in the adiabatic region, we fix the plate boundary conditions as no-slip and isothermal. Figure 20(a) presents the vertical temperature profiles $T((H-z)/H)$ at several streamwise locations within the adiabatic region. The differences are concentrated in the near-surface transition layer, with little variation in the bulk. Under the free-slip condition, the temperature departs from its near-surface value and approaches the bulk over a comparatively short depth with a gradual curvature; the transition is sharper under the no-slip condition. In addition, the near-surface temperature profiles at different $x/l_0$ are more tightly clustered under free-slip than under no-slip, indicating that the free-slip condition tends to smooth the streamwise variation of the local temperature and renders the vertical profiles less sensitive to streamwise position. This behaviour is associated with the streamwise evolution term $\mathrm{d}(\Delta T_b)/\mathrm{d}x$ in the coefficient $d$ defined in (3.10).

Figure 20(b) shows the normalised temperature profiles within the near-surface region for both velocity boundary conditions. The adiabatic thermal boundary condition imposes zero flux, $\theta '(0)=0$ , which reduces the sensitivity of the temperature profile to the detailed value of the tangential velocity at the wall. Furthermore, the outer portions of the streamfunction retain a jet-like structure in both cases, contributing to the similarity of the normalised profiles. Nevertheless, the velocity boundary condition still affects the curvature of $\theta (\xi )$ across the transition from the near-surface region to the bulk and, in particular, the depth at which the profile departs from the surface value.

Compared with the isothermal temperature equation (4.3), the adiabatic governing equation (3.9) contains an additional term $\text{d}\,\psi '(1-\theta )$ , which determines how the profile approaches the bulk state. The two velocity boundary conditions therefore differ primarily in how the near-surface tangential motion is initiated, and this distinction enters the temperature equation through $\psi '(\xi )$ :

(5.10) \begin{align} \psi '_{FS}(\xi )\approx 1+\mathcal{O}(\xi ^2)\quad (\xi \to 0) \ , \qquad \psi '_{NS}(\xi )\approx k_1\xi +\mathcal{O}(\xi ^3)\quad (\xi \to 0) \ , \end{align}

where the subscripts FS and NS denote free-slip and no-slip, respectively. Combined with the adiabatic condition $\theta '(0)=0$ , a Taylor expansion of $\theta$ about $\xi =0$ reveals that the near-surface departure of $\theta$ begins at second order under free-slip, $\theta \sim \xi ^2$ , and at third order under no-slip, $\theta \sim \xi ^3$ , consistent with the inset of figure 20(b).

From the perspective of transport, the contrast is set by the near-surface onset of advection. Under no-slip conditions, the tangential velocity vanishes at the boundary, so the advective contribution is negligible near the surface and becomes significant only beyond a certain depth. Under free-slip conditions, the tangential velocity is unconstrained at the boundary, so advective driving is effective immediately, causing the temperature profile to depart from its surface value at a shallower depth. The distinct velocity structures are also reflected in the functional form of the analytical solutions. In the no-slip case, the additional linear factor $k_1\xi$ breaks the algebraic closure between $\tanh (\xi )$ and $\operatorname {sech}^2(\xi )$ , leading to a more general singularity structure after transformation and to a confluent Heun form. Physically, this is consistent with the no-slip condition introducing an additional distance factor that delays effective thermal–momentum coupling away from the boundary, even though the normalised profiles appear similar over the range of $\xi$ of primary interest.