Boundary layers in turbulent vertical convection at high Prandtl number

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 $Ra$ and Prandtl number $Pr$. For the considered parameter range $1\leq Pr \leq 100$, $10^6 \leq 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.


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 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, 2 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 ≈ 0.7, whereas for salt diffusion in cold water the relevant parameter is Pr ≈ 2000. Numerical simulations are often restricted to Pr = (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 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 ≈ 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 setups have been used to investigate vertical convection boundary layers in numerical studies. Wang et al. (2021) recently simulated vertical convection in a closed box, but the presence of walls in that domain means that turbulent boundary layers are only observed at very high Ra, at which only 2-D simulations are computationally feasible. We instead simulate the flow in a vertical channel with periodic boundary conditions in the wallparallel 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 setup 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 = (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 setup 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.

Dynamical equations and control parameters
We consider the Navier-Stokes equations subject to the Oberbeck-Boussinesq approximation, where changes in density are only relevant in the buoyancy and a linear equation of state relates the density changes to temperature . These equations read ∇ · = 0 and where = ( , , ) is the velocity field, the kinematic pressure, kinematic viscosity, the molecular diffusivity of heat, gravitational acceleration, the thermal expansion coefficient, and 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 and the temperature difference between them is ∆. As in Ng et al. (2015) and shown in figure 1, we consider a domain of length 8 in the vertical ( ) and length 4 in the spanwise ( ) direction, and impose periodic boundary conditions on , , and in these directions, and . In a convective system, we can scale the velocity by the free-fall velocity = √︁ ∆ so that the dynamics of the system are solely determined by the Rayleigh and Prandtl numbers These are the only control parameters of the system, aside from parameters characterising the geometry of the flow domain. Their ratio Gr = Ra/Pr = 3 ∆/ 2 is also called the Grashof number.
The governing equations (2.1)-(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 timescale of the order of Pr −1 compared to the velocity field. We therefore also use the multiple-resolution technique of Ostilla-Monico et al. (2015) to evolve the scalar 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 ( ) direction using a clipped Chebyshev-type clustering. Uniform grid spacing is used in the and directions, and the base grid of all simulations are resolved down to Overview of the dimensionless parameters and grid resolutions used in the numerical simulations. Grid resolutions are listed here for the cases at highest Ra, and we distinguish between the base grid used to evolve the velocity and the refined grid used to evolve the temperature field. a factor of 2 times the Kolmogorov scale. The refined grid is such that the wall-normal grid spacing satisfies ∆ < 0.5 at the boundaries, and the grid spacing in the bulk satisfies ∆ , , < 4.5 , where = ( 2 / ) 1/4 is the Batchelor scale. The range of dimensionless control parameters simulated is shown in table 1. Simulations at 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 Ra are initialised using the final state of the simulation at Ra = 10 6 and Pr = 1, interpolated onto a new grid. Each computation is performed for at least 300 free-fall times, where / is the free-fall time unit. We average statistics over the last 250 time units once the system has reached a statistically steady state.

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.3) 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 is the horizontal heat flux, max is the peak value of the time-and spatially-averaged vertical velocity ( ), and * = √︁ / 0 is the friction velocity associated with the mean wall shear stress = / wall = 0 * 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 . Assuming a power-law relation between Nu and the Rayleigh number, dimensional analysis (e.g. Turner 1979) then requires the scaling Nu ∼ Ra 1/3 (Pr). This has been consistent with various experiments up to = 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 Ra, at which precise, controlled experiments and numerical Focus on Fluids articles must not exceed this page length 5 simulations are extremely difficult to perform. Finally, the Prandtl number dependence has hardly been addressed.
One important application for boundary layers in vertical convection 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 ∼ velocity × temperature change. Following Holland & Jenkins (1999), the parameterisation for the horizontal heat flux takes the form where is the vertical velocity of the rising plume, and − is the temperature difference between the ocean and the ice boundary. Taking = max and − = ∆/2, we note that the drag coefficient and 'transfer coefficient' from (2.5) are fully determined by the response parameters of (2.4) through (2.6) Accurate scaling laws for the quantities in (2.4) are thus crucial for determining and . The transfer coefficient is equivalent to a modified Stanton number where * is used for the velocity scale. In the ice-ocean literature, the transfer coefficient is often denoted Γ although we use here to avoid confusion with the aspect ratio Γ used throughout literature on convection. Both and 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 keeps the same value, but is reduced to reflect its dependence on the Schmidt number.
For , being constant, the dimensionless form of (2.5) is Nu ∼ RePr. Such a scaling is reminiscent of the 'ultimate' or 'diffusion-free' scaling hypothesised for Rayleigh-Bénard convection (RBC) at very high Ra (e.g. Kraichnan 1962;Spiegel 1971;Lohse & Toschi 2003;Ahlers et al. 2009). In that case, the heat flux is assumed independent of the molecular diffusivities and , such that dimensional analysis implies Nu ∼ ( ) 1/2 . In physical terms, this regime is associated with a dominant large-scale circulation that leads to sheardriven 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 to 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 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 on how the Prandtl number affects (a) the scaling of the above response parameters in the currently accessible range of Ra, and (b) any subsequent transition in the nature of the boundary layers.

Flow visualisation
To illustrate how the boundary layers in the flow change with Pr and Ra, we present a snapshot of the local dimensionless vertical shear stress and heat flux at the heated wall = 0 in figure 2. These quantities are defined as We note that averaging over the plane and over time gives the Nusselt number ( , , ) , , = Nu. In this sense, 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 ) are frequently co-located with instantaneous flow reversals (evidenced by white and blue patches appearing in the panels for ).
Panels ( )-( ) of figure 2 highlight the effect of increasing Pr in this setup while keeping Ra fixed (in this case at 10 8 ). As seen from the colour scales, although the range of the local Nusselt number 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 Gr as Pr is increased for fixed 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 field for the high Gr snapshot of ( ) compared to the lower Gr snapshot of ( ). By contrast, comparing panels ( ) and ( ) allows us to visualise the effect of changing Pr while keeping Gr fixed. Qualitatively the structures in both the and 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 Gr at the walls, although this similarity does not extend outside of the viscous boundary layer.
Compared to Rayleigh-Bénard convection (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 ( ) direction are prominent in figure 2, similar to those seen in the sheared RBC setup of Blass et al. (2021). Furthermore, in that study, an increased Pr (for fixed 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 vertical convection setup. Indeed, the simple analogy between the Grashof number and the square of the Reynolds number should not be overstated. As shown in figure 1( ), the peak value of the time-averaged vertical velocity does not simply scale with the free-fall velocity , 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.

Heat flux and Reynolds number parameterisation
In table 2 we report the observed Raand Pr-dependence of the response parameters from (2.4) and (2.6) in our simulations. An effective power-law dependence is assumed and twoparameter linear regression is used to obtain the effective scaling exponents. Precisely, we compute b = −1 y, where b = ( 1 , 2 , 3 ) and are constructed from the simulations for each response parameter , giving a linear fit = Ra 1 Pr 2 10 3 . We calculate the uncertainty of the power law exponents 1 and 2 through the variance matrix of b given by = 2 ( ) −1 , where 2 is the variance of y − b. The standard deviations of the slopes, given by √ 11 and √ 22 are presented in table 2.
The Nusselt number is consistent with the theoretical scaling relation Nu ∼ Ra 1/3 (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 ≈ 1, a regime transition to a shear-dominated boundary layer is underway at Ra = 10 9 , but following Grossmann & Lohse (2000), this transitional Ra can be expected to increase with Pr, as the smaller Reynolds number stabilizes the flow. Our results contrast with the effective scaling laws for laminar vertical convection derived by Shishkina (2016) difference is to be expected since our setup is far from the laminar state for which the scaling laws have been observed to hold (e.g. by Wang et al. 2021).
In figure 3 we plot Nu against both Ra and the shear Reynolds number Re . Figure 3a highlights the weak dependence of Nu on Pr, with higher Pr typically reducing Nu for a fixed value of 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 Ra 1/3 line on figure 3a at higher values of 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.
Against Re in figure 3b, we obtain a reasonable collapse for Nu by scaling with Pr 1/3 and observe a scaling close to Nu ∼ Re 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. For example, a recent study by Yerragolam et al. (2021) proposed a scaling theory for passive scalar transport in plane Couette flow where Nu ∼ Re 6/7 Pr 1/2 . This somewhat contrasts with the Pr 1/3 collapse observed in figure 3b, although the higher Re 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 Rayleigh-Bénard convection with a range of large Prandtl numbers. Lam et al. (2002) suggested that their results were consistent with the theoretical scaling relation Re ∼ Ra 4/9 Pr −2/3 proposed for the regime ( ) associated with Nu ∼ 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 ∼ Ra 4/9 Pr −2/3 scaling can also be derived from dimensional analysis by assuming that the vertical velocity max is solely determined by the buoyancy flux per unit area Φ = and the plate separation (as in the 'outer' scaling of George & Capp 1979), and also assuming the Malkus (1954) scaling Nu ∼ 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, like 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 Nu ∼ Ra 1/3 (Pr) implies a dimensional form for the heat flux that scales as ∼ ∆ 4/3 for fixed fluid properties. The heat flux is therefore independent of the bulk velocity 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 with both Ra and Pr. In all cases we find a value much larger than the high-Re limit of = 2.5 × 10 −3 , as used by Holland & Jenkins (1999). However, the scaling observed for the transfer coefficient ≈ 0.1Pr −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.6), we can express this result in terms of the Nusselt number as Nu ∼ Re Pr 1/3 or with dimensional quantities as ∼ Pr −2/3 * ∆. It may be tempting to associate the scaling Nu ∼ Re 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 Ra = 10 8 , 10 9 with a logarithmic -axis. From figure 5a, 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 sheardriven 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( ) 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.
Holland & Jenkins (1999) associate the scaling relation ∼ −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 vertical convection system.

Conductive thermal boundary layer
In a statistically steady state, the mean velocity and temperature profiles of the system satisfy where an overbar denotes an average in , , and time. Incompressibility ensures that ≡ 0, so the mean velocity only has components in the wall-parallel directions. The second equation of (5.1) implies that the heat flux at any wall-normal location must be constant, or in dimensionless terms 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 as the wall-normal location where the conductive heat flux − / is equal to the turbulent heat flux . In Rayleigh-Bénard convection at moderate 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 ∼ Re −1/2 (Pr). From our new simulations, we find a collapse of the data such that / ∼ Pr −1/3 (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 is more consistent with ∼ Re −2/3 Pr −1/3 , (5.4) as shown in figure 6( ).
For the case where max ∼ , the scaling of (5.4) is equivalent to / ∼ 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 Rayleigh-Bénard convection. However, as we already mentioned earlier, max does not simply scale with in our simulations. In figure 6( ), we plot the 'sublayer Rayleigh number' Ra = 3 ∆/ , 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 setup, would allow Re and Ra to be decoupled, and reveal the inherent parameter dependence of the boundary layer.

Conclusions
Through three-dimensional direct numerical simulations, we have investigated the multiparameter dependence of convection in a vertical channel for Prandtl 1 Pr 100 and Rayleigh numbers 10 6 Ra 10 9 . We observe Nusselt numbers consistent with the classical 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 Ra 0.491 Pr −0.735 , similar to that measured by Lam et al. (2002) in experiments of Rayleigh-Bénard convection. The discrepancy between the observed scaling and the theoretical prediction of Re ∼ 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 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 vertical convection that depends on both of the control parameters 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 Ra > 4 × 10 11 × Pr 1.89 . However, more work is needed to understand how this transition occurs, and whether local scaling exponents for 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 ≈ 0.1Pr −2/3 , matching values used in ice-ocean parameterisations. In other words, the friction velocity * in this flow seems to adjust such that the heat flux scales as ∼ * ∆ for each given value of Pr. The strong dependence of 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 / ∼ Re −2/3 Pr −1/3 to be consistent with our data. The emergent Rayleigh number Ra associated with this sublayer is found to depend strongly on Prandtl number, questioning the notion of marginal stability at a critical value of Ra . 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 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.
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.

Appendix A. Boundary layer transition prediction
In this appendix, we provide an estimate for the Pr-dependence of the transition to a sheardriven turbulent boundary layer, based on the critical Reynolds number criterion of Landau & Lifshitz (1987). From each simulation, we can calculate a Reynolds number Re * based on the displacement thickness * by * = ∫ where max is the maximum vertical velocity and 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 Re * = 0.159Ra 0.294 Pr −0.557 . (A 2) Assuming (somewhat ambitiously) that this scaling remains valid up to a critical Reynolds number of Re * = Re = 420, we infer a Pr-dependent critical Rayleigh number of = 4.27 × 10 11 × Pr 1.89 .
(A 3) For Pr = 1, this gives a value within the transition range of 3.8×10 10 10 12 predicted by Ng et al. (2017) in figure 10 of that paper.
In the context of a melting vertical ice face in the ocean, we can use (A 3) 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 saltfree, at a water temperature of 2°C the interfacial concentration of salinity is approximately 15 g kg −1 (see e.g. Kerr & McConnochie 2015). Combined with an ambient ocean salinity of 35 g kg −1 , a haline contraction coefficient of = 7.86 × 10 −4 (g kg −1 ) −1 , a kinematic viscosity of = 1 × 10 −6 m 2 s −1 , and a Schmidt number = / = 2600, we find = 3 ∆ ≈ 10 18 ≈ 4 × 10 14 3 , implying that ≈ 13.5 m. (A 4) Note that 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 (2008), the plume width can be linearly related to the height from the base of the ice by ≈ 0.1 . This relation is based on the constant entrainment rate assumption of classical plume theory as developed by Morton et al. (1956). The critical vertical position for a shear-driven boundary layer is then = 135 m, associated with a Rayleigh number of = 10 21 . This matches the prediction of Kerr & McConnochie (2015) 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 2016) or the pressure-dependence of the melt condition at the boundary of the ice (Hewitt 2020). 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 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 , the turbulent kinetic energy , and the equivalent quantities for the temperature field where the shear production P , which transfers energy between turbulence and the mean flow, is defined as Here P is an analogous term to the shear production described above, and quantifies the interaction between the mean temperature profile and the turbulent fluctuations: The thermal dissipation rates are given by In the statistically steady states reached by our simulations, the energies become constant in time, such that we get the following relations from (B 2) and (B 6) = P + , = P + , = P + , P = .
(B 10) These equations highlight how the total vertical heat flux = + can be related to the kinetic energy dissipation rate, and how the horizontal heat flux can be related to the thermal dissipation rate: = + = + , = + . (B 11) Figure 7 plots the relative contributions of each of these budget terms to the heat fluxes as a function of Ra and Pr. For the kinetic energy budget terms (shown in panels and ), we observe that the relative contributions of and increase with Ra and decrease with Pr. This also coincides with an increase in the relative magnitude of the shear production P . Since P 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( , ) 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( ) show very weak dependence on 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 Ra. As Pr increases the relative contribution of becomes greater. This highlights once again the key role that the thin, conductive boundary layers, whose strong gradients contribute to , have on the heat flux in vertical convection at high Pr.