Boundary layer asymmetry in turbulent spherical Rayleigh–Bénard convection: combined dependence on Prandtl number and radius ratio

Abstract

Direct numerical simulations (DNS) are performed to investigate the dependence of the Prandtl number ($\textit{Pr}$) and radius ratio ($\eta =r_{i}/r_{o}$) on the asymmetry of the mean temperature radial profiles in turbulent Rayleigh–Bénard convection (RBC) within spherical shells. Unlike planar RBC, the temperature drop, and the thermal and viscous boundary layer thicknesses, at the inner and outer boundaries are not identical in spherical shells. These differences in the boundary layer properties in spherical RBC contribute to the observed asymmetry in the radial profiles of temperature and velocity. The asymmetry originates from the differences in curvature and gravity at the two boundaries, and in addition, is influenced by $\textit{Pr}$. To investigate the $\eta$ and $\textit{Pr}$ dependence of these asymmetries, we perform simulations of Oberbeck–Boussinesq convection for $\eta = 0.2,0.6$ and $0.1 \leqslant Pr \leqslant 50$, and for a range of Rayleigh numbers ($Ra$) varying between $5 \times 10^{6}$ and $5 \times 10^{7}$. The Prandtl numbers that we choose cover a broad range of geophysical relevance, from low-$\textit{Pr}$ regimes ($\textit{Pr}=0.1$) representative of gas giants such as Jupiter and Saturn, to high-$\textit{Pr}$ regimes characteristic of organic flows used in the convection experiments ($\textit{Pr}=50$). A centrally condensed mass, with the gravity profile $g \sim 1/r^{2}$, is employed in this study. Our results show that the asymmetry at smaller $\eta$ exhibits a stronger $\textit{Pr}$ dependence than at larger $\eta$. Various assumptions for quantifying this asymmetry are evaluated, revealing that different assumptions are valid in different $\textit{Pr}$ regimes. It is shown that the assumption of the equal characteristic plume separation at the inner and outer boundaries, as well as the assumption of the identical thermal fluctuation scales between the two boundary layers, is valid only for $0.2 \lesssim Pr \lesssim 1$. In contrast, assumptions based on the equivalency of the local thermal boundary layer Rayleigh numbers and laminar natural-convective boundary layers are validated at $\textit{Pr}=50$ for the explored parameter space. Furthermore, new assumptions based on the statistical analysis of the inter-plume islands are proposed for $\textit{Pr}=0.1$ and $50$, and these are validated against the DNS data. These findings provide insights into the $(Pr,\eta)$ dependence of asymmetry in spherical RBC, and offer a framework for studying similar systems in geophysical and astrophysical contexts.

1. Introduction

Thermal convection in spherical shells is of vital importance in geophysical and astrophysical systems, including Earth’s mantle convection (Davies & Richards 1992), planetary interiors (Aurnou et al. 2008; Wicht & Sanchez 2019), and stellar convection zones (Hanasoge, Gizon & Sreenivasan 2016). Compared to Rayleigh–Bénard convection (RBC) in planar geometries, spherical RBC exhibits two key differences in convection dynamics: (i) the scaling relationships between the system response parameters – such as the dimensionless heat transfer (Nusselt number, $\textit{Nu}$ ) and dimensionless momentum transfer (Reynolds number, $Re$ ) – and the input parameters (Rayleigh number, $Ra$ , and Prandtl number, $\textit{Pr}$ ) differ from those observed in planar RBC; (ii) the differences in the radius of curvature at the inner and outer boundaries, and the radially varying gravity, lead to the asymmetric thermal and velocity boundary layers (Gastine, Wicht & Aurnou 2015; Fu et al. 2024). This asymmetry is also influenced by $\textit{Pr}$ (Tilgner 1996). To accurately model the heat and momentum transport in spherical convection, it is imperative to gain an understanding of the asymmetry, and its dependence on $\eta$ and $\textit{Pr}$ . In this work, we make an effort in this direction.

Experimental and numerical studies have provided significant insights into the dynamics of convection in spherical shells. Experimental studies, such as the ‘GeoFlow’ experiments on the International Space Station (ISS) (Egbers et al. 2003; Futterer et al. 2010; Zaussinger et al. 2020), use electrically induced radial body forces to mimic gravity, and focus mainly on mantle convection with high $\textit{Pr}$ ( $40 \lt Pr \lt 200$ ). However, these experiments are limited to relatively low Rayleigh numbers ( $Ra \lt 10^{7}$ ) due to the constraints on the electric field strength and fluid properties. Numerical studies, on the other hand, have explored a broader parameter space. For example, Tilgner (1996) investigated the scaling of global heat and momentum transport in spherical RBC, examining $\textit{Nu}(Ra, Pr)$ and $Re(Ra, Pr)$ relations in the ranges $0.06 \leqslant Pr \leqslant 10$ and $4 \times 10^{3} \leqslant Ra \leqslant 8 \times 10^{5}$ at a fixed radius ratio $\eta =0.4$ with a linear gravity profile ( $g(r) \sim r$ ). The author reports that the scalings $\textit{Nu} \sim Ra^{0.24}$ and $Re \sim Ra^{0.5}$ are valid within the explored parameter space. Similarly, Gastine et al. (2015) performed direct numerical simulations (DNS) to study spherical RBC with various gravity profiles ( $g(r) \in \left \{r/r_{o},1,(r_{o}/r)^{2}, (r_{o}/r)^{5} \right \}$ ) and radius ratios ( $0.2 \leqslant \eta \leqslant 0.95$ ) at $\textit{Pr}=1$ . The scaling behaviours of $\textit{Nu}(Ra)$ and $Re(Ra)$ were studied exclusively for the simulations with $Ra \leqslant 10^9$ and at a specific radius ratio $\eta = 0.6$ . Recently, Fu et al. (2024) systematically investigated the $\eta$ dependence of $\textit{Nu}(Ra)$ and $Re(Ra)$ scalings for $0.2 \leqslant \eta \leqslant 0.8$ at $\textit{Pr}=1$ . They observed larger scaling exponents for smaller $\eta$ cases, suggesting that the simulations with smaller $\eta$ reach the classical $\textit{Nu}\sim Ra^{1/3}$ regime at a relatively lower $Ra$ . The $\eta$ dependence in $\textit{Nu}(Ra)$ scaling relations was also quantified. These results highlight the impact of geometric curvature on the global transport properties in spherical RBC.

Boundary layer asymmetry studies in spherical RBC have primarily focused on two specific $\textit{Pr}$ regimes: $\textit{Pr} = 1$ , and the asymptotic $\textit{Pr} \rightarrow \infty$ limit. The asymmetry of thermal boundary layers in spherical RBC is typically characterised by the two key parameters: the mid-layer bulk temperature $\vartheta _{mid}$ , and the thermal boundary layer thickness ratio $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ , where $\lambda _{\vartheta }^{i}$ and $\lambda _{\vartheta }^{o}$ represent the thicknesses of the inner and outer thermal boundary layers, respectively.

For mantle convection at $\textit{Pr} \rightarrow \infty$ , Jarvis (1993) proposed that the thermal boundary layers are equally stable, leading to equal thermal boundary layer Rayleigh numbers at both the boundaries, where the thermal boundary layer Rayleigh number is defined based on the thermal boundary layer thickness as the length scale, and the temperature drop across the thermal boundary layer as the temperature scale. This assumption was extended by Vangelov & Jarvis (1994) and Jarvis, Glatzmaierand & Vangelov (1995) to study the $\eta$ dependence of boundary layer asymmetry in spherical RBC, which was further validated through their numerical simulations. Moreover, Shahnas et al. (2008) and Deschamps, Tackley & Nakagawa (2010) extended these models to systems with internal heating, highlighting the additional complexities introduced by internal heat sources.

For $\textit{Pr} = 1$ , Gastine et al. (2015) assumed equivalent mean plume separation at the edges of the inner and outer thermal boundary layers, and derived thermal boundary layer asymmetry as a function of $\eta$ , validating their predictions with DNS results. They also examined the assumption proposed by Vangelov & Jarvis (1994), and found that it did not align with their DNS data at $\textit{Pr} = 1$ . Furthermore, they tested the assumption of the equivalency of the temperature fluctuation scales at both boundary layers (Wu & Libchaber 1991), and observed that it fits DNS data only for a gravity profile $g(r) \sim 1/r^{2}$ at $\textit{Pr} = 1$ .

Despite these efforts, existing models have been validated only at specific $\textit{Pr}$ values, either $\textit{Pr} = 1$ or $\textit{Pr} \rightarrow \infty$ . Whether these models remain applicable to intermediate $\textit{Pr}$ regimes or require modifications remains an open question. This knowledge gap motivates the present study, which systematically investigates the $\textit{Pr}$ and $\eta$ dependence of thermal boundary layer asymmetry across a broader parameter space.

In astrophysical and geophysical contexts, $\textit{Pr}$ exhibits significant variations due to differences in the material properties of the flowing medium. For example, in the Earth’s outer core, $\textit{Pr}$ is estimated to lie between values $0.1$ and $1$ (Wicht & Sanchez 2019), while in the convective molecular layers of Jupiter and Saturn, the typical values are $\textit{Pr} \approx 0.1$ (Yadav, Heimpel & Bloxham 2020). The Earth’s near-surface atmosphere has Prandtl number approximately $0.7$ , whereas seawater has $\textit{Pr} \approx 7$ (at 20 $^\circ$ C) and $\textit{Pr} \approx 13$ (at 0 $^\circ$ C). Light organic liquids span a broader range, with $5 \leqslant Pr \leqslant 50$ (Shah & London 1978). Silicone oil with $\textit{Pr} \approx 64.64$ was used in the GeoFlow experiment aboard the ISS to study mantle-like convection in a spherical shell geometry (Futterer et al. 2010). In addition to a broad range of $\textit{Pr}$ values, the radius ratio ( $\eta = r_{i}/r_{o}$ ), which characterises the geometric configuration of spherical shells, also varies significantly among planetary bodies. Table 1 summarises the radius ratios and Prandtl numbers of selected planetary interiors. Since the asymmetry of the thermal boundary layers depends strongly on both $\textit{Pr}$ and $\eta$ , understanding how this asymmetry evolves with these parameters is essential for characterising the dynamics of planetary convection. Despite its significance, the combined dependence of thermal boundary layer asymmetry on $\eta$ and $\textit{Pr}$ regimes has not been investigated systematically, leaving ample room for further research.

Table 1. Radius ratios for the selected planetary interiors in the solar system. The radius ratios for Jupiter (Guillot et al. 2004) and Saturn (Christensen & Wicht 2008; Vazan et al. 2016) are chosen as the ratio between the radii of the inner metallic-core boundary and the outer upper atmosphere boundary.

This study investigates the dependence of thermal boundary layer asymmetry on $\textit{Pr}$ and $\eta$ through DNS. The simulations consider radius ratios $\eta = 0.2$ and $0.6$ , spanning $0.1 \leqslant Pr \leqslant 50$ and $Ra \leqslant 5 \times 10^{7}$ . Given the ranges of $\textit{Pr}$ and $\eta$ observed in natural systems, as discussed earlier, the chosen parameter space provides a reasonable representation of diverse convection scenarios. We assess the validity of the existing models, and propose new assumptions based on the statistical analysis of the inter-plume islands at the edges of the thermal boundary layers to explain the combined dependence of thermal boundary layer asymmetry on $\textit{Pr}$ and $\eta$ . Specifically, we develop new models for $\textit{Pr}=0.1$ and $50$ , and validate them using our simulation data.

The remainder of this paper is organised as follows. Section 2 describes the governing equations, diagnostic parameters and numerical methods used in this study. Section 3 gives an overview of the boundary layer asymmetry for different $\eta$ and $\textit{Pr}$ regimes. In § 4, we evaluate several existing models for quantifying thermal boundary layer asymmetry over $0.1 \leqslant Pr \leqslant 50$ , identifying their respective validation ranges. Section 5 introduces new thermal boundary layer asymmetry models for $\textit{Pr}=0.1$ and $50$ , developed through statistical analysis of the inter-plume islands near the inner and outer boundary layer edges. This section also explores the dependence of the inter-plume island statistics on $\textit{Pr}$ and $\eta$ . Finally, § 6 summarises our findings and discusses future research directions.

2. Model description

2.1. Governing equations

We study Oberbeck–Boussinesq convection in spherical shells with inner radius $r_{i}$ and outer radius $r_{o}$ . A no-slip boundary condition is employed at both boundaries for the velocity field. The temperature boundary condition is isothermal, with the inner boundary at $T_{i}$ , and the outer boundary at $T_{o}$ . The governing equations are non-dimensionalised by the temperature scale $\Delta T = T_{i}-T_{o}$ , the length scale $d=r_{o}-r_{i}$ , the viscous dissipation time scale $d^{2}/\nu$ , the momentum diffusive velocity scale $\nu /d$ , and the pressure scale $\rho \nu ^{2} / d^{2}$ . Gravity is non-dimensionalised by its value at the outer boundary $g_{o}=g(r_{o})$ . The dimensionless governing equations read

(2.1) \begin{equation} \boldsymbol{\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{u} =0, \end{equation}

(2.2) \begin{equation} \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{ \nabla } \boldsymbol{u} = -\boldsymbol{\nabla} p + \frac {Ra}{Pr}gT\boldsymbol{e}_r+\nabla^2 \boldsymbol{u}, \end{equation}

(2.3) \begin{equation} \frac {\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla} T = \frac {1}{Pr} \nabla^2 T. \end{equation}

In the above equations, $\boldsymbol{u}, T,p$ represent the velocity, temperature and pressure, respectively. The equations are controlled by three input parameters,

(2.4) \begin{equation} Ra=\frac {\alpha g_o\, \Delta T\, d^3}{\nu \kappa }, \quad Pr=\frac {\nu }{\kappa }, \quad \eta =\frac {r_i}{r_o}, \end{equation}

where $\alpha$ is the thermal expansion coefficient, and $\nu$ and $\kappa$ are the kinematic viscosity and thermal diffusivity, respectively. In our simulations, we assume a centrally condensed mass, with gravity profile $g(r)=(r_{o}/r)^{2}$ (Chandrasekhar 1981). A schematic of the computational domain and the coordinate system is provided in Appendix A.

2.2. Response parameters

We use the following notations throughout the paper for different averaging procedures. An overbar ( $\overline {{\cdot}}$ ) represents the time average of a dimensionless variable $f$ , while $\left \langle {\cdot }\right \rangle _s$ and $\left \langle {\cdot }\right \rangle$ represent the horizontal and volume averages, respectively. The definitions of these averaging procedures are

(2.5) \begin{eqnarray} \overline {f}=\frac {1}{\tau } \int _{t_0}^{t_0+\tau } f\mathrm{d} t, \end{eqnarray}

(2.6) \begin{eqnarray} \left \langle f \right \rangle _{s}=\frac {1}{4 \unicode{x03C0} } \int _{0}^{\unicode{x03C0} } \int _{0}^{2 \unicode{x03C0} } f \sin\theta\mathrm{d}\theta\mathrm{d}\phi, \end{eqnarray}

(2.7) \begin{eqnarray} \left \langle f \right \rangle=\frac {1}{V} \int _{r_{i}}^{r_{o}} \int _{0}^{\unicode{x03C0} } \int _{0}^{2 \unicode{x03C0} } f\sin\theta r^{2}\mathrm{d}r\mathrm{d}\theta\mathrm{d}\phi. \end{eqnarray}

In the above equations, $\theta$ is co-latitude and $\phi$ is longitude, and $V=4\unicode{x03C0} (r_{o}^3-r_{i}^3)/3$ is the volume of the spherical shell.

There are two key response parameters in RBC: the Nusselt number ( $\textit{Nu}$ ), which represents the dimensionless heat flux, and the Reynolds number ( $Re$ ), which is a dimensionless measure of flow velocity.

In Oberbeck–Boussinesq approximation, $\textit{Nu}$ for isothermal boundaries in spherical RBC reads

(2.8) \begin{equation} \textit{Nu} = \frac {\overline {\left \langle u_rT \right \rangle _s}-\dfrac{1}{Pr}\dfrac{{\rm d}\vartheta }{{\rm d}r}}{-\dfrac{1}{Pr}\dfrac{{\rm d}T_c}{{\rm d}r}} = \left. -\eta \frac {{\rm d} \vartheta }{{\rm d}r} \right |_{r=r_i} = \left. -\frac {1}{\eta } \frac {{\rm d} \vartheta }{{\rm d}r} \right |_{r=r_o}, \end{equation}

where

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

is the conductive temperature profile.

The Reynolds number is defined as

(2.10) \begin{equation} Re = \overline {\sqrt {\langle u^2 \rangle }} = \overline {\sqrt {\big\langle u_r^2+u_{\theta }^2+u_{\phi }^2 \big\rangle }}. \end{equation}

By adopting the notations (2.5) and (2.6), the time and horizontally averaged radial profiles of the mean temperature and the mean horizontal velocity are given, respectively, by

(2.11) \begin{equation} \vartheta (r) = \overline {\left \langle T \right \rangle _{s}}, \quad Re_h(r)=\overline {\sqrt {\big\langle u_{\theta }^2+u_{\phi }^2 \big\rangle_s }}. \end{equation}

2.3. Numerical settings and parameter space

In this work, all the simulations are conducted by using MagIC (https://magic-sph.github.io/), which is a pseudo-spectral solver in which all the unknown variables are expanded into complete sets of functions in the radial and horizontal directions. Chebyshev polynomials are applied in the radial direction, while spherical harmonic functions are used in the azimuthal and latitudinal directions. The equations are time-stepped by advancing the nonlinear terms using an explicit second-order Adams–Bashforth scheme, while the linear terms are time-advanced using an implicit Crank–Nicolson algorithm. At each time step, the linear terms are calculated in the spectral space, while the nonlinear terms are calculated in the physical space. For more information about MagIC, the interested reader is referred to Wicht (2002), Christensen & Wicht (2007) and Lago et al. (2021).

For $g(r)=(r_{o}/r)^2$ , if we dot (2.2) with $\boldsymbol{u}$ , and (2.3) with $T$ , and average the results temporally and volumetrically, then the following exact relations can be derived (Gastine et al. 2015):

(2.12) \begin{eqnarray} \epsilon _u &=& \frac {3}{1+\eta +\eta ^2} \frac {Ra}{Pr^2} (\textit{Nu}-1), \end{eqnarray}

(2.13) \begin{eqnarray} \epsilon _{\vartheta } &=& \frac {3\eta }{1+\eta +\eta ^2}\, \textit{Nu}. \end{eqnarray}

Here, $\epsilon _u$ and $\epsilon _{\vartheta }$ represent the time- and volume-averaged kinetic energy dissipation rate and thermal dissipation rate, respectively. Following our non-dimensionalisation, $\epsilon _u$ and $\epsilon _{\vartheta }$ are calculated as

(2.14) \begin{equation} \epsilon _u = \overline {\bigg\langle \left (\frac {\partial u_i}{\partial x_j} \right)^2 \bigg\rangle }, \quad \epsilon _{\vartheta } = \overline {\bigg\langle \left (\frac {\partial T}{\partial x_i} \right)^2 \bigg\rangle }. \end{equation}

The Nusselt number is computed based on the gradient of temperature at the boundaries as well as the dissipation rates:

(2.15) \begin{equation} \left. \begin{array}{ll} \displaystyle \textit{Nu}_{\epsilon _u} &\!\!\!\!\!\,= \dfrac {1+\eta +\eta ^2}{3} \dfrac{Pr^2}{Ra} \epsilon _u +1, \quad \textit{Nu}_{\epsilon _{\vartheta }} = \dfrac{1+\eta +\eta ^2}{3\eta } \epsilon _{\vartheta },\\ \\ \displaystyle \textit{Nu}_{b} &\!\!\!\!\!\,= \left. -\eta \dfrac {{\rm d} \vartheta }{{\rm d}r} \right |_{r=r_i}, \quad \textit{Nu}_{t} = \left. -\dfrac {1}{\eta } \dfrac {{\rm d} \vartheta }{{\rm d}r} \right |_{r=r_o}, \end{array}\right \} \end{equation}

where $\textit{Nu}_{b}$ and $\textit{Nu}_{t}$ follow from the definition of $\textit{Nu}$ (see (2.8)). The relative error between these $\textit{Nu}$ values is used as a measure of resolution in our simulations. A simulation is considered converged when the relative difference between these independent Nusselt number values is consistently $\lesssim 1\, \%$ . Additionally, we also ensured convergence by examining the balance of the kinetic energy budget, which is the balance between the volume-averaged input buoyancy power and the volume-averaged kinetic energy dissipation rate. See Appendix A.

In addition, to validate the numerical accuracy of MagIC, we conducted two benchmark simulations using the open-source spectral code Dedalus (https://dedalus-project.org/). In the first case ( $\textit{Pr} = 1$ , $\eta = 0.6$ , $Ra = 10^6$ ), MagIC produced $\textit{Nu} = 8.94$ and $Re = 260.46$ , while Dedalus gave $\textit{Nu} = 8.92$ and $Re = 260.36$ . In the second case ( $\textit{Pr} = 10$ , $\eta = 0.6$ , $Ra = 10^6$ ), MagIC yielded $\textit{Nu} = 9.17$ and $Re = 34.49$ , compared to $\textit{Nu} = 9.05$ and $Re = 34.14$ from Dedalus. These results show excellent agreement between the two codes, further confirming the reliability of MagIC for the problem considered in this study.

We conduct seven sets of simulations at $\textit{Pr}=0.1, 0.2, 0.5, 1, 5, 10,50$ with $5 \times 10^{6} \leqslant Ra \leqslant 5 \times 10^{7}$ , and $\eta =0.2$ and $0.6$ . For more details about the simulation parameters and grid resolution, refer to Table 3 in Appendix A.

3. Velocity and temperature asymmetry

Spherical shell RBC exhibits notable differences from its planar counterpart, primarily due to the curvature of the bounding surfaces, and the radial variation of gravity. This results in asymmetry in the radial profiles of temperature and velocity fluctuations.

Figure 1. Radial profiles of (a) temperature $\vartheta (r)$ , and (b) velocity $Re(r)$ , at different $\eta$ and $\textit{Pr}$ . Here, $Ra=10^{7}$ for all the cases shown.

Figure 1 shows the time and horizontally averaged radial profiles of temperature and velocity at $Ra=10^{7}$ , $\eta =0.2,0.6$ , at different $\textit{Pr}$ . As shown in figure 1(a), the majority of the temperature drop occurs near the boundaries. This behaviour can be attributed to the thermal shortcut (Ahlers, Grossmann & Lohse 2009) in the bulk flow. In turbulent RBC, temperature is well mixed in the bulk region, leading to a minimal temperature drop within the bulk, facilitating the majority of the temperature drop to occur within the thermal boundary layers. The mean temperature drop near the inner boundary is larger than that near the outer boundary, consistent with the observations reported by Gastine et al. (2015) and Fu et al. (2024). This phenomenon arises due to the conservation of heat flux, as well as the smaller surface area of the inner boundary compared to the outer boundary. Additionally, larger differences in curvature and gravity between the inner and outer boundaries yield larger deviations of the bulk temperature $\vartheta _{mid}=\vartheta ((r_{i}+r_{o})/2)$ , from $(T_{i}+T_{o})/2$ . These differences along with the corresponding deviations in $\vartheta _{mid}$ become more pronounced at smaller $\eta$ . This effect has also been observed in previous studies (Gastine et al. 2015; Fu et al. 2024). From figure 1(a), it can be seen clearly that the bulk temperature increases with increasing $\textit{Pr}$ . For example, at $Ra=10^{7}$ , $\eta =0.2$ , $\textit{Pr}=0.1$ , the temperature drop across the outer thermal boundary layer is $9.7\, \%$ of the total temperature drop across the shell. On the other hand, at $\eta =0.2,\ Pr=50$ , the temperature drop across the outer thermal boundary layer becomes $16.2\, \%$ of the total temperature drop. As $\textit{Pr}$ is varied over two and a half orders of magnitude, the temperature drop across the outer thermal boundary layer varies by approximately $67\, \%$ .

The asymmetry is also present in the time and horizontally averaged radial velocity profile $Re(r)$ . As shown in figure 1(b), two peaks are observed in the radial velocity profile. The peak value near the inner boundary is higher than that near the outer boundary, indicating greater turbulence intensity in the vicinity of the inner boundary. At fixed $Ra$ and $\textit{Pr}$ , the difference between the inner and outer peak values is higher at smaller $\eta$ , which is again due to the larger differences in curvature and gravity between the inner and outer boundaries at smaller $\eta$ . It is also observed that $Re$ is higher at lower $\textit{Pr}$ for fixed $Ra$ and $\eta$ . This occurs as the viscous effects are less significant compared to the thermal diffusion at lower $\textit{Pr}$ .

The asymmetry can also be quantified by the ratios of the inner and outer thermal and viscous boundary layer thicknesses. These are denoted as $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ and $\lambda _{u}^{i}/\lambda _{u}^{o}$ , where $\lambda _{\vartheta }$ and $\lambda _{u}$ represent the thermal and viscous boundary layer thicknesses, respectively. The superscripts $i$ and $o$ refer to the boundary layer thicknesses calculated at the inner and outer boundaries, respectively. We define the thermal and viscous boundary layer thicknesses using the slope method (Shishkina et al. 2010; Gastine et al. 2015). The thermal boundary layer thickness is defined as the distance where the linear fit of $\vartheta (r)$ at the boundary intersects the horizontal line $r=r_{mid}$ . The time and horizontally averaged radial profile of the horizontal velocity, $Re_{h}(r)$ , is used to characterise the viscous boundary layer. The viscous boundary layer is defined as the region between the boundary and the point where the linear fit of $Re_{h}(r)$ near the boundary intersects the horizontal line drawn from the local peak of the curve.

Since the temperature drop is chiefly concentrated near the boundaries, it is reasonable to assume that the entire temperature drop occurs within the thermal boundary layers, leading to

(3.1) \begin{equation} \Delta \vartheta _{i} + \Delta \vartheta _{o} = 1 \end{equation}

and

(3.2) \begin{equation} \Delta \vartheta _{i} = 1-\vartheta _{mid}, \quad \Delta \vartheta _{o} = \vartheta _{mid}. \end{equation}

Figure 2. Bulk temperature $\vartheta _{mid}=\vartheta ((r_{i}+r_{o})/2)$ with respect to $\textit{Pr}$ at different $\eta$ , for $Ra \geqslant 5 \times 10^{6}$ .

Figure 3. Ratio of the inner and the outer (a) thermal boundary layer thickness $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ , and (b) viscous boundary layer thickness $\lambda _{u}^{i}/\lambda _{u}^{o}$ , at different $\eta$ and $\textit{Pr}$ for $Ra \geqslant 5 \times 10^{6}$ .

Figure 2 shows $\vartheta _{mid}$ as a function of $\textit{Pr}$ at $\eta =0.2, 0.6$ for $Ra \geqslant 5 \times 10^{6}$ . As discussed earlier in this section, $\vartheta _{mid}$ increases with increasing $\textit{Pr}$ . For the case $\eta =0.2$ , as $\textit{Pr}$ increases, $\vartheta _{mid}$ rises from $0.097$ at $\textit{Pr}=0.1$ to $0.162$ at $\textit{Pr}=50$ , with absolute difference $0.065$ , corresponding to an approximately $67\, \%$ increase. On the other hand, for $\eta =0.6$ , $\vartheta _{mid}$ changes from $0.326$ at $\textit{Pr}=0.1$ to $0.373$ at $\textit{Pr}=50$ , resulting in absolute difference $0.047$ , which corresponds to an increase of approximately $14.4\, \%$ . This illustrates that $\textit{Pr}$ dependence of $\vartheta _{mid}$ changes with $\eta$ . This result aligns with our expectations, as when $\eta \rightarrow 1$ , $\vartheta _{mid} \rightarrow 0.5$ , leading to a reduced $\textit{Pr}$ dependence of $\vartheta _{mid}$ .

The boundary layer thickness ratios $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ and $\lambda _{u}^{i}/\lambda _{u}^{o}$ are shown in figure 3. As mentioned above, the chosen Rayleigh numbers are $Ra \geqslant 5 \times 10^{6}$ . At the same $\textit{Pr}$ and $\eta$ , the values of $\lambda _{u}^{i}/\lambda _{u}^{o}$ are almost identical to those of $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ . As shown in the figure, both $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ and $\lambda _{u}^{i}/\lambda _{u}^{o}$ decrease with increasing $\textit{Pr}$ , indicating that the asymmetry in boundary layer thicknesses becomes more pronounced at higher $\textit{Pr}$ . As an example, $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ decreases from $0.387$ at $\textit{Pr}=0.1$ to $0.207$ at $\textit{Pr}=50$ for $\eta =0.2$ , with absolute difference $0.18$ and relative change $86.9\, \%$ . In the case $\eta =0.6$ , $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ decreases from $0.722$ at $\textit{Pr}=0.1$ to $0.612$ at $\textit{Pr}=50$ , with absolute difference $0.11$ and relative change $15.2\, \%$ . The differences in $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ and $\lambda _{u}^{i}/\lambda _{u}^{o}$ with respect to $\textit{Pr}$ are smaller for $\eta =0.6$ compared to $\eta =0.2$ . This indicates a weaker $\textit{Pr}$ dependence of the boundary layer thickness ratios at higher $\eta$ . This behaviour is expected, as $\eta \rightarrow 1$ , $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o} \rightarrow 1$ and $\lambda _{u}^{i}/\lambda _{u}^{o} \rightarrow 1$ , indicating weaker $\textit{Pr}$ dependence in the $\eta \rightarrow 1$ limit.

4. Evaluation of models

In the following sections, we focus on the thermal boundary layer asymmetry, characterised by $\Delta \vartheta _{i},\ \Delta \vartheta _{o}$ and $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ , as the primary metrics.

In order to determine the three unknowns $\Delta \vartheta _{i},\ \Delta \vartheta _{o}$ and $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ , (3.1)–(3.2) can be utilised. However, an additional relation is still required. Several assumptions have been proposed to derive the third relation. For instance, Gastine et al. (2015) derive the third relation based on the assumption of equivalence of the mean plume density within the inner and outer thermal boundary layers. This assumption was validated in their simulations at $\textit{Pr}=1$ . Similarly, Wu & Libchaber (1991), in their experimental study on non-Boussinesq effects with low-temperature helium gas, found that the two thermal boundary layers adjust their temperature drops and thicknesses such that their thermal fluctuations are equal. This finding was later confirmed by Zhang, Childress & Libchaber (1997) in a follow-up experimental study using glycerol. Therefore, the equivalence of thermal fluctuation scales can also provide the third relation. Puthenveettil & Arakeri (2005) derived expressions for the heat flux through the thermal boundary layers based on the laminar natural convection boundary layer model of Rotem & Claassen (1969). By applying the principle of the conservation of heat flux, the third relation can be obtained by matching the heat fluxes at the inner and outer thermal boundary layers. In a two-dimensional cylindrical geometry with an infinite $\textit{Pr}$ , Jarvis (1993) proposed that the inner and outer thermal boundary layers are equally stable, resulting in equal thermal boundary layer Rayleigh numbers at both the surfaces. In a follow-up study, Vangelov & Jarvis (1994) applied this assumption to a two-dimensional asymmetric spherical model to obtain the third relation. This assumption was shown to work well in characterising the thermal boundary layer asymmetry in their simulations with an infinite $\textit{Pr}$ .

However, none of these assumptions has been validated across different $\textit{Pr}$ regimes. The range of applicability of the aforementioned assumptions remains unclear. In this section, we will examine the validity of these assumptions at different $\textit{Pr}$ .

4.1. Inter-plume density equivalence model

Proposed by Gastine et al. (2015), the inter-plume density equivalence model assumes that the characteristic plume density is equivalent at the edges of both thermal boundary layers. They estimate the characteristic plume separation based on the peak values of the probability distribution function (PDF) of the inter-plume areas at the edges of thermal boundary layers. The characteristic inter-plume density is defined as the number of plumes per unit surface area. Statistics were analysed at a fixed $\textit{Pr}$ of unity, with different radius ratios and gravity profiles. Their findings indicate that the characteristic inter-plume area is identical at the edges of the inner and outer thermal boundary layers, suggesting that the characteristic plume density is conserved across these surfaces. For the expressions of the characteristic inter-plume distance at the inner and outer surfaces, they use the relations developed by Puthenveettil et al. (2011) based on the laminar natural convection boundary layer model of Rotem & Claassen (1969). This approach leads to (Gastine et al. 2015)

(4.1) \begin{equation} \rho _{p}^{i} = \rho _{p}^{o} \quad \rightarrow \quad \frac {\alpha g_{i}\, \Delta \vartheta _{i}\, ({\lambda _{\vartheta }^{i}})^{5}}{\nu \kappa } = \frac {\alpha g_{o}\, \Delta \vartheta _{o}\, ({\lambda _{\vartheta }^{o}})^{5}}{\nu \kappa }, \end{equation}

where $\rho _{p}^{i},\ \rho _{p}^{o}$ are plume densities at the edges of the inner and outer thermal boundary layers. Combining (3.1)–(3.2) and (4.1), together with $g_{i}/g_{o}=(r_{o}/r_{i})^{2}=1/\eta ^{2}$ , we obtain the relations (Gastine et al. 2015)

(4.2) \begin{equation} \Delta \vartheta _{i} = \frac {1}{1+\eta ^{4/3}}, \quad \Delta \vartheta _{o}= \vartheta _{mid} = \frac {\eta ^{4/3}}{1+\eta ^{4/3}}, \quad \frac {\lambda _{\vartheta }^{i}}{\lambda _{\vartheta }^{o}}=\eta ^{2/3}. \end{equation}

These relations quantify the thermal boundary layer asymmetry derived from this model. The validation of this model will be discussed in § 4.5. The assumption on which this model rests, namely the equivalency of the characteristic plume density at the two boundaries, will also be examined carefully in § 5.

4.2. Thermal fluctuation equivalence model

The thermal fluctuation equivalence model was initially proposed by Castaing et al. (1989). They divided the flow domain into three regions: the boundary layers adjacent to the boundary, the mixing zone where plumes detach from the boundary layer and accelerate towards the bulk, and the turbulent bulk. Based on their theoretical arguments, the characteristic temperature fluctuation scale in the mixing zone, $\vartheta ^{mix}$ , was derived. They further assumed that $\vartheta ^{mix}$ is equivalent to the characteristic temperature fluctuation scale in the turbulent bulk, $\vartheta ^{bulk}$ , resulting in the relation

(4.3) \begin{equation} \vartheta ^{bulk} = \vartheta ^{mix} = \frac {\kappa \nu }{\alpha g \lambda _{\vartheta }^{3}}, \end{equation}

where symbols have their usual meanings. For the first time, Wu & Libchaber (1991) applied and successfully validated this relation in their studies of non-Boussinesq convection in helium. Zhang et al. (1997) later confirmed the validity of this assumption in the experiments using glycerol. Gastine et al. (2015) applied this assumption to determine the boundary layer temperature drop in spherical RBC with different gravity profiles ( $g \sim r, 1, 1/r^{2},1/r^{5}$ ) and $0.2 \leqslant \eta \leqslant 0.95$ , at $\textit{Pr}=1$ . It was observed that this assumption held true only when $g \sim 1/r^{2}$ . Wang et al. (2022) applied this assumption in axially co-rotating centrifugal convection system with cold inner and hot outer cylinders. The temperature drop and boundary layer thickness ratio derived from this assumption showed good agreement with their DNS data. In this study, we investigate the range of $\textit{Pr}$ for which this assumption is applicable in spherical RBC with $g \sim 1/r^{2}$ .

Invoking (4.3) at both the inner and outer thermal boundary layers yields

(4.4) \begin{equation} \vartheta ^{i}=\vartheta ^{bulk}=\vartheta ^{o}\quad \rightarrow\quad \frac {\kappa \nu }{\alpha g_{i} (\lambda _{\vartheta }^{i})^{3}} = \frac {\kappa \nu }{\alpha g_{o} (\lambda _{\vartheta }^{o})^{3}}. \end{equation}

Combining (4.4) with (3.1)–(3.2) leads to the same relations as those given by (4.2), with $g_{i}/g_{o}=1/\eta ^{2}$ .

4.3. Inter-plume distance $Ra$ equivalence model

Using the similarity solutions for a laminar natural-convection boundary layer in the limit $\textit{Pr} \rightarrow \infty$ (Rotem & Claassen 1969), Puthenveettil & Arakeri (2005) derived the following expression for $\textit{Nu}$ at the boundaries:

(4.5) \begin{equation} \textit{Nu} = -\frac {5 \times 2^{-14/15}}{3}\, Ra_{l}^{-2/15}\, Ra^{1/3}\, C(Pr), \end{equation}

where $Ra_{l}=\alpha g\,\Delta \vartheta _{w}\, l^{3} / (\nu \kappa)$ is the Rayleigh number based on the mean inter-plume distance near the thermal boundary layer, and $\Delta \vartheta _{w}$ is the temperature drop across the thermal boundary layer. Note that (4.5) applies only for large values of $\textit{Pr}$ . The coefficient $C(Pr)$ is a function of $\textit{Pr}$ but remains nearly constant when $\textit{Pr} \gt 5$ (Rotem & Claassen 1969; Puthenveettil & Arakeri 2005), with asymptotic value $C(Pr \rightarrow \infty) = -0.4601$ .

The conservation of $\textit{Nu}$ at the inner and outer boundaries yields

(4.6) \begin{equation} \textit{Nu}_{i}=\textit{Nu}_{o} \quad \implies \quad Ra_{l}^{i} = Ra_{l}^{o} \quad \implies \quad \frac {\alpha g_{i}\, \Delta \vartheta _{i}\, \overline {l}_{i}^{3}}{\nu \kappa } = \frac {\alpha g_{o}\, \Delta \vartheta _{o}\, \overline {l}_{o}^{3}}{\nu \kappa }, \end{equation}

where

(4.7) \begin{equation} \overline {l}_{i} \sim \sqrt {\frac {\alpha g_{i}\, \Delta \vartheta _{i}\, (\lambda _{\vartheta }^{i})^{5}}{\nu \kappa }}, \quad \overline {l}_{o} \sim \sqrt {\frac {\alpha g_{o}\, \Delta \vartheta _{o}\, (\lambda _{\vartheta }^{o})^{5}}{\nu \kappa }} \end{equation}

represent the mean inter-plume distances at the inner and outer thermal boundary layers (Puthenveettil et al. 2011).

By combining (3.1)–(3.2) with (4.6)–(4.7), we obtain the relations

(4.8) \begin{equation} \Delta \vartheta _{i} = \frac {1}{1+\eta }, \quad \Delta \vartheta _{o}=\vartheta _{mid}=\frac {\eta }{1+\eta }, \quad \frac {\lambda _{\vartheta }^{i}}{\lambda _{\vartheta }^{o}}=\eta. \end{equation}

Equation (4.8) describes the asymmetry of the thermal boundary layer.

4.4. Boundary layer $Ra$ equivalence model

In the mantle convection with an infinite $\textit{Pr}$ , it is hypothesised that the temperature field adjusts so that the local thermal boundary layer Rayleigh numbers are equal at both boundaries. The thermal boundary layer Rayleigh numbers are defined by using the thermal boundary layer thickness as the length scale, and the thermal boundary layer temperature drop as the temperature scale. This assumption implies that neither boundary layer is more unstable than the other, leading to the relation

(4.9) \begin{equation} Ra_{\lambda _{\vartheta }}^{i} = Ra_{\lambda _{\vartheta }}^{o} \quad \rightarrow \quad \frac {\alpha g_{i}\, \Delta \vartheta _{i}\, (\lambda _{\vartheta }^{i})^{3}}{\nu \kappa } = \frac {\alpha g_{o}\, \Delta \vartheta _{o}\, (\lambda _{\vartheta }^{o})^{3}}{\nu \kappa }. \end{equation}

Jarvis (1993) first applied this assumption in a two-dimensional, cylindrical shell model of mantle convection with uniform gravity. He derived the thermal boundary layer asymmetry and found a good agreement with his data. Vangelov & Jarvis (1994) further applied this assumption in an axisymmetric spherical model of mantle convection, achieving satisfactory results. Combining (3.1)–(3.2) and (4.9) for a spherical shell with $g \sim 1/r^{2}$ , we recover (4.8), implying that the assumptions in this subsection and § 4.3 lead to identical results for thermal boundary layer asymmetry.

The validation range (with respect to $\textit{Pr}$ ) of all the models presented above will be discussed in the next subsection.

Figure 4. (a) Compensated plot of $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ and (b) compensated plot of $\vartheta _{mid}$ , based on (4.2).

4.5. Model validation

In this study, we examine four widely used assumptions in spherical RBC to quantify the asymmetry of the thermal boundary layers. These four assumptions lead to two different relations as presented in (4.2) and (4.8). These relations can be expressed generally as

(4.10) \begin{equation} \frac {\lambda _{\vartheta }^{i}}{\lambda _{\vartheta }^{o}} = \eta ^{\gamma }, \quad \Delta \vartheta _{i} = \frac {1}{1+\eta ^{(2-\gamma)}}, \quad \Delta \vartheta _{o} = \vartheta _{mid} = \frac {\eta ^{(2-\gamma)}}{1+\eta ^{(2-\gamma)}}, \end{equation}

where $\gamma = \gamma _{1} = 2/3$ for the inter-plume density equivalence and thermal fluctuation equivalence models, and $\gamma = \gamma _{2} = 1$ for the inter-plume distance $Ra$ equivalence and boundary layer $Ra$ equivalence models.

Figure 4 shows the compensated plots of $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ and $\vartheta _{{mid}}$ , as calculated using (4.2), with respect to $\textit{Pr}$ . It is observed that the thermal boundary layer asymmetry relations derived from the inter-plume density equivalence model and the thermal fluctuation equivalence model show good agreement with the DNS data within the range $0.2 \lesssim Pr \lesssim 1$ . The deviation becomes increasingly significant outside this range. As we will demonstrate later, in § 5, this breakdown arises because the fundamental assumption underlying the model – the equivalence of the characteristic inter-plume distance at the inner and outer boundaries – no longer holds outside this $\textit{Pr}$ regime. Gastine et al. (2015) similarly reported that the inter-plume density equivalence model performs well at $\textit{Pr} = 1$ . In § 5, we further evaluate this underlying assumption across a range of $\textit{Pr}$ , and demonstrate that it does not hold for $1 \lt Pr \leqslant 50$ and $\textit{Pr} = 0.1$ .

Figure 5(a) shows $\vartheta ^{i}/\vartheta ^{o}$ as a function of $\textit{Pr}$ . The results indicate that the inner and outer thermal fluctuation scales are nearly equivalent for $0.2 \lesssim Pr \lesssim 1$ , but exhibit increasing deviations outside this range. Overall, $0.2 \lesssim Pr \lesssim 1$ is identified as a suitable validation range for both the inter-plume density equivalence model and the thermal fluctuation equivalency model. The reason why the thermal fluctuation equivalence model fails outside this $\textit{Pr}$ regime can be understood by examining the plume dynamics and boundary layer structures. For very low $\textit{Pr}$ (e.g. $\textit{Pr}=0.1$ ), thermal diffusivity dominates, causing the plumes to broaden significantly. This makes the boundary layers more susceptible to the large-scale disturbances, and thus disrupts the symmetry of the thermal fluctuation scales. At very high $\textit{Pr}$ (e.g. $\textit{Pr}=50$ ), viscosity dominates and suppresses heat transport, leading to thinner plumes and reduced thermal fluctuations. Consequently, the equivalence between the thermal fluctuation scales of the inner and outer boundary layers no longer holds. This mechanism is consistent with previous findings on boundary layer asymmetry caused by extreme variations in fluid properties and plume dynamics (Castaing et al. 1989; Wu & Libchaber 1991; Zhang et al. 1997), and explains why the assumption underlying the thermal fluctuation equivalence model ceases to be valid outside the intermediate $\textit{Pr}$ number range.

Figure 5. Ratios of the inner and outer (a) temperature fluctuation scales, $\vartheta ^{i}/\vartheta ^{o}$ , based on (4.4), and (b) boundary layer Rayleigh numbers, $Ra_{\lambda _{\vartheta }}^{i}/Ra_{\lambda _{\vartheta }}^{o}$ , based on (4.9).

Figure 6 presents the compensated plot of thermal boundary layer asymmetry using (4.8). It is evident that within the range of $\textit{Pr}$ considered, the relations derived from the inter-plume distance $Ra$ equivalence model and the boundary layer $Ra$ equivalence model show good agreement with the DNS data only at the high $\textit{Pr}$ value of $50$ .

The reason why these models are valid exclusively in the high- $\textit{Pr}$ regime can be understood through the plume and boundary layer dynamics inherent to high-viscosity fluids. According to Puthenveettil & Arakeri (2005), (4.5) was derived under the assumption that viscous effects dominate, a condition typically satisfied at large $\textit{Pr}$ . At large $\textit{Pr}$ , viscosity stabilises the boundary layers, resulting in thinner and more well-defined plume structures. This ensures that the plume separation distance, or equivalently, the thermal boundary layer thickness, is determined primarily by boundary layer instabilities. Thus the core assumption $Ra_{\lambda _{\vartheta }}^{i} = Ra_{\lambda _{\vartheta }}^{o}$ underlying these models holds robustly at high $\textit{Pr}$ , as illustrated directly in figure 5(b).

Moreover, the inter-plume distance $Ra$ equivalence model was originally formulated for mantle convection scenarios where $\textit{Pr} \rightarrow \infty$ (Jarvis 1993; Vangelov & Jarvis 1994), further highlighting the fundamental dependence of the model assumptions on extremely viscous regimes. Our simulations clearly demonstrate that the equivalence of boundary layer Rayleigh numbers is fulfilled only at $\textit{Pr} = 50$ , confirming that both the inter-plume distance $Ra$ equivalence model and the boundary layer $Ra$ equivalence model are inherently applicable solely in high- $\textit{Pr}$ environments.

Figure 6. Compensated plot of (a) $\vartheta _{mid}$ and (b) $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ , based on (4.8).

5. Characteristic inter-plume distance

In the models discussed in § 4, two among them rely on estimating the characteristic inter-plume distance. The inter-plume density equivalence model assumes that the characteristic inter-plume distance is identical at both the inner and outer boundaries, while the inter-plume distance $Ra$ equivalence model assumes the equivalency of $Ra_{l}$ at the boundaries. Accurate estimation of the characteristic inter-plume distance is therefore essential. In this section, the method for estimating the characteristic inter-plume distance will be discussed. Additionally, the dependence of the inter-plume distance on $\textit{Pr}$ will also be analysed. We also examine the basic assumptions underlying both the inter-plume density equivalence model and the inter-plume distance $Ra$ equivalence model. Finally, for $\textit{Pr}=0.1$ and $50$ , we propose and validate two new models based on the characteristic inter-plume distance.

5.1. Inter-plume distance estimation

To estimate the characteristic plume distance, we first need to delineate and extract the thermal plume from the surrounding turbulent background. Various methods exist for defining a thermal plume. For example, Zhou & Xia (2002) used the local temperature fluctuations $T^{\prime}$ to define a thermal plume. On the other hand, Shishkina & Wagner (2006, 2008) applied the local thermal dissipation rate as a criterion for thermal plume extraction. Additionally, vertical velocity has also been used for this purpose (Julien et al. 1999). Vipin & Puthenveettil (2013) found that the divergence of the horizontal velocity field could also effectively identify thermal plumes, as regions with negative divergence coincide with the areas of positive vertical acceleration in fluid parcels – a key feature of hot plumes. Gastine et al. (2015) compared the plume structures detected by the aforementioned methods, and reported that using the temperature fluctuations yields relatively better results in terms of the shape of the plume contours. We adopt the criterion used by Gastine et al. (2015). Specifically, the plumes are defined as

(5.1) \begin{equation} \left \{ \begin{array}{ll} \displaystyle T - \overline {\left \langle T \right \rangle _{s}} \geqslant \tfrac {1}{2} \sigma _{sd} & \text{at } r=r_{i} + \lambda _{\vartheta }^{i},\\ &\\ \displaystyle T - \overline {\left \langle T \right \rangle _{s}} \leqslant {-\tfrac {1}{2}} \sigma _{sd} & \text{at } r=r_{o} - \lambda _{\vartheta }^{o}, \end{array} \right. \end{equation}

where $\sigma _{sd}$ is the time and horizontally averaged standard deviation of temperature, defined as

(5.2) \begin{equation} \sigma _{sd}(r) = \overline {\sqrt {\big\langle \left ( T- \overline {\left \langle T \right \rangle _{s}} \right)^2 \big\rangle _s}}. \end{equation}

Figure 7. (a) Temperature fluctuation ( $T^{\prime}$ ) contours at the edge of the inner thermal boundary layer ( $r=r_{i}+0.0138$ ) for $\textit{Pr}=50,\ \eta =0.6,\ Ra=5 \times 10^{7}$ . The colour scale ranges from $T^{\prime}=-0.2$ (dark blue) to $T^{\prime}=0.2$ (dark red). (b) Binarised version of (a), following (5.1). The white regions represent plumes, and the black regions represent inter-plume islands.

Figure 7(a) shows the instantaneous temperature field at the edge of the inner thermal boundary layer, with the red regions representing the up-rising hot plumes, and the blue regions indicating the cold inter-plume regions. The plumes connect to form sheet-like structures (Shishkina & Wagner 2008; Stevens et al. 2018). Following (5.1), figure 7(a) is binarised, and the results are shown in figure 7(b). The white regions represent hot plumes, while the black regions represent inter-plume regions. The sheet-like plumes isolate the inter-plume areas into separate ‘islands’.

Puthenveettil & Arakeri (2005) and Zhou & Xia (2010) demonstrated that the statistics of the thermal plumes in turbulent RBC follow a log-normal distribution. Gastine et al. (2015) investigated the statistics of the areas of the inter-plume islands in spherical RBC, and observed a similar distribution. They used the mean plume distances ( $\overline {l}_i,\overline {l}_o$ ) derived from the convective boundary layer model (4.7) to estimate the mean plume separation, and further related it to the characteristic inter-plume areas ( $\overline {A}_{i},\overline {A}_{o}$ ) as

(5.3) \begin{equation} \overline {A}_{i} \sim \frac {\unicode{x03C0} }{4} \overline {l}_{i}^{2}, \quad \overline {A}_{o} \sim \frac {\unicode{x03C0} }{4} \overline {l}_{o}^2, \end{equation}

where the subscripts $i$ and $o$ have their usual meanings. They further used the area corresponding to the peak of the PDF of the inter-plume areas to estimate the characteristic inter-plume area from their simulations. Additionally, they observed that $\unicode{x03C0} \overline {l}_{i}^{2}/4,\ \unicode{x03C0} \overline {l}_{o}^{2}/4$ closely align with the peak values of the PDFs in their simulations, concluding that (5.3) provides a reliable estimate of the characteristic inter-plume area.

The area of an inter-plume island, $A(r)$ , at radial position $r=r_{o}$ is calculated by

(5.4) \begin{equation} A(r=r_{o}) = r_{o}^2 \iint _{\varSigma } \sin {\theta }\, \mathrm{d}\theta\, \mathrm{d}\phi, \end{equation}

where $\varSigma$ represents the boundary of the inter-plume island, and $\theta$ and $\phi$ are the co-latitude and longitude, respectively. We identify all inter-plume islands located at the radial positions $r=r_i + \lambda _{\vartheta }^{i}$ and $r = r_{o} - \lambda _{\vartheta }^{o}$ , and calculate the area of each island using (5.4). The PDFs of the inter-plume areas are averaged over a total time interval of at least 100 convective time units. The resulting PDFs are consistent with those of Gastine et al. (2015).

In the present work, we define the area corresponding to the peak of the PDF as the characteristic inter-plume area. Using $\mathcal{A}_{i}$ and $\mathcal{A}_{o}$ , the characteristic inter-plume areas at the edges of the inner and outer thermal boundary layers, we define characteristic inter-plume distances $L_{i}$ and $L_{o}$ as

(5.5) \begin{equation} \mathcal{A}_{i} = \frac {\unicode{x03C0} }{4} L_{i}^{2}, \quad \mathcal{A}_{o} = \frac {\unicode{x03C0} }{4} L_{o}^{2}. \end{equation}

Gastine et al. (2015) demonstrated that the characteristic inter-plume areas are equivalent at the inner and outer thermal boundary layers at $\textit{Pr}=1$ . They developed a model to quantify thermal boundary layer asymmetry based on this observation, as introduced in § 4.1. However, as already discussed, we observe significant deviations from this assumption when $\textit{Pr}$ is not unity. In the next subsection, we examine how the inter-plume area varies with $\textit{Pr}$ .

5.2. Inter-plume area: $\textit{Pr}$ dependence

In this subsection, we explore the $\textit{Pr}$ dependence of the characteristic inter-plume areas at the inner and outer thermal boundary layers.

Figure 8. The PDFs of the inter-plume area for (a) the inner thermal boundary layer, and (b) the outer thermal boundary layer, at $\eta =0.2,\ Ra=10^{7}$ .

Figure 8 presents the time-averaged PDFs of the inter-plume areas for both thermal boundary layers at $\eta =0.2,\ Ra=10^{7}$ across different $\textit{Pr}$ values. As shown in figure 8(a), the characteristic inter-plume area for the inner thermal boundary layer remains unchanged for $0.1 \leqslant Pr \leqslant 50$ . In contrast, the characteristic inter-plume area for the outer thermal boundary layer increases with $\textit{Pr}$ . This difference highlights that the $\textit{Pr}$ dependence of the characteristic inter-plume area varies between the inner and outer boundary layers.

Figure 9 provides qualitative insight into this phenomenon. The instantaneous contour plots of the binarised temperature fluctuations near the inner and outer boundary layers at $\eta =0.2$ for different $\textit{Pr}$ are shown. The morphology and distribution of inter-plume islands change significantly with increasing $\textit{Pr}$ ; they become more fragmented at both boundary layers. However, this fragmentation is more pronounced near the outer boundary layer. For example, at $\textit{Pr}=0.1$ , a ‘super-island’ exists near the outer boundary layer, where the thermal plumes (white regions in the binarised plot) are thick, leaving the inter-plume areas scattered and small. As $\textit{Pr}$ increases, smaller islands emerge and the plumes become thinner, leading to a larger number of inter-plume islands. The plumes become thinner also at the inner boundary layer with increasing $\textit{Pr}$ , but the characteristic inter-plume area remains constant. Moreover, the distribution of inter-plume islands becomes increasingly similar between the inner and outer thermal boundary layers as $\textit{Pr}$ increases.

Figure 9. Contour plots of the binarised temperature fluctuations for (a,c,e) the inner thermal boundary layer, and (b,d,f) the outer thermal boundary layer, at $\eta =0.2,\ Ra=10^{7}$ , for (a,b) $\textit{Pr}=0.1$ , (c,d) $\textit{Pr}=1$ , (e,f) $\textit{Pr}=50$ .

Figure 10 shows the time-averaged PDFs of the inter-plume areas near the inner and outer thermal boundary layers for $\eta =0.6,\ Ra=10^{7}$ at different $\textit{Pr}$ . The PDFs for both the inner and outer thermal boundary layers exhibit a similar trend in their $\textit{Pr}$ dependence, with the characteristic inter-plume area increasing as $\textit{Pr}$ increases. However, while the $\textit{Pr}$ dependence of the characteristic inter-plume area is comparable between the two boundary layers, the increase is more pronounced near the outer thermal boundary layer than near the inner one. As $\textit{Pr}$ increases, the thermal plumes become thinner, resulting in a larger characteristic inter-plume area. At $\textit{Pr}=50$ , the morphology and the distribution of the inter-plume islands become notably similar between the inner and outer thermal boundary layers, as shown in figure 11.

Figure 10. The PDFs of the inter-plume area for (a) the inner thermal boundary layer, and (b) the outer thermal boundary layer, at $\eta =0.6,\ Ra=10^{7}$ .

Figure 11. Contour plots of the binarised temperature fluctuations for (a,c,e) the inner thermal boundary layer, and (b,d,f) the outer thermal boundary layer, at $\eta =0.6,\ Ra=10^{7}$ , for (a,b) $\textit{Pr}=0.1$ , (c,d) $\textit{Pr}=1$ , (e,f) $\textit{Pr}=50$ .

The $\textit{Pr}$ dependence of the inter-plume area also helps to explain our earlier observation in § 3 that the asymmetry exhibits a weaker dependence on $\textit{Pr}$ at higher values of $\eta$ . Qualitatively, as $\eta$ approaches unity, the boundary layers become increasingly symmetric, significantly reducing asymmetries induced by variations in $\textit{Pr}$ . Furthermore, plume morphology and the statistical analysis of characteristic inter-plume areas also support this understanding. Specifically, at smaller radius ratios (e.g. $\eta =0.2$ ), radial differences in curvature and gravity ( $g \sim 1/r^2$ ) between the inner and outer boundaries are pronounced, causing significant variations in plume dynamics as $\textit{Pr}$ changes. In contrast, at higher radius ratios (e.g. $\eta =0.6$ ), these radial differences diminish, resulting in more symmetric boundary layer conditions and consequently a weaker sensitivity of the asymmetry to changes in $\textit{Pr}$ .

We now examine the basic assumption of the inter-plume density equivalence model, which posits that the characteristic plume separations are identical near the inner and outer thermal boundary layers. In Gastine et al. (2015), this equivalence of characteristic plume separation was reformulated as the equivalence of characteristic inter-plume area. Based on the PDF of the inter-plume areas, they demonstrated that the characteristic inter-plume area is nearly identical near the inner and outer thermal boundary layers for $\textit{Pr}=1$ . Following their methodology, we estimate the variation in characteristic plume separation by analysing the variation in the characteristic inter-plume area.

Figure 12 presents the time-averaged PDFs of the inter-plume area near the inner and outer thermal boundary layers for $\eta =0.2,\ Ra=10^{7}$ at various $\textit{Pr}$ . As shown in the figure, the characteristic inter-plume area for the outer thermal boundary layer is smaller than that for the inner boundary layer at $\textit{Pr}=0.1$ . Moreover, the characteristic inter-plume area for the outer thermal boundary layer increases with $\textit{Pr}$ , whereas it remains unchanged for the inner boundary layer, as previously illustrated in figure 8. The characteristic inter-plume area for the outer thermal boundary layer surpasses that of the inner boundary layer at approximately $\textit{Pr}=0.5$ , and becomes significantly larger for $\textit{Pr} \geqslant 5$ .

Figure 13 shows the corresponding PDFs for $\eta =0.6,\ Ra=10^{7}$ . Similar to the $\eta =0.2$ case, the characteristic inter-plume area for the outer thermal boundary layer is initially smaller than that for the inner boundary layer, at $\textit{Pr}=0.1$ . It surpasses the characteristic inter-plume area of the inner boundary layer at approximately $\textit{Pr}=5$ (though this specific case is not explicitly shown in figure 13). Furthermore, as the asymmetry decreases with increasing $\eta$ , PDFs for the inner and outer thermal boundary layers are similar at $\eta =0.6$ in comparison to the PDFs at $\eta =0.2$ .

From the above discussion, it follows that the characteristic inter-plume areas near the inner and outer boundary layers are not equal for either $\textit{Pr}=0.1$ or $5 \leqslant Pr \leqslant 50$ . Consequently, $0.2 \lesssim Pr \lesssim 1$ is identified as the appropriate validation range for the basic assumption of the inter-plume density equivalence model, which is in line with the range determined by comparing the thermal boundary layer asymmetry derived from this model, as discussed in § 4.5.

Figure 12. Time-averaged PDFs of the inter-plume areas for the inner (blue) and outer (green) thermal boundary layers at $\eta =0.2,\ Ra=10^{7}$ , at different $\textit{Pr}$ values: (a) $\textit{Pr}=0.1$ , (b) $\textit{Pr}=0.2$ , (c) $\textit{Pr}=0.5$ , (d) $\textit{Pr}=1$ , (e) $\textit{Pr}=10$ , (f) $\textit{Pr}=50$ .

Figure 13. Time-averaged PDFs of the inter-plume areas for the inner (blue) and outer (green) thermal boundary layers at $\eta =0.6,\ Ra=10^{7}$ , at different $\textit{Pr}$ values: (a) $\textit{Pr}=0.1$ , (b) $\textit{Pr}=0.2$ , (c) $\textit{Pr}=0.5$ , (d) $\textit{Pr}=1$ , (e) $\textit{Pr}=10$ , (f) $\textit{Pr}=50$ .

5.3. Proposed model for small $\textit{Pr}$

As demonstrated in the previous subsection, the PDFs of the inter-plume areas exhibit a dependence on $\textit{Pr}$ . At $\textit{Pr}=0.1$ , the characteristic inter-plume area near the inner thermal boundary layer ( $\mathcal{A}_{i}$ ) is larger than that near the outer thermal boundary layer ( $\mathcal{A}_{o}$ ). By analysing the PDFs of the inter-plume areas near both boundaries, we identify the following relation between $\mathcal{A}_{i}$ and $\mathcal{A}_{o}$ :

(5.6) \begin{equation} \mathcal{A}_{i} \lambda _{\vartheta }^{i} = \mathcal{A}_{o} \lambda _{\vartheta }^{o}, \end{equation}

where $\lambda _{\vartheta }^{i}$ and $\lambda _{\vartheta }^{o}$ represent the thicknesses of the inner and outer thermal boundary layers, respectively.

This model can be regarded as an extension of the inter-plume density equivalence model, which is based on the assumption that the characteristic inter-plume areas are identical near the inner and outer boundaries. This assumption can be interpreted as stating that the average surface area occupied by a single plume is the same at both boundaries. However, we find that this assumption does not hold at $\textit{Pr} = 0.1$ . Instead, (5.6) holds, which can be interpreted as the average volume occupied by a plume being equal between the two boundaries. This revised assumption shares a similar physical picture with the original inter-plume density equivalence model, but is better suited to describing the flow characteristics near $\textit{Pr}=0.1$ .

Figures 14 and 15 display the time-averaged PDFs of the inter-plume areas on a logarithmic scale ( $\log_{10}A$ ) and the PDFs of the inter-plume area multiplied by the corresponding thermal boundary layer thickness ( $\log_{10}(\lambda _{\vartheta }A)$ ). As shown in figures 14(b) and 15(b), the peak values of the PDFs for the inner and outer thermal boundary layers align, validating (5.6), as both $\lambda _{\vartheta }^{i}$ and $\lambda _{\vartheta }^{o}$ are constant.

Figure 14. Time-averaged PDFs of the inter-plume areas for the inner (blue) and outer (green) thermal boundary layers at $\eta =0.2,\ Ra=10^{7}$ : (a) original PDF and (b) scaled PDF. In (b), the blue symbols represent the scaled PDF of $\lambda _{\vartheta }^{i} A_{i}$ , and the green symbols denote the scaled PDF of $\lambda _{\vartheta }^{o} A_{o}$ .

Figure 15. Time-averaged PDFs of the inter-plume areas for the inner (blue) and outer (green) thermal boundary layers at $\eta =0.6,\ Ra=10^{7}$ : (a) original PDF and (b) scaled PDF. In (b), the blue symbols represent the scaled PDF of $\lambda _{\vartheta }^{i} A_{i}$ , and the green symbols denote the scaled PDF of $\lambda _{\vartheta }^{o} A_{o}$ .

Figure 16. Compensated plot of (a) $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ , and (b) $\vartheta _{mid}$ , based on (5.9).

Using the mean plume density defined in (4.7) to estimate the characteristic inter-plume area, we have

(5.7) \begin{equation} \mathcal{A}_{i} \approx \frac {\unicode{x03C0} }{4} \overline {l}_{i}^{2}, \quad \mathcal{A}_{o} \approx \frac {\unicode{x03C0} }{4} \overline {l}_{o}^{2}. \end{equation}

Substituting this into (5.6) gives

(5.8) \begin{equation} \overline {l}_{i}^{2} \lambda _{\vartheta }^{i} = \overline {l}_{o}^{2} \lambda _{\vartheta }^{o} \quad \Rightarrow \quad \frac {\alpha g_{i}\, \Delta \vartheta _{i}\, (\lambda _{\vartheta }^{i})^{6}}{\nu \kappa } = \frac {\alpha g_{o}\, \Delta \vartheta _{o}\, (\lambda _{\vartheta }^{o})^{6}}{\nu \kappa }, \end{equation}

which provides a new relation for deriving the thermal boundary layer asymmetry. Combining this with (3.1) and (3.2) yields

(5.9) \begin{equation} \Delta \vartheta _{i} = \frac {1}{1+\eta ^{10/7}}, \quad \Delta \vartheta _{o}= \vartheta _{mid} = \frac {\eta ^{10/7}}{1+\eta ^{10/7}}, \quad \frac {\lambda _{\vartheta }^{i}}{\lambda _{\vartheta }^{o}}=\eta ^{4/7}. \end{equation}

Figure 16 demonstrates the validity of (5.9) at the smallest Prandtl number ( $\textit{Pr}=0.1$ ) considered in our simulations. The results presented in figure 16 demonstrate excellent agreement with the scaling relation given by (5.9) for both $\eta =0.2$ and $\eta =0.6$ . Given that our model shares a similar physical picture with the inter-plume density equivalence model, and considering that the latter has already been validated over a broad range of $\eta$ values in the studies by Gastine et al. (2015) and Fu et al. (2024), we anticipate that our proposed model can be reasonably extrapolated to other $\eta$ values as well. Furthermore, our model is consistent in the limit $\eta \rightarrow 1$ , as it correctly predicts symmetric boundary layer temperature drops ( $\vartheta _{mid} = 1/2$ ) and equal boundary layer thicknesses ( $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}=1$ ), thus reinforcing its robustness and physical reliability.

5.4. Proposed model for large $\textit{Pr}$

In § 5.2, we demonstrated that the morphology and distribution of inter-plume areas become increasingly similar between the inner and outer thermal boundary layers as $\textit{Pr}$ increases. For instance, at $\textit{Pr}=50,\ \eta =0.6$ , the morphologies of inter-plume islands near the inner and outer thermal boundary layers appear quite similar, as shown in figures 11(e,f). However, despite the similar appearance of the Hammer projections at high $\textit{Pr}$ , the physical scales remain different. The total surface area of the edge of the inner thermal boundary layer is approximately $(r_{i}/r_{o})^{2}$ times smaller than that of the outer thermal boundary layer. Motivated by this observation, we propose a new assumption for the large $\textit{Pr}$ regime. To begin, we define the normalised inter-plume area $\tilde {A}(r)$ at a radial location $r=r_{o}$ as

(5.10) \begin{equation} \tilde {A}(r=r_{o}) = \frac {A(r=r_{o})}{r_{o}^{2}} = \iint _{\varSigma } \sin\theta\mathrm{d}\theta\mathrm{d}\phi, \end{equation}

where $\varSigma$ represents the boundary of an inter-plume island. We compute the PDF of the normalised inter-plume area, and define the characteristic normalised inter-plume area $\tilde {\mathcal{A}}$ as the value corresponding to the peak value of this PDF.

Our assumption is that as $\textit{Pr}$ increases, the characteristic normalised inter-plume areas near the inner and outer thermal boundary layers converge and eventually become identical at sufficiently large $\textit{Pr}$ . This is from the observation that the inter-plume islands become increasingly similar in shape at the inner and outer boundaries, which then, upon normalising with the surface area ( $4 \unicode{x03C0} r^2$ ) term, gives equivalent PDF peaks. Figures 17 and 18 display the PDFs of the normalised inter-plume areas at $Ra=10^{7}$ for various $\textit{Pr}$ values, with $\eta =0.2$ and $\eta =0.6$ , respectively. It is evident that the characteristic normalised inter-plume areas near the inner and outer thermal boundary layers approach each other as $\textit{Pr}$ increases, and become nearly identical at $\textit{Pr}=50$ .

Figure 17. The PDFs of the normalised inter-plume areas from (5.10) at $Ra=10^{7}$ , $\eta =0.2$ , for (a) $\textit{Pr}=0.1$ , (b) $\textit{Pr}=1$ , (c) $\textit{Pr}=10$ , (d) $\textit{Pr}=50$ .

Figure 18. The PDFs of the normalised inter-plume areas from (5.10) at $Ra=10^{7}$ , $\eta =0.6$ , for (a) $\textit{Pr}=0.1$ , (b) $\textit{Pr}=1$ , (c) $\textit{Pr}=10$ , (d) $\textit{Pr}=50$ .

Based on the above finding, we establish the following relation at $\textit{Pr}=50$ :

(5.11) \begin{equation} \tilde {\mathcal{A}}_{i} = \tilde {\mathcal{A}}_{o} \quad \Rightarrow \quad \frac {\mathcal{A}_{i}}{r_{i}^{2}} = \frac {\mathcal{A}_{o}}{r_{o}^{2}}. \end{equation}

Substituting (4.7) and (5.7) into (5.11), we can write

(5.12) \begin{equation} \frac {\mathcal{A}_{i}}{\mathcal{A}_{o}}=\frac {\overline {l}_{i}^{2}}{\overline {l}_{o}^{2}} = \eta ^{2} \quad \Rightarrow \quad \frac {\alpha g_{i}\, \Delta \vartheta _{i}\, (\lambda _{\vartheta }^{i})^{5}}{\nu \kappa } = \frac {\alpha g_{o}\, \Delta \vartheta _{o}\, (\lambda _{\vartheta }^{o})^{5}}{\nu \kappa } \eta ^{2}. \end{equation}

Equation (5.12) serves as the third relation for describing thermal boundary layer asymmetry. Combining this equation with (3.1) and (3.2) leads to the same results as in (4.8). As discussed in § 4.5, (4.8) provides an accurate description at $\textit{Pr}=50$ . Furthermore, our model for the high- $\textit{Pr}$ regime also remains consistent in the limiting case as $\eta \rightarrow 1$ . In this planar limit, the temperature drops across the inner and outer boundary layers become symmetric, satisfying $\vartheta _{mid} = 1/2$ , and the boundary layer thickness ratio approaches unity, i.e. $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}=1$ . These conditions align precisely with the expected symmetry and known results from planar RBC. Thus this consistency further supports the robustness and physical validity of our proposed high- $\textit{Pr}$ model.

The details of all thermal boundary layer asymmetry models examined in this study are summarised in table 2. Additionally, the assumptions underlying each model and their respective validation ranges in $\textit{Pr}$ space are also listed.

Table 2. Comparison of thermal boundary layer asymmetry models analysed in this study. The parameter $\gamma$ represents the exponent in the thermal boundary layer asymmetry equation (4.10).

6. Conclusion and outlook

In this study, we investigated the thermal boundary layer asymmetry of turbulent Rayleigh–Bénard convection (RBC) in spherical shells across a range of Prandtl numbers ( $0.1 \leqslant Pr \leqslant 50$ ) and radius ratios ( $\eta = 0.2, 0.6$ ) using direct numerical simulations (DNS). A centrally condensed mass with gravity profile $g(r)=(r_{o}/r)^{2}$ was considered. Two key parameters were employed to quantify the asymmetry of the thermal boundary layers: the time and horizontally averaged bulk temperature $\vartheta _{mid}$ , and the thermal boundary layer thickness ratio $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$ . The dependence of these parameters on $\eta$ was investigated for different $\textit{Pr}$ regimes.

We analysed four widely used thermal boundary layer asymmetry models using our DNS data, and examined their underlying assumptions. It was found that the inter-plume density equivalence model and the thermal fluctuation equivalence model are valid for $0.2 \lesssim Pr \lesssim 1$ , whereas the inter-plume distance $Ra$ equivalence model and the boundary layer $Ra$ equivalence model are valid only at $\textit{Pr}=50$ . These findings demonstrate that each model has a distinct validation range in $\textit{Pr}$ space.

We examined the $\textit{Pr}$ dependence of the PDF of the inter-plume area. The $\textit{Pr}$ dependence of the characteristic inter-plume area differed between the two boundaries. For $\eta =0.2$ , the characteristic inter-plume area increased with $\textit{Pr}$ for the outer thermal boundary layer, while it remained unchanged for the inner thermal boundary layer within the entire range $0.1 \leqslant Pr \leqslant 50$ . For $\eta =0.6$ , the characteristic inter-plume area increased for both boundary layers, but the increase was more pronounced for the outer boundary layer. At $\textit{Pr}=0.1$ , the characteristic inter-plume area was larger for the inner boundary layer than the outer one, but this trend reversed for $5 \leqslant Pr \leqslant 50$ .

By analysing our simulation data, we found that the relation $\mathcal{A}_{i} \lambda _{\vartheta }^{i} = \mathcal{A}_{o} \lambda _{\vartheta }^{o}$ holds at $\textit{Pr}=0.1$ . Based on this relation, a new model for the thermal boundary layer asymmetry at $\textit{Pr}=0.1$ was proposed and validated against DNS data. Furthermore, we observed that the morphologies of inter-plume areas become increasingly similar as $\textit{Pr}$ increases. Building upon this observation, we proposed that $\mathcal{A}_{i}/\mathcal{A}_{o}=\eta ^{2}$ holds for large values of $\textit{Pr}$ . This relation was validated at $\textit{Pr}=50$ . The proposed model yielded results consistent with the inter-plume density equivalence model and the boundary layer $Ra$ equivalence model, and showed excellent agreement with DNS data at $\textit{Pr}=50$ .

Table 3. Simulation details: $Ra$ is the Rayleigh number, $\textit{Nu}$ is the Nusselt number, $Err(\textit{Nu})$ represents the relative error of the four Nusselt numbers, as described in § 2.3, $Re$ denotes the global Reynolds number, $\vartheta _{mid}$ is the mean temperature at mid-depth, $\lambda _{\vartheta }^{i}$ and $\lambda _{\vartheta }^{o}$ are the inner and outer thermal boundary layer thicknesses, respectively, $\lambda _{u}^{i}$ and $\lambda _{u}^{o}$ correspond to the inner and outer viscous boundary layer thicknesses, $\epsilon _{\vartheta }^{bulk}$ represents the proportion of the bulk contribution to the thermal energy dissipation rate, $\epsilon _{u}^{bulk}$ denotes the proportion of the bulk contribution to the kinetic energy dissipation rate, $N_{\lambda _{\vartheta }}^{i}$ and $N_{\lambda _{\vartheta }}^{o}$ are the numbers of grid points within the inner and outer thermal boundary layers, respectively, and $N_{\lambda _{u}}^{i}$ and $N_{\lambda _{u}}^{o}$ represent the numbers of grid points within the inner and outer viscous boundary layers, respectively.

It is important to emphasise that all our analyses and conclusions are obtained specifically for the centrally condensed mass with the gravity profile $g(r)\sim 1/r^{2}$ , which is representative of convection within spherical shells characteristic of planetary atmospheres, such as gas giants. The applicability of our results to other gravity profiles, such as the linear profile $g(r)\sim r$ relevant to terrestrial mantle convection, is not straightforward. Indeed, different gravity distributions will alter buoyancy forcing and boundary layer instabilities significantly, particularly at low $\textit{Pr}$ where thermal diffusivity dominates the plume dynamics. Thus further targeted studies are necessary to quantitatively evaluate the robustness of our proposed scaling relations under alternative gravity conditions.

Moreover, an important aspect for future research is to examine the validity and robustness of the proposed models under more extreme parameter regimes. Specifically, extending our models to significantly higher Rayleigh numbers would be valuable, as the current scaling relationships are based on local balances and statistical properties verified primarily within the range $5 \times 10^6 \leqslant Ra \leqslant 5 \times 10^7$ . While our results demonstrate consistency across the considered $Ra$ range, some deviations might arise at much higher $Ra$ , potentially due to transitions to different turbulent regimes. Additionally, investigating the applicability of our models at smaller radius ratios ( $\eta \ll 1$ ) remains an open question. Such small- $\eta$ scenarios introduce increasingly pronounced curvature effects, which may fundamentally alter plume dynamics and boundary layer structures, potentially challenging the underlying assumptions of our models. Exploring these regimes would therefore provide valuable insights into the broader applicability and limitations of our theoretical framework.

Figure 19. Schematic of the computational domain and coordinate system.

Figure 20. Time series of the volume-averaged buoyancy power $\mathcal{P}_{buoy}$ , the volume-averaged kinetic energy dissipation rate $\epsilon _{u}$ , the total energy balance $\mathcal{P}_{buoy} + \epsilon _{u}$ , and the volume-averaged kinetic energy variation rate ${\rm d}E_{{kin}}/{\rm d}t$ at (a) $\eta =0.2,\ Ra=5 \times 10^{7}$ and (b) $\eta =0.6,\ Ra=5 \times 10^{7}$ . For both cases, $\textit{Pr}=50$ .