  Convection in the metallic cores of terrestrial planets is likely to be
subjected to lateral variations in heat flux through the outer boundary imposed
by creeping flow in the overlying silicate mantles. Boundary anomalies can
significantly influence global diagnostics of core convection when the Rayleigh
number, $Ra$, is weakly supercritical; however, little is known about the
strongly supercritical regime appropriate for planets. We perform numerical
simulations of rapidly rotating convection in a spherical shell geometry and
impose two patterns of boundary heat flow heterogeneity: a hemispherical
$Y_1^1$ spherical harmonic pattern; and one derived from seismic tomography of
Earth's lower mantle. We consider Ekman numbers $10^{-4} \le E \le 10^{-6}$,
flux-based Rayleigh numbers up to $\sim800$ times critical, and Prandtl number
unity. The amplitude of the lateral variation in heat flux is characterised by
$q_{\rm L}^*=0, 2.3, 5.0$, the peak-to-peak amplitude of the outer boundary
heat flux divided by its mean. We find that the Nusselt number, $Nu$, can be
increased by up to $\sim25$\% relative to the equivalent homogeneous case due
to boundary-induced correlations between the radial velocity and temperature
anomalies near the top of the shell. The $Nu$ enhancement tends to become
greater as the amplitude and length scale of the boundary heterogeneity are
increased and as the system becomes more supercritical. This $Ra$ dependence
can steepen the $Nu \propto Ra^\gamma$ scaling in the rotationally dominated
regime, with $\gamma$ for our most extreme case approximately $20\%$ greater
than the equivalent homogeneous scaling. Therefore, it may be important to
consider boundary heterogeneity when extrapolating numerical results to
planetary conditions.
