Mantle-induced temperature anomalies do not reach the inner core boundary

SUMMARY Temperature anomalies in Earth’s liquid core reﬂect the vigour of convection and the nature and extent of thermal core–mantle coupling. Numerical simulations suggest that longitudinal temperature anomalies forced by lateral heat ﬂow variations at the core–mantle boundary (CMB) can greatly exceed the anomalies that arise in homogeneous convection (i.e. with no boundary forcing) and may even penetrate all the way to the inner core boundary. However, it is not clear whether these simulations access the relevant regime for convection in Earth’s core, which is characterized by rapid rotation (low Ekman number E ) and strong driving (high Rayleigh number Ra ). We access this regime using numerical simulations of nonmagnetic rotating convection with imposed heat ﬂow variations at the outer boundary (OB) and investigate the amplitude and spatial pattern of thermal anomalies, focusing on the inner and outer boundaries. The 108 simulations cover the parameter range 10 − 4 ≤ E ≤ 10 − 6 and Ra = 1 − 800 times the critical value. At each Ra and E we consider two heat ﬂow patterns—one derived from seismic tomography and the hemispheric Y 11 spherical harmonic pattern—with amplitudes measured by the parameter q (cid:2) = 2.3, 5 as well as the case of homogeneous convection. At the OB the forcing produces strong longitudinal temperature variations that peak in the equatorial region. Scaling relations suggest that the longitudinal variations are weakly dependent on E and Ra and are much stronger than in homogeneous convection, reaching O (1) K at core conditions if q (cid:2) ≈ 35. At the inner boundary, latitudinal and longitudinal temperature variations depend weakly on Ra and q (cid:2) and decrease strongly with E , becoming practically indistinguishable between homogeneous and heterogeneous cases at E = 10 − 6 . Interpreted at core conditions our results suggest that heat ﬂow variations on the CMB are unlikely to explain the large-scale variations observed by seismology at the top of the inner core.


I N T RO D U C T I O N
Convection in Earth's liquid core sustains Earth's magnetic field through a dynamo process that converts kinetic energy into magnetic energy. The convection is driven by thermal and chemical buoyancy forces and the associated thermochemical anomalies contain important information regarding the operation of the dynamo and the dynamics of the core. Radial variations reflect the vigour of core convection and the relative strength of thermal and chemical driving (Lister & Buffett 1995). Seismically slow regions observed at the top and bottom of the core (Souriau & Calvet 2015) have been argued to reflect departures from uniform composition Helffrich & Kaneshima 2010;Gubbins & Davies 2013;Brodholt & Badro 2017;Wong et al. 2018), though thermal effects inevitably explain some of the signal. Lateral variations, particularly in longitude, may reflect coupling of core convection to temperature variations in the lowermost mantle (Buffett 2007). A key question is whether lateral variations in temperature imposed by mantle structure at the core-mantle boundary (CMB) persist to the inner core boundary (ICB) where a striking hemispheric variation in seismic properties is observed (Aubert et al. 2008;Monnereau et al. 2010;Gubbins et al. 2011;Souriau & Calvet 2015). The main obstacle in assessing this hypothesis is that the thermal structure of the core cannot be imaged directly and is difficult to disentangle from chemical variations. In this paper we use numerical simulations to investigate the magnitude and pattern of thermal anomalies in the core induced by heat flow variations at the CMB.
Radial temperature variations in the core are predominantly due to adiabatic compression (Jones 2015). The superadiabatic temperature difference T s between the CMB and ICB is usually assumed to be accommodated by thermal boundary layers at the top and bottom of the core that are a matter of metres thick, with an adiabatic bulk interior (Jones 2015). Taking a superadiabatic temperature gradient across these layers comparable to the adiabatic gradient of 1 K km −1  and assuming an adiabatic bulk gives T s ∼ 10 −2 −10 −3 K. Numerical simulations of spherical shell Mound & Davies 2017) and plane layer (Julien et al. 2012a,b) rotating convection display non-zero interior temperature gradients, which can account for as much as half of the total superadiabatic temperature drop in the plane layer case, though this does not significantly affect the estimate of T s . Variations in composition with depth are harder to estimate and are usually neglected altogether (Jones 2015).
With uniform thermal boundary conditions (referred to here as 'homogeneous' convection) lateral variations in temperature and composition within the core are expected to be tiny. Stevenson (1987) estimated that the density fluctuations associated with core convection are 9 orders of magnitude smaller than the mean density. Bloxham & Gubbins (1987) assumed a thermal wind balance just below the CMB and found that temperature anomalies of O(10 −3 ) K could account for core flows of O(10) km yr −1 (see also Bloxham & Jackson 1990). These thermal anomalies are 6-7 orders of magnitude smaller than estimates of the adiabatic temperature at the CMB  and comparable to the estimate of T s . Based on estimates of the respective buoyancy fluxes, compositional buoyancy can exceed the thermal buoyancy that drives the core flow (Lister & Buffett 1995), though the magnitude of the chemical anomalies will be about 5 orders of magnitude smaller than their thermal counterparts owing to the much larger compositional expansivity (Gubbins et al. 2004).
Heterogeneous convection arising due to lateral variations in heat flow imposed on the core by the mantle may significantly alter the estimates above. Seismic tomographic images of the lowermost mantle clearly identify two regions beneath Africa and the central Pacific characterized by seismic velocity anomalies that are several per cent slower than predicted from 1-D models (e.g. Masters et al. 1996;Grand et al. 1997;Garnero & McNamara 2008;Hernlund & McNamara 2015). The precise nature of these anomalies is still debated, but there is general consensus that they at least partly reflect thermal variations in the lowermost mantle (Gubbins 2003;Garnero et al. 2016). Since thermal anomalies in the mantle greatly exceed those in the core the former produce lateral variations in the heat flow across the CMB with seismically slow regions corresponding to locally depressed heat flow and fast regions corresponding to locally elevated heat flow (e.g. Gubbins 2003;Olson 2003). The amplitude of the heat flow anomalies cannot be observed directly, but some numerical simulations suggest that the lateral variations exceed the average CMB heat flow (Nakagawa & Tackley 2007). However, unlike homogeneous convection, even infinitesimal boundary forcing will drive thermal wind flows that transport heat both laterally and into the interior (e.g. Zhang & Gubbins 1992;Shishkina et al. 2016).
In non-magnetic rotating convection the effect of imposed boundary forcing depends on the forcing pattern and amplitude, the Rayleigh number Ra measuring the vigour of homogeneous convection, the Ekman number E measuring the ratio of viscous to rotational effects, and the Prandtl number Pr measuring the ratio of viscous and thermal diffusion. In this paper we consider two largescale heat flow patterns: (1) the pattern inferred by assuming that shear wave velocity variations in the lowermost mantle reflect thermal heterogeneity (the 'tomographic' pattern), which is dominated by spherical harmonic degree and order 2 (Masters et al. 1996); (2) the hemispheric spherical harmonic Y 1 1 pattern, which has been advocated as the large-scale pattern during supercontinent formation (Zhong et al. 2007). At low Ra, moderate E ∼ 10 −3 − 10 −5 and Pr 1 these forcings produce quasi-periodic solutions with convection rolls clustered under high boundary heat flow regions instead of the usual pattern of drifting uniformly spaced rolls (Davies et al. 2009). Time-dependent solutions emerge with decreasing E (Davies et al. 2009), decreasing Pr (Zhang & Gubbins 1996), increasing Ra (Sun et al. 1994) or the addition of a magnetic field Sreenivasan 2009). As Ra increases the boundary effects become harder to identify in snapshots, though they can be clearly seen in time-averages (Sun et al. 1994;Olson & Christensen 2002;Aubert et al. 2007) and we will hence use averages throughout this paper. The main effect of large-scale boundary forcing in the parameter regime studied to date is to induce longitudinal variations in the time-averaged core surface flow and magnetic field (e.g. Bloxham 2000;Christensen & Olson 2003;Olson 2003;Davies et al. 2008;Mound et al. 2015;Olson et al. 2017). Longitudinal thermal structure, the focus of this paper, has received less attention though a recent study predicts temperature variations of O(1) K at the CMB in geodynamo simulations with strong boundary forcing and imposed stable stratification (Christensen 2018), much larger than estimates based on homogeneous convection.
The influence of lateral CMB heat flow variations may not be confined to the uppermost regions of the core. Aubert et al. (2008), Gubbins et al. (2011) and  found that columnar boundary-induced flows would leave an imprint on the inner core, producing localized melting and freezing at the ICB. This mechanism has been invoked as a possible explanation for the laterally heterogeneous structure observed at the top of the inner core (Niu & Wen 2001;Waszek et al. 2011;Souriau & Calvet 2015). In this paper, we will quantify the strength and structure of thermal variations at the ICB across a broad range of model parameters.
Numerical simulations of core convection and dynamo action cannot model the small diffusivities that describe molecular transport processes in Earth's core. Systematic studies, rather than individual model runs, are needed to establish the effect of changing control parameters on the time-averaged properties. This is particularly challenging when considering heterogeneous boundary conditions as extra parameters (the pattern and amplitude of the lateral variations) must be varied compared to homogeneous models. The problem is further exacerbated since recent work on non-magnetic rotating convection shows that models with E ≥ 10 −4 , which represents the majority undertaken to date, do not access the rapidly rotating regime thought relevant to Earth's core Mound & Davies 2017). It is not currently practical to comprehensively study this regime with geodynamo simulations owing to the enormous computational costs (e.g. Matsui et al. 2016), but it can be accessed with non-magnetic models. A systematic study of nonmagnetic boundary-forced convection has been achieved at E = 5 × 10 −5 (Dietrich et al. 2016), but using a hemispheric heat flow pattern that is not directly related to present-day Earth.
We have recently produced a suite of 108 non-magnetic rotating convection simulations spanning the parameter range 10 −4 ≤ E ≤ 10 −6 , Ra c ≤ Ra < 800Ra c and Pr = 1. Here Ra c is the critical value of Ra for the onset of non-magnetic convection. At each E and Ra we have conducted simulations with hemispheric, tomographic and homogeneous outer boundary (OB) heat flow patterns. Of the 108 simulations, 106 are listed in Appendix B of Mound & Davies (2017) while new models with E = 10 −6 , Ra = 18 000, and a Y 1 heat flow pattern with two different forcing amplitudes have been conducted for this study. Heat transfer data suggest that the models with E = 10 −4 transition with increasing Ra from the weakly nonlinear regime involving a balance between viscous, Coriolis and buoyancy forces (Gillet & Jones 2006), directly to a regime where rotation plays a subdominant role Mound & Davies 2017). However, crucially, models at E = 10 −5 and 10 −6 reach the rapidly rotating and strongly driven regime thought to describe core dynamics.
In this paper, we use the suite of models in Mound & Davies (2017) to quantify the magnitude and spatial structure of the timeaveraged temperature, focusing on conditions at the top and bottom of the domain. We quantify latitudinal and longitudinal thermal variations and develop empirical scaling relations to gain insight into the thermal conditions that may pertain at the top and bottom of Earth's core. Model setup and outputs are described in Section 2 and results are presented in Section 3. A discussion of the limitations of our model and considerations required to apply the results to Earth are given in Section 4. Conclusions are presented in Section 5.

Numerical model and parameters
Complete details of the numerical simulations used in this study can be found in Mound & Davies (2017) and so only a brief description is given here. The numerical code ) employs spherical polar coordinates (r, θ , φ) and solves the equations governing conservation of mass, momentum and energy in a spherical shell that rotates about the verticalẑ direction with constant angular frequency . The dimensionless forms of these equations can be written, respectively, as where u is the fluid velocity,P is the modified pressure, and T = T C + T is the total temperature, where T C is the temperature in the absence of motion and T is the temperature perturbation. In writing the dimensionless equations we have assumed a constant kinematic viscosity ν, thermal diffusivity κ and thermal expansivity α, and scaled length by the shell thickness h, time by the thermal diffusion time τ d = h 2 /κ, and temperature by β/h. T C satisfies the equation ∂T C /∂r = −β/r 2 and hence the parameter β is related to q ave , the mean heat flow per unit area at r o , by q ave = −k∂ T C /∂r | ro = kβ/r 2 o where k is the thermal conductivity. The dimensionless parameters in eqs (1)-(3) are the Prandtl number Pr, the Ekman number E and the Rayleigh number Ra, which are defined as In this study all models have Pr = 1, a ratio of inner boundary (IB) and OB radii r i /r o = 0.35 and a gravity profile that varies linearly with r such that g = −(g/r o )r. The critical Rayleigh number depends on E: Ra c = 16.4 at E = 10 −4 , 24.7 at E = 10 −5 and 41.0 at E = 10 −6 (Mound & Davies 2017).
In all simulations the velocity boundary conditions are no-slip and fixed heat flux at r i and r o . There are three groups of simulations, defined by the pattern of heat flow heterogeneity imposed on the OB. The first, denoted by the letter 'H', have homogeneous heat flow at r o . The second, denoted by 'Y', have an imposed pattern corresponding to the spherical harmonic degree and order 1, Y 1 1 . The final group, denoted 'T', have an imposed heat flow pattern derived from the seismic tomography model of Masters et al. (1996) assuming that the lateral shear-wave velocity variations reflect thermal heterogeneity in the lowermost mantle. The tomographic pattern is dominated by the spherical harmonic Y 2 2 component, but also contains other significant contributions such as Y 0 2 . The amplitude of the lateral variations is defined by the parameter q : where q max , q min and q ave are, respectively, the maximum, minimum and average heat flow through the OB. Mound & Davies (2017) considered values of q = 2.3 and 5.0 (homogeneous cases correspond to q = 0). We use the shorthand notation Xq = C to distinguish different model groups, where X ∈ {H, Y, T} and C = 0, 2.3, 5.0. Our numerical code utilizes a toroidal-poloidal decomposition to represent the vector velocity field. Scalar fields are expanded in spherical harmonics on spherical surfaces and finite-differences in radius. Fields are evolved in time using a predictor-corrector scheme with adaptive timestepping. The code successfully reproduces the dynamo benchmark solutions Davies et al. 2011;Matsui et al. 2016). All solutions used in this study have been checked for spatial convergence, and run until robust time-averaged diagnostics of heat transfer behaviour were obtained (Mound & Davies 2017).

Diagnostics of temperature variations
We use several diagnostics to quantify latitudinal and longitudinal temperature variations at the IB and OB. Strong latitudinal variations can occur without heterogeneous boundary forcing, reflecting the efficiency of heat transport as well as the comparative vigour of convection inside and outside the tangent cylinder (Jones 2000). Longitudinal variations are expected to be weak in (time-averaged) homogeneous convection, but may be promoted by an imposed pattern of heat flow. Any influence of heterogeneous boundary conditions is expected to be most evident in time-averages, which we define by an overbar: where t 0 and t 1 are the start and end times, which usually span over 100 advection times (see Mound & Davies 2017, for details).
Horizontal averagesT h are defined as We characterize latitudinal variations inT at each r and θ bȳ The peak-to-peak amplitude at any radius is This is often, but not always, the pole-to-equator difference.
Downloaded from https://academic.oup.com/gji/article-abstract/218/3/2054/5513450 by University of Leeds user on 05 July 2019 We characterize longitudinal variations inT at each r and θ bȳ The largest longitudinal variation at any radius is When calculating differences between homogeneous and heterogeneous models the superscripts 'hom' and 'het' will be used, respectively; these will generally be omitted when there is no possibility of confusion.
Estimating the dimensional strength of temperature anomalies requires a temperature scale that can be observationally constrained.
Here we use T, the difference betweenT h on the OB and IB. We define as a dimensionless measure of the boundary temperature anomalies at the OB. An analogous expression is used at the IB.
In spherical shell convection with fixed flux boundaries T is an output and varies significantly between models. In this configuration T is related to the Nusselt number Nu by  (Mound & Davies 2017) where T c is the temperature difference between r i and r o due to conduction alone, that is when the velocity is zero. Therefore, an increase in the efficiency of convective heat transport corresponds to a decrease in T.

R E S U LT S
In this section, we consider properties of the time-averaged tem-peratureT at the IB and OB. We compare the latitudinal and longitudinal variations inT between models with homogeneous (H), tomographic (T) and Y 1 1 (Y) boundary patterns. We begin by presenting example cases with high Ra and low E before considering systematic behaviour as a function of Ra and E. Fig. 1 shows examples of the time-averaged temperature pertur-bationT −T h at the IB and OB for H, T and Y boundary patterns at E = 10 −5 and Ra = 526Ra c . Considering first the OB, the poles are anomalously hot compared to the mid-latitudes in all models since convection is suppressed inside the tangent cylinder (the imaginary cylinder aligned with the rotation axis and tangent to the inner core equator, Jones 2000), though the effect weakens as q increases. In the homogeneous modelT is dominantly axisymmetric and the equatorial region is anomalously cool, while in the heterogeneous modelsT clearly reflects the pattern of the boundary heat flow, with stronger anomalies at higher values of q as expected. At the IB T in the homogeneous model is similar to that at the OB, while in the tomographic modelsT is dominantly axisymmetric with little  evidence of the Y 2 2 component seen at the OB. In these models the southern hemisphere is slightly hotter than the northern hemisphere (see also Figs 2 and 3), which reflects the slight displacement of the low heat flow regions on the OB towards the southern hemisphere. In the Y 1 1 models the pattern of OB heat flow can be distinguished in the equatorial region of the IB at these parameter values. Overall, the effect of the heterogeneous boundary conditions is more apparent at the OB than the IB.
To quantify latitudinal and longitudinal temperature variations at high Ra we plotT θ andT φ for models with different E and similar Ra/Ra c . Figs 2 and 3 show results at the OB and IB, respectively, for H, Tq = 5 and Yq = 5 models. Solutions with E = 10 −4 display several features not seen at lower E. First, homogeneous E = 10 −4 solutions display strong equatorial asymmetry, particularly at the OB, whereas at lower E,T θ is strongly equatorially symmetric (Fig. 2a). Secondly, homogeneous E = 10 −4 models display much stronger longitudinal variations inT φ at all latitudes on the OB than the lower E models (Fig. 2d). Thirdly, heterogeneous E = 10 −4 solutions display at equatorial maximum inT θ at the IB rather than a minimum seen at lower E (Figs 3a-c). These differences may reflect the fact that our E = 10 −4 models mainly access the weakly nonlinear and the non-rotating convective regimes, while the E = 10 −5 and 10 −6 models mainly access the rapidly rotating regime (cf. Gastine et al. 2016;Mound & Davies 2017).
At the OB, higher values of q yield greater departures from the homogeneous case, with tomographic cases exhibiting the highest degree of equatorial asymmetry (Fig. 2), reflecting the imposed heat flow pattern. Thermal anomalies are enhanced in the equatorial region relative to the homogeneous case as expected. Latitudinal variations inT θ show little dependence on E, while the longitudinal variations are strongly dependent on E. The main effect of the heterogeneous boundary condition is to produce a hotter equator and colder poles at the OB compared to the homogeneous case. In the tomographic cases this arises because of the anomalously low heat flow regions under Africa and the Pacific, which suppress convection, and higher than average heat flux at higher latitudes. In the Y 1 1 cases this occurs because downwellings induced under the high heat flow hemisphere are narrower than upwellings induced under the low heat flow hemisphere (Mound & Davies 2017). At the IB, several contrasting features to the OB results can be noted. First, departures from the homogeneous case are greatly reduced in all models and particularly in the E = 10 −6 cases where the heterogeneous and homogeneous models are barely distinguishable (Fig. 3). Second, the heterogeneous solutions are dominantly equatorially symmetric, even at E = 10 −4 , for all forcing patterns. Third, latitudinal and longitudinal temperature variations depend strongly on E. In the higher E cases the heterogeneous boundary condition elevatesT φ in the equatorial region, but even this signature is strongly reduced at E = 10 −6 . In our low E and high Ra runs there are no morphological characteristics of the IB temperature that allow the homogeneous and heterogeneous cases to be meaningfully distinguished.
Temperature differences between the IB and OB reflect the efficiency of convective heat transfer and the thermal homogenization of the system. Fig. 4 showsT θ (r i ) −T θ (r o ) as a function of latitude for the low E models in Figs 2 and 3. The mean temperature at each radius has been subtracted and soT θ (r i ) −T θ (r o ) < 0 shows enhanced heat transport (reduced temperature difference). In the homogeneous modelsT θ (r i ) −T θ (r o ) in the polar regions decreases with increasing Ra/Ra c , which probably reflects the steep increase in heat transfer with Ra associated with convection inside the tangent cylinder (Yadav et al. 2015). Compared to the homogeneous case, heterogeneous models have lowerT θ (r i ) −T θ (r o ) in the equatorial region, reflecting enhanced heat transfer, and higherT θ (r i ) −T θ (r o ) in the mid-and high-latitude regions, reflecting reduced heat transfer. Increased heat transfer at low latitudes in heterogeneous simulations arises because the boundary heat flow variations promote strong correlations between the radial velocity and temperature (and hence the advective heat transport) which reduces the temperature gradient and hence T (Mound & Davies 2017). At high-latitudes the imposed heat flow patterns contain much weaker anomalies and their effect is correspondingly reduced.
The maximum variation ofT with longitude, δT φ , is compared to the maximum variation with latitude, δT θ , at the IB and OB in Fig. 5. Considering first the OB, δT φ /δT θ is larger by at least an order of magnitude in heterogeneous compared to homogeneous cases and exceeds 1 in all heterogeneous cases when Ra/Ra c ≥ 100 (Fig. 5a). In general δT φ /δT θ decreases with decreasing E, increases with q and increases weakly with Ra at high Ra (Fig. 5a). The variability in δT φ /δT θ between the different imposed patterns and forcing amplitudes does not seem to depend significantly on E (Fig. 5c).
Significant differences are seen in δT φ /δT θ at the IB compared to the OB. First, δT φ /δT θ differs by only a factor of 2-5 between homogeneous and heterogeneous cases and is smaller than 1 in the majority of cases (Fig. 5b). The E = 10 −4 models are again exceptional, producing larger longitudinal temperature variations across a broad range of Ra/Ra c . Secondly, and most importantly, δT φ /δT θ strongly decreases with decreasing E as does the variability between different imposed patterns (Fig. 5d). At E = 10 −6 latitudinal variations at the IB are stronger than longitudinal variations by an order of magnitude and do not depend significantly on the pattern or amplitude of the boundary forcing.
Another measure of the thermal heterogeneity induced by the OB anomalies is to compare the maximum longitudinal variation inT , δT φ , for heterogeneous and homogeneous cases. At the OB the ratio δT het φ /δT hom φ ≈ 2 − 6 ( Fig. 6a), while at the IB δT het φ /δT hom φ ≈ 1 − 2 independent of Ra/Ra c (Fig. 6b). Crucially, longitudinal temperature variations at the IB are almost identical between homogeneous and heterogeneous models at E = 10 −6 independent of the pattern and amplitude of boundary heterogeneity within the ranges considered (Fig. 6d).

D I S C U S S I O N
Our We have noted that homogeneous solutions with E = 10 −4 display much stronger latitudinal and longitudinal variations than models at lower E and heterogeneous E = 10 −4 models display maximum IB temperature at the equator rather than a minimum as seen at lower E. These differences are important when interpreting morphological characteristics of the temperature fields in these models since such features evidently do not persist into the lower E regime that is more appropriate to planetary cores.
Before considering the potential consequences of our results for Earth's core we must evaluate limitations in the modelling approach. Our suite of models is one of the largest created to study boundary heat flow variations in spherical shell rotating convection that systematically spans (E, Ra)-space, reaching values near the cutting edge of what is presently possible with direct numerical simulations. On the other hand it is still impossible to access the values of E ∼ 10 −13 − 10 −16 and Ra ∼ 10 17 estimated for the core (e.g. Mound et al. 2019) based on molecular properties of iron at high pressure and temperature (e.g. Davies et al. 2015) and so extrapolation from simulations is inevitable.
The value of q is also highly uncertain. Estimates of mantle thermal conductivity coupled with seismic tomography suggest that the CMB heat flow ranges from 0 to 140 mW m −2 (Stackhouse et al. 2015), which is consistent with predictions from some global mantle circulation models (Nakagawa & Tackley 2007;Olson et al. 2015). The mean superadiabatic heat flow at the CMB q ave depends on the difference between the total CMB heat flow, which is estimated to lie between 5 and 17 TW (Lay et al. 2009;Nimmo 2015;Jaupart et al. 2015), and the adiabatic heat flow, which is 4−16 TW depending on the assumed value of the core thermal conductivity . Estimates of q ave , and hence q , can therefore take either sign, and as a consequence of its definition q becomes infinite if the top of the core is neutrally stable. An alternative estimate of q ave is obtained by balancing inertial, Coriolis and buoyancy terms in the vorticity equation. The resulting estimates of the flow speed are supported by dynamo simulations (Christensen & Aubert 2006), which translate into a value of q ave ∼ 4 mW m −2 (Jones 2011). Using this estimate together with the values of q max − q min discussed above implies that q could be as large as 35. By comparison, numerical studies have tended to focus on lower values: the regime q ≤ 1 has been explored in detail (e.g. Zhang & Gubbins 1993;Aubert et al. 2007Aubert et al. , 2008Davies et al. 2009;Dietrich et al. 2016), while some studies have considered values of q up to approximately 2 (Olson & Christensen 2002;Sreenivasan 2009;Sahoo & Sreenivasan 2017). Laboratory experiments have studied strong boundary forcing with q max /q ave ≤ 95 (Sumita & Olson 2002).
At the OB, Figs 5 and 6 show that δT φ /δT θ and δT het φ /δT hom φ display a weak dependence on Ra and E and increase with q and the lengthscale L of the imposed heat flow pattern. This is consistent with previous studies that have found the effect of boundary anomalies tends to strengthen as q and L increase (Zhang & Gubbins 1993;Davies et al. 2009;Calkins et al. 2015). Since we expect q to be larger in Earth than in our models, the results suggest that longitudinal temperature variations would exceed latitudinal temperature variations and be clearly distinguishable from temperature anomalies due to homogeneous convection at core conditions.
By contrast, at the IB, Figs 5 and 6 show that δT φ /δT θ and δT het φ /δT hom φ are almost independent of Ra and L at high Ra, increase with q (though to a lesser extent than at the OB) and decrease strongly as E decreases. At low E and high Ra it is clear that latitudinal temperature variations dominate longitudinal variations at the IB and that the longitudinal variations induced by the imposed heat flow are indistinguishable from those that arise due to homogeneous convection. Since the E effect dominates in our simulations we would expect these conclusions to be reinforced upon moving to the lower E regime that characterizes Earth's core. Sumita & Olson (2002) observed large-scale spiralling jets that reached the IB at somewhat similar conditions to those used here. They conducted rotating convection experiments in a hemispherical shell with E = 2.4 × 10 −6 , Ra T /Ra c ≤ 62 and q max /q ave ≤ 95 and applied heterogeneous forcing through localized heat sources on the OB. Here Ra T is a Rayleigh number based on an imposed temperature difference, rather than the flux-based Rayleigh number used in the present simulations. Sumita & Olson (2002) found that the radial scale of the jets was sensitive to the location and pattern of boundary heating; in particular, they observed that jets do not penetrate to the IB when an isolated boundary heat source is placed at low latitudes or when two heat sources are placed at antipodal locations. This may explain why such jets are not observed in our simulations, which have the strongest heat flow anomalies in the equatorial region and contain regions of anomalously high and anomalously low heat flow.
The strong reduction of persistent IB thermal anomalies with E in our models arises because of the disparity between the characteristic lengthscales of the flow, l, and the boundary forcing, L. It is well known that the influence of thermal boundary forcing on the global dynamics is most significant when L = l, in which case the whole flow pattern can become locked to the pattern of boundary heterogeneity (Zhang & Gubbins 1993). When L ≈ l quasi-locked solutions can be found at low Ra (e.g. Willis et al. 2007;Sreenivasan 2009), but as E (and therefore l) decreases all solutions become time-dependent (Davies et al. 2009). The heat that is transported near the boundary by the large-scale thermal wind flow (Sreenivasan 2009) is advected into the deeper interior by flow of much smaller scale and correspondingly shorter diffusion time, thereby reducing the probability that advective flow will retain its morphological identity all the way to the IB.
Another effect that stifles the transfer of the large-scale OB pattern to the IB is the emergence of a zonal flow at high Ra. This flow is aligned with the rotation axis and, in nonmagnetic convection, is strong and retrograde outside the inner core (Christensen 2002;Aubert 2005). The zonal flow shears convective anomalies and sweeps them around the inner core, smearing out any longitudinal variations that would arise in its absence.
It is worth considering effects not included in our simulations that might alter the lengthscale of the deep flow. Our simulations do not include a self-generated magnetic field, which reduces computational costs and allows systematic exploration of the low E, high Ra regime. A body of work now suggests that the appropriate dynamical balance at core conditions is between the pressure gradient and Coriolis effect at zeroth order (geostrophy) and between magnetic, Archimedian (buoyancy) and Coriolis terms at first order (MAC balance) (Davidson 2013;Aubert et al. 2017; Aurnou   (Dormy 2016;Yadav et al. 2016), which might hinder the propagation of OB effects to depth. Future dynamo simulations with heterogeneous boundary conditions are needed to quantify the relative importance of these effects. Magnetic fields also significantly alter the zonal flow. Aubert (2005) found prograde zonal flows in his simulations, which would still act to smear out the effect of OB anomalies as in our models. In non-magnetic convection the zonal flow is dominantly columnar and draws its power from the non-axisymmetric flow via the Reynolds stress; however, in dynamo simulations the magnetic field can break the axial alignment, in which case the zonal flow is powered by the thermal wind (Aubert 2005). Different scaling behaviour for the zonal flow speed is predicted in the two cases. In the simulations the zonal flow is weakened by the magnetic field (Aubert 2005;Schaeffer et al. 2017), and the scaling behaviour predicts that the dynamo-generated zonal flow would be slower at core conditions than the non-magnetic zonal flow [see eqs (3.1) and (3.2) of Aubert  Aubert (2005) predicts a zonal flow velocity of ∼10 −4 m s −1 in the core and so the time taken for one revolution of the zonal flow around the inner core is comparable to the time for fluid parcels to descend from the CMB to the ICB. Therefore it is possible that the smearing effect could be important in the core.
Figs 5 and 6 do not allow a simple interpretation of the amplitude of lateral thermal variations because δT θ and δT hom φ are poorly known for the core. To estimate the amplitude we derive empirical scalings based on T at the OB and IB (Fig. 7). At the OB T increases with increasing q and is approximately constant with E and Ra. At the IB T is much less sensitive to the OB heat flow pattern and q and decreases with Ra/Ra c at high Ra. These observations suggest fitting T to a function of the form (Ra/Ra c ) a 10 b . The dependence of T on the heat flow pattern is not simple, while the three values of E in our suite of simulations do not allow a robust estimate of the E-dependence of T . Therefore we empirically fit T for each q , E and heat flow pattern. Results are given in Table 1, which also provides values of T extrapolated to Ra/Ra c = 10 3 and 10 8 , which span the range of plausible estimates provided by Gubbins (2001). Note that the extrapolation to Ra/Ra c = 10 8 at higher E should be interpreted with caution since such high values of Ra/Ra c may not lie in the rapidly rotating regime thought appropriate for core dynamics.
As expected, the results in Table 1 show weak Ra dependence at both the OB and IB, while the E-dependence is weaker at the OB than at the IB. The extrapolated values of T can exceed 10 at the OB for the most extreme cases and are always <1 at the IB. To estimate max (T ) − min (T ) for the Earth we assume T = 10 −3 K, the lower end of estimates given in Section 1. In this case the lateral thermal anomalies in the core are estimated to be a few mK at the core-mantle boundary and around one tenth of a mK at the inner core boundary. However, at the OB, changing q from 2.3 to 5.0 increases the predicted value of T by an order of magnitude or more. The trend with increasing q is difficult to determine from our dataset. If the change in T scales like (q ) 2 then extrapolating to q = 35 would suggest T ∼ 10 2 − 10 3 and lateral temperature anomalies of ∼0.01 − 0.1 K; if T scales like q then T ∼ 10 3 − 10 4 at q = 35. These values are comparable with estimates from geodynamo simulations for the temperature anomalies arising from penetration of boundary heat flow anomalies into a stable layer below the OB (Christensen 2018). At the IB, at sufficiently low E, the amplitude of lateral variations is insensitive to the pattern or amplitude of OB heat flow. Therefore temperature anomalies at the ICB are expected to be comparable to those arising in homogeneous convection, nominally O(10 −4 ) K.

C O N C L U S I O N S
We have studied non-magnetic convection in a rotating spherical shell with lateral heat flow variations imposed at the OB. Our study covers the parameter range 10 −4 ≤ E ≤ 10 −6 , 1 ≤ Ra/Ra c ≤ 800 and considers Y 1 1 and tomographic heat flow patterns with amplitudes measured by the parameter q = 2.3, 5 as well as the homogeneous case with q = 0. We have focused on the time-averaged temperature anomalies produced at the IB and OB of the fluid shell. We find that: (i) At the OB, latitudinal and longitudinal temperature variations induced by the boundary forcing greatly exceed variations produced by homogeneous convection and are most prominent in the equatorial region. Longitudinal variations exceed latitudinal anomalies; based on empirical scaling estimates they are weakly dependent on E and Ra and may reach O(1) K if the forcing is strong. Such strong temperature variations may be seismically observable and searching for lateral variations in seismic velocity near the top of the core may provide new constraints on the amplitude of temperature variations induced by CMB heterogeneity.
(ii) At the IB, boundary-induced temperature variations become practically indistinguishable from temperature anomalies due to homogeneous convection as E decreases. In our most extreme models latitudinal variations exceed longitudinal variations. Induced anomalies decrease strongly with E and are weakly dependent on Ra and q . Extrapolated to core conditions our results suggest that heat flow variations on the CMB are unlikely to explain the largescale variations observed by seismology at the top of the inner core.

A C K N O W L E D G E M E N T S
CD acknowledges a Natural Environment Research Council personal fellowship, reference NE/L011328/1. This work used the ARCHER UK National Supercomputing Service (http://www.arch er.ac.uk) and ARC2, part of the High Performance Computing facilities at the University of Leeds, UK. We are grateful for the constructive comments provided by two reviewers.