A regime diagram for the slurry F-layer at the base of Earth's outer core
AA regime diagram for the slurry F-layer at the base ofEarth’s outer core
Jenny Wong a,b, ∗ , Christopher J. Davies c , Christopher A. Jones d a Universit´e de Paris, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France b EPSRC Centre for Doctoral Training in Fluid Dynamics, University of Leeds, Leeds, LS29JT, UK c School of Earth and Environment, University of Leeds, Leeds, LS2 9JT, UK d School of Mathematics, University of Leeds, Leeds, LS2 9JT, UK
Abstract
Seismic observations of a slowdown in P wave velocity at the base of Earth’souter core suggest the presence of a stably-stratified region known as the F-layer. This raises an important question: how can light elements that drive thegeodynamo pass through the stably-stratified layer without disturbing it? Weconsider the F-layer as a slurry containing solid particles dispersed within theliquid iron alloy that snow under gravity towards the inner core. We present aregime diagram showing how the dynamics of the slurry F-layer change uponvarying the key parameters: P´eclet number (
P e ), the ratio between advectionand chemical diffusion; Stefan number ( St ), the ratio between sensible andlatent heat; and Lewis number ( Le ), the ratio between thermal and chemicaldiffusivity. We obtain four regimes corresponding to stable, partially stable,unstable and no slurries. No slurry is found when the heat flow at the base ofthe layer exceeds the heat flow at the top, while a stably-stratified slurry ariseswhen advection overcomes thermal diffusion ( P e (cid:38) Le ) that exists over a widerange of parameters relevant to the Earth’s core. Our results estimate thata stably-stratified F-layer gives a maximum inner-core boundary (ICB) body ∗ Corresponding author
Email address: [email protected] (Jenny Wong)
Preprint submitted to EPSL December 15, 2020 a r X i v : . [ phy s i c s . g e o - ph ] D ec ave density jump of ∆ ρ bod ≤
534 kgm − which is compatible with the lowerend of the seismic observations where 280 ≤ ∆ ρ bod ≤ ,
100 kgm − is reportedin the literature. With high thermal conductivity the model predicts an innercore age between 0 . . Keywords: slurry, iron snow, inner core, crystallisation, F-layer
1. Introduction
Seismic, geomagnetic and dynamical explanations for and against the pres-ence of stably-stratified layers in the Earth’s liquid core is an active researchtopic of significant geophysical interest. Confirmation of their existence wouldwarrant a shift away from the usual approximation that the core is adiabaticallystratified with thin boundary layers. Stably-stratified layers should exhibit dy-namics that distinguish them from the turbulent bulk of the core; elucidatingthis behaviour will help to uncover their signature in seismic and geomagneticobservations [1]. In this paper we focus on the F-layer at the base of the liq-uid core to further understand the conditions where stable stratification can besustained.There is a strong consensus that a slowdown in the P wave velocity comparedwith the Preliminary Reference Earth Model (PREM) [2] is observed at thebase of the Earth’s outer core [3, 4, 5]. PREM follows the Adams-Williamsonequation that assumes the outer core is adiabatically stratified and homogeneousthroughout, therefore a P wave velocity lower than PREM is attributed to ananomalously higher density structure than expected. This departure in densityaway from neutral stability means that the seismic observations suggest a stably-stratified layer exists that cannot be explained by adiabatic compression alone.2stimates of the layer thickness vary from 150 km [3] to 400 km [4] thick, whichgreatly exceeds the thermal diffusion length scale, hence the layer cannot besimply explained by a thermal boundary layer alone and another mechanism isrequired to maintain stratification [6].The F-layer is intimately linked to the geodynamo process that generatesEarth’s magnetic field. Motion of the liquid iron core is powered by heat ex-tracted at the core-mantle boundary (CMB), which leads to freezing of the solidinner core from the centre of the Earth because the melting curve is steeper thanthe core adiabat. Freezing releases latent heat and light elements into the liquid;light elements drive compositional convection in the core and is thought to bethe main power source for the dynamo at the present day [7]. The key issueis buoyant light material excluded from the inner core must pass through theF-layer and into the overlying core while preserving stable stratification.Previous geodynamic models try to explain the F-layer by inner core transla-tion [8], a thermochemical layer [9], or a slurry (iron snow) layer [10, 11]. Innercore translation supposes melting occurs on the eastern hemisphere of the in-ner core while freezing occurs on the western hemisphere to form a dense layerabove the inner-core boundary (ICB), though recent upward revisions of thethermal conductivity of iron [12, 13] together with the likely presence of compo-sitionally stabilising conditions suggest that the instability cannot arise in thepresent-day [14, 15]. [9] proposed a thermochemical model with the F-layer onthe liquidus that succeeded in producing a stably-stratified layer, however, themodel effectively imposed a stable composition without a physical mechanismexplaining why. Wong et al. [11] explains this mechanism in a self-consistentway by proposing a slurry layer that captures the physics of how light elementand solid is transported. The authors used a simple box model to demonstratethat a slurry layer produced stably-stratified layers that match the seismic ob-3 CB ICB speedfixed heatflux fixed heatfluxfixed oxygenconcentration
CSB zero solidfluxsolid flux ∝ ICB speedliquidustemperature n e t i n w a r d t r a n s p o r t o f h e a vy i r o n n e t o u t w a r d t r a n s p o r t o f li g h t e l e m e n t fixed layer thickness Figure 1: A sketch of the slurry layer at the base of Earth’s outer core, with the boundary con-ditions imposed at the inner core boundary (ICB) and core slurry boundary (CSB) describedin Section 2. servations of density and layer thickness. In this work we focus on the slurryscenario and build upon Wong et al. [11], hereafter referred to as W18.In the slurry, pure solid iron particles crystallise throughout the entire F-layer while light elements remain in the liquid (see Figure 1 for a sketch). Heavygrains of iron fall under gravity to accumulate at the base of the layer, therebyproducing a net inward transport of dense solid and a net outward transportof light elements to give an overall stable density stratification. W18 simplifiedthe full slurry theory by Loper and Roberts [10, 16] and developed a reduced,thermodynamically self-consistent framework that accounts for solid and liquidphase, in addition to the influence of pressure, temperature and composition.The fraction of solid in a slurry is assumed to be small, and so cannot transmitshear waves to create an impedance contrast at the top of the F-layer that wouldotherwise have been seismically observable.The principal assumptions of the W18 slurry model are (1) the fast-meltinglimit and (2) a binary alloy. (1) supposes that an infinitesimal material volumecontains either solid or liquid phase exclusively so that the system is in phase4quilibrium. As a consequence, minimising the Gibbs free energy constrains theslurry temperature to the liquidus. (2) presumes an Fe-O composition sinceoxygen is able to explain the core density deficit [17].
Ab initio calculationsshow that oxygen almost entirely partitions into the liquid phase upon freezing[18] and cannot occur with silicon (another potential light element candidate) inlarge quantities [19]. The solid produced by the slurry is thereby assumed to becomposed of pure iron, which avoids a complicated particle history dependencethat would be difficult to model, where the growth of solid grains at differentpressure-temperature conditions affects the composition of light elements in thesolid.We generalise the W18 model by moving from a Cartesian to spherical ge-ometry and update the boundary conditions to a more geophysically realisticsetup. We suppose that the inner core is isothermal and not convecting due tothe upward revisions of the thermal conductivity [20, 21] of iron alloys at coreconditions, which means that the present-day inner core is unlikely to convect.We consider a thin compacting layer on the order of kilometres thick on thesolid side of the inner core [22], parameterised by interfacial freezing at the ICBwhich generates a flux of latent heat into the slurry. Accordingly, the ICB ad-vances at a rate composed of the interfacial freezing speed and the snow speedgiven by the accumulation of solid iron particles from the slurry. Our modelself-consistently determines the ICB speed and therefore allows the inner coregrowth rate to be independently determined, whereas this was fixed in W18 andyielded unrealistically high ICB heat flows. Following the fast-melting limit, weset the temperature at the top of the layer to the liquidus temperature for thebulk core composition obtained from the literature, so that the ICB temperatureis free to be self-determined by the slurry whereas this was fixed in W18.In this paper, we perform a systematic parameter search to elucidate the5onditions that promote stable stratification of the F-layer. The generalisedmodel of W18 and the derivation of the dimensionless system with its associateddimensionless control parameters is given in Section 2. In Section 3, we presentand discuss example solutions of the temperature, oxygen concentration, solidflux and density given by the slurry model with different layer thicknesses. Weproceed to test the sensitivity of the model to the boundary conditions on thetemperature and composition at the top of the layer, before mapping a regimediagram of the slurry based on the main dimensionless control parameters. Thisvastly improves upon the computed solution space of W18 and uncovers thepossible regimes occupied by the slurry. We demonstrate the predictive powerof our model for seismic observations and inner core age, before concluding witha summary of our findings in Section 4.
2. Model and Methods
Detailed development of the slurry model is given W18 so here we presentthe key details. We start with the general equations for the conservation ofoxygen, temperature and the liquidus constraint from W18 (equations (4),(6)and (13)), given by ρ sl D ξ D t = −∇ · i , (1) ρ sl c p D T D t = ∇ · ( k ∇ T + L j ) + ρ sl L D φ D t , (2) ∇ T = T ∆ V s,lF e L ∇ p − T ξ l (cid:0) ∂µ/∂ξ l (cid:1) L ∇ ξ l , (3)where p is the pressure, T is the temperature, ξ is the oxygen concentration, ξ l is the oxygen concentration in the liquid phase, φ is the solid fraction, ρ sl isthe reference density of the slurry taken to be the PREM value at the CSB, i
6s the oxygen flux vector, j is the solid flux vector, c p is the specific heat ca-pacity, k is the thermal conductivity, L is the latent heat of fusion, ∆ V s,lF e isthe change in specific volume between liquid iron and solid iron, and ∂µ/∂ξ l is the thermodynamic derivative of the chemical potential, µ , with respect to ξ l . Throughout this paper superscripts s , l and sl denote solid, liquid or slurryphases respectively, and subscripts F e and O denote the iron and oxygen com-ponents, respectively. For further reference, a table of symbols and values isprovided in Appendix A, Table A.2.We adopt a spherical, one-dimensional geometry where we assume a global F-layer in a non-convective steady state that depends on radius only and excludeslateral variations. The layer thickness, d , is fixed so that the slurry is definedover the interval [ r i , r sl ], where r i denotes the ICB radius and r sl ≡ r i + d isthe core-slurry boundary (CSB). We are interested in a timescale that is longcompared with the timescale of freezing and comparable to the evolution ofthe inner core [23], so we assume a Boussinesq slurry with a reference state inhydrostatic equilibrium, and resolve the slurry layer in a frame moving at therate of inner core growth, v ˆr , with v the ICB speed and ˆr the unit vector pointingoutwards and normal to the inner core surface. Mass conservation ∇ · v = 0implies that v ( r ) = v i ( r i /r ) , where the changes in ICB speed due to radius arenegligible and we therefore assume that v is constant. This speed is composedof two parts, v = v s + v f , where v s is the snow speed, in which all solid particlesfrom the slurry are assumed to accumulate at the base of the layer; and v f is thefreezing speed, which represents the growth of the inner core due to compactionon the solid side of the ICB, and assumed to occur quickly compared with theslurry timescale. The material derivative is hence D / D t −→ − v d / d r . As inW18, the fraction of solid, φ , in the slurry is small so that φ (cid:28)
1, and weassume that their variations, d φ , are negligible (see discussion in Appendix B).7e apply the above assumptions to equations (1 – 3), to obtain, − vρ sl d ξ l d r = 1 r dd r (cid:32) r ρ sl ∆ V s,lF e,O D O RT a O exp (cid:32) F (cid:0) r − r i (cid:1) d (cid:33) d p d r (cid:33) + ξ l d j d r + j d ξ l d r + 2 r ξ l j, (4) − vρ sl c p d T d r = k d T d r + 2 r k d T d r + L d j d r + 2 r Lj, (5)d T d r = − T ∆ V s,lF e L gρ − RT a O L d ξ l d r . (6)On the LHS of (4), we have the advection of oxygen, and on the RHS, thefirst term is the effect of barodiffusion and the last three terms describe thephysical displacement of oxygen as solid iron particles sediment. On the LHSof (5), we have the advection of heat, and on the RHS, the first two termscorrespond to thermal diffusion and the last two terms correspond to the latentheat release associated with phase change. The liquidus constraint (6) describesthe change in the liquidus temperature as a consequence of changes in pressure,given by the first term on the RHS, and changes in composition, given by thesecond term on the RHS.We have applied ideal solution theory [24] so that the thermodynamic deriva-tive of the chemical potential with respect to ξ l can be expressed as ξ l ∂µ/∂ξ l = RT /a O , where R is the gas constant and a O is the atomic weight of oxy-gen. We expect penetrative convection to occur at the top of the layer becauseof the velocity difference between the non-convective slurry and the overlyingconvective outer core. We have thus invoked a turbulent mixing sublayer belowthe CSB that promotes the transport of oxygen out of the layer and also allowsus to impose a vanishing solid flux boundary condition at the CSB (see equation(10) below). This effect is modelled by modifying the self-diffusion coefficient8f oxygen to take the exponential form¯ D = D O exp (cid:32) F (cid:0) r − r i (cid:1) d (cid:33) , (7)where F is the mixing parameter that appears in the first term on the RHS of(4), and D O is the self-diffusion coefficient of oxygen taken from the literature[25]. We have three equations (4),(5), and (6) for three output variables
T, ξ l , j ,and two output parameters (eigenvalues) F and v . The system is fourth order,therefore we require six boundary conditions to determine a unique solution.We suppose that the oxygen concentration is constant at the CSB and equalto the uniform value in the bulk of the core, ξ sl , since the liquid core is vigor-ously convecting and its composition does not change much over the timescalesconsidered [26], so that ξ l ( r sl ) = ξ sl . (8)At the ICB the solid flux is proportional to the ICB speed and at the CSB thesolid flux vanishes, so that j ( r i ) = − ρ s − v, (9) j ( r sl ) = 0 , (10)where ρ s − denotes the density on the solid side of the ICB.We fix the heat flux per unit area at the CSB to a constant value, q sl , andthis is a parameter to be varied since there are no independent estimates. By9ourier’s law, the boundary condition on the CSB temperature gradient isd T d r (cid:12)(cid:12)(cid:12)(cid:12) r = r sl = − q sl k . (11)In the following two boundary conditions, we digress from the conditions im-posed in W18. First, compaction below the F-layer releases latent heat atthe interface while there is no specific heat contribution from an isothermal IC,whereas W18 assumed specific heat loss from the inner core. Second, we assumethat the temperature at the top of the layer is coincident with the liquidus tem-perature of the bulk composition, T sl ( r sl ), known from estimates given in theliterature [25, 27], so thatd T d r (cid:12)(cid:12)(cid:12)(cid:12) r = r i = − ρ s − v f Lk = − q s k , (12) T ( r sl ) = T sl ( r sl ) . (13)where q s is the ICB heat flux per unit area. We derive the dimensionless equations using the following scalings, r = r sl ˆ r, T = q sl ρ sl c p v f ˆ T ,ξ l = ξ sl ˆ ξ, j = ρ s − v f ˆ ,g = g sl ˆ g, ρ = ρ sl ˆ ρ,v = v f ˆ v (14)where the hat symbol denotes a dimensionless quantity that depends on radius.The slurry layer is defined over the dimensionless interval [ r i /r sl , − ˆ v d ˆ ξ dˆ r = − r ddˆ r (cid:32) Li p R ρ Li ξ P eStR v ˆ g ˆ ρ ˆ r ˆ T exp (cid:34) F (cid:0) r sl ˆ r − r i (cid:1) d (cid:35)(cid:33) + ˆ ξ dˆ dˆ r + ˆ d ˆ ξ dˆ r + 2ˆ r ˆ ξ ˆ , (15) − ˆ v d ˆ T dˆ r = LeP e (cid:32) d ˆ T dˆ r + 2ˆ r d ˆ T dˆ r (cid:33) + 1 St (cid:18) dˆ dˆ r + 2ˆ r ˆ (cid:19) , (16)d ˆ T dˆ r = − Li p ˆ g ˆ ρ ˆ T − Li ξ StR ρ ˆ T d ˆ ξ dˆ r . (17)where the dimensionless numbers are defined as R ρ = ρ sl ρ s − , R v = ∆ V s,lF e ∆ V s,lF e,O , Li p ≡ ∆ V s,lF e g sl ρ sl r sl L , Li ξ ≡ Rξ sl a O c p ,P e ≡ v f r sl D O , St ≡ q sl ρ s − v f L , Le ≡ kρ sl c p D O . (18) R ρ is the ratio between the reference density and the density on the solid side ofthe ICB, and R v is the ratio between the change in specific volumes upon phasechange of pure iron and the iron alloy. The dimensionless numbers Li p and Li ξ arise from the pressure and compositional parts of the liquidus constraint(17), respectively. P e is the P´eclet number that measures the ratio betweenadvection and chemical diffusion. St is the Stefan number and gives the ratiobetween sensible and latent heat. Le is the Lewis number that describes theratio between thermal and chemical diffusivity.Inserting the scalings (14) into (8)–(13) yields the dimensionless boundary11onditions ˆ T (1) = T sl c p R ρ StL , (19)d ˆ T dˆ r (cid:12)(cid:12)(cid:12)(cid:12) ˆ r = rirsl = − P eStLe , (20)d ˆ T dˆ r (cid:12)(cid:12)(cid:12)(cid:12) ˆ r =1 = − P eLe , (21)ˆ ξ (1) = 1 , (22)ˆ (cid:18) r i r sl (cid:19) = − ˆ v, (23)ˆ (1) = 0 . (24) We solve the dimensionless system numerically using solve bvp , which isa boundary value problem solver included in scipy – an open source Pythonlibrary [28]. We rewrite the problem as a system of first-order ODEs, and varythe model parameters
P e , St , Le , Li p , Li ξ , R ρ , R v and T sl (see Table 1). Thenumerical code written to solve the slurry equations is freely available online . Unless otherwise specified, all other physical properties of the F-layer areassumed constant, as specified in Table A.2. P e varies because of the range of seismic estimates of the F-layer thickness, d , and the unknown freezing speed, v f . We obtain a proxy for v f through theICB heat flow, Q s , which we presume cannot be much greater than the adiabaticvalue estimated as Q a = 1 . Q s upto 3 . Q a in theliterature. St varies due to Q s and Q sl . While there is no geophysical constrainton Q sl , we test Q sl over a wide range up to 12 TW which is verified a posteriori https://github.com/jnywong/nondim-slurry
12o give 5 ≤ Q c ≤
17 TW [7, 29] (see Appendix C for the calculation of Q c ). Le depends on the thermal conductivity and we explore two values that representthe higher and the lower estimates given in the literature [13, 21, 30]. Li p and R ρ change according to the layer thickness and also the reference density, ρ sl ,though the density on the solid side of the ICB remains the same and is takenfrom PREM, ρ s − = 12 ,
764 kgm − . Li ξ depends on ξ sl , wherein the preciseamount of oxygen present in the bulk of the Earth’s core is difficult to constrainwith reported values varying between 0 . . ξ sl between 2 . . R v , which depends on both the layerthickness and the CSB oxygen concentration and in turn affects the change inspecific volume between solid iron and the liquid iron alloy.Input parameter Symbol Units F-layerLayer thickness d km 150 , , , , , k W m − K − , Q s TW 0 – 3 . Q sl TW 0 – 12Sedimentation prefactor k φ kgm − s 10 − – 10 − CSB liquidus temperature T sl K 4 ,
500 – 6 , ξ sl mol . % 2 – 12P´eclet number P e St Le
354 – 1196Liquidus number (pressure) Li p .
163 – 0 . Li ξ .
004 – 0 . R ρ .
940 – 0 . R v .
203 – 0 . Table 1: Relevant dimensional and dimensionless parameter range for the F-layer, taking otherphysical parameters as constant given in Table A.2.
We constrain T sl in boundary condition (19) from the melting curves of iron,which is usually reported in the range of 5 ,
500 K [32] to 6 ,
350 K [33]. Light13lements present in the iron alloy depress the liquidus temperature of pure iron,and a typical value of ∆ T ξ = 700 K is used, though this can also vary between500 and 1 ,
000 K [31]. We therefore opt to vary T sl between 4 ,
500 and 6 ,
000 K.We constrain the results of the parameter study in three ways:1. by the seismically determined density jump at the ICB [9]. The jumpobtained from normal modes, ∆ ρ mod , has a long wavelength on the orderof hundreds of kilometres and represents the difference between the av-erage densities of the top of the inner core and the bottom of the outercore. The jump obtained from body waves, ∆ ρ bod , has a short wavelengthon the order of several kilometres, and therefore represents the differencein densities either side of the ICB. The difference between the normalmode and body wave estimates therefore points to a density anomalycaused by the F-layer. Normal mode studies suggest 600 kg m − ≤ ∆ ρ mod ≤ ±
180 kg m − [2, 34], whereas body wave studies suggest520 ±
240 kg m − ≤ ∆ ρ bod ≤ ,
100 kg m − [35, 36]. Hence solutionsobtained from the slurry model should satisfy a maximum density jumpacross the layer of max(∆ ρ mod − ∆ ρ bod ) = 1 , −
280 = 720 kg m − ,and a minimum bound min(∆ ρ mod − ∆ ρ bod ) < ρ s + − ρ sl , where ρ s + is the density on the slurry side of the ICB. The14otal density is ρ = ρ H + ρ (cid:48) , (25)where ρ H = (cid:18) gK ( r − r sl ) + 1 ρ sl (cid:19) − , (26)is the hydrostatic part, in which K denotes the bulk modulus, and ρ (cid:48) = ρ sl (cid:104) − αT (cid:48) − α ξ ξ l (cid:48) + ( α φ + α ξ ξ l ) φ (cid:48) (cid:105) . (27)is the density perturbation, where T (cid:48) , ξ l (cid:48) and φ (cid:48) are perturbations in the de-composition T = T sl + T (cid:48) ,ξ l = ξ sl + ξ l (cid:48) ,φ = φ sl + φ (cid:48) , with φ sl = 0 at the CSB. In the equation of state (27), the first two terms arewell known from double diffusive/thermochemical convection theory, whereasthe last term is unique to the slurry. The expansion coefficient of the solid isdefined as α φ = ρ sl ∆ V s,lF e,O , (see derivation in Appendix D). To determine φ (cid:48) , the dimensional solid fluxyielded from the solution of the slurry equations is related to the solid fractionby j = b ( φ )∆ V s,lF e,O ∇ p. (28)15e employ a Stokes’ flow model of mobility [37] where we assume that thesolid iron particles created in the slurry layer drift towards the ICB under theinfluence of gravity. Consequently the sedimentation coefficient is given by b ( φ ) = k φ φ / = (cid:18) ρ s − ( ρ sl ) π ν N (cid:19) / φ / , (29)where N is the number of particles in a unit volume. We write k φ to groupthe prefactors that multiply the solid fraction. This k φ sedimentation prefactorintroduces a degree of freedom since N and the value of the kinematic viscosity, ν , in the slurry is unknown. We further consider this issue in Section 3.2, wherewe conduct a sensitivity analysis on k φ .We evaluate the gradient of the density perturbations to assess the stabilityof the layer. Solving the slurry equations yields T , ξ l and φ and subtractingthe reference values from these quantities gives the perturbations T (cid:48) , ξ (cid:48) and φ (cid:48) ,which allows us to evaluate the density perturbation (27). If the layer is stablethen the change in density relative to the hydrostatic reference should decreaseas the radius increases, i.e. d ρ (cid:48) / d r <
0, and vice versa if the layer is unstable,i.e. d ρ (cid:48) / d r >
0. We relax the condition for stability slightly by stipulating thatthe layer is partially stable if d ρ (cid:48) / d r <
3. Results and discussion
Figure 2 shows example solutions for a wide range of layer thicknesses be-tween 150 and 400 km, where Q s = 2 . v f = 0 .
44 mmyr − ), Q sl = 5 TW, k = 100 Wm − K − and T sl = 5 ,
547 K taken from [25, 27], for an 82%Fe-8%O-10%Si mixture. For all layer thicknesses the temperature gradient is strictly16 T e m p e r a t u r e ( K ) (a) Liquidus (Davies et al. 2015) O xy g e n ( m o l . % ) (b) S o li d f l u x ( k g m s )
1e 7 (c) D e n s i t y ( k g m ) (d)
150 km (box)150 km200 km250 km300 km350 km400 kmPREM
Figure 2: Profiles of (a) temperature, (b) oxygen concentration, (c) solid flux and (d) densityacross the slurry with layer thicknesses between 150 and 400 km. The dashed line in (a)is the uniform composition liquidus determined from ab initio calculations [25, 27], and thedash-dotted line in (d) is the PREM density [2]. The blue dotted line is the equivalent boxmodel case from W18 for d = 150 km (see text). Control parameters are Q s = 2 . Q sl = 5 . k = 100 Wm − K − and T sl = 5 ,
547 K (colour online). negative so the slurry is thermally destabilising, whereas the compositional gra-dient may either be stable or unstable. Solid flux is always negative in thedirection towards the ICB, achieving its largest magnitude at the ICB itselfand then vanishing to zero at the top of the layer as imposed by the boundarycondition (24). Temperature, oxygen and solid flux all contribute to the overalldensity across the layer and the layer is stably-stratified when 150 ≤ d ≤
300 km,however for d ≥
350 km the layer is unstable.We observe in Figure 2(c) that the solid flux eventually increases exponen-tially with radius under the influence of turbulent mixing from the bulk of theliquid outer core. As the solid flux diminishes its gradient sharply increases,17quivalent to precipitating more iron particles. The latent heat release associ-ated with the phase change in the sublayer prompts the slurry to depress thetemperature in Figure 2(a) to stay on the liquidus and the oxygen concentrationgradient increases in response to this as seen in Figure 2(b).For completeness, we compare the solution in spherical geometry for d =150 km with the equivalent boundary conditions from the Cartesian geometry ofW18, except for the box model we have shifted the CSB temperature to coincidewith the liquidus value from [25, 27]. The main difference is that the effect ofspherical geometry suppresses the effect of the mixing sublayer, since Figure2(c) shows that the solid flux gradient in the upper part of the layer is shallowerthan the box model case. Melting at the base of the layer in the box model ishighlighted by the negative gradient in the solid flux, which is not present in thespherical model, and contributes positively to the density variations at the baseof the layer in the box model. Despite this, turbulent mixing in the sublayerdominates to depress the temperature and oxygen variations, producing a muchsmaller density difference overall in the box model.The control parameters P e and Li p depend linearly on d and parameters R ρ , R v and Le depend indirectly on d through ρ sl . For a sense of scale in thevariation of parameters from d = 150 km to 400 km for the reference solutiongiven in Figure 2, P e increases by 18%, Li p increases by 35%, R v increasesby 14% and the changes in R ρ and Le are less than 1 . F and the strength of the turbulent mixing layerthrough the barodiffusion term in (15), as seen in Figure 2(c), which combineto reduce the magnitude of the density variations in the layer overall. Figure3 shows more clearly that contributions to dρ (cid:48) / d r from the oxygen gradient for18
300 1400 1500 1600Radius (km)1.251.000.750.500.250.000.25 d d r = sl ( T d T d r d d r + ( + ) d )d rsl T d T d rsl d d rsl ( + ) d d r
150 400Layer thickness (km)
Figure 3: Gradients of the density variation (grey), d ρ (cid:48) / d r , separated into contributionsfrom the temperature (red), − ρ sl α T d T (cid:48) / d r , oxygen (blue), − ρ sl α ξ d ξ (cid:48) / d r , and solid fraction(orange), ρ sl ( α φ + α ξ ξ )d φ (cid:48) / d r , across the slurry for layer thicknesses of 150 (dotted) and400 km (solid). If d ρ (cid:48) / d r < ρ (cid:48) / d r > Q s = 2 . Q sl = 5 . k = 100 Wm − K − (colour online). d = 400 km can become positive in the mid-depths, since the pressure part of theliquidus relation (17) outweighs the temperature part, before becoming negativeagain under the stabilising influence of the mixing sublayer. Oxygen variationspredominantly control the gradient of the density variations. When d = 400 kmthe sign of d ρ (cid:48) / d r changes from negative (stable), to positive (unstable), andback to negative again, producing an “S”-shaped density profile seen in Figure2(d). On the other hand, Figure 3 shows that for d = 150 km the layer is stablesince d ρ (cid:48) / d r < φ (cid:48) , are negligible over the majority of the layer, perhaps apart from a very thinregion at the top (see Appendix B for further discussion on d φ ).19 .2. Sensitivity analysis We evaluate the sensitivity of the slurry system to the sedimentation pref-actor, k φ , CSB temperature, T sl , and CSB oxygen concentration, ξ sl . We varythe value of k φ by a large range between 10 − and 10 − kgm − s and recomputethe example given in Figure 2(a) where the layer thickness is fixed at 150 km.Figure 4(a) shows that density stratification is increased with smaller values of k φ and the solutions converge as k φ exceeds 10 − kgm − s. Figure 4(b) showsthe corresponding solid fraction, φ , profiles within the layer. It can be seenthat for the smallest value of k φ , the solid fraction is close to the rheologicaltransition at φ m = 0 . φ (cid:28) k φ = 10 − kgm − s that is high enough so that thedensity becomes independent of k φ and consistent with φ (cid:28) d = 150 kmand vary T sl between 4 ,
500 and 6 ,
000 K with intervals of 100 K and for theCSB oxygen concentration we vary its value between 2 . . . T sl and ξ sl , respectively. We find that no solutions wereobtained for values below ξ sl = 2 . ξ sl , there appears to be a limited effect on the resulting density pro-files, whereas Figure 5(a) shows that changing T sl introduces a greater spreadin the solutions obtained. In terms of the dimensionless parameters, varying2 . . % < ξ sl < . . % through 0 . < Li ξ < .
028 does not signifi-cantly affect the system, where Li ξ enters into the pressure part of the liquidusrelation and the barodiffusion term of equation (15). Changing T sl in boundarycondition (19) shifts the anchor point of the temperature solution and will notaffect the curvature of the temperature profile since the heat fluxes into and outof the slurry remain fixed, so the temperature perturbations are the same for20 D e n s i t y ( k g m ) (a) S o li d f r a c t i o n (b) k = 10 kgm sk = 10 kgm sk = 10 kgm sk = 10 kgm sk = 10 kgm s Figure 4: (a) Density and (b) solid fraction profiles across the layer for different values of k φ .The dashed line in (a) is PREM and the dashed in (b) is φ m = 0 . all T sl . However, T sl does affect the oxygen concentration through the liquidusrelation, where a higher T sl will decrease the oxygen concentration gradient andcreate smaller perturbations in the oxygen concentration, therefore generatingsmaller density perturbations to produce a lower slurry density overall. A lower T sl has the opposite effect to produce a higher slurry density. The spread in theICB density with different T sl is roughly 60 kgm − , and there is up to a 40%change in the density jump across the whole layer with respect to the examplecase. This sensitivity may have an impact on the overall stability of the layer,however for the regime diagram we shall continue to use the default value fromthe literature of T sl = 5 ,
547 K for a layer 150 km thick.21
225 1250 1275 1300 1325 1350Radius (km)12100121501220012250 D e n s i t y ( k g m ) (a) T sl = 4500 K T sl = 5457 K T sl = 6000 K sl = 2.0 mol.% sl = 8.0 mol.% sl = 12.0 mol.%PREM Figure 5: Density profiles across the layer for (a) 4 , < T sl < ,
000 K and (b) 2 . < ξ sl < . . %. Default values of T sl = 5 ,
547 K from [25, 27] in (a) and ξ sl = 8 . . % in (b)are given by the grey lines. All other input parameters are the same as the example case, andthe layer thickness is fixed at 150 km (colour online). We present a regime diagram of the
P e, St –space by varying Q s and Q sl .There are two cases: high thermal conductivity, k = 100 Wm − K − , corre-sponding to Le = 1181 and low thermal conductivity, k = 30 Wm − K − ,corresponding to Le = 354. By fixing d = 150 km and ξ sl = 8 mol . %, then Li p = 0 . Li ξ = 0 . R ρ = 0 .
952 and R v = 0 . ρ (cid:48) / d r < ρ (cid:48) / d r < ρ (cid:48) / d r >
0) slurries, andalso areas with no slurry where no solution is found.No slurry can exist when the latent heat release at the ICB interface exceedsthe CSB heat flux, so thatif Q sl < Q s then St = q sl q s < (cid:18) r i r sl (cid:19) ≡ St ∗ which corresponds to a critical Stefan number of 0 . ≤ St ∗ ≤ . ≤ .00.51.01.52.02.53.0 S t NO SLURRYUNSTABLE STABLEPS (a) Le =1181 NO SLURRY U N S T A B L E STABLE PSPS (b) Le =354 Pe S t NO SLURRYUNSTABLE STABLEPS (c) Pe NO SLURRY U N S T A B L E STABLE PSPS (d) C M B h e a t f l u x ( T W ) d e n s i t y j u m p ( k g m ) Figure 6: (Left column) Regime diagrams for the high Lewis number case, (right column)and the low Lewis number case with d = 150 km. (Top row) Contours of the CMB heat fluxconstrained by 5 < Q c <
17 TW and (bottom row) contours of the density jump, ρ s + − ρ sl .The phase space is divided into stable (contour fill), partially stable, denoted PS (contour fill,dashed line), unstable (grey) and no slurry (white) regions. The yellow star represents theexample case from Figure 2 (colour online). d ≤
400 km. When the layer thickness is 150 km, above St ∗ = 0 . P e and stable solutions at higher
P e , andthis transition generally occurs when
P e T ≡ P eLe = v f r sl κ (cid:39) , (30)where P e T is the thermal P´eclet number and κ = k/ρ sl c p is the thermal dif-fusivity. P e T controls the boundary condition on the temperature gradient in2320) and (21). An unstable slurry develops when P e < Le because more heatis conducted through the slurry which is unavailable for equilibration throughthe liquidus constraint, hence creating a smaller oxygen concentration differencethat produces a smaller density contrast that is insufficient to stabilise the layer.A stable slurry develops when
P e > Le as the rate of advection from the com-pacting layer overcomes the thermal diffusion rate. Increasing
P e T stabilises theslurry layer since the temperature gradient steepens at the boundaries, thereforegenerating greater positive variations in the oxygen gradient that enhance thedensity anomaly. If d = 150 km then the transition at P e (cid:39) Le correspondsto Q s = 1 . Le and Q s = 0 . Le . Figure 6(c) and(d) show that the density jump across the layer in all cases is well below themaximum limit of 720 kgm − and also suggests that the density jump roughlyscales with P e , with higher values encountered at low Le .A higher St steepens the CSB temperature gradient relative to the ICB, andcondition (30) is met by strengthening the turbulent sublayer at the top of theslurry by increasing the mixing parameter F . Within the sublayer the rate ofcrystallisation is increased, therefore the density jump across the layer increaseswith St as seen in Figure 6(c) and (d). The critical P e T where the slurrytransitions from unstable to stable slightly decreases as St increases and thelayer is more likely to become partially stable. We find that the dimensionlessICB speed, and thereby the magnitude of the solid flux at the ICB, is roughlyproportional to St . The ICB speed affects contributions to the secular cooling,gravitational power and latent heat release in the core and therefore influencesthe total CMB heat flow. This is reflected in Figure 6(a) and (b) on the totalCMB heat flux, which scales linearly with P eSt .24 .4. Seismic properties of the F-layer and the inner core age P w a v e s p ee d ( k m / s ) slurryOhtaki et al. (2015)PREMak135 Figure 7: P wave speed in the F-layer derived from the high Le example solution with d =150 km , k = 100 Wm − K − , Q s = 2 . Q sl = 5 . et al. [4](blue, dotted) (colour online). Motivated by the seismic observations, we have demonstrated that a slurry isable to produce stable layers with varying degrees of stratification depending onthe parameters selected. If we invert the process, then what constraints can themodel provide on observations? From our slurry model we can determine thedensity at every point across the layer and hence the P wave speed, v p = K/ρ .The bulk modulus, K , depends on the core chemistry and therefore an appro-priate equation of state should be applied [19, 39]. For the sake of simplicity weapproximate the bulk modulus by PREM in our calculations.Figure 7 illustrates the v p profile of the example solution from Section 3.1compared with the seismic models of PREM [2], ak135 [40] and the FVW modelof Ohtaki et al. [4]. The P wave speed from the slurry model is reduced relativeto PREM by up to 0 .
13% and there is a difference of 0 .
1% compared withOhtaki et al. [4], which reported that v p ( r i ) = 10 . − and d v p / d r =25 . × − s − . For this example slurry solution we obtain a density jumpacross the layer of ρ s + − ρ sl = 112 kgm − and calculate that Q c = 10 . ρ s + = 12 ,
198 kgm − , therefore we predict that the density jump frombody waves for this particular slurry is ∆ ρ bod = ρ s − − ρ s + = 502 kgm − . Thedensity jump from normal modes is ∆ ρ mod = ∆ ρ bod + ρ s + − ρ sl = 628 kgm − ,and this is fixed in all solutions of the slurry by choosing ρ s − and ρ sl from PREM.Note that the example solution given above is one of many compatiblesolutions given by the regime diagram. By considering the entire range ofparameters permitted by the regime diagram that fit with the geophysicalconstraints, the high Le case gives 81 ≤ ρ s + − ρ sl ≤
140 kg m − and thelow Le case gives 101 ≤ ρ s + − ρ sl ≤
331 kg m − . This leads to bounds of475 ≤ ∆ ρ bod ≤
534 kg m − and 283 ≤ ∆ ρ bod ≤
513 kg m − for high and low Le respectively. The observations available give 520 ± ≤ ∆ ρ bod ≤ ,
100 kg m − ,therefore our model suggests that a stably-stratified F-layer is opposed to thehigher values of ∆ ρ bod recommended by the observations and that the slurrymodel gives an upper bound of ∆ ρ bod ≤
534 kg m − for both Le cases.We can approximate the inner core age using the model output of the ICBspeed, v , by assuming that inner core growth is proportional to the square-rootof time [41], r i ( τ ) = r i (cid:16) ττ i (cid:17) , (31)where τ is time normalised by the age of the inner core, τ i , and r i is the present-day IC radius. By differentiating (31) and considering the present-day where τ = 0, we have τ i = r i v , (32)26here we recall that v ≡ v f + v s = d r i / d τ is the inner core growth rate.Estimating τ i is of significant geophysical interest and its value is subject todebate since τ i is directly related to the evolution of the core global heat balance,which is significantly affected by the thermal conductivity [21]. For the samerange of solutions considered above to estimate ∆ ρ bod , our results from theregime diagram for high Le give 0 . ≤ v ≤ .
98 mm yr − , which is relativelyfast compared to the low Le solutions where 0 . ≤ v ≤ .
31 mm yr − . Thisyields an inner core age of 0 . ≤ τ i ≤ . Le , which is in linewith other estimates from the literature [7], and 2 . ≤ τ i ≤ . Le ,which is mostly beyond oldest estimates of 2 .
4. Conclusions
In this study, we have investigated the conditions that produce a stably-stratified slurry layer at the base of the Earth’s outer core. We non-dimensionalisethe governing equations given by W18 and elucidate the key dimensionless pa-rameters that control the system behaviour. By varying the P´eclet and Stefannumbers, we map a regime diagram of the slurry demarcating the conditionsthat favour stable, partially stable and unstable density configurations, as wellas no slurry. Solutions are obtained for high and low Le numbers, which reflecthigh or low thermal conductivity values of the core. We constrain the resultsby evaluating the density jump across the layer and the total CMB heat flux.Our main result is that stably-stratified solutions can be produced by a slurryfor a wide range of parameters that span plausible values for Earth’s core. Wehave identified regions of the parameter space containing many solutions that27re compatible with observations of the F-layer, and we have also establishedconditions that are not suitable. We find that • P e T = P eLe (cid:39) • higher P´eclet number slurries generally facilitate stable density stratifica-tion because this controls the steepness of the temperature gradient at theboundaries, which in turn increases the magnitude of density variations inthe slurry, • the Stefan number has a stabilising role through crystallising more parti-cles in the turbulent mixing sublayer, • no slurry can exist with St < St ∗ = ( r i /r sl ) since the heat flow at thebottom of the layer is greater than at the top of the layer • the slurry model suggests that ∆ ρ bod ≤
534 kg m − for high and low Le • estimates of the inner core age suggests a slurry with high core conduc-tivity is compatible with independent estimates from the literatureDensity perturbations of compositional origin dominate contributions to theoverall stability of the layer, and can be destabilised through increasing the layerthickness. We have also investigated the sensitivity to the CSB temperature andoxygen concentration and deduced a limit for the prefactor k φ = 10 − kgm − sthat defines the sedimentation coefficient, b ( φ ), which relates the solid flux, j , tothe solid fraction, φ , through Stokes’ flow. We find that the model is insensitiveto the concentration of oxygen at the CSB, however, the effect of the liquidustemperature is significant to the overall stability of the layer.From mineral physics or seismology of model uncertainties, improved esti-mates such as the layer thickness, CSB temperature and the CSB oxygen con-centration, will benefit the slurry model greatly. Further progress determining28he liquidus curve for an iron mixture, using experiments or first-principle cal-culations, could improve the temperature condition at the CSB, and improvedestimates of the density jump across the layer from normal and body wavestudies would help constrain the space of geophysically consistent solutions.We examined the P wave speed in the slurry layer to show that the modelcan produce a profile that compares well with other seismic models. This couldbe further verified comprehensively so that the model may provide a useful toolfor corroborating observations, such as constraining the seismic density jumpat the ICB derived from body waves, ∆ ρ bod . There is also the potential forthe slurry to provide improved and independent estimates of the inner core agestipulated by having a stable F-layer that agrees with the seismic and heat flowrequirements.We implemented a simple compacting layer on the solid side of the IC wheresolid particles produced by the slurry accumulate and instantly compact. Re-alistically the physics of this process is extremely complicated. Studies suggestthat a mush solidification regime can occur [6] where liquid channels permeate amatrix of solid with a high solid fraction. A more advanced model may seek toincorporate this process. The timescale of interest in our slurry is comparableto the slow growth of the inner core, where we have applied the fast meltinglimit. On short timescales at the microscopic level, the effect of supercoolingon the nucleation of solid iron at core conditions may come into play and thistopic is an active area of interest [42, 43, 44].Future work could focus on assessing the dominant balances between termsin the slurry equations that are responsible for stable density stratification sothat the physics controlling the transitions between different regimes can beelucidated. For example, the importance of the mixing sublayer and its influenceon the production of solid phase needs to be quantitatively evaluated. A more29ophisticated model coupling the slurry to outer core convection, perhaps similarto a recent study conducted by Bouffard et al. [45] for the stably-stratified layerat the top of the core, may fully capture the physics of this process and shedlight on the effect of entrainment.Our study also opted to keep the layer thickness constant over time in orderto ascertain solutions that are concomitant with the present-day seismic obser-vations. Relaxing the layer thickness requires an extra constraint on the modelto be developed in its place [46]. Constructing a fully time-dependent frame-work could provide insight into the factors controlling the growth and declineof a slurry layer, which may shed light on how the F-layer came into existenceover the core’s history.
5. Acknowledgements
The numerical code used to solve the slurry equations in this paper is freelyavailable at https://github.com/jnywong/nondim-slurry . We thank MarineLasbleis and an anonymous reviewer for their constructive comments that helpedimprove this paper. JW acknowledges support from the Fondation Simone andCino Del Duca of Institut de France and the Engineering and Physical SciencesResearch Council (EPSRC) Centre for Doctoral Training in Fluid Dynamics(EP/L01615X/1). CJD is supported by the Natural Environment ResearchCouncil (NERC) Independent Research Fellowship (NE/L011328/1). Figureswere produced using Matplotlib [47]. 30 ppendix A. Table of valuesSymbol Definition Value Units Source p Pressure kgm − s − T Temperature K ξ Light element concen-tration Mass frac-tion ξ l Light element concen-tration in the liquidphase Mass frac-tion φ Solid fraction Mass frac-tion j Solid flux kgm − s − ρ Density kgm − F Mixing parameter v ICB speed ms − v f Freezing speed ms − v s Snow speed ms − µ Chemical potential Jkg − ρ H Hydrostatic density kgm − ρ (cid:48) Density perturbation kgm − T (cid:48) Temperature perturba-tion K ξ l (cid:48) Light element perturba-tion in the liquid phase Mass frac-tion φ (cid:48) Solid fraction perturba-tion Mass frac-tion31 ( φ ) Sedimentation coeffi-cient kgm − s k φ Sedimentation coeffi-cient prefactor kgm − s ν Kinematic viscosity m s − N No. of solid iron parti-cles per unit volumeΦ Gibbs free energy
P e
P´eclet number St Stefan number Le Lewis number Li p Liquidus number (pres-sure) Li ξ Liquidus number (com-position) R ρ Ratio between liquidand solid iron R v Ratio between thechange in specificvolumes upon phasechange from pure ironand iron alloy R Ideal gas constant 8 .
31 Jkg − mol . − a O Atomic weight of oxy-gen 16 Da r i ICB radius 1221 × m PREM [2]32 Layer thickness (150 , , , , × m [3], [5] r sl CSB radius 1371–1621 × m [3], [5] ρ s − Density of iron on thesolid side of the ICB 12 . × kgm − PREM [2] ρ s + Density of iron on theliquid side of the ICB kgm − ρ sl Density of liquid iron at r = r sl kg m − PREM [2] g Gravity ms − PREM [2] K Bulk modulus kgm − s − PREM [2] ξ sl Oxygen concentrationin the bulk of the liquidcore 2–12 mol . % [31] T sl Liquidus temperatureat the CSB 4 , ,
000 K [27, 25] c p Specific heat capacity 715 Jkg − K − [23] α Thermal expansion co-efficient 1 × − K − [23] α ξ Compositional ex-pansion coefficient ofoxygen 1 . L Latent heat of fusion 0 . × Jkg − [23] k Thermal conductivity 30, 100 Wm − K − [21], [27]33 D Modified self–diffusioncoefficient of oxygen m s − D O Self–diffusion coefficientof oxygen 0 . × − m s − [25] V sF e Specific volume of solidiron m kg − Ideal solutiontheory V lF e,O Specific volume of liq-uid iron and oxygen m kg − Ideal solutiontheory∆ V s,lF e,O Change in specific vol-ume between liquid andsolid phase m kg − Ideal solutiontheory∆ V s,lF e Change in specific vol-ume between liquid ironand solid iron m kg − Ideal solutiontheory α φ Expansion coefficient ofsolid Ideal solutiontheory q s ICB heat flow per unitarea Wm − q sl CSB heat flow per unitarea Wm − Q s ICB heat flow TW Q sl CSB heat flow TW Q c CMB heat flow 5–17 TW [29], [7] Q s Secular cooling TW Q g Gravitational power TW Q l Latent heat power TW34 able A.2: Symbols and values (if applicable) used in the slurry model. Group 1: fundamentalthermodynamic variables. Group 2: slurry notation. Group 3: Dimensionless numbers. Group4: Physical constants. Group 5: Seismic properties. Group 6: Properties determined from ab initio calculations and high-pressure experiments. Group 7: Properties that assume idealsolution theory. Group 8: Heat flow estimates.
Appendix B. Variations in the solid fraction, d φ We suppose that the variations in the solid fraction, d φ , are negligible. Herewe assess the validity of this assumption a posteriori . The variation in solidfraction enters the energy equation (2) and the equation of state (27). In thefirst instance, including these variations adds an extra term to the latent heatrelease due to phase change. The dimensionless energy equation (16) becomes − ˆ v d ˆ T dˆ r = LeP e (cid:32) d ˆ T dˆ r + 2ˆ r d ˆ T dˆ r (cid:33) + 1 St (cid:18) dˆ dˆ r + 2ˆ r ˆ − R ρ d φ dˆ r (cid:19) . Using the control parameters from Figure 3, we compare the order of magnitudeof each part of the last term in the expression above where the d φ/ dˆ r term ap-pears. In Figure B.8 we can see that with r i /r sl = 0 .
89 and 0 .
75 (correspondingto d = 150 and 400 km, respectively), the extra term due to the variations in φ is almost three orders of magnitude smaller compared with the other two termsin the majority of the layer. We observe that | r ˆ | ∼ | R ρ d φ dˆ r | approaching the topof the layer, which is a consequence of ˆ tending to 0 as dictated by boundarycondition ˆ (1) = 0 in equation (24). This occurs within a very thin region atthe top of the layer and so | R ρ d φ dˆ r | makes a negligible difference to the overallenergy balance in the slurry.Turning to the equation of state we observe that the contributions from thevariations in solid also arise by imposing boundary condition (24), which can35 .75 0.80 0.85 0.90 0.95 1.00 r | dj / dr | |2 j / r | | R d / dr | Figure B.8: Comparing magnitudes | dˆ dˆ r | , | r ˆ | and | R ρ d φ dˆ r | , across the slurry for layer thick-nesses of 150 (dashed, blue fill) and 400 km (solid). Control parameters are Q s = 2 . Q sl = 5 . k = 100 Wm − K − (colour online). be evaluated analytically by assuming Stokes’ flow and writinglim ˆ r → (cid:18) d φ dˆ r (cid:19) = lim ˆ r → (cid:32) − (cid:18) ρ s − v f K φ (cid:19) ( − ˆ ) − dˆ dˆ r (cid:33) = −∞ On a physical basis Stokes’ flow is unlikely to hold in the turbulent mixingregion close to the CSB, and the contributions from the variations in φ in thisthin region are considered to be a numerical artefact of maintaining ˆ (1) = 0.This thin region comprises less than 5% of the layer and has a limited impacton assessing the overall stability of the slurry. Appendix C. CMB heat flow, Q c In general contributions to the CMB heat flux can be separated into thesecular cooling, Q s , gravitational power, Q g , and the latent heat Q l , whilepressure freezing is neglected and radiogenic heating is ignored for the sake of36implicity [24]. This gives a core heat balance of Q c = Q ls + Q lg + Q sl , (C.1)where Q ls = − c p T c d T c d t (cid:90) V l ρT a d V l (C.2)is the secular cooling in the liquid volume, V l , and Q lg = (cid:90) V l ρψα ξ d ξ d t d V l − π ( r sl ) ρ sl ψ sl α ξ ξ sl v (C.3)is the gravitational power in the liquid volume, and Q sl is the heat flux imposedat the CSB that also contains Q l . In (C.2), the adiabatic temperature is givenby [23] T a ( r ) = T sl exp (cid:18) − (cid:90) rr sl gγφ d r (cid:19) . (C.4)We calculate the core cooling rate, d T c / d t , by constructing an adiabat anchoredat the present day CSB temperature, followed by constructing a new adiabatanchored at a new CSB temperature after time ∆ t due to the advancing layer,and then finding the resulting decrease in the CMB temperature, T c = T a ( r c ).This is given byd T c d t ≈ ∆ T c ∆ t = T sl ( r sl + v ∆ t ) exp (cid:16) − (cid:82) r c r sl + v ∆ t gγφ d r (cid:17) − T sl ( r sl ) exp (cid:16) − (cid:82) r c r sl gγφ d r (cid:17) ∆ t . The gravitational power (C.3), is composed of two parts: the first part is pro-portional to the change in oxygen concentration in the bulk of the volume, V l ,and the second part is from the motion of the CSB, where there is no additional37ontribution from the motion of the ICB since the CSB and ICB move at thesame rate. By mass conservation, the change in oxygen concentration in thebulk is given by (cid:90) V l ∂ξ∂t d V l = − π ( r sl ) ρ sl ξ sl vM l , where M l = (cid:82) V l ρ d V l is the mass of the liquid outer core. Appendix D. Definition of the solid expansion coefficient, α φ The coefficient, α φ , is a dimensionless expansion coefficient that influencesthe contribution of the solid fraction to the slurry density. From equation (A3)of W18, we have that α φ ≡ − ρ sl (cid:18) ∂V∂φ (cid:19) p,T,ξ . (D.1)From (A1) of W18, the expression for the Gibbs free energy, d Φ, defines thespecific volume as V ≡ (cid:18) ∂ Φ ∂p (cid:19) T,ξ,φ . (D.2)The lever rule in equation (A9) of W18 gives (cid:18) ∂ Φ ∂φ (cid:19) p,T,ξ = Φ s − Φ l (D.3)and (cid:18) ∂ Φ s ∂p (cid:19) T,ξ,φ = V sF e , (cid:18) ∂ Φ l ∂p (cid:19) T,ξ,φ = V lF e,O , (D.4)38here Φ s is the solid part of the Gibbs free energy and Φ l is the liquid part.Substituting (D.2), (D.3) and (D.4) into (D.1) gives the final result α φ = − ρ sl (cid:0) V sF e − V lF e,O (cid:1) = ρ sl ∆ V s,lF e,O . (D.5) References [1] C. M. Hardy, J. Wong, Stably stratified layers within Earth’s core, Astron-omy & Geophysics 60 (3) (2019) 3–30.[2] A. Dziewonski, D. Anderson, Preliminary Reference Earth Model, Phys.Earth Planet. Int. 25 (1981) 297–356.[3] A. Souriau, G. Poupinet, The velocity profile at the base of the liquid corefrom PKP(BC+Cdiff) data: an argument in favor of radial inhomogeneity,Geophys. Res. Lett. 18 (1991) 2023–2026.[4] T. Ohtaki, S. Kaneshima, Independent estimate of velocity structure ofEarth’s lowermost outer core beneath the northeast Pacific from PKiKP–PKPbc differential traveltime and dispersion in PKPbc, Journal of Geo-physical Research: Solid Earth 120 (11) (2015) 7572–7586.[5] Z. Zou, K. Koper, V. Cormier, The structure of the base of the outer coreinferred from seismic waves diffracted around the inner core, J. Geophys.Res. 113 (2008) B05314.[6] R. Deguen, Structure and dynamics of Earth’s inner core, Earth Planet.Sci. Lett. 333–334 (2012) 211–225.[7] F. Nimmo, Energetics of the core, in: G. Schubert (Ed.), Treatise on Geo-physics 2nd Edn, Vol. 9, Elsevier, Amsterdam, 2015, pp. 31–65.398] T. Alboussi`ere, R. Deguen, Asymmetric dynamics of the inner core andimpact on the outer core, J. Geodyn. 61 (2012) 172–182.[9] D. Gubbins, G. Masters, F. Nimmo, A thermochemical boundary layer atthe base of Earth’s outer core and independent estimate of core heat flux,Geophys. J. Int. 174 (2008) 1007–1018.[10] D. Loper, P. Roberts, On the motion of an iron-alloy core containing aslurry: I. General theory, Geophys. Astrophys. Fluid Dyn. 9 (1977) 289–321.[11] J. Wong, C. J. Davies, C. A. Jones, A Boussinesq slurry model of the F–layer at the base of Earth’s outer core, Geophysical Journal International.[12] N. de Koker, G. Steinle-Neumann, V. Vojtech, Electrical resistivity andthermal conductivity of liquid Fe alloys at high P and T and heat flux inEarth’s core, Proc. Natl. Acad. Sci. 109 (2012) 4070–4073.[13] M. Pozzo, C. Davies, D. Gubbins, D. Alf`e, Thermal and electrical conduc-tivity of iron at Earth’s core conditions, Nature 485 (2012) 355–358.[14] M. Pozzo, C. Davies, D. Gubbins, D. Alf`e, Thermal and electrical con-ductivity of solid iron and iron-silicon mixtures at Earth’s core conditions,Earth Planet. Sci. Lett. 393 (2014) 159–164.[15] R. Deguen, T. Alboussi´ere, S. Labrosse, Double-diffusive translation ofEarth’s inner core, Geophysical Journal International 214 (1) (2018) 88–107.[16] D. Loper, P. Roberts, A Boussinesq model of a slurry, Structure and Dy-namics of Partially Solidified Systems (1987) 291–323.[17] F. Birch, Elasticity and the constitution of Earth’s interior, J. Geophys.Res. 66 (1952) 227–286. 4018] D. Alf`e, M. Gillan, G. Price, Composition and temperature of the Earth’score constrained by combining ab initio calculations and seismic data,Earth Planet. Sci. Lett. 195 (2002) 91–98.[19] J. Badro, A. Cˆot´e, J. Brodholt, A seismologically consistent compositionalmodel of Earth’s core, Proc. Natl. Acad. Sci. 111 (2014) 7542–7545.[20] H. Gomi, K. Ohta, K. Hirose, S. Labrosse, R. Caracas, V. Verstraete,J. Hernlund, The high conductivity of iron and thermal evolution of theEarth’s core, Phys. Earth Planet. Int. 224 (2013) 88–103.[21] Q. Williams, The thermal conductivity of Earth’s core: A key geophysicalparameter’s constraints and uncertainties, Annual Review of Earth andPlanetary Sciences 46 (1).[22] R. Deguen, T. Alboussi`ere, D. Brito, On the existence and structure of amush at the inner core boundary of the Earth, Phys. Earth Planet. Int. 164(2007) 36–49.[23] D. Gubbins, D. Alfe, G. Masters, G. Price, M. Gillan, Can the Earth’sdynamo run on heat alone?, Geophys. J. Int. 155 (2003) 609–622.[24] D. Gubbins, D. Alf`e, G. Masters, G. Price, M. Gillan, Gross thermodynam-ics of two-component core convection, Geophys. J. Int. 157 (2004) 1407–1414.[25] M. Pozzo, C. Davies, D. Gubbins, D. Alf`e, Transport properties for liquidsilicon-oxygen-iron mixtures at Earth’s core conditions, Phys. Rev. B 87(2013) 014110.[26] C. Davies, Cooling history of Earth’s core with high thermal conductivity,Phys. Earth Planet. Int. 4127] C. Davies, M. Pozzo, D. Gubbins, D. Alf`e, Constraints from material prop-erties on the dynamics and evolution of Earth’s core, Nat. Geosci. 8 (2015)678–687.[28] P. Virtanen, R. Gommers, T. Oliphant, M. Haberland, T. Reddy, D. Cour-napeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al., Scipy1.0: fundamental algorithms for scientific computing in Python, Naturemethods 17 (3) (2020) 261–272.[29] T. Lay, J. Hernlund, B. Buffett, Core-mantle boundary heat flow, Nat.Geosci. 1 (2008) 25–32.[30] F. Stacey, D. Loper, A revised estimate of the conductivity of iron alloyat high pressure and implications for the core energy balance, Phys. EarthPlanet. Int. 161 (2007) 13–18.[31] K. Hirose, S. Labrosse, J. Hernlund, Compositional state of Earth’s core,Annual Review of Earth and Planetary Sciences 41 (2013) 657–691.[32] R. Sinmyo, K. Hirose, Y. Ohishi, Melting curve of iron to 290 GPa de-termined in a resistance-heated diamond-anvil cell, Earth and PlanetaryScience Letters 510 (2019) 45–52.[33] J. Jackson, W. Sturhahn, M. Lerche, J. Zhao, T. Toellner, E. Ercan Alp,S. Sinogeikin, J. Bass, C. Murphy, J. Wicks, Melting of compressed iron bymonitoring atomic dynamics, Earth Planet. Sci. Lett. 362 (2013) 143–150.[34] G. Masters, D. Gubbins, On the resolution of density within the Earth,Phys. Earth Planet. Int. 140 (2003) 159–167.[35] K. D. Koper, M. Dombrovskaya, Seismic properties of the inner core bound-ary from PKiKP/P amplitude ratios, Earth Planet. Sci. Lett. 237 (2005)680–694. 4236] H. Tkal˘ci´c, B. Kennett, V. Cormier, On the inner–outer core density con-trast from PKiKP/PcP amplitude ratios and uncertainties caused by seis-mic noise, Geophys. J. Int. 179 (2009) 425–443.[37] D. Loper, P. Roberts, On the motion of an iron-alloy core containing aslurry: II. A simple model, Geophys. Astrophys. Fluid Dyn. 16 (1980) 83–127.[38] V. Solomatov, Magma oceans and primordial mantle differentiation., in:G. Schubert, B. Romanowicz, A. Dziewonski (Eds.), Treatise on geophysics,Elsevier, Amsterdam, 2015, pp. 655–693.[39] S. Cottaar, T. Heister, I. Rose, C. Unterborn, Burnman: A lower man-tle mineral physics toolkit, Geochemistry, Geophysics, Geosystems 15 (4)(2014) 1164–1179.[40] B. Kennett, E. Engdahl, R. Buland, Constraints on seismic velocities inthe Earth from traveltimes, Geophys. J. Int. 122 (1995) 108–124.[41] S. Labrosse, Thermal and compositional stratification of the inner core, C.R. Geosci. 346 (2014) 119–129.[42] L. Huguet, J. Van Orman, S. Hauck II, M. Willard, Earth’s inner corenucleation paradox, Earth Planet. Sci. Lett. 487 (2018) 9–20.[43] C. J. Davies, M. Pozzo, D. Alf`e, Assessing the inner core nucleation paradoxwith atomic-scale simulations, Earth and Planetary Science Letters 507(2019) 1–9.[44] M. Lasbleis, M. Kervazo, G. Choblet, The fate of liquids trapped duringthe Earth’s inner core growth, arXiv preprint arXiv:1912.12258.[45] M. Bouffard, M. Landeau, A. Goument, Convective erosion of a primordialstratification atop Earth’s core, Geophysical Research Letters.4346] J. Wong, A slurry model of the f-layer in the earth’s core, Ph.D. thesis,University of Leeds (2018).[47] J. D. Hunter, Matplotlib: A 2D graphics environment, Computing in Sci-ence & Engineering 9 (3) (2007) 90–95. doi:10.1109/MCSE.2007.55doi:10.1109/MCSE.2007.55