Boundary layers in turbulent vertical convection at high Prandtl number

Abstract

Many environmental flows arise due to natural convection at a vertical surface, from flows in buildings to dissolving ice faces at marine-terminating glaciers. We use three-dimensional direct numerical simulations of a vertical channel with differentially heated walls to investigate such convective, turbulent boundary layers. Through the implementation of a multiple-resolution technique, we are able to perform simulations at a wide range of Prandtl numbers ${Pr}$. This allows us to distinguish the parameter dependences of the horizontal heat flux and the boundary layer widths in terms of the Rayleigh number $\mbox {{Ra}}$ and Prandtl number ${Pr}$. For the considered parameter range $1\leq {Pr} \leq 100$, $10^{6} \leq \mbox {{Ra}} \leq 10^{9}$, we find the flow to be consistent with a ‘buoyancy-controlled’ regime where the heat flux is independent of the wall separation. For given ${Pr}$, the heat flux is found to scale linearly with the friction velocity $V_\ast$. Finally, we discuss the implications of our results for the parameterisation of heat and salt fluxes at vertical ice–ocean interfaces.

1. Introduction

When a fluid is heated from a side boundary, buoyancy drives a flow up the boundary via convection. The laminar flow along a heated surface has long been understood (Batchelor 1954; Kuiken 1968; Shishkina 2016) but there is no formal solution for the case where the flow becomes turbulent. This occurs when the Rayleigh number Ra of the flow is sufficiently high. In many environmental applications of this so-called vertical convection (VC), such as the flow in a cavity wall, high Rayleigh numbers imply that an accurate understanding of the turbulent flow is needed to describe the heat transfer to the environment.

Such convective boundary layers are not only generated by surface heating. For example, a vertical ice face submerged in salty water will drive convection due to the generation of fresh meltwater at the ice–water interface as it melts or dissolves (McConnochie & Kerr 2015; Malyarenko et al. 2020). In this case, the buoyancy driving the flow is primarily due to the salinity difference between the meltwater and the ambient water. One key difference between the two applications mentioned so far is the ratio of the diffusivities of momentum and heat (or salt), known as the Prandtl (or Schmidt) number ${Pr}$. In air the Prandtl number is ${Pr}\approx 0.7$, whereas for salt diffusion in cold water the relevant parameter is ${Pr}\approx 2000$.

Numerical simulations are often restricted to ${Pr}=O(1)$ because high spatial resolution is needed at high ${Pr}$ to resolve sharp scalar gradients that diffuse more slowly than the velocity gradients. However, understanding the role of the Prandtl number is vital for interpreting the results of such research for environmental or geophysical applications. We shall therefore investigate boundary layers in turbulent vertical convection at ${Pr}\gg 1$ with the aim of bridging the gap from classical studies of convection to geophysical applications.

In this study, we use direct numerical simulations to investigate turbulent convective boundary layers for a range of Rayleigh and Prandtl numbers. By using the multiple-resolution technique of Ostilla-Monico et al. (2015), we can efficiently simulate flows at high ${Pr}$, and we vary ${Pr}$ from 1 to 100. Although this is still considerably lower than the ${Pr}\approx 2000$ applicable to salt diffusion in the ocean, it is large enough to extract scaling laws in the large ${Pr}$ regime, which we expect to also hold in oceanographic flows.

Many different set-ups have been used to investigate VC boundary layers in numerical studies. Wang et al. (2021) recently simulated VC in a closed box, but the presence of walls in that domain means that turbulent boundary layers are only observed at very high $\mbox {{Ra}}$, at which only two-dimensional simulations are computationally feasible. We instead simulate the flow in a vertical channel with periodic boundary conditions in the wall-parallel directions. As originally described by Batchelor (1954), this domain approximates the flow at mid-heights in a tall vertical cell. A recent study by Ke et al. (2020) used this domain to simulate the temporally evolving boundary layer at a single heated wall, but in order to obtain converged statistics for a wide range of parameters, we instead consider the vertical channel set-up where one wall is heated and the other is cooled. This flow configuration achieves a statistically steady state with an anti-symmetric velocity profile and has been the subject of numerous numerical studies at ${Pr}=O(1)$ (e.g. Versteegh & Nieuwstadt 1999; Pallares et al. 2010; Ng et al. 2015).

The remainder of this paper is organised as follows. In § 2 we outline the numerical model and the set-up of the simulations. This is followed by flow visualisations in § 3 and a qualitative discussion of ${Pr}$-dependence of this flow. In § 4 we describe how various parameterisations for turbulent heat flux perform when applied to our simulations, and in § 5 we identify appropriate scaling laws for the boundary layer thicknesses. Finally, we conclude and discuss important remaining open questions for convective boundary layers in VC in § 6. The paper is supplemented by a concrete translation of our results into the geophysical context, focusing on the transition from laminar-type to turbulent-type boundary layers (Appendix A) and a detailed analysis of the energy dissipation and thermal dissipation budgets.

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. (2015) 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 (1996) and van der Poel et al. (2015). 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. (2015) 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 (1954) 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 1979) 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 1968; Tsuji & Nagano 1988) and is often provided in engineering reference texts such as Holman (2010).

However, recent analysis of numerical simulations by Ng et al. (2017) 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 (1999), 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. 2020) 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 (1972). However, recent analysis by Malyarenko et al. (2020) 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 1962; Spiegel 1971; Lohse & Toschi 2003; Ahlers, Grossmann & Lohse 2009). 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 2016). 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 (2000, 2001). This theory has shown excellent agreement with subsequent experimental and numerical investigations over a large range of $\mbox {{Ra}}$ and ${Pr}$ (Ahlers et al. 2009; Stevens et al. 2013). 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 (2008) applied ideas from the GL theory about boundary layer transition to geophysical-scale convection at a vertical wall, and Ng et al. (2015, 2017) 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. (2010), 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. (2021). 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 (2016) for VC and Grossmann & Lohse (2000) 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 (2016) 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 1954). Ng et al. (2017) 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 (2000), 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 (2016), 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. 2021).

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. (2017), 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 (1972), 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. (2002) from experiments of RBC with a range of large Prandtl numbers. Lam et al. (2002) 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 (2000, 2001), although Lam et al. (2002) 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 1979), and also assuming the Malkus (1954) 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 (2017), 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 (1999). 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. (2012). 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 (1999) 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 (2008) and in the spirit of Grossmann & Lohse (2000), 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. 2009; Ching et al. 2019) 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. (2015) 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 (1968), 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 (2016) 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 2016), applied to the regimes of Grossmann & Lohse (2000) 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 (1968) and Shishkina (2016). The ${Pr}^{-1/3}$ scaling is often also observed in empirical data for turbulent flows (e.g. Kader 1981). 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 (2008), similar to the marginally stable boundary layer argument of Malkus (1954) 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.

6. Conclusions

Through three-dimensional direct numerical simulations, we have investigated the multi-parameter dependence of convection in a vertical channel for Prandtl $1\leq {Pr} \leq 100$ and Rayleigh numbers $10^{6} \leq \mbox {{Ra}} \leq 10^{9}$. We observe Nusselt numbers consistent with the classical $\mbox {{Ra}}^{1/3}$ scaling combined with some weak but non-trivial dependence on the Prandtl number. The Reynolds number associated with the large-scale ‘wind’ exhibits a scaling of $\mbox {{Ra}}^{0.491} {Pr}^{-0.735}$, similar to that measured by Lam et al. (2002) in experiments of RBC. The discrepancy between the observed scaling and the theoretical prediction of ${Re}\sim \mbox {{Ra}}^{4/9}{Pr}^{-2/3}$ from Grossmann & Lohse (2000, 2001) for RBC, however, suggests there is more work to be done to build a theoretical understanding for the behaviour of turbulent VC. We cannot rule out the possibility that our observations arise due to a mixed scaling with contributions from multiple flow regimes.

As previously highlighted by McConnochie & Kerr (2017), such a scaling for $\mbox {{Nu}}$ is inconsistent with the commonly used heat flux parameterisation of Holland & Jenkins (1999). Our simulations highlight that this discrepancy is due to a highly variable drag coefficient in VC that depends on both of the control parameters $\mbox {{Ra}}$ and ${Pr}$. The absence of logarithmic velocity profiles suggests that the lack of shear-driven turbulent boundary layers is to blame for the large variation in the drag coefficient. By considering the critical Reynolds number of Landau & Lifshitz (1987) in Appendix A, we infer that transition to such turbulent boundary layers will only occur for $\mbox {{Ra}} > 4\times 10^{11} \times {Pr}^{1.89}$. However, more work is needed to understand how this transition occurs, and whether local scaling exponents for $\mbox {{Nu}}$ become impacted by multiple regimes and logarithmic corrections, as is the case for RBC (Grossmann & Lohse 2011) and convection from rough walls (MacDonald et al. 2019).

In contrast to the variation in the drag coefficient, the transfer coefficient (or modified Stanton number) satisfies $C_T\approx 0.1 {Pr}^{-2/3}$, matching values used in ice–ocean parameterisations. In other words, the friction velocity $V_\ast$ in this flow seems to adjust such that the heat flux scales as $q_T\sim V_\ast \varDelta$ for each given value of ${Pr}$. The strong dependence of $C_T$ on ${Pr}$ suggests that the conductive sublayer at the wall plays an important role in the total heat flux. We diagnose the width of this sublayer from the simulations and find the scaling $\delta _T/H \sim {{Re}}^{-2/3}{Pr}^{-1/3}$ to be consistent with our data. The emergent Rayleigh number $\mbox {{Ra}}_\delta$ associated with this sublayer is found to depend strongly on Prandtl number, questioning the notion of marginal stability at a critical value of $\mbox {{Ra}}_\delta$. This is similar to RBC, where the marginal stability theory of Malkus (1954) is also insufficient to fully describe the control parameter dependence of the heat flux (Ahlers et al. 2009).

Understanding how generic these results are will be vital for environmental applications. For example, Jackson et al. (2020) recently highlighted the role of mean horizontal flows in enhancing heat and salt transport at melting ice faces. In such a mixed convection scenario, ${Re}$ is not necessarily coupled to $\mbox {{Ra}}$ as it is in vertical convection. Thus understanding the underlying parameter dependence is an important topic for future research. As reviewed by Malyarenko et al. (2020), many factors not considered here can also be important for the ablation of ice in the ocean. In particular, the presence of both temperature and salinity variations and the dynamic melting condition may modify the nature of the boundary layers in this geophysical setting.