2. Numerical set-up, simulations and control and response parameters

2.1. Dynamical equations and control parameters

We consider the Navier–Stokes equations subject to the Oberbeck–Boussinesq approximation, where changes in density $\rho$ are only relevant in the buoyancy and a linear equation of state relates the density changes to temperature $T$. These equations read $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0$ and

(2.1)\begin{gather} \partial_t \boldsymbol{u} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} ={-}\frac{1}{\rho_0}\boldsymbol{\nabla} p + \nu \nabla^{2} \boldsymbol{u} + g \alpha T\hat{\boldsymbol{y}}, \end{gather}

(2.2)\begin{gather}\partial_t T + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \kappa \nabla^{2} T, \end{gather}

where $\boldsymbol {u}=(u,v,w)$ is the velocity field, $p$ the kinematic pressure, $\nu$ kinematic viscosity, $\kappa$ the molecular diffusivity of heat, $g$ gravitational acceleration, $\alpha$ the thermal expansion coefficient and $\rho _0$ a reference density. We solve these equations in a vertical channel domain between two no-slip, impermeable, isothermal walls. These walls are separated by a distance $H$ and the temperature difference between them is $\varDelta$. As in Ng et al. (Reference Ng, Ooi, Lohse and Chung2015) and shown in figure 1, we consider a domain of length $8H$ in the vertical ($y$) and length $4H$ in the spanwise ($z$) direction, and impose periodic boundary conditions on $\boldsymbol {u}$, $p$ and $T$ in these directions, $y$ and $z$. In a convective system, we can scale the velocity by the free-fall velocity $U_T=\sqrt {g\alpha {\rm \Delta} H}$ so that the dynamics of the system is solely determined by the Rayleigh and Prandtl numbers

(2.3a,b)\begin{equation} \mbox{{Ra}} = \frac{g\alpha H^{3} \varDelta}{\nu\kappa},\quad {Pr} = \frac{\nu}{\kappa}. \end{equation}

Figure 1. (a) A schematic of the simulation domain. (b,c) Mean profiles of the vertical velocity and the temperature for $\mbox {{Ra}}=10^{8}$ and a range of ${Pr}$. Recall that the mean profiles are anti-symmetric such that $\bar {v}(x) =-\bar {v}(H-x)$.

These are the only control parameters of the system, aside from parameters characterising the geometry of the flow domain. Their ratio $\mbox {{Gr}}=\mbox {{Ra}}/{Pr}=g\alpha H^{3} \varDelta /\nu ^{2}$ is also called the Grashof number.

The governing equations (2.1) and (2.2) are solved numerically using a second-order finite difference scheme for spatial derivatives and a third-order Runge–Kutta scheme for time stepping, as described in Verzicco & Orlandi (Reference Verzicco and Orlandi1996) and van der Poel et al. (Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015). For high values of ${Pr}$, the temperature field must be resolved at smaller scales than the velocity field because the temperature field diffuses on the time scale of the order of ${Pr}^{-1}$ compared with the velocity field. We therefore also use the multiple-resolution technique of Ostilla-Monico et al. (Reference Ostilla-Monico, Yang, van der Poel, Lohse and Verzicco2015) to evolve the scalar $T$ on a refined grid. Interpolation between the two grids is achieved through a four-point Hermitian method. Grid stretching is also implemented in the wall-normal ($x$) direction using a clipped Chebyshev-type clustering. Uniform grid spacing is used in the $y$ and $z$ directions, and the base grid of all simulations are resolved down to a factor of 2 times the Kolmogorov scale. The refined grid is such that the wall-normal grid spacing satisfies $\varDelta _x <0.5L_B$ at the boundaries, and the grid spacing in the bulk satisfies $\varDelta _{x,y,z} < 4.5L_B$, where $L_B = (\nu \kappa ^{2}/\varepsilon )^{1/4}$ is the Batchelor scale.

The range of dimensionless control parameters simulated is shown in table 1. Simulations at $\mbox {{Ra}}=10^{6}$ are initialised using the laminar, purely conductive solution of Batchelor (Reference Batchelor1954) with the addition of small amplitude random noise to trigger a transition to turbulence. Simulations at higher $\mbox {{Ra}}$ are initialised using the final state of the simulation at $\mbox {{Ra}}=10^{6}$ and ${Pr}=1$, interpolated onto a new grid. Each computation is performed for at least $300$ free-fall times, where $H/U_T$ is the free-fall time unit. We average statistics over the last $250$ time units once the system has reached a statistically steady state.

Table 1. Overview of the dimensionless parameters and grid resolutions used in the numerical simulations. Grid resolutions are listed here for the cases at highest $\mbox {{Ra}}$, and we distinguish between the base grid used to evolve the velocity and the refined grid used to evolve the temperature field.

2.2. Response parameters and theoretical scaling laws

Before presenting the results of the simulations, we now provide an overview of the key quantities of interest and existing theoretical frameworks used for their prediction.

Understanding how the global horizontal heat transport in vertical convection depends on the control parameters of (2.3a,b) is vital for many applications. Varying the control parameters also leads to changes in the peak velocity of the rising flow and the mean shear stress on the boundary. These can be quantified through the following dimensionless response parameters: the Nusselt number, the Reynolds number and the shear Reynolds number

(2.4a–c)\begin{equation} \mbox{{Nu}} = \frac{H q_T}{\kappa \varDelta},\quad {Re} = \frac{V_{max} H}{\nu},\quad {Re}_\tau = \frac{V_\ast H}{\nu}, \end{equation}

where $q_T=\kappa |{\textrm {d}}\overline {T}/\textrm {d}x|_{wall}$ is the horizontal heat flux, $V_{max}$ is the peak value of the time- and spatially averaged vertical velocity $\overline {v}(x)$ and $V_\ast = \sqrt {\tau _w/\rho _0}$ is the friction velocity associated with the mean wall shear stress $\tau _w=\mu \,{\textrm {d}}\overline {v}/\textrm {d}x|_{wall}=\rho _0 {V_\ast }^{2}$.

In turbulent convection, many studies follow the so-called ‘classical’ regime as a theoretical starting point. This regime relies on the assumption that the thermal driving is sufficiently strong such that the heat flux becomes independent of the plate separation $H$. Assuming a power-law relation between $\mbox {{Nu}}$ and the Rayleigh number, dimensional analysis (e.g. Turner Reference Turner1979) then requires the scaling $\mbox {{Nu}} \sim \mbox {{Ra}}^{1/3}f({Pr})$. This has been consistent with various experiments up to $Ra = 10^{12}$ (Warner & Arpaci Reference Warner and Arpaci1968; Tsuji & Nagano Reference Tsuji and Nagano1988) and is often provided in engineering reference texts such as Holman (Reference Holman2010).

However, recent analysis of numerical simulations by Ng et al. (Reference Ng, Ooi, Lohse and Chung2017) suggests that a power-law description may be insufficient and that the ‘classical’ scaling does not accurately describe the data even in this range. Furthermore, there are open questions regarding the relevant scaling at even higher $\mbox {{Ra}}$, at which precise, controlled experiments and numerical simulations are extremely difficult to perform. Finally, the Prandtl number dependence has hardly been addressed.

One important application for boundary layers in VC is to predict the dissolution or melting of a vertical ice face in the ocean. This is why parameterising the heat and salt fluxes is crucial. In regional ocean models, the heat flux through the turbulent boundary layer at such locations is often parameterised by invoking the heat flux balance $q_T \sim \textrm {velocity} \times \textrm {temperature change}$. Following Holland & Jenkins (Reference Holland and Jenkins1999), the parameterisation for the horizontal heat flux takes the form

(2.5)\begin{equation} q_T = C_T C_D^{1/2} U (T - T_b), \end{equation}

where $U$ is the vertical velocity of the rising plume, and $T-T_b$ is the temperature difference between the ocean and the ice boundary.

Taking $U=V_{max}$ and $T-T_b = \varDelta /2$, we note that the drag coefficient $C_D$ and ‘transfer coefficient’ $C_T$ from (2.5) are fully determined by the response parameters of (2.4a–c) through

(2.6a,b)\begin{equation} C_D = \left(\frac{V_\ast}{V_{max}}\right)^{2} = \frac{{{Re}_\tau}^{2}}{{{Re}}^{2}}, \quad C_T = \frac{2q_T}{V_\ast \varDelta} = \frac{2\mbox{{Nu}}}{{Re}_\tau {Pr}}. \end{equation}

Accurate scaling laws for the quantities in (2.4a–c) are thus crucial for determining $C_D$ and $C_T$. The transfer coefficient $C_T$ is equivalent to a modified Stanton number where $V_\ast$ is used for the velocity scale. In the ice–ocean literature, the transfer coefficient is often denoted $\varGamma _T$ although we use $C_T$ here to avoid confusion with the aspect ratio $\varGamma$ used throughout literature on convection. Both $C_D$ and $C_T$ are typically set to constant values in melt parameterisations (see e.g. Jackson et al. Reference Jackson, Nash, Kienholz, Sutherland, Amundson, Motyka, Winters, Skyllingstad and Pettit2020) based on the reasoning that the boundary layers in ice–ocean applications are strongly shear driven, and are in accordance with the classical results of Kader & Yaglom (Reference Kader and Yaglom1972). However, recent analysis by Malyarenko et al. (Reference Malyarenko, Wells, Langhorne, Robinson, Williams and Nicholls2020) of ice shelf observations suggests that the Reynolds numbers may not always be large enough to justify this shear-driven boundary layer assumption. An equivalent equation to (2.5) is often used to parameterise the salt flux, where $C_D$ keeps the same value, but $C_T$ is reduced to reflect its dependence on the Schmidt number.

For $C_T$, $C_D$ being constant, the dimensionless form of (2.5) is $\mbox {{Nu}}\sim {Re}{Pr}$. Such a scaling is reminiscent of the ‘ultimate’ or ‘diffusion-free’ scaling hypothesised for Rayleigh–Bénard convection (RBC) at very high $\mbox {{Ra}}$ (e.g. Kraichnan Reference Kraichnan1962; Spiegel Reference Spiegel1971; Lohse & Toschi Reference Lohse and Toschi2003; Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009). In that case, the heat flux is assumed independent of the molecular diffusivities $\nu$ and $\kappa$, such that dimensional analysis implies $\mbox {{Nu}} \sim (RaPr)^{1/2}$. In physical terms, this regime is associated with a dominant large-scale circulation that leads to shear-driven turbulent boundary layers. The dominant mean flow arising in VC is analogous to such a coherent large-scale circulation (Shishkina & Horn Reference Shishkina and Horn2016). RBC provides a useful comparison with VC thanks to its identical geometry (except for the direction of gravity) and its dependence on the same control parameters.

In the case of RBC, a unifying theory describing the transitions between various regimes in RBC was proposed by Grossmann & Lohse (Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001). This theory has shown excellent agreement with subsequent experimental and numerical investigations over a large range of $\mbox {{Ra}}$ and ${Pr}$ (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Stevens et al. Reference Stevens, van der Poel, Grossmann and Lohse2013). Although VC lacks the global relation between the Nusselt number and the mean dissipation rate of kinetic energy required to close the equations corresponding to those of the Grossmann–Lohse (GL) theory, it remains appealing to search for parallels between RBC and VC to understand how the heat flux can be parameterised as the boundary layers evolve. Wells & Worster (Reference Wells and Worster2008) applied ideas from the GL theory about boundary layer transition to geophysical-scale convection at a vertical wall, and Ng et al. (Reference Ng, Ooi, Lohse and Chung2015, Reference Ng, Ooi, Lohse and Chung2017) considered how changes in the boundary layer structure relate to an increased bulk contribution to turbulent dissipation. However, these studies left the issue of ${Pr}$-dependence largely unresolved. In this study, we aim to gain insight into how the Prandtl number affects (a) the scaling of the above response parameters in the currently accessible range of $\mbox {{Ra}}$, and (b) any subsequent transition in the nature of the boundary layers.

3. Flow visualisation

To illustrate how the boundary layers in the flow change with ${Pr}$ and $\mbox {{Ra}}$, we present a snapshot of the local dimensionless vertical shear stress $\hat {\tau }$ and heat flux $\hat {q}$ at the heated wall $x=0$ in figure 2. These quantities are defined as

(3.1a,b)\begin{equation} \hat{\tau}(y,z) = \frac{H}{U_T} \left.\frac{\partial v}{\partial x}\right|_{x=0}, \quad \hat{q}(y,z) ={-}\frac{H}{\varDelta} \left.\frac{\partial T}{\partial x}\right|_{x=0}. \end{equation}

We note that averaging $\hat {q}$ over the plane and over time gives the Nusselt number $\langle \hat {q}(y,z,t)\rangle _{y,z,t}=\mbox {{Nu}}$. In this sense, $\hat {q}$ can be thought of as a ‘local and instantaneous Nusselt number’. The snapshots, taken at the end time of each simulation, highlight the striking localisation of the heat flux at the wall. Consistent with the analysis of Pallares et al. (Reference Pallares, Vernet, Ferre and Grau2010), the regions of strongest heat flux (being the dark patches in the panels of $\hat {q}$) are frequently co-located with instantaneous flow reversals (evidenced by white and blue patches appearing in the panels for $\hat {\tau }$).

Figure 2. Final-time snapshots of the dimensionless horizontal shear stress $\hat {\tau }$ and the local dimensionless heat flux $\hat {q}$ at the heated wall $x=0$ for a range of $\mbox {{Ra}}$ and ${Pr}$. It can be seen how large ${Pr}$ smooths the fields, even at a large $\mbox {{Ra}}=10^{8}$.

Figure 2(a–c) highlights the effect of increasing ${Pr}$ in this set-up while keeping $\mbox {{Ra}}$ fixed (in this case at $10^{8}$). As seen from the colour scales, although the range of the local Nusselt number $\hat {q}$ remains similar as ${Pr}$ increases, a significant decrease in the mean dimensionless shear stress is observed at high ${Pr}$. This is due to a drop in the Grashof number $\mbox {{Gr}}$ as ${Pr}$ is increased for fixed $\mbox {{Ra}}$. The Grashof number quantifies the ratio of buoyancy effects to viscosity, and is analogous to a squared Reynolds number based on the free-fall velocity.

This analogy with the Reynolds number provides some further intuition for the snapshots of figure 2, where a much wider range of length scales can be observed in the $\hat {\tau }$ field for the high $\mbox {{Gr}}$ snapshot of (a) compared with the lower $\mbox {{Gr}}$ snapshot of (c). By contrast, comparing panels (b,d) allows us to visualise the effect of changing ${Pr}$ while keeping $\mbox {{Gr}}$ fixed. Qualitatively the structures in both the $\hat {\tau }$ and $\hat {q}$ snapshots appear similar. However, the mean values of both quantities vary as ${Pr}$ increases. To obtain a more quantitative evaluation of the boundary layer structures and to more quantitatively extract length scales, we have also calculated the relevant power spectra for each of the simulations shown in figure 2. These results (not shown here) emphasise the similarity of structures for constant $\mbox {{Gr}}$ at the walls, although this similarity does not extend outside of the viscous boundary layer.

Compared with RBC, where large-scale thermal structures do not exhibit a preferred direction, the mean shear at the wall in VC introduces significant anisotropy to the wall structures. Streaky structures elongated in the vertical ($y$) direction are prominent in figure 2, similar to those seen in the sheared RBC set-up of Blass et al. (Reference Blass, Tabak, Verzicco, Stevens and Lohse2021). Furthermore, in that study, an increased ${Pr}$ (for fixed $\mbox {{Ra}}$ and ${Re}$) was found to enhance momentum transport from the walls, allowing the wall shear to affect the flow structures in the bulk more easily. However, as mentioned in § 2, the wall shear in VC is not pre-determined and instead arises as a response parameter of the system.

The snapshots of figure 2 highlight the complex multi-parameter dependence in the VC set-up. Indeed, the simple analogy between the Grashof number and the square of the Reynolds number should not be overstated. As shown in figure 1(b), the peak value of the time-averaged vertical velocity does not simply scale with the free-fall velocity $U_T$, but varies depending on ${Pr}$. In the following section, we shall investigate the multi-parameter dependence more quantitatively by identifying scaling relations for key response parameters of the system.

4. Heat flux and Reynolds number parameterisation

In table 2 we report the observed $\mbox {{Ra}}$- and ${Pr}$-dependence of the response parameters from (2.4a–c) and (2.6a,b) in our simulations. An effective power-law dependence is assumed and two-parameter linear regression is used to obtain the effective scaling exponents. Precisely, we compute $\boldsymbol {b} = X^{-1}\boldsymbol {y}$, where $\boldsymbol {b}=(b_1, b_2, b_3)^\textrm {T}$ and

(4.1a–d)\begin{equation} x_{i 1} = \log \mbox{{Ra}}_i, \quad x_{i 2} = \log {Pr}_i, \quad x_{i 3} = 1, \quad y_i = \log R_i, \quad i=1,\ldots,n \end{equation}

are constructed from the $n$ simulations for each response parameter $R$, giving a linear fit $R=\mbox {{Ra}}^{b_1} {Pr}^{b_2} 10^{b_3}$. We calculate the uncertainty of the power-law exponents $b_1$ and $b_2$ through the variance matrix of $\boldsymbol {b}$ given by $V=\sigma ^{2} (X^\textrm {T} X)^{-1}$, where $\sigma ^{2}$ is the variance of $\boldsymbol {y} - X \boldsymbol {b}$. The standard deviations of the slopes, given by $\sqrt {v_{11}}$ and $\sqrt {v_{22}}$ are presented in table 2.

Table 2. Observed effective scalings laws for various dimensionless response parameters. Only simulations with ${Re}>150$ are included in the linear regression. The uncertainty shown is the standard deviation of the estimated slopes, as described in the text of § 4. Theoretical scaling relations for laminar VC and turbulent RBC from Shishkina (Reference Shishkina2016) for VC and Grossmann & Lohse (Reference Grossmann and Lohse2000) for RBC in the so-called $\mathrm {IV}_u$ are provided for comparison. ${Re}_\tau$ is calculated for these scaling relations using the similarity variable of Shishkina (Reference Shishkina2016) and using the Blasius drag law $C_D\sim Re^{-1/2}$ for the GL theory.

The Nusselt number is consistent with the theoretical scaling relation $\mbox {{Nu}}\sim \mbox {{Ra}}^{1/3}f({Pr})$ that arises when the heat flux is assumed to be independent of the plate separation (Malkus Reference Malkus1954). Ng et al. (Reference Ng, Ooi, Lohse and Chung2017) suggested that for ${Pr}\approx 1$, a regime transition to a shear-dominated boundary layer is underway at $\mbox {{Ra}}=10^{9}$, but following Grossmann & Lohse (Reference Grossmann and Lohse2000), this transitional $\mbox {{Ra}}$ can be expected to increase with ${Pr}$, as the smaller Reynolds number stabilises the flow. Our results contrast with the effective scaling laws for laminar VC derived by Shishkina (Reference Shishkina2016), where $\mbox {{Nu}}\sim \mbox {{Ra}}^{1/4}$ and ${Re}\sim \mbox {{Ra}}^{1/2}{Pr}^{-1}$ for ${Pr}\gg 1$. This difference is to be expected since our set-up is far from the laminar state for which the scaling laws have been observed to hold (e.g. Wang et al. Reference Wang, Liu, Verzicco, Shishkina and Lohse2021).

In figure 3 we plot $\mbox {{Nu}}$ against both $\mbox {{Ra}}$ and the shear Reynolds number ${Re}_\tau$. Figure 3(a) highlights the weak dependence of $\mbox {{Nu}}$ on ${Pr}$, with higher ${Pr}$ typically reducing $\mbox {{Nu}}$ for a fixed value of $\mbox {{Ra}}$. Note that a simple, single power-law fit is unlikely to adequately describe the heat transfer outside of the currently accessible parameter range. Even within the data presented here, the ${Pr}=1$ cases appear to trend downwards relative to the $\mbox {{Ra}}^{1/3}$ line on figure 3(a) at higher values of $\mbox {{Ra}}$. This observation is consistent with Ng et al. (Reference Ng, Ooi, Lohse and Chung2017), who attribute the decrease to a lower heat flux contribution from regions of weak shear. Later in this section, and in Appendix A, we shall discuss at which parameter values we may expect a transition to shear-driven turbulent boundary layers and how this would affect the scaling of the Nusselt number.

Figure 3. Nusselt number against (a) Rayleigh number (compensated by $\mbox {{Ra}}^{1/3}$), and (b) against shear Reynolds number (compensated by ${Pr}^{1/3}$).

Against ${Re}_\tau$ in figure 3(b), we obtain a reasonable collapse for $\mbox {{Nu}}$ by scaling with ${Pr}^{1/3}$ and observe a scaling close to $\mbox {{Nu}}\sim {Re}_\tau {Pr}^{1/3}$. Since this is consistent with the high ${Pr}$ limit of passive heat transport in turbulent boundary layers from Kader & Yaglom (Reference Kader and Yaglom1972), we are motivated to compare with passive scalar transport in other turbulent flows. A recently proposed a scaling theory for passive scalar transport in plane Couette flow suggests that $\mbox {{Nu}}\sim {Re}_\tau ^{6/7}{Pr}^{1/2}$ (Yerragolam, Stevens, Verzicco, Lohse & Shishkina, personal communication). This somewhat contrasts with the ${Pr}^{1/3}$ collapse observed in figure 3(b), although the higher ${Re}_\tau$ values of our data do exhibit a local scaling exponent less than one and close to $6/7$.

We note that the Reynolds number scaling in table 2 is close to that reported by Lam et al. (Reference Lam, Shang, Zhou and Xia2002) from experiments of RBC with a range of large Prandtl numbers. Lam et al. (Reference Lam, Shang, Zhou and Xia2002) suggested that their results were consistent with the theoretical scaling relation ${Re}\sim \mbox {{Ra}}^{4/9}{Pr}^{-2/3}$ proposed for the regime ($IV_u$) associated with ${\mbox {{Nu}}\sim \mbox {{Ra}}^{1/3}}$ in the ‘GL theory’ of Grossmann & Lohse (Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001), although Lam et al. (Reference Lam, Shang, Zhou and Xia2002) acknowledge that this measured effective ${Pr}$ exponent shows a relatively large deviation from the theory. Furthermore, these deviations varied depending on the definition of the Reynolds number inferred from their experiments. We note that the ${Re}\sim \mbox {{Ra}}^{4/9}{Pr}^{-2/3}$ scaling can also be derived from dimensional analysis by assuming that the vertical velocity $V_{max}$ is solely determined by the buoyancy flux per unit area $\varPhi =g\alpha q_T$ and the plate separation $H$ (as in the ‘outer’ scaling of George & Capp Reference George and Capp1979), and also assuming the Malkus (Reference Malkus1954) scaling $\mbox {{Nu}}\sim \mbox {{Ra}}^{1/3}$. As seen from table 2, this ${Re}$ scaling does not perfectly capture the observed data, and we cannot rule out the effect of multiple regimes on the effective scaling exponent, as in the GL theory for RBC. More work is needed to provide a theoretical understanding for these results.

As highlighted by McConnochie & Kerr (Reference McConnochie and Kerr2017), the scaling relation $\mbox {{Nu}} \sim \mbox {{Ra}}^{1/3}f({Pr})$ implies a dimensional form for the heat flux that scales as $F_T\sim {\rm \Delta} T^{4/3}$ for fixed fluid properties. The heat flux is therefore independent of the bulk velocity $V_{max}$, making the shear-based model of (2.5) an inappropriate parameterisation for this regime. Indeed, as shown in figure 4, we observe significant variation in the drag coefficient $C_D$ with both $\mbox {{Ra}}$ and ${Pr}$. In all cases we find a value much larger than the high-${Re}$ limit of $C_D=2.5\times 10^{-3}$, as used by Holland & Jenkins (Reference Holland and Jenkins1999). However, the scaling observed for the transfer coefficient $C_T \approx 0.1 {Pr}^{-2/3}$ is consistent with the values used for parameterising heat and salt fluxes in that work and subsequent melting studies. Using the definition from e.g. (2.6a,b), we can express this result in terms of the Nusselt number as $\mbox {{Nu}} \sim {Re}_\tau {Pr}^{1/3}$ or with dimensional quantities as $q_T \sim {Pr}^{-2/3} V_\ast \varDelta$.

Figure 4. Drag coefficient $C_D$ and transfer coefficient $C_T$ defined in (2.6a,b). Dashed lines in (a) use the exponents obtained from the linear regression in table 2.

It may be tempting to associate the scaling $\mbox {{Nu}}\sim {Re}_\tau {Pr}^{1/3}$ with the appearance of turbulent boundary layers in the sense of Prandtl and von Kàrmàn, where log-law profiles appear in the mean velocity and temperature profiles. However, this is not the case for our simulations. In figure 5 we plot these mean profiles from the simulations at $\mbox {{Ra}}=10^{8}, 10^{9}$ with a logarithmic $x$-axis. From figure 5(a), it is clear that log layers are absent from the velocity profile. Indeed, we are far from the critical Reynolds number for transition to such a shear-driven boundary layer. As we explore in Appendix A, Rayleigh numbers above $10^{11}$ are likely to be necessary for this transition and such critical values only increase with ${Pr}$. By contrast, the temperature profiles of figure 5(b) appear consistent with logarithmic profiles. This observation is somewhat unsurprising, given the appearance of such profiles in the ‘classical’ regime of RBC by Ahlers et al. (Reference Ahlers, Bodenschatz, Funfschilling, Grossmann, He, Lohse, Stevens and Verzicco2012). A logarithmic profile in the temperature field does not imply the presence of a shear-driven turbulent boundary layer.

Figure 5. Mean profiles of (a) vertical velocity and (b) temperature on a logarithmic $x$ axis. The $x$-axis is scaled in terms of viscous wall units, such that $x^{+} = xV_\ast /\nu$. Vertical velocity is scaled by the friction velocity $V_\ast$, and temperature is scaled by the equivalent ‘friction temperature’ scale $T_\ast = q_T/V_\ast$. Solid lines denote simulations at $\mbox {{Ra}}=10^{8}$, whereas dotted lines represent the two simulations at $\mbox {{Ra}}=10^{9}$. The dashed black lines denote the linear profiles $\bar {v}=V_\ast x^{+}$ and $\bar {T} = T_\ast {Pr}x^{+}$ in (a) and (b), respectively. The inset in (b) is a zoom out of the main figure highlighting the results for ${Pr}=100$.

Holland & Jenkins (Reference Holland and Jenkins1999) associate the scaling relation $C_T \sim Pr^{-2/3}$ with the strong influence of a molecular sublayer where conduction is the dominant mechanism of heat transport. Motivated by this result, we proceed by investigating how the width of this boundary layer depends on the control parameters of the VC system.

5. Conductive thermal boundary layer

In a statistically steady state, the mean velocity and temperature profiles of the system satisfy

(5.1a,b)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}\kern0.7pt x} \overline{u^{\prime} \boldsymbol{u}^{\prime}} = \nu \frac{\mathrm{d}^{2} \bar{\boldsymbol{u}}}{\mathrm{d}\kern0.7pt x^{2}} + \bar{T} \hat{\boldsymbol{y}},\quad \frac{\mathrm{d}}{\mathrm{d}\kern0.7pt x} \overline{u^{\prime} T^{\prime}} = \kappa \frac{\mathrm{d}^{2} \bar{T}}{\mathrm{d}\kern0.7pt x^{2}}, \end{equation}

where an overbar denotes an average in $y$, $z$ and time. Incompressibility ensures that $\bar {u}\equiv 0$, so the mean velocity $\bar {\boldsymbol {u}}$ only has components in the wall-parallel directions. The second equation of (5.1a,b) implies that the heat flux at any wall-normal location must be constant, or in dimensionless terms

(5.2)\begin{equation} \mbox{{Nu}} = \frac{H}{\kappa\varDelta} \left(\overline{u^{\prime} T^{\prime}} - \kappa \frac{\mathrm{d}\bar{T}}{\mathrm{d}\kern0.7pt x}\right) = \mathrm{constant}. \end{equation}

Following Wells & Worster (Reference Wells and Worster2008) and in the spirit of Grossmann & Lohse (Reference Grossmann and Lohse2000), we divide the flow into thermal boundary layers, where the heat flux is dominated by molecular diffusion of the mean, and bulk regions, where the heat flux is due to the ‘wind’ of turbulence. Precisely, we define the conductive thermal boundary layer width $\delta _T$ as the wall-normal location where the conductive heat flux $-\kappa \,{\textrm {d}}\bar {T}/{\textrm {d}x}$ is equal to the turbulent heat flux $\overline {u^{\prime } T^{\prime }}$.

In RBC at moderate $\mbox {{Ra}}$, there is a general consensus from existing literature (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Ching et al. Reference Ching, Leung, Zwirner and Shishkina2019) that, scaling-wise, the thickness of the boundary layers follows a laminar-like scaling according to Prandtl, Blasius and Pohlhausen, that is

(5.3)\begin{equation} \frac{\delta_T}{H} \sim {Re}^{{-}1/2} f({Pr}) . \end{equation}

For VC, Ng et al. (Reference Ng, Ooi, Lohse and Chung2015) suggested the application of the same form as (5.3) at moderate Reynolds number, although only cited ‘fair’ agreement with their direct numerical simulations reporting an effective ${Re}$-exponent of $-0.60$. The scaling (5.3) is applicable in the case of a fully laminar flow as studied by Kuiken (Reference Kuiken1968), who derived an equivalent scaling of $\delta _T/H \sim \mbox {{Gr}}^{-1/4}$ in the limit of high ${Pr}$. The scaling laws $\mbox {{Nu}}\sim \mbox {{Ra}}^{1/4}$, ${Re}\sim \mbox {{Ra}}^{1/2}Pr^{-1}$ proposed by Shishkina (Reference Shishkina2016) are also consistent with (5.3).

From our new simulations, we find a collapse of the data such that $\delta _T/H \sim {Pr}^{-1/3} f({Re})$, as shown in figure 6. This ${Pr}$-dependence is well known from the similarity scaling of a laminar boundary layer at a horizontal wall (e.g. Schlichting & Gersten Reference Schlichting and Gersten2016), applied to the regimes of Grossmann & Lohse (Reference Grossmann and Lohse2000) where the thermal dissipation rate is dominated by boundary layer contributions. However, the ${Pr}^{-1/3}$ factor does not arise in the laminar solutions for VC considered by Kuiken (Reference Kuiken1968) and Shishkina (Reference Shishkina2016). The ${Pr}^{-1/3}$ scaling is often also observed in empirical data for turbulent flows (e.g. Kader Reference Kader1981). Indeed, rather than observing a laminar-like ${Re}^{-1/2}$ scaling, our data are more consistent with

(5.4)\begin{equation} \frac{\delta_T}{H} \sim {Re}^{{-}2/3} {Pr}^{{-}1/3} , \end{equation}

as shown in figure 6(b).

Figure 6. (a) Dimensionless conductive thermal boundary layer width $\delta _T/H$ against Reynolds number. (b) Plot of the same data compensated by ${Re}^{2/3} {Pr}^{1/3}$. (c) Measured sublayer Rayleigh number $\mbox {{Ra}}_\delta$ as a function of Prandtl number. Dashed lines in panel (a) mark the suggested ${Re}^{-2/3} {Pr}^{-1/3}$ scaling.

For the case where $V_{max}\sim U_T$, the scaling of (5.4) is equivalent to $\delta _T/H \sim \mbox {{Ra}}^{-1/3}$ and one can interpret the boundary layer width as being set by a critical Rayleigh number. This is the ‘buoyancy-controlled sublayer’ scaling as described by Wells & Worster (Reference Wells and Worster2008), similar to the marginally stable boundary layer argument of Malkus (Reference Malkus1954) for RBC. However, as we already mentioned earlier, $V_{max}$ does not simply scale with $U_T$ in our simulations. In figure 6(c), we plot the ‘sublayer Rayleigh number’ $\mbox {{Ra}}_\delta = g\alpha \delta _T^{3} \varDelta /\nu \kappa$, and find that this value is not constant, but instead strongly depends on the Prandtl number. Further studies are certainly needed to understand how to interpret these results. It remains an open question whether a ${Pr}$-dependent critical Rayleigh number is appropriate for limiting the conductive boundary layer width or whether the Reynolds number plays a more significant role. The addition of a spanwise mean flow to the system, forming a three-dimensional mixed convection set-up, would allow ${Re}$ and $\mbox {{Ra}}$ to be decoupled, and reveal the inherent parameter dependence of the boundary layer.