## 1. Introduction

The effect of an off-centre gravity location on the flow structures in turbulent thermal convection has not been studied previously. Here, we use spherical Rayleigh–Bénard (RB) convection as an idealized model system to study the effect of an off-centre gravity location on thermal convection in a spherical system. The system consists of a fluid layer enclosed between two spherical shells, heated from the inner sphere and cooled from the outer one. Convection in spherical shells differs from the classical RB convection (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012) owing to the curvature of the boundaries, non-uniform gravity, and the geometrical asymmetry between the boundary layer at the inner and outer spheres (Busse Reference Busse1970; Spiegel Reference Spiegel1971; Jarvis, Glatzmaierand & Vangelov Reference Jarvis, Glatzmaierand and Vangelov1995; Tilgner Reference Tilgner1996; Shahnas *et al.* Reference Shahnas, Lowman, Jarvis and Bunge2008; Deschamps, Tackley & Nakagawa Reference Deschamps, Tackley and Nakagawa2010; O'Farrell, Lowman & Bunge Reference O'Farrell, Lowman and Bunge2013; Gastine, Wicht & Aurnou Reference Gastine, Wicht and Aurnou2015).

Most studies on spherical RB convection focus on the case in which the gravity centre coincides with the geometric one. Busse (Reference Busse1975) used perturbation analysis to show the qualitative difference between convective flow patterns of odd and even spherical harmonic order just above the onset of convection. Subsequently, in an infinite Prandtl number ($Pr$) fluid with application to mantle convection, Bercovici, Schubert & Glatzmaier (Reference Bercovici, Schubert and Glatzmaier1989) showed that these convective patterns persist for Rayleigh numbers ($Ra$) up to 100 times larger than the critical value. For $Ra>10^5$, the axisymmetric convective patterns break down when the flow starts to show time-dependent behaviour (Iwase & Honda Reference Iwase and Honda1997; Deschamps *et al.* Reference Deschamps, Tackley and Nakagawa2010; Futterer *et al.* Reference Futterer, Krebs, Plesa, Zaussinger, Hollerbach, Breuer and Egbers2013; Yanagisawa, Kameyama & Ogawa Reference Yanagisawa, Kameyama and Ogawa2016). In spherical RB convection, the influence of curvature is important. For constant $Ra$, the length scale of the preferentially convective patterns decreases with increasing radius ratio (Deschamps *et al.* Reference Deschamps, Tackley and Nakagawa2010). Yanagisawa & Yamagishi (Reference Yanagisawa and Yamagishi2005) find that the convective rolls become smaller with increasing $Ra$ when the radius ratio is kept constant. Gastine *et al.* (Reference Gastine, Wicht and Aurnou2015) show that for $Pr=1$, sheet-like thermal plumes are formed near the inner and outer spheres, and these plumes undergo morphological changes into mushroom-like plumes when they eject from the boundary layers to the bulk. These plume dynamics are similar to the situation for classic RB convection (Puthenveettil & Arakeri Reference Puthenveettil and Arakeri2005; Shishkina & Wagner Reference Shishkina and Wagner2008; Zhou & Xia Reference Zhou and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012). However, due to the asymmetries between the hot and cold surfaces in spherical shells, the mushroom-like plumes emitted from the outer sphere are thicker than those emitted from the inner sphere (Gastine *et al.* Reference Gastine, Wicht and Aurnou2015). Nevertheless, studies have shown that there are similarities between spherical and classical planar RB convection. For example, Gastine *et al.* (Reference Gastine, Wicht and Aurnou2015) showed that the effective heat transport (characterized by the Nusselt number $Nu$) scaling exponent $\alpha$ in $Nu \sim Ra^\alpha$ in spherical RB convection has a similar $Ra$ number dependence as in planar RB convection, where for $Ra\lesssim 10^{11}$, typically $0.28 \lesssim \alpha \lesssim 0.31$ (Ahlers *et al.* Reference Ahlers, Grossmann and Lohse2009).

However, the misalignment between the gravity centre and the geometric one in the spherical RB convection set-up has been overlooked to date. In this study, we use three-dimensional direct numerical simulations to investigate the effect of the gravity centre offset on the flow structures and heat transfer in spherical RB convection. The paper is organized as follows. In § 2, we explain the numerical method and discuss the considered parameter regime. In §§ 3.1–3.3, we discuss the large-scale flow pattern induced by the off-centre gravity. The effect of the large-scale structure on heat transfer is shown in § 3.4. In § 3.5, we study the effect of the gravity profile and $Ra$ on flow structures and heat transfer. The paper ends with a summary of our findings and conclusions in § 4.

## 2. Numerical method and parameters

### 2.1. Set-up of convection in spherical shells

The spherical RB geometry is illustrated schematically in figure 1. Fluid fills a spherical shell between the inner sphere of radius $r_i$ and the outer sphere of radius $r_o$. The radius ratio between the inner and outer spheres is given by

In this study, we consider a fixed aspect ratio $\eta =0.3$. The surface temperatures of the inner and outer spheres are kept constant at $T_i$ and $T_o$, respectively, with $T_i > T_o$. No-slip boundary conditions are imposed on both boundaries. Point $G$ indicates the gravity centre, which is offset from the geometric centre $O$. The gravity shift is defined as

where $|\boldsymbol {OG}|$ indicates the displacement of the gravity compared to the geometrical centre; see figure 1. In the absence of rotation (and magnetic field), for a sphere, there are no north and south directions, therefore any diametrically opposite symmetry is equivalent to another. Nevertheless, since the spherical coordinates are naturally clustered around the polar axis, it is beneficial to consider the jet aligned along this axis, referred to as north–south, from a computation point of view. For convenience, we always shift the gravity centre $G$ to the south.

The dynamics of spherical RB convection is controlled by the Rayleigh and Prandtl numbers

*a*,

*b*)\begin{equation} Ra=\frac{\beta g_o\,\Delta T\,d^3}{\kappa \nu}, \quad Pr=\frac{\nu}{\kappa}, \end{equation}

where $\beta$ is the thermal expansion coefficient, $g_o$ is the surface-averaged gravity at the outer sphere, $\nu$ is the kinematic viscosity, and $\kappa$ is the thermal diffusivity. We normalize the fields and distances by the length scale $d=r_o-r_i$, the temperature difference $\Delta T$ between the inner and outer spheres, and the free-fall velocity $U=\sqrt {\beta g_o\,\Delta T\,d}$.

### 2.2. Underlying dynamical equation and numerical method

We solve the Navier–Stokes equations under the Boussinesq approximation in spherical coordinates, which in dimensionless form read

where $\boldsymbol {u}$, $p$, $T$ and $\boldsymbol {C}_{g}g$ denote the fluid velocity, pressure, temperature and gravitational acceleration acting in the three directions. The coefficients $C_{g,i}$ (where $i=1,2,3$ correspond to the $\hat {\theta }$, $\hat {r}$ and $\hat {\varphi }$ directions, respectively) denote the decomposition of $g$, which is detailed in Appendix A. Since we always shift the gravity centre $G$ to the south, the three components of the buoyancy coefficient can be written as

in which $L_{GP}=|\boldsymbol {GP}|$ indicates the distance between a fluid point and the gravity centre; see figure 1. We consider the gravity profile

and we use $n_g=-2$, $-1$, $0$ and $1$. The value of $n_g$ can be inferred by physical considerations on the mass ratio between the nucleus ($r< r_i$) and the spherical shell ($r_i< r< r_o$). The gravity is constant ($n_g=0$) when the mass of the spherical shell is negligible. This assumption is typically used to model the Earth's mantle (Bercovici *et al.* Reference Bercovici, Schubert and Glatzmaier1989). If the mass is condensed centrally, then one obtains $n_g=-2$. Gastine *et al.* (Reference Gastine, Wicht and Aurnou2015) showed that only for this gravity distribution function is there an analytical relation between the viscous dissipation rate ($\epsilon _U$) and $Nu$, i.e. $\epsilon _U = ({3}/({1+\eta +\eta ^2})) ({1}/{Pr})(Nu-1)$. When the density is constant within the fluid layer, and no mass is contained in the inner sphere, gravity is directly proportional to the radial coordinate $r$ ($n_g=1$) (Tilgner Reference Tilgner1996).

The governing equations (2.4) and (2.5) are discretized by a staggered central second-order finite-difference scheme in spherical coordinates. The numerical scheme is based on the method by Verzicco & Orlandi (Reference Verzicco and Orlandi1996), which has been extended recently to spherical coordinates by Santelli, Orlandi & Verzicco (Reference Santelli, Orlandi and Verzicco2021). The advantage of this scheme is that it allows arbitrary non-uniform grids in the radial and co-latitudinal directions. Here, it is worthwhile to mention that singularities at the poles are prevented by introducing a new set of quantities $(u_\theta,u_r r^2, u_\varphi \sin \varphi )$, which results in trivial boundary conditions $u_\varphi \sin \varphi =0$ at the north pole ($\varphi =0$) and south pole ($\varphi ={\rm \pi}$) in the co-latitudinal direction. The code has been validated carefully by Santelli *et al.* (Reference Santelli, Orlandi and Verzicco2021), and we refer to that work for details on the method. For additional validation results relevant to the flows considered here, we refer the reader to Appendix B.

### 2.3. Choice of parameters

To study the effect of the gravity centre location on the flow dynamics in spherical RB convection, we considered different gravity distributions ($n_g \in \{-2,-1,0,1\}$, see (2.7)) and gravity centre offset ($\varepsilon$ from $0$ to $0.8$, see (2.2)) for $Ra=10^7$ and $10^8$. All simulations in this study are for $Pr=1$. To ensure that the flow is resolved fully, we place a sufficient number of computational grid points in the bulk and the boundary layers (Stevens, Verzicco & Lohse Reference Stevens, Verzicco and Lohse2010). As we will use spectral analysis to study the flow structures, we use a uniform grid in the longitudinal and co-latitudinal directions. In the radial direction, the grid cells are clustered towards the inner and outer spheres to ensure that the boundary layers are resolved adequately (Shishkina *et al.* Reference Shishkina, Stevens, Grossmann and Lohse2010).

We calculate the Nusselt number $Nu$ from the normalized averaged temperature gradients at the inner and outer spheres as follows:

*a*,

*b*)\begin{equation} Nu_i = \left. -\eta\,\frac{{\rm d} \overline{\langle T \rangle_{\theta,\varphi}}}{{\rm d}r} \right| _{r_i},\quad Nu_o= \left. -\frac{1}{\eta}\,\frac{{\rm d} \overline{\langle T \rangle_{\theta,\varphi}}}{{\rm d}r} \right| _{r_o}. \end{equation}

Here, $\langle \cdots \rangle _{\theta,\varphi }$ represents the average over a spherical surface with fixed $r$, and $\overline {\cdots }$ indicates a time average. Also, $Nu_i$ and $Nu_o$ as functions of the co-latitudinal direction are defined as

*a*,

*b*)\begin{equation} Nu_i(\varphi) = \left. -\eta\,\frac{{\rm d} \overline{\langle T \rangle_\theta}}{{\rm d}r} \right| _{r_i},\quad Nu_o(\varphi)= \left. -\frac{1}{\eta}\,\frac{{\rm d} \overline{\langle T \rangle_\theta}}{{\rm d}r} \right| _{r_o}, \end{equation}

where $\langle \cdots \rangle _{\theta }$ represents the average over the azimuthal direction. The difference in $Nu$ values obtained at the inner and outer spheres is always less than $0.2\,\%$. Also, we verify that $Nu$ calculated at the spheres is the same as the value obtained from the volume-averaged Nusselt number $Nu_v$, which is calculated from the surface-averaged Nusselt number $Nu_h$ averaged over concentric spherical surfaces with fixed $r$, with

where $T_c(r)=\eta /[(1-\eta )^2r]-\eta /(1-\eta )$ is the purely conductive temperature profile for constant temperature boundary conditions. Subsequently, we perform radial integration and time averaging to obtain $Nu_v$ as

The volume-averaged root-mean-square (r.m.s.) Reynolds number is given by

Here, $\langle {\cdots }\rangle$ denotes the average over the whole volume of the spherical shell. The details of the simulations considered in this study are summarized in table 1.

## 3. Results

Before we discuss the turbulent flows in §§ 3.2–3.5, we first focus on the large-scale circulation formation induced by off-centre gravity in a steady axisymmetric flow obtained at $Ra=10^3$. As indicated in (2.6), the radial and co-latitudinal buoyancy coefficients depend on the gravity centre offset $\varepsilon$ and the gravity distribution $n_g$. To disentangle these effects, for turbulent flows we first examine the effect of the gravity centre offset for constant gravity ($n_g=0$) at $Ra=10^8$ in §§ 3.2–3.4. The effect of $n_g$ and $Ra$ is studied in § 3.5.

### 3.1. Dynamics of steady flow

When the gravity centre offsets with the geometric one in spherical RB convection, the crossing between temperature gradient ($\boldsymbol {\nabla } T$) and gravitational acceleration ($\boldsymbol {g}$) results in the baroclinic torque, which plays an important role in the generation of vorticity by density variations. Figure 2(*a*) shows gravitational potential isolines to intersect the isotherms in an axisymmetric flow. The curl of the buoyancy term ($\boldsymbol {\nabla } \times \boldsymbol {g} T$) denotes the Boussinesq baroclinic torque under the Boussinesq approximation (Schladow, Patterson & Street Reference Schladow, Patterson and Street1989; Heifetz, Maor & Guha Reference Heifetz, Maor and Guha2020). As shown in figure 2(*b*), the baroclinic torque is strongest in the vicinity of the inner sphere where the angle between $\boldsymbol {\nabla } T$ and $\boldsymbol {g}$ is the largest. The baroclinic torque drives a northward-directed flow along the inner sphere. This results in a large-scale circulation pattern as shown in figure 2(*c*). In steady flows, large-scale circulations induced by baroclinic torque are also reported for side-heated cavity (Jiang, Sun & Calzavarini Reference Jiang, Sun and Calzavarini2019), stratified flow in a semicylindrical cavity (Kumar & Pothérat Reference Kumar and Pothérat2020), and natural convection in spherical annuli (Scurtu, Futterer & Egbers Reference Scurtu, Futterer and Egbers2008).

For different $n_g$, figure 3(*a*) shows the radial ($C_{g,r}g \hat {r}$) and co-latitudinal ($C_{g,\varphi }g \hat {\varphi }$) buoyancy coefficients for $\varepsilon =0.8$ on the mid-plane of the spherical shell. It shows that $g$ increases monotonically from the north to the south pole when $n_g<0$. For $n_g>0$, the opposite trend is observed, while for $n_g=0$, the gravity is constant as function of the co-latitude location. We keep $Ra=10^3$ and vary $n_g$ from $-2$ to $1$. The baroclinic torque ($\boldsymbol {\nabla } \times \boldsymbol {g} T$) for $\varepsilon =0.8$ is shown in figure 3(*b*): $\boldsymbol {\nabla } \times \boldsymbol {g} T$ concentrates in the south polar region when $n_g=-2$, whereas the strongest $\boldsymbol {\nabla } \times \boldsymbol {g}T$ gradually moves to the equator with increasing $n_g$. The corresponding large-scale circulation is shown in figure 3(*c*); the velocity field strengthens close to the inner sphere and its highest value moves from the south to the equator, which is the same trend as for the baroclinic torque. In addition, the magnitude of the velocity field decreases with increasing $n_g$ due to the decreasing mean buoyancy power ($\overline {\langle g u_r T \rangle }$, see Gastine *et al.* Reference Gastine, Wicht and Aurnou2015) averaged over the spherical shell.

### 3.2. Generation of convective jet and large-scale circulation

Figure 4 shows the temperature field on the meridional cut in figures 4(*a*–*d*), the equatorial cut in figures 4(*e*–*h*, and the surfaces corresponding to the locations of the outer and inner thermal boundary layer heights in figures 4(*i*–*p*). Figures 4(*b*–*d*) show that for $\varepsilon \ge 0.05$, the jet ejects hot plumes from the inner sphere towards the outer sphere. When $\varepsilon$ is larger, e.g. $\varepsilon \ge 0.4$, thermal plumes preferably ascend northwards along the inner sphere, as shown in figures 4(*o*,*p*). This happens because the plumes that form along the southern part of the inner sphere move upwards, following the large-scale flow pattern induced by the baroclinic torque; see figure 2(*c*). Consequently, the plumes combine at the northern part of the inner sphere to form the northbound plume towards the outer sphere (figures 4*c*,*d*). This creates a convective jet impinging on the north pole of the outer sphere, after which thermal plumes descend southwards along the outer sphere, as shown in figures 4(*k*,*l*). As a result, a limited number of plumes are ejected radially from the hot sphere to the cold sphere in the equatorial region, as shown in figures 4(*g*,*h*). In addition, we visualize the advection of temperature fluctuation by baroclinic flow $\boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla } T'$ as shown in figures 4(*q*–*t*). We find that $\boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla } T'$ is intense where $T$ is high, which indicates that baroclinic flow works on the temperature perturbation to promote the generation of thermal plumes on the northern side of the sphere. This confirms that the plumes combine at the northern part of the inner sphere to form the northbound plume towards the outer sphere.

We compare the unsteady flow fields for $\varepsilon =0.01$ and $0.05$, and different $Ra$, in figures 5(*a*) and 5(*b*), respectively. The $Ra=10^3$ cases are shown as reference in the leftmost panel. For the very small $\varepsilon =0.01$ in figure 5(*a*), the upward convective jet is the strongest for $Ra=10^7$. However, the convective jet and meridional circulation are weakened significantly for large $Ra=10^8$. For $\varepsilon =0.05$, figure 5(*b*) shows that the convective jet is the strongest for $Ra=10^6$. Again, the convective jet weakens with increasing $Ra$. In figure 5(*c*), we compare different $\varepsilon$ for $Ra=10^8$. As expected, the convective jet is enhanced significantly with increasing $\varepsilon$. However, the shape of the large-scale circulation is almost unchanged for $\varepsilon \le 0.4$. When $\varepsilon$ increases to 0.8, the large-scale circulation becomes a crescent shaped circulation, which results in its centre moving northwards.

### 3.3. Modal analysis to identify large-scale structures

When the gravity centre is moved towards the south pole (increasing $\varepsilon$), the total turbulent kinetic energy ($\mathrm {TKE}$), which is defined as

first increases slightly ($\varepsilon \leq 0.4$) and then decreases considerably ($\varepsilon = 0.8$) (see figure 6*a*). To study the effect of the off-centre gravity on the distribution of $\mathrm {TKE}$, we perform a modal analysis of the flow. In spherical coordinates, a surface spherical harmonic function is the expansion of plane waves, where the order $m$ and degree $\ell$ represent the wavenumbers along the longitudinal and co-latitudinal directions, respectively. Here, we define the power spectrum on a spherical surface as

where $a_{\ell, m}$ and $b_{\ell, m}$ are the spherical harmonics expansion coefficients, which are detailed in Appendix C. We define the power spectrum in the longitudinal direction as

where $N_\varphi$ and $N_\theta$ represent the numbers of grid points in the co-latitudinal and longitudinal directions, respectively. Radial integration and temporal averaging are used to obtain $\overline {\langle \cdots \rangle }_r$. The large-scale circulation is represented by $(m=0,\ell >0)$ modes as demonstrated in Appendix C.

Figure 6(*a*) shows that the average $\mathrm {TKE}$ is almost unchanged for $\varepsilon \leq 0.4$, but the distribution among the different modes is changed. As we have seen in figure 5, a large-scale circulation is formed when $\varepsilon \geq 0.05$ for $Ra=10^8$. Figure 6(*b*) shows that the axisymmetric ($m=0$) mode, which indicates the large-scale mode, increases strongly with $\varepsilon$. For $\varepsilon = 0.8$, both the $m=0$ mode and the total $\mathrm {TKE}$ decrease sharply, while it is still enhanced compared with $\varepsilon =0$. Interestingly, only a small gravity centre offset ($\varepsilon \geq 0.05$) is required to form the convective jet and the corresponding meridional circulation for $Ra=10^8$. The inset of figure 6(*b*) shows the distribution of the $m=0$ mode over the different co-latitudinal structures ($\overline {\langle C_{\ell, m=0} \rangle _r}$), which shows that the large-scale circulation is dominated by the co-latitudinal modes with $\ell \leq 3$. With increasing $\varepsilon$, $\overline {\langle C_{\ell, m=0} \rangle _r}$ is enhanced for low ($\ell \leq 3$) and high ($\ell \geq 70$) wavenumbers, compared with $\varepsilon =0$.

### 3.4. Effect of large-scale flow on heat transfer distribution

As we have shown above, the gravity centre offset influences strongly the flow structures. This raises the question of how the local temperature profiles and heat fluxes are affected. Figure 7(*a*) shows the mean temperature in the meridian plane. For $\varepsilon =0$, the temperature is uniform in the co-latitude direction. However, with increasing $\varepsilon$, the effect of the north-pole convective jet, which brings hot fluid from the inner to the outer sphere, becomes more pronounced. The effect of the jet is more visible in figure 7(*b*), which shows the temperature profiles at selected angular positions in the co-latitudinal direction. Most of the temperature drop occurs in the thermal boundary layer regardless of $\varepsilon$ and $\varphi$. However, for $\varepsilon =0.4$ and $0.8$, the positive temperature gradient (${\rm d}\overline {\langle T \rangle _\theta }/{\rm d}r$), observed at $\varphi ={\rm \pi} /2$ and shown in the inset, indicates some thermal stratification in the bulk region.

Next, we examine how the local heat fluxes depend on $\varepsilon$. Figure 8 shows the Nusselt number along the inner $Nu_i(\varphi )$ and outer $Nu_o(\varphi )$ spheres. It is observed that the local heat fluxes on the outer sphere decrease from $\varphi =0$ (north pole) to $\varphi ={\rm \pi}$ (south pole), whereas the opposite trend is observed on the inner sphere. This happens because the hot convective jet impinges on the outer boundary layer, resulting in a thinner thermal boundary layer and larger heat flux than without the jet. On the inner sphere, the formation of the convective jet reduces temperature gradients, and therefore the heat transport, at the north pole of the inner sphere, as shown in figure 7(*a*). Surprisingly, the formation of the large-scale circulation and corresponding non-uniform heat flux distribution has a limited influence on the overall heat flux, which globally is not sensitive to the large-scale dynamics in the flow (Ahlers *et al.* Reference Ahlers, Grossmann and Lohse2009) (see inset in figure 8*b*), although some reduction in the heat transport is observed for $\varepsilon =0.8$ when the flow becomes stratified more thermally.

### 3.5. Effect of $Ra$ and gravity profile

So far, we have focused on the effect of the gravity centre offset on the flow structures. In this subsection, instead, we show that convective jet and large-scale flow structures are formed for a wide range of $n_g$ and $Ra$, and that the control parameters modify significantly the intensity of the large-wavelength structures. However, as we have seen above, the global heat flux depends only weakly on $\varepsilon$ in the considered parameter regime.

Figures 9(*a*–*d*) show visualizations of the convective jet for $Ra=10^7$, $\varepsilon =0.4$, and different gravity profiles. The figures show that the convective jet forms with $u_r>0$ in the north pole region when the gravity centre shifts to the south. As expected, the jet becomes more turbulent when $n_g$ is decreased, as this increases the mean buoyancy power over the spherical shell. As shown previously in figure 3(*a*), the co-latitudinal buoyancy component $C_{g,\varphi }g \hat {\varphi }$ points to the north pole in all circumstances. However, the radial buoyancy coefficient $C_{g,r}g \hat {r}$ is strongly dependent on the co-latitude. It is well accepted that the stronger radial velocity component tends to induce more intense thermal plumes. However, the convective jet is still observed in the north pole region rather than the south pole region when $n_g=-2$, as shown in figure 9(*d*). This confirms that the large-scale circulation is a result of the baroclinic imbalance due to the misalignment of the gravity isopotential and isothermal lines, which is indicated by steady flows shown in figures 3(*b*,*c*).

The corresponding $Re_{rms}(r)$, presented in figure 10(*a*), shows that the turbulent intensity is enhanced greatly with decreasing $n_g$. This reveals that the flow becomes more turbulent when $n_g$ decreases, due to the increased buoyancy power over the volume. Correspondingly, figure 10(*b*) shows that the $\mathrm {TKE}$ spectrum of $\overline {\langle C_{m} \rangle _r}$ increases for all wavenumbers. We note that $\overline {\langle C_{m=0} \rangle _r}$, which dominates $\mathrm {TKE}$ and indicates the strength of the large-scale circulation, increases with decreasing $n_g$. With increasing $Ra$ from $10^7$ to $10^8$, the normalized $Re_{rms}(r)$ is enhanced as shown in figure 10(*c*). As shown in figure 10(*d*), the enhancement happens for all wavenumbers, especially for the high wavenumber structures. Thus for a constant $\varepsilon =0.4$, the large-scale circulation is more energetic for higher $Ra$ and smaller $n_g$, due to the increased mean buoyancy power over the volume and the relatively stronger baroclinic torque as indicated in figures 3(*b*,*c*).

The generation of large-scale structures can induce highly non-uniform co-latitudinal heat flux distribution. We showed previously that the global heat transfer is relatively insensitive to $\varepsilon$ for $n_g=0$ and $Ra=10^8$. Figure 11 shows that a similar trend is observed for different $n_g$ and $Ra$. In particular, for $\varepsilon \leq 0.8$, the difference between $Nu(\varepsilon >0)$ and $Nu(\varepsilon =0)$ is less than $8\,\%$, which indicates that the effect of the large-scale flow structures on the overall heat transport is relatively limited.

## 4. Summary and conclusions

We demonstrated using direct numerical simulations that even a small shift of the gravity centre location can introduce pronounced changes in the flow organization and local heat fluxes in spherical RB convection. In the Rayleigh number ($Ra$) space explored in the present study, we observe the presence of the baroclinic flow, which induces a large-scale circulation and jet for $\varepsilon \ge 0.05$. However, for very small $\varepsilon =0.01$ and large $Ra=10^8$, the baroclinic flow is weakened significantly due to the dominance of the strong convective flow. In turbulent flows with $\varepsilon \geq 0.05$, a strong convective jet on the northern side of the inner sphere is formed, which leads to hot/cold fluid going towards the north/south along the inner/outer sphere. The formation of the large-scale circulations has been confirmed by modal analysis. Due to the large-scale flow organization, the heat transfer at the inner and outer sphere is highly non-uniform. In particular, around the north pole, the heat flux is enhanced on the outer sphere due to the impingement of the convective jet, while the formation of the convective jet at the north pole of the inner sphere reduces temperature gradients, and therefore the local heat transport. However, surprisingly, the global heat flux is relatively insensitive to the gravity centre shift. These findings do not seem to depend much on $Ra$ (over our explored range) or the particular gravity profile.

## Supplementary movies

Supplementary movies are available at https://doi.org/10.1017/jfm.2022.360.

## Acknowledgements

The authors thank the three anonymous referees for constructive comments that improved the paper. We also acknowledge the Irene at Très Grand Centre de Calcul du CEA (TGCC) under PRACE project 2019215098, and the national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organization for Dutch education and research. We acknowledge PRACE for awarding us access to MareNostrum at Barcelona Supercomputing Center (BSC), Spain (projects 2020225335 and 2020235589).

## Funding

G.W. and R.J.A.M.S. acknowledge financial support from ERC (the European Research Council) Starting grant no. 804283 UltimateRB. This work was sponsored by NWO Science for the use of supercomputer facilities.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Buoyancy components of off-centre gravity

The two-step transformation of a vector $\boldsymbol {GP}$ in Cartesian coordinates $xyz$ to spherical coordinates $\varphi \theta r$ ($x''y''z''$) is shown in figure 12. A random fixed point $G(x_G,y_G,z_G)$ is representing the gravity centre location. The other random point, $P(r\sin {\varphi } \cos {\theta },r\sin {\varphi } \sin {\theta }, r\cos {\varphi })$, is a moving point representing a fluid parcel. Initially, both $P$ and $G$ are in Cartesian coordinates $xyz$. Using two steps, we will transfer $\boldsymbol {GP}$ to spherical coordinates $\varphi \theta r$ ($x''y''z''$).

(i) Rotating $xyz$ to $x'y'z'$ about $z$ by an angle $\theta$ using the right-hand side rule. The corresponding rotation matrix is

(A 1) \begin{equation} R_{z}(\theta)=\left[\begin{array}{@{}ccc@{}} {\cos \theta} & {-\sin \theta} & {0}\\ {\sin \theta} & {\cos \theta} & {0} \\ {0} & {0} & {1} \end{array}\right]. \end{equation}(ii) Rotating $x'y'z'$ to $x''y''z''$ about $y'$ by an angle $\varphi$ using the right-hand side rule. The corresponding rotation matrix is

(A 2) \begin{equation} R_{y}(\varphi)=\left[\begin{array}{@{}ccc@{}} {\cos \varphi} & {0} & {\sin \varphi} \\ {0} & {1} & {0} \\ {-\sin \varphi} & {0} & {\cos \varphi} \end{array}\right]. \end{equation}

After the above two steps we obtain that $x''y''z''$ overlaps with $\varphi \theta r$. The transformation matrix $\mathcal {M}_{\theta \varphi }$ from $\boldsymbol {GP}_{xyz}$ to $\boldsymbol {GP}_{x''y''z''}$ is defined by

In the Cartesian coordinate frame, $\boldsymbol {GP}_{xyz}=(r\sin {\varphi } \cos {\theta }-x_G,\, r\sin {\varphi } \sin {\theta }-y_G,\, r\cos {\varphi }-z_G)$. Using the above transformations, we find that $\boldsymbol {GP}$ in spherical coordinates $x''y''z''$ is

In the dimensionless equations for the velocity $\boldsymbol {u}$, the buoyancy term is given by $gT\,\boldsymbol {GP}/L_{GP}$, where

The three components in the $\hat {\varphi }$, $\hat {\theta }$ and $\hat {r}$ directions can be decomposed as

## Appendix B. Code validations

We validated extensively our code against the results presented by Gastine *et al.* (Reference Gastine, Wicht and Aurnou2015), who performed a systematic parameter study for spherical RB convection with $Pr=1$ covering radius ratios $0.2 \le \eta \le 0.95$ for different gravity profiles in table 2. The table shows that $Nu$ and the volume-averaged Reynolds number $Re_{rms}$ agree within $2\,\%$ with the results from Gastine *et al.* (Reference Gastine, Wicht and Aurnou2015) for all cases. To ensure that the off-centre gravity is implemented correctly, we verified that the flow dynamics for the same $L_{\boldsymbol {OG}}$ (see figure 1) are identical for five different gravity centre locations. The last five rows of table 2 confirm that both $Nu$ and $Re_{rms}$ obtained from these simulations are consistent.

## Appendix C. Spectral analysis in spherical coordinates

This appendix describes the spectral analysis procedure used in this study. For further reading on this topic, we refer to pp. 122–145 of MacRobert (Reference MacRobert1947). In planar configurations, Fourier transformations can be applied when performing spectral analysis. However, for a spherical configuration, the basis wave functions are found by solving the Laplace equation. Analogous to the complex exponential in planar configuration, the spherical harmonics ${Y}_{\ell }^{m}({\theta }, {\varphi })$ read

where $\ell$ and $m$ ($-\ell \le m \le \ell$) represent the wavenumbers along a meridian and the equatorial plane, respectively. The polar angle $\varphi$ ranges from $0$ to ${\rm \pi}$, and $\theta$ is the azimuthal angle, in the range $0 \leq \theta \leq 2{\rm \pi}$. Also, $P_{\ell }^{m} (\varphi )$ are the associated Legendre functions. For computational purposes, the normalized associated Legendre functions in (C2) are more attractive since $P_{\ell }^{m}(\varphi )$ can overflow the computer (Swarztrauber Reference Swarztrauber1979):

The use of the spherical harmonics for approximating functions on a sphere is motivated by the fact that they form a complete system of orthogonal functions. Therefore, any function $f(\theta, \varphi )$ that is continuous and has continuous derivatives up to second-order may be expanded in an absolutely and uniformly convergent series

where the expansion coefficients are given by

The spherical harmonics ${Y}_{\ell }^{m}({\theta }, {\varphi })$ are orthonormal, so

where $\delta _{i,j}$ is the Kronecker delta, and $\textrm {d} \varOmega =\sin \varphi \,\textrm {d} \theta \,\textrm {d} \varphi$ is the differential solid angle in spherical coordinates. Coefficients (C4) and (C5) are obtained by solving (C3) on a uniform spherical grid by using the SPHEREPACK library (Adams & Swarztrauber Reference Adams and Swarztrauber1999).

The power spectrum on a spherical surface $C_{\ell, m}$ and in the longitudinal direction $C_{m}$ is defined in (3.2) and (3.3), respectively. To demonstrate the formation of $(m=0, \ell >0)$ modes corresponding to the meridional large-scale circulation, we perform a simulation for $n_g=-2$ and $\varepsilon =0.4$ with $Ra=10^3$ (much lower than the critical Rayleigh number when $\varepsilon =0$). The advantage of using this steady flow is a regular large-scale circulation generated without longitudinal dependent structures ($m \neq 0$), which can be seen from figures 13(*a*,*b*). The spectrum of the kinetic energy, i.e. $\overline {\langle C_{\ell,m=0} \rangle }_r$, is shown in figure 13(*c*). Except for the mode of $(m=0,\ell =0)$ representing the mean flow, apparently the large-scale circulation is represented by $(m=0,\ell >0)$ modes.