Hostname: page-component-7c8c6479df-hgkh8 Total loading time: 0 Render date: 2024-03-29T06:27:46.834Z Has data issue: false hasContentIssue false

Boundary layers in turbulent vertical convection at high Prandtl number

Published online by Cambridge University Press:  16 November 2021

Christopher J. Howland*
Affiliation:
Physics of Fluids Group, Max Planck Center for Complex Fluid Dynamics, MESA+ Institute and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, Netherlands
Chong Shen Ng
Affiliation:
Physics of Fluids Group, Max Planck Center for Complex Fluid Dynamics, MESA+ Institute and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, Netherlands
Roberto Verzicco
Affiliation:
Physics of Fluids Group, Max Planck Center for Complex Fluid Dynamics, MESA+ Institute and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, Netherlands Dipartimento di Ingegneria Industriale, University of Rome ‘Tor Vergata’, Via del Politecnico 1, Roma 00133, Italy Gran Sasso Science Institute, Viale F. Crispi, 7, 67100 L'Aquila, Italy
Detlef Lohse
Affiliation:
Physics of Fluids Group, Max Planck Center for Complex Fluid Dynamics, MESA+ Institute and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, Netherlands Max Planck Institute for Dynamics and Self-Organization, Am Fassberg 17, 37077 Göttingen, Germany
*
Email address for correspondence: c.j.howland@outlook.com

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.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© The Author(s), 2021. Published by Cambridge University Press

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 Reference Batchelor1954; Kuiken Reference Kuiken1968; Shishkina Reference Shishkina2016) 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 Reference McConnochie and Kerr2015; Malyarenko et al. Reference Malyarenko, Wells, Langhorne, Robinson, Williams and Nicholls2020). 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. (Reference Ostilla-Monico, Yang, van der Poel, Lohse and Verzicco2015), 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. (Reference Wang, Liu, Verzicco, Shishkina and Lohse2021) 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 (Reference Batchelor1954), this domain approximates the flow at mid-heights in a tall vertical cell. A recent study by Ke et al. (Reference Ke, Williamson, Armfield, Norris and Komiya2020) 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 Reference Versteegh and Nieuwstadt1999; Pallares et al. Reference Pallares, Vernet, Ferre and Grau2010; Ng et al. Reference Ng, Ooi, Lohse and Chung2015).

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. (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.4ac)\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.4ac) 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.4ac) 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(ac) 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.4ac) 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.1ad)\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.

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. (Reference Lam, Shang, Zhou and Xia2002) 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 (Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001) 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 (Reference McConnochie and Kerr2017), such a scaling for $\mbox {{Nu}}$ is inconsistent with the commonly used heat flux parameterisation of Holland & Jenkins (Reference Holland and Jenkins1999). 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 (Reference Landau and Lifshitz1987) 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 Reference Grossmann and Lohse2011) and convection from rough walls (MacDonald et al. Reference MacDonald, Hutchins, Lohse and Chung2019).

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 (Reference Malkus1954) is also insufficient to fully describe the control parameter dependence of the heat flux (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009).

Understanding how generic these results are will be vital for environmental applications. For example, Jackson et al. (Reference Jackson, Nash, Kienholz, Sutherland, Amundson, Motyka, Winters, Skyllingstad and Pettit2020) 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. (Reference Malyarenko, Wells, Langhorne, Robinson, Williams and Nicholls2020), 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.

Acknowledgements

We thank three anonymous referees for their helpful and insightful comments.

Funding

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 804283). We acknowledge PRACE for awarding us access to Joliot-Curie at GENCI@CEA, France, and this work was also sponsored by NWO Science for the use of supercomputer facilities.

Declaration of interests

The authors report no conflict of interest.

Data availability statement

The data that support the findings of this study are openly available in 4TU.ResearchData at http://doi.org/10.4121/16836874. The software used to perform the simulations described in this study can be freely obtained at https://github.com/chowland/AFiD-MuRPhFi.

Appendix A. Boundary layer transition prediction

In this appendix, we provide an estimate for the ${Pr}$-dependence of the transition to a shear-driven turbulent boundary layer, based on the critical Reynolds number criterion of Landau & Lifshitz (Reference Landau and Lifshitz1987). From each simulation, we can calculate a Reynolds number ${Re}_{\delta ^{\ast }}$ based on the displacement thickness $\delta ^{\ast }$ by

(A1a,b)\begin{equation} \delta^{{\ast}} = \int_0^{x_{max}} 1 - \frac{\bar{v}(x)}{V_{max}}\,\mathrm{d}\kern0.7pt x,\quad Re_{\delta^{{\ast}}} = \frac{V_{max}\delta^{{\ast}}}{\nu}, \end{equation}

where $V_{max}$ is the maximum vertical velocity and $x_{max}$ is the wall-normal location of this maximum. Performing the same linear regression as described in § 4 on this data, we obtain the power-law relation

(A2)\begin{equation} {Re}_{\delta^{{\ast}}} = 0.159 \mbox{{Ra}}^{0.294} {Pr}^{{-}0.557} . \end{equation}

Assuming (somewhat ambitiously) that this scaling remains valid up to a critical Reynolds number of ${Re}_{\delta ^{\ast }}={Re}_c = 420$, we infer a ${Pr}$-dependent critical Rayleigh number of

(A3)\begin{equation} Ra_c = 4.27\times 10^{11} \times {Pr}^{1.89} . \end{equation}

For ${Pr}=1$, this gives a value within the transition range of $3.8\times 10^{10}\lesssim Ra_c \lesssim 10^{12}$ predicted by Ng et al. (Reference Ng, Ooi, Lohse and Chung2017) in figure 10 of that paper.

In the context of a melting vertical ice face in the ocean, we can use (A3) to estimate the length scales at which a shear-driven boundary layer may be relevant in describing the salt flux towards the ice due to natural convection. Although the ice can be considered salt free, at a water temperature of $2\,^{\circ }\textrm {C}$ the interfacial concentration of salinity is approximately $15\,\textrm {g}\,\textrm {kg}^{-1}$ (see e.g. Kerr & McConnochie Reference Kerr and McConnochie2015). Combined with an ambient ocean salinity of $35\,\textrm {g}\,\textrm {kg}^{-1}$, a haline contraction coefficient of $\beta = 7.86\times 10^{-4}\,{(\textrm {g}\,\textrm {kg}^{-1})^{-1}}$, a kinematic viscosity of $\nu =1 \times 10^{-6}\,\textrm {m}^{2}\,\textrm {s}^{-1}$ and a Schmidt number $Sc=\nu /\kappa _S = 2600$, we find

(A4)\begin{equation} Ra_c = \frac{g\beta H_c^{3} \varDelta_S}{\nu\kappa_S} \approx 10^{18} \approx 4\times 10^{14} H_c^{3},\quad \textrm{implying that}\ H_c \approx 13.5\,\textrm{m}. \end{equation}

Note that $H_c$ is the critical horizontal length scale. In the context of convection at an ice face, where the domain is essentially unbounded in one direction, this is best compared with the local plume width. Following Wells & Worster (Reference Wells and Worster2008), the plume width $H$ can be linearly related to the height $Z$ from the base of the ice by $H\approx 0.1 Z$. This relation is based on the constant entrainment rate assumption of classical plume theory as developed by Morton, Taylor & Turner (Reference Morton, Taylor and Turner1956). The critical vertical position for a shear-driven boundary layer is then $Z_c=135\,\textrm {m}$, associated with a Rayleigh number of $Ra_z=10^{21}$. This matches the prediction of Kerr & McConnochie (Reference Kerr and McConnochie2015) who used GL theory to estimate the transition. Over such large vertical distances, other physical phenomena are likely to play an important role in the dynamics, such as ambient stratification (McConnochie & Kerr Reference McConnochie and Kerr2016) or the pressure dependence of the melt condition at the boundary of the ice (Hewitt Reference Hewitt2020). It is therefore unlikely that a shear-driven boundary layer would develop at an ice face solely due to natural convection, without some external forcing such as subglacial discharge or a mean horizontal current.

Appendix B. Turbulence budgets

Finally, to gain more insight into the nature of the flow as $\mbox {{Ra}}$ and ${Pr}$ vary, we present results from the turbulence budgets of our simulations and describe how the turbulent kinetic and thermal dissipation rates are related to the heat flux in the system. From the governing equations, we can construct evolution equations for the kinetic energy of the mean flow $\overline {E_K}$, the turbulent kinetic energy $E_K'$, and the equivalent quantities for the temperature field

(B1ad)\begin{equation} \overline{E_K} = \tfrac{1}{2}\left\langle\left|\bar{\boldsymbol{u}}\right|^{2} \right\rangle, \quad {E_K}^{\prime} = \tfrac{1}{2} \left\langle\left|\boldsymbol{u}^{\prime}\right|^{2}\right\rangle, \quad \overline{E_T} = \tfrac{1}{2}\left\langle |\bar{T}|^{2} \right\rangle, \quad {E_T}^{\prime} = \tfrac{1}{2}\left\langle |T^{\prime}|^{2} \right\rangle, \end{equation}

whereas in the main text an overbar denotes an average over $y$ and $z$, and angle brackets denote a domain average. The evolution equations for the kinetic energies read

(B2a,b)\begin{equation} \frac{{\textrm{d}}\overline{E_K}}{{\textrm{d}}t} ={-}\mathcal{P}_S - \bar{\varepsilon} + \bar{q}, \quad \frac{{\textrm{d}}{E_K}^{\prime}}{{\textrm{d}}t} = \mathcal{P}_S - \varepsilon' + q^{\prime} , \end{equation}

where the shear production $\mathcal {P}_S$, which transfers energy between turbulence and the mean flow, is defined as

(B3)\begin{equation} \mathcal{P}_S ={-}\overline{u'\boldsymbol{u}'\boldsymbol{\cdot} \frac{\partial \bar{\boldsymbol{u}}}{\partial x}}, \end{equation}

the dissipation rates of mean KE and TKE are

(B4a,b)\begin{equation} \bar{\varepsilon} = \nu \left\langle \left| \frac{\partial \bar{\boldsymbol{u}}}{\partial x}\right|^{2} \right\rangle , \quad \varepsilon' = \nu \left\langle \frac{\partial u_i}{\partial x_j} \frac{\partial u_i}{\partial x_j} \right\rangle, \end{equation}

and the vertical heat fluxes due to the mean and turbulent profiles are given by

(B5a,b)\begin{equation} \bar{q} = g\alpha \left\langle \bar{v}\bar{T}\right\rangle,\quad q' = g\alpha \left\langle v'T'\right\rangle . \end{equation}

The mean square temperature and the temperature variance evolve according to similar equations, namely

(B6a,b)\begin{equation} \frac{{\textrm{d}}\overline{E_T}}{{\textrm{d}}t} ={-}\mathcal{P}_T - \bar{\chi} + q_T, \quad \frac{{\textrm{d}}{E_T}^{\prime}}{{\textrm{d}}t} = \mathcal{P}_T - \chi' . \end{equation}

Here, $\mathcal {P}_T$ is an analogous term to the shear production described above, and quantifies the interaction between the mean temperature profile and the turbulent fluctuations

(B7)\begin{equation} \mathcal{P}_T ={-}\left\langle u'T' \frac{\partial \bar{T}}{\partial x} \right\rangle. \end{equation}

The thermal dissipation rates are given by

(B8a,b)\begin{equation} \overline\chi = \kappa \left\langle \left(\frac{\partial \bar{T}}{\partial x} \right)^{2} \right\rangle ,\quad \chi' = \kappa \left\langle \left|\boldsymbol{\nabla} T\right|^{2}\right\rangle , \end{equation}

and $q_b$ is the mean horizontal heat flux through the boundaries

(B9)\begin{equation} q_T = \frac{\kappa}{2}\left(\left.\frac{\partial \bar{T}}{\partial x}\right|_{x=0} + \left.\frac{\partial \bar{T}}{\partial x}\right|_{x=H}\right) = \frac{\mbox{{Nu}} H}{\kappa {\rm \Delta} T}. \end{equation}

In the statistically steady states reached by our simulations, the energies become constant in time, such that we get the following relations: from (B2a,b) and (B6a,b)

(B10ad)\begin{equation} \bar{q} = \mathcal{P}_S + \bar{\varepsilon}, \quad \varepsilon' = \mathcal{P}_S + q', \quad q_T = \mathcal{P}_T + \bar{\chi}, \quad \mathcal{P}_T = \chi'. \end{equation}

These equations highlight how the total vertical heat flux $q_v = \bar {q} + q^{\prime }$ can be related to the kinetic energy dissipation rate, and how the horizontal heat flux $q_T$ can be related to the thermal dissipation rate

(B11a,b)\begin{equation} q_v = \bar{q} + q^{\prime} = \bar{\varepsilon} + \varepsilon^{\prime}, \quad q_T = \bar{\chi} + \chi^{\prime} . \end{equation}

Figure 7 plots the relative contributions of each of these budget terms to the heat fluxes as a function of $\mbox {{Ra}}$ and ${Pr}$. For the kinetic energy budget terms (shown in (a,b)), we observe that the relative contributions of $\varepsilon '$ and $\bar {q}$ increase with $\mbox {{Ra}}$ and decrease with ${Pr}$. This also coincides with an increase in the relative magnitude of the shear production $\mathcal {P}_S$. Since $\mathcal {P}_S$ is positive in all our simulations, this means that energy is always (on average) transferred from the mean flow to the turbulent perturbations. The trends observed in figure 7(a,b) suggest that the kinetic energy budget terms may be most sensitive to the Reynolds number of the flow. By contrast, the relative contributions of the thermal dissipation rates plotted in figure 7(c) show very weak dependence on $\mbox {{Ra}}$. For ${Pr}$ fixed at 10, the dissipation of the mean temperature accounts for 60 % of the horizontal heat flux, and this fraction changes by less than 3 % over three decades of $\mbox {{Ra}}$. As ${Pr}$ increases the relative contribution of $\bar {\chi }$ becomes greater. This highlights once again the key role that the thin, conductive boundary layers, whose strong gradients contribute to $\bar {\chi }$, have on the heat flux in VC at high ${Pr}$.

Figure 7. Relative contributions to the heat flux due to the energy budget terms. (a,b) Plot of the kinetic energy budget terms as a fraction of the total vertical heat flux; (c,d) plot of the thermal dissipation rates as a fraction of the total horizontal heat flux. (a,c) Show variation with Rayleigh number for simulations at fixed ${Pr} = 10$, (b,d) show variation with Prandtl number for simulations at fixed $\mbox {{Ra}} = 10^{8}$.

References

REFERENCES

Ahlers, G., Bodenschatz, E., Funfschilling, D., Grossmann, S., He, X., Lohse, D., Stevens, R.J.A.M. & Verzicco, R. 2012 Logarithmic temperature profiles in turbulent Rayleigh–Bénard convection. Phys. Rev. Lett. 109 (11), 114501.CrossRefGoogle ScholarPubMed
Ahlers, G., Grossmann, S. & Lohse, D. 2009 Heat transfer and large scale dynamics in turbulent Rayleigh–Bénard convection. Rev. Mod. Phys. 81 (2), 503537.CrossRefGoogle Scholar
Batchelor, G.K. 1954 Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures. Q. Appl. Maths 12 (3), 209233.CrossRefGoogle Scholar
Blass, A., Tabak, P., Verzicco, R., Stevens, R.J.A.M. & Lohse, D. 2021 The effect of Prandtl number on turbulent sheared thermal convection. J. Fluid Mech. 910, A37.CrossRefGoogle Scholar
Ching, E.S.C., Leung, H.S., Zwirner, L. & Shishkina, O. 2019 Velocity and thermal boundary layer equations for turbulent Rayleigh–Bénard convection. Phys. Rev. Res. 1 (3), 033037.CrossRefGoogle Scholar
George, W.K. & Capp, S.P. 1979 A theory for natural convection turbulent boundary layers next to heated vertical surfaces. Intl J. Heat Mass Transfer 22 (6), 813826.CrossRefGoogle Scholar
Grossmann, S. & Lohse, D. 2000 Scaling in thermal convection: A unifying theory. J. Fluid Mech. 407, 2756.CrossRefGoogle Scholar
Grossmann, S. & Lohse, D. 2001 Thermal convection for large Prandtl numbers. Phys. Rev. Lett. 86 (15), 33163319.CrossRefGoogle ScholarPubMed
Grossmann, S. & Lohse, D. 2011 Multiple scaling in the ultimate regime of thermal convection. Phys. Fluids 23 (4), 045108.CrossRefGoogle Scholar
Hewitt, I.J. 2020 Subglacial plumes. Annu. Rev. Fluid Mech. 52 (1), 145169.CrossRefGoogle Scholar
Holland, D.M. & Jenkins, A. 1999 Modeling thermodynamic ice–ocean interactions at the base of an ice shelf. J. Phys. Oceanogr. 29 (8), 17871800.2.0.CO;2>CrossRefGoogle Scholar
Holman, J.P. 2010 Heat Transfer, 10th edn. Boston: McGraw Hill.Google Scholar
Jackson, R.H., Nash, J.D., Kienholz, C., Sutherland, D.A., Amundson, J.M., Motyka, R.J., Winters, D., Skyllingstad, E. & Pettit, E.C. 2020 Meltwater intrusions reveal mechanisms for rapid submarine melt at a tidewater glacier. Geophys. Res. Lett. 47 (2), e2019GL085335.CrossRefGoogle Scholar
Kader, B. 1981 Temperature and concentration profiles in fully turbulent boundary layers. Intl J. Heat Mass Transfer 24 (9), 15411544.CrossRefGoogle Scholar
Kader, B.A. & Yaglom, A.M. 1972 Heat and mass transfer laws for fully turbulent wall flows. Intl J. Heat Mass Transfer 15 (12), 23292351.CrossRefGoogle Scholar
Ke, J., Williamson, N., Armfield, S.W., Norris, S.E. & Komiya, A. 2020 Law of the wall for a temporally evolving vertical natural convection boundary layer. J. Fluid Mech. 902, A31.CrossRefGoogle Scholar
Kerr, R.C. & McConnochie, C.D. 2015 Dissolution of a vertical solid surface by turbulent compositional convection. J. Fluid Mech. 765, 211228.CrossRefGoogle Scholar
Kraichnan, R.H. 1962 Turbulent thermal convection at arbitrary Prandtl number. Phys. Fluids 5 (11), 13741389.CrossRefGoogle Scholar
Kuiken, H.K. 1968 An asymptotic solution for large Prandtl number free convection. J. Engng Maths 2 (4), 355371.CrossRefGoogle Scholar
Lam, S., Shang, X.-D., Zhou, S.-Q. & Xia, K.-Q. 2002 Prandtl number dependence of the viscous boundary layer and the Reynolds numbers in Rayleigh–Bénard convection. Phys. Rev. E 65 (6), 066306.CrossRefGoogle ScholarPubMed
Landau, L.D. & Lifshitz, E.M. 1987 Fluid Mechanics, 2nd edn, Course of Theoretical Physics, vol. 6. Pergamon Press.Google Scholar
Lohse, D. & Toschi, F. 2003 Ultimate state of thermal convection. Phys. Rev. Lett. 90 (3), 034502.CrossRefGoogle ScholarPubMed
MacDonald, M., Hutchins, N., Lohse, D. & Chung, D. 2019 Heat transfer in rough-wall turbulent thermal convection in the ultimate regime. Phys. Rev. Fluids 4 (7), 071501.CrossRefGoogle Scholar
Malkus, W.V.R. 1954 The heat transport and spectrum of thermal turbulence. Proc. R. Soc. Lond. A 225 (1161), 196212.Google Scholar
Malyarenko, A., Wells, A.J., Langhorne, P.J., Robinson, N.J., Williams, M.J.M. & Nicholls, K.W. 2020 A synthesis of thermodynamic ablation at ice–ocean interfaces from theory, observations and models. Ocean Model. 154, 101692.CrossRefGoogle Scholar
McConnochie, C.D. & Kerr, R.C. 2015 The turbulent wall plume from a vertically distributed source of buoyancy. J. Fluid Mech. 787, 237253.CrossRefGoogle Scholar
McConnochie, C.D. & Kerr, R.C. 2016 The effect of a salinity gradient on the dissolution of a vertical ice face. J. Fluid Mech. 791, 589607.CrossRefGoogle Scholar
McConnochie, C.D. & Kerr, R.C. 2017 Testing a common ice-ocean parameterization with laboratory experiments. J. Geophys. Res. 122 (7), 59055915.CrossRefGoogle Scholar
Morton, B.R., Taylor, G. & Turner, J.S. 1956 Turbulent gravitational convection from maintained and instantaneous sources. Proc. R. Soc. A 234 (1196), 123.Google Scholar
Ng, C.S., Ooi, A., Lohse, D. & Chung, D. 2015 Vertical natural convection: application of the unifying theory of thermal convection. J. Fluid Mech. 764, 349361.CrossRefGoogle Scholar
Ng, C.S., Ooi, A., Lohse, D. & Chung, D. 2017 Changes in the boundary-layer structure at the edge of the ultimate regime in vertical natural convection. J. Fluid Mech. 825, 550572.CrossRefGoogle Scholar
Ostilla-Monico, R., Yang, Y., van der Poel, E.P., Lohse, D. & Verzicco, R. 2015 A multiple-resolution strategy for direct numerical simulation of scalar turbulence. J. Comput. Phys. 301, 308321.CrossRefGoogle Scholar
Pallares, J., Vernet, A., Ferre, J.A. & Grau, F.X. 2010 Turbulent large-scale structures in natural convection vertical channel flow. Intl J. Heat Mass Transfer 53 (19), 41684175.CrossRefGoogle Scholar
van der Poel, E.P., Ostilla-Mónico, R., Donners, J. & Verzicco, R. 2015 A pencil distributed finite difference code for strongly turbulent wall-bounded flows. Comput. Fluids 116, 1016.CrossRefGoogle Scholar
Schlichting, H. & Gersten, K. 2016 Boundary-Layer Theory, 9th edn. Springer, Berlin, Heidelberg.Google Scholar
Shishkina, O. 2016 Momentum and heat transport scalings in laminar vertical convection. Phys. Rev. E 93 (5), 051102.CrossRefGoogle ScholarPubMed
Shishkina, O. & Horn, S. 2016 Thermal convection in inclined cylindrical containers. J. Fluid Mech. 790, R3.CrossRefGoogle Scholar
Spiegel, E.A. 1971 Convection in stars. I. Basic Boussinesq convection. Annu. Rev. Astron. Astrophys. 9, 323352.CrossRefGoogle Scholar
Stevens, R.J.A.M., van der Poel, E.P., Grossmann, S. & Lohse, D. 2013 The unifying theory of scaling in thermal convection: the updated prefactors. J. Fluid Mech. 730, 295308.CrossRefGoogle Scholar
Tsuji, T. & Nagano, Y. 1988 Characteristics of a turbulent natural convection boundary layer along a vertical flat plate. Intl J. Heat Mass Transfer 31 (8), 17231734.CrossRefGoogle Scholar
Turner, J.S. 1979 Buoyancy Effects in Fluids, 1st edn. Cambridge: Cambridge University Press.Google Scholar
Versteegh, T.A.M. & Nieuwstadt, F.T.M. 1999 A direct numerical simulation of natural convection between two infinite vertical differentially heated walls scaling laws and wall functions. Intl J. Heat Mass Transfer 42 (19), 36733693.CrossRefGoogle Scholar
Verzicco, R. & Orlandi, P. 1996 A finite-difference scheme for three-dimensional incompressible flows in cylindrical coordinates. J. Comput. Phys. 123 (2), 402414.CrossRefGoogle Scholar
Wang, Q., Liu, H.-R., Verzicco, R., Shishkina, O. & Lohse, D. 2021 Regime transitions in thermally driven high-Rayleigh number vertical convection. J. Fluid Mech. 917, A6.CrossRefGoogle Scholar
Warner, C.Y. & Arpaci, V.S. 1968 An experimental investigation of turbulent natural convection in air at low pressure along a vertical heated flat plate. Intl J. Heat Mass Transfer 11 (3), 397406.CrossRefGoogle Scholar
Wells, A.J. & Worster, M.G. 2008 A geophysical-scale model of vertical natural convection boundary layers. J. Fluid Mech. 609, 111137.CrossRefGoogle Scholar
Figure 0

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)$.

Figure 1

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.

Figure 2

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 3

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.

Figure 4

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}$).

Figure 5

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.

Figure 6

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$.

Figure 7

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.

Figure 8

Figure 7. Relative contributions to the heat flux due to the energy budget terms. (a,b) Plot of the kinetic energy budget terms as a fraction of the total vertical heat flux; (c,d) plot of the thermal dissipation rates as a fraction of the total horizontal heat flux. (a,c) Show variation with Rayleigh number for simulations at fixed ${Pr} = 10$, (b,d) show variation with Prandtl number for simulations at fixed $\mbox {{Ra}} = 10^{8}$.