2.1. Governing equations

We consider Boussinesq convection of homogeneous fluid in a spherical shell of inner radius $r_i$ and outer radius $r_o$ that uniformly rotates at $\boldsymbol \varOmega =\varOmega \boldsymbol {\hat z}$. Convection is driven by a fixed temperature difference $\Delta T=T_{i}-T_{o}$ between the inner and outer boundaries, under the gravity $\boldsymbol {g}=-g_o \boldsymbol {r}/r_o$. Using the shell thickness $D = r_{o} - r_{i}$ as the length scale, $\varOmega ^{-1}$ as the time scale, $\Delta T$ as the temperature scale and gravity at the outer boundary $g_o$ as the reference value, the dimensionless governing equations can be expressed as

(2.1)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t}+ \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}+2 \widehat{\boldsymbol{z}} \times \boldsymbol{u} =-\boldsymbol{\nabla} P+E \boldsymbol{\nabla}^{2} \boldsymbol{u}+ \frac {RaE^2}{Pr} T \boldsymbol{r}, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial T}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T =\frac{E}{Pr} \boldsymbol{\nabla} ^{2} T, \end{gather}$$

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

where $\boldsymbol {u}$ is the velocity, $P$ is the reduced pressure and $T$ is the temperature. The system is defined by three dimensionless control parameters, the Ekman number $E$, the Rayleigh number $Ra$ and the Prandtl number $Pr$:

(2.4a–c)\begin{equation} E=\frac{\nu}{\varOmega D^2},\quad Ra=\frac{\alpha g_{0}\Delta T D^3}{\nu \kappa}, \quad Pr=\frac{\nu}{\kappa}. \end{equation}

Here $\nu$ is the kinematic viscosity, $\kappa$ is the thermal diffusivity and $\alpha$ is the thermal expansion coefficient. The rotational modified Rayleigh number $Ra^*$ is often used in the rotating convection and $Ra^*$ is also the squared convective Rossby number, which is used in many rotating convection studies (e.g. Gilman Reference Gilman1977; Gastine et al. Reference Gastine, Yadav, Morin, Reiners and Wicht2013). We have

(2.5)\begin{equation} Ra^*=\frac{\alpha_{T} g_{0}\Delta T}{\varOmega^2 D}=\frac{RaE^2}{Pr}={Ro_c}^2, \end{equation}

which provides a clear relationship between thermal buoyancy and the Coriolis force in rotating convection.

In this study we fix the radius ratio of the shell $\eta =r_i/r_o=0.35$. Both boundaries are impermeable, no-slip and held at constant temperatures.

2.2. Numerical technique

We use the open-source code XSHELLS (https://www.bitbucket.org/nschaeff/xshells/) to solve the governing equations (2.1)–(2.3) subjected to the boundary conditions. The fluid is assumed to be incompressible, and the velocity $\boldsymbol u$ can be decomposed into toroidal and poloidal components:

(2.6)\begin{equation} \boldsymbol u=\boldsymbol{\nabla} \times(\mathcal{T} \boldsymbol r)+\boldsymbol{\nabla} \times \boldsymbol{\nabla} \times(\mathcal{P} \boldsymbol r). \end{equation}

The toroidal $\mathcal {T}$ and poloidal $\mathcal {P}$ scalar fields and the temperature field $T$ are expanded in terms of spherical harmonic expansion on spherical surfaces. The code XSHELLS uses a second-order finite differences method in the radial direction and pseudo-spectral spherical harmonic expansion. The spectral expansion is truncated up to spherical harmonics of degree $L_{max}$ and $Nr$ denotes the number of radial gird points. The time-stepping scheme is second order, and treats the diffusive terms implicitly, while the nonlinear and Coriolis terms are handled explicitly. This code also uses the SHTns library to speed up the spherical harmonic transformations (Schaeffer Reference Schaeffer2013).

2.3. Diagnostics

We analyse several diagnostics properties to quantify the influence of different control parameters on heat and momentum transports. We adopt several notations regarding averaging procedures. Overbars ${\overline {\cdots }}$ correspond to temporal averaging, angular brackets $\langle \cdots \rangle$ to spatial averaging over the entire spherical shell volume and $\langle \cdots \rangle _s$ to an average over a spherical surface:

(2.7a–c)\begin{align} \bar{f}=\frac{1}{\tau} \int_{t_{0}}^{t_{0}+\tau}f\, \mathrm{d}t, \quad \langle \,f\rangle=\frac{1}{V} \int_{V} f(r, \theta, \phi) \,\mathrm{d} V, \quad \langle \,f\rangle_s=\frac{1}{4 {\rm \pi}} \int_{0}^{2 {\rm \pi}} \int_{0}^{\rm \pi} f(r, \theta, \phi) \sin \theta\, \mathrm{d} \theta \,\mathrm{d} \phi, \end{align}

where $\tau$ is the time averaging interval, $V$ is the volume of the spherical shell, $r$ is the radius, $\theta$ is the colatitude and $\phi$ is the longitude.

The Nusselt number $Nu$ denotes heat transport, the ratio of the total heat flux to the conduction heat flux. In spherical shells the conductive temperature profile $T_{c}$ is the solution of

(2.8a–c)\begin{equation} \frac{\mathrm{d} T_{c}}{\mathrm{d} r} =-\frac{r_{i}r_{o}}{r^2}, \quad T_{c}(r_{i}) = 1, \quad T_{c}(r_{o}) = 0. \end{equation}

Following Gastine, Wicht & Aurnou (Reference Gastine, Wicht and Aurnou2015) yields

(2.9)\begin{equation} T_{c}(r)=\frac{\eta}{(1-\eta)^2}\frac{1}{r}-\frac{\eta}{1-\eta},\end{equation}

where $\eta =r_{i}/r_{o}$ is the radius ratio. The notation $\vartheta$ is introduced to define the time and horizontally averaged radial dimensionless temperature profile

(2.10)\begin{equation} \vartheta(r)=\overline{\langle T\rangle_s}.\end{equation}

Then the Nusselt number is

(2.11)\begin{equation} Nu=-\eta\frac{\mathrm{d}\vartheta}{\mathrm{d}r}(r=r_{i}) =-\frac{1}{\eta}\frac{\mathrm{d}\vartheta}{\mathrm{d}r}(r=r_{o}).\end{equation}

Here we only consider the global averaged $Nu$, but the efficiency of convective heat transfer slightly varies as latitude in rotating convection (Wang et al. Reference Wang, Santelli, Lohse, Verzicco and Stevens2021; Gastine & Aurnou Reference Gastine and Aurnou2023). The total kinetic energy is given by

(2.12)\begin{equation} E_{k}=\frac{1}{2} \int_{v} u^2 \,\mathrm{d}V=\sum_{l = 1}^{l_{max}} \sum_{m = 0}^{l} \mathcal{E}_{l}^{m}(t), \end{equation}

where $\mathcal {E}_{l}^{m}$ is the dimensionless kinetic energy at a spherical harmonic degree $l$ and order $m$. The total kinetic energy can be decomposed into zonal and non-zonal parts:

(2.13a,b)\begin{equation} E_{zon}=\frac{1}{2} \int_{v}{(u_{\phi}^{0})}^2 \mathrm{d}V, \quad E_{non}=E_{k}-E_{zon}.\end{equation}

We define the Rossby number as

(2.14)\begin{equation} Ro=\frac{U_{rms}}{\varOmega D},\end{equation}

where $U_{rms}$ is the dimensional root-mean-square velocity. Based on the non-dimensionalisation we used, the Rossby number can be determined through kinetic energy as

(2.15)\begin{equation} Ro=\sqrt{\frac{2E_{k}}{V}}.\end{equation}

Accordingly, we have the non-zonal Rossby number

(2.16)\begin{equation} Ro_{non}=\sqrt{\frac{2E_{non}}{V}},\end{equation}

and the zonal Rossby number is

(2.17)\begin{equation} Ro_{zon}=\sqrt{\frac{2E_{zon}}{V}}.\end{equation}

We also define the local Rossby number $Ro_{\ell }$ as

(2.18)\begin{equation} Ro_{\ell}=\frac{U_{rms}}{\varOmega \ell}=Ro\frac{D}{\ell},\end{equation}

where $\ell$ is the typical flow length scale and is determined from the time-averaged kinetic energy spectrum following Christensen & Aubert (Reference Christensen and Aubert2006),

(2.19)\begin{equation} \ell^{-1} = \overline{\left(\frac{\displaystyle D \sum_{l = 1}^{l_{max }} \sum_{m = 0}^{l} l \mathcal{E}_{l}^{m}(t)}{\displaystyle {\rm \pi}\sum_{l = 1}^{l_{max }} \sum_{m = 0}^{l} \mathcal{E}_{l}^{m}(t)}\right)}. \end{equation}

In rotating spherical shell convection the zonal flow component does not contribute to the heat transport from the inner to the outer shell, so we focus on the non-zonal component and it can be determined through the non-zonal Rossby number as

(2.20)\begin{equation} Re_{non}=\frac{Ro_{non}}{E},\end{equation}

and the zonal Reynolds number is

(2.21)\begin{equation} Re_{zon}=\frac{Ro_{zon}}{E}.\end{equation}

3.1. Overview

We have simulated a total of 211 cases in the parameter space of $1.79\times 10^{5} \le Ra \le 5\times 10^{9}$, $1\times 10^{-6} \le E \le 1\times 10^{-4}$ and $0.01 \le Pr \le 100$. Details of input and diagnostic parameters for all models are listed in table 1 in Appendix A. Figure 1 shows $Nu-1$ as a function of $Ra$ for all the numerical models. We use different colours and shapes to distinguish different values of $Pr$ and $E$, respectively. Broadly speaking, the ($Nu-1, Ra$) diagram shows that $Nu$ is insensitive to $Pr$ when $Pr \ge 1$ (black, blue and green symbols) but highly dependent on $Pr$ when $Pr < 1$ (orange and red symbols). This is in line with previous studies on the $Pr$ effects in Rayleigh–Bénard convection (Verzicco & Camussi Reference Verzicco and Camussi1999), and suggests that we need to consider different scaling behaviours for small $Pr$ and large $Pr$.

Figure 1. Regime diagram of all numerical models in the parameter space $(Nu-1, Ra)$. Different colours and shapes correspond to different $Pr$ and $E$. Open symbols denote the onset (ON) regime; half-filled symbols denote the multiple modes (MM) regime; filled symbols denote the geostrophic turbulence (GT) regime; central dot symbols denote the weakly rotating (WR) regime. The dashed line corresponds to the non-rotating scaling $(Nu-1)\propto Ra^{1/3}$.

We also see from figure 1 that slopes of $Nu-1$ change as $Ra$ increases because of different convective regimes. Our numerical models can be separated into four different regimes depending on the flow morphology and heat transport efficiency. When the Rayleigh number is slightly above the critical value $Ra_c$, convection is characterised by a single ON mode, which is referred to as the ON regime (open symbols in figure 1). The critical $Ra_c$ and the structure of the ON mode depend on $E$ and $Pr$ (see table 2 in Appendix B). In the ON regime, weak convection has little contribution to the heat transfer, so we can see that $(Nu-1)\ll 1$ but increases steeply as a function of $Ra$. As we increase $Ra$, the convective flow remains laminar, but several convective modes can coexist, which is referred to as the multiple modes (MM) regime (half-filled symbols in figure 1). In the MM regime the Nusselt number remains small, i.e. $(Nu-1)<1$, meaning that heat is mainly transferred by conduction. Further increasing $Ra$, convection becomes more complex and even turbulent but exhibits columnar structures along the rotation axis, which is referred to as the GT regime (filled symbols in figure 1). As the transition from laminar to turbulent flow is not always well defined, it is not straightforward to define the boundary between the MM and GT regimes. There exists sudden jumps of $Nu-1$ as increasing $Ra$ for $Pr<1$ cases (see red and orange symbols in figure 1), so we define such jumps as the boundary between the MM and GT regimes when $Pr<1$. However, $Nu-1$ smoothly increases as $Ra$ increases when $Pr\ge 1$, so we simply define $Nu>2$ as the criterion to enter the GT regime for $Pr\ge 1$ cases following Gastine et al. (Reference Gastine, Wicht and Aubert2016). This criterion amounts to saying that the heat transported by convection overtakes the conductive heat transfer. When $Ra$ is sufficiently large for fixed $E$ and $Pr$, convection becomes less geostrophic by breaking the rotational constraint and eventually approaches non-rotating convection. Again it is not straightforward to define the transition from rotating to non-rotating turbulence. Some scalings and criteria have been proposed to characterise the transition from rotating to non-rotating convection (Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012; Gastine et al. Reference Gastine, Wicht and Aubert2016; Long et al. Reference Long, Mound, Davies and Tobias2020), but it is difficult to reconcile numerical simulations with different $Pr$ as we shall show. Here we simply use the local Rossby number $Ro_\ell >0.1$ as the criterion for the weakly rotating (WR) regime (central dot symbols in figure 1). The local Rossby number is the ratio of advection to Coriolis force and $Ro_\ell \lesssim 0.1$ is often seen as an indicator of rotating flows (e.g. Davidson Reference Davidson2014). This criterion is a bit subjective but we do see that the columnar structures tend to break when $Ro_\ell > 0.1$. We define different flow regimes to facilitate our following discussions, but determining the regime boundaries is not the focus of this study.

Figures 2 and 3 show typical flow structures of four different regimes for cases with $Pr=7$ and $Pr=0.01$, respectively. We have mentioned that the heat transfer exhibits different behaviours between models with $Pr\ge 1$ and $Pr< 1$. Such differences are also reflected in the flow structures in all regimes. In the ON regime the convective mode at $Pr=7$ takes the form of columnar rolls in the vicinity of the tangent cylinder (Dormy et al. Reference Dormy, Soward, Jones, Jault and Cardin2004; Barik et al. Reference Barik, Triana, Calkins, Stanley and Aurnou2023), while the convection at $Pr=0.01$ shows spiralling columnar structures that occur in the whole domain outside the tangent cylinder (Zhang Reference Zhang1992). The spiralling structure corresponds to the transitional onset between the viscous convection and inertial convection (Zhang & Liao Reference Zhang and Liao2017). In the MM regime the convective flows are similar to the ON modes but apparently show other unstable modes coexisting. We will show in § 3.3 that the interaction of MM may take place in triadic resonances (Lin Reference Lin2021). In the GT regime, convection flows are chaotic in the equatorial plane but are organized along the rotation axis. The horizontal length scale is dramatically different between cases with $Pr=7$ and $Pr=0.01$ in the same $E$. In the case of $E=1\times 10^{-5}$ and $Pr=7$, the range of $\ell$ is $0.064\unicode{x2013}0.089$. However, when $Pr=0.01$, the range of $\ell$ is $0.488\unicode{x2013}0.529$. In the WR regime, columnar convective structures are broken as the rotational constraint becomes less important. Figures 2 and 3 provide an overview of the flow morphology in different regimes. In the following we focus on the scaling behaviours of the heat transfer and typical convective flow speeds in terms of control parameters $Ra$, $E$ and $Pr$.

Figure 2. Contours of velocity magnitude in units of the Rossby number at meridional and equatorial planes, and spherical surfaces showing typical flow structures of four different regimes with $Pr=7$. The inner (outer) surface corresponds to a spherical surface of radius $r_i+0.1D$ ($r_o-0.1D$). The colour bar represents Rossby number $Ro$. (a) The ON regime; (b) MM regime; (c) GT regime; (d) WR regime.

Figure 3. Same as figure 2 but for cases with $Pr=0.01$. (a) The ON regime; (b) MM regime; (c) GT regime; (d) WR regime.

3.2. Onset regime

We have shown the flow structures of convection ON modes are different depending on $Pr$. In this subsection we show scaling behaviours of the heat transfer and convective flow speed near the onset, with particular attention to $Pr$ dependence. The critical Rayleigh number $Ra_{c}$ and azimuthal wavenumber $m_{c}$ of the ON mode are given in table 2 for different $E$ and $Pr$. It has been show that the convective heat transport increases linearly with $Ra/Ra_c$ near the onset (Busse & Or Reference Busse and Or1986; Gillet & Jones Reference Gillet and Jones2006),

(3.1)\begin{equation} (Nu-1) \propto \left(\frac{Ra}{Ra_{c}}-1\right).\end{equation}

Figure 4(a) shows $Nu-1$ as a function of $Ra/Ra_{c}-1$ for the cases in the ON regime. We see that $Nu-1$ is linearly proportional to $Ra/Ra_c-1$ as predicted, but the prefactor of the scaling (3.1) depends on $Pr$, which reflects the $Pr$ dependence of the convection onset. For numerical models with $Pr \ge 1$, the prefactor is similar for different $Pr$ but weakly depends on $E$, corresponding to the viscous convection mode. At low $Pr<1$, the ON modes become spiralling columnar structures that fill the whole domain outside the tangent cylinder (figure 3a), corresponding to the transitional mode between viscous convection and inertial convection (Zhang Reference Zhang1992). As the ON mode and critical $Ra_c$ depends on $Pr$ in the transitional convection, the prefactor in the scaling (3.1) also highly depends on $Pr$ when $Pr<1$. Our numerical models do not reach sufficiently low $Pr$ to show the purely inertial convection in which the onset of convection is in the form of inertial modes (Zhang Reference Zhang1994).

Figure 4. (a) Nusselt number $Nu-1$ as a function of $Ra/Ra_{c}-1$ in the ON regime; (b) $Re_{non}$ as a function of $Ra_{Q}Pr^{-2}E^{2/3}$ in the ON regime. Different colours and shapes represent different values of $Pr$ and $E$.

For the convective flow speed $Re_{non}$ near the onset, a scaling based on the visco-Archimedean–Coriolis (VAC) force balance has been proposed (Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001; King & Buffett Reference King and Buffett2013; King, Stellmach & Buffett Reference King, Stellmach and Buffett2013):

(3.2)\begin{equation} Re_{non}\sim Ra_{Q}^{1/2}Pr^{-1}E^{1/3}.\end{equation}

Here $Ra_Q=(Nu-1)Ra$ is the flux-based Rayleigh number. The above VAC scaling can also be derived from the balance between the viscous dissipation rate and work done by the buoyancy force (Gastine et al. Reference Gastine, Wicht and Aubert2016). Figure 4(b) shows $Re_{non}$ vs $Ra_{Q}Pr^{-2}E^{2/3}$ for all cases near the onset. We can see that numerical results with different $Pr$ fit well a unified VAC scaling (3.2). We note that the VAC scaling is given in terms of $Ra_Q$ that implicitly carries the $Pr$ dependence of $Nu$. The VAC scaling has been confirmed by previous numerical simulations at $Pr$ around unity (Gastine et al. Reference Gastine, Wicht and Aubert2016; Long et al. Reference Long, Mound, Davies and Tobias2020). At low $Pr$ in our numerical model, convection near the onset presumably should be the transitional convection (Zhang Reference Zhang1992), in which the effect of inertial force is gradually becoming apparent. However, the VAC scaling tends to be still valid for our numerical models at low $Pr$. This means that the $Pr$ we calculated is still not small enough to reach the inertial ON regime where the inertial force plays a dominant role.

In summary for the ON regime, both the heat transfer and typical convective flow speed can be explained by previous theoretical predictions. The prefactor of the heat transfer scaling depends on $Pr$, which reflects the $Pr$ dependence of the convection ON mode.

3.3. Multiple modes regime

As we increase $Ra$, the convective flow remains laminar but exhibits several mode interactions (figures 2b and 3b). We refer to such flow patterns as the MM regime, which corresponds to a transition from simple convection ON modes to more complex convection modes. As this regime exists only in a narrow parameter range (see figure 1), we do not analyse systematic scaling behaviours but focus on the characteristics of MM interaction in the MM regime.

Previous numerical simulations (Horn & Schmid Reference Horn and Schmid2017; Lam, Kong & Zhang Reference Lam, Kong and Zhang2018) and laboratory experiments (Aurnou et al. Reference Aurnou, Bertin, Grannan, Horn and Vogt2018) using liquid gallium ($Pr \approx 0.025$) observed multiple mode interactions in rotating convection. More recently, Lin (Reference Lin2021) revealed that the MM interaction take place in the form of triadic resonances in spherical rotating convection at low $Pr\le 0.01$. Triadic resonance is a generic instability mechanism in rotating fluids (Kerswell Reference Kerswell2002; Le Bars, Cébron & Le Gal Reference Le Bars, Cébron and Le Gal2015), in which a primary inertial mode with azimuthal wavenumber $m_0$ and frequency $\omega _0$ can excite a pair of unstable inertial modes with wavenumbers $m_1, m_2$ and frequencies $\omega _1$, $\omega _2$ matching the resonance conditions:

(3.3a,b)\begin{equation} \omega_0=\omega_1 \pm \omega_2, \quad m_0=m_1 \pm m_2.\end{equation}

In this study we find that triadic resonances can take place at both small and moderate $Pr$. Figure 5 shows time-averaged energy spectra (a,c) and time evolution of the kinetic energy contained in the four dominant components (b,d) for two cases with $Pr=0.1$ (a,b) and $Pr=7$ (c,d). For the case with $Pr=0.1$, the four highest peaks in the $m$ spectrum correspond to $m=11, 0, 9, 2$, respectively. The $m=11$ component represents the primary convective mode, while the $m=0$ component can be attributed to the mean flow generated by the nonlinear interaction of the primary mode (Zhang & Liao Reference Zhang and Liao2017). As the energy in the primary mode saturated, two other components with $m=2$ and $m=9$ that satisfy the triadic resonance condition start to exponentially grow and then saturate. This is a typical behaviour of the triadic resonance (Lin Reference Lin2021) and suggests that two secondary modes are excited through the triadic resonance with the primary mode. The case with $Pr=7$ in the bottom panel exhibits similar behaviours but with the primary mode of $m=12$ and two secondary modes of $m=4$ and $m=8$. We should note that triadic resonances at moderate $Pr$ should be due to nonlinear interactions of three thermal Rossby waves rather than purely inertial waves at very low $Pr$ rotating convection (Lin Reference Lin2021) or in the mechanical driven rotating flows (Le Bars et al. Reference Le Bars, Cébron and Le Gal2015).

Figure 5. Energy spectra of two cases in the MM regime with possible triadic resonances. (a,b) The case with $E=1\times 10^{-5}$, $Pr=0.1$ and $Ra=5.1\times 10^{6}$; (c,d) the case with $E=3\times 10^{-5}$, $Pr=7$ and $Ra=4.67\times 10^{6}$. (a,c) Time-averaged kinetic energy spectrum as a function of the azimuthal wavenumber $m$. (b,d) Time evolution of the kinetic energy contained in different $m$ components indicated by circles in the left panel.

Our numerical results, in line with the early work of Lin (Reference Lin2021), suggest that the triadic resonance may provide a generic mechanism of the transition from the single ON mode to MM coexisting in rotating convection with both small and moderate $Pr$. Further increasing $Ra$ would excite more and more unstable modes and eventually lead to turbulent convection that will be discussed in the following subsection.

3.4. Geostrophic turbulence regime

Geostrophic turbulence is thought to be the most relevant regime to convection in rapidly rotating stars and planets. Both the flow morphology and scaling behaviours of heat transfer and momentum transport have been extensively discussed in previous studies (e.g. Aurnou et al. Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015). In this section we focus on the scaling behaviours of heat transfer $Nu$ and the convective flow speed $Re_{non}$ at different $Pr$.

In the GT regime the well-known diffusion-free scalings have been widely discussed and compared with numerical simulations and laboratory experiments (Cheng et al. Reference Cheng, Aurnou, Julien and Kunnen2018; Hawkins et al. Reference Hawkins, Cheng, Abbate, Pilegard, Stellmach, Julien and Aurnou2023). The diffusion-free scalings (also called inertial scalings or CIA scalings) were derived using different approaches and arguments in the literature (e.g. Stevenson Reference Stevenson1979; Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001; Gillet & Jones Reference Gillet and Jones2006; Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012; Barker et al. Reference Barker, Dempsey and Lithwick2014; Aurnou et al. Reference Aurnou, Horn and Julien2020), which are essentially based on the CIA force balance and the mixing length theory, leading to scalings independent of the fluid viscosity and thermal diffusivity. These scalings predict that the efficiency of convective heat transfer follows

(3.4)\begin{equation} (Nu-1)\sim Ra^{3/2}E^{2}Pr^{-1/2},\end{equation}

and the convective flow speeds follow

(3.5)\begin{equation} Re_{non}\sim (Ra_{Q}Pr^{-2}E^{1/2})^{2/5},\end{equation}

where $Ra_Q=(Nu-1)Ra$ is the flux-based Rayleigh number. Previous studies found that the scalings ((3.4) and (3.5)) can fit numerical and experimental data in a certain parameter regime, but most numerical simulations set $Pr$ around unity and laboratory experiments usually use water ($Pr\approx 7$) as the working fluid (but see other rotating convection experiments with different $Pr$, e.g. Aurnou et al. Reference Aurnou, Bertin, Grannan, Horn and Vogt2018; Abbate & Aurnou Reference Abbate and Aurnou2023). Systematic tests of the diffusion-free scalings over a wide range of $Pr$ seem to still be lacking. Here we compare the diffusion-free scalings with our numerical models over a wide range of $Pr$ ($10^{-2}\le Pr \le 10^2$).

Figure 6 shows $Nu-1$ as a function of the control parameters $Ra$, $E$ and $Pr$ in the GT regime. As already noted from figure 1, a single scaling is not able to reconcile numerical results with different $Pr$ because $Nu-1$ highly depends on $Pr$ when $Pr \le 1$ and becomes nearly independent of $Pr$ when $Pr \ge 1$. Therefore, in figure 6 we plot $Nu-1$ separately for cases with $Pr \le 1$ and $Pr \ge 1$. We can see from figure 6(a) that the diffusion-free scaling (3.4) matches numerical data reasonably well with $Pr \le 1$ and $(Nu-1)\ge 1$, except for some data points in the left-bottom corner that correspond to cases with $(Nu-1)<1$. In these cases, the convective flow becomes chaotic but $Nu$ remains small because of efficient thermal conduction at small $Pr$. The diffusion-free scaling is expected to be valid in the asymptotic regime of $Nu\gg 1$. However, it is a huge numerical challenge to achieve large $Nu$ but remain in the GT regime with small $Pr$, as this would require reducing $E$ and increasing $Ra$ meanwhile. So we find that the diffusion-free scaling (3.4) for the convective heat transfer tends to be valid when $Pr \le 1$ and $(Nu-1)\ge 1$ based on our numerically accessible models.

Figure 6. Nusselt number versus control parameters in the GT regime: (a) $Nu-1$ as a function of $Ra^{3/2} E^2Pr^{-1/2}$ for cases with $Pr\le 1$, (b) $Nu-1$ as a function of $Ra^{3/2} E^2$ for cases with $Pr\ge 1$. Different shapes and colours represent different $E$ and $Pr$ as in figure 1.

As $Nu$ becomes nearly independent of $Pr$ when $Pr\ge 1$ (in figure 1), figure 6(b) plots $Nu-1$ as a function of $Ra^{3/2}E^2$ without $Pr$ dependence for cases with $Pr\ge 1$ in the GT regime. Here we basically drop the $Pr$ dependency and keep the same power law for $Ra$ and $E$ as in the diffusion-free scaling (3.4) for reference. We can see that $Nu-1$ does not show obvious dependence on $Pr$ for fixed $Ra$ and $E$. The almost independence of $Nu$ on $Pr$ at moderate and large $Pr$ is in line with previous studies on non-rotating convection (Verzicco & Camussi Reference Verzicco and Camussi1999; Li et al. Reference Li, He, Tian, Hao and Huang2021). The scattering of data points from the solid line in figure 6(b) is mainly due to different Ekman numbers (different shapes in the plot), suggesting that the power law $E^2$ is also not suitable for the $E$ dependence. For fixed $E$ and $Pr$, we see that $Nu-1$ is proportional to $Ra^{3/2}$ mostly but tends to deviate from the power law $Ra^{3/2}$ and approaches the non-rotating power law $Ra^{1/3}$ at large $Ra$, although the convection remains rotationally dominated based on the criterion of $Ro_{\ell } <0.1$. We discuss the transition from rotating turbulent convection to non-rotating convection in § 3.5. In short, figure 6(b) shows that the widely used diffusion-free scaling (3.4) for the convective heat transfer does not fit our numerical models with $Pr \ge 1$ in the GT regime, which is also shown by recent laboratory experiments with moderate to high $Pr$ (Abbate & Aurnou Reference Abbate and Aurnou2023).

We now turn to examine the scaling behaviour of the convective flow speed measured by the non-zonal Reynolds number $Re_{non}$. Figures 7(a) and 7(b) show comparisons of $Re_{non}$ from numerical models in the GT regime with the VAC scaling and the CIA scaling, respectively. By comparing the results of least-squares fitting, we find that our numerical results more closely collapse on the VAC scaling. However, we note that the flow speed tends to approach the CIA scaling at low $Pr$ and large $Re_{non}$ ($\gtrsim 10^3$). This implies that the viscosity plays a non-negligible role in most of our numerical models. It is interesting that convective velocities roughly follow a unified scaling, despite the very different scaling behaviours for the heat transfer. However, we should mention that both the VAC and CIA scalings are given in terms of the flux-based Rayleigh number $Ra_Q=(Nu-1)Ra$, which already takes into account different scaling behaviours of $Nu-1$. The scaling behaviour of the convective velocities in our numerical simulations is also in agreement with recent RRBC laboratory experiments at different $Pr$ (Abbate & Aurnou Reference Abbate and Aurnou2023).

Figure 7. Non-zonal Reynolds number $Re_{non}$ versus control parameters in the GT regime: (a) $Re_{non}$ as a function of $Ra_{Q}Pr^{-2}E^{2/3}$, where the solid black line corresponds to the VAC scaling (3.2); (b) $Re_{non}$ as a function of $Ra_{Q}Pr^{-2}E^{1/2}$, where the solid black line corresponds to the CIA scaling (3.5). The dashed lines corresponds to the least-square fit to the data.

In summary for the GT regime, both heat transfer and convective velocities asymptotically approach the diffusion-free scalings at low $Pr\le 1$, suggesting that both the global heat transfer and convective flows are controlled by inviscid dynamics in the bulk. At large $Pr>1$, the efficiency of heat transfer becomes nearly independent of $Pr$ and approaches the non-rotating scaling at large $Ra$, while the convective velocities closely follow the scaling based on the VAC force balance. The scaling behaviours at large $Pr$, in line with experimental results of Abbate & Aurnou (Reference Abbate and Aurnou2023), suggest that the heat transfer is controlled by the boundary layers, whereas the typical flow speeds are controlled by the interior force balance in currently accessible numerical models.

3.5. Weakly rotating regime

Further increasing $Ra$, the convective flows become less geostrophic due to the weakening rotational constraint (figures 2d and 3d). At sufficiently large $Ra$, a strong buoyancy force would break rotational constraints and lead to non-rotating convection. However, the transition from rotating to non-rotating convection does not mean an abrupt regime change. It is not straightforward to quantitatively define when the convection is rotationally dominated. The regime changes from rotating to non-rotating convection have been extensively discussed and several transition criteria were proposed and tested (e.g. King et al. Reference King, Stellmach, Noir, Hansen and Aurnou2009; Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012; Gastine et al. Reference Gastine, Wicht and Aubert2016; Long et al. Reference Long, Mound, Davies and Tobias2020), yet no consensus has been reached. As introduced in § 3.1, we use the local Rossby number $Ro_{\ell } \le 0.1$ as a tentative criterion for the rotation-dominated convection and, thus, set cases with $Ro_{\ell } > 0.1$ as the WR regime. In this section we examine changes of the Nusselt number and local Rossby number as increasing $Ra$ from GT to WR regimes. As we shall show the transitional behaviours are rather complicated and it is impractical to define a single transition criteria that can reconcile both the heat transfer and flow morphology with different $Pr$.

In order to monitor the behaviour changes of the heat transfer, we define a rescaled Nusselt number $\widetilde {Nu}$, i.e.

(3.6)\begin{equation} \widetilde{Nu} =\left\{ \begin{array}{@{}ll} \dfrac{Nu-1}{Ra^{3/2}E^{2}} & \mathrm{when}\ Pr \ge 1, \\ \dfrac{Nu-1}{Ra^{3/2}E^{2}Pr^{-1/2}} & \mathrm{when}\ Pr \le 1, \end{array} \right. \end{equation}

which takes into account different scalings of $Nu-1$ at different $Pr$ in the GT regime as shown in figure 6. The rescaled $\widetilde {Nu}$ should be flat in the rotation-dominated convection and would start to drop as $Ra$ increases to the WR regime. The turning point is usually seen as a sign for the transition from rotation to non-rotating convection (Long et al. Reference Long, Mound, Davies and Tobias2020).

Figures 8(a) and 8(b) show the rescaled Nusselt number $\widetilde {Nu}$ and the local Rossby number $Ro_\ell$ as a function of $RaE^{8/5}$, which was proposed to be a key control parameter for the transition (Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012). We can see that the heat transfer behaviours indeed change at around $RaE^{8/5}\approx 5$ for data points with $Pr\ge 1$, but this criteria does not reconcile data points with low $Pr$. It is interesting that $RaE^{8/5}\approx 5$ also corresponds to $Ro_{\ell }\approx 0.1$ at $Pr=1$. However, the $Ro_{\ell }$ at high $Pr$ (blue and black symbols in figure 8b) is still much smaller than 0.1 even when $RaE^{8/5}>5$. This means that the convective flow remains quite geostrophic in the bulk but the heat transfer already approaches the non-rotating scaling at high $Pr$. Therefore, the transitional criterion based on the heat transfer and based on the flow morphology would be quite different. This again points to the scenario that the heat transfer is controlled by the boundary layers but the convective flows are controlled by the force balance in the bulk.

Figure 8. Rescaled Nusselt number $\widetilde {Nu}$ and local Rossby number $Ro_\ell$ in the GT and WR regimes: (a) $\widetilde {Nu}$ and (b) $Ro_\ell$ as a function of $RaE^{8/5}$; (c) $\widetilde {Nu}$ and (d) $Ro_\ell$ as a function of $RaE^{8/5}Pr^{-4/5}$. Data points with $Nu<2$ are excluded in (a) and (c). Symbols have the same definitions as those in figure 1.

Figure 8(b) clearly shows that the $Ro_{\ell }$ depends on $Pr$. By an empirical fitting of data, we find that the $Pr$ dependence can be approximated by a power exponent of around 4/5. Therefore, we define a control parameter

(3.7)\begin{equation} Ra_G=RaE^{8/5}Pr^{-4/5}, \end{equation}

and plot $Nu$ and $Ro_{\ell }$ as a function of $Ra_G$ in figure 8(c,d). Gastine et al. (Reference Gastine, Wicht and Aubert2016) derived a similar control parameter of $RaE^{8/5}Pr^{-3/5}$ based on the convective Rossby number $Ro_c$ in the boundary layers. Here we find that $Ra_G=RaE^{8/5}Pr^{-4/5}$ can better fit our numerical data of the $Ro_{\ell }$. Figure 8(d) shows that $Ra_G\approx 5$ corresponds to $Ro_{\ell }=0.1$ for all different $Pr$. We also see from figure 8(c) that $Ra_G\approx 5$ corresponds to the turning point of the rescaled Nusselt number $\widetilde {Nu}$ for data points with $Pr\le 1$. This implies that the transition criteria can be unified for the heat transfer and flow morphology at low $Pr$. The control parameter $Ra_G$ is not suitable to describe the transition of the heat transfer for the model with $Pr>1$ as expected.

Combining the scaling behaviours shown in § 3.4 and the transitional behaviours in this section, we suggest that the heat transfer is controlled by the boundary layers at large $Pr$ and by the interior dynamics at small $Pr$. We should mention that $Pr=1$ corresponds to a quite special parameter regime because data points can fit both $Pr$-dependent and $Pr$-independent scalings. The convective flow velocities and flow morphology are mainly controlled by the force balance in the interior. If we characterise the transition from GT to WR based on the $Ro_{\ell }$ in the bulk, we find that $Ra_G=RaE^{8/5}Pr^{-4/5}$ provides a unified control parameter for the transition with all different $Pr$ based on our numerical data.

3.6. Mean zonal flows

The above analysis has focused on the scaling behaviour of non-zonal flows, which are related to the heat transfer. Mean zonal flows are spontaneously generated in rotating convection systems due to the nonlinear effects (e.g. Christensen Reference Christensen2002; Miyagoshi, Kageyama & Sato Reference Miyagoshi, Kageyama and Sato2010). In this section we show the magnitude of zonal flows and scaling behaviours with different $Pr$ in the GT regime where the zonal flow may be significant. Previous studies have shown that zonal flows are more readily developed at low Prandtl numbers $Pr<1$ (Zhang Reference Zhang1992; Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001). The zonal flow amplitude is measured by the zonal Reynolds number $Re_{zon}$ as defined in § 2.3. Figure 9(a) shows the ratio of $Re_{zon}/Re_{non}$ as a function of $Ra$ with different $E$ and $Pr$ in the GT regime. We can see that the ratio is around or larger than unity for most cases at $Pr<1$, suggesting strong zonal flows are generated at low Prandtl numbers. In contrast, most of the cases show the ratio $Re_{zon}/Re_{non} <1$ when $Pr >1$. These results are in line with previous numerical and experimental studies on zonal flows driven by rotating convection (Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001; Christensen Reference Christensen2002). It is also found that as $Ra$ increases, $Re_{zon}/Re_{non}$ in both $Pr=1$ and $Pr = 7$ cases shows a significant increase, especially at relatively low Ekman number $E=3 \times 10^{-6}$. This demonstrates that at sufficiently small $E$ and large $Ra$, the zonal flows can become dominant over the non-zonal flows even at high $Pr$. However, within the parameter interval we calculated, dominant zonal flows are more likely to be generated at low $Pr$.

Figure 9. (a) Reynolds number ratio $Re_{zon}/Re_{non}$ as a function of $Ra$ in the GT regime; (b) $Re_{zon}$ as a function of $Ra^{1.29}E^{1.15}Pr^{-1.42}$ in the GT regime. Symbols have the same definitions as those in figure 1.

It is of great interest to see if there is any systematic dependence of the zonal flows on the control parameters. As there is no existing theoretical prediction on scaling of zonal flows in rotating turbulent convection, we made a least-square fitting of $Re_{zon}$ in the form of the power law in control parameters from numerical data with $Re_{zon} \ge 100$. As shown in figure 9(b), the power law scaling $Re_{zon} \sim Ra^{1.29}E^{1.1}Pr^{-1.42}$ (solid black line) can well describe the data. It is found that the zonal flow strength is stronger for $Pr \le 1$ and, for $Pr \ge 1$, the $Re_{zon}$ hardly exceeds 100 in our simulation interval. The scaling of $Pr^{-1.42}$ demonstrates that the strong zonal flow is more likely to be generated at $Pr \le 1$ with fixed $Ra$ and $E$. For $Pr \ge 1$, the strong zonal flow can also be generated by increasing $Ra$ with fixed $E$.