Mean velocity and temperature profiles in turbulent Rayleigh–Bénard convection at low Prandtl numbers

Abstract We report a direct numerical simulation (DNS) study of the mean velocity and temperature profiles in turbulent Rayleigh–Bénard convection (RBC) at low Prandtl numbers ($Pr$). The numerical study is conducted in a vertical thin disk with $Pr$ varied in the range $0.17\leq Pr\leq 4.4$ and the Rayleigh number ($Ra$) varied in the range $5\times 10^8\leq Ra \leq 1\times 10^{10}$. By varying $Pr$ from 4.4 to 0.17, we find a sharp change of flow patterns for the large-scale circulation (LSC) from a rigid-body rotation to a near-wall turbulent jet. We numerically examine the mean velocity equation in the bulk region and find that the mean horizontal velocity profile $u(z)$ can be determined by a balance equation between the mean convection and turbulent diffusion with a constant turbulent viscosity $\nu _t$. This balance equation admits a self-similarity jet solution, which fits the DNS data well. In the boundary-layer region, we find that both the mean temperature profile $T(z)$ and $u(z)$ can be determined by a balance equation between the molecular diffusion and turbulent diffusion. Within the viscous boundary layer, both $u(z)$ and $T(z)$ can be solved analytically and the analytical results agree well with the DNS data. Our careful characterisation of the mean velocity and temperature profiles in low-$Pr$ RBC provides a further understanding of the intricate interplay between the LSC, plume emission and boundary-layer dynamics, and pinpoints the physical mechanism for the emergence of a pronounced LSC in low-$Pr$ RBC.

We report a direct numerical simulation (DNS) study of the mean velocity and temperature profiles in turbulent Rayleigh-Bénard convection (RBC) at low Prandtl numbers (Pr). The numerical study is conducted in a vertical thin disk with Pr varied in the range 0.17 ≤ Pr ≤ 4.4 and the Rayleigh number (Ra) varied in the range 5 × 10 8 ≤ Ra ≤ 1 × 10 10 . By varying Pr from 4.4 to 0.17, we find a sharp change of flow patterns for the large-scale circulation (LSC) from a rigid-body rotation to a near-wall turbulent jet. We numerically examine the mean velocity equation in the bulk region and find that the mean horizontal velocity profile u(z) can be determined by a balance equation between the mean convection and turbulent diffusion with a constant turbulent viscosity ν t . This balance equation admits a self-similarity jet solution, which fits the DNS data well. In the boundary-layer region, we find that both the mean temperature profile T(z) and u(z) can be determined by a balance equation between the molecular diffusion and turbulent diffusion. Within the viscous boundary layer, both u(z) and T(z) can be solved analytically and the analytical results agree well with the DNS data. Our careful characterisation of the mean velocity

Introduction
Thermal convection at low Prandtl numbers is a common phenomenon both in nature and in many engineering applications. Here, the Prandtl number is defined as Pr = ν/κ, it represents the relative importance of the kinematic viscosity ν to the thermal diffusivity κ of the convecting fluid. For example, astronomical observations showed that the Sun has an outer convecting layer of low-density plasma, where the value of Pr is in the range of 10 −7 -10 −4 (Hanasoge, Gizon & Sreenivasan 2016;Schumacher & Sreenivasan 2020). Seismology measurements showed that the outer core of the Earth is a fluid layer made of a liquid iron-nickel alloy with Pr ∼ 10 −2 (Stevenson 1981). The convection of the Earth's outer core is linked to the generation and reversal of the Earth's magnetic field (Glatzmaier & Roberts 1995). The convection of the Earth's atmosphere with Pr 0.7 drives the large-scale atmospheric circulation and plays a critical role in the global climate and water balance (Hartmann, Moy & Fu 2001). Liquid-metal convection is also used in material processing (Brodova, Popel & Eskin 2001), nuclear engineering (Ihli et al. 2008), and liquid-metal batteries for renewable energy storage (Wang et al. 2014).
In the laboratory, controlled Rayleigh-Bénard convection (RBC) experiments are conducted in a closed cell, in which a fluid layer is heated from below and cooled from above. In addition to the Prandtl number Pr, the Rayleigh number Ra is another control parameter of thermal convection, which is defined as Ra = gα TH 3 /(νκ), where g is the gravitational acceleration, α is the thermal expansion coefficient of the fluid and T is the temperature difference across the fluid layer of height H. When the Rayleigh number is sufficiently large (e.g. Ra 10 8 for Pr 4.4), the bulk fluid becomes turbulent and a large-scale circulation (LSC) is formed across the convection cell (Krishnamurti & Howard 1981;Zocchi, Moses & Libchaber 1990). The LSC is driven by the warm and cold plumes emitted from the unstable thermal boundary layers near the bottom and top conducting plates and is maintained in a turbulent environment. This large-scale flow with Pr > 1 has been studied extensively in upright cylindrical cells of aspect ratio unity (Du & Tong 2000;Qiu & Tong 2001;Xi, Lam & Xia 2004;Sun, Xia & Tong 2005), in which the LSC has a single roll structure, with its size comparable to the cell height.
Another intriguing feature of turbulent RBC is that its thermal boundary layer has significant fluctuations resulting from intermittent eruption of thermal plumes from the boundary layer, even when the boundary layer is not fully turbulent (Du & Tong 2000;du Puits, Resagk & Thess 2010;Zhou & Xia 2010). The structure and dynamics of the thermal boundary layer and its interaction with the LSC are of great importance in determining the global heat transport of the system (Kadanoff 2001;Ahlers, Grossmann & Lohse 2009). For these reasons, the past decade has witnessed a continuing growth of experimental and theoretical efforts aimed at understanding the boundary-layer dynamics in turbulent RBC. Recent studies of the thermal boundary layer for Ra 10 12 and Pr 1 (Belmonte, Tilgner & Libchaber 1993, 1994Lui & Xia 1998;Du & Tong 2000 Shishkina, Horn & Wagner 2013;Zhou & Xia 2013;Scheel & Schumacher 2014) showed that the measured (and numerically calculated) normalised mean temperature profile θ(z) has a universal form θ(ξ ), where ξ ≡ z/δ T is the vertical distance from the conducting plate normalised by the thermal boundary-layer thickness δ T . The measured θ(ξ ) was found to be invariant with Ra and has the Prandtl-Blasius-Pohlhausen (PBP) form (Landau & Lifshitz 1987;Schlichting & Gersten 2016) for a laminar boundary layer (Grossmann & Lohse 2000) only when ξ is in the region ξ 0.6 (Zhou & Xia 2013). Deviations of θ(ξ ) from the PBP form were found when 0.6 ξ 4 (Scheel et al. 2012;Shi et al. 2012;Wagner et al. 2012;du Puits et al. 2013). More recently, Shishkina et al. (2015) considered the effect of boundary-layer fluctuations and obtained an analytical form of θ(ξ ) for the thermal boundary layers with Pr > 1, which explained the observed deviations of θ(ξ ) from the PBP form. The theoretical prediction made by Shishkina et al. (2015) was verified by the convection experiments conducted in a quasi-two-dimensional (quasi-2-D) thin-disk cell (Wang, He & Tong 2016). Wang et al. (2018) further extended the boundary-layer theory to the temperature variance profile.
Compared with the large number of investigations on the dynamics of the LSC and boundary layers for Pr 1, our understanding on the LSC and boundary-layer dynamics for fluids with Pr < 1 is rather limited. This is partially caused by the fact that laboratory experiments and numerical simulations of turbulent RBC at low Pr are quite challenging. On the experimental side, liquid metals, which are often used as a working fluid to obtain a sufficiently low value of Pr, are opaque and thus exclude optical imaging or particle tracking of the velocity field. Owing to their high thermal conductivity, it is also quite difficult to drive the low-Pr convection to reach sufficiently high Rayleigh numbers. On the numerical simulation side, because the Kolmogorov scale in the low-Pr fluid is smaller than the Batchelor scale by a factor of Pr 1/2 (Grötzbach 1983;Shishkina et al. 2010), a higher spatial resolution is needed to resolve the smaller Kolmogorov scale for low-Pr convection at a comparable value of Ra.
Early low-Pr experiments using liquid mercury (Pr 0.024) (Takeshita et al. 1996;Cioni, Ciliberto & Sommeria 1997;Mashiko et al. 2004;Tsuji et al. 2005) or gases (such as nitrogen and sulfur hexafluoride) and gas mixtures (0.18 ≤ Pr ≤ 0.88) (Hogg & Ahlers 2013) studied the Nusselt number (normalised heat flux) scaling, thermal boundary-layer profile, LSC dynamics and single-point velocity and temperature statistics in the bulk region. Recent low-Pr experiments using liquid gallium and its alloys (Pr 0.027) (Vogt et al. 2018;Zürner et al. 2019) investigated the three-dimensional (3-D) structure and dynamics of LSC in the Ra range 10 5 Ra 6 × 10 7 . These experiments found that the LSC at low Pr is more coherent compared with that at high Pr (> 1). Recent advances in computational power and numerical techniques allow the study of turbulent RBC at low Pr by direct numerical simulation (DNS) (Schumacher, Götzfried & Scheel 2015;, 2017Schumacher et al. 2016). It was found that the mean streamwise velocity is like a near-wall jet (Scheel & Schumacher 2017). Many of the experimental and DNS studies were conducted in upright cylindrical cells with an aspect ratio close to unity. In the cylindrical cells, the large-scale flow has several 3-D flow modes, such as the torsional and sloshing modes (Funfschilling & Ahlers 2004;Brown & Ahlers 2008Xi et al. 2009;Ji & Brown 2020), which may cause additional complications to the study of the LSC and boundary-layer dynamics. There are also corner flows in the closed cylinder (Sun et al. 2005), which may destabilise the large-scale flow (Sugiyama et al. 2010). The strong coupling between the boundary-layer dynamics and complex 3-D large-scale flow in a closed cylinder, which has been studied in recent numerical simulations (van Reeuwijk et al. 2008;Scheel et al. 2012;Shi et al. 2012;Stevens et al. 2012;Wagner et al. 2012;van der Poel et al. 2013;Shishkina et al. 2013;Scheel & Schumacher 2014), makes a quantitative comparison between the experiment and 2-D theory difficult.
In this paper, we report a systematic numerical study of the mean horizontal velocity profile and mean temperature profile for turbulent RBC in the low-Pr regime. DNS runs were performed in a vertical thin disk with its circular cross-section aligned parallel to gravity (see figure 1 for details). This specially designed quasi-2-D convection cell has two unique features for the DNS study attempted here. First, because the cell shape matches the single-roll structure of the LSC perfectly, there is no corner flow inside the cell. The LSC in the circular cross-section has a steady rotation along a fixed orientation. Owing to the axial symmetry of the thin disk cell, the initial orientation of the LSC is random for different initial conditions and different values of Ra. Second, because the flow is confined in a thin circular disk, no 3-D flow modes can be excited in this system. The quasi-2-D flow in the thin disk cell, therefore, has a better geometry satisfying the assumption of the boundary-layer theory for a 2-D flow. Even with these simplifications, the quasi-2-D system retains the key features of turbulent convection, which have been observed in the upright cylinders, and has been used in recent experimental and numerical studies of the LSC dynamics Song, Villermaux & Tong 2011;Song et al. 2014) and thermal boundary-layer profiles with Pr > 1 (Wang et al. 2016(Wang et al. , 2018. This is a 'simple but not simpler' convection system, which offers a natural platform to study the dynamics of the LSC and boundary layers, and their interactions in the low-Pr regime, and to test different theoretical ideas. In this DNS study, we vary the Rayleigh number in the range 5 × 10 8 ≤ Ra ≤ 1 × 10 10 and focus our attention on Pr = 0.17. This is the lowest value of Pr that can be reached by using a gas mixture of H 2 and Xe (Bajaj, Ahlers & Pesch 2002), and we choose this value of Pr so that the DNS results can be used for comparison with future experiment. This value of Pr is small enough so that the viscous and thermal boundary layers are well separated with the viscous boundary layer being nested within the thermal boundary layer. In this case, the LSC can produce a strong shear to the thermal boundary layer, which suppresses the emission of thermal plumes from the thermal boundary layer. On the other hand, the value of Pr used is not too small so that we can carry out the DNS study at high Rayleigh numbers with an adequate spatial resolution using the available computational resources.
The remainder of the paper is organised as follows. We first present the numerical method used and DNS set-up in § 2. A comparison of the temperature and velocity fields between the low-Pr and high-Pr regimes is given in § 3. The DNS results of the mean horizontal velocity and temperature profiles in the boundary-layer region are presented in § 4. The DNS results of the mean horizontal velocity profile in the bulk region are presented in § 5. Finally, the findings of this study are summarised in § 6.

Direct numerical simulation
The governing equations for turbulent RBC are the incompressible Navier-Stokes equations and the convective heat equation under the Boussinesq approximation. The dimensionless forms of these equations are given bŷ The length, time, velocity, pressure and temperature are made dimensionless by the cell height H, the free-fall time T f = √ H/(gα T), the free-fall velocity U f = √ gα TH, the free-fall pressure p f = ρgα TH and the temperature difference T across the cell, respectively. Figure 1(a) shows a sketch of the convection cell used for numerical simulations. The dimensionless boundary conditions are given bŷ Another control parameter is the aspect ratio Γ = d/H, where d is the thickness of the thin disk cell. The governing equations are solved numerically using the open-source code Nek5000 (Fischer 1997), which uses a spectral element method to accurately resolve the gradients in the velocity fieldû(r, t) and temperature fieldT(r, t). In the simulation, the time-derivative terms are discretised by backward-differentiation formula, the nonlinear convective terms are treated explicitly and the linear diffusive terms are approximated implicitly. This scheme leads to a Poisson equation for pressure and Helmholtz equations for velocity components and temperature. These equations are written in a weak formulation and discretised by the Galerkin method using the Nth-order Lagrangian interpolation polynomials as the basis functions on Gauss-Lobatto-Legendre (GLL) collocation points (Deville, Fischer & Mund 2002). More details of the numerical scheme and mesh resolution requirements can be found in Fischer (1997), Deville et al. (2002) and Scheel, Emran & Schumacher (2013). All the gradients in the post-processing are also calculated on the GLL collocation points with spectral accuracy.
The computational mesh is designed for the highest value of Ra = 1 × 10 10 with Pr = 0.17 and the minimal primary mesh size near the boundary is set at 1.944 × 10 −3 H, which is closed to the estimated viscous boundary-layer thickness. All the simulations are performed with the polynomial order N = 7, so that we have 8 3 = 512 grid points within each primary element. We verify that the mesh resolution at Ra = 1 × 10 10 and Pr = 0.17 satisfies the Grötzbach's criterion (Grötzbach 1983;Scheel et al. 2013). We use the same mesh for other simulations at lower values of Ra and with Pr in the range 0.1-4.4. This is done to save the programming time during the post-processing of the DNS results, and at the same time we have a sufficient mesh resolution for all simulations. An adaptive time step is applied to ensure that the Courant number is always below 0.5 during the simulations. The viscous and thermal boundary layers are well resolved by using the polynomial order N = 7. We run each simulation for at least 100T f to reach the steady state, followed by a continuous running for at least another 200T f to conduct time averaging. As shown in figure 1(b), the vertical profile of the local properties is computed along a thin column surrounding the vertical z-axis of the cell with x = y = 0 and is averaged over the cross-section of the thin column with a small area of 0.027 thickness square. Other details about the numerical simulations can be found in Wang et al. (2018). Table 1 lists a summary of the parameters used for simulations at different Rayleigh numbers and Prandtl numbers. In table 1, we also include the numerically calculated values of the local Nusselt number Nu, the Reynolds number Re, the normalised viscous boundary-layer thickness δ v /H and the normalised thermal boundary-layer thickness δ T /H. Here the local Nusselt number is defined as Nu = −∂ z T| z=0 H/ T and the Reynolds number is defined as Re = U m H/ν. The definition of the viscous and thermal boundary-layer thicknesses, δ v and δ T , are given in (4.4) and (4.18), respectively. For all the DNS runs, we use the same aspect ratio Γ = 0.2 and the same total number of spectral elements, N e = 96768.
3. Comparison of the temperature and velocity fields between the low-and high-Pr fluids Figure 2 shows the contour plots of the instantaneous 2-D temperature field on the middle cross-section for two Prandtl numbers: Pr = 0.17 and Pr = 4.4. It is found that warm rising plumes (in red) are generated from the bottom conducting plate and cold falling plumes (in blue) are generated from the top plate. These plumes drive the LSC. Because the low-Pr fluid has a relatively larger thermal diffusivity, the thermal plumes in the low-Pr fluid have a shorter lifetime. As a result, they have a lower chance to move into the bulk region and most of them are concentrated in the narrow region near the curved sidewall.
On the other hand, the thermal plumes in the high-Pr fluid have a longer lifetime and thus they have a higher chance to be found in the bulk region.
The spatial distribution of the thermal plumes has a profound influence on the velocity field. Figure 3 shows a comparison of the in-plane mean velocity field and corresponding mean pressure field across the middle cross-section of the cell between two Prandtl numbers: (a) Pr = 0.17 and (b) Pr = 4.4. It is seen that the mean flow field in the low-Pr fluid is more coherent and has a mean flow concentrated mainly in the narrow region near the circular sidewall. The pressure contour in the bulk region shows a homogeneous  circular shape, which fits nicely to the circular sidewall. This pressure field provides a long-range homogeneous shear flow along the conducting plates. Owing to the rotational symmetry of the LSC, the pressure gradient along the conducting plates is negligibly small. The thin disk cell thus offers a simple flow structure for the study of the intrinsic properties of the LSC and boundary layers at low Pr. The mean flow field in the high-Pr fluid, however, has more spatial variations and spreads deep inside the bulk region. The pressure contour in the bulk region shows a slightly tilted elliptical shape, which is caused by the accumulation of rising warm plumes on the lower right region and falling cold plumes on the upper left region. On the other hand, the mean velocity and pressure fields in convection cells of other shapes, such as thin square cells and cylinders of aspect ratio unity, are very different from those shown in figure 3. Previous numerical studies  (Xu 2014;Chong et al. 2018) have revealed that the mean velocity field in thin square cells has a large tilted LSC in the bulk region, which coexists with two small rolls of opposite rotation located in two opposing corners of the cell. Similar corner flows were also observed experimentally by particle image velocimetry in cylindrical cells (Sun et al. 2005). The interaction between the corner flow and LSC generates a complex structure near the conducting plates, which is different from the simple shear flow assumed by most boundary-layer theories. As shown in figure 3, no corner flow is observed in the thin disk cell.
To describe the velocity field more quantitatively, we show, in figure 4, the normalised mean horizontal velocity profile u(z)/U f as a function of the normalised vertical distance z/H for five different values of Pr. At Pr = 4.4, the bulk flow behaves like a rigid-body rotation with a zero mean velocity at the cell centre. The mean velocity increases linearly with the radial distance away from the cell centre. Similar single-roll structures of the LSC were also observed in the upright cylindrical cells of aspect ratio unity (Du & Tong 2000;Qiu & Tong 2001;Xi et al. 2004;Sun et al. 2005;. As Pr decreases, the effect of the fluid viscosity decreases and the rotation speed of the LSC increases. The faster rotation of the LSC gives the thermal plumes less time to penetrate into the bulk region. As a result, u(z)/U f moves to the region near the circular sidewall. At Pr = 0.17, the obtained u(z)/U f behaves like a near-wall jet with a sharp velocity peak near the bottom plate. The strong near-wall flow, in turn, produces a large entrainment effect carrying the thermal plumes in the narrow near-wall region. The strong coupling between the LSC and thermal plumes is self-organised so that a continuous LSC is maintained even when the thermal plumes have a short lifetime. The shape of the near-wall jet remains approximately the same when Pr is further reduced to 0.1. This result suggests that the flow has reached a new steady-state regime at low Pr.
In the following, we investigate the functional form of the normalised mean horizontal velocity profileũ(z) ≡ u(z)/U m and the normalised mean temperature profile θ(z) ≡ [T b − T(z)]/ b , where T(z) is the actual mean temperature profile, T b is the bottom plate temperature and b is the temperature difference across the bottom thermal boundary layer. This study is carried out at a fixed Pr = 0.17 and a fixed aspect ratio Γ = 0.2. As mentioned previously, Pr = 0.17 is chosen because the flow has reached a new steady-state regime at low Pr and the DNS results can be obtained more efficiently and be used for comparison with future experiments using a gas mixture of H 2 and Xe (Bajaj et al. 2002).

Mean horizontal velocity and temperature profiles in the boundary-layer region
To find the functional form of the mean horizontal velocity and temperature profiles in the boundary-layer region, we consider a 2-D convective flow over a long horizontal heating plate with a constant horizontal mean velocity U m . Using the coordinate system shown in figure 1(a) and taking the Reynolds decomposition of the velocity field u(r, t) = u (r) + u (r, t) and temperature field T(r, t) = T (r) + T (r, t), we obtain the viscous and thermal boundary-layer equations where all the time-average notation · on the first-order terms are omitted for conciseness. These equations are obtained by applying the boundary-layer approximations |w| |u|, |∂/∂x| |∂/∂z| and |∂ 2 /∂x 2 | |∂ 2 /∂z 2 |, and a long time average to the dimensional form of (2.2) and (2.3). Based on the boundary-layer approximations and dimensional analysis, we find that the convection and diffusion terms in the vertical direction are much smaller than those in the horizontal direction, as shown in (4.2). As a result, only the pressure gradient and buoyancy terms remain in the viscous boundary-layer equation in the vertical direction, namely, ρ −1 ∂ z p = gα(T − T 0 ) (Ching et al. 2019). As shown in figure 3(a) and more quantitatively in figure 5(a), the pressure gradient term ρ −1 ∂ x p in (4.2) is negligibly small, so that the vertical pressure-buoyancy balance equation does not  Figure 5. (a) Contributions of the mean convection term u∂ x u + w∂ z u, molecular momentum diffusion components −ν∂ 2 x u and −ν∂ 2 z u, Reynolds stress gradient components ∂ x u u and ∂ z w u , and pressure gradient 1/ρ∂ x p, as a function of the normalised vertical distance z/δ v . The unit of the terms is U 2 f /H. The two solid lines show the calculated terms of the velocity equation using the analytical solutionũ(ξ ) in (4.14) with a = 1.19. (b) Contributions of the mean convection term u∂ x T + w∂ z T, molecular thermal diffusion components −κ∂ 2 x T and −κ∂ 2 z T, and turbulent heat flux gradient components ∂ x u T and ∂ z w T , as a function of z/δ v . The unit of the terms is U f T/H. The DNS data used for the calculations shown in (a) and (b) are obtained at Pr = 0.17 and Ra = 1 × 10 9 in a vertical thin disk with Γ = 0.2. affect the above viscous and thermal boundary-layer equations explicitly. For low-Pr RBC, we still expect the classical boundary-layer assumptions to be valid for both the viscous and thermal boundary layers. This is because the thicknesses of the two boundary layers are comparable at Pr = 0.17, with the thickness ratio between the viscous and thermal boundary layers being Pr 1/2 0.4.
With the DNS data, we numerically calculate the contributions of each term in the mean velocity equation and the mean temperature equation. Figures 5(a) and 5(b) show, respectively, the terms of the velocity and temperature equations as a function of the normalised vertical distance z/δ v away from the bottom conducting plate. Here the viscous boundary-layer thickness δ v is defined as a vertical distance, at which the tangent of the mean velocity profile at the bottom plate (z = 0) intersects the maximum velocity U m , namely It can be seen from figure 5(a) that the boundary-layer approximations hold for the mean velocity equations in the near-wall region. We verify numerically that the horizontal pressure gradient ρ −1 ∂ x p is negligible at low Pr so that the viscous boundary-layer equation (4.2) is decoupled from the thermal boundary-layer equation (4.3) (Ching et al. 2019). In addition, we find numerically that the mean convection term, u∂ x u + w∂ z u, in (4.2) is negligibly small for this low-Pr RBC system. As a result, only two dominant terms remain in the mean velocity equation (4.2), namely, the molecular momentum diffusion ν∂ 2 z u and turbulent momentum diffusion ∂ z w u , which balance each other in the near-wall region, Similarly, the mean temperature equation (4.3) also has two dominant terms, namely, the molecular thermal diffusion κ∂ 2 z T and turbulent thermal diffusion ∂ z w T . There is a small contribution from the mean convection term, u∂ x T + w∂ z T, for Ra = 5 × 10 8 . From the DNS results, we find the contribution of the mean convection decreases with increasing Ra. As shown in figure 5(b), the mean convection can be approximately ignored for Ra = 1 × 10 9 . Therefore, in the Ra range 1 × 10 9 ≤ Ra ≤ 1 × 10 10 , we have (4.6)

Mean horizontal velocity profile in the boundary-layer region
With the turbulent viscosity ν t (z) defined as equation (4.5) becomes an ordinary differential equation Integrating both sides and applying the definition (4.4) as the boundary condition, we obtain The formal solution of (4.9) is given bỹ ds, (4.10) whereũ = u/U m and ξ = z/δ v . With the no-slip boundary condition at the bottom conducting plate, we find the Reynolds stress w u and its derivatives satisfy the following boundary conditions at From the definition of ν t in (4.7) and the linear velocity profile at the wall, i.e. du/dz ∝ z 0 , we obtain the boundary conditions for ν t at ξ = 0, Therefore, the leading order of ν t has the form where a is a constant. A similar boundary condition was also obtained for κ t /κ in the thermal boundary-layer equation (Shishkina et al. 2015). Substituting (4.13) into (4.10), we obtain the analytical form of the normalised horizontal velocity profilẽ .
Figure 6(a) shows the normalised turbulent viscosity ν t /ν as a function of the normalised vertical distance z/δ v . Using (4.7), we numerically calculate ν t = − w u /(du/dz). For all three values of Ra, the obtained ν t /ν is well described by (4.13) with a = 1.19 (green dashed line) for small values of z up to z/δ v 2. Deviations from the green dashed line are observed when z/δ v 2. These deviations exhibit an interesting dependence on Ra, with the obtained ν t /ν at Ra = 1 × 10 9 as an exceptional case, which has a wider power-law range up to z/δ v 6. With the obtained ν t /ν in figure 6(a), we numerically compute the solution of (4.10), as shown by the three solid lines in figure 6(b) for the three different values of Ra. It is seen that the agreement between the numerical solutions of (4.10) and the DNS data (open symbols) improves significantly with increasing Ra. A small deviation between the green solid line and black circles is observed at larger values of z/δ v for Ra = 5 × 10 8 . This is caused by the fact that for small values of Ra, a small contribution from the mean convection term becomes noticeable, which affects the accuracy of (4.5). Because the obtained ν t /ν at Ra = 1 × 10 9 is well described by (4.13) with a = 1.19 over a wider range of z/δ v , the corresponding velocity profile u/U m as a function of z/δ v (red triangles in figure 6b) is equally well described by (4.14) with the same fitting parameter a = 1.19. By assuming (4.13) holds for distances far enough from the wall, one may use the asymptotic boundary condition,ũ(∞) = 1, to determine the value of a = 2 √ 3π/9 1.2. The fitted value of a = 1.19 is very close to this limiting value. Figure 6 thus confirms that the viscous boundary-layer model discussed earlier captures the essential physics.

Mean temperature profile in the boundary-layer region
Using the turbulent thermal diffusivity κ t (z) defined as (4.15) equation (4.6) becomes an ordinary differential equation Integrating both sides and applying the relation at the boundary, which defines the thermal boundary-layer thickness δ T , we obtain The formal solution of (4.18) is given by For high-Pr RBC, the thermal boundary layer is nested within the viscous boundary layer, so that κ t (z) obeys a cubic power law similar to (4.13) across the entire thermal boundary layer. In this case, the mean temperature profile θ(ξ ) takes the same functional form as shown in (4.14). This was shown previously by Shishkina et al. (2015). For low-Pr RBC, however, the thermal boundary layer is thicker than the viscous boundary layer, so that the LSC may produce a strong shear on the outer portion of the thermal boundary layer, as illustrated in figure 7(a). The numerically calculated κ t (z) for Pr = 0.17 is found to obey the cubic power law in (4.13) only up to z/δ v 1. This is shown in figure 7(b). Outside the viscous boundary layer (but still within the thermal boundary layer), the obtained κ t (z) goes as z 2.2 with the power-law exponent being smaller than 3. In this case, no analytical solution of (4.19) is available and one needs to numerically integrate (4.19). The situation is somewhat similar to a boundary layer near a partial-slip wall, as illustrated in figure 7(a). Owing to the large streamwise velocity at the edge of the viscous boundary layer, ∂ 2 z w T is no longer guaranteed to be zero, so that the leading order of κ t is reduced from z 3 to z 2 (authors' unpublished observations).
With the obtained κ t /κ in figure 7(b), we numerically compute the solution of (4.19), as shown by the three solid lines in figure 8 for the three different values of Ra. It is seen that the agreement between the numerical solutions of (4.19) and the DNS data (open symbols) improves significantly with increasing Ra. Small deviations are observed at larger values of z/δ v for the two lower values of Ra. This is caused by the fact that for small values of Ra, a small contribution from the mean convection term becomes noticeable (see figure 5b), which affects the accuracy of (4.6). The good agreement between the calculated solution of (4.19) and the DNS result at Ra = 1 × 10 10 further confirms that the thermal boundary-layer model discussed previously captures the essential physics.

Mean horizontal velocity profile in the bulk region
Figure 9(a) shows the terms of the velocity equation as a function of the normalised vertical distance z/H. Compared with those terms in the boundary-layer region, all of the terms are small in the bulk region. Among them, the mean convection and the vertical Ra = 5 × 10 8 Ra = 1 × 10 9 Ra = 1 × 10 10 . Normalised mean temperature profile θ(z) as a function of z/δ v for Ra = 5 × 10 8 (black circles), 1 × 10 9 (red triangles) and 1 × 10 10 (blue diamonds). The green, blue and brown solid lines show, respectively, the numerical solutions of (4.19) using the numerically calculated κ t /κ shown in figure 7(b) for Ra = 5 × 10 8 , 1 × 10 9 and 1 × 10 10 . The DNS data used for the calculations are obtained at Pr = 0.17 in a vertical thin disk with Γ = 0.2.
Reynolds stress gradient decay with z/H much slower than the other terms. By balancing the mean convection with the vertical Reynolds stress gradient, we have As shown in figure 9(b), the turbulent viscosity ν t (x, z) in the bulk region is approximately independent of z and, therefore, we have  Figure 9. (a) Contributions of mean convection u∂ x u + w∂ z u (black circles), vertical molecular momentum diffusion −ν∂ 2 z u (blue diamonds) and vertical Reynolds stress gradient ∂ z w u (brown triangles) as a function of the normalised vertical distance z/H for Ra = 1 × 10 9 . The unit of the terms is U 2 f /H. (b) Normalised turbulent viscosity ν t (z)/(U f H) as a function of z/H for three values of Ra: 5 × 10 8 (black circles), 1 × 10 9 (red triangles) and 1 × 10 10 (blue diamonds). The DNS data used for the calculations shown in (a) and (b) are obtained at Pr = 0.17 in a vertical thin disk with Γ = 0.2. Substituting (5.2) into (5.1), we obtain Equation (5.3) admits a self-similarity solution for a free turbulent jet (Pope 2000). By assuming the large-scale flow in the bulk region is a half turbulent jet, we find u = 1 − tanh 2 z − z 0 η , (5.4) whereũ = u/U m , z 0 is the starting position of the half-jet, and η represents the jet width, (5.5) Figure 10(a) shows the normalised mean horizontal velocity profile u(z)/U f as a function of z/H. It is seen that all of the numerically calculated u(z)/U f for different values of Ra collapse approximately on to a single master curve. This suggests that the obtained u(z)/U f in the bulk region is approximately a scaling function of z/H, which is invariant with Ra. From figures 9(b) and 10(a), we find ν t 1.0 × 10 −4 U f H and U m 0.42U f . By assuming the streamwise length x scales with the circumference πH of the circular disk, i.e.x πH, we obtain the jet width η 0.055H from (5.5). This result thus explains why the numerically calculated u(z)/U f in the bulk region is a scaling function of z/H. In addition, figure 10(b) shows that the obtained u(z)/U f is adequately described by (5.4) and the fitted values of η is 0.08H, which also agrees with the estimated value of η 0.055H by (5.5). The good agreement between the DNS results and the analytical solution confirms that the turbulent jet model discussed above captures the essential physics. We also compute u(z)/U f in the Γ = 0.1 disk cell and the resulting u(z)/U f exhibits a similar shape of near-wall jet, as shown in figure 10 Figure 10. (a) Normalised mean horizontal velocity profiles u(z)/U f as a function of the normalised vertical distance z/H for three different values of Ra: 5 × 10 8 (black circles), 1 × 10 9 (red triangles) and 1 × 10 10 (blue diamonds). (b) Replot of the mean horizontal velocity profiles u/U m normalised by its maximum value U m as a function of z/H for Ra = 5 × 10 8 (black circles), 1 × 10 9 (red triangles) and 1 × 10 10 (blue diamonds). For clarity, the origin of the red and blue curves is shift to the right by 0.1 and 0.2 normalised distances, respectively. The green, blue and brown solid lines show the same plot of (5.4) with η = 0.08H and z 0 = 0.02H. The DNS data used for the calculations shown in (a) and (b) are obtained at Pr = 0.17 in a vertical thin disk with Γ = 0.2. is comparable to the disk thickness, we suspect that this weak Ra-dependence is caused by the confinement effect of the thin disk. Indeed, when the disk thickness is increased to Γ = 0.2, the obtained u(z)/U f is found to have an approximate universal form, as shown in figure 10(a), which does not change significantly with Ra.

Summary
We have carried out a systematic numerical study of the mean velocity and temperature profiles in turbulent Rayleigh-Bénard convection at low Prandtl numbers. The DNS runs were conducted in a vertical thin disk with the Prandtl number (Pr) varied in the range 0.17 ≤ Pr ≤ 4.4 and the Rayleigh number (Ra) varied between 5 × 10 8 and 1 × 10 10 . For a fixed value of Ra = 5 × 10 9 and with varying values of Pr from 4.4 to 0.17, we find a sharp change of flow patterns for the LSC. For high-Pr RBC, the lifetime of thermal plumes is longer so that they have a better chance to enter the bulk region and drive the flow globally. As a result, the LSC in high-Pr RBC rotates like a rigid body. For low-Pr RBC, however, the thermal plumes have a shorter lifetime so that they have less chance of surviving in the bulk region. Consequently, the thermal plumes are concentrated in the near-wall region and the LSC in low-Pr RBC behaves like a near-wall turbulent jet.
For low-Pr RBC, the viscous boundary layer is nested inside the thermal boundary layer. Thus, the strong shearing effect by the viscous boundary layer causes the thermal boundary layer to be separated into two sub-layers. This effect makes the thermal boundary layer at low Pr more complicated than that at high Pr. We numerically examine the mean velocity and temperature equations in the boundary-layer region and find that both the mean horizontal velocity profile u(z) and temperature profile T(z) can be determined by a balance equation between the molecular diffusion and turbulent diffusion. Furthermore, the no-slip boundary condition at the bottom conducting wall dictates that both the turbulent viscosity ν t (z) and turbulent thermal diffusivity κ t (z) as a function of distance z away from the bottom conducting wall obey a cubic power law (∼ z 3 ) for small values of z within the viscous boundary layer. In this case, both u(z) and T(z) can be solved analytically and the solutions agree well with the DNS data. Owing to the large streamwise velocity at the edge of the viscous boundary layer, the leading order of the obtained κ t (z) at larger values of z outside the viscous boundary layer (but still within the thermal boundary layer) is reduced from z 3 to z 2 . In this case, the mean temperature profile T(z) can be solved numerically using the calculated κ t (z).
For a fixed value of Pr = 0.17 and with varying Ra in the range 5 × 10 8 ≤ Ra ≤ 1 × 10 10 , we find that the mean horizontal velocity profile u(z) can be determined by an equation with the mean convection balanced by turbulent diffusion with a z-independent turbulent viscosity ν t . This balance equation admits a self-similarity jet solution, which fits the DNS data well. The width η of the near-wall jet is determined by the diffusion equation η 2 = 4ν t τ , where τ πH/U m is the turnover time of the LSC with πH being the circumference of the circular disk and U m is the maximum velocity of the LSC. Our DNS data also reveal that the azimuthal velocity profiles at different polar angles θ have a similar near-wall jet shape, as shown in figure 10. However, the peak velocity U m has a weak angular dependence. We believe that this is caused by the angular dependence of the buoyancy term, g sin θ, where g is the gravitational acceleration. Owing to this buoyancy term, the momentum equation is coupled to the thermal equation, which makes the problem difficult to solve mathematically. Here we choose the vertical z-direction with θ = 0, so that the two equations are decoupled. In this case, we find a simple mathematical solution to describe the vertical profile of the mean horizontal velocity u(z), which captures the essential feature of the LSC in low-Pr RBC.
Our work thus provides a full characterisation of the mean velocity and temperature profiles in low-Pr RBC. This characterisation allows us to further understand the intricate interplay between the LSC, plume emission and boundary-layer dynamics, and pinpoint the physical mechanism for the emergence of a pronounced LSC in low-Pr RBC. This study lays down a foundation for further study of the effects of cell shape and confinement by cell end walls, as reported by Chong et al. (2018) for a different cell geometry.