Core flows and heat transfer induced by inhomogeneous cooling with sub- and supercritical convection

The amount and spatial pattern of heat extracted from cores of terrestrial planets is ultimately controlled by the thermal structure of the lower rocky mantle. Using the most common model to tackle this problem, a rapidly rotating and differentially cooled spherical shell containing an incompressible and viscous liquid is numerically investigated. To gain the physical basics, we consider a simple, equatorial symmetric perturbation of the CMB heat flux shaped as a spherical harmonic $Y_{11}$. The thermodynamic properties of the induced flows mainly depend on the degree of nonlinearity parametrised by a horizontal Rayleigh number $Ra_h=q^\ast Ra$, where $q^\ast$ is the relative CMB heat flux anomaly amplitude and $Ra$ is the Rayleigh number which controls radial buoyancy-driven convection. Depending on $Ra_h$ we characterise three flow regimes through their spatial patterns, heat transport and flow speed scalings: in the linear conductive regime the radial inward flow is found to be phase shifted $90^\circ$ eastwards from the maximal heat flux as predicted by a linear quasi-geostrophic model for rapidly rotating spherical systems. The advective regime is characterised by an increased $Ra_h$ where nonlinearities become significant, but is still subcritical to radial convection. There the upwelling is dispersed and the downwelling is compressed by the thermal advection into a spiralling jet-like structure. As $Ra_h$ becomes large enough for the radial convection to set in, the jet remains identifiable on time-average and significantly alters the global heat budget in the convective regime. Our results suggest, that the boundary forcing not only introduces a net horizontal heat transport but also suppresses the convection locally to such an extent, that the net Nusselt number is reduced by up to 50\%, even though the mean CMB heat flux is conserved.


Introduction
The cooling of the liquid iron cores of terrestrial planets is due to radial heat transport towards the core mantle boundary (CMB) via heat conduction and in case the entropy gradient is sufficiently negative supported by buoyancy driven convection. The lateral variation of heat conducted out of the core and hence through the CMB q cmb is mainly controlled by the lower mantle temperature pattern T lm ðh; /Þ, such that q cmb ðh; /Þ ¼ k T lm ðh; /Þ À T core where k is the thermal conductivity and d cmb the thickness of the thermal boundary layer at the bottom of the mantle, respectively. To first order, the core temperature is rather uniform due to a much more efficient conductive and convective heat transport therein. If the heat transport is only via conduction, thermal inhomogeneities at the CMB are thought to drive baroclinic flows (Zhang and Gubbins, 1992), whereas in a convecting core lateral variations of convective vigour, the dynamo process and the stimulation of mean horizontal flows are expected. The importance of thermal coupling between Earth's mantle and core due to inhomogeneities of the lower mantle temperature was first suggested by Bloxham and Gubbins (1987). Seismological evidence for the non-homogeneous CMB heat flux came from the mantle tomography (e.g. Masters et al., 2000) and the detection of LLSVPs (large low shear velocity provinces) (Yuen et al., 1993) which are associated with lower local mantle temperatures, but it was also reported that only for an isochemical mantle the variation of shear wave velocity and the temperature is linear (Nakagawa and Tackley, 2008). Such a mantle control has been suggested to influence the geomagnetic secular variation (Bloxham, 2000;Davies et al., 2008;Olson et al., 2010;Aubert et al., 2013) and field strength (Takahashi et al., 2008), concentrate magnetic flux patches (Olson and Christensen, 2002;Amit et al., 2010), lock the (usually drifting) core convection (Davies et al., 2009) and dynamo action (Willis et al., 2007) but also introduce mean large scale flows (Gibbons et al., 2007). Furthermore the CMB heat flux anomalies are also reported to affect the buoyancy flux from the inner core (Aubert et al., 2008;Amit and Choblet, 2009;Gubbins et al., 2011). Note, in contrast to most of these dynamo studies, we focus on the hydrodynamic aspects and hence exclude magnetic fields.
Mantle induced CMB heat flux variations are also expected to influence core flows and the magnetic field generation process of other terrestrial planets. For example the rather localised present-day magnetisation on Mars was targeted to be explained by an equatorial antisymmetric mantle-induced CMB heat flux anomaly of variable strength (Stanley et al., 2008;Amit et al., 2011;Dietrich and Wicht, 2013;Monteux et al., 2015). These studies report the induction of a global magnetic field with strong equatorially asymmetric intensity reminiscent of the recently measured distribution of magnetised crust on the surface of Mars (Acuña et al., 1999). Also the equatorial asymmetry of Mercury's magnetic field (Cao et al., 2014;Wicht and Heyner, 2014) and the axisymmetry of Saturn's magnetic field (Stanley, 2010) were investigated in the framework of a non-homogeneous CMB heat flux.
In addition, terrestrial exoplanets orbiting their host star in close proximity and typically in a synchronous orbit will receive strong stellar irradiation at the near side with a latitudinal maximum at the equator. Most likely, the heating-cooling dichotomy at the surface will result in a smooth azimuthal varying, equatorial symmetric thermal forcing pattern at the CMB reminiscent of which we use here (Y 11 ). There have been attempts to model mantle convection under such a specific heating pattern (Gelman et al., 2011) suggesting the development of a single-plume mantle convection mode. Hence the core flows and the probable induction of a core dynamo in a tidally locked terrestrial exoplanet might be significantly influenced by the enormous difference of stellar irradiation between the near and far side.
The problem of a differentially cooled and rapidly rotating fluid shell has received much attention during the last decades. Zhang and Gubbins (1992) first investigated the pure effect of laterally varying temperature at the outer boundary in the absence of radial convection. The numerical results revealed that the radial inflows are not found where the local outer boundary temperature is lowest hence cooling most efficient, but they are phase shifted by a quarter of the azimuthal wavelength of the thermal anomaly to where the azimuthal gradient of the temperature is maximal. This results from the vorticity balance between Coriolis and buoyancy terms frequently referred to as a thermal wind balance, but more specifically is a Sverdrup balance in the geophysical fluid dynamics (Pedlosky, 1979, see also below). This implies, that in a rotation dominated system thermal anomalies induce vorticity rather than flows directly. This holds as long as the flows are assumed to be rotation dominated, inviscid and nonlinearities due to temperature advection or inertia are negligible (Gibbons et al., 2007) and is hence not found in models with small rotation rates or infinite Prandtl number (Sun et al., 1994;Zhang and Gubbins, 1996;Davies et al., 2009). A mathematically more straight-forward analysis of the linear quasi-geostrophic model (Busse annulus) by Yoshida and Hamano (1993) including the effect of the magnetic fields, confirmed the azimuthal phase shift when there is no magnetic field altering the leading order force balance.
For a physically more realistic model of the Earth's core, experiments by Sumita and Olson (1999) perturbed the outer boundary of a vigorously convecting and rapidly rotating spherical shell with a local anomalous heat flux. There was also a phase shift between local minimum temperature and position of the inward flow reported, however due to the strongly nonlinear driving it was deformed in a jet-like structure spiralling inwards. As the characteristic hydrodynamic numbers (e.g. Ekman or Reynolds number) in experiments are closer to real planets than those accessible by numerical models, this might reflect better the situation in a realistic planetary core. In a follow-up study (Sumita and Olson, 2002) a broad parameter survey was reported, featuring a detailed description of the flow and scaling relations of how the heat flux anomaly, pattern and strength affects the induced flow amplitudes in azimuthal and radial direction. Typically these experiments are set in a strongly convective regime, where a rather localised and very strong heat flux anomaly induces a sharp front separating the cold east from the hot west and strong azimuthal flows connecting them.
As some of the reported effects are due to complex interactions between core convection, boundary forcing and magnetic fields, a clear physical interpretation might not be always possible. Thus to gather clearer insights we limit this study to the hydrodynamic aspects for the sake of a more coherent physical description. Starting from an analytical formulation of the linear theory, we numerically model core flows induced by thermal CMB inhomogeneities for cores subcritical to buoyancy driven core convection (baroclinic flows) and compare them to models featuring radial convection. Note, also the baroclinic flows obey strong nonlinear features, such as bending or compression of the emerging inward jet. As the various linear and nonlinear flow regimes have been studied only individually, our main focus is to distinguish them and discuss their properties.
More precisely we modify the outer boundary heat flux of rapidly rotating spherical shell with a smoothly varying pattern along azimuthal and latitudinal direction. For simplicity the heating of the shell is exclusively supplied by a constant internal heat source modelling secular cooling of the shell. It was shown that internal heated systems are most sensitive to inhomogeneities at the outer boundary (Hori et al., 2014), as the strongest temperature gradients are typically found close the outer boundary. The specific shape is an anomaly of the outer boundary heat flux of spherical harmonic degree and order one (Y 11 ). This purely equatorially symmetric and non-axisymmetric pattern was taken mainly for application to terrestrial exoplanets orbiting their host star in synchronous rotation. Note, today's Earth CMB heat flux variation is dominated by a sectorial spectral mode Y 22 , but the core convection is mainly driven from compositional buoyancy release at the inner core boundary. Hence if applied to real planetary systems, our models best describe cores of terrestrial planets before inner core nucleation. Apart from the different azimuthal length scales, our linear results might be applicable to a general sectorial anomaly pattern as shown below on Section 3.1. However, for the nonlinear results the interaction between several sectorial heat flux anomalies are beyond this study. Emphasis is put on a clear physical description of the induced flow structures and the radial (and horizontal) heat transport.
We also vary the amplitude of the CMB heat flux variation q Ã relative to the mean value as it controls the strength of the boundary forcing to first order and is not even well known for the Earth's core. We are interested in how sensitive convection is to a thermal boundary anomaly of variable strength. Studying q Ã might provide a smooth transition between models with homogeneous and heterogeneous CMB heat flux. Note, we limit q Ã to unity to avoid any heat flux into the core which might not be consistent with our model assumptions. A fundamental property of the here applied Boussinesq-approximation is that all temperatures and heat fluxes are understood as superadiabatic fluctuations on top of the (excluded) adiabatic background state of heat conduction. Hence a (q Ã > 1)-heat anomaly describes locally a superadiabatic heat flux into the core in a Boussinesq-model, whereas the physical meaning is a subadiabatic heat flux out of the core.
In a stably stratified, rapidly rotating fluid shell motions are introduced due to baroclinic torques emerging from the misalignment of surfaces of constant pressure and constant density. As we consider a fluid in Boussinesq approximation, the equation of state simplifies to a linear relation between density and temperature due to thermal expansion. Hence the leading order vorticity balance is given by the curl of the Navier-Stokes equation

À2X
@ṽ @z where X is the rotation vector,ṽ the (dimensional) flow velocity, q the reference density, a the thermal expansivity, T the temperature andg the linear increasing gravity. Physically this implies, that any horizontal temperature variation introduces a finite flow (baroclinicity). However, if the temperature varies only along radius, pressure and density are aligned and the flow comes to a rest (hydrostatic equilibrium). The modelling strategy of this paper attempts to cover all possible linear and nonlinear flow states. We start off with deriving the analytical description like Yoshida and Hamano (1993), but for spherical shells rather than a cylindrical annulus, adding the ageostrophic part of the flow and without magnetic fields. These results are verified by numerical simulations set in the almost linear (conductive) regime (Section 3.1). Then we increase the driving to reach a nonlinear regime, that is yet not supercritical to radial convection and discuss the properties of the emerging inward jet (Section 3.2). To understand how radial convection and the inward jet interact, we finally analyse the convective regime where radial convection is strongly participating. Furthermore it is an open question, how and to which extent boundary forcing alters the global heat transfer budget. Therefore an investigation of the threedimensional redistribution of heat and the mean heat transport properties such as the Nusselt number is described in Section 4. We also parametrise the influence of the vigour of convection and the various relative forcing amplitudes.

Methods and model
The outer core spherical shell with outer and inner radius, r o and r i , is assumed to be filled with an incompressible and viscous liquid. Within the Boussinesq-approximation the evolution of the fluid flow is governed by the dimensionless Navier-Stokes equation: whereũ is the incompressible velocity field (r Áũ ¼ 0), P the nonhydrostatic pressure,ê z the direction of the rotation axis and T the super-adiabatic temperature fluctuation on top of a hydrostatic equilibrium state. We use the shell thickness D ¼ r o À r i as length The evolution of the thermal energy is affected by thermal diffusion and advection along with the flow, such that @T @t þũ ÁrT ¼ 1 Pr where is a uniform heat source density.
We assume non-penetrative and no-slip mechanical boundary conditions because the mantle and inner core are solid. When core convection is driven by secular cooling only, we set the heat flux at the outer boundary to unity (q o ¼ 1) and apply zero flux on the inner core boundary (q i ¼ 0). From the modelling perspective, the system is then powered exclusively by the internal heat source and a simple balance between and the mean CMB heat flux q o is given by where b ¼ r i =r o is the radius ratio of the spherical shell. Such an internal heating setup models a shell undergoing secular cooling, but neglects compositional and thermal buoyancy arising from inner core solidification. The thermal effect of inhomogeneous cooling is modelled with a spherical harmonic of degree and order unity (Y 11 ) with variable amplitude q Ã ¼ ½0 . . . 1:0, such that The boundary anomaly proportional to Y 11 provides a smaller than average heat flux in the sectors I and IV (see Fig. 1). Actually, for q Ã ¼ 1 it is exactly zero at / ¼ 0. In the sectors II and III the heat flux is enhanced, with a maximum of 2q o at / ¼ p. In between the spherical harmonic description guarantees a smooth variation (see also Fig. 1b). The emerging horizontal difference in convective intensity is measured in terms of a horizontal Rayleigh number Ra h (Willis et al., 2007;Monteux et al., 2015), such that As we use the Boussinesq-approximation, the horizontal Rayleigh number shall not exceed the (radial) reference Rayleigh number of the homogeneous model. This is assured when keeping q Ã 6 1. The combination of radial and horizontal temperature gradients will lead to a complex superposition of various flows. Fig. 1 illustrates the equatorial plane of the three-dimensional model and clarifies the nomenclature used throughout the paper. X gives the sense of rotation, the azimuthal angle / is counted in prograde direction starting from the right hand side. Hence the upper half (0 6 / 6 p) is termed eastern, the lower western hemisphere. For convenience the plane is further subdivided into the four azimuthal sectors (I to IV). The heat flux anomaly is positioned such that, the hottest mantle feature is on the right at / ¼ 0, where then also the smallest amount of heat is extracted from the core as indicated by the small red arrow in the figure. Thus the largest CMB heat flux is at the opposite point (at / ¼ p). We expect that in the absence of radial convection, the boundary anomaly should drive a mean flow consisting of two cells, a (eastern) cold cyclone and a western hot anticyclone. The cells converge into an upwelling (downwelling) where the CMB heat flux is small (large). However, the dominant Coriolis force will turn the entire flow system. Indeed, as we shall show later, for models featuring radial convection the downwelling is typically shifted by an angle 0 6 U 6 90 somewhere in sector III. For clarity, we measure the phase shift U in degrees and generic azimuthal angle / in radians. Earlier weakly nonlinear numerical (Zhang and Gubbins, 1992) and linear analytical studies (Yoshida and Hamano, 1993) suggested a phase shift of U ¼ 90 between the heat flux maximum and the radial inflow, whereas for a non-rotating model U ¼ 0 would be expected.
We numerically solve the governing equations (Eqs. (3) and (4)) with the pseudo-spectral MagIC3 code (Wicht, 2002). The radial and horizontal resolution is mainly 49 and 288, respectively, but is increased up to 61 and 320 for the computational most demanding cases with high Rayleigh numbers. The radial resolution is not equidistant featuring higher grid point density close to the boundaries. E.g., the thickness of the Ekman layer in the vicinity of the boundaries is given by d=D % E 1=2 ¼ 0:01 and is resolved with several radial grid points. In addition the spectral decay in the flow spectrum from the largest to the smallest length scales was checked to consistently span three orders of magnitude. In total 80 numerical models with different setups are solved, where each individual case is time integrated from an initial condition taken from a convecting or randomised solution until a steady or statistically stationary state is reached and time-averaged over a viscous diffusion time. A steady state is characterised by an exactly vanishing time derivative, whereas a statistically stationary state is reached when the kinetic energy fluctuates around a constant mean value. For investigating qualitatively the various linear and nonlinear flow regimes, 28 different models covering 8 orders of magnitude in Rayleigh number are calculated, where all other parameters are kept fixed (q Ã ¼ 1: not stated otherwise. For reference we ran 26 equivalent homogeneous models (q Ã ¼ 0:0). For the combination of Ekman number, Prandtl number, aspect ratio and the internal heating mode with flux boundaries, an equivalent homogeneously cooled model would require a Rayleigh number of at least Ra c ¼ 1:6 Á 10 6 for the onset of convection. Note, this value was determined numerically by monitoring the growth or decay of random perturbations from the basic state. For analysing the transition from homogeneous to inhomogeneous boundary heat flux 8 different amplitudes (q Ã ) are tested as well. To get an idea of the Ekman number dependence there are another 16 models with E ¼ 5 Á 10 À5 .

Results
In contrast to the flows introduced by the boundary anomaly, buoyancy driven convective instabilities are opposed by the stabilizing pressure gradient and require a driving stronger than a critical value (Ra > Ra c ). According to the vorticity equation (Eq. (2)) there is always a finite flow for inhomogeneous outer boundary heat flux independent of the presence of radial convection. The horizontal Rayleigh number Ra h serves here as the main parameter to find different regimes. Verifying our analytical results of the linear, spherical problem we start the numerical survey at Ra h ¼ 10, which is strongly subcritical (Ra h =Ra c ¼ 6:25 Á 10 À6 ). As long as the resulting flow amplitudes in terms of the Reynolds or Péclet number, Re or Pe, are less than unity the nonlinear temperature advection and inertia terms are negligible, hence this is the linear (conductive) regime. For the nondimensional viscous time scale chosen here, Re is equivalent to the rms flow speed and Pe ¼ RePr. Note, as the Prandtl number is fixed throughout the paper to Pr ¼ 1:0, Reynolds and Péclet number give identical measures of the flow amplitudes. Increasing Ra h will increase Re as well, and the thermal advection is found to strongly modify the flow. Hence the nonlinear (advective) regime is defined by nonlinear flows, which are still subcritical to convection. We set the regime boundary to the nonlinear convective regime, at the critical Rayleigh number Ra c for the homogeneous system (Ra h =Ra c ¼ 1). Hence a suite of numerical simulations starting with Ra h ¼ 10 and reaching up to Ra h ¼ 6:4 Á 10 8 (Ra h =Ra c ¼ 400) is expected to cover all possible flow states and transitions. As the Ekman number is fixed at E ¼ 10 À4 , a further enhancement of Ra or Ra h might lead to an inertia-dominated regime, which is unrealistic for planetary applications. Zhang and Gubbins (1992) predicted for the rapidly rotating, non-magnetic and inviscid limit (E ! 0), that radial flows are determined by the azimuthal temperature gradients rather than the outer boundary radial heat flux. Hence up-and downwellings are not expected to be found where the temperature is largest and smallest, but where the azimuthal temperature gradient is maximised. A mathematical more straight-forward analytical study by Yoshida and Hamano (1993) including magnetic fields confirmed this result within the linear theory of the Busse annulus (Busse, 1970). To understand the physical origin of these flows, we follow the approach of Yoshida and Hamano (1993) but generalise it to spherical shells, add the ageostrophic flows while neglecting the effect of the magnetic field. Numerical simulations with very small Ra h are used to test the analytical results.

Linear conductive regime
For rapidly rotating systems the thermodynamic balance between Coriolis and buoyancy force can be found by taking the curl of the Navier-Stokes equation (Eq. (3)). For the inertia-less, inviscid and steady limit, the leading order vorticity balance (Eq. (2)) is then given in dimensionless form by We separate the flow into its geostrophic (ũ g ) and ageostrophic part (ũ a ). As geostrophic flows are by definition invariant along the axis of rotation (ê z ), it is convenient to separateũ g andũ a using cylindrical coordinates (s; z; /): u ¼ũ g ðs; /Þ þũ a ðs; /; zÞ: The geostrophic part is represented by the z-averaged flow (ũ g ¼ hũi z ), where the averaging operator is defined by for a generic function f and H ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi is the s-dependent height of a spherical shell measured from the equator outside the tangent cylinder (TC). The geostrophic flow can be found by investigating the z-averaged z-component of Eq. (8): The left hand side of the equation can be directly integrated. Incorporating the sloping boundaries conditions (e.g. Busse, 1970;Gillet and Jones, 2006;Hori et al., 2015) gives For the right hand side, the balance (Eq. (11)) has in general a latitudinal and an azimuthal component, where the latter is antisymmetric with respect to the equator for any equatorial symmetric boundary anomaly and averages out (h @T @h i z =0). Hence z-averaging the RHS of Eq. (11) yields 1 2H Finally evaluating the equation in the equatorial plane, where s ¼ r, u g s ¼ u g r ,ê h Áê z ¼ À1 gives the solution for the radial flow outside TC u g s;r ¼ À From this equation it can be seen that the radial geostrophic flow is amplified where the azimuthal gradients of the temperature are large, i.e. a phase shift between the two. Note, this relation for the baroclinic driving of geostrophic radial flows holds for any azimuthal variation of the temperature. This is equivalent to the Sverdrup relation, which is commonly known in oceanography (Pedlosky, 1979). There usually the thermal anomaly is replaced by the curl of the wind stress exerted onto the free surface of an ocean. The geostrophic azimuthal flow (u g / ) will connect the up-and downwellings and is found from the z-averaged incompressibility condition in cylindrical coordinates. The solenoidal condition for the geostrophic flow simplifies to: Hence the azimuthal, geostrophic flow u g / is not phase shifted relative to, but obeys a more complex radial structure than the temperature field. As can be seen, azimuthal flows converge (diverge) where radial flows are increasing (decreasing) along radius. We therefore expect a prograde (retrograde) azimuthal flow in sectors I and III (sectors I and IV).
We discuss the temperature and the ageostrophic part of the flow as well. As the heat flux anomaly is given in spherical harmonics, switching back to spherical coordinates (r; h; /) is physically and mathematically more suitable, however on the expense of a mild inconsistency with the cylindrical treatment of the geostrophic flow. If temperature advection is negligible, Eq. (4) simplifies to a diffusion equation subject to the outer boundary heat flux variation, which shall be shaped like a general sectorial spherical harmonic Y mm ¼ sinðm/Þ P m m ðcos hÞ with amplitude unity: As there is no temperature advection considered, the boundary anomaly will diffuse as a temperature anomaly through the entire shell, hence where RðrÞ is a undetermined function describing the radial dependence only. The azimuthal and latitudinal gradients can be given as: As a consequence only non-axisymmetric (m ! 1) disturbances of the CMB heat flux drive radial geostrophic flows, where then the azimuthal displacement between the temperature and radial flow is a quarter of the azimuthal variational wave length. Further we note, the z-average of Eq. (21) disappears always for sectorial patterns as the integrand (cot h P m m ðcos hÞ) is antisymmetric with respect to the equator.
By definition geostrophic flows are invariant along z, hence the ageostrophic flow is found in the spherical components of the Eq.
Hence the ageostrophic part of the flow is found by the temperature gradient along colatitude for the axial variation of the azimuthal flow and along azimuthal angle for the latitudinal flow variation. The azimuthal component (Eq. (24)) is more commonly known as the thermal wind balance in thick spherical shells (Jones, 2007). The ageostrophic flows are entirely determined by Eqs. (22)-(24) hence act only along azimuth / and latitude h, but are non-radial. As the heat flux anomaly gets weaker towards the poles, the latitudinal temperature gradient is either negative where the equator is hotter than the pole (/ ¼ 0) or positive where the equator is colder than the pole (at / ¼ p). A schematic flow is shown in Fig. 2, panel b. That implies (Eq. (24)) a prograde (retrograde) azimuthal flow at the equator in sectors II and III (IV and I), hence the far side (near side) of the planet. The ageostrophic azimuthal flow u a / is thus consistent with its geostrophic counterpart u g / and retrograde (prograde) at the front (far) side of the shell in the equatorial plane.
Furthermore, as the temperature decreases (increases) in the eastern (western) hemisphere, the latitudinal flows are expected to converge (diverge) at the equator. Note both flow parts, geostrophic u g r ðr; /Þ; u g / ðr; /Þ and ageostrophic u a r ðh; /Þ; u a / ðh; /Þ, provide a unique two dimensional closed circulation seeking to equilibrate the heat flux anomaly. Fig. 2 schematically depicts the expected flow structures and compares to the simulations. Subpanel (a) shows that the numerical results are consistent with the theoretical treatment for the geostrophic flow, where (b) shows the temperature distribution and the horizontal flow components at a spherical surface for the ageostrophic flow. The reversed u / at higher latitudes is the ageostrophic equivalent of the geostrophic return flows closer to the inner core. Whereas in the geostrophic flow an azimuthal temperature gradient drives a z-independent radial flow, for the ageostrophic flow this results in an equatorially antisymmetric convergence or divergence of latitudinal flow. Hence equatorward flows are expected in the eastern hemisphere, whereas polewards flows are driven in the western hemisphere. Note, that the ageostrophic flow does not require any radial flows in order to equilibrate the heating anomaly, however the geostrophic part cannot be neglected as it exceeds the ageostrophic part in amplitude and will naturally emerge in any rapidly rotating, boundary forced system. In addition from the Eqs. (14) and (23) it can be seen, that the (purely geostrophic) radial flow (u r ¼ u g r ) and the (purely ageostrophic) latitudinal flow gradient (@u h =@z ¼ @u a h =@z) are both proportional to @T=@/ and hence are governed by the azimuthal temperature gradients rather than the temperature itself. For the azimuthal flow (u / ), the geostrophic part found from the Eq. (15) and the ageostrophic flow gradients defined by Eq. (24) are proportional to the temperature rather than the azimuthal gradient.
In addition we numerically isolate the geostrophic (z-independent) part of the flow resulting from the numerical models by transforming from spherical to cylindrical coordinate representation and measuring the rms flow amplitude by Re where the index i stands for either radial (r), latitudinal (h) or azimuthal flow (/). Here we use the time-averaged flow assuming only the boundary induced flow structures are time-persistent. The geostrophy as the ratio Re g i =Re a i is calculated from the simulations and listed in the Table 1 for several Ra h . For the latitudinal direction, indeed the geostrophic contribution is much smaller than the ageostrophic as expected from the analytical results. Interestingly this is less clear for the (as expected purely geostrophic) radial flow, where the ageostrophic part is probably larger than expected due to boundary or viscous effects. As evident from Table 1 the mean azimuthal Reynolds number of the geostrophic flow is seven times larger than the ageostrophic.  Table 1 Geostrophy of the radial (Rer), latitudinal (Reh) and azimuthal (Re/) calculated as the ratio of the geostrophic and ageostrophic Reynolds number (Re g i and Re a i , respectively) for the six study cases from Fig. 4 To visualise the ratio of geostrophic and ageostrophic flow contributions the bending of lines of constant u / is plotted. Fig. 3 depicts the azimuthal flows u / , temperature T and the left-and right-hand side of the thermal wind (Eq. (24)) in meridional cuts along / ¼ 0 and / ¼ p. The rather weak bending confirms the geostrophic flow (at least the azimuthal) exceeds the ageostrophic flow by far. However it can be seen, the weak variation of u / along z is given by the azimuthal thermal wind. Interestingly, an additional shear flow feature not related to the thermal wind appears at the TC. This was suggested to be caused by viscous effects on the geometric discontinuity across the TC creating azimuthal shear flows (Livermore and Hollerbach, 2012).

Nonlinear regimes: advective and convective
As convincing and clear the analytical results are, the simple linear treatment oversimplifies a more realistic scenario. The nonlinear advection terms in the temperature and momentum equation will significantly contribute and complicate flow and temperature once their nondimensional measure (Reynolds and Péclet number) reaches unity. We therefore numerically investigate the effect of an increasing (horizontal) Rayleigh number and characterise the different regimes.
As an overview, Fig. 4 shows the time-averaged flow and temperature distribution found in the equatorial plane for several Ra h ranging from Ra h ¼ 10 to 3:2 Á 10 8 . As extensively discussed above, for the smallest Ra h (top row) the flow structure is in agreement with numerical and analytical results (Zhang and Gubbins, 1992;Yoshida and Hamano, 1993). Note, a few plots from Fig. 2 are repeated to complete the overview. As introduced in the Section 3.1 an equatorial stream function yields a clearer characterisation. Radial and azimuthal flows can then be deduced by u r ¼ 1 r @W @/ and u / ¼ À @W @r ; ð27a; bÞ in spherical coordinates (cf. Eq. (16a,b)). The anticylconic structure further east of the jet appears negative (red), whereas the cyclonic structure is positive (blue).
To better quantify the different phase relations, azimuthal profiles of mean temperature T, azimuthal temperature gradient @T=@/ and radial flows u r are defined by averages with respect to time, radius and colatitude: where hf ðr; h; /Þi t ¼ 1=Dt R f ðh; /; tÞdt is the time average operator. Note, ju 0 r jð/Þ is derived as the time averaged intensity of nonaxisymmetric radial flow in individual snapshots, whereas u r ð/Þ is the remaining radial flow after time-averaging. Fig. 5 plots the azimuthal profiles of T; dT=d/; u r and ju 0 r j. The black curves correspond to the upper case (Ra h ¼ 10, conducting regime) from Fig. 4 and show nicely the anti-correlation between the radial flow and azimuthal temperature gradient (Eq. (8)) and the p=2 phase shift between the temperature T and u r . Note there is a small offset of ca. 15 between u r and @T=@/, that we attribute to the smallness of the flow Reynolds-number (Re ¼ Pe ¼ 6:7 Á 10 À4 ). Hence the flow is affected by viscosity and tends slightly to a viscous-buoyancy balance, where u r $ T. At smaller Ekman numbers this offset is expected to vanish as viscosity is not important in such a limit (Zhang and Gubbins, 1992).  Table 1. Parameters: E ¼ 10 À4 ; q Ã ¼ 1:0.
If the driving is increased to Ra h ¼ 3 Á 10 4 , the symmetry between the eastern and western hemisphere is broken (Fig. 4b).
We picked this case, because the mean flow Reynolds-number firstly exceeds unity here (Re ¼ Pr % 1:8), hence the effect of temperature advection becomes significant. As shown in the plot the downwelling is stronger and more concentrated than the upwelling. This is caused by the azimuthal temperature advection Àu / =s @T=@/, which is added on the right-hand side of the figure. If this term is locally negative (positive) the temperature is increased (decreased). As can be seen from its structure the nonlinear advection provides heating (cooling) in sector I and III at larger (smaller) radii, and vice versa in sector II and IV. Due to the incompressibility, the flows at the larger radii are stronger and global effect is then an increase (decrease) of the azimuthal temperature gradient in the western (eastern) hemisphere. In response, the downwelling (upwelling) in the western (eastern) hemisphere is enhanced (suppressed) and the azimuthal flows are concentrated in the west. At a slightly higher Ra h ¼ 3 Á 10 5 (Fig. 4c), the upwelling has almost ceased whereas the downwelling developed a strong and confined jet-like structure, which is inclined in prograde direction while descent.
The azimuthal profiles in Fig. 5 (red and brown) show that temperature flattens out along azimuth, but obeys a strong slope in the vicinity of the jet. This clearly indicates that the jet is separating the cold east from the hot western hemisphere. As discussed before, the local azimuthal temperature gradient is directly responsible and hence indicative for the inward traveling jet. Remarkably, the temperature T looses its phase relation with Fig. 4. Flow regimes and transitions. Each row shows (from left to right) the time average of temperature, azimuthal gradient of temperature, equatorial stream function, the nonlinear temperature advection term along azimuth in the equatorial plane and radial (upper) and azimuthal flow (lower) on a spherical surface at mid-depth. The (horizontal) Rayleigh number is increased from (a) Ra h ¼ 10, (b) Ra h ¼ 3 Á 10 4 , (c) Ra h ¼ 3 Á 10 5 , (d) Ra h ¼ 5 Á 10 6 , (e) Ra h ¼ 4 Á 10 7 to (f) Ra h ¼ 3:2 Á 10 8 . Contour maxima and minima are indicated at the bottom right in each panel. Note, Table 1 lists ratios of geostrophic and ageostrophic flow amplitudes for the depicted cases. Parameters: E ¼ 10 À4 ; q Ã ¼ 1:0. respect to the radial flows, but the strong anti-correlation between u r and @T=@/ still holds. Hence also for nonlinear regimes, where advection or convection is significant, on time average the boundary induced flow is still fully determined by the temperature structure. Note, this was also successfully tested for the experimental setup of Sumita and Olson (1999), where a model thermal anomaly was applied to find the mean flow structure analytically (Sumita and Yoshida, 2003). In this regime radial temperature (density) variations are enhanced by Ra h , but not sufficiently to reach unstable stratification. As the major nonlinearity stems from the temperature advection this regime might be termed as 'advective'.
As the jet sinks in and gains in amplitude with increasing Ra h , the double cyclone structure is deformed and the strong Coriolis force might deflect any inward radial flow into prograde azimuthal direction. Alongside with the jet compression, Fig. 4 shows that a radial dependence of the jet phase shift UðrÞ is emerging with increasing Ra h for the advective regime. The phase shift close to the outer and inner boundary, U o ¼ Uðr o Þ and U i ¼ Uðr i Þ, respectively is plotted in Fig. 6, whereas the difference j U o À U i j estimates the bending. It can be seen, from Ra h ¼ 3 Á 10 5 onwards significant jet bending occurs, as the jet anchor points move prograde in the interior and retrograde close to the outer boundary. The reason for a more retrograde position of the upper boundary anchor point seems to be related to the asymmetry between the cold cyclone and the hot anticyclone. The former features a prograde flow close to the outer boundary, whereas the latter is retrograde. Both flows advect the jet, but the anticyclone is increasingly stronger with higher Ra h and hence the jet sinks down at smaller /. On the other hand, close to the inner boundary the azimuthal flow in the dominating anticyclone is prograde yielding a westward shift of U i , which can be seen in Fig. 6, red line. Note, this trend was also observed in the experiments by Sumita and Olson (2002). The total jet bending hence increases with Ra h until the onset of convection is reached.
Radial convective instabilities driven by unstable buoyancy stratification might participate when the buoyancy scaling factor Ra h supersedes the critical value Ra c for the homogeneous model (given in Section 2). Note, it is not clear how the boundary forcing affects the onset of convection. However, we use Ra c as the regime boundary between advective and convective regime as there are several key properties coinciding. Besides being the critical Rayleigh number in the homogeneous case at Ra ¼ Ra c , we also observe steady flows for Ra h < Ra c and fluctuating flows for Ra h > Ra c which can be evident when taking snapshots, as shown below in Section 4. The time dependence might be related to both, emergence of time-dependent convective flows and a wiggling of the jet. For strongly convecting models (Ra h ) Ra c ) at any instance of time, the jet is hardly identifiable but remains persistent after time-averaging. It is thus not the same periodic fluctuations caused by periodic formation and destruction of the front observed in the similar experimental setup by Sumita and Olson (2002). π/2 π 3π/2 2π azimuth φ |u' r | tively. The jet bending in the advective regime stemmed from the dominant anticyclone yielding a deformation of both cyclonic structures. As the radial jet defines the temperature front, azimuthal bending implies strong radial gradients through the jet. It is reasonable that the radial convective motions setting at Ra h ¼ Ra c homogenise radial temperature gradients and only allow for azimuthal gradients to persist. As a consequence the jet is dominantly radial in the convective regime. However for higher Ra h ) Ra c the radial jet is azimuthally deflected by the Coriolis force yielding a weak, but consistent bending.
Throughout the convective regime with Ra h > Ra c , the jet separates the cold cyclone and the hot anti-cyclone and the strongest radial temperature gradients develop (as expected) at the prograde edge of the eastern cyclone. In opposite to the retrograde azimuthal flow in the hot anticyclonic western cell, the prograde azimuthal flow in the eastern cyclone advects hot fluid along a strong azimuthal temperature gradient and hence provides cooling. If buoyancy is allowed and a fluid parcel is sufficiently cooled it can become negatively buoyant due to unstable stratification. Hence, once cooled enough it is the jet itself, that becomes negatively buoyant (heavy) and sinks down at smaller U (further east).
Interestingly the buoyancy is not only responsible for the jet to sink in but also cools the eastern cyclone by radial convection. Even for the time-dependent convective regime ( Fig. 4e and f), amplitude and position of the jet are perfectly anti-correlated to the azimuthal temperature gradients (Fig. 5, orange and green).
For Ra h ¼ 4 Á 10 7 (Fig. 4e) the flow has reached full nonlinearity in the sense, that both radial and azimuthal temperature variations drive flows as the critical Rayleigh number for radial convection is exceeded by far (Ra h =Ra c % 25). However, the azimuthal profile of the intensity of nonaxisymmetric flows ju 0 r j (Fig. 5, bottom plot) shows the dominance of the jet and that the radial convection in the colder, eastern cyclone is much more energetic than in the hotter, western anticyclone. In addition it can be seen, that the stronger the radial convection gets, the closer it is to the dominant jet. As the jet is eastward phase-shifted, the azimuthal profile of ju 0 r j does not follow the locally prescribed boundary heat flux (as introduced by the anomaly), but does follow the temperature distribution created by the mean (boundary forced) flow. In other words, the phase-shifting effects of the Sverdrup balance due to rotation also translate into the 'regular' radial convection.

Flow speed scaling
An interesting question is whether the jet amplitudes can be related to the system parameters, such as the boundary forcing amplitude or the rotation rate, and how dominant they are in the presence of additional radial convective instabilities. For the conductive linear regime, that is when the flow is controlled by the boundary anomalies and advective terms are not important, the axial vorticity equation (Eq. (8)) should predict the amplitude of the global circulation. Assuming that the azimuthal temperature gradient is well represented by the boundary forcing Eq. (8) can be approximated by where ' is the variational length scale for flow and temperature. This leads to the flow scaling: where the temperature scale is approximated by the boundary anomaly strength (DT % q Ã ). As a measure of the kinetic energy the global flow kinetic Reynolds number of the time averaged flow is used. Making use of the time-persistence of the boundary induced flows, Re takes radial, azimuthal and longitudinal flows into account but might ignore rapid flow fluctuations such as the convective flow. Fig. 7 tests this scaling for various combinations of Ra (red), E (blue) and q Ã (green). The dot-dashed grey line indicates the linear slope and resembles the mean flow amplitude in terms of Reynolds number surprisingly well up to a certain point equivalent to Ra h ¼ 3 Á 10 4 (see also Fig. 4, 2nd row). All red cases are calculated with fixed E ¼ 10 À4 ; q Ã ¼ 1 and variable Rayleigh number ranging from 10 to 6:4 Á 10 8 . For the blue cases in the figure the Ekman number was halved, and the resulting Reynolds numbers seem to be halved as well. As the third factor of the thermal wind scaling (Eq. (33)), two test cases with halved q Ã ¼ 0:5 and doubled Ra fall back to mean trend.
It can be seen from the Fig. 7, when reaching Ra h E=2 % 1 the linear relation fails and overestimates the flow amplitude. At this point the nonlinear temperature advection term (ũ Á rT) overcomes the thermal diffusion as indicated by flow the Reynolds or Péclet number (as Pr = 1) reaching Re % 1. We set this as the boundary between linear conductive and nonlinear advective regime. A further increase in Rayleigh number leads to the previously discussed jet focusing and amplification between Ra h ¼ 3 Á 10 4 (advective regime) and the onset of radial convection at Ra c % 1:6 Á 10 6 (convective). For cases with supercritical convection, the flow amplitude seem to be much better resembled by a

Ra 1=2
h -scaling. Such an exponent is theoretically expected for a buoyancy-viscosity balance ) and well confirmed by numerical models and experiments (Christensen and Aubert, 2006;. From the axial vorticity equation (Eq. (8)) Sumita and Olson (2002) also suggested a flow scaling / Ra 1=2 h , but this was deduced for the boundary induced azimuthal flow only. We hypothesise, that our global measure is in agreement with their findings as the azimuthal flow dominates the timeaveraged kinetic energy. However for the nonlinear advective regime this simple scaling law fails to predict the correct mean flow amplitude.
The kinetic energy corresponding to the two cell flow might be spectrally strongly dominated by kinetic energy contributions obeying a spherical harmonic order of m ¼ 1. A measure of the relative strength of (m ¼ 1)-flows is given by Hori et al. (2014) as E=1e-4 q*=0.5,2Ra E=5e-5 1-Em1 As a consequence of the jet focusing, Em1 is expected to decay fast when Ra h is increased and the spectrum is not as monochromatic (with m ¼ 1), but involves higher modes. Fig. 7 shows the deviation from a monochromatic (m ¼ 1)-flow. In the plot this quantity increases with a slope of 2 for Ra h E=2 < 1 (conductive regime), undergoes a transition (advective) and saturates to a constant value corresponding to Em1 % 0:24 for Ra h E=2 > 100 in the convective regime, where radial convection is active. Note, at the first transition between linear conductive and nonlinear advective regime, Em1 has only marginally decreased, e.g. at Ra h ¼ 3 Á 10 4 we find Em1 ¼ 0:961.

Heat transport in the convective regime
This section contains a comprehensive description of how the radial convection is affected and how the temperature anomaly introduced by the boundary anomaly is extracted from the core. Radial heat transport in a convective flow is given by the radial component of the heat advection term and should usually direct radially outwards. This implies that a positive (outward) radial flow provides cooling if the temperature decreases with radius hence hot material will be advected towards the cooling surface. Note, these flows will be accompanied with the same amount of negative radial flow pumping cool material to the hotter interior (cooling as well). However, if the flow is time-dependent and fluctuating there will be also inward heat flux contributions. We therefore investigate the radial heat transport in a series of individual snapshots of a time-variable flow, collect the contributions of the radial heat transport into outward (j o ) and inward (j i ) and timeaverage their intensities. Afterwards the sum of the two might give an idea of the net (outward) radial heat transport. Locally the direction of the radial heat transport is determined by the sign of the radial temperature gradient, hence @T @r where both quantities will be present in a fluctuating flow. However, the net radial heat flux will be positive as the fluid shell is cooling. The classic Nusselt number as the nondimensional ratio between convective and conductive heat flux is related to the global average of the radial heat transport j r (e.g. Otero et al., 2002) by As we use fixed flux boundary conditions, the normalisation temperature scale DTðRaÞ is an output quantity. The temperature scale usually drops with increasing convective vigour as stirring and mixing is more efficient. In accordance to the radial heat transport, we define similarly azimuthal (east/west) and longitudinal (poleward/equatorial) heat transport by Presumably for a model with homogeneous outer boundary heat flux there will be no net azimuthal or longitudinal heat transport. Although the individual directions will be significant, on time average they will cancel out each other.

Convective structures
We plot snapshots of radial (u r ) and azimuthal flow (u / ) and heat transport along radius (j r ) and azimuth (j a ) in the equatorial plane for several combinations of convective driving Ra and forcing amplitudes q Ã in Fig. 8. For the homogeneous case the convective instability sets in at Ra c ¼ 1:6 Á 10 6 with an azimuthal wavenumber of m c ¼ 13 (Fig. 8, panel a). This is in line with the results of Hori et al. (2012), where a similar system was studied. Note, at Ra ¼ Ra c , the system is steady and hence there is only a positive contribution to the radial heat flux (j r ). The convection sets in at mid-depth. Note for this special case, the very weak azimuthal contribution j a was enhanced by a factor of 50 relative to the radial transport to visualise the structure better. For all other cases in Fig. 8, in each row j r and j a share the same contour levels.
Applying a boundary forcing with q Ã ¼ 1:0 changes the flow drastically, where an azimuthal wavenumber of m ¼ 1 is dominant. The radial and azimuthal flows (Fig. 8b) resemble the time average structures discussed before (Section 3.2). The compressed inward jet harbours already both, inward and outward radial flow. The azimuthal flows feeding into the jet start to contribute to the redistribution of heat. Especially the hot anticyclonic flow transports a large amount of heat towards the jet. For this case, the boundary induced flows are much stronger than the weak radial convection.
Increasing the convective supercriticality slightly to Ra ¼ 5 Á 10 6 (Ra=Ra c ¼ 3:13, Fig. 8c; cf. Fig. 4d) and keeping q Ã ¼ 1:0 gives rise to radial convective flows significantly participating in the flow. Especially stronger convection and hence more efficient cooling develops west of the jet, but is clearly suppressed in the opposing hemisphere. The radial outward heat transport supported by the convective flows is now strongest close to the outer boundary. The azimuthal transport is equally partitioned between the eastward advection (positive, red colours) west of the jet and westward advection (blue colours) east of the jet. That implies both azimuthal flows advect heat towards the jet.
At an even higher Ra ¼ 4 Á 10 7 (Ra=Ra c ¼ 25, Fig. 8d; cf. Fig. 4e) and q Ã ¼ 1:0 the flow has developed full time-dependence. Radial convection fills most of the left half (sectors II and III) of the shell, but seems most energetic west of the jet where the temperatures are colder and hence stronger convection can occur. The radial heat flux is outwards everywhere, but the sinking jet clearly deposits heat into the interior as well (blue) and hence heats up the areas beneath the anticyclone, which suppresses convection. The azimuthal direction is more mixed between eastward and westward, but is still mainly eastward (westward) in front of (behind) the jet.
In Fig. 8(e), we reduced the forcing amplitude to q Ã ¼ 0:5, where the jet is still clearly identifiable, though its position and most likely its power have been altered. It is quite remarkable, that the convective vigour and hence the net radial heat transport remain a strong function of azimuth, even though the flow is clearly supercritical everywhere.
As a reference case, we also add Ra=Ra c ¼ 25 with no anomalous heat flux in Fig. 8(f). Here the entire shell is filled homogeneously with turbulent convection. There is no net azimuthal transport, whereas the radial one shows the expected behaviour of dominantly outward radial heat transport (net cooling) with weak azimuthal variation.

Radial, azimuthal and longitudinal heat transport
To show the concept of the three components of j we compare two almost identical cases characterised by fixed parameters (Ra ¼ 4 Á 10 7 ), but one with homogeneous (q Ã ¼ 0) outer boundary heat flux and one with boundary forcing (q Ã ¼ 1:0). The convective supercriticality for this setup is Ra=Ra c ¼ 25 hence we expect rich nonlinear dynamics with strong time variability. Correlating the heat transport over a set of 30 individual snapshots might provide a confident mean value. Fig. 9 plots the radial, azimuthal and longitudinal (upper, middle, bottom) heat transport averaged over r and h as function of azimuthal angle.
For the homogeneous case (light blue, blue and dark-blue) the outward and inward heat transport profiles are rather flat and sum up to a net outward radial transport. This reflects the action of a fluctuating convective flow. There is no net horizontal heat flux. However, if the variable outer boundary heat flux is applied, a strong azimuthal dependence and net horizontal heat transport develops. As a reminder, the minimum (maximum) heat flux at the outer boundary is at / ¼ 0 (/ ¼ p). The large peaks in all curves for q Ã ¼ 1:0 (orange, red and dak-red) align well with the phase shifted inward jet and are not anchored where the heat flux is maximal but are azimuthally displaced by ca. 30 further east. Also, it can be clearly seen, that each of the profiles is strongly suppressed where the heat flux is smaller. The net radial heat transport features a maximum slightly westward of the jet and minimum slightly eastward of the jet. This represents nicely the enhanced cooling efficiency in the colder eastern cyclone and thermal blanketing in the hot western hemisphere. The jet diverts the azimuthal flows downwards and hence transports both hot and cold fluid to the interior (see also Fig. 8d). Note, on an azimuthal average the amount of heat transported outwards is less in the boundary forced than in the homogeneous case and hence the core shell remains hotter. The boundary forcing stimulates strong azimuthal flows towards the jet, hence we expect a significant contribution of azimuthal heat transport when q Ã > 0. The middle plot of Fig. 9 proves that both, the cyclone and the anticyclone, participate into a net azimuthal flux. It can be seen, that even in the net azimuthal heat transport (j a ), the transport is (weakly) eastwards further west of the jet (smaller /) and westward further east of the jet. Both converge into the jet. However, as the eastern cyclone is more isothermal due to the enhanced radial convective flows the anticyclonic westward transport is stronger and hence the net azimuthal flow is westward and peaks at the eastern edge of the jet. Interestingly also a net longitudinal heat transport is created by the boundary anomaly. This transport is aligned with the jet and directed polewards. As at any radial level the ageostrophic part (Eqs. (23) and (24)) dictate the strong radial inward flow of the jet to coincide with strong polarward flows, this suggests a decreasing temperature towards the poles which is found in our models (not shown).

Nusselt number vs. Rayleigh number
We expect that when convection becomes more vigourous at higher Rayleigh number the efficiency of stirring and mixing is strongly enhanced, and hence the influence of any thermal boundary inhomogeneity might diminish. We therefore calculate the net radial heat transport in terms of the Nusselt number (Eq. (37)) for cases between weakly nonlinear with Ra ¼ 10 7 up to strongly turbulent Ra ¼ 6:4 Á 10 8 , ranging in terms of convective supercriticality from 6:25 6 Ra=Ra c 6 400. Averages are taken over radius and longitude hence the j r ð/Þ-curves are similar to the top plot in Fig. 9. The Rayleigh numbers used for Fig. 10 are stepwise increased by factor four and q Ã ¼ 1 is fixed. Note, for visualisation purposes the thick, bright curves are smoothed using gnuplot's standard Bezier smoothing function, whereas the darker and thinner curves give the real data. Interestingly, the enhanced radial convection between / ¼ ½p=2 . . . 3p=2 (or in sectors II and III) and the suppression on the opposite half seems quite independent of the Rayleigh number. Also for all cases the convection is stronger (weaker) west of the jet, and weaker further east. The independence on the convective vigour Ra h is a quite remarkable result, as for rather high Ra=Ra c the flow is almost supercritical everywhere.
The suppression of the net radial heat transport in sectors I and IV should be visible as well in the Nusselt numbers. We calculate Nu with Eq. (37) for a set of 13 different Rayleigh numbers for the homogeneous reference cases Nu 0 (Fig. 11, red), for the equivalent boundary forced models Nu 1 (black) and their ratio (blue). The smallest is the critical Rayleigh number for the homogeneous reference case. By definition Nu 0 should be marginal there (black curve), but the boundary forcing actually enhances the heat transport (red). It is expected by linear theory that at the onset any boundary anomaly will enhance Nu (Kelly and Pal, 1978). However, increasing the Rayleigh number towards more nonlinear systems shows that at any supercritical test case the homogeneous case (Nu 0 ) features a higher Nusselt number or more efficient heat transport than the boundary forced case (Nu 1 ). The suppression by the thermal boundary anomaly seem to settle at ca. 50% for strongly driven models. Given the sparse data, we do not attempt π/2 π 3π/2 2π azimuth φ j r Ra=1e7 4e7 1.6e8 6.4e8 to derive a scaling relation like Nu $ Ra a . However, it is well accepted that a is smaller than unity for weakly rotating systems and depends in the model setup etc. E.g., King et al. (2010) suggest a ¼ 2=7 for a weakly rotating spherical shell with fixed temperature contrast and a dynamo process, whereas a comparable model but with rapid rotation follows a ¼ 6=5. We find the rapid rotating behaviour only for models very close to the onset of convection.
For the majority of our models a % 0:45 seems reasonable, when Ra h ) Ra c as indicated in Fig. 11. The reason for the suppression of convective efficiency by the boundary forcing is then directly related to that exponent. As a toy model of the boundary forced cases, we assume the Rayleigh number, which is determined by the radial temperature gradient, shall be doubled in one half of the core and set to zero in the other one. This leads to different values of the globally-averaged radial heat transport, i.e. the Nusselt number. Compared to a model with homogeneous boundary heat flux, the heat transport can be given as yielding a reduction of heat transport in the boundary forced case for 0 < a < 1. Note, the suppression measured in the numerical models is much stronger than this simplistic estimate (% 1:46 for a ¼ 0:45) as the ratio Nu 0 =Nu 1 seems to be bordered by 1:8 and 2:5. Even though areas with reduced local q are exactly compensated with enhanced q areas and hence on zonal average the total heat flux is identical between a ðq Ã ¼ 0:0Þ-and a (q Ã ¼ 1:0)-model, the boundary forcing reduces the net radial heat transport independently of the convective vigour. Physically this implies, that the additional azimuthal heat transport introduced by the boundary anomaly consumes a significant fraction of the available convective kinetic energy, hence the mean radial heat transport is reduced.

Dependence on heat flux anomaly amplitude
As a final set of models, we keep all parameters fixed (Ra ¼ 4 Á 10 7 ; Ra=Ra c ¼ 25; E ¼ 10 À4 ) and vary q Ã . Such a procedure will test the sensitivity of a vigourously convecting system regarding inhomogeneities at the CMB heat flux. Starting with the azimuthal profiles of the net radial heat flux j r in Fig. 12, we conduct numerical experiments for a few other q Ã -values. Note, only the model related to q Ã ¼ 1:0 (black curve) features a neutrally critical heat flux at / ¼ 0. However, the smooth transition between q Ã ¼ 0 (grey) and q Ã ¼ 1:0 (black) suggests a rather linear dependence on q Ã . Note, that the position of the jet indicated by the rough decay of j r is shifting towards larger U and the amplitude is weaker with decreasing q Ã . For measuring the azimuthal asymmetry of the radial heat transport we define a hemisphericity like where the maximal j r max ¼ max½j r is taken between p=2 6 / 6 3p=2 (sectors II and III) and the minimal j r min ¼ min½j r from the other half. Hence we suggest that the effects on the convection do not depend at all on the convective vigour and only linearly on the boundary forcing amplitude.

Conclusions
We numerically investigated core flows under inhomogeneous cooling at the outer boundary. When the core is subcritical to radial thermal convection, these thermal anomalies induce core flows whose linear and nonlinear properties are discussed in detail. We also study the situation when, in contrast, the core is unstable to thermal convection and a non-homogeneous CMB heat flux will strongly alter the convective structures and the heat transport. To model the horizontal variation of the lower mantle temperature in a terrestrial planet, we apply a perturbation of the CMB heat flux shaped as a spherical harmonic Y 11 and perform a comprehensive study with a numerical model of an incompressible and viscous fluid heated from within and cooled from above, which is contained in a rapidly rotating spherical shell. Such a simple heating setup takes only secular cooling into account and hence represents best the liquid cores of the Early Earth, Mars or any terrestrial exoplanet before inner core solidification. As the main control parameters we conduct a parameter survey regarding the Rayleigh number Ra and forcing amplitude q Ã . In the absence of radial convection (Ra h ¼ Raq Ã < Ra c ), a linear conductive regime (determined by the vorticity balance between heat conduction and Coriolis effects and the ageostrophic thermal wind) and a nonlinear advective regime controlled by thermal advection are identified. If Ra h > Ra c radial convection participates strongly in the nonlinear dynamics of a convective regime.
In general such a heat flux anomaly introduces a global azimuthal and latitudinal temperature gradient, which tends to be equilibrated by a dominantly geostrophic mean two-cell circulation given by a broad cyclonic structure in the eastern hemisphere and an anticyclone in the west, which converge into a radial inward flow. As it was previously suggested and in agreement with the linear theory, see Section 3.1 in this paper (or Zhang and Gubbins, 1992;Yoshida and Hamano, 1993) the radial flows are not located where the outer boundary heat is maximal, but where the azimuthal temperature gradient is maximised. For the linear model this phase shift is 90 eastward or a quarter of the azimuthal wavelength of the anomaly. This can be understood as the Sverdrup relation, which is wellknown in the context of geophysical fluid dynamics.
Starting from a non-convecting, and almost linear numerical model (Ra h ¼ 10), the numerical findings confirm the analytical predictions of the induced geostrophic and ageostrophic flows. For higher Ra and once the advective terms become dynamically important, the radial inflow is compressed into a jet, bended in prograde direction and moves retrograde to smaller phase shift angles. The nonlinearities in the system heat up the western and cool the eastern side of the core. It might be stated, in the absence of convection any inhomogeneously cooled system drives baroclinic flows, whose nonlinear properties, such as the effect of advection or heat transport were analysed here.
If the driving is enhanced beyond the critical value for the onset of convection in the equivalent homogeneous model (Ra h > Ra c ), radial convective flows appear and are enhanced where the heat flux is increased, and strongly suppressed where the heat flux is weaker than average. This convective asymmetry is a consistent property, emerges independent of the convective vigour and increases linearly with the heat flux anomaly amplitude q Ã . It could be also shown, that the convection is most energetic west of the inward jet where cooling is most efficient.
We also investigated the properties of heat transport, and find that the boundary anomaly introduces a net azimuthal heat trans-port towards the descending jet. As a consequence of the laterally variable cooling efficiency, the Nusselt number representing the global radial heat transport is always smaller for boundary forced than for homogeneously cooled models. Implying that the mean radial convective heat transport is globally reduced by at least 50%, although the azimuthal averaged heat flux is the same in both models. The amount of kinetic energy available for radial convection is reduced by the boundary driven azimuthal flows and the inward jet. As a consequence, the liquid iron cores of terrestrial planets facing inhomogeneous cooling might be hotter and much less well mixed. In addition core convection and the potential induction of global magnetic field might last longer. Remarkably, the local suppression of core convection is rather independent of convective vigour (Rayleigh number), but is linearly amplified by the relative strength of the thermal CMB inhomogeneity.
It seems contradictory, that in a series of comparable studies (Sreenivasan, 2009;Aurnou and Aubert, 2011;Dietrich and Wicht, 2013), the magnetic field induction process was amplified or triggered by CMB heat flux anomalies, whereas our results suggest a strong suppression of core convection for any anomalous heat flux. It has been shown in our models that the boundary forcing drives additional flows favouring a supercritical dynamo, where it would fail for a homogeneous system. The local suppression and amplification of convective flows in response to the weaker or stronger CMB heat flux will only lead to an azimuthal variation of the magnetic field intensity, but not terminate the dynamo. Indeed the local enhancement makes a dynamo more likely. It might be also said, that equatorially antisymmetric and axisymmetric heat flux anomalies drive fierce axisymmetric flows which provide shearing hence strengthen the dynamo process.
Our models are applicable to terrestrial exoplanets orbiting their host star in synchronous rotation. Assuming a similar frontto-back thermal anomaly at the CMB as imprinted on the surface by the irradiation of the host star, the Sverdrup balance yields significant azimuthal flows on the far side converging into a prograde spiralling inward jet. If the liquid core is further convecting, a potential magnetic field induction process might be strongly concentrated on the far side of the planet.
Compared to experimental works by Sumita and Olson (1999) our models reproduce several key features, such as the development of a sharp temperature front located east of the heat flux anomaly and hosting the radial inward spiralling jet and separating the cold east from the hot west. Furthermore the experiments reported a scaling relation suggesting the amplitude boundary induced flows increases like the square root of the total anomalous heat flux through the anomaly (Sumita and Olson, 2002), which is confirmed within our numerical results. This implies the very efficient turbulent stirring and mixing due to the regular radial convection in the core might not suppress the influence of boundary anomalies.