Constraining mantle convection models with palaeomagnetic reversals record and numerical dynamos

S U M M A R Y We present numerical models of mantle dynamics forced by plate velocities history in the last 450 Ma. The lower-mantle rheology and the thickness of a dense basal layer are systematically varied and several initial procedures are considered for each case. For some cases, the dependence on the mantle convection vigour is also examined. The resulting evolution of the CMB heat flux is analysed in terms of criteria to promote or inhibit reversals inferred from numerical dynamos. Most models present a rather dynamic lower mantle with the emergence of two thermochemical piles towards present-day. Only a small minority of models present two stationary piles over the last 450 Myr. At present-day, the composition field obtained in our models is found to correlate better with tomography than the temperature field. In addition, the temperature field immediately at the CMB (and thus the heat flux pattern) slightly differs from the average temperature field over the 100-km thick mantle layer above it. The evolution of the mean CMB heat flux or of the amplitude of heterogeneity seldom presents the expected correlation with the evolution of the palaeomagnetic reversal frequency suggesting these effects cannot explain the observations. In contrast, our analysis favours ‘inertial control’ on the geodynamo associated with polar cooling and in some cases break of Taylor columns in the outer core as sources of increased reversal frequency. Overall, the most likely candidates among our mantle dynamics models involve a viscosity increase in the mantle equal or smaller than 30: models with a discontinuous viscosity increase at the transition zone tend to agree better at present-day with observations of seismic tomography, but models with a gradual viscosity increase agree better with some of the criteria proposed to affect reversal frequency.


I N T RO D U C T I O N
The present-day long wavelength pattern of the lower mantle has emerged consistently from global models of shear-wave velocity tomography with various parameterizations (Lekic et al. 2012).It involves two large low shear wave velocity provinces (LLSVPs) beneath the Pacific and Africa.The dynamical interpretation of these seismic structures is still under debate.Some authors favour a purely (or dominantly) thermal origin (isochemical models; Schuberth et al. 2009a,b;Davies et al. 2012) with LLSVPs corresponding to clusters of thermal plumes.Since tomography models also report a negative correlation between shear wave and bulk sound velocities (e.g.Masters et al. 2000) suggesting compositional heterogeneity in the deepest part of the lower mantle, a growing number of dynamical models interpreted LLSVPs as piles of intrinsically denser material (thermochemical models, cf.e.g.McNamara & Zhong 2005;Deschamps & Tackley 2008;Tackley 2011;Deschamps et al. 2012) either of primordial origin or resulting from the recycling of lithospheric material.
A second debate is whether thermochemical piles correspond to relatively stationary features in the last hundreds of million years or these structures evolve significantly.A spatial correlation between the margins of LLSVPs deep in the mantle and the locations at Earth's surface of present-day hotspots (Thorne et al. 2004) as well as the reconstructed eruption sites of Large Igneous Provinces (LIPs) and kimberlites in the last 100 s of Myr (Burke & Torsvik 2004;Burke et al. 2008;Torsvik et al. 2008Torsvik et al. , 2010;;Burke 2011) has been interpreted as fixed birth sites of hot convective instabilities and a global morphology of LLSVPs that remained unchanged during this period of time.The hypothesis of fixed LLSVPs is at odds with some modelling efforts for mantle convection with prescribed plate motion histories that report dynamic thermochemical piles (Zhang et al. 2010;Zhang & Zhong 2011), with the African superplume possibly forming as a consequence of Pangea assembly as in the alternating degree-1/degree-2 scenario initially proposed by Zhong et al. (2007b).Further numerical tests (Bower et al. 2013;Bull et al. 2014;Rudolph & Zhong 2014;Zhong & Rudolph 2015) did not lead yet to a consensus.For syntheses on this debate, we refer to Hernlund & McNamara (2015) and Zhong & Liu (2016).
Such a discussion on the evolution of mantle dynamics is emblematic.While one can expect that further numerical tests and improvements of numerical tools will contribute to clarify apparent conflicts, it is also pertinent to consider external constraints on such models.Concerning mantle evolution in the last 100 s of Myr, such constraints include the geological evidence for the assembly and dispersal of supercontinents (cf.e.g.Li & Zhong 2009), as exemplified by the reconstruction of LIPs/kimberlites locations (Torsvik et al. 2010).Other observations such as True Polar Wander (TPW; cf.Besse & Courtillot 2002) or global sea-level changes (cf.e.g.Conrad 2013) are also typically deciphered in relation to past mantle dynamics.However, the evolution of deep mantle structures is hardly resolved by the evidence of TPW.In addition, although global sea-level change is clearly sensitive to upper-mantle features, its relation to the number of LLSVPs and their locations through possibly associated dynamic topography is unclear (Conrad 2013).Another constraint could arise directly from plate motions.Conrad et al. (2013) analysed plate motion models in terms of the decomposition between degree-1 and degree-2 of the divergence of velocities since 250 Ma: These models were compatible with the presence of an upwelling beneath Africa during the whole period.Rudolph & Zhong (2013) argued nevertheless that this does not necessarily imply a fixed African pile.
Here, we take advantage of palaeomagnetic data, more specifically of long term variations of the frequency of reversals (Geomagnetic Polarity Time Scale, GPTS) in an attempt to constrain mantle dynamics scenarios.While the chaotic core convection has been proposed as an explanation for the very irregular variations in GPTS (Ryan & Sarson 2007;Wicht et al. 2009), mantle control is a much more probable mechanism due to the difference between GPTS timescales (e.g.superchron duration) and the much shorter core convection timescales as well as due to the similarity between GPTS timescales and mantle convection timescales (Glatzmaier et al. 1999;Kutzner & Christensen 2004;Driscoll & Olson 2011;Olson et al. 2013;Amit & Olson 2015).Recent attempts to interpret mantle evolution models as potential causes for GPTS favour a dynamical mantle (Zhang & Zhong 2011;Olson et al. 2013).Here, we propose to systematically evaluate the history of CMB heat flux inferred from such models in light of their potential to inhibit or promote reversals.
As in Zhang & Zhong (2011), our approach is to prescribe plate velocities as an upper kinematic boundary condition to mantle convection models (Section 2).Two main aspects are varied in the mantle model: the global rheology and the thickness of a dense basal layer (the isochemical hypothesis is not investigated).The resulting CMB heat flux is analysed according to specific criteria for reversal frequency inferred from numerical dynamos (Section 3).Results are described in Section 4 and discussed in Section 5. Perspectives are detailed in Section 6.
The philosophy of this paper is similar to that of Zhang & Zhong (2011): we assume that the relations between CMB heat flux aspects and reversal frequency inferred from the numerical dynamo literature are correct.Because hypothetically each one of the criteria may be much more influential than the rest, we do not assign (a generally unknown) relative weight to each criterion.We therefore consider a mantle model adequate if its time-dependent CMB heat flux complies with the evolution of reversal frequency according to at least one of the various dynamo criteria.However, we do not pretend to propose a best candidate mantle model.Instead, we examine the effect of each ingredient in the mantle model on the agreement between CMB heat flux and reversal frequency evolutions according to the dynamo criteria.

Principle
The link between plate velocities and mantle flow is inherent to the plate tectonics theory and already appears in the convection hypothesis proposed by Holmes (1931) to account for the then newly proposed idea of continental drift.Early attempts to couple instantaneous mantle flow to present-day plate velocities date from Hager & O'Connell (1979).Bunge et al. (1998) proposed the first numerical model where a plate motion history since 120 Ma is prescribed as a kinematic boundary condition at the surface of a mantle convection model.Since then, increasingly sophisticated models and plate motion histories of growing durations and improved precision have been proposed- Zhang & Christensen (1993) termed this group of numerical models 'semi-dynamic', a term we use in the following.For a clear synthesis, we refer to Zhong & Liu (2016) where such semi-dynamic models are described and their predictions are compared with various observations, first and foremost, global tomography models.
Upper-mantle dynamics are extremely simplified in most of the semi-dynamic models introduced so far, most importantly because the dichotomy between oceans and continents is not introduced and because the rheology is simplified.In the most realistic case (McNamara & Zhong 2005;Zhang et al. 2010;Zhang & Zhong 2011;Bower et al. 2013;Bull et al. 2014;Rudolph & Zhong 2014;Zhong & Liu 2016), as in this study, viscosity is temperatureand pressure-dependent.As a consequence, a cold and viscous lid should develop naturally beneath the surface (Solomatov 1995) if the 'Earth-like' plate regime is not forced by prescribed plate velocities.Note also that the temperature-dependence of viscosity is reduced for numerical reasons -possible consequences are discussed below.Possible ways to overcome this problem include the prescription of the shallow thermal structure of slabs as proposed by Bower et al. (2013) and Flament et al. (2014) or the use of self-consistent plate rheology (Bello et al. 2015).We consider in this study that the role of upper-mantle dynamics is only to provide reasonable locations and timing for the penetration of slabs in the lower mantle.It should be recalled however that a larger temperature-dependence as well as the use of pseudo-plasticity can affect the morphology of slabs even at large depths (Bello et al. 2015).

Plate velocity model
Defining plate motion in the past is anything but trivial, and the further back in time kinematic models go, the more the uncertainty increases.In order to reconstruct past plate motion seafloor age maps are the most robust data sets that can be used.As the current seafloor is not older than 180 Ma, reconstructions that go beyond this maximal age clearly rely on assumptions that are more speculative.The geological record provides indirect arguments to define a plausible temporal evolution of surface plate motion (Seton et al. 2012).Available palaeomagnetic data prior to 180 Ma constrain the palaeolatitudes of continental areas and allow to reconstruct plate motions at earlier times (Zhang et al. 2010).Strictly remaining immune from any considerations on the relationships with the underlying mantle structure is impossible because an absolute reference frame needs to be set, which is by essence based on conjectures on the dynamic interplay between the mantle and the lithosphere.The hotspot reference frame-either fixed (e.g.Morgan 1972) or moving (e.g.O'Neill et al. 2005)-certainly provides the most commonly used basis, but other observables such as seafloor spreading directions (e.g.Williams et al. 2016) or TPW (e.g.Steinberger & Torsvik 2008) can also be used.Seton et al. (2012), whose reconstructions we implement in the following, use a hybrid reference frame that combines a moving hotspot reference frame (O'Neill et al. 2005) back to 100 Ma and a TPW based reference frame (Steinberger & Torsvik 2008) for earlier times.Overall, it should be kept in mind that such inferences are based on premises that become increasingly debatable in the distant past so that conjecture gradually takes over certainty, to an extent that is hardly quantifiable.Alternatively, based on the hypothesis of fixed LLSVPs, a plate motion model that allows going back deep in Palaeozoic times has been proposed (Domeier & Torsvik 2014;Torsvik et al. 2014).Because our objective is to explore the relationships between plate motion, mantle convection, and core dynamics, hypotheses on the past structure of the mantle are discarded.Other assumptions may help to reconstruct further backwards in time, including minimizing the longitudinal drift (Torsvik et al. 2008), linking slab remnants in the lower mantle to geological structures at the surface (Van Der Meer et al. 2010), considering subduction kinematics (Williams et al. 2015) or relying on a reference frame based on hotspot tracks in both hemispheres (Conrad et al. 2013).
We assembled a comprehensive plate motion history within the past 450 Ma based on the compilations of Zhang et al. (2010) and Seton et al. (2012).23 plate motion snapshots are used between 450 Ma and 120 Ma corresponding to intervals of varying durations ranging between 5 and 30 Myr (Zhang et al. 2010).11 plate motion snapshots are used for the last 120 Ma, again corresponding to intervals with varying durations ranging between 6 and 18 Myr (Seton et al. 2012).The average plate velocities for these models are shown in Fig. 1.

Mantle convection model
Global scale mantle dynamics are described by solving the conservation equations for thermochemical convection in the infinite Prandtl number approximation.Since plate velocities are prescribed as a surface mechanical boundary condition in our models, external work is supplied to the system and the Boussinesq approximation is preferred to one that would involve dissipation.Boundary conditions are as follows: time-dependent horizontal velocities are prescribed at the surface (see Section 2.2) and free-slip at the CMB; dimensionless temperature is set to 0 at the surface and 1 at the CMB.Initial conditions are described in Section 2.3.3.Aside from these, the viscosity law and four dimensionless numbers (Ra, Bu, h, f) suffice to determine the various numerical models.The basal Rayleigh number is Figure 1.Time evolution of the surface average prescribed plate velocity model (thick dashed line) as well as the volume average velocity for 54 numerical mantle convection models (coloured lines, colour code defined in Fig. 2).The two light blue rectangles indicate the KRS and the CNS.
with thermal expansivity α, reference density ρ, gravity g, temperature difference T, mantle thickness d, thermal diffusivity κ and basal reference viscosity η b = η(r = r b ).In the Boussinesq approximation, α, ρ, g and κ are uniform.Note that in reality these parameters vary: for example, the thermal expansivity may be reduced by a factor of two to three in the lower mantle (Tackley 2012).
The buoyancy number is with intrinsic density increase associated with the dense material ρ C .The dimensionless (homogeneous) volumetric heating rate is with dimensional volumetric heating rate h * and thermal conductivity k. h includes a radiogenic component but also mimics secular cooling in the mantle in this 'steady-state' model (cf.e.g.Choblet & Sotin 2000).Finally the ratio of inner radius to outer radius is with the core radius r b and the outer radius of the Earth r s .The values of these dimensionless numbers are identical for all numerical models.These dimensionless numbers as well as specific values for some variables that are used in the following to present results with physical dimensions, are given in Table 1.They are comparable to the values used by Zhang & Zhong (2011) in their reference model.In this study, a given mantle model can be distinguished by only two main aspects: the viscosity law and the structure of the dense layer at the base of the mantle.These aspects are described in the next sections.

Viscosity
Viscosity depends on temperature and pressure (or radius): where E a is the dimensionless activation energy and T the dimensionless temperature.The function η r mimicking variations with radius r is inspired by the law proposed by Zhang & Zhong (2011).η r = 1 in the upper 100 km of the model is mimicking the viscous lithosphere, η r = η um in the upper mantle below the lithosphere, η r = η lm exp (V a r ) in the lower mantle.Note that an abrupt viscosity jump η um /η lm is introduced in some cases across the transition zone (located at a depth of r s − r tz = 670 km).In the lower mantle, viscosity increases exponentially with depth r = (r tz − r)/(r tz − r b ).Finally, a constant η * = exp(E a − V a )/η lm ensures that the bottom viscosity reference value is set to 1.In our models, the lithosphere is intrinsically more viscous than the upper mantle 100 km below the surface (η um ≤ 1, see Table 2).In the real mantle, causes for such intrinsic variations might involve the dehydration of the oceanic lithosphere (Hirth & Kohlstedt 1996) and the presence of partial melt and strain weakening beneath the lithosphere.Such effects are shown to favour plate-like behaviour in self-consistent models of plate generation (Tackley 2000).These are simplified in the present models-note, for example, that the distinction between oceans and continents is not included.Practically, this choice for η r also enables our models to handle large viscosity gradients associated with surface plates as the effect of temperature is moderate.
In the absence of strong laboratory constraints on the viscosity of lower-mantle materials, geophysical observations (geoid, postglacial rebound) favour a viscosity increase factor of ∼10-100 in the lower mantle (cf.e.g.Hager et al. 1985;Ricard et al. 1993;Mitrovica & Forte 1997).It is not clear nevertheless whether such an increase occurs abruptly at the transition zone (due to changes in mineralogy and possibly composition) or more gradually in the Figure 2. Radial profiles for laterally averaged viscosity at the end of numerical models (corresponding to present-day plate velocities) for the six viscosity laws: V1 (red), V2 (blue), V3 (green), V4 (cyan), V5 (yellow) and V6 (purple).The same colours are used throughout the paper to refer to these various viscosity laws.The reference viscosity at the CMB is 5 × 10 20 Pa s. lower mantle (due for example to the effect of pressure for a given composition/mineralogy, as measured by the activation volume), even though efforts at modelling the sinking speed of subducting slabs tend to favour an almost isoviscous lower mantle unless postperovskite strongly weakens the lowermost mantle ( Čížková et al. 2012).Here, we investigate both endmembers (cf.Table 2): viscosity models V1, V2 and V3 prescribe a sharp increase η um /η lm of 100, 30, 10, respectively, while the pressure effect is null in the lower mantle.In contrast, viscosity models V4, V5 and V6 prescribe nonzero activation volumes V a leading to overall increases of 100, 30, 10 in the lower mantle, respectively, while there is no discontinuity in viscosity profiles across the transition zone.
Fig. 2 displays the radially averaged viscosity profiles for all numerical models.A given viscosity model corresponding to a given function η r is associated with a specific colour.Differences among profiles of the same colour thus only result from different temperature fields among the models.For numerical reasons, models V4 and V5 are associated with a larger value of η um than V6.Hence, the evolution from V4 to V5 and to V6 is not as straightforward as the one from V1 to V2 and to V3.However, this will be shown to have only moderate consequences in the following.

Dense basal layer
As indicated above, compositional gradients are likely present in the deepest mantle although the detailed structure and its dynamical effect is still a matter of debate.Here we consider a rudimentary model for thermochemical piles as a dynamical equivalent to LLSVPs.In the context of our study, such piles are a major ingredient of large scale mantle dynamics whose interplay with subducting slabs shapes heat flux heterogeneities at the CMB.Whether these piles are of primordial origin or caused by continuous recycling of oceanic lithosphere (Deschamps et al. 2012) is only of minor Unless stated otherwise (see Section 2.3.3), a dense basal layer is initially prescribed above the CMB with a flat interface at radius r co (initial thickness d co = r co − r b ).The amount of material included in LLSVPs might correspond to ∼2 per cent of the mantle volume (Hernlund & Tackley 2008).We thus consider three values for r co (see Table 3) corresponding to 64 km (i.e. a volume V co 1 per cent), 127 km (V co 2 per cent) and 191 km (V co 3.5 per cent), which are denoted A, B, and C, respectively.For simplicity, the dense material is associated with a single value of the buoyancy number, Bu = 0.5.We performed several numerical tests indicating that, to a first order, the long term dynamics of the dense basal layer is dictated by its global buoyancy measured by the product (Bu × d co ) so that, for example, cases with a 127 km-thick layer and a twice smaller value of Bu behave similarly to case A. We do not introduce any specific viscosity contrast associated with the dense material, nor extra heat sources possibly caused by an enrichment in heat-producing elements.Tests indicated that the latter effect is of secondary importance.Though similar tests indicated that the viscosity effect is significant, it is not investigated here, partly due to the absence of physical constraints.

Numerical set-up
We use the OEDIPUS numerical tool (Choblet 2005;Choblet et al. 2007).The compositional field is advected using a high-resolution scheme based on the Superbee flux limiter method also used for conservation of energy.The Lenardic filter treatment (Lenardic & Kaula 1993) is used to further reduce the effect of numerical diffusion which results in good agreement with the 3D tests presented by Tackley & King (2003).We acknowledge that the use of tracers would ultimately be the most appropriate method.Nevertheless, we consider that our approach is sufficient as only the evolution of the shape of the dense basal layer is of interest in our study (i.e.subtle chemical exchanges that might be biased with a field method are not considered here).In all models, the grid mesh involves 128 cells in the radial direction and 6 blocks of 64 × 64 cells in lateral directions.
The models discussed above for the viscosity law and the dense basal layer define 18 groups of models denoted by a code combining the two parameters, for example, V6-B denotes viscosity law V6 and dense basal layer B. In order to test the sensitivity to ini-tial conditions, three numerical models are performed for each of the 18 groups: two of these, noted with indices i and ii correspond to 900 Myr long models, while the third one, noted with index s, correspond to a shorter, 450 Myr long model.As indicated above we use plates velocity models for the last 450 Myr.In the case of longer models (i and ii), we initially run a 450 Myr long preliminary model with stationary plate velocities corresponding to the oldest plate velocity model (t = 450 Ma).Plate velocities are then varied in the remaining time of the models (t ∈ [450 Ma,present]).For the short model (s), only the second phase is performed.Models i and s are initialized from a flat interface for the dense basal layer and a radial thermal structure is prescribed corresponding to the average steady-state solution of a similar model without the dense basal layer.Models ii start from the final state of models i (t = 0 Ma).

T H E R E L AT I O N B E T W E E N C M B H E AT
Numerical dynamo simulations with imposed outer boundary heat flux shed light on the relation between various aspects of the CMB heat flux and the reversal frequency (for a review see Amit et al. 2010).It was found that the reversal frequency may depend on the mean heat flux, the amplitude of its lateral heterogeneity, and the pattern of the heterogeneity.Some of these criteria (but not all) were already used by Zhang & Zhong (2011) to test the consistency of their mantle convection models with the palaeomagnetic reversal frequency.Below we review the results of numerical dynamo studies for an updated and more complete set of criteria on CMB heat flux concerning reversal frequency.

Mean CMB heat flux
Numerical dynamos with homogeneous outer boundary conditions (isothermal or fixed heat flux) generally find that for a given rotation rate and core fluid physical properties there is a critical threshold of convection vigour for the onset of reversals (Christensen & Aubert 2006;Olson & Christensen 2006;Aubert et al. 2009;Heimpel & Evans 2013).When all other parameters are fixed, stronger core convection driven by larger CMB heat flux results in stronger inertial effects (Christensen & Aubert 2006;Olson & Christensen 2006;Aubert et al. 2009) and reversing dynamos, whereas more moderate core convection associated with weaker CMB heat flux may result in non-reversing dynamos.Furthermore, studies of reversing numerical dynamos found that increasing CMB heat flux yields systematic increase in reversal frequency (Kutzner & Christensen 2004;Driscoll & Olson 2009;Olson & Amit 2014).Stronger core convection increases turbulence and mixing, rendering the dipole more vulnerable to collapses and eventual reversals.This criterion is generally considered as the most robust condition for reversal frequency.
Thermal convection in the core occurs when the heat flux exceeds the adiabat.High estimates of core thermal conductivity (Pozzo et al. 2012;Gomi et al. 2013) suggest that parts of the outer core are stably stratified, whereas lower estimates (Konôpková et al. 2016) favour whole core convection.In any case increased CMB heat flux (even within a subadiabatic regime) would enhance inner core growth rate and strengthen compositional core convection, hence increase reversal frequency.Monitoring the mean CMB heat flux is therefore useful even if convection in the core is not dominantly thermal.1170 G. Choblet, H. Amit and L. Husson We therefore compare time-series of the mean CMB heat flux q 0 with the palaeomagnetic reversal frequency.A model is considered successful if the correlation coefficient between these two properties is large and positive.

Amplitude of CMB heat flux heterogeneity
Numerical dynamos with various patterns of heterogeneous outer boundary heat flux found that in general increasing the amplitude of the heterogeneity (for a given pattern) increases reversal frequency (Olson et al. 2010;Heimpel & Evans 2013;Olson & Amit 2014).Larger CMB heat flux heterogeneity produces locally stronger heat flux and hence locally stronger core convection and more turbulent conditions.It is therefore possible that reversals may be triggered locally in locations where the CMB heat flux is anomalously large.
For each model we monitor the time evolution of the amplitude of the heat flux heterogeneity q * = (q max − q min )/2q 0 (e.g.Olson & Christensen 2002).A model is considered successful if the correlation coefficient between q * and the reversal frequency is large and positive.

Equatorial versus polar cooling
Several studies of reversal frequency in numerical dynamos found that large equatorial heat flux, most commonly represented by a positive spherical harmonic Y 0 2 (and in some cases Y 0 4 as well), increases reversal frequency, whereas large polar heat flux (i.e.negative Y 0 2 ) stabilizes the dipole and may even induce superchrons (Glatzmaier et al. 1999;Kutzner & Christensen 2004;Olson et al. 2010).Polar cooling may localize downwelling at high-latitudes and field concentration there may maintain the axial dipole; Conversely, equatorial cooling may attract magnetic flux and hence the dipole axis to the equator, thus inducing reversals (Amit et al. 2010).Some exception to this rule was found by Olson et al. (2010).Although their Y 0 2 heat flux pattern dynamo models produce by far the highest reversal frequency, somewhat surprisingly their models with a −Y 0 2 pattern reversed more frequently than their uniform reference case.Olson & Amit (2014) found that right beyond the onset of reversals equatorial cooling indeed enhances reversal frequency, in agreement with previous studies (Glatzmaier et al. 1999;Kutzner & Christensen 2004), a state they termed 'geographic control'; however, well within the reversing regime, polar cooling enhances reversal frequency, a state they termed 'inertial control'.In the latter scenario, polar cooling drives a meridional circulation that is compatible with that of the homogeneous dynamo (Aubert 2005;Amit & Olson 2006), hence enhancing core convection, turbulence and reversal frequency (Olson & Amit 2014).
Some studies reported dynamo failure when the CMB heat flux heterogeneity is increased (Olson & Christensen 2002;Sreenivasan 2009).However, in these studies the heterogeneity amplitude is large-in some places heat flows from the mantle to the core, potentially violating the Boussinesq approximation on which the dynamo models rely (e.g.Dietrich et al. 2016).Even with such large heterogeneity amplitudes not all dynamos fail (Stanley et al. 2008;Dietrich & Wicht 2013).
For each model we calculate the time evolution of q 0 2 .A model is considered successful if the correlation coefficient between q 0 2 and the reversal frequency is large whether the value is positive or negative, corresponding to the 'geographic control' and 'inertial control' scenarios, respectively.

Equatorial symmetry
Rapid rotation effects govern the fluid dynamics in Earth's atmosphere, ocean and outer core.These rotational effects are expected to produce a flow that is invariant in the direction parallel to the rotation axis and thus symmetric about the equator (e.g.Busse 1975;Jault 2008).Such equatorially symmetric columnar flow may concentrate opposite polarity magnetic field structures at equatorially symmetric locations on the opposite hemispheres, and by that maintain the dominant axial dipole that characterizes the magnetic fields of Earth, Jupiter, Saturn and possibly Mercury.Breaking this equatorial symmetry by convection driven by the lower-mantle heterogeneity may therefore break these fluid columns, diminish the axial dipole and result in reversals (Berhanu et al. 2007;Pétrélis et al. 2009Pétrélis et al. , 2011;;Biggin et al. 2012).
For each model we compute the time evolution of the level of equatorial symmetry given by the power in the symmetric heat flux normalized by the total power.Note that here a model is considered successful if the correlation coefficient between the level of equatorial symmetry and the reversal frequency is large and negative.

R E S U LT S 4.1 A typical numerical model
We select model V2-Bi (sharp viscosity increase factor of 30 at the transition zone and initially flat dense basal layer of intermediate thickness) as an example to describe the outcome of the numerical models.Fig. 3 displays radial profiles of quantities characterizing heat transfer and flow at the end of the model.The average temperature profile (black, left panel) has the typical structure obtained for a convective layer.The temperature value in the bulk of the layer ( 0.6) essentially results here from the interplay between sphericity that tends to lower it when compared to the Cartesian case, and volumetric heating that results in an increase in temperature.A simple scaling relationship for this set up (ratio of inner to outer radii f = 0.55 and basal power corresponding approximately to one third of the power through the surface) would lead to a 10 per cent smaller value for the internal temperature in the case of an isoviscous fluid (Choblet 2012), showing that the complex viscosity profile does not affect dramatically the average temperature.This is due to the non-monotonic radial dependency of the viscosity (function η r ) as a purely temperature-dependent viscosity would cause a much larger departure from the isoviscous case (Solomatov 1995).
Maximum and minimum temperatures at a given radius (red and blue curves, left panel, Fig. 3) reflect the instantaneous dynamics.Small bumps denote the radii reached by ascending or descending instabilities in this specific snapshot.Longer-term features such as the shape of the maximum temperature profile between r − r b = 0 and r − r b = 0.5 denote the influence of the dense basal region (compare with the green curve in the left panel, Fig. 3).A similar signature is observed in the viscosity profiles (red and black curves, centre panel, Fig. 3).Internal convection in this compositional layer can even be detected in the radial velocity profile as a local maximum in ascending velocities (red curve, right panel, Fig. 3) at r − r b 0.15.
Other notable features are the abrupt viscosity changes involved in function η r (eq.5), that is, independent of the actual temperature profile, especially at the base of the lithosphere (centre panel, Fig. 3): due to simplifications in our rheological model, cold subducting slabs become much less viscous once they reach the bottom of the lithosphere arbitrarily fixed at r − r b 0.965, which is a non-realistic feature.Beneath the lithosphere, the upper mantle is a low viscosity region.As a consequence strong local maxima in velocities are always located there (blue and red curves, right panel, Fig. 3).Finally, it must be emphasized that, with the possible exception of the average composition at a given depth (green curve, left panel, Fig. 3), in the case of an initially flat interface (models i and s), the global shape of these profiles remain similar throughout the numerical models.Fig. 4 displays the temperature and composition fields at the end of model V2-B i for radii in the lower part of the mantle (corresponding to discrete indices 1, 2, 4, 8 and 16, see the associated values of r − r b in Fig. 4).Focusing on the compositional field, a structure with thermochemical piles is clearly identified, displaying a decent agreement with the observed locations of LLSVPs.This example is actually associated with one of the best correlations with tomography: a linear correlation coefficient between the global tomography model of Masters et al. (2000) and the composition field averaged over the lower 180 km, both truncated at max = 8, provides a significant value of 0.73 and will be discussed in detail below (cf.Section 4.3).Nevertheless, although a global degree 2 is observed, the Pacific pile is split in two parts: while the possibility of such a bimodal structure for the Pacific pile is suggested by global tomography (cf.e.g.Lekic et al. 2012), it is much more pronounced in this model.The northeastern extent of the Pacific pile below Japan in this model also differs from tomographic observations.In this example where the amount of dense material is moderate (case B where the initial condition is a flat upper interface 127 km above the CMB), the material corresponding to the African pile does not reach elevations above 500 km (Fig. 4d).Material from the Pacific pile culminates approximately 900 km above the CMB (Fig. 4d).Extremely small fractions of dense material can still be observed up to the middle of the layer, but such features probably correspond to the resolution limit for the entrainment of dense material by our numerical approach.
The global features of the temperature field (Fig. 4, left) reflect the same structure although less clearly as both the dense basal layer and the mantle above convect efficiently so that significant lateral temperature gradients can be observed.In this example the coldest region in the dense basal material approximately has the same temperature as the hottest region of the normal mantle, especially in the lowermost mantle.Compare for example, the two light red regions adjacent to the North-South oriented, dark red, linear feature in the Southern Atlantic (Fig. 4a): this ridge corresponding to the hottest material is associated with the boundary of the African pile, with normal mantle material on the West and denser material on the East.As a consequence, the agreement in global pattern with the tomography is not remarkable (Pearson coefficient: 0.44, cf.Section 4.3).The margins of the dense piles often correspond to these hottest regions but some of the hottest material also forms linear features within the piles.These sheet-like hot upwellings develop in the upper part of the lower mantle, above the upper interface of the bulk dense basal regions (Fig. 4d, above the African pile).In the upper mantle (not shown here), the hot upwellings form plumes rather than linear features.These plumes are not as distinct since the average temperature is closer to the maximum temperature.The temperature pattern there mostly reflects the location of downwelling cold slabs and this extends in the lower mantle (the transition zone effect on viscosity is included in these models, but not on buoyancy, Fig. 4e).The distribution of the coldest regions in the lowermost part of the mantle (Figs 4a-c) reflects the interplay between earlier plate motions and the inherent dynamics of dense piles and normal mantle.
The time evolution of the criteria we propose for the CMB heat flux in terms of promoting higher or lower reversal frequency is exemplified in Fig. 5 (red curves).As these time-series depend stronger on the initial configuration of the numerical model (i.e.either group i, ii or s, for a given rheology and a given structure of the dense basal layer) at earlier stages than at later stages, and because the palaeomagnetic record is most reliable from the ocean floor, we decided to only focus on the last 160 Ma, a period that includes the Cretaceous Normal Superchron (CNS) but not the Kiaman Reversed Superchron (KRS).We systematically compared the evolutions of these criteria with the observed variations in reversal frequency [a 5 Myr running average window has been performed on the GPTS polarity record (Gradstein et al. 2012), see blue curve in Fig. 5].In the example of numerical model V2-B i , the correlations between the two curves for each panel are respectively: −0.06 (average q 0 ), −0.27 (amplitude of heterogeneity q * ), −0.56 (degree 2, order 0, q 20 ), 0.16 (equatorial symmetry).In this case, a successful correlation therefore appears only for the polar cooling criterion: the CMB heat flux is larger at high latitudes throughout (q 0 2 < 0) and minimal polar cooling during CNS (especially when compared to later times) is in agreement with 'inertial control' (Olson & Amit 2014).Other criteria behave in contrast to inferences from numerical dynamos.The q 0 criterion, expected to reach a minimum during the CNS, increases with a recent peak.The q * criterion, also expected to be minimal during the CNS, peaks slightly before the CNS and dips to present-day hyper-reversing state.Equatorial symmetry, expected to peak during the CNS, increases steadily during the whole period.

Mantle dynamics in the last 450 Myr
Overall, 54 numerical models are performed.Fig. 1 displays the time evolution of the volume average velocity for all these numerical models, corresponding to the active phase of long models (groups i and ii) (i.e. when plate velocities evolve with time) and to the whole duration of short models (group s, plate velocities evolve throughout the entire simulation).The prescribed plate velocities (thick dashed curve) affect internal dynamics.Peak surface velocities during CNS or before the KRS around 350 Ma result in peak volumetric velocities at the same periods.However, while changes in surface velocities are clearly visible on the whole mantle velocity field, the latter evolves on longer timescales, and some high frequency temporal fluctuations in the plates do not appear as clearly in all the mantle models (see e.g.peak at 240 Ma).Dif- ferences in averaged volumetric velocity amplitudes among the models of a factor smaller than 2 can be attributed to the rheological structure of the mantle.Globally, these differences constitute a monotonic relationship between average viscosity at the end of the model in the lower mantle and the average amplitude of the rms-velocity (compare Figs 1 and 2).Models of group s display initially (i.e. at t 450 Ma) a larger slope corresponding to the initiation of mantle dynamics (Fig. 1).Models i and ii also present such evolutions at t 900 Ma.In the case of group ii it is a large scale reorganization of mantle dynamics rather than its initiation; this results from the hiatus introduced by the restart procedure as the state corresponding to the present-day plate velocities (t = 0 Ma from model i) is injected as an initial condition while the initial plate velocities (t = 450 Ma) are prescribed.
Such a strong effect of surface plate velocities is observed in the evolution of the surface heat flux but not on the CMB heat flux (Fig. 6): the amplitude of temporal variability at the CMB is much smaller than at the surface, with some surface peaks extremely attenuated at the CMB.While the range of surface heat flux variations among the models is larger at 450 Ma (typically 15 TW) than at the end of models (typically 7 TW), reflecting the different initial conditions, it remains almost constant in time for the CMB heat flux (less than 10 TW, leading to large relative differences of a factor four for this smaller value).The absolute surface heat flux values are in agreement with expected heat power out of the convecting mantle (Jaupart et al. 2015) although this shall by no means be considered as an indication that the global mantle dynamics in these models reflect the real mantle dynamics.Values for the CMB heat power ranging from 3.5 to 12 TW at the end of models are at the lower end of what might be the real heat flux out of the core at present (cf.e.g.Lay et al. 2008;Labrosse 2014;Hernlund & McNamara 2015).In a significant fraction of the models, the CMB heat flux might be smaller than the flux along the adiabatic gradient in the core which would prevent thermal convection but nevertheless might not cancel the mantle control on the dynamo (see the discussion in Section 3.1).As will be emphasized below (see also Section 3), we focus here on the time evolution of either the average heat flux or of the pattern of lateral heterogeneities rather than their absolute values which are roughly consistent with the expected order of magnitude.
Figure 7. Time evolution of the long-wavelength structure of thermochemical piles for all numerical models.The compositional field is averaged for the 32 discrete layers above the CMB (bottom quarter of the mantle).The ratio of degree-2, order-2 power (S 22 ) to degree-1 power (S 1 ) is reported.Models are grouped according to the initial conditions.Top: models • s , middle: models • ii , bottom: models • i .The two light blue rectangles indicate the CNS and the KRS.The black half-disks on the right side denote the corresponding ratios in the seismic shear wave velocity model of et al. (2000) at the lowermost mantle.
We wish to stress that the long-wavelength structure of the lower mantle is far from stationary in most of our numerical models.Fig. 7 reports the time evolution of the long-wavelength compositional field averaged in the deepest region of our models (extending from the CMB to one fourth of the mantle radius).We consider the ratio of the degree-2 order-2 power (S 22 ) to the degree-1 power (S 1 ) as an indication of the presence or absence of two thermochemical piles centred at low latitudes.This measurement provides a test for the scenario proposed by Zhong et al. (2007a) in which the degree-1 structure prevails up to the formation of Pangea (up to 330 Ma) and later the present-day degree-2 structure emerges as a consequence of circum-Pangea subduction.The time evolution of this ratio is found to highly depend on the initial condition considered among various models.
For models • i (cf.Fig. 7, bottom), where the dense layer with an initially flat interface is first subjected to a 450 Myr long phase with velocities at 450 Ma prescribed at the surface, a global increase is obvious with the emergence of the present-day two piles structure between 200 Ma and 120 Ma, depending on the viscosity model.Larger viscosity contrasts between the lower and upper mantles lead to a slower formation of the two piles.Note that the present-day ratio is often smaller than the one observed by seismic velocity models (black half-disks in Fig. 7).As will be discussed later on, some of the models actually differ significantly at present from the structure in global tomography models [viscosity models V1 (red) and V4 (cyan), see Section 4.3].But the two-piles structure is visible in the lower mantle, with maybe the exception of model V6-B i leading to the smallest value of S 22 /S 1 at present.Common trends in all models • i denote the influence of surface velocities.Such a coherent envelope is not observed for models • ii (cf.Fig. 7, middle) where the fields obtained at present for models • i are prescribed as an initial condition and then subjected to a similar 450 Myr long phase with stationary surface velocities.While the overall two piles structure is injected in the models (with again the exception of model V6-B i ), these remain globally stationary only in the case of viscosity model V3 (with a modest viscosity jump at the transition zone) and in the specific case of model V2-A ii .In all other cases the piles either drift significantly while preserving a degree-2 order-2 pattern (viscosity models V2 and V1) or give rise to a degree 1 structure before 450 Ma which then keeps this spectral signature, although the pattern varies geographically (viscosity models V4, V5 and V6 with a gradual increase in the lower mantle).In the case of models • s (cf.Fig. 7, top), an initially flat interface is prescribed for the dense basal layer at 450 Ma, instantaneously subjected to time varying surface velocities.An initial adjustment stage is witnessed in these models that sometimes lasts up to 150 Myr, until the KRS.The subsequent evolution points to a clear amplification of the degree-2 order-2 structure until present-day.While the evolution in models • s (Fig. 7, top) apparently leads to long-wavelength present-day structures comparable to the one witnessed by global tomography, the global correlation will be shown to be less satisfactory than for models • i (cf.Section 4.3).
In conclusion, most models suggest that the deep mantle structure is quite dynamic in the last 450 Myr.This is essentially in agreement with the recent appraisal of the time evolution of the mantle at long wavelengths by Rudolph & Zhong (2014) or Zhong & Rudolph (2015) proposing that the present-day structure was formed around 200 Ma.In contrast to the results obtained by Bull et al. (2014), the dominantly degree-2 structure involving two piles still observed at present remains globally stationary only in four numerical models out of 17 where it is injected as an initial condition (Fig. 7, middle).It is interesting to note that the range obtained for the present-day structure in our models includes the values in global tomography models (Masters et al. 2000, half-disks in Fig. 7b).However, significant variations are observed among models, possibly because of the delay induced by the viscosity law.Some models should be rejected on this basis: this is discussed in the next section.

Correlation of lower-mantle present-day fields with global tomography
The individual linear correlation coefficients r between numerical results corresponding to the present-day and the global tomography model of Masters et al. (2000) are compiled in Fig. 8.The temperature (Fig. 8a) and the composition (Fig. 8b) fields are  (Masters et al. 2000) and either (a) the temperature field or (b) the composition field at the base of the mantle at present-day.Spherical harmonics decompositions of the three quantities have been truncated at degree max = 8.Colours indicate the viscosity model (as in Fig. 2): V1 (red), V2 (blue), V3 (green), V4 (cyan), V5 (yellow) and V6 (purple) (see Table 2).Suffix A, B or C (on the x-axis) indicates the model for the dense basal layer (cf.Table 3).Filled circles denote models of group i, empty circles denote models of group ii and empty stars denote models of group s.The five solid horizontal lines denote confidence limits above which there is (from bottom to top) a 50 per cent, 40 per cent, 30 per cent, 20 per cent or 10 per cent chance that the correlation occurs by chance.considered individually.These fields were averaged over the 8 discrete layers above the CMB (lowest 180 km) because the lowermost layer sometimes differ from the discrete layers above it, especially for the temperature field (see e.g.Fig. 4).These are compared to the lowermost layer of the tomography model, roughly corresponding to a similar thickness (Masters et al. 2000).For comparison of large-scale structures, both the numerical results and the tomography model are expanded in spherical harmonics until degree max = 8.The appraisal of the significance of the correlation coefficient r is then performed assuming that the statistic t = r (N − 2)/(1 − r 2 ) follows a Student's t-distribution (Press et al. 1992) in the case of no correlation, with N the number of independent coefficients (N = 45 for max = 8).The six solid horizontal lines in Fig. 8 then denote values above which the correlation is not fortuitous with confidence levels of 50 per cent, 60 per cent, 70 per cent, 80 per cent and 90 per cent.
A first observation is that the composition field correlates better than the temperature field with the pattern of seismic wave anomalies: for cases above the confidence limit of 90 per cent, the value of the correlation coefficient for composition always exceeds the value for temperature.This is caused by the clearer signature of thermochemical piles as compositionally distinct material than as hotter material (see Fig. 4, cf.Section 4.1).Second, as noted above, the present-day results strongly depend on the initial configuration and the duration of the numerical models.The values within a column (cases i, ii, or s) for a given viscosity structure (colour) and a given structure of the dense basal layer (suffix A, B or C) often differ significantly, with models i typically presenting a larger correlation coefficient.Third, the viscosity model (V j ) j = 1, 6 clearly affects the agreement with tomography: models associated with a 100 increase in the lower mantle (either discontinuous, V1, or gradual, V4) globally correspond to a poorer correlation.The increase of the amount of dense basal material (from A to B to C) is not as significant.
The statistically significant correlations in Fig. 8(b) indicate that the prescribed plate velocities tend to favour a shape for dense piles roughly consistent with what is observed at present by global tomography.Reasons for a clear departure from the tomographic pattern are multiple as can be illuminated by considering variables describing the interaction of cold slabs with the thermochemical piles (Fig. 9).The maximum sinking velocity of slabs (V slab ) ranges ∼1-2 cm yr −1 .Accounting for the variability of the sinking velocity with depth in the lower mantle of a given slab (Fig. 3) as well as possible variations from one slab to another, the mean sinking velocity (which is technically more difficult to obtain) may range ∼0.5-1.6 cm yr −1 .This range overlaps (though slightly underestimates) inferences from matching predicted slabs with tomographic features (Van Der Meer et al. 2010;Čížková et al. 2012;Steinberger et al. 2012).Calculations with a shorter duration s naturally involve a shorter interaction of cold slabs with the dense basal layer.As a consequence, the formation of piles is less developed which induces large values of variable C b corresponding to the surface area of the CMB covered by dense material at present.This also explains a less pronounced viscosity contrast in the lowermost mantle (ratio between maximum and minimum viscosity at a given depth, averaged over 180 km, denoted η in Fig. 9): while η provides a measurement of the typical viscosity contrast between cold sinking slabs and the hot dense basal material, its value is naturally decreased when slabs do not penetrate significantly the dense layer Figure 9. Three variables measuring the interaction between cold slabs and thermochemical piles: maximum viscosity contrast between slabs and piles η, maximum sinking velocity of slabs in the lower mantle V slab , and average value of the compositional field at the CMB C b .A dimensionless value of 10 3 for V slab corresponds approximately to 1 cm yr −1 .Same notation as for Fig. 8. down to the CMB.As expected, C b also increases with the volume of dense material (models A, B and C, see Table 3) although the effect on the correlation coefficient is minor.A value of C b = 0.21 was reported by Burke et al. (2008).If that is the case, all the s models should be discarded as they correspond to values larger than 0.5.Furthermore, among the remaining models, those involving a gradual viscosity increase (V4, V5, V6) are also problematic as these either induce too large values of C b (Fig. 9) or a small correlation coefficient with tomography (Fig. 8)-in the following, we will only consider the latter criterion.Finally, while two piles are observed at present in models V2-A ii and V2-C ii (as witnessed by their S 22 /S 1 ratios, cf.Fig. 7, middle), their locations do not correspond to the seismic observations which causes poor correlations.

Correlation of the time evolution of CMB heat flux criteria with observed reversal frequency
While the comparison of our models with tomography involved a region of 180 km above the CMB (cf.Section 4.3), we focus now on the heat flux precisely at the CMB (i.e. the locus of mantle control on the dynamo).In addition, while the composition field was earlier found to present a better agreement with tomography, heat flux variations in our models rely only on temperature anomalies (e.g.possible compositional effects on thermal conductivity are not taken into account in the framework of the Boussinesq approximation).We emphasize that, as a consequence, the present-day pattern for CMB heat flux shall not be confused with the pattern revealed by seismic tomography.
Fig. 10(a) compiles for all numerical models the linear correlation coefficient between the time evolution of each quantity indicative of the criteria listed above for CMB heat flux and the reversal frequency.Fig. 10(b) only displays models that provide present-day agreement with global tomography (i.e. a 90 per cent confidence that the composition field correlates with tomography, see Fig. 8): 33 out of the 54 numerical models conducted in the present study satisfy this criterion.In addition to model V2-B i shown in Fig. 5 and guided by the results in Fig. 10, we present in Fig. 11 four other examples of the time evolution of the various criteria.Corresponding temperature fields above the CMB are displayed in Fig. 12.
These results reveal first how, practically, a given criterion can be used to discriminate among mantle models.A first quality is that general trends can be observed in the results, that is, that the effect of the viscosity model or of the structure of the dense layer induce observable tendencies.Overall, we emphasize that while one single model presenting a good correlation for only one of the criteria could correspond to the real mantle and as such, may be satisfactory, we are here interested first in the lessons learnt from the whole ensemble of results.Which model or group of models might appear as the best candidates to explain reversal frequency variability is merely a secondary conclusion.
The criterion on the average CMB heat flux (q 0 in Figs 10 and 11) does not lead to the discrimination of successful mantle models: Figure 10.Linear correlation coefficient between the time evolution of reversal frequency and the time evolution of criteria for the CMB heat flux: average power q 0 , amplitude of heterogeneity q * , relative amplitude of degree 2, order 0 coefficient q 20 , equatorial symmetry.Same notation as for Fig. 8.Only the last 160 Myr have been considered, see Fig. 5. Shaded rectangles denote statistically significant correlation/anticorrelation (a 95 per cent confidence for the number of degrees of freedom considered, approx.20).(a) All models, (b) only models whose present-day composition field correlates sufficiently well with tomography (with at least a 90 per cent confidence, see Fig. 8).
only one model presents a significantly positive correlation (Fig. 10) as is expected if the influence of the average power were preponderant.About two thirds of the models actually lead to a negative correlation (as already noted by Olson et al. 2013).The reason is that fast plate motions during the CNS drive vigorous mantle convection (Fig. 1) which is opposite to the expectation for a superchron.Fig. 11(a) indeed shows a peak of q 0 at the end of the CNS, leading to anticorrelation for this criterion (see also Fig. 12).Furthermore, in all models, the time variability of global CMB heat flux is low (cf.Fig. 6).Because reversal frequency varies as the squared root of the mean CMB heat flux (Olson & Amit 2014;Amit & Olson 2015), this casts doubts on the likelihood that such weak CMB thermal evolution may cause the strong variations of observed reversal frequency during this period.As for the average CMB heat flux, only one of the mantle models (V4-B s ) satisfies the criterion measuring the amplitude of CMB heat flux heterogeneity (q * ) (Fig. 10a).In addition, this specific model is associated with a poor correlation at present-day with tomography, thus it is discarded (Fig. 10b).We note that for this criterion, the influence of initial conditions on the correlation coefficients is modest for models with a sharp viscosity increase at the transition zone (V1, V2, V3) as well as a moderate effect of the initial thickness of the basal layer.A clear result is that 43 out of the 54 models are associated with an anti-correlated evolution for this criterion, sometimes with a significant amplitude of temporal variations (obvious in the example of V2-B s or V5-A i , cf.Figs 11c and d, and 12c and   d).All these results tend to strongly disfavour the likelihood that variations in the amplitude of heterogeneity in CMB heat flux are responsible for the temporal variations in reversal frequency.
The criterion measuring the latitudinal distribution of CMB heat flux (variable q 20 ) also displays relatively moderate differences in correlation attributable to the initial conditions and length of numerical models.Among the models presenting a correct agreement at present with tomography, none leads to a clearly positive correlation.A clear majority of candidate models (25 out of the 33 satisfying the tomography agreement criterion) provide however a clearly anti-correlated evolution: models V1-B i , V2-B s or V5-A i are three examples (Figs 11a, c and d, and 12a, c and d).In addition, note that the time variability of this criterion (typically several tens of per cents variations, see Figs 11a, c, and d) is more significant than for the average heat flux for example.These results suggest that polar cooling may explain the palaeomagnetic reversal frequency with boundary driven core flow reinforcing the homogeneous dynamo meridional circulation leading to strongly turbulent conditions in the core (Olson & Amit 2014) while superchrons are periods of relatively weak polar cooling.Overall, this criterion is the one associated with the largest number of successful models.
Finally, for the criterion measuring the level of equatorial symmetry, only 5 models (out of the 33 satisfying the tomography agreement criterion) present significantly large correlation.As a consequence, equatorial symmetry is a less likely criterion for mantle control on reversal frequency.

Influence of mantle convection vigour
In the absence of precise viscosity estimates, constraints from the sinking speed of slabs ( Čížková et al. 2012) can be considered.All the models presented above roughly agree with (though slightly underestimate) this magnitude (see V slab in Fig. 9).To truly match the observed estimations of mean slab velocities, values in our models should be increased by a factor of 1.3-2 (Section 4.2).According to simple scaling laws, this would be obtained for Rayleigh numbers 1.5-3 times larger than the ones prescribed here.Furthermore, if the constraint on slab sinking velocities is relaxed and larger values of sinking velocities are permitted, even smaller viscosities (i.e.larger Rayleigh numbers) can be introduced for the lower mantle.
An assessment of this effect is shown in Fig. 13 where 3 larger values of the Rayleigh number were considered for two specific models.While the largest value considered (Ra = 4 × 10 8 ) most probably corresponds to an extreme endmember (maximum sinking velocities roughly one order of magnitude larger, ∼10 cm yr −1 ), we note that significant changes appear already at intermediate Ra values.First, as could be expected, the amplitude of time variability is in general enhanced when Ra is increased (Fig. 13a).Temporal variations in average CMB heat flux q 0 are ∼4-7 times larger for the endmember value Ra = 4 × 10 8 .While the likeli-hood that this specific criterion controls reversal frequency is low for the set of models described above due to very modest amplitude of time variability, such a conclusion does not hold if Ra is increased.The examination of the correlation between the time evolution of the mean CMB heat flux and the reversal frequency (Fig. 13b, top) indicates that no significantly positive coefficients are observed.However, as noted above, the introduction of other ingredients in the mantle could potentially shift these more pronounced peaks.The variable measuring the CMB heat flux heterogeneity, q * , does not present a larger amplitude of time variability as Ra is increased (but instead a globally flat evolution).Note however that q * is a relative quantity so that the absolute amplitude of lateral heterogeneity actually scales as q 0 and thus present a similar increase with Ra.Although variations in amplitude associated with latitudinal heterogeneity (q 20 ) and symmetry of the CMB heat flux do not display obvious trends (Fig. 13a, bottom two panels), some of the larger Ra models lead to interesting correlation results (Fig. 13b).A clearly positive correlation is obtained for the evolution of q 20 in two larger Ra models (corresponding to Ra = 4 × 10 7 ) , that is, polar cooling is maximal during CNS, in contrast to the conclusions obtained for Ra = 10 7 .This suggests that for larger values of Ra, it is conceivable that 'geographic control' (Olson & Amit 2014) can be considered as a viable cause to explain variations in reversal frequency.In addition, some of the larger Ra models yield acceptable symmetry correlations while the corresponding Ra = 10 7 models did not (Fig. 13b, bottom).In summary, in our models, if slab velocities are larger than those inferred by Van Der Meer et al. (2010) (∼12 mm yr −1 ), conclusions on the admissible criteria for mantle control on reversal frequency may differ significantly, especially with regard to the mean CMB heat flux, often considered as the most important criterion.

D I S C U S S I O N A N D C O N C L U S I O N S
We considered numerical models of mantle evolution in the last 450 Ma, forced by plate velocities history.Three main ingredients are varied: the lower-mantle rheology, the thickness of a dense basal layer and the initialization scheme.Most models tend to present a rather dynamic lower mantle for this period with the emergence of two thermochemical piles towards present-day (Fig. 7), although in a small minority of models two piles remain globally stationary during the last 100 s of Myr.Overall, initial conditions and the duration of numerical models is found to strongly affect the temporal evolution of the lower mantle in the last 100 s of Myr.At present-day, the composition field obtained in our models is found to correlate better with tomography in the lowermost mantle than the temperature field.Correlation coefficients between composition and tomography are statistically significant if the viscosity increase in the lower mantle is 30 or smaller.Criteria on the morphology of thermochemical piles (thickness of the piles, surface area of the CMB covered by anomalously dense mantle) tend to favour models with a discontinuous viscosity increase at the transition zone when compared to models with a gradual viscosity increase in the lower mantle.
In most of our models, the temporal variability of mean CMB heat flux is weak.When increasing the mantle convection vigour, Fig. 13 demonstrates that this variability increases.This may have important consequences for the likelihood of the reversals criteria.
Focusing on the CMB heat flux, we note that the pattern sometimes differs from the average temperature field in the 100 km mantle layer above.Furthermore, because the present-day composition better correlates with tomography than temperature does, the CMB heat flux pattern might differ significantly from the tomographic pattern.These findings motivate constraining geodynamo models with alternative CMB heat flux patterns which deviate from tomography based on lower-mantle dynamical scenarios (Amit et al. 2015).
Five criteria are proposed for CMB heat flux that can promote or inhibit reversals.We did not attempt to assign weights to these criteria.Even if such weights were combined to appraise the mantle models, the same overall score could be reached by different combinations of the individual criteria.We only address here whether such criteria can be used to discriminate among mantle models as all might inherently be pertinent for the real Earth.Perhaps more problematic, the correlation coefficients do not take into account critical thresholds below which there are no reversals at all and hence changes in CMB heat flux would have no effect on the zero reversal frequency, which is likely the case at the height of a superchron (Olson & Amit 2015).Again we do not introduce such thresholds simply because it is difficult to estimate them for the Earth's core.Nevertheless, the observation of specific trends among models raises the likelihood that a given criterion might indeed be relevant.
In our models, the evolution of the mean CMB heat flux very rarely correlates with the evolution of reversal frequency.Furthermore, initial conditions are found to strongly affect this correlation due to the small absolute values of time variations of the power out of the core.As a consequence, definite conclusions cannot be drawn for this criterion.These results oppose the general notion that reversal frequency is governed mostly by core convection vigour.
The amplitude of CMB heat flux heterogeneity is a well constrained criterion (i.e.weakly dependent on initial conditions) in the case of models involving a discontinuous viscosity increase at the transition zone.However generally, its expected correlation with reversal frequency is not observed.As with the mean heat flux, our results suggest that reversal frequency variability cannot be explained by changes in the amplitude of heat flux heterogeneity.
The evolution of the latitudinal distribution of CMB heat flux often anticorrelates with reversal frequency, especially for models involving a gradual viscosity increase in the lower mantle.This may suggest 'inertial control' on the geodynamo with polar cooling enhancing turbulent core conditions (Olson & Amit 2014).During superchrons the polar cooling is relatively weak.
The level of equatorial symmetry of the CMB heat flux pattern is found to provide in a few models a viable explanation for variations in the reversal frequency.A plausible symmetry correlation coefficient is found if the viscosity increase in the lower mantle is gradual.In these few models the break of Taylor columns in the outer core may provide a source of reversals.During superchrons, in these models, the level of equatorial symmetry is highest (as already noted by Biggin et al. 2012, in their models).
Overall, the most likely candidates among our preliminary models involve a large viscosity increase in the mantle.Models with a discontinuous viscosity increase at the transition zone tend to agree better at present with observations of seismic tomography, but models with a gradual viscosity increase provide better correlation coefficients with two criteria proposed to affect reversal frequency (latitudinal distribution of the CMB heat flux, and to a lesser extent level of equatorial symmetry).

P E R S P E C T I V E S
The choice of a specific model for plate motion history could affect the results.The global variance of the divergence for the period 0-120 Ma as well as the peak in plate speed during the CNS are comparable between the plate history models of Lithgow-Bertelloni & Silver (1998) and Seton et al. (2012).However, at earlier times, the plate history models of Zhang et al. (2010) and Furthermore, while Bull et al. (2014) conclude that the model of Domeier & Torsvik (2014) tends to favour the stability of the deep mantle African and Pacific thermochemical piles, this result is not reproduced with a similar plates set-up by Zhong & Rudolph (2015) where, instead, a dominant degree-1 structure is observed in the lowermost mantle for the early Pangea (∼300 Ma).In this context, an interesting path to investigate may be to systematically vary the plate motion histories prescribed as surface conditions of mantle models.Ideally, the palaeomagnetic frame should be entirely corrected based on TPW (or better, TPW could be evaluated from the mantle models depending on their structure in terms of density and rheology) and the solid rotation induced by TPW should then be applied to the CMB heat flux pattern (Biggin et al. 2012).Such efforts shall also focus on two characteristics of plate motion histories, namely divergence and net rotation, as these affect mantle dynamics in a substantially different manner.
While mantle models presented here aim at incorporating the effect of plate geometry and relative motion on mantle flow, one missing ingredient concerns the specific role of continents.These have been shown to impose a spatial heterogeneity in the heat flux boundary condition at the upper surface of the mantle (Jaupart & Mareschal 1999) that tends to increase the temperature (Grigné & Labrosse 2001) and to promote hot upwellings below them (while cold downwellings develop at their edges) (see Ricard 2015, for a brief synthesis).The locations of slabs dictated by plate motion is included the present models via imposed mechanical boundary conditions.Nevertheless, the effects on hot structures possibly developing very deep in the mantle are worth investigating.
Concerning the bulk mantle, while the radially averaged viscosity profiles introduced here overall encompass a reasonable range of variations with depth, other rheological aspects are simplified, partly owing to the semi-dynamic nature of the model where lithospheric plates do not emerge self-consistently from lithospheric rheology (cf.Section 2.1).Although still preliminary, attempts at constraining mantle models including self-consistent plate rheology (Bello et al. 2015) with plate motion models through a data assimilation method (Bocher et al. 2016) offer a very promising option to obtain a more realistic mantle model.
A potentially major ingredient missing in our mantle models concerns the influence of phase changes on radial transport.Most importantly, the phase change at 660 km induces lateral density gradients since the precise depth of the transition is modulated by temperature which might inhibit the penetration of slabs that are locally less dense than hotter ambient lower mantle (Ringwood & Irifune 1988).Early numerical models confirmed this proposition of layered convection with the possibility of intermittent global 'avalanches' in the lower mantle (e.g.Machetel & Weber 1991).More sophisticated descriptions (more realistic values for the Clapeyron slope and other thermodynamic and rheological properties) proposed more localized avalanches (Deschamps & Tackley 2009).Nevertheless, in the context of the present study, this effect could clearly add a specific temporal buffer between the plate motion history and the resulting CMB heat flux which has the potential to shift the peak of the latter from the CNS and thus revolutionize the correlation between CMB heat flux and reversal frequency.This potential temporal buffer is especially relevant to criteria presenting a significant amplitude of time variability, such as the lateral heterogeneity of CMB heat flux where notable peaks are observed although not coinciding correctly with the CNS in our models (see e.g.q * in Fig. 11).
Another simplification relates to the structural complexity of the deep lower mantle as revealed by seismology (again, see Lay 2015, for a synthesis).While LLSVPs are incorporated in our models as compositionally denser reservoirs, other aspects are simply omitted.These include the perovskite to post-perovskite phase change possibly associated with subducted lithosphere as the origin of the D discontinuity (cf.e.g.Nakagawa & Tackley 2008).Another aspect not considered is extremely thin regions with a distinct composition, possibly partially molten remnants of a basal magma ocean (Labrosse et al. 2007), as the origin of ULVZs (see Hernlund & McNamara 2015, for a synthesis).Such complexities are local/regional features but their effect on CMB heat flux and the geodynamo could be significant (Amit et al. 2015).While our results suggest that very large scale structures in the lower mantle control the geomagnetic dipole's stability (Fig. 7), the question remains to what extent such smaller scale structures also affect the reversal frequency.Finally, as indicated above, we have not considered here isochemical models that interpret LLSVPs in terms of large clusters of relatively thin thermal plumes (cf.e.g.Schuberth et al. 2009a).We even found higher correlations between the composition field and tomography than between temperature and tomography.Compositional gradients in the lower mantle might not lead to dynamical consequences if the associated density effect is too small, in spite of being detectable by seismology.Testing both groups of models (isochemical versus thermochemical) on the basis of the tomographic signature they should produce did not permit to decide whether the isochemical hypothesis can be ruled out (Bull et al. 2009;Schuberth et al. 2009a;Davies et al. 2012).
Clearly the set of criteria for reversal frequency that we gathered from the numerical dynamos literature (cf.Section 3) may be extended or refined if new studies unravel different aspects of the CMB heat flux ingredients that trigger reversals.Moreover, even with the existing criteria, here we relied solely on temporal trend agreement between reversal frequency and each criterion via a correlation coefficient, while amplitude variability may also be incorporated.Considering the amplitude of the mean CMB heat flux is plausible because the reversal frequency is expected to vary as √ q 0 (Olson & Amit 2014;Amit & Olson 2015).However, the relation between the reversal frequency and the amplitudes of the other criteria is not known, hence incorporating amplitude considerations for these criteria is premature.Finally, correlation in one criterion might be cancelled by anti-correlation in another.Ideally a quality factor that combines all criteria with appropriate weights should be defined, as was done for morphological criteria of the geomagnetic field (Christensen et al. 2010).Unfortunately, at the moment we do not know how to assign such weights so we avoid defining a quality factor.
Other ingredients apart from CMB heat flux may also affect reversal frequency, for example core convection type and inner core size.Olson et al. (2013) combined core convective and rotational evolution (including a growing inner core size) together with CMB heat flux pattern evolution to set up their dynamo models.They found that the CMB conditions affect reversal frequency much more than the core evolution.Nevertheless, it is worth considering such ingredients as additional criteria for reversal frequency.
Overall, this paper introduces a framework to independently test mantle convection models with palaeomagnetic observations and inferences from numerical dynamos.It is a first systematic attempt and improvements on the three building blocks in our approach (plate motion history models, mantle convection models, criteria on CMB heat flux inferred from numerical dynamos) can be envisioned.

Figure 3 .
Figure3.Radial profiles at the end (t = 0 Ma) of model V2-Bi: left-temperature as well as lateral average of composition (green); centre-viscosity; right-maximum ascending (red) and descending (blue) velocity.In the left and centre panels, minima are in blue, average in black and maxima in red.

Figure 4 .
Figure 4. Temperature (left) and composition (right) at various radii above the CMB at the end of model V2-B i (corresponding to present-day).Note that the colour scale is identical for all cross-sections in the case of composition while it varies for temperature.

Figure 5 .
Figure 5.Time evolution of the values associated with the criteria on CMB heat flux for model V2-B i (red) versus reversal frequency (blue).The CNS is highlighted in light blue.

Figure 6 .
Figure 6.Surface (dashed curves) and CMB (solid curves) heat power for the 54 numerical models.The two light blue rectangles indicate the CNS and the KRS.

Figure 8 .
Figure8.Linear correlation coefficient between a global tomography model for seismic shear wave velocity at the lowermost mantle(Masters et al. 2000) and either (a) the temperature field or (b) the composition field at the base of the mantle at present-day.Spherical harmonics decompositions of the three quantities have been truncated at degree max = 8.Colours indicate the viscosity model (as in Fig.2): V1 (red), V2 (blue), V3 (green), V4 (cyan), V5 (yellow) and V6 (purple) (see Table2).Suffix A, B or C (on the x-axis) indicates the model for the dense basal layer (cf.Table3).Filled circles denote models of group i, empty circles denote models of group ii and empty stars denote models of group s.The five solid horizontal lines denote confidence limits above which there is (from bottom to top) a 50 per cent, 40 per cent, 30 per cent, 20 per cent or 10 per cent chance that the correlation occurs by chance.

Figure 11 .
Figure 11.Time evolution of criteria on CMB heat flux (red) compared to the reversal frequency (blue) in the case of four selected numerical models: (a) V1-B i , (b) V2-B ii , (c) V2-B s , (d) V5-A i .The number in each panel (right upper corner) denotes the correlation coefficient between the two curves.The light blue rectangle denotes the CNS.

Figure 12 .
Figure 12.Four snapshots of the temperature right above the CMB (inversely proportional to the CMB heat flux) for the same numerical models as in Fig. 11.The temperature field has been expanded in spherical harmonics until degree max = 10.

Figure 13 .
Figure13.Influence of Rayleigh number Ra (i.e.average viscosity).Three additional models are performed for cases V2-B i (blue) and V5-B i (yellow): while the value of Ra for all models presented above is 10 7 , the new models correspond to 4 × 10 7 , 10 8 and 4 × 10 8 .(a) Absolute time variability of the four criteria in the period 160 Ma to present as a function of Ra.(b) Linear correlation coefficient between the time evolution of reversal frequency and the time evolution of criteria for the CMB heat flux (as in Fig.10) as a function of Ra for the same period: left-all models (two corresponding to results already presented above for V2-B i and V5-B i and six larger Ra models); right-only models whose present-day composition field correlates sufficiently well with tomography (with at least a 90 per cent confidence).

Table 1 .
Parameters for the mantle convection models (top) and for dimensionalization of the results (bottom).

Table 2 .
Parameter values for the six viscosity laws V • .

Table 3 .
McNamara et al. 2010ree models A, B and C for the dense basal layer.Each model involves only two parameters corresponding to the buoyancy number Bu and the initial dimensionless thickness d co (the value in km is indicated in parentheses).giventhetime scale of our numerical models that focus on the most recent geological history (<500 Ma).Observations of ULVZs (cf.e.g.McNamara et al. 2010) might reflect the presence of partial melts or indicate that compositional heterogeneities cannot be restricted to the 'normal mantle' versus piles dichotomy.Such important details as well as the presence of post-perovskite might affect the mantle dynamics in this region, and thereby the CMB heat flux.Nevertheless, as a rudimentary attempt to describe the largest scale structure, we introduce a single compositional field in our models accounting for denser material. importance