Heat transport and flow morphology of geostrophic rotating Rayleigh–Bénard convection in the presence of boundary flow

Abstract

Using direct numerical simulations, we investigate the heat transport in bulk and boundary flows separately in rotating Rayleigh–Bénard convection in cylindrical cells. In the bulk we observe a steep scaling relationship between the Nusselt number ($Nu$) and the Rayleigh number ($Ra$), which is consistent with the results from simulations using periodic boundary conditions. For the boundary flow, we observe a power law $Nu_{BF}\sim (Ra/Ra_w)^1$ at the leading order, where $Nu_{BF}$ is the local Nusselt number of the boundary flow and $Ra_w$ is the onset Rayleigh number of the wall mode. We develop a model using the boundary layer marginal stability theory to explain this power law, and further show that a more precise description of the data can be obtained if a higher-order correction is introduced. A striking finding of our study is the observation of a sharp transition in flow state, manifested by a sudden drop in $Nu_{BF}$ with a corresponding collapse of the boundary flow coherency. After the transition, the boundary flow breaks into vortices, leading to a reduction in flow coherency and heat transport efficiency. As the physical properties of the vortices should not depend on the aspect ratio, $Nu_{BF}$ for all aspect ratios collapse together after the transition. Moreover, the centrifugal force helps trigger the breakdown of the coherent boundary flow state. For this reason, $Nu_{BF}$ for the cases with non-zero centrifugal force collapse together. We further develop a method that enables us to separate the contributions from the bulk and boundary flows in the global Nusselt number using only the global $Nu$ and it does not require the centrifugal force to be absent.

1. Introduction

Natural convection under the influence rotation occurs widely in many geo- and astrophyiscal environments (Chen & Guo 2003; King et al. 2009; Cheng et al. 2015, 2018; Plumley & Julien 2019). For example, the large-scale vortices observed in the giant planets, such as the ‘Great Red Spot’ at the atmosphere of Jupiter, are sustained by convection under a strong Coriolis force (Guervilly, Hughes & Jones 2014; Favier, Guervilly & Knobloch 2019; Cai 2021). Moreover, moist convection is considered an important factor in explaining the formation of a polar vortex in giant planets (O'Neill, Emanuel & Flierl 2015) and the upscale energy transfer strengthening the large cyclones at Jovian high latitudes (Siegelman et al. 2022). Understanding the dynamics of these systems relies deeply on the observations. However, owing to the difficulty in data acquisition and the complexity of these systems, idealized models are very useful complementary approaches used to investigate the dynamics of these systems. A classical physical model for rotating convection is the rotating Rayleigh–Bénard convection (RBC). When the system is rapidly rotating, the horizontal pressure gradient is balanced by the Coriolis force at the leading order, which is known as the geostrophic balance. Rayleigh–Bénard convection under geostrophic balance possesses features relevant to many important issues of the geo- and astrophysical systems, thus, it has attracted considerable interest in recent decades (Knobloch 1998; King et al. 2009; Stevens et al. 2009; Zhong et al. 2009; Weiss et al. 2010; King, Stellmach & Aurnou 2012; Ecke & Niemela 2014; Stellmach et al. 2014; Kunnen 2021). In this system, the rotation rate can be characterized by the Ekman number $Ek\equiv \nu /(2\varOmega H^2)$, and the strength of buoyancy forcing is represented by the Rayleigh number $Ra\equiv (\alpha g\Delta H^3)/(\kappa \nu )$. Here $\alpha$ is the thermal expansion coefficient of the fluid, $g$ is the gravitational acceleration, $\varDelta =T_{hot}-T_{cold}$ is the temperature difference between the hot and cold plates, $H$ is the height of convection cell, $\nu$ is the kinetic viscosity, $\kappa$ is the thermal diffusivity and $\varOmega$ is the rotating rate. Rotation suppresses the vertical fluid motions as a result of the Taylor–Proudman constraint. For a laterally unlimited domain, bulk convection emerges when $Ra$ is beyond the onset Rayleigh number $Ra_c$ (Niiler & Bisshopp 1965)

(1.1)\begin{equation} Ra_c=(8.7-9.63Ek^{1/6})Ek^{{-}4/3}. \end{equation}

For the systems with lateral boundaries, however, distinct flow structures can be found in the region adjacent to the sidewalls. These structures are called wall mode, which is azimuthally periodic and precesses in the retrograde direction (Zhong, Ecke & Steinberg 1991; Ecke, Zhong & Knobloch 1992; Favier & Knobloch 2020). Different from the direction of precession, the mean azimuthal velocity of the outer part of the wall mode is positive (prograde). For a cylindrical system with non-slip boundaries, the onset Rayleigh number of the wall mode is given by (Zhang & Liao 2009)

(1.2)\begin{equation} Ra_w={\rm \pi}^2\sqrt{6\sqrt{3}}Ek^{{-}1}+46.6Ek^{{-}2/3}. \end{equation}

For a large $Pr$ number, the wall mode emerges at a smaller $Ra$ than the bulk convection. When the bulk flow emerges, the wall-localized flow structure is also understood as the boundary zonal flow (Zhang et al. 2019; Zhang, Ecke & Shishkina 2021; Ecke, Zhang & Shishkina 2021; Wedi et al. 2022). To avoid confusion, hereafter we refer to this flow structure adjacent to the sidewall as the boundary flow, whether or not the bulk convection sets in.

Among the various topics in geostrophic RBC, determination of the scaling relationship of heat transport efficiency is one of the key issues. This is not only attributed to the fact that heat transport itself is crucial in many physical problems, heat transport may also reflect the transition in flow morphology and the dominating force balance in the system. For this reason, the scaling relationship of the heat transport efficiency can help one to better understand the essential physics of the system, which may also enable simplification of the physical models or parameterization of the main physical quantities of the system. Moreover, the scaling relationship of heat transport efficiency provides the possibility to estimate the physical quantities or explain the features of the geo- and astrophysical environments that are beyond the achievable parameter range in experiments or simulations. The heat transport efficiency is quantified by the Nusselt number $Nu\equiv (\langle u_zT-\kappa \partial T/\partial z\rangle _{x,y})/(\kappa \varDelta /H)$, where $\langle {\cdot }\rangle _{x,y}$ denotes the average over the horizontal plane. For laterally unbounded cases, a steep scaling relationship $Nu\sim (Ra/Ra_c)^3$ can be found numerically (King et al. 2012; Cheng et al. 2015). However, recent experiments and numerical simulations using cylindrical cells observed much smaller scaling exponents than 3 (Lu et al. 2021). Such a difference suggests that the boundary flow has a distinct $Ra$ dependence of the heat transport efficiency from the bulk. In practical situations, the existence of lateral boundaries is inevitable. Moreover, in order to reach the geostrophic regime, many experiments are chasing for even smaller $Ek$ and higher $Ra$. One can enhance $Ra$ by using the so-called annular centrifugal Rayleigh–Bénard convection (ACRBC) system through the generation of supergravity induced by rapid rotation and the strong centrifugal force (Jiang et al. 2020; Wang et al. 2022), or simply by increasing the system height $H$. However, for the ACRBC system, the rotating axis is perpendicular to the effective gravity (centrifugal force), which is different from the rotating RBC system (rotating axis is parallel to the gravity). Thus, increasing $H$ is a more reasonable method to increase $Ra$ in a rotating RBC system. On the other hand, increasing $H$ or the rotation rate $\varOmega$ can both effectively extend the lower limit of the achievable $Ek$, but for the latter, the centrifugal effect could be strong so that the flow structures are significantly changed (Horn & Aurnou 2018; Hu et al. 2021; Hu, Xie & Xia 2022; Hu & Xia 2023). For this reason, in a rotating RBC system, it is believed to be advantageous to extend the parameter range by utilizing slender cells with a large system height (Cheng et al. 2015, 2020; de Wit et al. 2020). In this case, the influence of boundary flows can be significant. Thus, it is crucial to quantitatively understand the influence of boundary flow, so that one can assess and even extract the contributions of the boundary flow from the total heat flux.

In this study we use direct numerical simulations (DNS) to investigate rotating RBC in cylindrical cells with various aspect ratios, aiming to obtain a systematic understanding of the heat transport in rotating RBC in the presence of the boundary flows. To achieve this objective, we examine the bulk and boundary flows separately. The remainder of this paper is organized as follows. In § 2.1 we briefly introduce the numerical method used in this study. In § 2.2 we explain how the bulk and boundary flows are decomposed in our study, and demonstrate the reliability of our method. In §§ 3.1 and 3.2 we discuss the heat transport efficiency of the bulk and boundary flows separately. In § 3.3 we report a sharp transition in flow morphology, which is manifested as the collapse of coherency of the boundary flow structure and accompanied with a sudden drop in heat transport efficiency associated with the boundary flow. In § 3.4 we take a whole-system approach and provide a unifying understanding of the global heat transport efficiency.

2. Methodology

2.1. Numerical set-up

We consider a rotating cylindrical RBC cell with an upward rotational vector. The top and bottom plates are non-slip and with constant temperature $T_{cold}$ and $T_{hot}$, respectively. The sidewall is non-slip and adiabatic. The geometrical property of the cell is described by the aspect ratio $\varGamma =D/H$, where $D$ is the diameter of the convection cell. This system is governed by the three-dimensional Navier–Stokes equations with Oberbeck–Boussinesq approximation. Since the centrifugal effect is usually unavoidable in many of the rapid rotating experiments, we have conducted a series of simulations that consider the centrifugal effect to examine the influence of centrifugal force. The centrifugal force is quantified by the Froude number $Fr\equiv (R\varOmega ^2)/g$, which is the ratio of centrifugal acceleration over the gravitational one. Here $R\equiv D/2$ is the radius of the cell. Length, velocity, time and temperature are respectively non-dimensionalized by the height of the convection cell $x_{ref}=H$, the convective free-fall velocity $u_{ref}=(\alpha gH\varDelta )^{1/2}$, the free-fall time $t_{ref}=x_{ref}/u_{ref}$ and the temperature difference between the hot and cold plates $T_{ref}=\varDelta$. The dimensionless governing equations with centrifugal effect are

(2.1)\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}={-}\boldsymbol{\nabla} p+\sqrt{\frac{Pr}{Ra}}\nabla^2\boldsymbol{u}-\frac{1}{Ek}\sqrt{\frac{Pr}{Ra}}\hat{e}_z\times\boldsymbol{u}+\theta \hat{e}_z-Fr\frac{2r}{\varGamma}\theta\hat{e}_r, \end{gather}

(2.2)\begin{gather}\frac{\partial \theta}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \theta=\frac{1}{\sqrt{RaPr}}\nabla^2\theta, \end{gather}

(2.3)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0. \end{gather}

Here $Pr\equiv \nu /\kappa$ is the Prandtl number, which is the ratio between momentum and thermal diffusivity of the working fluid. All our simulations are done at $Pr=4.38$ (corresponding to water at about ${40}\,^{\circ }{\rm C}$). The dimensionless temperature variation $\theta$ is defined as $\theta \equiv (T-T_m)/\varDelta$, where $T_m\equiv (T_{hot}+T_{cold})/2$. For the sake of simplicity, symbols without special notation stand for the non-dimensional quantities in the rest of this paper.

The governing equations are solved in the cylindrical coordinates using a well-tested DNS code called CUPS, which is a fully parallelized DNS code based on the finite volume method with fourth-order precision (Kaczorowski & Xia 2013; Kaczorowski, Chong & Xia 2014; Chong, Ding & Xia 2018). In this study we focus on the systems with $Pr>1$, in which case the characteristic length scale of the temperature field (Bachelor length scale) is smaller than that of the velocity field (Kolmogorov length scale). In a traditional single-resolution scheme, the grid spacing is determined by Bachelor length scale, which is over-resolved for the velocity field. To improve computational efficiency without any sacrifice in precision, we used a multiple-resolution strategy. In this case, the momentum and temperature equations are solved in two grid sets. The momentum equations are solved in a coarser grid set than the temperature equation, which can still allow Batchelor and Kolmogorov length scales both being resolved. Using this algorithm, computational sources spent on the momentum solver is significantly reduced. We conducted three series of simulations: set I varies $Ra$ and fixes $Ek$ at $Ek=1.85\times 10^{-6}$; set II has similar conditions as set I, but includes the centrifugal effect; and in set III, $Ra$ is fixed at $Ra=2\times 10^{7}$ and $Ek$ is varied. Three aspect ratios ($\varGamma$) are used in sets I and II, which are $\varGamma =0.5$, $1$ and $2$. In set III we fix $\varGamma =4$. Additionally, we also conduct a series of simulations in a cubic cell with aspect ratio 1 and lateral periodic boundaries. These simulations use the same Ekman number as sets I and II. All statistical quantities are taken over 400 convective free-fall time units and, for those with a high rotating rate, the statistical period is longer. In the main text of this paper, we focus on sets I and II with fixed $Ek$. In Appendix F the cases with fixed $Ra$ and the data of set III are discussed. Table 1 in the supplemental material available at https://doi.org/10.1017/jfm.2023.872 lists all the numerical set-ups for the simulations in this study, and the parameters for the cases in sets I and II are shown in figure 1.

Figure 1. Parameter space for the cases in set I and II. In these two data sets, $Ek$ is fixed to be $1.85\times 10^{-6}$ and $Pr=4.38$. The black dashed and grey dot-dashed lines, respectively, correspond to $Ra=Ra_c$ and $Ra=Ra_w$.

2.2. Decomposition of bulk and boundary flows

The key to quantitatively understand the global heat transport of this system is to decompose the corresponding contributions from the bulk and the boundary flows. In previous studies the decomposition was achieved using a series of different radii as the cutoff (Ecke, Zhang & Shishkina 2022). In this study we wish to decompose the respective contributions from these two regions according to the distinct features of these two regions. In order to develop a reliable and reasonable method for the decomposition across various situations with differing turbulence intensities, we first concentrate on the flow structures present within this system. In figure 2(a) we present the distributions of the azimuthal velocity $\langle u_\phi \rangle _{t,\phi }$ averaged over time and the azimuthal direction for $Ra=8.71\times 10^8$. From figure 2(a) one can see flow structures with positive azimuthal velocity near the sidewall. It would be intuitive to separate the boundary flow and the bulk flow according to the zero crossing of $u_\phi$ from the outer positive part close to the sidewall to the inner negative part. However, such wall-localized boundary flow is a two-layer structure in which the outer and inner layers have opposite azimuthal velocities. Thus, using the simple zero-crossing criterion cannot extract the boundary flow in its entirety. We remark that this two-layer structure is consistent with the analytical solution of the velocity component parallel to the sidewall for the wall mode given by Herrmann & Busse (1993). Such similarity suggests the connection between the wall mode and the boundary flow. This two-layer structure can also be found in the averaged convective heat flux $\langle J_z\rangle _{t,\phi }$ shown in figure 2(c), where $J_z$ is the local convective heat flux defined as $J_z\equiv \sqrt {RaPr}u_z\theta$. We see that the region adjacent to the sidewall and with positive azimuthal velocity has positive convective heat flux; and the region next to it and with negative $u_\phi$ has negative $J_z$.

Figure 2. (a) Mean azimuthal velocity. (b) Vertical averaged r.m.s. of the azimuthal velocity $\langle u^2_\phi \rangle ^{1/2}_{t,z,\phi }$. (c) Mean convective vertical heat flux. (d) The convective vertical heat flux averaged temporally and in the vertical and azimuthal directions $\langle J_z\rangle _{t,z,\phi }$. The insets of (b,d) show the corresponding data from $r/R=0$ to $1$, respectively. The control parameters for (a–d) are $Ra=8.71\times 10^8$, $Ek=1.85\times 10^{-6}$, $Fr=0$ and $\varGamma =1$. The dot-dashed and dashed lines in (a–d) respectively correspond to $r_{BF}$ and $r_0$. (e) Width of the boundary flow $\delta _{BF}/H$ (close symbols) and $\delta _{0}/H$ (open symbols) as a function of $Ra$. The blue dot-dashed and red dashed lines in (e) represent the best fits of the data for $\delta _{BF}$ and $\delta _0$, respectively.

To better identify the two-layer structure of the boundary flow and explain the method decomposing the bulk and the boundary flows, we present the root-mean-square (r.m.s.) value of the azimuthal velocity $\langle u^2_\phi \rangle ^{1/2}_{t,z,\phi }$ in figure 2(b) and the mean convective heat flux $\langle J_z\rangle _{t,z,\phi }$ in figure 2(d). The subscript $\langle {\cdot }\rangle _{t,z,\phi }$ denotes the average over time and the vertical and azimuthal directions. We remark that the azimuthal component of the boundary flow is azimuthally periodic and its flow strength cannot be correctly reflected from $\langle u_\phi \rangle _{t,z,\phi }$, thus decomposing the bulk and the boundary flow according to $\langle u_\phi \rangle _{t,z,\phi }$ can be problematic especially at high $Ra$ when the bulk fluctuations are strong. For this reason, we instead use the r.m.s. $\langle u^2_\phi \rangle ^{1/2}_{t,z,\phi }$ in the bulk-boundary-flow decomposition. From the profile of $\langle u^2_\phi \rangle ^{1/2}_{t,z,\phi }$ presented in figure 2(b), one can clearly see a two-layer structure near the boundary and a bulk with a relatively weak flow strength. We fit the data near the two peaks using parabolic functions, which are shown as the red dot-dashed and blue dashed curves in figure 2(b), respectively. We fit the bulk part of $\langle u^2_\phi \rangle ^{1/2}_{t,z,\phi }$ ($r/R<0.7$) using linear functions, which are denoted by the grey dotted lines in figure 2(b). We then define the left intersection between the parabolic fitting of the inner peak and the linear fitting as $r_{BF}$ (indicated by the red triangles), which is the edge separating the bulk and the boundary flow, as shown by the vertical black dot-dashed lines in figure 2. The blue and red shaded regions in figure 2(b,d) respectively represent the bulk and the boundary flows. We then define $r_0$ as the zero crossing of the outer parabolic fitting (indicated by the blue stars), which separates the prograde (outer) and retrograde (inner) parts of the boundary flow. As heat transport is one of the main topics of this study, we need to ensure that the method we used for the decomposition can effectively identify the different domains with distinguishing features especially in terms of the heat transfer. To examine our method, we also present the mean convective heat flux $\langle J_z\rangle _{t,z,\phi }$ and the result of the decomposition in the same plot, as shown in figure 2(d). We see that the boundary flow with distinguishable features in the convective heat flux can be properly identified, which provides certain validation for our method. More examples of the decomposition for various $Ra$ and different flow states are presented in Appendix A.

The width of the entire boundary flow $\delta _{BF}\equiv \varGamma /2-r_{BF}$ and the outer boundary flow $\delta _{0}\equiv \varGamma /2-r_0$ can then be determined, respectively. We present the results of $\delta _{BF}$ and $\delta _0$ for sets I and II in figure 2(e). We find that the results of different $\varGamma$ collapse together for both $\delta _{BF}$ and $\delta _0$. Interestingly, it seems that $\delta _{BF}$ and $\delta _0$ exhibit similar $Ra$ dependence but only differ by a constant, which is found to be $\delta _{BF}\approx 2.7\delta _0$ in our study. Zhang et al. (2021) observed a scaling relationship $\delta _0\sim Ra^{1/4}$ using $Pr=0.8$. More recently, the experiments conducted by Wedi et al. (2022) suggests that for $Pr=6.55$, $\delta _0$ is insensitive to $Ra$. The numerical study conducted by Ecke et al. (2022) also has observed a relatively weak $Ra$ dependence of $\delta _0$, which is given by $\delta _0/H\sim (Ra-Ra_w)^{1/6}Ek^{2/3}$. We plot the black and red dashed lines respectively referring to $\delta _{BF,0}\sim Ra^{1/4}$ and $\delta _{BF,0}\sim Ra^{0}$ in figure 2(e) for comparison. For $Ra\lesssim 1\times 10^9$, both $\delta _{BF}$ and $\delta _{0}$ are insensitive to $Ra$, which agrees with $\delta _{BF,0}\sim Ra^{0}$ indicated by Wedi et al. (2022). When $Ra$ approaches to about $1\times 10^9$, $\delta _{BF}$ and $\delta _0$ slightly increases with $Ra$. Even though the data are scattered for $Ra\gtrsim 1\times 10^9$, it seems that $\delta _{BF,0}\sim Ra^{1/4}$ deviates from the data of either $\delta _{BF}$ or $\delta _0$. As pointed out by Wedi et al. (2022), the exponent of the power law for $\delta_0$ with $Ra$ could depend on $Pr$. Thus, such a discrepancy between our study and Zhang et al. (2021) may be attributed to the difference in Prandtl numbers. Nevertheless, for $Ra\lesssim 1\times 10^9$, it is reasonable to consider that $\delta _{BF}$ and $\delta _0$ are not sensitive to $Ra$.

3. Results

3.1. Heat transport properties of the bulk

Using the criterion discussed in § 2.2, we proceed to investigate the heat transport in the bulk and boundary flows separately. We define $Nu_{bulk}$ as the Nusselt number of the bulk

(3.1)\begin{equation} Nu_{bulk}\equiv\left\langle \sqrt{RaPr}u_z\theta-\frac{\partial\theta}{\partial z}\right\rangle_{t,z,\phi,r\leqslant r_{BF}}. \end{equation}

We plot $Nu_{bulk}$ as a function of $Ra/Ra_c$ in figure 3(a). For comparison, we also present $Nu$ from simulations using lateral periodic boundary conditions. For $Ra$ smaller than the bulk onset $Ra_c$, all simulations yield $Nu_{bulk}=1$, meaning that heat is transferred conductively in the bulk. For $Ra>Ra_c$, steep scaling relationships can be observed for all cases. The periodic results can be described by $Nu\sim Ra^{3.3\pm 0.2}$, which agree with the steep scaling relationships observed in previous studies with a small Ekman number and using non-slip top and bottom plates and lateral periodic boundary conditions (King et al. 2009, 2012; Stellmach et al. 2014; Cheng et al. 2015) and the polar region of a spherical shell convection system (Gastine & Aurnou 2023). This steep scaling is attributed to the nonlinear effect arising from the Ekman pumping near the non-slip top and bottom plates (Stellmach et al. 2014; Julien et al. 2016; Plumley et al. 2016, 2017). Despite the minor data scattering, all $Nu_{bulk}$ for cylindrical cases have roughly the same scaling relationships $Nu_{bulk}\sim Ra^{3.3\pm 0.1}$, which is consistent with that of the periodic cases. For the cases close to the bulk onset, one can also describe the data using $Nu_{bulk}-1\sim (Ra/Ra_c-1)^{\gamma ^*}$. We find that the exponent $\gamma ^*$ for this scaling relationship is roughly $1.39\pm 0.07$, as we show in figure 3(b). The $Ra_c^*$ in figure 3(b) is the onset of bulk convection determined by extrapolating the best fit power law shown in figure 3(a) to $Nu_{bulk}=1$. Although the magnitude of $\gamma ^*$ is different from the exponent $\gamma$ in the scaling relationship $Nu_{bulk}\sim (Ra/Ra_c)^{\gamma }$, we find that the result for the bulk is consistent with the periodic case, as shown in figure 3(b). This result demonstrates the similar feature in heat transport between the bulk of a system with lateral boundaries and the case with periodic boundary conditions.

Figure 3. (a) Plot of $Nu_{bulk}$ as a function of $Ra/Ra_c$. The grey and black dashed lines respectively represent the power-law fittings of the periodic and all cylindrical results. The vertical black dot-dashed line indicates $Ra_w$. (b) Plot of $Nu_{bulk}-1$ as a function of $Ra/Ra^*_c-1$. The orange dotted curves refer to $Nu_{bulk}-1\sim (Ra/Ra_c^*-1)^{1.39}$. The inset in (b) presents the log–log plot.

Additionally, we find that the onset of bulk convection in the cylindrical cells is slightly larger than the linear stability analysis of the unbounded problem as given by (1.1) ($Ra_c^*\approx 1.40Ra_c$). We also note that although all the cylindrical data are very close to each other, those of $\varGamma =2,Fr=0$ seem to have an onset value closer to $Ra_c$ than the others. This suggests that the difference in $Ra_c$ may be related to the aspect ratio. Moreover, when examining data near $Ra/Ra_c\approx 1$, we observe that the $Nu_{bulk}$ values for these cases are actually greater than 1, implying the presence of convection in the bulk. However, as the $Ra$ dependence for these cases is significantly weaker than the so-called ‘steep scaling’, we deduce that this vertical convective heat transport should arise from plume-like structures emanating from the boundary flow rather than the bulk convective instability. Such flow emission from the wall-localized structures has been reported by Favier & Knobloch (2020) and will be further discussed in § 3.2. As the onset of bulk convection is affected by both bulk convective instability and perturbations from the boundary flow, the influence of aspect ratio to the onset of bulk convection cannot be sufficiently explained by onset theories for either the unbounded system (Niiler & Bisshopp 1965) or the wall mode (Herrmann & Busse 1993; Zhang & Liao 2009). Further theoretical study is required to understand the onset of bulk convection in a rotating system with a lateral boundary. Additionally, the flow strength (quantified by the Reynolds number) for the bulk is also independent of the aspect ratio when bulk convection emerges. Details of the Reynolds number for the bulk can be found in Appendix E.

3.2. Heat transport properties of the boundary flow

Since the bulk of the cylindrical cell reproduces the heat transport properties of the domains with lateral periodic boundary conditions, the difference in global heat transport between these two systems must then be attributed to the boundary flow. Thus, a quantitative investigation of the heat transport in the boundary flow is the key for understanding the global heat transport in the system with lateral boundaries. Analogous to (3.1), the Nusselt number for the boundary flow $Nu_{BF}$ can be defined as

(3.2)\begin{equation} Nu_{BF}\equiv\left\langle \sqrt{RaPr}u_z\theta-\frac{\partial\theta}{\partial z}\right\rangle_{t,z,\phi,r\geqslant r_{BF}}. \end{equation}

We plot $Nu_{BF}$ as a function of $Ra/Ra_w$ in figure 4(a). Compared with $Nu_{bulk}$, as shown in figure 4(a), the behaviour of $Nu_{BF}$ appears to be more complicated. For $Ra/Ra_w\lesssim 3$, $Nu_{BF}$ for all aspect ratios can be described by a linear relationship $Nu_{BF}-1\approx 6(Ra/Ra_w-1)$, which agrees with the reported supercritical behaviour of the global Nusselt number beyond the onset of convection (Zhong et al. 1991; Ecke et al. 1992; Cross & Hohenberg 1993; Zhong, Ecke & Steinberg 1993; Julien et al. 2012; Ecke et al. 2021). In figure 4(b) we plot $Nu_{BF}-1$ as a function of $Ra/Ra_w-1$, which better demonstrates the linear dependence of convective heat flux near the onset of the wall mode.

Figure 4. (a) Plot of $Nu_{BF}$ as a function of $Ra/Ra_w$. The vertical grey dot-dashed line indicates $Ra_c$. The inset of (a) is the zoomed-in view of the results around $Ra/Ra_w=50$. (b) Plot of $Nu_{BF}-1$ as a function of $Ra/Ra_w-1$. The grey dot-dashed curves correspond to $Nu_{BF}-1=6(Ra/Ra_w-1)$. The dashed lines refer to the best fit of $Nu_{BF}=C(Ra/Ra_w)^{\gamma _{BF}}$, and the red dot-dashed line is the guide for the eye indicating $Nu_{BF}\sim (Ra/Ra_w)^1$. The grey dotted curves correspond to the relationship (B3) with higher-order correction (HOC).

As $Ra/Ra_w$ further increases, results for different aspect ratios diverge. When $Ra/Ra_w\gtrsim 10$, $Nu_{BF}$ roughly follows a similar scaling relationship for all cases. The qualitative behaviour is the same whether or not the centrifugal force persists. Comparing the data with the red dot-dashed line in figure 4(a) indicating $Nu_{BF}\sim (Ra/Ra_w)^1$, we find that the observed scaling exponent is close to 1.

Here we propose a theoretical model explaining the leading order of this scaling relationship, following the idea of the marginal boundary layer stability analysis (Malkus 1954; King et al. 2012). For $Ek\ll 10^{-1}$, the second term in (1.2) is negligible, thus, it can be simplified to

(3.3)\begin{equation} Ra_w\approx{\rm \pi}^2\sqrt{6\sqrt{3}}Ek^{{-}1}. \end{equation}

We remark that (3.3) is identical to the onset Rayleigh number of the wall mode in a semi-infinite system with a shear free top and bottom (Herrmann & Busse 1993). Within the thermal boundary layer (with thickness $\delta$), the buoyancy unstablizes the flow, while the Coriolis force acts as a dominating stablizer. The local stability criterion at the edge of the thermal boundary layer is given by $Ra^\delta /Ra_w^\delta \approx 1$, where $Ra^\delta =\alpha g\varDelta _\delta \delta ^3/(\kappa \nu )$, $Ra_w^\delta =Ra_w(Ek^\delta )$, $Ek^\delta =\nu /(2\varOmega \delta ^2)$. Here $\varDelta _\delta$ is the temperature difference across the thermal boundary layer. Substituting $Ek^\delta$ into (3.3) we obtain $Ra_w^\delta \sim A(Ek^\delta )^{-1}$, where $A$ is a prefactor. The marginal stability criterion can then be written as

(3.4)\begin{equation} \frac{Ra^\delta}{Ra_w^\delta}\approx\frac{\alpha g\varDelta_\delta\delta^3Ek^\delta}{A\kappa\nu}=\frac{\alpha g\varDelta_\delta\delta}{2A\varOmega\kappa}\approx1, \end{equation}

which gives

(3.5)\begin{equation} \delta=\frac{2A\varOmega\kappa}{\alpha g\varDelta_\delta}. \end{equation}

Since the heat transport of the boundary flow is dominated by the thermal boundary layer, we can assume that

(3.6)\begin{equation} Nu\approx\frac{\kappa\varDelta_\delta/\delta}{\kappa\varDelta/H}=\frac{\varDelta_\delta}{\varDelta}\frac{H}{\delta}. \end{equation}

Assuming

(3.7)\begin{equation} \varDelta_\delta\sim\varDelta/2, \end{equation}

we can have

(3.8)\begin{equation} Nu\sim H/\delta. \end{equation}

Using this assumption, we obtain the following relationship:

(3.9)\begin{equation} Nu\sim\frac{\alpha g\Delta H}{4A\varOmega\kappa}\sim RaEk\sim(Ra/Ra_w)^1. \end{equation}

The data presented in figure 4(a) exhibit scaling relationships close to (3.9).

We further fit the data using the power-law relationship $Nu_{BF}\approx C(Ra/Ra_w)^\gamma _{BF}$, with $C$ being an aspect-ratio dependent prefactor and $\gamma _{BF}$ the exponent; we list the results in table 1. From table 1 we can see that the obtained scaling exponents are all smaller than unity, suggesting that the above theory may not be able to fully explain the data, and other effects should be introduced. A plausible reason for such discrepancy is the neglect of the higher-order term in the critical Rayleigh number $Ra_w$. The critical Rayleigh number of the wall mode given by Zhang & Liao (2009) consists of two terms, as shown in (1.2). We drop the second term in (1.2) when deriving the scaling relationship of $Nu_{BF}$ (3.9). However, since $Ek^\delta \sim 1/\delta ^2$, for a large $Nu$ number, the thermal boundary layer could be sufficiently thin so that the second term in (1.2) is not negligibly small. By including the second term in (1.2), one can then obtain a relationship for $Nu_{BF}$ with high-order correction,

(3.10)\begin{equation} Nu_{BF}\bigg({\rm \pi}^2\sqrt{6\sqrt{3}}+46.55Ek^{1/3}\left(\frac{\varDelta}{\varDelta^\delta}Nu_{BF}\right)^{2/3}\bigg)\sim RaEk. \end{equation}

Detailed derivations of (3.10) are presented in Appendix B.

Table 1. The results of the best fits of $Nu_{BF}=C(Ra/Ra_w)^{\gamma _{BF}}$.

In figure 4(a) we plot the solutions of (3.10) as grey dashed curves. By introducing the second term of (1.2) in the derivation, (3.10) results in a close agreement with the data, suggesting that neglecting the higher-order term in (1.2) could be one of the reasons leading to the discrepancy between the obtained exponents and (3.9). However, one can see that (3.10) is not a power-law relationship but instead an equation for $Nu_{BF}$, which is mathematically more complicated than (3.9). As (3.9) describes the leading-order behaviour of $Nu_{BF}$, in the latter discussions we use (3.9) instead of (3.10) describing $Nu_{BF}$, but it is important to note that the higher-order correction is required in a strict manner.

3.3. Transition in boundary flow morphology

In figure 4(b), when $Ra$ increases and is beyond the range obeying (3.9), we can see a sudden drop in $Nu_{BF}$ at $Ra/Ra_w\approx 50$. After the sudden drop, $Nu_{BF}$ becomes independent of $\varGamma$ and all data collapse together. To understand this sudden drop in $Nu_{BF}$, we need to examine the flow structure before and after this transition.

We present the temperature fields at the mid-height for different $Ra$ numbers close to the drop of $Nu_{BF}$ in figure 5. Here figure 5(a) corresponds to $Ra=8.71\times 10^8$ (before the sudden drop in $Nu_{BF}$, marked by ‘Before’ in figure 4a) and figure 5(b) $Ra=1.02\times 10^9$ (after the sudden drop, marked by ‘After’ in figure 4a). One can find a coherent boundary flow structure with clear azimuthal periodicity adjacent to the sidewall in the case ‘Before’. The flow in the azimuthal direction is constrained within the boundary flow region ($r>r_{BF}$). One can also observe the emissions of plume-like structures at the contacts of the hot and cold parts. Such plume-like structures move radially to the bulk. As bulk convection is insufficiently strong for mixing, one can still observe an azimuthally periodic pattern away from the boundary flow in the temperature field. When the transition occurs, the boundary flow breaks into fragmented vortices leading to a collapse of the coherent boundary flow state. In this case, the convective mixing is sufficiently strong and the periodic pattern away from the boundary flow does not exist anymore. These plume-like structures emitting from the bulk boundary flow region act as an important source for vortex production, leading to the result that vortex density at the edge of the boundary flow is higher than the uniform distribution. Details of the vortex distribution can be found in Appendix D.

Figure 5. (a,b) Instantaneous temperature field $\theta$ at the mid-height and (c,d) the space–time plot of the sidewall temperature at mid-height for (a,c) $Ra=8.71\times 10^8$ (case ‘Before’) and (b,d) $1.02\times 10^9$ (case ‘After’), respectively. All plots are for $\varGamma =1$, $Ek=1.85\times 10^{-6}$ and $Fr=0$. The grey dashed and black dot-dashed circles in (a,b) respectively correspond to $r_0$ and $r_{BF}$.

We denote the state with coherent boundary flow as the coherent boundary flow state and the state after transition as the vortical boundary flow state. One of the consequences of this transition is a sharp decrease of the size of coherent structures: from the wavelength $\lambda$ of the boundary flow unit to the diameter of vortices. Once the coherent structure is fragmented, heat will be dissipated more easily during the transport and, hence, becomes less efficient in heat transport. Such transition in flow morphology is consistent with the sudden drop in $Nu_{BF}$ observed in figure 4(a). Moreover, as the size of the boundary vortices should be independent of $\varGamma$, $Nu_{BF}$ in the vortical boundary flow state will also be independent of $\varGamma$, which can be found from the data for $Ra$ beyond the transition shown in figure 4.

However, we wish to remark that although after the transition occurs the boundary flow breaks into vortices, one can still observe some heritage features of the coherent boundary flow state, such as a certain integer mode number and retrograde precession of the flow mode. In figure 5(c,d) we present the space–time plot of the sidewall temperature at mid-height for the cases ‘Before’ and ‘After’, respectively. In figure 5(d), even if the temperature signal is highly fluctuating, one can still observe mode number N switching between $N = 3$ and 4. The boundary flow structure precesses in the retrograde direction, which is similar to the coherent state. A similar phenomenon has also been observed in Favier & Knobloch (2020).

This transition also diminishes the difference between the bulk and the boundary flow in flow morphology. We show the streamlines for four cases in figure 6, where (a,b) corresponds to the coherent boundary flow state, and (c,d) the vortical boundary flow state. Comparing figures 6(a) and 6(b), one can see that the boundary flow structure for $Ra=8.71\times 10^8$ is not much different from that for $Ra=8\times 10^7$, even though in the former case bulk convection has already emerged and the latter is below the onset of bulk convection. In figure 6(b) both the coherent boundary flow and bulk convective Taylor columns coexist. These two structures can be easily distinguished from the streamlines. For $Ra=1.02\times 10^9$, the system enters the vortical boundary flow state. As shown in figures 6(c) and 6(d), the boundary flow breaks into columns. In this case, it is difficult to distinguish the bulk and the boundary flows according to flow morphology. However, we remark that even in the vortical boundary flow state, the boundary flow and the bulk flow are still statistically different regions.

Figure 6. Streamlines coloured by the temperature fluctuation $\theta -\langle \theta \rangle _{r,\phi }$ for (a) $Ra=8.0\times 10^7$, (b) $8.71\times 10^8$ (case ‘Before’), (c) $1.02\times 10^9$ (case ‘After’) and (d) $1.30\times 10^9$. All panels are with $\varGamma =1$, $Ek=1.85\times 10^{-6}$ and $Fr=0$.

Such transition in flow morphology also help to understand the influence of centrifugal force to $Nu_{BF}$. For the cases with centrifugal force (set II), one can neither observe a sudden drop nor the $\varGamma$ dependence in $Nu_{BF}$ in figure 4(a). To understand such phenomenon, we present the vertical snapshots of the temperature field along the azimuthal direction at $r/R=0.95$ (within the boundary flow) in figures 7. In figures 7(a) and 7(c) we respectively present the temperature snapshots for the cases ‘Before’ and ‘After’ the transition, clearly depicting the change from the coherent boundary flow state to the vortical one. In figure 7(b) we present the snapshot of the case with parameters similar to the case ‘Before’ ($Ra=8.71\times 10^8, \varGamma =1$), except that it has non-zero centrifugal force $Fr=0.12$. The temperature distribution in figure 7(b) exhibits noticeable differences from figure 7(a) (corresponding to case ‘Before’), but looks very similar to figure 7(c) (corresponding to case ‘After’). Instead of a coherent boundary flow structure, vortical columns are observed in figure 7(b), suggesting that the centrifugal force helps trigger the breakdown of the boundary flow. We also present the temperature field for $Fr=0.12$ and $Ra=1.02\times 10^9$ in figure 7(d), which has a similar $Ra$ to case ‘After’. Similar vortical flow structures can be observed in both figures 7(b) and 7(d), suggesting that the presence of the centrifugal force strongly suppresses the transition in boundary flow morphology. The breakdown of boundary flow coherency induced by the centrifugal force may then lead to a decrease in heat transport efficiency, as shown in figure 4. This could be the reason why no significant drop in $Nu_{BF}$ can be found for the cases with non-zero $Fr$, and their magnitudes of $Nu_{BF}$ are overall smaller than those with zero $Fr$.

Figure 7. Vertical snapshots of the temperature field along the azimuthal direction at $r/R=0.95$. The left panels refer to cases with $Fr=0$ and the right panels with $Fr=0.12$. The Rayleigh number for (a,b) is $Ra=8.71\times 10^8$ and for (c,d) is $1.02\times 10^9$. The aspect ratio $\varGamma$ equals to 1 for all cases. Snapshots (a,c) correspond to the cases ‘Before’ and ‘After’, respectively.

Additionally, we also observe a weak $\varGamma$ dependence in the Reynolds number ($Re$) of the boundary flow before the breakdown of boundary flow coherency. A drop in $Re$ can be found when the transition occurs, although the drop in $Re$ in much less significant than that of $Nu_{BF}$. Details of the Reynolds number for the boundary flow can be found in Appendix E.

3.4. Global heat transport in a cylindrical convection cell

With quantitative understandings of the heat transport in the bulk and the boundary flow, as respectively discussed in §§ 3.1 and 3.2, we are now in a position to explain the global Nusselt number. In this study we have two variables to indicate the flow state: $Ra/Ra_w$ for the properties of the boundary flow and $Ra/Ra_c$ for the bulk. When we discuss the global heat transport, we first need to choose which variable would be a better choice. If we focus on the parameter range in which only geostrophic convection is present, the boundary flow always persists owing to the fact that $Ra_w$ is smaller than $Ra_c$. On the other hand, the bulk region can be either conductive or convective in this parameter range, depending on whether it is beyond the onset of bulk convection. In these two different bulk states, the heat transport of the bulk and, hence, the global one can be significantly different. To better distinguish conductive and convective bulk states, it will be convenient to use $Ra/Ra_c$ as the variable of choice for examining the global heat transport efficiency. Additionally, as discussed in 3.1, the scaling exponent $\gamma$ for $Nu_{bulk}\sim (Ra/Ra_c)^\gamma$ and $\gamma ^*$ for $Nu_{bulk}-1\sim (Ra/Ra_c-1)^{\gamma ^*}$ are in fact different when $Ra$ is close to onset of bulk convection. Such a difference may introduce inaccuracy when reconstructing the global $Nu$. Nevertheless, for $Ra/Ra_c\approx 1$, the contributions from the boundary flow are significant and the difference between the descriptions of $Nu_{bulk}$ may have a minor influence. For reason of simplicity, we hereby use $Nu_{bulk}\sim (Ra/Ra_c)^\gamma$ to describe the heat transport efficiency for the bulk. As for the boundary flow, (3.9) only captures the leading-order effect and (3.10) obviously provides a more precise description. Also for the purpose of a concise final result that could capture the dominating heat transport properties of the system, we use (3.9) for $Nu_{BF}$ in the following discussion. We also wish to remark that one can obtain a more precise description if more accurate (but probably more complicated) relationships for the $Nu_{bulk}$ and $Nu_{BF}$ are applied.

The global Nusselt number can be viewed as a superposition of the contributions from the bulk and the boundary flow. When the width of the boundary flow $\delta _{BF}$ is much smaller than the cell radius $R$, we can assume that $r_{BF}/R\approx 1$. In this case, $Nu$ can be written as

(3.11)\begin{equation} Nu=\frac{r_{BF}^2}{R^2}Nu_{bulk}+\frac{R^2-r_{BF}^2}{R^2}Nu_{BF}\approx Nu_{bulk}+\frac{4\delta_{BF}/H}{\varGamma}Nu_{BF}. \end{equation}

For fixed $Ek$ and $Ra\gtrsim Ra_c$, the heat transport efficiency of the boundary flow follows the power-law relation $Nu_{BF}\sim (Ra/Ra_w)^1$. We define $Nu_{BF}^c\equiv Nu_{BF}|_{Ra/Ra_c=1}$. We note that $Ek$ is fixed for sets I and II; thus, the local Nusselt number of the boundary flow can be written as

(3.12)\begin{equation} Nu_{BF}=Nu_{BF}^c (Ra/Ra_c)^1. \end{equation}

On the other hand, $Nu_{bulk}$ can be described by the steep scaling $Nu_{bulk}\approx Ra/Ra_c^\gamma$. Thus, (3.11) is written as

(3.13)\begin{equation} Nu\approx (Ra/Ra_c)^\gamma+B (Ra/Ra_c)^1\ (Ra\gtrsim Ra_c), \end{equation}

where the prefactor $B$ is defined as

(3.14)\begin{equation} B\equiv\frac{4\delta_{BF}/H}{\varGamma}Nu_{BF}^c. \end{equation}

For the coherent boundary flow state ($Ra\leqslant Ra_t$), when the data obey (3.9), the width of boundary flow $\delta _{BF}$ is independent of $Ra$. Thus, for a given $\varGamma$ and $Ek$, the prefactor $B$ is roughly a constant. At $Ra=Ra_c$, one can obtain $B\approx Nu|_{Ra/Ra_c=1}-1$, wherein $B$ can be determined using $Nu$ at $Ra=Ra_c$. However, the existence of sidewall may slightly alter the onset of bulk flow, as shown in figure 3(a). Thus, the exact onset Rayleigh number of the bulk may need to be determined in experiments. We note that, for $Ra\lesssim Ra_c$, the $Ra$ dependence of $Nu$ is given by

(3.15)\begin{equation} Nu\approx1+B(Ra/Ra_c)^1\ (Ra\lesssim Ra_c), \end{equation}

which is noticeably different from (3.13). Thus, the exact onset Rayleigh number of the bulk convection and $Nu|_{Ra=Ra_c}$ can be determined according to transition from (3.15) to (3.13). We plot $Nu$ vs $Ra$ in figure 8, and also (3.13) for $Ra>Ra_c^*$ and (3.15) for $Ra< Ra_c^*$ as dashed curves in the same figure. Here $Ra_c^*$ is the onset Rayleigh number determined according to the transition from (3.15) to (3.13). We remark that (3.15) is valid only when $Nu_{BF}$ follows (3.9). Due to the existence of the transitional regime, the Rayleigh number that (3.9) becomes valid is close to $Ra_c$ for sets I and II, as shown in figure 4(a). Thus, not all data for $Ra< Ra_c$ obey (3.15). Nevertheless, for cases beyond and below the onset of bulk convection, their $Ra$ dependences of $Nu$ are significantly different, so that $Nu|_{Ra=Ra_c}$ can still be determined. From figure 8, one can see that (3.13) can properly describe the obtained $Nu$ in the geostrophic regime. The red vertical dotted line corresponds to the transition in flow morphology in the boundary flow ($Ra=Ra_t$). When the transition occurs, the scaling relation for the boundary flow (3.9) becomes invalid and $Nu$ gradually deviates from (3.13). As $Ra$ further increases, buoyancy gradually becomes dominant and $Nu$ asymptotically approaches the results of a non-rotating RBC.

Figure 8. Plot of $Nu$ as a function of $Ra$ for sets I and II. The grey upper triangles are data collected from Funfschilling et al. (2005). The dashed curves refer to (3.13) for $Ra>Ra_c^*$ and (3.15) for $Ra< Ra_c^*$. The red dotted line corresponds to $Ra=Ra_t$.

4. Conclusion

In this study we use a DNS approach to investigate the geostrophic rotating RBC in cylindrical cells. We explore the rotating RBC system in cylindrical cells with different aspect ratios and over a wide range of $Ra$ from the onset of boundary flow to the rotation-unaffected regime. Our simulation yielded three data sets. Set I has $Ek$ fixed and $Ra$ varying and in set III we have $Ra$ fixed and $Ek$ varying. We compare these two data sets in order to reveal the $Ra$ and $Ek$ dependences of the heat transport and flow morphology of this system. As centrifugal force is inevitable in experiments, to ensure that our findings can be applied to the practical situations, set II is produced to study the effects of the centrifugal force. To understand the heat transport of this system with lateral boundaries, we decompose the bulk and the boundary flow and investigate these two domains separately.

We first focus on sets I and II with $Ek$ fixed. For the bulk region, the $Ra$ dependence of $Nu_{bulk}$ is consistent with the numerical simulations using periodic boundary conditions. As for the boundary flow, after an intermediate regime, we observed a similar scaling relationship close to $Nu_{BF}\sim (Ra/Ra_w)^1$ for all cases. For set I, such a scaling relationship can be observed at $Ra/Ra_w\gtrsim 10$ ($Ra\gtrsim \times 10^8$), and deviates from it at $Ra/Ra_w\approx 60$ ($Ra\approx 1\times 10^9$). For set II, data for $Ra/Ra_w\lesssim 60$ ($Ra\lesssim 1\times 10^9$) obey the above scaling relationship. By introducing the idea of boundary layer marginal stability theory, one can obtain the leading-order scaling relationship $Nu_{BF}\sim (Ra/Ra_w)^1$ if only the dominating term for $Ra_w$ is used, and a more precise description of the data can be obtained if the higher-order term in $Ra_w$ is considered.

As $Ra$ increases, we observe a sharp transition in boundary flow morphology from the coherent state to the vortical state, which is accompanied by a sudden drop in $Nu_{BF}$. When the transition occurs, the scaling relationship $Nu_{BF}\sim (Ra/Ra_w)^1$ breaks down. The transitional Rayleigh number $Ra_t$ for $Ek=1.85\times 10^{-6}$ is found approximately equaling $8.71\times 10^8$ in our study. After the transition, the boundary flow breaks into vortices, manifesting the breakdown of flow coherency. Consequently, a drop in $Nu_{BF}$ can be observed as a result of the transition in flow morphology that the boundary flow becomes fragmented and heat exchange with the surrounding becomes more significant. After the transition, it becomes difficult to identify the bulk and the boundary flows based on flow morphology alone. Although this transition is difficult to discern from the global $Nu$, we remark that such a transition is well defined according to the sharp change in both $Nu_{BF}$ and the flow morphology. Additionally, comparing the temperature field for sets I and II, we find that the centrifugal force helps trigger the breakdown of the boundary flow and results in a weaker flow coherency. This probably is the reason why one cannot observe such a transition for the cases with non-zero $Fr$, and why $Nu_{BF}$ for the cases with zero $Fr$ is larger than those with non-zero ones. However, the current study is not sufficient to fully understand the influence of centrifugal force on the boundary flow. This would be an interesting topic for future studies.

Finally, we then revisit the global $Nu$ and try to provide a unifying understanding. For the geostrophic rotating RBC, the global $Nu$ can be described by (3.13). When centrifugal force exists but does not significantly alter flow structures, its influence is reflected on the magnitude of the prefactor $B$ only and will not affect the scaling exponent. We further show that the prefactor $B$ can be determined using the global $Nu$. Thus, (3.13) can still be used to describe $Nu$ for experiments whether or not the centrifugal force is present. For set III with fixed $Ra$ and $Ek$ varying, the main discoveries are analogous to those for sets I and II, but with some variations in the width of the boundary flow and also the global $Nu$. Details of the discussions can be found in Appendix F.