Geomagnetic semblance and dipolar-multipolar transition in top-heavy double-diffusive geodynamo models
ssubmitted to
Geophys. J. Int.
Geomagnetic semblance and dipolar-multipolar transition intop-heavy double-diffusive geodynamo models
Th´eo Tassin , Thomas Gastine and Alexandre Fournier Universit´e de Paris, Institut de Physique du Globe de Paris, CNRS, F-75005 Paris, France
12 January 2021
SUMMARY
Convection in the liquid outer core of the Earth is expected to be driven by thermal and chemi-cal density perturbations. For the sake of parsimony and in order to mitigate the computationalcost, the difference between the two buoyancy sources has been simply ignored in the majorityof geodynamo models published to date. The most common approach assumes that moleculardiffusivities should be replaced by eddy diffusivities in the turbulent liquid core of terrestrialplanets. Thermal and chemical perturbations are then aggregated into one single variable termedthe codensity. The main purpose of this study is to examine the impact of double-diffusive con-vection on magnetic field generation by means of three-dimensional global geodynamo models.We focus here on the so-called “top-heavy” regime of double-diffusive convection, when boththermal and compositional background gradients are destabilizing. We carry out a systematicparameter survey by computing numerical dynamo models spanning various fractionning ofbuoyancy sources. Using a linear eigensolver, we show that the onset of convection is facili-tated by the addition of a second buoyancy source. The critical onset mode is similar to classi-cal thermal Rossby waves observed with the codensity formalism. Using a rating parameter toquantify the morphological semblance of the models magnetic field with the geomagnetic field,we show that a good agreement can be attained for any partitioning of the convective inputpower. Next, we show that the transition between dipole-dominated and multipolar dynamosstrongly depends on the nature of the buoyancy forcing. Classical parameters expected to gov-ern this transition, such as the local Rossby number -a proxy of inertia over Coriolis force- orthe degree of equatorial symmetry of the flow, fail to capture the transition observed in our set ofdynamo models. A scale-dependent analysis of the force balance at work instead reveals that thedipole breakdown occurs when the ratio of inertia to Lorentz force at the dominant convectiveflow lengthscale reaches . , independently of the distribution of input power between ther-mal buoyancy and compositional buoyancy. The ratio of integrated kinetic to magnetic energy E k /E m appears to be a reasonable proxy of this force ratio. Given that E k /E m ≈ − − − in the Earth’s core, the geodynamo should operate far from the dipole-multipole transition. Ithence appears unlikely that the occurrence of geomagnetic reversals is related to dramatic andpunctual changes of the amplitude of inertial forces in the Earth’s core. Key words:
Dynamo: theories and simulations – Core – Numerical modelling – Compositionand structure of the core – Magnetic field variations through time
The exact composition of Earth’s core remains unclear but it is ad-mitted that it is mainly composed of iron and nickel with a mixtureof lighter elements in liquid state, such as silicon or oxygen (seeHirose et al. 2013, for a review). The ongoing crystallization of theinner core releases light elements and latent heat at the inner-coreboundary (ICB), while the mantle extracts thermal energy from theouter core at the core-mantle boundary (CMB). This combinationof processes is responsible for the joint presence of thermal andchemical inhomogeneities within the outer core.Convection with two distinct sources of mass anomaly is termed double-diffusive convection (e.g. Radko 2013). A key phys-ical parameter of double-diffusive convection is the Lewis num-ber Le , defined as the ratio of the thermal diffusivity κ T to thechemical diffusivity κ ξ . In the liquid core of terrestrial planets, Le reaches at least − (e.g. Loper & Roberts 1981; Li et al.2000).Under double-diffusive conditions, the increase of the back-ground density with depth does not necessarily imply the stabilityof the fluid in response to perturbations. As shown in Fig. 1, threeconfigurations can be considered: (i) the salt fingering regime, when the background thermal gradient a r X i v : . [ phy s i c s . g e o - ph ] J a n T. Tassin, T. Gastine and A. Fournier
Destabilizing background chemicalgradient Destabilizingbackground thermal gradient Non-linear hydrodynamicLinear onsetand weaklynon-linearDynamoPure thermalconvectionPure chemicalconvection
Salt fingersTop-heavySemi-convection
Manglik et al. (2010)Takahashi et al. (2019)Monville et al.(2019)Silva et al. (2019)Monville et al.(2019)Net et al. (2012)Silva et al. (2019) Breuer et al. (2010)Busse (2002)Manglik et al. (2010)Takahashi (2014)Net et al. (2012)Takahashi et al. (2019)Trümper et al. (2012)This studySilva et al. (2019)
Stabilizing background chemical gradientStabilizingbackground thermal gradient
Glatzmaier and Roberts (1996)Bouffard(2017)
No convection
Simitev et al. (2011)
Figure 1.
Schematic overview of previous studies on double-diffusive con-vection and dynamo action in terrestrial interiors, inspired by the regime di-agram of Ruddick (1983). The parameter space has two independent direc-tions, defined by the prescribed background temperature and compositiongradients, respectively. It is divided into four quadrants, since each gradi-ent can have a stabilizing or destabilizing effect on fluid motion. The blueand red half straight lines correspond to purely thermal and purely chemicalconvection, respectively (in which case there exists a unique destabilizingbackground profile). The publications that appear in the three quadrants ofdynamical interest are the following: Glatzmaier & Roberts (1996); Busse(2002); Breuer et al. (2010); Manglik et al. (2010); Simitev (2011); Net et al.(2012); Tr¨umper et al. (2012); Takahashi (2014); Bouffard (2017); Monvilleet al. (2019); Silva et al. (2019); Takahashi et al. (2019). Crosses denote dy-namo studies, triangles linear and weakly non-linear hydrodynamic studies,and squares non-linear hydrodynamic studies. ∇ T is stabilizing and the mean compositional gradient ∇ ξ isdestabilizing; (ii) the semi-convection regime, when ∇ T is destabilizing and ∇ ξ is stabilizing; (iii) the top heavy convection regime (also known as double-buoyant),when both ∇ T and and ∇ ξ are destabilizing.These three cases correspond to three quadrants in Fig. 1, whichis inspired by Ruddick (1983). Note that being located inside oneof these three quadrants does not necessarily guarantee that convec-tion occurs, since, for all configurations, the onset does not coincidewith the origin of the diagram: a critical contrast in temperature, orcomposition, or both, is required to trigger convective motions. Inaddition, this diagram does not account for the influence of back-ground rotation, or a magnetic field, on the onset (see Monvilleet al. 2019, for the impact of rotation on the onset in the salt fingersregime).Double-diffusive convection has been extensively studied inphysical oceanography (e.g. Radko 2013), as heat and salinity pro-vide two different sources of buoyancy for seawater, as well as instellar interiors (e.g. Spiegel 1972), where the molecular diffusivityof composition is only a tiny fraction of the diffusivity of temper-ature (e.g. Moll et al. 2016, and references therein). To our knowl-edge, a limited number of studies have been devoted to the anal-ysis of double-diffusive convection in the context of Earth’s core, or more generally in the context of the metallic core of terrestrialplanets. They are listed in Fig. 1, with possibly several occurrencesof a given study when several of the aforementioned regimes wereconsidered. It appears that the top-heavy configuration has been in-vestigated the most; this is also the configuration we shall focus onin this work.Previous studies devoted to the top-heavy regime are listed inthe top-right quadrant of Fig. 1, starting with hydrodynamic, non-magnetic studies. Considering the Le (cid:29) limit for a fluid fill-ing a rapidly-rotating annulus, Busse (2002) theoretically predictedthat the addition of chemical buoyancy facilitated the onset of stan-dard thermal convection, in two ways. First by lowering the criticalvalue of the Rayleigh number required to trigger the classical spi-ralling thermal Rossby waves of size and frequency proportionalto E / and E − / , respectively, where E = ν/ Ω D is the non-dimensional Ekman number, ν being the kinematic viscosity, Ω therotation rate and D the size of the fluid domain. Second, and moreimportantly, by enabling a second class of instability. The latter ischaracterized by a much lower onset, independent of E , a criti-cal length scale of the order of the size of the fluid domain, anda very small frequency, proportional to E . Busse (2002) referredto this class of instability as “nearly steady” (i. e. slow) convec-tion, and he argued that it could facilitate “immensely” convectionin Earth’s core, though with several caveats (see Busse 2002, fordetails). Simitev (2011) pushed the analysis further for various Le ,and stressed that the critical onset curve for convective instabil-ity in a rotating annulus forms disconnected regions of instabilitiesin the parameter space. In particular, the second family of modes(i.e. the slow modes) are stable whenever the compositional gra-dient is destabilizing. Tr¨umper et al. (2012), Net et al. (2012) andSilva et al. (2019) studied the onset of double-diffusive convec-tion in a spherical shell geometry. In the top-heavy regime with ≤ Le ≤ , Tr¨umper et al. (2012) confirmed that the additionof a secondary buoyancy source facilitates the convective onset.Net et al. (2012) showed that the properties of the critical onsetmode, such as its drift frequency or its azimuthal wave number,strongly depends on the fractioning between thermal and composi-tional buoyancy. Using a linear eigensolver, Silva et al. (2019) car-ried out a systematic survey of the onset of convection in sphericalshells for the different double-diffusive regimes. In the top-heavyconfiguration with Le = 25 , they showed that the convection on-set is characterized by an abrupt change between the purely ther-mal and the purely compositional eigenmodes depending on therelative proportion of the two buoyancy sources. In addition, theydemonstrated that the onset mode features the same asymptotic de-pendence on the Ekman number as classical thermal Rossby wavesover the entire top-heavy regime (top-right quadrant of Figure 1),thereby casting some doubt on the likelihood of the occurrence ofslow modes.Tr¨umper et al. (2012) performed a series of non-linear, mod-erately supercritical, rotating convection calculations at constant Le = 10 , for different proportions of chemical and thermal driving.To that end, they conducted a parameter survey, varying the chemi-cal and thermal Rayleigh numbers, to be defined below, while keep-ing their sum constant. This way of sampling the parameter space,however, does not guarantee that the total buoyancy input powerstays constant. This complicates the interpretation of their results,for instance regarding the influence of compositional and thermalforcings on the convective flow properties.Let us now turn our attention to the few self-consistentdynamo calculations in the top-heavy regime published to date(marked with a cross in the top-right quadrant of Fig. 1). The first ouble-diffusive geodynamo models integration was reported by Glatzmaier & Roberts (1996). This cal-culation was an anelastic, double-diffusive extension of the cele-brated Boussinesq simulation of the geodynamo by Glatzmaier &Roberts (1995). Glatzmaier & Roberts (1996) assumed enhancedand equal values of the diffusivities, i.e. Le = 1 . Consequently, thatsimulation did not exhibit stark differences with the purely thermalconvective model, except that the dipole did not reverse over thecourse of the simulated , years. The second, and last to date,numerical investigation of geodynamo models driven by top-heavyconvection has been conducted by Takahashi (2014). His modelswere based on the Boussinesq approximation with Le = 10 andrelatively large Ekman numbers ( E ≥ × − ). Top-heavy dy-namo simulations were also performed by Manglik et al. (2010)and Takahashi et al. (2019) in the context of modelling Mercury’sdynamo.Double-diffusive models of the geodynamo are the exceptionrather than the rule, essentially on the account of Occam’s razor. Ef-forts carried out in the community since the mid nineteen-ninetieshave been towards understanding the most salient properties of thegeomagnetic field using a minimum number of ingredients (e.g.Wicht & Sanchez 2019, for a review). To that end, the codensity for-malism introduced by Braginsky & Roberts (1995) is particularlyattractive: it assumes that the molecular values of κ ξ and κ T can bereplaced by a single turbulent transport property. Consequently, themass anomaly field can be described by a single scalar, termed thecodensity, that aggregates the two sources of mass anomaly. Thisapproach has the benefit of ( i ) removing one degree of freedomand ( ii ) mitigating the numerical cost by suppressing the scale sep-aration between chemical and thermal fields when Le (cid:29) . Withregard to Fig. 1, this amounts to restricting the diagram to either theblue or the red half straight line.The codensity formalism was quite successful in reproduc-ing some of the best constrained features of the geomagnetic fieldand its secular variation (e.g. Christensen et al. 2010; Aubert et al.2013; Schaeffer et al. 2017; Wicht & Sanchez 2019). Most geody-namo models actually assume that the diffusivity of the codensityfield equals the kinematic viscosity, yielding a Prandtl number ofunity. A remarkable property of the geodynamo that remains to beexplained satisfactorily from the numerical modelling standpoint,is its ability to reverse its polarity every once in a while, that isto go from a dipole-dominated state to another dipole-dominatedstate through a transient multipolar state (see e.g. Valet & Fournier2016, for a recent review of the relevant paleomagnetic data). Apossibility is that the geodynamo has been, at least punctually inits history, in a dynamical state that can enable the switch betweendipole-dominated and multipolar states to occur. A key questionthat follows is therefore: what are the physical processes that con-trol the transition between dipole-dominated dynamos and multipo-lar dynamos? This has been analyzed intensively numerically, start-ing with the systematic approach of Kutzner & Christensen (2002),who demonstrated that a stronger convective driving led to a dipolebreakdown, and that for intermediate values of the forcing, the sim-ulated field could oscillate between dipolar and multipolar states.Sreenivasan & Jones (2006) showed that an increasing role ofinertia (through stronger driving) perturbed the dominant Magneto-Archimedean-Coriolis force balance to the point that it led to a lessstructured and less dipole-dominated magnetic field. In the samevein, Christensen & Aubert (2006) assumed that the transition isdue to a competition between inertia and Coriolis force. They in-troduced a diagnostic quantity termed the local Rossby number Ro L = u rms / Ω L , as a proxy of this force ratio. Here u rms denotesthe average flow speed and L is an integral measure of the con- vective flow scale. Based on their ensemble of simulations, Chris-tensen & Aubert (2006) concluded that the breakdown of the dipoleoccurred above a critical value of about Ro L (cid:39) . , in what ap-peared a relatively sharp transition (see also Christensen 2010). Ifthis reasoning gives a satisfactory account of the numerical dataset,its extrapolation to Earth’s core regime raises questions (e.g. Oruba& Dormy 2014). Since the geomagnetic dipole reversed in the past,the numerical evidence collected so far (see Wicht & Tilgner 2010,for a review) suggests that the geodynamo could lie close to thetransition between dipolar and multipolar states. This implies that Ro L could be of the order of . for the Earth’s core. Geomagneticreversals should then reflect the action of a convective feature ofscale L of about
50 m . It is very unlikely that such a small-scaleflow could significantly alter a dipole-dominated magnetic field.According to Soderlund et al. (2012), the breakdown of thedipole is rather due to a decrease of the relative helicity of the flow.In numerical dynamo simulations, coherent helicity favors large-scale poloidal magnetic field through the α -effect (see Parker 1955)at work in convection columns and it therefore contributes activelyto the production and the maintenance of dipolar field (e.g. Olsonet al. 1999). Conversely, the dipolar field can promote a more heli-cal flow, with the Lorentz force enhancing the flow along the axis ofconvection columns, as shown by Sreenivasan & Jones (2011). Bymeasuring the integral force balance for dipolar and multipolar nu-merical dynamos, Soderlund et al. (2012) noticed that the Coriolisforce remained dominant, even in multipolar models. Accordingly,they suggested that, in their models, the modification of the flowstructure was rather controlled by a competition of second-orderforces, with the ratio of inertia to viscous forces as the parametercontrolling the transition.The role played by viscous effects was further stressed byOruba & Dormy (2014), who proposed that the transition fromdipolar to multipolar dynamos in the numerical dataset was con-trolled by a triple force balance between Coriolis force, viscosity,and inertia. A local Rossby number constructed using the viscouslengthscale, E / D , as opposed to the integral scale L discussedabove, carries the same predictive power in separating dipolar frommultipolar dynamo models (their Fig. 4). This argument was sub-sequently refined by Garcia et al. (2017), who included the Prandtlnumber dependence of the critical convective length scale whendefining the local Rossby number. From a mechanistic point ofview, Garcia et al. (2017) showed that the dipole breakdown wasnot necessarily correlated to a decrease of the relative helicity, butrather to a weakening of the equatorial symmetry of the flow. In-troducing the proportion of kinetic energy contained in this equato-rially symmetric component, they demonstrated that the transitioncould be satisfactorily explained by a sharp decrease of that quan-tity when the refined local Rossby number exceeded a value of . (their Fig. 3d). This would imply that the transition from a dipo-lar state to a multipolar state would essentially be a hydrodynamictransition. Discussing the implication of their results for the geo-dynamo, they clearly stated that the role of inertia was presumablyoverestimated in the numerical dataset that had been investigated sofar, stressing the need for stronger field dynamos, where the mag-netic field could possibly have an active role.The early numerical dataset admittedly contained a majorityof dynamos operating at large Ekman numbers, and relatively lowmagnetic Prandtl number P m . In the dynamos studied by Soder-lund et al. (2012), the convective flow was not dramatically alteredby the presence of a self-sustained magnetic field. Since then, alarge number of simulations have been published with lower valuesof E and comparatively larger values of P m (Yadav et al. 2016;
T. Tassin, T. Gastine and A. Fournier
Schaeffer et al. 2017; Schwaiger et al. 2019; Menu et al. 2020). Forthose strong-field dynamos, in the sense of a ratio of bulk magneticenergy to bulk kinetic energy larger than one, the magnetic fieldhas a significant impact on the flow, and on the dipolar-multipolartransition. Menu et al. (2020) reported simulations with a prevail-ing Lorentz force, that remain dipolar way beyond the supposedlycritical Ro L (cid:39) . value. Given that Earth hosts a strong-fielddynamo, it is then worth investigating whether the competition be-tween Lorentz force and inertia may actually lead to the transition.In order to shed light on this particular issue, and to fur-ther strengthen our understanding of double-diffusive dynamos, wehave performed a suite of novel dynamo models, comprising , , simulations with top-heavy, purely thermal, purely chemi-cal driving, respectively. The goal of this work is twofold: on theone hand, we will follow an approach similar to that of Takahashi(2014) and study to which extent the relative proportion of chemi-cal to thermal driving impacts the earth-likeness (in a morphologi-cal sense) of the simulated magnetic fields. On the other hand, wewill examine whether this proportion has an impact on the dipo-lar to multipolar transition. Takahashi (2014) reported a drop ofthe dipolarity when the relative contribution of thermal convectionto the total input power exceeds
65 % . We will analyze whetherthis was a fortuitous consequence of the cases he considered, andwhether we can find a more general rationale to explain the transi-tion, by careful inspection of the force balance at work.This paper is organized as follows: the derivation of the gov-erning equations and their numerical approximations are presentedin Section 2. The results are described in Section 3. We proceedby firstly investigating the onset of convection for top-heavy con-vection. The impact of the input power distribution on the Earth-likeness is then explored. Finally we examine the transition bew-teen dipolar and multipolar dynamos. Section 4 discusses the re-sults and their geophysical implications.
We operate in spherical coordinates ( r, θ, ϕ ) and consider a spher-ical shell of volume V o filled with a fluid delimited by the innercore boundary (ICB), located at the radius r i , on one side and bythe core mantle boundary (CMB), located at the radius r o , on theother side with r i /r o = 0 . . The shell rotates about the ˆz -axiswith a constant rotation rate Ω , where ˆz is the unit vector in thedirection of rotation. The equation of state ρ = ρ [1 − α T ( T − T ) − α ξ ( ξ − ξ )] , (1)describes how the density of the fluid ρ varies with temperature T and composition ξ . In this equation, α T and α ξ are the coefficientsof thermal and chemical expansion, T , ρ and ξ the average tem-perature, density and composition of lighter elements in the outercore.The properties of the fluid, including its kinematic viscosity ν , its magnetic diffusivity η , its specific heat C p , its chemical andthermal diffusivity ( κ ξ , κ T ), its coefficients of thermal and chem-ical expansion ( α T , α ξ ) are assumed to be spatially uniform andconstant in time. Due to the almost uniform density in the outercore, we also assume a linear variation of the acceleration of grav-ity g with radius. The physical and thermodynamical properties ofthe Earth’s core relevant for this study are given in Table 1. Convection of an electrically-conducting fluid gives rise to a mag-netic field B . The state of the fluid is then described by the ve-locity field u , the magnetic field B , the pressure p , the temper-ature T and the composition ξ . The equations governing the dy-namics of the flow under the Boussinesq approximation are castin a non-dimensional form. We adopt the thickness of the shell D = r o − r i as reference length scale and the viscous diffusion time D /ν as time scale. A velocity scale is then given by ν/ D . Com-position is scaled by ∂ξ = | ∂ r ξ ( r i ) |D , pressure by ρ ( ν/ D ) ,power by ν ρ / D and magnetic induction by √ ρ µ η Ω , where µ is the magnetic permeability of vacuum. Temperature unit isbased on the temperature gradient at r i , ∂T = | ∂ r T ( r i ) |D , or at r o , ∂T = | ∂ r T ( r o ) |D , depending on the thermal boundary con-ditions. Under the Boussinesq approximation, the equation for theconservation of mass is ∇ · u = 0 . (2)The dynamics of the flow is described by the Navier-Stokes equa-tion, expressed in the frame rotating with the mantle ∂ t u + u · ∇ u = − E ˆz × u + (cid:18) Ra T P r T + Ra ξ Sc ξ (cid:19) rr o ˆr − ∇ p + 1 EP m [( ∇ × B ) × B ] + ∇ u , (3)where ˆr is the unit vector in the radial direction. The time evolutionof the magnetic field under the magnetohydrodynamics approxima-tion is given by the induction equation ∂ t B = ∇ × ( u × B ) + 1 P m ∇ B with ∇ · B = 0 . (4)Finally, the evolution of entropy and composition are governed bythe similar transport equations ∂ t T + u · ∇ T = 1 P r ∇ T + h T , (5)and ∂ t ξ + u · ∇ ξ = 1 Sc ∇ ξ + h ξ , (6)where h T is a volumetric internal heating and h ξ a chemical volu-metric source.The set of equations (2 – 6) is governed by 6 dimensionlessnumbers: The Ekman number E expresses the ratio between vis-cous and Coriolis forces E = ν Ω D , the Prandtl (Schmidt) number between viscous and thermal (chem-ical) diffusivities P r = νκ T and Sc = νκ ξ and the magnetic Prandtl number between viscous and magneticdiffusivities P m = νη . The thermal and chemical Rayleigh numbers Ra T = α T g D ∂Tνκ T and Ra ξ = α ξ g D ∂ξνκ ξ , where g o is the gravitational acceleration at the CMB, measure the ouble-diffusive geodynamo models Table 1.
Physical and thermodynamical parameters of Earth’s outer core relevant for this study. The corresponding references are listed in the rightmostcolumn, which may comprise several entries if a bracket of values are provided.
Definition Symbol Value Reference (Lower bound - Upper bound)Inner radius r i . Dziewonski & Anderson (1981)Outer radius r o Dziewonski & Anderson (1981)Earth angular velocity
Ω 7 . × − rad · s − Gravitational acceleration at CMB g o .
68 m · s − Dziewonski & Anderson (1981)Core density at CMB ρ o × kg · m − Dziewonski & Anderson (1981)Specific heat C p (850 ±
80) J · kg − · K − Labrosse (2003)Heating power from core Q −
16 TW
Buffett (2015)Thermal conductivity at CMB k T −
100 W · m − · K − Konˆopkov´a et al. (2016) - Pozzo et al. (2013); Zhang et al. (2020)Thermal diffusivity κ T ( . − . ) × − m · s − Estimated using values of k T , ρ o and C p .Coefficient of thermal expansion α T (1 . ± . × − K − Labrosse (2003)Kinematic viscosity ν − m · s − Roberts & King (2013)Magnetic diffusivity η . − . · s − Pozzo et al. (2013) - Konˆopkov´a et al. (2016)Estimated magnetic field strength B rms × − T Gillet et al. (2010)Superadiabatic composition contrast ∆ ξ κ ξ ( × − – . × − ) m · s − Loper & Roberts (1981) - Li et al. (2000)Coefficient of chemical expansion α ξ . − . Braginsky & Roberts (1995) - Labrosse (2015)Estimated flow velocity u rms (0 . − . × − m · s − Finlay & Amit (2011)
Table 2.
Dimensionless control parameters. The two rightmost columns provide estimates of these parameters for Earth’s core and the values spanned by thesimulations computed in this study. Earth’s core values were estimated thanks to Tab. 1.
Name Symbol Definition Core This study
Ekman
E ν/ Ω D − − – − Thermal Rayleigh Ra T α T g D F Ti /νκ T − − Chemical Rayleigh Ra ξ α ξ g D F ξi /νκ ξ − − Magnetic Prandtl
P m ν/η (3 . − · − P r ν/κ T . – . Sc ν/κ ξ – Le κ T /κ ξ – Table 3.
Characteristic time scales for the Earth’s core. The values were estimated thanks to Tab. 1.
Name Symbol Definition Core
Rotation period τ Ω / Ω 4 h
Turnover time τ adv D /u rms −
250 yr
Magnetic diffusion time τ η D /η − yr Thermal diffusion time τ T D /κ T − yr Viscous diffusion time τ ν D /ν yr Chemical diffusion time τ ξ D /κ ξ – yr Table 4.
Output parameters of the numerical simulations and their estimates for Earth’s core
Name Symbol Definition Earth’s core This study Reference
Relative thermal convective power P % T see Eq. 15 −
70 % 0 −
100 %
Lister & Buffett (1995) -Takahashi (2014)Rossby
Ro u rms / Ω D (1 . − × − (0 . − × − Table 1Local Rossby Ro L u rms / Ω L . × − − .
09 0 . − . Davidson (2013) - Olson& Christensen (2006)Relative equat. symmetric kinetic energy ζ . − . . − . Aubert et al. (2017)Magnetic Reynolds
Rm u rms D /η (0 . − × − × Table 1Elsasser Λ B /µ ηρ o Ω 6 . −
39 0 . − × Table 1Dipolarity parameter f dip . − . . − Gillet et al. (2015) vigour of thermal and chemical convection. Note that the Lewisnumber Le discussed in the Introduction is the ratio of the Schmidtnumber to the Prandtl number, Le = κ T κ ξ = ScP r .
Table 2 provides estimates of these control parameters for theEarth’s core. These nondimensional number express the ratio of characteristic physical time scales E = τ Ω τ ν , P m = τ η τ ν and Le = τ ξ τ T , where τ Ω = 1 / Ω is the rotation period, τ ν = D /ν the viscous dif-fusion time, τ η = D /η the magnetic diffusion time, τ T = D /κ T the thermal diffusion time and τ ξ = D /κ ξ the chemical diffusiontime.Earth’s outer core evolves on a broad range of time scales, as T. Tassin, T. Gastine and A. Fournier one can glean from the inspection of Table 3. In particular, evenif the evolution of temperature and composition are governed bysimilar transport equations (see Eq. 5 and 6), thermal diffusionis much more efficient than chemical diffusion, which causes theLewis number to greatly exceeds unity (see Tab. 2). This impliesthat the typical lengthscale of chemical heterogeneities is possiblyseveral orders of magnitude smaller than the length scale of thermalheterogeneities. The values of the control parameters spanned bythe simulations presented in this work are discussed in section 2.6.
The inner core is growing and is ejecting light elements becauseof its crystallization, which in principle yields coupled boundaryconditions for temperature and composition (e.g. Glatzmaier &Roberts 1996; Anufriev et al. 2005). For the sake of simplicity, ther-mal and chemical boundary conditions are considered as decoupledhere, as in e.g. Takahashi (2014). We consider two different setupsto explore the impact of varying the thermal boundary conditions: (i)
Fixed fluxes: thermal and composition fluxes are imposed at bothboundaries with (cid:26) ∂ r T ( r i ) = − ,∂ r ξ ( r i ) = − , ∂ r T ( r o ) = 0 ,∂ r ξ ( r o ) = 0 . (7) (ii) Hybrid: temperature is fixed at the ICB while the temperature fluxis imposed at the CMB, the boundary conditions on chemical com-position are the same as in the previous setup (cid:26) T ( r i ) = 0 ,∂ r ξ ( r i ) = − , ∂ r T ( r o ) = − ,∂ r ξ ( r o ) = 0 . (8)No-slip mechanical boundary conditions are used at both bound-aries. The mantle is assumed to be insulating, such that the mag-netic field at the CMB has to match a source-free potential field.The inner core is treated as a rigid electrically-conducting spherewhich can freely rotate around the ˆz -axis (e.g. Hollerbach 2000;Wicht 2002). Its rotation is a response to the viscous and magnetictorques exerted by the outer core on the inner core. The conductiv-ity of the outer core is assumed to be equal to that of the outer core,and its moment of inertia is calculated by using the same density asthe liquid outer core. We solve the system of equations (2– 6) using the open-source geo-dynamo code MagIC (cid:63) (Wicht 2002; Gastine et al. 2016). This codehas been validated against a benchmark for double-diffusive con-vection initiated by Breuer et al. (2010).The solenoidal vectors u and B are decomposed in poloidaland toroidal potentials (cid:26) u ( r , t ) = ∇ × ∇ × [ W ( r , t ) ˆr ] + ∇ × [ Z ( r , t ) ˆr ] , B ( r , t ) = ∇ × ∇ × [ G ( r , t ) ˆr ] + ∇ × [ H ( r , t ) ˆr ] , where r is the radius vector. The new unknowns are then W , Z , G , H , T , ξ and p .Each of these scalar fields is expanded in spherical harmonicsto maximum degree and order (cid:96) max in the horizontal direction. The (cid:63) https://magic-sph.github.io/ spherical harmonic representation of the magnetic poloidal poten-tial G reads G ( r, θ, ϕ, t ) (cid:39) (cid:96) max (cid:88) (cid:96) =0 (cid:96) (cid:88) m = − (cid:96) G (cid:96)m ( r, t ) Y m(cid:96) ( θ, ϕ ) where G (cid:96)m ( r, t ) is the coefficient associated to Y m(cid:96) , the spheri-cal harmonic of degree (cid:96) and order m . The non-linear terms arecalculated in physical space. The open-source SHTns library † (seeSchaeffer 2013) is used to compute the forward and inverse spectraltransforms on the unit-sphere.In the radial direction, MagIC uses either a finite differencescheme, or a Chebyshev collocation method (see Boyd 2001). Thefinite difference grid, whose number of points is denoted by N r , isregularly spaced in the bulk of the domain, and follows a geometricprogression near the boundaries (see Dormy et al. 2004). Whenthe collocation approach is selected, Chebyshev polynomials aretruncated at degree N and the N r collocation points r k are definedby ∀ k ∈ [[1 , N r ]] , r k = 12 [( r o − r i ) x k + r o + r i ] x k = cos (cid:20) ( k − πN r − (cid:21) . (9)Due to the particular choice of spatial grid given by the equationabove, the transforms between physical grid and Chebyshev repre-sentation are carried out by fast discrete cosine transform (see Presset al. 1992, chapter 12). This discretisation yields a point densifica-tion close to the boundaries which could impose severe time steprestrictions when the magnetic field is strong (see Christensen et al.1999). To mitigate this effect, we adopt the mapping proposed byKosloff & Tal-Ezer (1993) and replace x k in Eq. (9) by ∀ k ∈ [[1 , N r ]] , X k = arcsin( αx k )arcsin( α ) with α ∈ ]0 , , where α is the mapping coefficient. To maintain the spectral con-vergence of the simulation α has to verify α ≤ α max = (cid:26) cosh (cid:20) | ln( (cid:15) ) | N r − (cid:21)(cid:27) − , (10)where (cid:15) is the machine precision.Figure 2 shows the minimum grid spacing ∆ r min as a functionof N r for a regular grid, the finite difference grid with geometricalclustering near the boundaries, and two collocation grids, with α → (Gauss-Lobatto grid) and α = α max . Because ∆ r min ∼ N − r when using the classical Gauss-Lobatto grid, adopting the mappingby Kosloff & Tal-Ezer (1993) yields a possible increase of ∆ r min by a factor − when N r (cid:38) . The time step size could inprinciple rise by a comparable amount, should it be controlled bythe propagation of Alfv´en waves in the vicinity of the boundaries(Christensen et al. 1999). Using the finite difference method en-ables even larger grid spacing and hence possible additional gainin the time step size. This speed-up shall however be mitigated bythe fact that N r has to be increased by a factor . to when usingfinite differences to achieve an accuracy comparable to that of theChebyshev collocation method (see Christensen et al. 2001; Matsuiet al. 2016). G (cid:96)m is then expanded in truncated Chebyshev series G (cid:96)m ( r k , t ) (cid:39) (cid:114) N r − N − (cid:88) (cid:48)(cid:48) n =0 G (cid:96)mn ( t ) T n ( r k ) , (11) † https://bitbucket.org/nschaeff/shtns ouble-diffusive geodynamo models N r − − − − ∆ r m i n ∼ N − r ∼ N − r ∼ N − r ∼ N − r Regular grid ( α → α = α max Chebyshev-Gauss-Lobatto ( α → Figure 2.
Minimal grid spacing ∆ r min between two radial points as a func-tion of the number of collocation points N r for a regularly-spaced grid, thegrid used when finite differences are employed, the collocation grid with themapping by Kosloff & Tal-Ezer (1993) with a mapping coefficient α max (Eq. 10) and the standard Gauss-Lobatto grid. with G (cid:96)mn ( t ) (cid:39) (cid:114) N r − N r (cid:88) (cid:48)(cid:48) n =1 G (cid:96)m ( x k , t ) T n ( x k ) , where the double primes on the summations indicate that the firstand the last indices are multiplied by one half (see Glatzmaier1984). In the above equations, T n ( x k ) is the n -th order Chebyshevpolynomial at the collocation point x k defined by T n ( x k ) = cos[ n arccos( x k )] = cos (cid:20) πn ( k − N r − (cid:21) . Further details on spherical harmonics and Chebyshev polynomialsexpansion can be found in Tilgner (1999) and Christensen & Wicht(2015).Once the spatial discretisation has been specified, the set ofequations (2–6) complemented by the boundary conditions (seeEq. 7 and 8) forms a semi-discrete system where only the timediscretisation remains to be expressed. As an example, the timeevolution of the poloidal potential for the magnetic field G (cid:96)m ( r k ) (see Eq. 11) can be written as an ordinary differential equation d G (cid:96)m d t ( r k , t ) = E [ u , B ] + I [ G (cid:96)m ] ,G (cid:96)m ( x k , t ) = G (cid:96)m ( x k ) , (12)where G (cid:96)m ( r k ) is the initial condition, E a non linear-function of u and B and I a linear function of G (cid:96)m . The above equation serves asa canonical example of the treatment of the different contributions:the non-linear terms are treated explicitly (function E ) while the re-maining linear terms are handled implicitly (function I ). In MagIC,several implicit/explicit (IMEX) time schemes are employed totime advance the set of equations (2–6) from t to t + δt : (i) a combination of a Crank-Nicolson for the implicit terms anda second-order Adams-Bashforth for the explicit terms calledCNAB2 (see Glatzmaier 1984), (ii) two IMEX Runge-Kutta : PC2 (see Jameson et al. 1981) andBPR353 (see Boscarino et al. 2013).CNAB2 has been commonly used in geodynamo models sincethe pioneering work of Glatzmaier (1984). IMEX Runge-Kutta schemes have been rarely employed in the context of geodynamomodels (Glatzmaier & Roberts 1996), rapidly-rotating convectionin spherical shells (Marti et al. 2016) or quasi-geostrophic mod-els of 2-D convection (Gastine 2019). For IMEX Runge-Kuttaschemes, s substages are solved to time-advance Eq. (12) from t to t + δtG i(cid:96)m ( r k ) = G (cid:96)m ( r k , t ) + δt i (cid:88) j =0 (cid:16) a E i,j E j + a I i,j I j (cid:17) , (13)where i ∈ [[0 , s ]] , G i(cid:96)m is the intermediate solution at substage i , E j = E [ B , u ]( t + c E j δt ) and I j = I [ G (cid:96)m ]( t + c I j δt ) . G (cid:96)m ( r k , t + δt ) is then given by G (cid:96)m ( r k , t + δt ) = G (cid:96)m ( r k , t ) + δt s (cid:88) j =0 (cid:16) b E j E j + b I j I j (cid:17) (14)In the above equations, the matrices a I and a E and the vectors b I , b E , c E and c I form the so-called Butcher tables of the IMEXRunge Kutta schemes given in Appendix A. Since the last lines of a E , I are equal to b E , I for PC2 and BPR353, the last operation toretrieve G (cid:96)m ( r k , t + δt ) (Eq. 14) is actually redundant with the lastsub-stage (see Ascher et al. 1997). IMEX Runge-Kutta schemesrequire more computational operations to time advance the set ofequations(2–6) from t to t + δt than CNAB2. However, they allowlarger time step sizes that compensate for this extra numerical cost(Marti et al. 2016). For each diagnostic quantity f , we adopt in the following overbarsfor time averaging and triangular brackets for spatial averaging ¯ f = 1 t avg (cid:90) t + t avg t f ( t ) d t and (cid:104) f (cid:105) V = 1 V (cid:90) V f ( r , t ) d V, where t avg corresponds to the averaging time. The magnetic E m and kinetic E k energies are given by E m ( t ) + E k ( t ) = 12 (cid:20)(cid:90) V o + V i B ( r , t ) EP m d V + (cid:90) V o u ( r , t ) d V (cid:21) , where V i is the inner core volume. By multiplying the Navier-Stokes equation (3) by u and the induction equation (4) by B , weobtain the following power balance dd t ( E k + E m ) ( t ) = P T ( t ) + P ξ ( t ) − D ν ( t ) − D η ( t ) . In the case of double-diffusive convection, the energy is providedby chemical and thermal buoyancy power P ξ and P T defined by P ξ ( t ) = V o (cid:28) Ra ξ Sc ξ ( r , t ) rr o u r ( r , t ) (cid:29) V o P T ( t ) = V o (cid:28) Ra T P r T ( r , t ) rr o u r ( r , t ) (cid:29) V o and dissipated by viscous and Ohmic dissipations D ν and D η givenby T. Tassin, T. Gastine and A. Fournier D ν ( t ) = V o (cid:10) [ ∇ × u ( r , t )] (cid:11) V o D η ( t ) = ( V o + V i ) (cid:28) [ ∇ × B ( r , t )] EP m (cid:29) V o + V i . Once a statistically steady state has been reached, the inputbuoyancy powers should compensate the Ohmic and viscous dissi-pations. Following King et al. (2012), to assess the consistency ofthe numerical computations, we measure the time-average differ-ence between input and output powers ∆ P ∆ P = 100 × P T + P ξ − D ν − D η P T + P ξ . We made sure that this difference remained lower than . forall the simulations reported in this study.For each simulation, the total convective power P tot and therelative thermal convective power P % T are defined by P tot = P T ( t ) + P ξ ( t ) and P % T = P T ( t ) P T ( t ) + P ξ ( t ) × . (15) P % T hence vanishes for a purely chemical forcing and is equal to
100 % for a purely thermal forcing.Following Christensen & Aubert (2006) and Schwaiger et al.(2019), we introduce two quantities to characterise the typical flowlengthscale. The integral scale L already discussed in the introduc-tion is obtained from the time-averaged kinetic energy spectrum L = π E k ( t ) (cid:88) (cid:96) (cid:96) u (cid:96) ( t ) · u (cid:96) ( t ) , where u (cid:96) · u (cid:96) / is the kinetic energy contained in spherical har-monic degree (cid:96) , while the dominant lengthscale (cid:98) (cid:96) is defined as thepeak of the poloidal kinetic energy spectra (Schwaiger et al. 2019,2020) (cid:98) (cid:96) = argmax (cid:96) (cid:104) E (cid:96)k,p ( t ) (cid:105) , where E (cid:96)k,p is the contribution of spherical harmonic degree (cid:96) tothe total poloidal kinetic energy.In order to explore the impact of the equatorial symmetry ofthe flow on the dipole-multipole transition, we consider the relativeequatorially-symmetric kinetic energy ζ introduced by Garcia et al.(2017) ζ = E s , NZ k E NZ k , (16)where E NZ k is the kinetic energy contained in the non-zonalflow and E s , NZ k the kinetic energy contained in the equatorially-symmetric part of the non-zonal flow.We measure the mean convective flow amplitude either by themagnetic Reynolds number Rm or by the Rossby number Ro de-fined by Rm = (cid:112) E k ( t ) P m = τ η τ adv , Ro = (cid:112) E k ( t ) E = τ Ω τ adv , where τ adv = D /u rms is the characteristic turnover time and . . . . . . f dip C oun t o f s i m u l a ti on Multipolar Dipolar
Figure 3.
Distribution of the dipolar fraction f dip for the numericalsimulations of this study. Multipolar simulations are defined as having f dip < . , and they will be marked by a cross in subsequent figures. Dipo-lar simulations ones will be displayed using a circle. The vertical dashedline marks the f dip = 0 . limit between dipolar and multipolar simula-tions. u rms the root mean square flow velocity. Following Christensen& Aubert (2006), we define a local Rossby number by Ro L = Ro D L .
The dipolar character of the CMB magnetic field is quantifiedby its dipolar fraction f dip , defined as the ratio of the axisymmetricdipole component to the total field strength at the CMB in sphericalharmonics up to degree and order .Figure 3 shows the statistical distribution of f dip for the sim-ulations reported in this study. Two distinct groups of numericalsimulations separated by a gap at f dip ≈ . are visible in this fig-ure. A magnetic field is considered as dipolar when f dip > . andmultipolar otherwise. This bound differs from the original thresh-old of f dip = 0 . considered by Christensen & Aubert (2006),but it is found to better separate the two types of dynamo modelscontained in our dataset. Note that the same bound of 0.5 was re-cently chosen by Menu et al. (2020) in their study. The magneticfield amplitude is measured by the Elsasser number ΛΛ = √ E m . Finally, transports of heat and chemical composition are quan-tified by using the Nusselt Nu and the Sherwood Sh number de-fined by (see Goluskin 2016, chapter 1) Nu = ∆ T ∆ T and Sh = ∆ ξ ∆ ξ , where ∆ T , ∆ ξ are the temperature and composition differencesbetween the ICB and and the CMB, and the subscript stands forthe background conducting state. For both the fixed fluxes and hy-brid configurations (recall section 2.3 above), we obtain a back-ground composition contrast given by ∆ ξ = η ( η + 2)2( η + η + 1) , (17)where η = r i /r o is the radius ratio. The background temperaturedrop ∆ T depends on the imposed boundary conditions. For thefixed flux setup, it reads ∆ T ff0 = η ( η + 2)2( η + η + 1) = ∆ ξ , (18) ouble-diffusive geodynamo models while it becomes ∆ T hyb0 = 1 η , (19)for the hybrid setup.Table 4 gives the definition of most of these integral diagnos-tics and provides estimates for Earth’s core, along with the bracketof values obtained in the numerical dataset presented here. We compute simulations varying the Ekman number, the mag-netic Prandtl number, the thermal and the chemical Rayleigh num-bers. The properties of this dataset are listed in Tab. A1. The lessturbulent simulations have been initialized with a strong dipolarfield and a random thermo-chemical perturbation. Their final stateshave been used as initial conditions for the more turbulent simu-lations, in order to shorten their transients. Three different Ekmannumbers are considered in this study : − , × − and − . Ra T has been varied between and × and Ra ξ between and . × to study the influence of the convective forcing andspan the transition between dipole-dominated and multipolar dy-namos. We adopt P r = 0 . and Sc = 3 (i.e. Le = 10 ) for a bettercomparison with previous studies (e.g. Takahashi 2014) and to mit-igate the computational cost associated with large Lewis numbers. P m varies between . and , depending on the Ekman number,in order to maintain Rm > . The numerical models were inte-grated for at least 20 % of a magnetic diffusion time τ η for the mostturbulent (and demanding) ones, and for more than one τ η for theothers, in order to ensure that a statistically steady state had beenreached. Figure 4 shows the location of the computed numerical simu-lations for the three considered Ekman numbers in the parameterspace ( G T , G ξ ) defined by G T = 1 + Gr T E / and G ξ = 1 + Gr ξ E / , where Gr T ( ξ ) corresponds to the thermal (chemical) Grasshofnumber, Gr T = Ra T P r and Gr ξ = Ra ξ Sc .
Adding to Gr T ( ξ ) E / in the above equations allows us to uselogarithmic scales in the top row of Fig. 4.We determine the onset of convection using the open-sourcesoftware SINGE ‡ which computes linear eigenmodes for incom-pressible, double-diffusive fluids enclosed in a spherical cavity(see Schaeffer 2013; Vidal & Schaeffer 2015; Kaplan et al. 2017;Monville et al. 2019). For a fixed Ra T ( Ra ξ ), the code solves thegeneralized eigenvalue problem formed by the linearized Navier-Stokes and transport equations. It seeks eigenmodes f of the form f ( t, r, θ, ϕ ) = ˆ f ( r, θ ) exp[ i ( mϕ − ωt )] , where ˆ f is a function of r and θ , m is the azimuthal wave num-ber and ω the complex angular frequency. Starting at a specific ‡ https://bitbucket.org/vidalje/singe ( Ra T , Ra ξ ), the critical mode is determined by varying one of theRayleigh numbers (keeping the other fixed) in order to obtain an ω with a vanishing imaginary part. The onset mode is then char-acterized by its critical chemical and thermal Rayleigh numbers( Ra cξ , Ra cT ), its critical azimuthal wave number m c and its (real)angular drift frequency ω c . Connecting the ( Ra cξ , Ra cT ) pairs thatwe collected with SINGE gives rise to the critical curves plotted inFig. 4. In each panel, the intersection of these curves with the x -axis(resp. y -axis) defines the critical Rayleigh number for purely ther-mal (resp. chemical) convection. Underneath the critical curves, thegrey-shaded areas are regions of the parameter space where thermaland chemical perturbations are unable to trigger a convective flow.We note that the shape delimited by the critical curves in thebottom row of Fig. 4 is not rectangular, as would be the case if thetwo sources of buoyancy were independent. Instead, we observe adecrease of the critical G cξ when G cT increases for the three differentEkman numbers. As previously reported by Busse (2002), Tr¨umperet al. (2012) and Net et al. (2012), and discussed in the Introduction,this demonstrates that the addition of a second buoyancy sourcefacilitates the onset of convection.Starting from G cT = 0 and following the critical curve for eachEkman number, m c grows until one reaches the upper right ”cor-ner” of the onset region and then decreases to a value comparableto the starting G cT = 0 m c value when G cξ tends to zero. We ob-serve that m c is nearly constant on the vertical branches, while itincreases much faster on the horizontal ones, as already reportedby Tr¨umper et al. (2012).As can been seen in Figs. 4(e) and 4(f), the shaded shape ismuch wider with fixed-flux boundary conditions than with hybridboundary conditions. This comes from the difference in the temper-ature contrasts of the background conducting states. An adequateway to compare both setups resorts to using diagnostic Grasshofnumbers (cid:103) Gr T = Gr T ∆ T , (cid:103) Gr ξ = Gr ξ ∆ ξ , (20)that are based on the temperature (composition) contrasts ∆ T ( ∆ ξ ) of the reference conducting state instead of the control ther-mal (chemical) Grasshof numbers Gr T ( Gr ξ ) (see Johnston & Do-ering 2009; Goluskin 2016). Using Eqs. (18) and (19), the temper-ature difference between ICB and CMB for the conducting statesfor both setups reads ∆ T hyb0 ∆ T ff0 = 2( γ + γ + 1) γ ( γ + 1) ≈ . , which is in good agreement with the actual ratio of critical G cT ob-served in panels (e) and (f) of Fig. 4.The analysis of the onset of double-diffusive convection be-comes even more straightforward if one adopts the formalism in-troduced by Silva et al. (2019). This framework rests on two pa-rameters: first, the diagnostic effective Grasshof number (cid:103) Gr c = (cid:114)(cid:16) (cid:103) Gr cT (cid:17) + (cid:16) (cid:103) Gr cξ (cid:17) , (21)and second the Grasshof mixing angle Θ such that cos Θ = (cid:103) Gr cT (cid:103) Gr c , sin Θ = (cid:103) Gr cξ (cid:103) Gr c . (22)The pair ( (cid:103) Gr c , Θ) can be interpreted as the polar coordinates of thecritical onset mode in the ( (cid:103) Gr T , (cid:103) Gr ξ ) Cartesian parameter space. T. Tassin, T. Gastine and A. Fournier G ξ a) E = − Fix. flux b) E = × − Hyb.
Multipolar Dipolar Hybrid Fixed flux Onset of convection10 c) E = − (*) (**)
10 20123456 G ξ d) G T e) f) 2345678910111213141516171819 m c Figure 4.
Linear onset of top-heavy convection in ( Gr T E / , Gr ξ E / ) parameter space, where Gr T , Gr ξ , and E are the thermal Grasshof,chemical Grasshof and Ekman numbers, respectively, for the three Ekman numbers E considered in this study: − (left column), × − (centercolumn), and − (right column). Critical curves correspond to the edges of the gray shaded areas. Dark gray areas were obtained for fixed-flux boundaryconditions and light gray areas for hybrid boundary conditions, the latter present only for E = 3 × − and E = 10 − . The bottom panels (d), (e) and (f)show zoomed-in insets of upper panels (a), (b) and (c). The edges of the gray areas, which define the critical curves, connect discs whose color defines thecritical azimuthal wavenumber m c . The top row (in logarithmic scale in both directions) also features the location in parameter space of the simulationscomputed in this study. Circles (resp. crosses) represent dipolar (resp. multipolar) simulations. Circles and crosses with gray (resp. black) edges correspond tofixed-flux (resp. hybrid) boundary conditions. Simulations (*) and (**) are reference simulations discussed in detail in the text (see also Tab. A1). A mixing angle
Θ = 0 hence corresponds to purely thermal con-vection, while
Θ = π/ corresponds to purely chemical convec-tion. Figure 5 shows the critical effective Grasshof number (cid:103) Gr c ,the critical azimuthal wave number m c and the critical angulardrift frequency | ω c | , multiplied in each instance by their expectedasymptotic dependence on the Ekman number for purely thermalconvection, as a function of the Grasshof mixing angle Θ . Adopt-ing a diagnostic effective Grasshof number conveniently enablesthe merging of the onset curves associated with the two sets of ther-mal boundary conditions (fixed-flux and hybrid).The onset curves can be separated into two branches: ( i ) From Θ = 0 up to Θ ≈ π/ , the onset mode almost behaves as apure low- P r thermal mode with little change in m c E / and alarge drift speed; ( ii ) a sharp transition to another kind of onsetmode, reminiscent of the P r (cid:38) convection onset, is observed for Θ (cid:38) π/ . The latter is characterized by a smaller drift frequency,a higher azimuthal wavenumber and a lower effective Grasshofnumber (cid:103) Gr c . The critical azimuthal wavenumber reaches its max-imum for Θ ≈ π/ before gradually decreasing to reach a valuecomparable to that expected for purely thermal convection towards Θ = π/ .The Ekman dependence is almost perfectly captured by theasymptotic scalings E − / and E − / for the critical wavenum-ber m c and the drift frequency ω c . The case of mixing angles Θ ≥ π/ with E = 10 − constitutes an exception to thisrule with a sharp drop to a constant critical wavenumber m c = 2 (see the small kink in the upper left part of Fig. 4d). The depen-dence of (cid:103) Gr c on the Ekman number shows a more pronounced de-parture from the leading-order asymptotic scaling Gr c ∼ E − / .As shown by Dormy et al. (2004) for differential heating (see e.g.their Fig. 5), the higher-order terms in the asymptotic expansionof Ra cT as a function of the Ekman number are still significantfor E > − (see also Schaeffer 2016). Given the range of Ek-man numbers considered for this study, it is not suprising that theasymptotic scaling for Gr c is not perfectly realized yet.In summary, the onset of double-diffusive convection in thetop-heavy regime takes the form of thermal-like drifting Rossby-waves, the nature of which strongly depends on the fraction be-tween chemical and thermal forcings. This confirms the results pre-viously obtained by Silva et al. (2019) (their Fig. 9). We will now focus on the simulation marked by an asterisk (*) inTable A1 and in Figure 4 to highlight specific double-diffusive con-vection features.This simulation corresponds to Ra T = 3 . × , Ra ξ = 4 . × , E = 10 − and P m = 0 . with hybrid bound-ary conditions. The thermal convective power amounts on average ouble-diffusive geodynamo models g G r c × E / a) E = − Hybrid E = × − Fixed flux E = − . . . . m c × E / b) π / π / π / π / π /
16 3 π / π / π / Θ . . . . . . . . | ω c | × E / c) Figure 5. (a) Critical effective Grasshof number multiplied by E / as afunction of the Grasshof mixing angle Θ . (b) Critical azimuthal wavenum-ber m c multiplied by E / as a function of Θ . (c) Critical drift frequency ω c multiplied by E / as a function of Θ . Symbol colors correspond tothe Ekman number, open-symbols to hybrid boundary conditions and filledsymbols to fixed-flux boundary conditions. for
46 % of the total input power. Since the local Rossby number Ro L reaches . , this simulation is expected to operate in a pa-rameter regime close to the transition between dipolar and multi-polar regimes ( Ro L ≈ . ) put forward by Christensen & Aubert(2006). This high value of Ro L also indicates the sizeable roleplayed by inertia in the force balance of this simulation. Figure 6shows 3-D renderings of several fields extracted from a snapshottaken over the course of the numerical integration: (a) temperatureperturbation, (b) composition, (c) zonal velocity and magnitude ofthe velocity vector, (d) radial magnetic field and magnitude of themagnetic field vector. We chose the radius of the inner and outerspheres of these renderings to place ourselves outside thermal andcompositional boundary layers. Convection is primarily driven bythe chemical composition flux at the ICB. Because of the contrastin diffusion coefficients ( Le = 10 ), compositional plumes developat a much smaller scale than that of thermal plumes (see Fig. 6a andFig. 6b). Having Le = 10 also induces a chemical boundary layermuch thinner than the thermal one, as illustrated by the Sherwoodnumber Sh being five times larger than the Nusselt number Nu inthis case (44.8 versus 8.0, see Tab. A1). The emission of a thermalplume is likely triggered by an impinging chemical plume. Accord- ingly, one would expect temperature fluctuations to be enslaved tocompositional fluctuations, which may explain the outstanding spa-tial correlation between temperature and composition noticeable inFig. 6(a) and Fig. 6(b). This correlation was already reported byTr¨umper et al. (2012) for predominantly thermal convection.In typical dipole-dominated dynamos, the velocity field fea-tures extended sheet-like structures that span most of the core vol-ume (see Yadav et al. 2013; Schwaiger et al. 2019). Because of thestrong forcing that characterizes the reference simulation, the con-vective flow is organised on smaller scales, which likely reflectsthe sizeable amplitude of inertia. The magnetic field is dominatedby its axisymmetric dipolar component (see Fig. 6d), despite thesupposed proximity of the simulation with the dipolar-multipolartransition zone. Inspection of Table A1 indicates that this snapshot-based observation can in fact be extended to the entire duration ofthe simulation (close to half a magnetic diffusion time), as the av-erage dipolar fraction f dip is equal to . .The radial magnetic field features intense localized fluxpatches of mostly normal polarity in each hemisphere, with a fewreverse flux patches paired with normal ones, mostly at high lati-tudes. Those fluid regions hosting a locally strong magnetic fieldare characterized by more quiescent flows, as can be seen in theequatorial planes of Fig. 6(c) and Fig. 6(d). In a global sense, ourreference simulation can be qualified as a strong-field dynamo sincethe ratio of the magnetic energy E m to kinetic energy E k is slightlylarger than unity, E m /E k = 1 . .Figure 7 shows a comparison between the radial component ofthe magnetic field at the CMB (left) and its representation truncatedat spherical harmonic degree (cid:96) max = 13 (right). The truncated fieldpresents several key similarities with the geomagnetic field at theCMB, such as a significant axial dipole and patches of reverse po-larity in both hemispheres. For a more quantitative assessment of the Earth-likeness of the dy-namo models, we employ the rating of compliance χ introducedby Christensen et al. (2010). This quantity is derived from four cri-teria based on the magnetic field at the CMB truncated at the degree (cid:96) = 8 (i) the relative axial dipole energy AD/NAD , which corresponds tothe ratio of the magnetic energy in the axial dipole field to that ofthe rest of the field up to degree and order eight, (ii) the equatorial symmetry
O/E corresponding to the ratio of themagnetic energy at the CMB of components that have odd valuesof ( (cid:96) + m ) for harmonic degrees between two and eight to its coun-terpart in components with (cid:96) + m even, (iii) the zonality Z/NZ , which corresponds to the ratio of the zonalto non-zonal magnetic energy for harmonic degrees two to eight atthe CMB, (iv) the flux concentration factor
F CF , defined by the variance in thesquared radial field.To evaluate these quantities for the Earth, Christensen et al. (2010)used different models based on direct measures such as gufm1 model by Jackson et al. (2000) and IGRF-11 model (from Fin-lay et al. 2010), as well as archeomagnetic and lake sediment data(model CALS7K.2 from Korte et al. 2005) and a statistical modelfor paleofield (see Quidelleur & Courtillot 1996). These models al-low to estimate the evolution of the mean value of the Gauss coef-ficients and their variances. Finally, they obtained the values givenin Table 5 for the four rating parameters. These values are used T. Tassin, T. Gastine and A. Fournier a) b)d)c)
Figure 6.
Three-dimensional renderings of a snapshot of simulation (*) (see Tab. A1). On each panel, the inner and outer spherical surfaces correspond todimensionless radii r = 0 . and r = 1 . , respectively. The large red arrow is the direction of global rotation. (a): Temperature perturbation. (b): Composition.(c): Velocity. The outer surface shows the zonal velocity u ϕ , while the inner surface, the equatorial cut and the two meridional cuts display the magnitude ofthe velocity field, | u | . (d): Magnetic field. The outer surface shows the radial field B r , while the inner surface, the equatorial cut and the two meridional cutsdisplay the magnitude of the magnetic field, | B | . All fields are dimensionless. Figure 7.
Hammer projection of the dimensionless radial magnetic field at the CMB truncated at the spherical harmonic degree (cid:96) max = 341 (a) and atspherical harmonic degree (cid:96) max = 13 (b) for the numerical simulation (*) shown in Figure 6. It corresponds to the same snapshot as in Figure 6. ouble-diffusive geodynamo models − − − − − − R m Earth-like P % T = Sc = Multipolar Dipolar Hybrid Fixed flux
Christensen et al. (2010) − − − − − − E η Earth-like0% < P % T < Pr = . Sc = − − − − − − Earth-like P % T = Pr = . Excellent Good Marginal Non compliant
Figure 8.
Compliance parameter χ as a function of the magnetic Ekman number E η and of the magnetic Rayleigh number Rm for three different setups:purely chemical convection (left), double-diffusive convection (center) and pure thermal convection (right). Red dashed lines mark the limits of the Earth-likedomain defined by Christensen et al. (2010) and the black symbol (cid:76) marks the approximate position of the Earth’s dynamo in this representation. Thesignificant size of the error bars is due to the wide range of estimates for u rms and η (see Tab. 2). The orange triangles correspond to the simulations providedby Christensen et al. (2010) with P r = 0 . or Sc = 3 . Table 5.
Time-average of the rating parameters defined by Christensen et al.(2010) for Earth and the simulation (*) (see Tab. A1).
Name Earth’s value Simulation (*)
AD/NAD
O/E
Z/NZ
F CF to determine the rating of compliance between numerical dynamomodels and the geomagnetic field χ expressed by χ = (cid:88) ψ k ln( ψ k ) − ln( ψ ⊕ k )ln (cid:20)(cid:114) Var (cid:16) ψ ⊕ k (cid:17) (cid:21) , where ψ k ∈ { AD/NAD, O/E, Z/NZ, F CF } , Var( ψ ⊕ k ) is thevariance of ψ ⊕ k and the exponent ⊕ stands for the Earth core. Theagreement between simulation and Earth is termed by Christensenet al. (2010) as excellent if χ < , as good when ≤ χ ≤ ,as marginal wen < χ ≤ and non-compliant when χ > .We adopt the same classification in the following. According to Ta-ble 5, the relative axial dipole power ( AD/NAD ) and the equato-rial symmetry (
O/E ) of the reference simulation (*) are too largein comparison with the reference geomagnetic values, which pe-nalises the overall compliance of the simulation. Nevertheless, thesimulation remains in excellent agreement with the geomagneticfield with χ = 1 . .Based on this rating of compliance χ , Christensen et al.(2010) propose a representation to classify the different numeri-cal dynamos according to the ratio of three different timescales:the rotation period τ Ω , the advection time τ adv and the magneticdiffusion time τ η . Those can be cast into two dimensionless num-bers: the magnetic Reynolds number Rm defined by τ η /τ adv andthe magnetic Ekman number E η defined by τ Ω /τ η = E/P m .Figure 8 shows our numerical simulations plotted in the param- eter space ( E η , Rm ). For comparison purposes, the simulationsby Christensen et al. (2010) with P r = 0 . or Sc = 3 (tri-angular markers) have been added to our dynamo models. Tosingle out the effect of the diffusivities, the purely chemical, thedouble-diffusive and the purely thermal simulations have been plot-ted separately. The black symbol marks the approximate position ofEarth’s dynamo in this representation (see Tab. 1 and Tab. 4). Chris-tensen et al. (2010) posited the existence, in this parameter space,of a triangular wedge (delimited by red dashed lines in Fig. 8),inside which the numerical dynamos yield Earth-like surface mag-netic fields (those from our dataset are shown in Fig. 8 with darkblue and blue disks for Excellent and Good ratings, respectively).Below the wedge, weakly super-critical dynamos feature too dipo-lar magnetic fields, on the account of the modest value of the mag-netic Reynolds number. Conversely, a significant increase of the in-put power at a given E η yields small-scale convective flows, whichpossibly lead to the breakdown of the axial dipole (multipolar sim-ulations displayed with crosses in Fig. 8). We observe in Fig. 8 thatthe upper boundary of the wedge actually depends on the value of P % T : while several purely chemical simulations are already multi-polar below Rm (cid:39) , the transition to multipolar dynamos isdelayed to Rm > for the purely thermal ones. Although allthose dynamo models that possess a good or excellent semblancewith the geomagnetic field at the CMB lie within the boundariesof the original wedge, having a pair ( E η , Rm ) that lies within thiswedge cannot be considered as a sufficient condition to producean Earth-like magnetic field, since the wedge includes a number ofsimulations with marginal or poor semblance.To further discuss the impact of P % T on the morphology ofthe magnetic field at the CMB, Fig. 9 shows χ in the parameterspace ( P tot , P % T ) for the three different Ekman numbers consid-ered here. The orange dashed lines mark the tentative boundariesbetween dipolar and multipolar simulations in term of P tot . Thedipole-multipole transition is delayed to larger input power P tot at lower Ekman numbers. Decreasing E indeed enables the explo-ration of a physical regime with lower Ro prone to sustain dipole- T. Tassin, T. Gastine and A. Fournier P % T E = − a) P tot E = × − b) (*) (**) E = − c) HybridFixed flux dipolarmultipolar
Excellent Good Marginal Non compliant
Figure 9.
Compliance parameter χ as a function of the input power P tot and of the relative thermal buoyancy power P % T for the three different Ekmannumbers considered in this study. The different colors correspond to the compliance quality defined by Christensen et al. (2010). The orange dashed lines marka tentative extrapolation of the transition between dipolar and multipolar dynamos in this parameter space. dominated dynamos (see Kutzner & Christensen 2002; Christensen& Aubert 2006). The input power required to obtain multipolardynamos is multiplied by roughly when decreasing E from − to − . For each Ekman number, the width of the dipo-lar window strongly depends on P % T since the actual input powerneeded to reach the transition is an order of magnitude lower forpure chemical convection than for pure thermal convection. Simu-lations with P % T = 40 % are found to behave similarly as purelythermal convection. We further observe that Earth-like dynamoscan be obtained for any partitioning of power injection with thebest agreement obtained close to the dipole-multipole transition.This is a consequence of the way we have sampled the parameterspace, mainly adopting one single P m value for each Ekman num-ber. Magnetic Reynolds numbers Rm ∼ O (1000) conducive toyield Earth-like fields are then attained at strong convective forc-ings. Adopting larger P m values at more moderate chemical andthermal Rayleigh numbers could hence produce Earth-like fieldsfurther away from the dipole-multipole transition (see Christensenet al. 2010).Additional diagnostics are hence required to better understandwhy the dipole-multipole transition depends so strongly on the na-ture of the convective forcing.
The physical reasons which cause the breakdown of the dipole innumerical models remain poorly known. Several previous studiessuggest that the dipole may collapse when inertia reaches a sizeablecontribution in the force balance (see Sreenivasan & Jones 2006;Christensen & Aubert 2006).Christensen & Aubert (2006) introduced the local Rossbynumber Ro L as a proxy of the ratio between inertia and Coriolisforce and found no dipole-dominated dynamos for Ro L > . (see Christensen 2010). Figure 10 shows f dip as a function of Ro L for the simulations computed in this study. The vertical linemarks the threshold value of Ro L = 0 . put forward by Chris-tensen & Aubert (2006), while the horizontal line corresponds to f dip = 0 . , the boundary between dipolar and multipolar dynamosadopted in this study. To single out the effect of partitioning theinput power between chemical and thermal forcings, the symbolshave been color coded according to P % T . Each subset of models − − Ro L . . . . . . f d i p EarthDipole-multipole transition
HybridFixed flux E = − E = × − E = − P % T Figure 10.
The dipolar fraction f dip as a function of the local Rossby num-ber Ro L . Geophysical range of f dip based on the COV-OBS.x1 model byGillet et al. (2015). Blue markers correspond to simulations with P % T =100 % , while red ones correspond to simulations with P % T = 0 % . Thehorizontal dashed orange line marks the limit between dipolar and multi-polar dynamos adopted in this study (see Sec. 2.5 for details). The verticaldashed green line marks the expected limit between dipolar and multipolardynamos according to Christensen & Aubert (2006). Vertical and horizon-tal black segments attached to each symbol represent one standard deviationabout the time-averaged values. with comparable P % T exhibits the same decrease of f dip with Ro L as reported by Christensen & Aubert (2006). However, the dipole-multipole transition occurs at lower Ro L when P % T decreases. For P % T ≥
40 % , the transition happens close to Ro L = 0 . while ithappens around Ro L = 0 . for P % T = 0 % . In addition, the dy-namo models with P % T ≥
40 % are clearly separated in two groupsof simulations with either f dip ≥ . or f dip < . , while thedipole-multipole transition is much more gradual for pure chemicalforcing. Ro L hence fails to capture the transition between dipolarand multipolar dynamos, independently of the transport propertiesof the convecting fluid.Following Christensen & Aubert (2006), Garcia et al. (2017)although envision that the increasing role of inertia would be re-sponsible for the dipole breakdown. They however define anotherparameter to characterize it. They suggest that the transition is re- ouble-diffusive geodynamo models .
65 0 .
70 0 .
75 0 .
80 0 .
85 0 .
90 0 .
95 1 . ζ . . . . . . f d i p Earth
HybridFixed flux E = − E = × − E = − P % T Figure 11.
The dipolar fraction f dip as a function of the relative portion ofequatorially-symmetric kinetic energy ζ for the simulations computedin this study. Geophysical range of f dip based on the COV-OBS.x1 modelby Gillet et al. (2015). Geophysical estimates of ζ are based on the study ofAubert (2014). Blue markers correspond to simulations with P % T = 100 % ,while red ones correspond to simulations with P % T = 0 % . The horizontaldashed orange line marks the limit between dipolar and multipolar dynamosadopted in this study (see Sec. 2.5 for details). Vertical and horizontal blacksegments attached to each symbol represent one standard deviation aboutthe time-averaged values. lated to a change in the equatorial symmetry properties of the con-vective flow. To quantify it, they introduce the ratio ζ previously de-fined in Eq. (16). Figure 11 shows f dip as a function of ζ . Geophys-ical estimates for ζ are based on the study of Aubert (2014). The de-crease of the relative equatorial symmetry ζ goes along with a grad-ual weakening of the axisymmetric dipolar field. Below ζ = 0 . ,no dipolar dynamos are obtained while conversely the models with ζ ≥ . are all dipolar. However, this parameter has little pre-dictive power to separate the dipolar solutions from the multipolarones over the intermediate range . ≤ ζ ≤ . .Another way to examine the dipole-multipole transition re-sorts to looking at the force balance governing the outer core flowdynamics (see Soderlund et al. 2012, 2015). To do so, we analyseforce balance spectra following Aubert et al. (2017) and Schwaigeret al. (2019). Each force entering the Navier-Stokes equation (3) ishence decomposed into spherical harmonic contributions. The spa-tial root mean square F RMS of a vector F reads F RMS ( t ) = (cid:113) (cid:104) F ( r , t ) (cid:105) V o \ δ , (23)where δ represents the thickness of the viscous boundary layer and V o \ δ the outer core volume that excludes those boundary layers. Byusing the decomposition in spherical harmonics, the above expres-sion can be rearranged as F ( t ) = 1 V o \ δ (cid:96) max (cid:88) (cid:96) =0 (cid:90) r o − δr i + δ (cid:96) (cid:88) m = − (cid:96) | F (cid:96)m ( r, t ) | r d r, We define the time-averaged spectrum F (cid:96) as a function of the har-monic degree (cid:96) for the force F by the relation F (cid:96) = (cid:118)(cid:117)(cid:117)(cid:116) V o \ δ (cid:90) r o − δr i + δ (cid:96) (cid:88) m = − (cid:96) | F (cid:96)m ( r, t ) | r d r . (24)Figure 12 shows the time-averaged force balance spectra for one dipolar and one multipolar dynamo with E = 10 − and P T % (cid:39) . For both panels, the spherical harmonic at whichthe poloidal kinetic energy peaks (cid:98) (cid:96) is indicated by filled markers.The left panel corresponds to the force balance of the refer-ence case (*) which is in excellent agreement with the geomag-netic field in terms of its low χ value (recall Sec.3.2). Its spectrafeature a dominant quasi-geostrophic balance between Coriolis andpressure forces up to (cid:96) ≈ accompanied by a magnetostrophicbalance at smaller scales. The difference between pressure andCoriolis forces, forming the so-called ageostrophic Coriolis force(red dashed line), is then balanced by the two buoyancy sources(black dashed line) at large scales and by Lorentz force (greenline) at small scales. This forms the quasi-geostrophic Magneto-Archimedean-Coriolis balance (QG-MAC) devised by Davidson(2013) and expected to control the outer core fluid dynamics (seeRoberts & King 2013). This hierarchy of forces is similar to the oneobserved in geodynamo models that use a codensity approach (e.g.Aubert et al. 2017; Schwaiger et al. 2019). The breakdown of buoy-ancy sources reveals a dominant contribution of chemical forcingswhich grows at small scales. The QG-MAC balance is perturbedby a sizeable inertia, which reaches almost a third of the amplitudeof Lorentz force below (cid:96) ≈ , while viscous effects are deferredto more than one order of magnitude below. Because of the strongconvective forcing, the force separation is hence not as pronouncedas in the exemplary dipolar cases by (Schwaiger et al. 2019, theirFig. 2).By increasing the input power by a factor while keeping P % T constant, we obtain a numerical model (**) (see Tab A1 andFig. 4c) where dynamo action yields a multipolar field. As com-pared to the dipole-dominated solution, the amplitude of each con-tribution is hence shifted to higher values. Most noticeable changesconcern the prominent contribution of the chemical buoyancy forthe degree (cid:96) = 1 and the ratio of inertia to Lorentz force. Theformer comes from a pronounced equatorial asymmetry of thechemical fluctuations. The development of strong equatorially-asymmetric convective motions has been observed by Landeau &Aubert (2011) and Dietrich & Wicht (2013) with a codensity ap-proach and flux boundary conditions. Below (cid:96) ≈ , inertia reachesa comparable amplitude to Lorentz force, while the smaller scalesare still controlled by magnetic effects. This differs from the mul-tipolar dynamo model described by Schwaiger et al. (2019), whereinertia was significantly larger than Lorentz force at all scales form-ing the so-called quasi-geostrophic Coriolis-Inertia-Archimedianbalance (e.g. Gillet & Jones 2006). Here the situation differs likelybecause of the larger P m , which enables a stronger magnetic field(see Menu et al. 2020). At the dominant lengthscale (cid:98) (cid:96) , the ratio F (cid:96)i /F (cid:96)L is around for the multipolar model, while it is less than . for the dipolar one. To examine whether the dipole-multipoletransition is controlled by the ratio of inertia over Lorentz forces,we hence focus on the force balance at the dominant lengthscale (cid:98) (cid:96) ,in contrast to previous studies, which analysed ratio of integratedforces (Soderlund et al. 2012, 2015; Yadav et al. 2016).Figure 13 shows the time-averaged force balance at (cid:98) (cid:96) for thesimulations with E = 10 − and
30 % ≤ P % T ≤
60 % . Thedynamics at (cid:98) (cid:96) is primarily controlled by the geostrophic balancebetween the Coriolis force and the pressure gradient. The othercontributions grow differently with P tot : viscosity and inertia in-crease continuously while Lorentz force at (cid:98) (cid:96) hardly increases be-yond P tot ≈ × . The dipole-multipole transition occurswhen inertia reaches a comparable amplitude to Lorentz force at (cid:98) (cid:96) (crosses).The relevance of this force ratio for sustaining the dipolar field T. Tassin, T. Gastine and A. Fournier (cid:96) F o r ce s ( R M S ) Ra T = . × , Ra ξ = . × (*)a) CoriolisPressure Coriolis-PressureTotal buo. Chemical buo.Thermal buo. InertiaLorentz Viscosity10 (cid:96) Ra T = . × , Ra ξ = . × (**)b) Figure 12.
Force balance spectra as a function of the spherical harmonic degree (cid:96) for a dipolar (a) and a multipolar (b) simulations with E = 10 − and P % T ≈
46 % . Solid and dashed lines represent time average of each force. The corresponding shaded regions represent one standard deviation about themean. The abscissa of the markers corresponds to the dominant lengthscale (cid:98) (cid:96) for each simulation. A circle corresponds to a dipolar simulation while a crosscorresponds to a multipolar one. Both simulations are referenced as simulations (*) and (**) in Tab. A1. P tot F o r ce s ( R M S ) (*) (**) E = − CoriolisPressureChemical buo. Thermal buo.LorentzInertia ViscosityMultipolarDipolar
Figure 13.
Time-averaged force balance spectra at the dominant lengthscale (cid:98) (cid:96) as a function the total buoyancy power P tot for numerical models with E = 10 − and
30 % < P % T <
60 % . The dynamo models (*) and (**)correspond to simulations referenced in Tab. A1. The horizontal and verti-cal segments attached to each symbol correspond to one standard deviationabout the time-averaged values. has already been put forward by Menu et al. (2020), using modelswith a purely thermal forcing and
P r = 1 . By considering tur-bulent simulations with large
P m , they show that strong Lorentzforces at large scale prevent the collapse of the dipole by inertia.As a result, they report dipole-dominated simulations with Ro L which exceeds the limit of . proposed by Christensen & Aubert (2006). Here, we quantify the contribution of inertia at the domi-nant convective lengthscale by dividing the amplitude of inertia F (cid:98) (cid:96)i by the amplitude of Coriolis F (cid:98) (cid:96)C and Lorentz F (cid:98) (cid:96)L forces.Figure 14(a) shows our simulations in the parameter spacedefined by the amplitude ratios ( F (cid:98) (cid:96)i /F (cid:98) (cid:96)L , F (cid:98) (cid:96)i /F (cid:98) (cid:96)C ). Each simula-tion is characterized by the proportion of thermal convection P % T and the nature of its thermal boundary conditions. Increasing theinput power of the dynamo leads to a growth of inertia, suchthat the strongly-driven cases all lie in the upper right quadrantof Fig. 14(a). The transition from dipolar to multipolar dynamosoccurs sharply when F (cid:98) (cid:96)i /F (cid:98) (cid:96)L exceeds . over a broad range of F (cid:98) (cid:96)i /F (cid:98) (cid:96)C ranging from . to . . This indicates that the transitionis much more sensitive to the ratio of inertia over Lorentz force thanto the ratio of inertia over Coriolis force. The transition for purelychemical simulations (dark red symbols) is reached at lower valuesof F (cid:98) (cid:96)i /F (cid:98) (cid:96)C and is more continuous than for the thermal ones (darkblue symbols). This confirms the trend already observed in Fig. 10,where f dip shows a much more gradual decreases with Ro L when P % T = 0% .Figure 14(b) shows f dip as a function of F (cid:98) (cid:96)i /F (cid:98) (cid:96)L . In con-trast with the previous criteria, the ratio F (cid:98) (cid:96)i /F (cid:98) (cid:96)L successfully cap-tures the transition between dipole-dominated and multipolar dy-namos which happens when F (cid:98) (cid:96)i /F (cid:98) (cid:96)L (cid:39) . , independently of thebuoyancy power fraction. The simulation, which singles out in theupper right quadrant of Fig. 14(b), is an exception to this cri-terion. This numerical model corresponds to P % T ≈
47 % with Ra T = 6 . × , Ra ξ = 9 . × with hybrid boundary con-ditions. Although it features F (cid:98) (cid:96)i /F (cid:98) (cid:96)L > . , its magnetic field is ontime-average dominated by an axisymmetric dipole ( f dip = 0 . ).This dynamo however strongly varies with time with several drops ouble-diffusive geodynamo models − F (cid:98) (cid:96) i / F (cid:98) (cid:96) L − − F (cid:98) (cid:96) i / F (cid:98) (cid:96) C a) − F (cid:98) (cid:96) i / F (cid:98) (cid:96) L . . . . . . f d i p b) HybridFixed fluxMultipolarDipolar 052040608095100 P % T Figure 14.
Ratio of inertia to Coriolis force at the dominant lengthscale F (cid:98) (cid:96)i /F (cid:98) (cid:96)C (left) and the dipolar fraction f dip (right) as a function of the ratio of inertiato the Lorentz force at the dominant lengthscale F (cid:98) (cid:96)i /F (cid:98) (cid:96)L . Blue markers correspond to simulations with P % T = 100 % , while red markers correspond tosimulations with P % T = 0 % . Horizontal orange dashed line marks the limit between dipolar and multipolar dynamo according to Sec. 2.5. Vertical greendashed line marks the limit between dipolar and multipolar dynamo in term of F (cid:98) (cid:96)i /F (cid:98) (cid:96)L ratio. Vertical and horizontal black segments attached to the symbolcorrespond to one standard deviation about the time-averaged values. − E k / E m . . . . . . f d i p HybridFixed flux E = − E = × − E = − Christensen et al. (2010) P % T Figure 15.
Dipolarity parameter f dip as a function of the kinetic to magnetic energy ratio E k /E m . The orange triangles correspond to those simulations byChristensen (2010) with P r (cid:54) = 1 . The horizontal dashed orange line marks the f dip = 0 . limit between dipolar and multipolar dynamos. The vertical dashedorange line corresponds to E k /E m = 0 . . Vertical and horizontal black segments attached to the symbols correspond to one standard deviation about thetime-averaged values for f dip and E k /E m , respectively. T. Tassin, T. Gastine and A. Fournier of the dipolar component below f dip = 0 . (see Fig. A1 in theappendix). Although the numerical model has been integrated formore than two magnetic diffusion times, the stability of the dipolecannot be granted for certain. Convection in the liquid outer core of the Earth is thought to bedriven by density perturbations from both thermal and chemicalorigins. In the vast majority of geodynamo models, the differencebetween the two buoyancy sources is simply ignored. In planetaryinteriors with huge Reynolds numbers, diffusion processes associ-ated with molecular diffusivities could indeed possibly be super-seded by turbulent eddy diffusion (see Braginsky & Roberts 1995).This hypothesis forms the backbone of the so-called “codensity”approach which assumes that both thermal and compositional dif-fusivities are effectively equal. This approach suppresses some dy-namical regimes intrinsic to double-diffusive convection (Radko2013).The main goal of this study is to examine the impact of double-diffusive convection on the magnetic field generation when boththermal and compositional gradients are destabilizing (the so-calledtop-heavy regime, see Takahashi 2014). To do so we have com-puted global dynamo models, varying the fraction between ther-mal and compositional buoyancy sources P % T , the Ekman number E and the vigor of the convective forcing using a Prandtl number P r = 0 . and a Schmidt number Sc = 3 . We have explored theinfluence of the thermal boundary conditions by considering twosets of boundary conditions for temperature and composition.Using a generalised eigenvalue solver, we have first investi-gated the onset of thermo-solutal convection. In agreement withprevious studies (Busse 2002; Net et al. 2012; Tr¨umper et al. 2012;Silva et al. 2019), we have shown that the incorporation of a desta-bilizing compositional gradient actually facilitates the onset of con-vection by reducing the critical thermal Rayleigh number. The criti-cal onset mode in the top-heavy regime of rotating double-diffusiveconvection is otherwise similar to classical thermal Rossby wavesobtained in purely thermal convection (Busse 1970).To quantify the Earth-likeness of the magnetic fields producedby the non-linear dynamo models, we have used the rating param-eters introduced by Christensen et al. (2010). Using geodynamomodels with a codensity approach, Christensen et al. (2010) sug-gested that the Earth-like dynamo models are located in a wedge-like shape in the 2-D parameter space constructed from the ratioof three typical timescales, namely the rotation rate, the turnovertime and the magnetic diffusion time. Here, we have shown thatthe physical parameters at which the best morphological agree-ment with the geomagnetic field is attained strongly depend on theratio of thermal and compositional input power. In particular, weobtain purely-compositional multipolar dynamo models that liewithin the wedge region of Earth-like dynamos (recall Figure 8a).This questions the relevance of the regime boundaries proposed byChristensen et al. (2010).We have then used our set of double-diffusive dynamos toexamine the transition between dipolar and multipolar dynamos.We have assessed the robustness of several criteria controlling thistransition that had been proposed in previous studies. Sreenivasan& Jones (2006) suggested that the dipole breakdown results froman increasing role played by inertia at strong convective forcings.Christensen & Aubert (2006) then introduced the local Rossbynumber Ro L as a proxy of the ratio of inertial to Coriolis forces. They suggested that Ro L ≈ . marks the boundary betweendipole-dominated and multipolar dynamos over a broad range ofcontrol parameters. Our numerical dynamo models with P % T ≥ follow a similar behaviour, while the transition between dipo-lar and multipolar dynamos occurs at lower Ro L ( ≈ . ) whenchemical forcing prevails. A breakdown of the dipole for dynamomodels with Ro L < . was already reported by Garcia et al.(2017) using P r > under the codensity hypothesis (their Fig. 1).Using non-magnetic numerical models, Garcia et al. (2017)further argued that the breakdown of the dipolar field is correlatedwith a change in the equatorial symmetry properties of the con-vective flow. They introduced the relative proportion of kinetic en-ergy contained in the equatorially-symmetric convective flow, ζ ,and suggested that multipolar dynamos would be associated with alower value of this quantity. However, our numerical dataset showsthat multipolar and dipolar dynamos coexist over a broad range of ζ ( . − . , recall Fig. 11), indicating that this ratio has littlepredictive power in separating dipolar from multipolar simulations.While neither Ro L nor ζ provide a satisfactory measure tocharacterise the dipole-multipole transition, the analysis of theforce balance governing the dynamo models has been found to be amore promising avenue to decipher the physical processes at stake(Soderlund et al. 2012, 2015). By considering a spectral decompo-sition of the different forces (e.g. Aubert et al. 2017; Schwaigeret al. 2019), we have shown that the transition between dipolarand multipolar dynamos goes along with an increase of inertia atlarge scales. The analysis of the force ratio at the dominant scaleof convection has revealed that the dipole-multipole transition ismuch more sensitive to the ratio of inertia to Lorentz force than tothe ratio of inertia to Coriolis force. The transition from dipolar tomultipolar dynamos robustly happens when the ratio of inertial tomagnetic forces at the dominant lengthscale of convection exceeds . , independently of P % T and the Ekman number. This confirmsthe results by Menu et al. (2020) who argued that a strong Lorentzforce prevents the demise of the axial dipole, delaying its break-down beyond Ro L ≈ . (their Figs. 3 and 4).Providing a geophysical estimate of the ratio of inertial tomagnetic forces at the dominant scale of convection in the Earth’score is not an easy task. Recent work by Schwaiger et al. (2020)suggests that the dominant scale of convection should be that atwhich the Lorentz force and the buoyancy force, both second-orderactors in the force balance, equilibrate. Extrapolation of this findingto Earth’s core yields a scale of approximately km, that corre-sponds to spherical harmonic degree . This is far beyond whatcan be constrained through the analysis of the geomagnetic secu-lar variation. Estimating the strength of both Lorentz and inertialforces at that scale is hence out of reach.We can however try to approximate the ratio of these twoforces by a simpler proxy, namely the ratio of the total kinetic en-ergy E k to total magnetic energy E m . To examine the validity ofthis approximation, Fig. 15 shows f dip as a function of the ratio E k /E m for our numerical simulations complemented with the co-density simulations of Christensen (2010) that have P r (cid:54) = 1 . Weobserve that the dipolar fraction f dip exhibits a variation similar tothat shown in Fig. 14(b) for the actual force ratio. The transitionbetween dipolar and multipolar dynamos is hence adequately cap-tured by the ratio E k /E m (Kutzner & Christensen 2002). All butone of the numerical dynamos of our dataset become multipolar for E k /E m > . , independently of E , P % T , and the type of thermalboundary conditions prescribed. Using the physical properties from ouble-diffusive geodynamo models Tab. 1 leads to the following estimate for the Earth’s core E k E m = µ ρ o u B ≈ − − − (cid:28) O (1) . Christensen (2010) and Wicht & Tilgner (2010) argued that convec-tion in the Earth’s core should operate in the vicinity of the transi-tion between dipolar and multipolar dynamos in order to explain thereversals of the geomagnetic field. This statement, however, postu-lates that reversals and the dipole breakdown are governed by thesame physical mechanism. The smallness of E k /E m for Earth’score indicates that it should operate far from the dipole-multipoletransition, contrary to the numerical evidence accumulated up tonow. The paleomagnetic record indicates that during a reversal oran excursion, the intensity of the field is remarkably low, whichsuggests that the strength of the geomagnetic field could decreaseby about an order of magnitude. This state of affairs admittedlybrings the ratio E k /E m closer to unity, yet without reaching it.So it seems that the occurrence of geomagnetic reversals is notdirectly related to an increase of the relative amplitude of inertia.Other mechanisms proposed to explain geomagnetic reversals relyon the interaction of a limited number of magnetic modes, whosenonlinear evolution is further subject to random fluctuations (e.g.Schmitt et al. 2001; P´etr´elis et al. 2009). In this useful conceptualframework, the large-scale dynamics of Earth’s magnetic field isgoverned by the induction equation alone. The origin of the fluc-tuations that can potentially lead to a reversal of polarity is not ex-plicited and it remains to be found, but evidence from mean fieldsmodels by Stefani et al. (2006) suggests that the likelihood of re-versals increases with the magnetic Reynolds number. In practice,these fluctuations could very well occur in the vicinity of the con-vective lengh scale, and have either an hydrodynamic, or a mag-netic, or an hydromagnetic origin, depending on the process driv-ing the instability. Shedding light on the origin of these fluctuationsconstitutes an interesting avenue for future numerical investigationsof geomagnetic reversals. We would like to thank Ulrich Christensen and Julien Aubert forsharing their data. Numerical computations were performed at S-CAPAD, IPGP and using HPC resources from GENCI-CINESand TGCC-Irene KNL (Grant 2020-A0070410095). Those com-putations required roughly . millions single core hours on In-tel Haswell CPUs. This amounts for
33 MWh , or equivalentlyto a bit more than tons of carbon dioxide emissions. We thankGabriel Hautreux, from GENCI-CINES, for providing us with adirect measure of the power consumption of MagIC. All the figureshave been generated using matplotlib (Hunter 2007) and par-aview ( ). The colorbar used inthis study were designed by Thyng et al. (2016) and Crameri (2018,2019). REFERENCES
Anufriev, A. P., Jones, C. A., & Soward, A. M., 2005. The Boussinesqand anelastic liquid approximations for convection in the Earth’s core,
Physics of the Earth and Planetary Interiors , , 163–190.Ascher, U. M., Ruuth, S. J., & Spiteri, R. J., 1997. Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, AppliedNumerical Mathematics , (2-3), 151–167.Aubert, J., 2014. Earth’s core internal dynamics 1840–2010 imaged by in-verse geodynamo modelling, Geophysical Journal International , (3),1321–1334.Aubert, J., Finlay, C. C., & Fournier, A., 2013. Bottom-up control ofgeomagnetic secular variation by the Earth’s inner core, Nature , ,219–223.Aubert, J., Gastine, T., & Fournier, A., 2017. Spherical convective dy-namos in the rapidly rotating asymptotic regime, Journal of Fluid Me-chanics , , 558–593.Badro, J., Fiquet, G., Guyot, F., Gregoryanz, E., Occelli, F., Antonangeli,D., & Astuto, M., 2007. Effect of light elements on the sound velocitiesin solid iron:Implications for the composition of Earth’s core, Earth andPlanetary Sciences Letters , , 233–238.Boscarino, S., Pareschi, L., & Russo, G., 2013. Implicit-Explicit Runge–Kutta Schemes for Hyperbolic Systems and Kinetic Equations in the Dif-fusion Limit, SIAM J. Sci. Comput. , (1), A22–A51.Bouffard, M., 2017. Double-Diffusive Thermochemical Convection in theLiquid Layers of Planetary Interiors : A First Numerical Explorationwith a Particle- in-Cell Method , Theses, Universit´e de Lyon.Boyd, J. P., 2001.
Chebyshev and Fourier Spectral Methods: Second Re-vised Edition , Courier Corporation.Braginsky, S. I. & Roberts, P. H., 1995. Equations governing convectionin earth’s core and the geodynamo,
Geophysical & Astrophysical FluidDynamics , , 1–97.Breuer, M., Manglik, A., Wicht, J., Tr¨umper, T., Harder, H., & Hansen, U.,2010. Thermochemically driven convection in a rotating spherical shell, Geophysical Journal International , , 150–162.Buffett, B., 2015. Core-Mantle Interactions, in Treatise on Geophysics ,vol. 8: Core Dynamics, pp. 345–358, Elsevier, 2nd edn.Busse, F. H., 1970. Thermal instabilities in rapidly rotating systems,
Jour-nal of Fluid Mechanics , (3), 441–460.Busse, F. H., 2002. Is low Rayleigh number convection possible in theEarth’s core?, Geophysical Research Letters , (7), 9.Butcher, J. C., 1964. On Runge-Kutta processes of high order, Journal ofthe Australian Mathematical Society , (2), 179–194.Christensen, U., Olson, P., & Glatzmaier, G. A., 1999. Numerical mod-elling of the geodynamo: A systematic parameter study, Geophys. J. Int. , (2), 393–409.Christensen, U., Aubert, J., Cardin, P., Dormy, E., Gibbons, S., Glatz-maier, G., Grote, E., Honkura, Y., Jones, C., Kono, M., Matsushima, M.,Sakuraba, A., Takahashi, F., Tilgner, A., Wicht, J., & Zhang, K., 2001.A numerical dynamo benchmark, Physics of the Earth and PlanetaryInteriors , (1-4), 25–34.Christensen, U. R., 2010. Dynamo Scaling Laws and Applications to thePlanets, Space Sci Rev , (1), 565–590.Christensen, U. R. & Aubert, J., 2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetarymagnetic fields, Geophysical Journal International , , 97–114.Christensen, U. R. & Wicht, J., 2015. Numerical Dynamo Simulations, in Treatise on Geophysics , vol. 8: Core Dynamics, pp. 245–282, Elsevier,2nd edn.Christensen, U. R., Aubert, J., & Hulot, G., 2010. Conditions for Earth-like geodynamo models,
Earth and Planetary Science Letters , , 487–496.Crameri, F., 2018. Geodynamic diagnostics, scientific visualisation andStagLab 3.0, Geosci. Model Dev. , (6), 2541–2562.Crameri, F., 2019. Scientific Colour Maps, Zenodo.Davidson, P., 2013. Scaling laws for planetary dynamos, GeophysicalJournal International , (1), 67–74.Dietrich, W. & Wicht, J., 2013. A hemispherical dynamo model: Impli- T. Tassin, T. Gastine and A. Fournier cations for the Martian crustal magnetization,
Physics of the Earth andPlanetary Interiors , , 10–21.Dormy, E., Soward, A. M., Jones, C. A., Jault, D., & Cardin, P., 2004. Theonset of thermal convection in rotating spherical shells, J. Fluid Mech. , , 43–70.Dziewonski, A. M. & Anderson, D. L., 1981. Preliminary reference earthmodel, Physics of the Earth and Planetary Interiors , , 297–356.Finlay, C. C. & Amit, H., 2011. On flow magnitude and field-flow align-ment at Earth’s core surface: Core flow magnitude and field-flow align-ment, Geophysical Journal International , (1), 175–192.Finlay, C. C., Maus, S., Beggan, C. D., Bondar, T. N., Chambodut, A.,Chernova, T. A., Chulliat, A., Golovkov, V. P., Hamilton, B., Hamoudi,M., Holme, R., Hulot, G., Kuang, W., Langlais, B., Lesur, V., Lowes,F. J., L¨uhr, H., Macmillan, S., Mandea, M., McLean, S., Manoj, C., Men-vielle, M., Michaelis, I., Olsen, N., Rauberg, J., Rother, M., Sabaka, T. J.,Tangborn, A., Tøffner-Clausen, L., Th´ebault, E., Thomson, A. W. P.,Wardinski, I., Wei, Z., & Zvereva, T. I., 2010. International Geomag-netic Reference Field: The eleventh generation, Geophys J Int , (3),1216–1230.Garcia, F., Oruba, L., & Dormy, E., 2017. Equatorial symmetry breakingand the loss of dipolarity in rapidly rotating dynamos, Geophysical &Astrophysical Fluid Dynamics , (5), 380–393.Gastine, T., 2019. Pizza: An open-source pseudo-spectral code for spher-ical quasi-geostrophic convection, Geophysical Journal International , (3), 1558–1576.Gastine, T., Wicht, J., & Aubert, J., 2016. Scaling regimes in sphericalshell rotating convection, J. Fluid Mech. , , 690–732.Gillet, N. & Jones, C. A., 2006. The quasi-geostrophic model for rapidlyrotating spherical convection outside the tangent cylinder, J. Fluid Mech. , (-1), 343.Gillet, N., Jault, D., Canet, E., & Fournier, A., 2010. Fast torsional wavesand strong magnetic field within the Earth’s core, Nature , (7294),74–77.Gillet, N., Barrois, O., & Finlay, C. C., 2015. Stochastic forecasting of thegeomagnetic field from the COV-OBS. x1 geomagnetic field model, andcandidate models for IGRF-12, Earth, Planets and Space , (1), 71.Glatzmaier, G. A., 1984. Numerical simulations of stellar convective dy-namos. I - The model and method, Journal of Computational Physics , , 461–484.Glatzmaier, G. A. & Roberts, P. H., 1995. A three-dimensional convectivedynamo solution with rotating and finitely conducting inner core andmantle, Physics of the Earth and Planetary Interiors , , 63–75.Glatzmaier, G. A. & Roberts, P. H., 1996. An anelastic evolutionarygeodynamo simulation driven by compositional and thermal convection, Physica D Nonlinear Phenomena , , 81–94.Goluskin, D., 2016. Internally Heated Convection and Rayleigh-B´enardConvection , Springer Briefs in Applied Sciences and Technology,Springer.Hirose, K., Labrosse, S., & Hernlund, J., 2013. Composition and state ofthe core,
Annual Review of Earth and Planetary Sciences , , 657–691.Hollerbach, R., 2000. A spectral solution of the magneto-convection equa-tionsin spherical geometry, International journal for numerical methodsin fluids , (7), 773–797.Hunter, J. D., 2007. Matplotlib: A 2D Graphics Environment, Comput.Sci. Eng. , (3), 90–95.Jackson, A., Jonkers, A. R. T., & Walker, M. R., 2000. Four centuriesof geomagnetic secular variation from historical records, PhilosophicalTransactions of the Royal Society of London. Series A: Mathematical,Physical and Engineering Sciences , (1768), 957–990.Jameson, A., Schmidt, W., & Turkel, E., 1981. Numerical solution of theEuler equations by finite volume methods using Runge Kutta time step-ping schemes, in , Ameri-can Institute of Aeronautics and Astronautics, Palo Alto,CA,U.S.A.Johnston, H. & Doering, C. R., 2009. Comparison of Turbulent ThermalConvection between Conditions of Constant Temperature and ConstantFlux, Phys. Rev. Lett. , (6), 064501.Kaplan, E. J., Schaeffer, N., Vidal, J., & Cardin, P., 2017. Subcritical Ther-mal Convection of Liquid Metals in a Rapidly Rotating Sphere, Phys. Rev. Lett. , (9), 094501.King, E. M., Stellmach, S., & Aurnou, J. M., 2012. Heat transfer byrapidly rotating Rayleigh–B´enard convection, J. Fluid Mech. , , 568–582.Konˆopkov´a, Z., McWilliams, S. R., G´omez-P´erez, N., & Goncharov, A. F.,2016. Direct measurements of thermal conductivity in solid iron at plan-etary core conditions, Nature , (7605), 99–101.Korte, M., Genevey, A., Constable, C. G., Frank, U., & Schnepp, E., 2005.Continuous geomagnetic field models for the past 7 millennia: 1. A newglobal data compilation, Geochemistry, Geophysics, Geosystems , (2).Kosloff, D. & Tal-Ezer, H., 1993. A Modified Chebyshev PseudospectralMethod with an O(N-1) Time Step Restriction, Journal of ComputationalPhysics , (2), 457–469.Kutzner, C. & Christensen, U., 2002. From stable dipolar towards revers-ing numerical dynamos, Physics of the Earth and Planetary Interiors , (1), 29–45.Labrosse, S., 2003. Thermal and magnetic evolution of the Earth’s core, Physics of the Earth and Planetary Interiors , (1-3), 127–143.Labrosse, S., 2015. Thermal evolution of the core with a high thermalconductivity, Physics of the Earth and Planetary Interiors , , 36–55.Landeau, M. & Aubert, J., 2011. Equatorially asymmetric convection in-ducing a hemispherical magnetic field in rotating spheres and implica-tions for the past martian dynamo, Physics of the Earth and PlanetaryInteriors , (3 – 4), 61–73.Li, Y., Fruehan, R. J., Lucas, J. A., & Belton, G. R., 2000. The chemicaldiffusivity of oxygen in liquid iron oxide and a calcium ferrite, Metalland Materi Trans B , (5), 1059–1068.Lister, J. R. & Buffett, B. A., 1995. The strength and efficiency of thermaland compositional convection in the geodynamo, Physics of the Earthand Planetary Interiors , , 17–30.Loper, D. E. & Roberts, P. H., 1981. A study of conditions at the innercore boundary of the earth, Physics of the Earth and Planetary Interiors , (4), 302–307.Manglik, A., Wicht, J., & Christensen, U. R., 2010. A dynamo modelwith double diffusive convection for Mercury’s core, Earth and Plane-tary Science Letters , , 619–628.Marti, P., Calkins, M. A., & Julien, K., 2016. A computationally efficientspectral method for modeling core dynamics, Geochemistry, Geophysics,Geosystems , (8), 3031–3053.Matsui, H., Heien, E., Aubert, J., Aurnou, J. M., Avery, M., Brown, B.,Buffett, B. A., Busse, F., Christensen, U. R., Davies, C. J., Feather-stone, N., Gastine, T., Glatzmaier, G. A., Gubbins, D., Guermond, J.-L.,Hayashi, Y.-Y., Hollerbach, R., Hwang, L. J., Jackson, A., Jones, C. A.,Jiang, W., Kellogg, L. H., Kuang, W., Landeau, M., Marti, P., Olson, P.,Ribeiro, A., Sasaki, Y., Schaeffer, N., Simitev, R. D., Sheyko, A., Silva,L., Stanley, S., Takahashi, F., Takehiro, S.-i., Wicht, J., & Willis, A. P.,2016. Performance benchmarks for a next generation numerical dynamomodel, Geochemistry, Geophysics, Geosystems , (5), 1586–1607.Menu, M. D., Petitdemange, L., & Galtier, S., 2020. Magnetic effects onfields morphologies and reversals in geodynamo simulations, Physics ofthe Earth and Planetary Interiors , p. 106542.Moll, R., Garaud, P., & Stellmach, S., 2016. A new model for mixing bydouble-diffusive convection (semi-convection). III. Thermal and compo-sitional transport through non-layered ODDC,
The Astrophysical Jour-nal , , 33.Monville, R., Vidal, J., C´ebron, D., & Schaeffer, N., 2019. Rotating con-vection in stably-stratified planetary cores, Geophysical Journal Interna-tional , (Supplement 1), S195–S218.Net, M., Garcia, F., & S´anchez, J., 2012. Numerical study of the onset ofthermosolutal convection in rotating spherical shells, Physics of Fluids , (6), 064101.Olson, P. & Christensen, U. R., 2006. Dipole moment scaling forconvection-driven planetary dynamos, Earth and Planetary Science Let-ters , (3-4), 561–571.Olson, P., Christensen, U. R., & Glatzmaier, G. A., 1999. Numerical mod-elling of the geodynamo: Mechanism of field generation and equilibra-tion, Journal of Geophysical Research: Solid Earth , (B5).Oruba, L. & Dormy, E., 2014. Transition between viscous dipolar and in- ouble-diffusive geodynamo models ertial multipolar dynamos, Geophysical Research Letters , (20), 7115–7120.Parker, E. N., 1955. Hydromagnetic dynamo models, The AstrophysicalJournal , , 293.P´etr´elis, F., Fauve, S., Dormy, E., & Valet, J.-P., 2009. Simple mechanismfor reversals of earth’s magnetic field, Physical Review Letters , (14),144503.Pozzo, M., Davies, C., Gubbins, D., & Alf`e, D., 2013. Transport prop-erties for liquid silicon-oxygen-iron mixtures at Earth’s core conditions, Physical Review B , .Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T., 1992. Numerical Recipes in Fortran 77 , vol. 1 of
Fortran Numerical Recipes ,Cambridge University Press, 2nd edn.Quidelleur, X. & Courtillot, V., 1996. On low-degree spherical harmonicmodels of paleosecular variation,
Physics of the earth and planetary in-teriors , (1-2), 55–77.Radko, T., 2013. Double-Diffusive Convection , Cambridge UniversityPress.Roberts, P. H. & King, E. M., 2013. On the genesis of the Earth’s mag-netism,
Reports on Progress in Physics , .Ruddick, B., 1983. A practical indicator of the stability of the water col-umn to double-diffusive activity, Deep Sea Research Part A. Oceano-graphic Research Papers , (10), 1105–1107.Schaeffer, N., 2013. Efficient spherical harmonic transforms aimedat pseudospectral numerical simulations, Geochemistry, Geophysics,Geosystems , (3), 751–758.Schaeffer, N., 2016. Effective scaling for the onset ofthermal convection in rotating planetary cores, Figshare:doi:10.6084/m9.figshare.4231376.v2 .Schaeffer, N., Jault, D., Nataf, H. C., & Fournier, A., 2017. Turbulent geo-dynamo simulations: A leap towards Earth’s core,
Geophysical JournalInternational , (1), 1–29.Schmitt, D., Ossendrijver, M. A. J. H., & Hoyng, P., 2001. Magnetic fieldreversals and secular variation in a bistable geodynamo model, Physicsof the Earth and Planetary Interiors , (1-4), 119–124.Schwaiger, T., Gastine, T., & Aubert, J., 2019. Force balance in numer-ical geodynamo simulations: A systematic study, Geophysical JournalInternational , (Supplement 1), S101–S114.Schwaiger, T., Gastine, T., & Aubert, J., 2020. Relating force balancesand flow length scales in geodynamo simulations, Geophys J Int .Silva, L., Mather, J. F., & Simitev, R. D., 2019. The onset of thermo-compositional convection in rotating spherical shells,
Geophysical & As-trophysical Fluid Dynamics , (4), 377–404.Simitev, R. D., 2011. Double-diffusive convection in a rotating cylindricalannulus with conical caps, Physics of the Earth and Planetary Interiors , , 183–190.Soderlund, K. M., King, E. M., & Aurnou, J. M., 2012. The influenceof magnetic fields in planetary dynamo models, Earth and PlanetaryScience Letters , , 9–20.Soderlund, K. M., Sheyko, A., King, E. M., & Aurnou, J. M., 2015. Thecompetition between Lorentz and Coriolis forces in planetary dynamos, Progress in Earth and Planetary Science , (1), 24.Spiegel, E. A., 1972. Convection in Stars: II. Special effects, Annual Re-view of Astronomy and Astrophysics , , 261.Sreenivasan, B. & Jones, C. A., 2006. The role of inertia in the evolutionof spherical dynamos, Geophysical Journal International , (2), 467–476.Sreenivasan, B. & Jones, C. A., 2011. Helicity generation and subcriticalbehaviour in rapidly rotating dynamos, Journal of Fluid Mechanics , ,5–30.Stefani, F., Gerbeth, G., G¨unther, U., & Xu, M., 2006. Why dynamosare prone to reversals, Earth and Planetary Science Letters , (3-4),828–840.Takahashi, F., 2014. Double diffusive convection in the Earth’s core andthe morphology of the geomagnetic field, Physics of the Earth and Plan-etary Interiors , , 83–87.Takahashi, F., Shimizu, H., & Tsunakawa, H., 2019. Mercury’s anomalousmagnetic field caused by a symmetry-breaking self-regulating dynamo, Nature Communications , , 208.Thyng, K., Greene, C., Hetland, R., Zimmerle, H., & DiMarco, S., 2016.True Colors of Oceanography: Guidelines for Effective and AccurateColormap Selection, Oceanog , (3), 9–13.Tilgner, A., 1999. Spectral methods for the simulation of incompressibleflows in spherical shells, International Journal for Numerical Methodsin Fluids , (6), 713–724.Tr¨umper, T., Breuer, M., & Hansen, U., 2012. Numerical study on double-diffusive convection in the Earth’s core, Physics of the Earth and Plane-tary Interiors , , 55–63.Valet, J.-P. & Fournier, A., 2016. Deciphering records of geomagneticreversals, Rev. Geophys. , (2), 410–446.Vidal, J. & Schaeffer, N., 2015. Quasi-geostrophic modes in the Earth’sfluid core with an outer stably stratified layer, Geophys. J. Int. , (3),2182–2193.Wicht, J., 2002. Inner-core conductivity in numerical dynamo simulations, Physics of the Earth and Planetary Interiors , (4), 281–302.Wicht, J. & Sanchez, S., 2019. Advances in geodynamo modelling, Geo-physical & Astrophysical Fluid Dynamics , (1-2), 2–50.Wicht, J. & Tilgner, A., 2010. Theory and Modeling of Planetary Dy-namos, Space Sci Rev , (1), 501–542.Yadav, R. K., Gastine, T., & Christensen, U. R., 2013. Scaling laws inspherical shell dynamos with free-slip boundaries, Icarus , (1), 185–193.Yadav, R. K., Gastine, T., Christensen, U. R., Wolk, S. J., & Poppenhaeger,K., 2016. Approaching a realistic force balance in geodynamo simula-tions, PNAS , (43), 12065–12070.Zhang, Y., Hou, M., Liu, G., Zhang, C., Prakapenka, V. B., Greenberg, E.,Fei, Y., Cohen, R. E., & Lin, J.-F., 2020. Reconciliation of Experimentsand Theory on Transport Properties of Iron and the Geodynamo, Phys.Rev. Lett. , (7), 078501. T. Tassin, T. Gastine and A. Fournier
APPENDIX A: TIME SCHEME
We provide in this appendix the matrices a I and a E and the vectors b I , b E , c E and c I of the two IMEX Runge-Kutta schemes that weresorted to for this study. These vectors and matrices can conveniently be represented using Butcher tables (Butcher 1964). BPR353: (Boscarino et al. 2013)Implicit component c I A I b I = 0 01 1 / / / / − / /
21 1 / /
21 1 / / − / / / / − / / Explicit component c E A E b E = 0 01 1 02 / / / / / / / / / PC2: (Jameson et al. 1981)Implicit component c I A I b I = 0 01 1 / /
21 1 / /
21 1 / / / / Explicit component c E A E b E = 0 01 1 01 1 / / / / / / APPENDIX B: NUMERICAL SIMULATIONS
Table A1: Control parameters and simulation diagnostics for the numerical simulations computed for this study. Simulations computed using the finitedifference method in radius are marked with a superscript f (the others were computed using the Chebyshev collocation method in radius). Simulations withhybrid boundary conditions are marked by an H in the first column. Simulations are sorted by growing Ekman number and then by growing magnetic Reynoldsnumber. The averaging and running times t avg and t run are expressed in units of magnetic diffusion time τ η . Ra T Ra ξ ( n r , (cid:96) max ) P m α t scheme t avg t run Rm Ro L Λ f ohm χ ζ f dip P % T Nu Sh ( × )( × ) ( × − ) E = 1 × − , P r = 0 . , Sc = 3 H . , . . BPR353 .
80 1 .
80 176 1 . . .
66 20 . .
80 0 .
91 32 1 . . . , . . BPR353 .
79 3 .
09 196 0 . . .
79 12 . .
92 0 .
94 100 2 . .
00 100 (193 , . . BPR353 .
29 0 .
29 444 4 . . .
51 3 . .
76 0 .
69 0 1 . . H . , . . BPR353 .
20 0 .
37 514 4 . . .
74 1 . .
69 0 .
80 40 4 . .
10 270 (320 , f . None BPR353 .
47 0 .
95 685 7 . . .
54 7 . .
73 0 .
33 0 1 . .
20 400 (380 , f . None BPR353 .
16 0 .
20 880 8 . . .
49 14 . .
74 0 .
24 0 1 . .
046 0 (109 , . . BPR353 .
38 0 .
90 763 5 . . .
75 2 . .
71 0 .
77 100 7 . . H . , . . BPR353 .
17 0 .
28 760 6 . . .
70 2 . .
70 0 .
75 44 6 . . H . , f . None BPR353 .
46 0 .
67 580 10 . . .
65 1 . .
71 0 .
77 46 8 . . (*)
100 0 (129 , . . BPR353 .
26 1 .
29 1213 10 21 . .
67 0 . .
72 0 .
71 100 10 . . H . , f . None BPR353 .
54 1 .
88 843 15 . . .
60 1 . .
73 0 .
69 47 11 . . (x)
300 0 (129 , . . BPR353 .
30 1 .
22 1235 17 . . .
58 19 . .
65 0 .
23 100 16 . . H
13 1900 (408 , f . None PC2 .
16 0 .
24 1252 19 . . .
57 21 . .
68 0 .
20 47 14 . . (**)
600 0 (129 , . . BPR353 .
22 0 .
43 1640 21 . . .
61 11 . .
65 0 .
32 100 20 . . E = 3 × − , P r = 0 . , Sc = 30 0 . , . . CNAB2 .
37 1 .
55 113 1 . . .
30 52 . .
89 0 .
86 0 1 . . H .
015 0 . , . . CNAB2 .
95 2 .
33 154 1 . . .
47 24 . .
93 0 .
90 44 1 . . H . . , . . CNAB2 .
15 4 .
39 98 2 0 . .
31 36 . .
88 0 .
88 21 1 . . H . . , . . BPR353 .
15 3 .
11 178 1 . . .
48 10 . .
85 0 .
86 22 1 . . . . , . . BPR353 .
35 1 .
79 193 1 . . .
51 13 . .
83 0 .
89 19 1 . .
20 2 (65 , . . CNAB2 .
81 2 .
26 204 2 . . .
43 8 . .
79 0 .
86 0 1 . . . , . . CNAB2 .
62 1 .
71 255 1 6 . .
68 8 . .
96 0 .
88 100 1 . .
00 4 (97 , . . CNAB2 .
24 1 .
42 288 3 . . .
38 2 . .
78 0 .
60 0 1 . . H . , . . CNAB2 .
04 1 .
32 307 2 . . .
54 3 . .
78 0 .
81 20 2 . .
51 0 (49 , . . CNAB2 .
96 6 .
41 180 1 . . .
69 10 . .
92 0 .
93 100 2 . .
01 0 (49 , . . CNAB2 .
79 1 .
18 281 1 . . .
81 3 . .
81 0 .
78 100 2 . . ouble-diffusive geodynamo models Table A1: Control parameters and simulation diagnostics for the 79 numerical simulations computed for this study (continued). Ra T Ra ξ ( n r , (cid:96) max ) P m α t scheme t avg t run Rm Ro L Λ f ohm χ ζ f dip P % T Nu Sh ( × )( × ) ( × − )0 8 (97 , . . CNAB2 .
91 5 .
45 385 4 . .
41 2 . .
77 0 .
67 0 1 . . . . , . . CNAB2 .
30 1 .
77 164 2 . . .
79 5 . .
77 0 .
80 86 2 . . H . . , . . CNAB2 .
21 1 .
21 373 3 . . .
62 2 . .
76 0 .
81 35 2 . . . . , . . CNAB2 .
61 1 .
89 317 1 . . .
79 5 . .
75 0 .
65 87 2 . . .
05 5 . , . . CNAB2 .
63 1 .
81 342 3 12 . .
69 2 . .
72 0 .
80 39 2 . .
91 11 (97 , . . CNAB2 .
10 1 .
30 297 5 . . .
51 4 . .
78 0 .
82 39 3 . .
91 11 (145 , . . CNAB2 .
08 1 .
23 541 4 . . .
59 1 . .
75 0 .
76 39 3 . .
40 20 (97 , . . CNAB2 .
04 1 .
06 581 6 . . .
43 5 . .
76 0 .
44 0 1 . . . , . . CNAB2 .
31 1 .
41 328 5 12 . .
71 2 . .
74 0 .
84 100 4 . . . , . . CNAB2 .
88 1 .
02 613 4 . . .
72 3 . .
73 0 .
70 100 4 . . . , . . CNAB2 .
83 1 .
08 738 7 . . .
54 1 . .
76 0 .
65 40 4 . .
35 0 (97 , . . CNAB2 .
89 1 .
06 726 5 . . .
69 2 . .
72 0 .
71 100 4 . . . , . . CNAB2 .
07 1 .
07 857 8 . . .
53 1 . .
76 0 .
63 41 4 . . . , . . BPR353 .
02 1 .
10 837 8 . .
57 1 . .
73 0 .
70 40 4 . . H . , . . BPR353 .
95 1 .
09 755 6 . . .
68 4 . .
69 0 .
63 84 4 . .
33 35 (129 , . . CNAB2 .
13 1 .
13 976 9 . . .
52 1 . .
76 0 .
59 42 5 . . H . , . . BPR353 .
43 0 .
51 949 9 . . .
57 1 . .
73 0 .
67 48 4 . .
20 100 (193 , . . BPR353 .
13 0 .
18 1378 12 . . .
40 26 . .
68 0 .
16 0 1 . .
512 0 (129 , . . CNAB2 .
94 1 .
02 1237 11 37 . .
59 2 . .
73 0 .
66 100 6 . .
07 70 (193 , . . CNAB2 .
87 1 .
22 1579 13 . .
44 27 . .
69 0 .
15 47 7 . .
19 100 (193 , . . PC2 .
10 0 .
30 1765 14 . .
49 12 . .
72 0 .
27 45 8 . . H . , . . PC2 .
12 0 .
16 1790 17 . . .
47 26 . .
71 0 .
15 51 7 . .
640 0 (97 , . . CNAB2 .
34 0 .
38 2521 19 . .
50 22 . .
71 0 .
17 100 11 . .
068 0 (129 , . . BPR353 .
59 0 .
59 1671 23 . .
53 7 . .
71 0 .
26 100 13 . . E = 1 × − , P r = 0 . , Sc = 30 .
013 0 .
045 (49 ,
85) 5 . None CNAB2 .
09 1 .
78 175 1 . . .
32 10 . .
98 0 .
82 40 1 . .
80 0 . , . None CNAB2 .
17 1 .
31 161 1 . . .
47 4 . .
81 0 .
88 0 1 . . .
02 0 .
09 (49 ,
85) 5 . None CNAB2 .
54 7 .
73 180 1 . . .
62 4 . .
86 0 .
86 36 1 . . . , . None CNAB2 .
86 2 .
16 363 1 . . .
67 11 . .
90 0 .
63 100 1 . .
00 0 . ,
85) 5 . None CNAB2 .
65 1 .
99 378 3 . .
32 1 . .
79 0 .
58 0 1 . . .
047 0 .
33 (65 ,
85) 5 . None BPR353 .
50 3 .
64 345 2 . . .
58 2 . .
78 0 .
75 37 1 . . .
08 0 .
19 (65 , . None CNAB2 .
10 1 .
18 334 1 . .
69 3 . .
81 0 .
71 67 1 . .
00 0 . ,
85) 5 . None CNAB2 .
02 2 .
43 438 4 . . .
31 2 . .
78 0 .
52 0 1 . . .
07 0 . , . None BPR353 .
35 1 .
70 438 3 . . .
56 3 . .
78 0 .
71 40 1 . .
10 1 (65 , . None CNAB2 .
06 1 .
33 522 5 . . .
31 8 . .
77 0 .
35 0 1 . . .
14 0 .
34 (65 , . None CNAB2 .
08 1 .
09 482 3 . . .
64 3 . .
78 0 .
68 70 2 . . .
19 0 (65 , . None CNAB2 .
14 1 .
36 469 2 . . .
69 5 . .
80 0 .
64 100 2 . . .
12 0 .
52 (81 , . None BPR353 .
89 2 .
24 516 3 . . .
58 3 . .
78 0 .
69 56 2 . . .
027 1 . , . None CNAB2 .
11 2 .
44 580 5 . . .
34 6 . .
78 0 .
36 10 1 . . .
027 1 . , . None CNAB2 .
14 1 .
30 570 5 . . .
37 2 . .
79 0 .
54 11 1 . . . .
69 (65 , . None CNAB2 .
12 1 .
45 533 4 . . .
55 3 . .
78 0 .
69 44 2 . . .
049 1 (65 , . None CNAB2 .
05 1 .
16 573 5 . . .
41 2 . .
79 0 .
60 20 1 . . .
049 1 (65 , . None BPR353 .
13 1 .
38 566 5 . . .
42 2 . .
78 0 .
61 20 1 . .
10 1 . , . None CNAB2 .
32 1 .
61 615 6 . .
31 16 . .
75 0 .
23 0 1 . . .
14 1 (65 , . None BPR353 .
33 2 .
10 412 5 . . .
50 2 . .
79 0 .
74 45 2 . . .
14 1 (81 , . None BPR353 .
37 1 .
72 655 5 . . .
51 3 . .
77 0 .
66 46 2 . . .
14 1 (81 , . . BPR353 .
42 2 .
36 896 5 . . .
52 4 . .
77 0 .
62 46 2 . .
40 2 (97 , . . BPR353 .
88 2 .
98 750 7 . . .
32 39 . .
72 0 .
10 0 1 . . . , . None CNAB2 .
16 1 .
24 616 3 . . .
66 4 . .
76 0 .
62 100 2 . . .
21 1 . , . . BPR353 .
32 2 .
74 813 6 . .
51 3 . .
77 0 .
63 49 2 . . . , . . BPR353 .
28 2 .
13 1144 9 . . .
44 6 . .
79 0 .
44 43 3 . . . , . . BPR353 .
34 2 .
08 1244 9 . . .
54 6 . .
74 0 .
56 100 3 . . . , . . BPR353 .
29 1 .
85 1576 12 . . .
38 42 . .
72 0 .
09 35 4 . .
45 0 (73 , . . BPR353 .
33 1 .
49 3363 23 . .
39 17 . .
75 0 .
25 100 7 . .
020 0 (129 , . . CNAB2 .
12 0 .
39 6160 45 . .
36 21 . .
76 0 .
21 100 10 . . T. Tassin, T. Gastine and A. Fournier .
00 0 .
25 0 .
50 0 .
75 1 .
00 1 .
25 1 .
50 1 . τ η )0 . . . . . f d i p . . . . . . . . E m / E k Figure A1.
Time evolution of the dipolar fraction f dip (green) and of the magnetic to kinetic energy ratio E m /E k (blue) for the simulation ( x ) of Table A1.The orange dashed line corresponds to the boundary between dipole-dominated and multipolar dynamos ( f dip = 0 . ). Time is scaled by the magneticdiffusion time. APPENDIX C: SIMULATION X In this appendix we provide in Figure A1 the detailed time evolution of the dipolar fraction f dip and of the magnetic to kinetic energy ratio E m /E k of the anomalous simulation that appears for instance in the top-right quadrant of Fig. 14(b) and Fig. 15. This simulation is markedwith a (x)(x)