Constraining the Si composition and thermal history of Earth's liquid core from ab initio calculations

22 Earth’s core has sustained a global magnetic field for much of the last 4 bil23 lion years which is, at present, sustained by the power associated with inner 24 core growth. High thermal conductivity of the core suggests the solid inner 25 core is young and models of this predict that there is insufficient power from 26 secular cooling to sustain a geodynamo prior to inner core formation. Pre27 cipitation of light elements dissolved into the liquid core offers an alternative 28 power source for the magnetic field in the absence of inner core growth. We 29 present the first ab initio calculations of the silicon partition coefficient at 30 core-mantle boundary conditions and a thermodynamic partitioning model 31 based on interaction parameters which captures previous experimental re32 ∗a.j.wilson1@leeds.ac.uk Preprint submitted to Earth and Planetary Science Letters January 26, 2022 sults. We report our model and its implications for the past and present 33 core composition as well as the effect of silicon precipitation on the early 34 geodynamo. Oxygen competes with silicon in the liquid metal, meaning for 35 one to be abundant, the other must be sparse. We calculate precipitation 36 rates of ∼ 10−4 to 10−6 wt % K for oxygen concentrations of 0.6 to 3.1 37 wt%. Incorporating our partitioning model into a classic thermal evolution 38 model of the core coupled to a parameterised model of the solid mantle, we 39 show that precipitation of Si can satisfy constraints of the present inner core 40 size, convective heat flux of the mantle and mantle temperature, all whilst 41 sustaining a magnetic field until inner core formation, but requires that the 42 initial oxygen content of the core was < 3 wt%. We find that the core inner 43 age is between 840 and 940 Myrs and that the ancient core was hot, with a 44 core mantle boundary temperature of ∼4700 K, 3.5 Ga. 45

Earth's core has sustained a global magnetic field for much of the last 4 bil-23 lion years which is, at present, sustained by the power associated with inner 24 core growth. High thermal conductivity of the core suggests the solid inner 25 core is young and models of this predict that there is insufficient power from 26 secular cooling to sustain a geodynamo prior to inner core formation. Pre-27 cipitation of light elements dissolved into the liquid core offers an alternative 28 power source for the magnetic field in the absence of inner core growth. We the inner core to be far less than 1 Gyrs old and requires the geodynamo to 71 be powered by heat loss from the core alone. This rapidly cooling scenario 72 means the mantle would have been subject to a super-solidus core-mantle 73 boundary (CMB) temperature for much of Earth history (Nimmo, 2015b; 74 Davies, 2015; Labrosse, 2015). Davies and Greenwood (2022) proposed that 75 the presence of a basal magma ocean (BMO) may provide a resolution, al-76 though this approach relies upon the uncertain evolution of the BMO as well 77 as requiring a conductivity at the lower limit of the recent high estimates. 78 In search of an alternate explanation for the long-lived geodynamo, prior implement an interaction parameter model (Ma, 2001) with differing num-98 bers of included elements and interactions. Davies and Greenwood (2022) 99 show that these differences change the onset time and power of precipitation

126
To evaluate the influence of Si precipitation on the geodynamo we use ab 127 initio molecular dynamic simulations of iron-rich liquids and silicate liquids 128 to calculate equilibrium constants at CMB conditions. We compare our ab 129 initio results to the results of a thermodynamic model fit to experimental 130 partitioning data of Si between silicate and metallic liquids using the inter-131 action parameter formulation of Ma (2001). This thermodynamic model is 132 then coupled to a core evolution model to describe Si solubility in the liq-133 uid core as it cools. We then evaluate the influence of precipitation on the 134 thermal history of the core. We conduct density functional theory (Hohenberg and Kohn, 1964;Kohn 137 and Sham, 1965) molecular dynamic simulations of silicate and iron-rich liq-138 uids to calculate the excess chemical potentials of individual chemical com-139 ponents. Chemical potentials (µ i ) can be described as the free-energy change 140 (∂F ) of a system when the quantity of a species is changed (2) Method 2 computes the change in free energy as a result of changing the Here µ i is dependent on v, T and composition. Separating out the configu-158 rational portion of µ i , which plays no role in partitioning, gives whilstμ SiO 2 =μ Si + 2μ O in the liquid, which when rearranged (for a disso-160 ciation reaction) becomes equal the distribution coefficient allowing us to validate our thermodynamic model. 162 We focus on pressures and temperatures most relevant to the CMB (124

163
GPa and 4500-5500 K), as these are the most crucial for the evolution of 164 the core, and also to avoid complications with changes in magnetic moment 165 at shallower conditions. Simulations were run using the VASP code (Kresse 166 and Furthmüller, 1996) in the canonical ensemble using a Nosé thermostat 167 (Nosé, 1984) and with the Brillioun Zone sampled at the Γ point. A timestep 168 of 1 fs was used and runs lasted between 10 and 100 ps. The plane wave 169 cutoff was set to 500 eV and the projector augmented wave method (Kresse 170 and Joubert, 1999) was used with the generalised gradient approximation 171 functional PW91 (Perdew et al., 1992 The equilibrium constant (K) describes a reaction at equilibrium where a i are activities, α i are reaction exponents, x i and x j are the mo- to K by for which K d is given respectively by where ∆F r is the Helmholtz free energy of reaction and ∆G r is the equivalent 195 Gibbs free energy change, representing conditions of constant volume (V ) and 196 pressure (P ), respectively. ∆H r , ∆S r and ∆V r are the changes in enthalpy, 197 entropy and volume with reaction, respectively, and k B is the Boltzmann 198 constant. Eq. 14 is often written where a, b and c describe the entropy (S), enthalpy (H) and volume (V ) 200 changes of reaction, respectively. This naming convention is adopted over 201 traditional thermodynamic notation because these quantities are not exclu-202 sively represented by a, b or c; entropy for example, will have a pressure and 203 temperature dependence which is not captured by a and so, in practice, these 204 effects will be absorbed into b and c. What cannot be absorbed into these 205 parameters is the compositional dependence of the reaction.

206
Combining Eq. 7 and 15 gives the model which is fit to calculated K d from previous experimental partitioning results 208 via a least squares approach. The best-fit parameters are used to calculate 209 the equilibrium concentration of Si for a dissociation reaction (which we later 210 show to be favourable) in the liquid metal using and T 0 = 1873 K. The activity coefficient of the solvent (γ Fe ) is described

244
In Table 1 we report our ab inito results of the calculated excess chemical 245 potential of SiO 2 at 5500 K and 4500 K, both at 124 GPa. We show consis-246 tency with experimental K d at comparable T and P in Fig. 1      including those below the solidus of SiO 2 (e.g. Usui and Tsuchiya (2010)). 338 We assume the core to be well mixed throughout for convenience, although with the overlying mantle. We assume that the core is mixed thoroughly 448 on timescales far shorter than the timestep of our simulation such that the 449 liquid core has no compositional variation nor stable layers.

450
To balance energies in the core, we follow Davies (2015) where, if small 451 terms are ignored, the heat flow across the CMB (Q cmb ) is where Q s is the secular heat stored in the core and Q L is the latent heat 453 release due to inner core growth. Q p is the gravitational energy from mixing 454 the dense, iron-rich residual liquids post precipitation across the outer core where ρ is density, α i ppt is expansivity, ψ is gravitational potential, C ppt is the 456 precipitation rate (see Fig. 5), t is time and V c is volume of the liquid core. core has no effect on the Si concentration of the liquid core. Gubbins et al.

461
(2004) show that the entropy budget of the core can be balanced by where E α is the entropy due to barodiffusion throughout the core which 463 is negligible (Gubbins et al., 2004;Davies, 2015) and so is ignored. E k is 464 the entropy from thermal conduction and the other terms follow the same 465 notation as their energy counterparts.

466
If when evaluated, these entropy sources produce a positive E j , the geo-467 dynamo can be sustained. This presents the difficulty in a high thermal 468 conductivity core, the entropy balance now has a far larger E k , taking power 469 from E j , and so to sustain a magnetic field before inner core growth, one or 470 more of the r.h.s terms must be increased. Because time before inner core 471 nucleation excludes the influence of E L and E g , a more rapidly cooling core 472 (E s ) or precipitation (E ppt ) are needed. 473 We vary the upper to lower mantle viscosity ratio (f visco ) and initial CMB 474 temperature (T t=0 cmb ) to regulate the core temperature such that the final state 475 of our models matches constraints of the present-day core. These constraints    is because for the core to cool sufficiently to produce the inner core of present 491 radius, they fail to consistently sustain a positive E j . We find that including 492 the energy and entropy effects of precipitation can maintain a positive E j for 493 the majority of Earth history preceding inner core formation, also found to be 494 the case by Hirose et al. (2017), producing an older inner core. We find that 495 high initial oxygen concentration cases require f visco < 1 in order to grow the 496 core to present-day size meaning the upper mantle is more viscous than the lower mantle. We do not expect this to be the case in reality (Rudolph et al., 2015), highlighting the requirement for modest O content for Si precipitation 499 to sustain a geodynamo whilst satisfying present-day constraints.  Figure 6: Thermal evolution of the Earth's core with an initial condition of low O concentration (2 mol%) and Si saturation with inclusion (teal) and exclusion (pink) of power and entropy of precipitation (cases D and D P from Table 2, respectively). Inner core radius (a), mid-mantle potential temperature (b, left, solid lines) and CMB temperature (b, right, dotted lines), convective mantle heat flux (c, left, solid lines) and CMB heat flow (c, right, dotted), and core entropy from ohmic dissipation (d). Black dashed lines show present-day target values.
We show the outcomes of thermal history cases from Table 2  to produce a magnetic field from ∼800 Myrs prior to inner core formation 518 (crossed symbols, Fig. 7). 519 Figure 7: Inner core age and core temperature at 3.5 Ga for our model (coloured symbols) with and without convective power from precipitation (connected by dashed lines). Initial Si saturation is shown as red and undersaturation as blue whilst O rich initial conditions are squares and O poor conditions are circles. Numbers within symbols give the CMB heat flow at 3.5 Ga and where models fail to maintain positive E j prior to inner core nucleation, symbols are crossed out. Also shown are the models of Davies and Greenwood (2022) (who examine MgO precipitation) where up, right and down triangles have a ppt. rate of 0, 0.3 and 1.5 ×10 −5 K -1 respectively and colours denote the core properties in terms of the density jump at the ICB (0.6 (light grey) and 1.0 g cm -3 (dark grey)) which represent bounding extremes of the density jump.

531
We find that the precipitation of Si allows the early core to cool more 532 slowly than it would otherwise and supplies power to the geodynamo through- transition of the magma ocean should occur between 40% to 60% melt frac-537 tion (Abe, 1997;Solomatov, 2015), rather than at the intersect of the liquidus 538 or solidus, which would correspond to the occurrence of complete freezing and 539 first partial melt, respectively. Our thermal histories for the low oxygen con- The dataset using in this study is available to download from the sup-  Badro, J., Aubert, J., Hirose, K., Nomura, R., Blanchard, I., Borensztajn, core and its potential to drive an early exsolution geodynamo. Geophysical