Examining the power supplied to Earth ’ s dynamo by magnesium precipitation and radiogenic heat production

We examine magnesium and potassium solubility in liquid Fe mixtures, representative of Earth ’ s core composition, in equilibrium with liquid silicate mixtures representative of an early magma ocean. Our study is based on the calculation of the chemical potentials of MgO and K 2 O in both phases, using density functional theory. For MgO, we also study stability against precipitation of the solid phase. We use thermal evolution models of the core and mantle to assess whether either radiogenic heating from 40 K decay or Mg precipitation from the liquid core can resolve the new core paradox by powering the geodynamo prior to inner core formation. Our results for K show that concentrations in the core are likely to be small and the effect of 40 K decay on the thermal evolution of the core is minimal, making it incapable of sustaining the early geodynamo alone. Our results also predict small concentrations of Mg in the core which might be sufficient to power the geodynamo prior to inner core formation, depending on the process by which it is transported across the core mantle boundary.


Introduction
Classical core evolution studies based on cooling and inner-core growth with high thermal conductivity (κ) values point to a potential power shortage for maintaining the Earth's magnetic field prior to innercore formation, believed to be around 0.5-1 Gyr ago (Davies, 2015;Labrosse, 2015;Nimmo, 2015a), which is inconsistent with paleomagnetic determinations of a field back to at least 3.5 Ga (e.g.Biggin et al., 2011) and potentially 4 Ga (Tarduno et al., 2015(Tarduno et al., , 2020;;Fu et al., 2021;Bono et al., 2022).This incongruity in the high κ scenario defines the new core paradox.The exact value of the core's thermal conductivity is debated.Original estimates suggested a low value (<30 Wm − 1 K − 1 , Stacey and Anderson, 2001;Stacey and Loper, 2007) with some experimental studies continuing to support this (Konôpková et al., 2016;Hasegawa et al., 2019;Saha et al., 2020).However, over the last decade a growing majority of theoretical (Pozzo et al., 2012;de Koker et al., 2012;Pozzo et al., 2013Pozzo et al., , 2014;;Pourovskii et al., 2020;Zhang et al., 2022;Pozzo et al., 2022) and experimental studies (Gomi et al., 2013;Ohta et al., 2016;Inoue et al., 2020) suggest that values of κ for the Earth's core are much higher (70-140 Wm − 1 K − 1 ).Significant attention has been focused on identifying alternative power sources that can help to sustain the geodynamo prior to inner-core formation.In this paper we consider two of the prime candidates: heat released from the decay of 40 K, and precipitation of magnesium.In both cases the key challenge is to determine elemental partitioning behaviour at conditions up to and above core-mantle boundary (CMB) pressure P ~ 135 GPa and temperature T ~ 4000 K (where shallower conditions are relevant to equilibration during core formation).
Radiogenic heating provides power for magnetic field generation, though it is thermodynamically inefficient compared to the release of latent heat and light elements that accompany inner core growth because the heat is released throughout the core (Nimmo, 2015a).Indeed the more significant effect of radiogenic heating is to reduce the core cooling rate for the same CMB heat flow, which slows inner core growth.The effect depends on the nature and abundance of radiogenic elements in the core.Both uranium and thorium have been proposed to enter the core during its formation (Wohlers and Wood, 2015;Chidester et al., 2022), but most thermal history studies have focused on potassium (Nimmo et al., 2004).Early experimental investigations at relatively low P (~1 − 24 GPa) and moderate T (~2000 K) found that up to a few hundred ppm K could enter the core during its formation (Gessmann and Wood, 2002;Murthy et al., 2003) depending on the abundance of O and S in the metal (Bouhifd et al., 2007).However, laser heated diamond anvil cell experiments (Hirao et al., 2006;Watanabe et al., 2014;Blanchard et al., 2017;Chidester et al., 2022) and ab initio calculations (Xiong et al., 2018) on molten iron alloys at high P (>50 GPa) and T (>3500 K) suggested small concentrations of only 8-40 ppm.Core-mantle evolution models with low κ included ~ 400 to 800 ppm to help satisfy constraints on mantle cooling, inner core size and continued dynamo generation (Nimmo et al., 2004;Nakagawa and Tackley, 2010).Other models have argued that at least 250 ppm is required to match the present inner core size and maintain dynamo action with high κ (Driscoll and Bercovici, 2014).Additionally, including small 40 K concentrations (30 ppm) has been shown to make little difference to the predicted inner core age and ancient core temperature (Pozzo et al., 2022).Whilst these studies cannot be directly compared, it is clear that the effect of heating from 40 K can be significant and the concentration in the core is not agreed upon.
A second proposal to address the reduced power supply to the geodynamo is ascribed to light elements such as Mg and Si precipitating out of solution early in the history of the core, releasing power by leaving behind a heavy liquid that sinks and mixes the bulk core (O' Rourke and Stevenson, 2016;Badro et al., 2016;Mittal et al., 2020).We have recently re-examined the case of Si precipitation (Wilson et al., 2022) and so here we focus on Mg.The power provided by Mg precipitation depends on 1) the amount of Mg dissolved in the core during its formation, c Mg i (where c is mass fraction of solute); 2) the equilibrium concentration of Mg in the core, c Mg C , and; 3) the rate of Mg precipitation once the equilibrium concentration falls below the concentration of Mg initially dissolved in the core.
The initial core Mg concentration c Mg i is difficult to estimate because it depends on the manner in which the core formed.A recent review (Davies and Greenwood, 2023) used the range 0.3-3.6 wt%, where the lower estimates come from single-stage core formation models (O'Rourke and Stevenson, 2016;Helffrich et al., 2020) while the upper bounds were obtained from formation models that included a late high T event such as a giant impact (O'Rourke and Stevenson, 2016;Badro et al., 2016).The equilibrium Mg concentration c Mg C has been estimated by modelling high T experiments of partitioning between metal and silicate melts, which show a small pressure effect (Badro et al., 2016(Badro et al., , 2018;;Du et al., 2019) and so MgO precipitation is expected to occur first at the CMB where the core temperature is lowest.However, the uncertainties on both c Mg i and c Mg C mean that the onset time for MgO precipitation is poorly known.Indeed, using the c Mg C from previous studies (Badro et al., 2018;Du et al., 2019), Davies and Greenwood (2023) showed high c Mg i (~3.6 wt%) would have allowed precipitation for all temperatures below 6000 K, i.e. over most of Earth's history, while low c Mg i (~0.3 wt%) implies that precipitation does not occur for T above the present-day CMB temperature of ~4000 K (Davies et al., 2015), i.e.Mg has never precipitated from the core.(Badro et al., 2018;Du et al., 2019).A recent review obtained a precipitation rate in the range 0.3-1.5×10− 5 K − 1 (Davies and Greenwood, 2023) based on the aforementioned thermodynamic models and a range of plausible core and lower mantle chemistry.Compared to the case with no precipitation, the lower rate produced a minor change in inner core age and early core temperature, while the upper rate could double the predicted inner core age and reduce early core temperatures.
In this paper we present new ab initio determinations of MgO and K 2 O partitioning between liquid metal and both solid and liquid silicate at CMB conditions, complementing experimental studies that generally access lower PT.We employ our recently developed methodology for computing chemical potentials, which gives good agreement with extrapolations based on experimental determinations of FeO (Pozzo et al., 2019) and SiO 2 (Wilson et al., 2022) partitioning.We model Fe-rich metallic liquid alloyed with O and Si, since these lighter elements are generally predicted to be incorporated into the early core (Rubie et al., 2015;Badro et al., 2015) and can satisfy the present-day core mass and inner core boundary density jump (Davies et al., 2015).Reasonable compositions which are consistent with seismic observation contain up to 15 mol% O (Badro et al., 2015;Davies and Greenwood, 2023).We include a silicate melt as representative of the early mantle, when thermal history models predict CMB temperatures far above the pyrolite solidus (Nimmo, 2015b;Davies et al., 2015).We compare our results to literature data and incorporate them into core evolution models to predict the viability of dynamo action over geological time and to constrain the age of the inner core.

Methods
Partitioning of elements between the core and mantle is represented here by the partition coefficient of a species between silicate and ironrich liquid.These coefficients are calculated from the difference of chemical potential of the species in each system which in turn are evaluated via free energies in ab initio molecular dynamic calculations.In this section we describe the theory of chemical equilibrium and the partition coefficients which form the basis of this work, as well as how these are calculated from first principles using chemical potentials.

Chemical equilibrium
Chemical equilibrium is reached when the chemical potentials μ of all species are equal in the liquid iron mixture and the liquid silicate.Experiments usually report the distribution of composite species in the silicate, such as FeO, SiO 2 , MgO and K 2 O, and so the relevant equations are, for example for K 2 O and MgO: where , … are the molar concentrations of elements i, j, … in the core (superscript C) and mantle (superscript M).For simplicity of notation in the following we will leave out explicitly writing the dependence of μ on p, T, x i , x j , ….If the composite species are dissolved into their respective systems, then Eq. 1 and 2 can be written in terms of the chemical potentials of the single elements: This is indeed the case for the liquid iron mixture, and to some extent also for the liquid silicate where individual elements can be present in multiple species.
To obtain the relation that governs partitions it is useful to re-write the chemical potential by separating the configurational part (leaving just the excess chemical potential μ) and so Eq. 1 for K 2 O becomes and similarly for MgO.Eq. 5 can be rearranged in terms of the partition coefficients: where μK2O = 2μ K + μO for the liquid metal and liquid silicate.For MgO we obtain Complete details of the approach we outline here are given by Pozzo et al. (2019).μ is calculated via several different computational methods, all based on thermodynamic integration.Here, we use two of these approaches, referring to them henceforth as Method 1 and Method 2. In Method 1, a system A is slowly transformed (that is, allowing the system to remain in thermodynamic equilibrium) into system B, and the reversible work performed in this alchemical transmutation is equal to the free energy difference between B and A. This transformation involves changing the number of solute units in the system, meaning the change in free energy is equal to μ of the solute.Method 2, also described in Pozzo et al. (2019), is to refer to an external potential of known free energy, both for system A and system B. The transmutations from the external potential to the ab-initio potential then gives access to the total free energies of A and B, and from their difference one can obtain once again the chemical potential of interest.These two approaches are completely independent from one another, and by applying them both we can double check the internal consistency of our results and quantify uncertainty of the overall method.
We also calculate K D Mg between the liquid core and the solid B1 structure of MgO at the centre of the Earth.This provides an additional test of our method and helps to determine the process of Mg exsolution from the liquid core.This requires a slightly different approach to the other solutes in the core and an adjustment to the methods laid out by Pozzo et al. (2019).We used Method 1, in which the reference potential is the harmonic system, obtained by expanding the DFT potential energy function as a function of atomic displacements from their equilibrium zero temperature positions, and including only the quadratic term in the expansion: where U h (R) is the total harmonic energy function of the system which depends on the positions (R = r 1 , …r N ) of all the atoms.u i = r i − r i 0 is the displacement of atom i from its zero temperature equilibrium position r i 0 , U 0 the value of the potential with zero displacements, and i ∂r 0 j the force constant matrix, with the derivatives calculated at the equilibrium positions.The force constant matrix is computed using the small displacement method, as implemented in the phon code (Alfè, 2009).
The free energy per formula unit of MgO of the harmonic system U ref is obtained by summing the contributions of each normal mode ω q, s : where ℏ is the reduced Plank's constant, ∑ s runs over the 6 phonon branches (3 acoustic and 3 optical in the MgO crystal), and 1

∑
q is used to approximate the integral over all wavevectors q in the Brillouin Zone.The latter sum usually converges very quickly w.r.t. the number of q included and it is straightforward to compute once the force constant matrix is known.In fact, since the chemical potential of MgO in the liquid core is calculated by assuming that the atomic nuclei behave classically, we use the classical approximation for the harmonic free energy, given by: We note that at the conditions of interest the difference between the classical approximation and the full quantum free energy is only ~ 4-5 meV/formula unit, which is negligible for all practical purposes.By using the classical approximation for the solid as well as the liquid errors are minimised.

Ab initio simulations
The calculations are based on density functional theory (Hohenberg and Kohn, 1964;Kohn and Sham, 1965), using the VASP code (Kresse and Furthmüller, 1996), with the projector augmented wave (PAW) method (Blöchl, 1994;Kresse and Joubert, 1999) and the generalised gradient corrected functional known as PW91 (Wang and Perdew, 1991).Single particle wavefunctions were expanded in plane waves, with an energy cutoff of 500 eV.The electronic configuration of the various elements and the core radii are detailed in Table 1.
Electronic levels were occupied according to Fermi-Dirac statistics, with an electronic temperature equal to the ionic temperature.An efficient extrapolation of the charge density was used to speed up the ab initio molecular dynamics (AIMD) simulations (Alfè, 1999), which were performed by sampling the Brillouin Zone (BZ) with the Γ point only.The temperature was controlled with a Nosé thermostat (Nosé, 1984) and the time step was set to 1 fs.The simulation cells contained between 148 and 160 atoms in total, depending on composition.

Results
We calculate the chemical potentials of MgO and K 2 O in the magma ocean and in the liquid core.For convenience in our calculations, we set the composition of the mixtures at the outset, and then subtract a number of molecules dN.As a result, the compositions of the liquids are slightly different in the various cases, which could result in small differences in the chemical potentials if they depend on concentration.However, within the statistical accuracy of our calculations we cannot detect any such dependency.These chemical potentials are used to calculate the partition coefficients of Mg and K at the CMB, which will set the equilibrium composition of the core at the CMB.
Our calculations are carried out at pressures relevant to the CMB (124 GPa) as well as at mid-mantle pressures to examine the effect of pressure on K D .We run the majority our simulations at 5500 K to emulate the conditions of the hot early core, as these are more important for the thermo-chemical evolution of the core than the lower temperatures of the present day CMB (Davies and Greenwood, 2023).K D at these temperatures can be extrapolated to lower temperatures using the heat of reaction, also calculated here.We also include a low temperature (3600 K) result to better examine the T dependence of K D .We study a silicate composition which is close to pyrolitic (43.75 mol% MgO, 6.25 mol% FeO, 50 mol% SiO 2 ) and an additional case representing a more reduced case (55 mol% MgO, 16 mol% FeO, 29 mol% SiO 2 ).The metal compositions are chosen to explore reasonable O (4-16 mol%) and Si (0-8 mol%) concentrations, and an extreme case to better understand the role of O in the metal.

Magnesium
In Table 2 we report excess chemical potential differences δμ MgO between the core and the silicate mixtures containing Mg, as well as the resulting partition coefficients K D Mg , and the differences in the heat of reactions δH MgO .We complete some calculations using both Method 1 and Method 2 described in section 2. The two approaches show good internal consistency, and therefore for each species we take the weighted averages as our final results, with weights given by the inverse of the squares of the standard deviations.Low values of K D Mg for MgO are consistent with previous works (Wahl and Militzer, 2015) and imply low solubility of Mg in iron rich alloys, consistent with experimental studies (e.g.Badro et al., 2016Badro et al., , 2018;;Chidester et al., 2017;Du et al., 2017;Jackson et al., 2018).Details of separate metal and silicate calculations can be found in the supplementary information.
We compare predicted magnesium partitioning at the CMB from our results with the experimental data.The transfer of Mg can be represented using two possible reactions that can represent the transfer of Mg are dissociation (dc) and dissolution (dl).These reactions are written respectively as The equations determining the partition coefficients K D i for reaction i are Here x i is the molar concentration of species i, γ i the activity coefficient, and a, b and c are coefficients that are fixed by fitting to experimental data.γ i account for compositional variation in K D meaning that for the ideal case, γ i = 1.Note that the activity coefficients in the dissolution reaction arise because it is assumed that dissolved MgO further breaks down into Mg and O (Badro et al., 2018), while silicate activities are set to 1. Figure 1 compares the Mg partition coefficients in Table 2 to literature data, where because MgO is expected to break down to ionic species in the dissolution case, our K D dl and K D ds are equal.We use the dataset and thermodynamic modelling approach from Badro et al. (2018).Briefly, the model uses the interaction parameter formulation of Ma (2001) to represent the compositional dependence of K D and considers interactions between Fe, O, Si, Mg, C, and S. In the figure, red points show the data assuming no compositional dependence (activity coefficients set to 1), white points show K D with γ i ∕ = 0, and squares show our data (lines show K D Mg projected by ±500 K based on the local derivatives of the chemical potential).We find a strong T and weak P dependence of K D , in agreement with previous studies (e.g.Fischer et al., 2015;Badro et al., 2018) 2019) is around 2 log units, similar to our findings.We find P dependence to be small, which is consistent with past studies (e.g.Badro et al., 2018;Du et al., 2019).
By comparing fits of partitioning data to the predictions from Eq.s 13-14 it is hopefully possible to elucidate a single underlying reaction and use this to estimate the equilibrium Mg concentration at the CMB.Badro et al. (2018) found that these two reactions gave comparable fits to their dataset, but favoured the dissolution model as it produced less scatter at high Mg and O metal concentrations.Badro et al. (2018) also examined the exchange reaction but found it inferior to the other reactions, so we do not consider it here.For the dissociation reaction our K D Mg values bracket the experimental dataset at high T, while for the dissolution reaction our x O C < 30 mol. % result lies below the experimental range at high T.We therefore favour the dissociation reaction over dissolution.
We also calculate K D Mg for MgO between the liquid core and solid B1 MgO to further test our ab initio methods and provide a more complete picture of how Mg is exsolved from the liquid core.Results, reported in Table 3, agree qualitatively with Wahl and Militzer (2015), showing that only a very small amount of B1 MgO would be stable in solution in the liquid core, against precipitation of its solid phase.K D is smaller here than for the lower P and T liquid silicate interaction.This is largely because of the limited configurational space in the solid B1 structure compared to the liquid silicate.The T dependence of solubility means that exsolution of Mg will occur at the coldest region of the core first, but these results also show that solid precipitate would not be stable deep in the core either.Precipitation of Mg must therefore occur at the CMB where metal-silicate interaction is present.

Potassium
Due to a lack of experimental partitioning data for K at high pressure and temperature, it is not possible to effectively compare different reaction types as we have done with Mg.We only consider the dissociation reaction, due to compatibility with the ionic nature of liquid metals, minimal assumptions compared to other reactions (exchange implies a certainty of FeO exchange coupled to K 2 O, excluding the possibility of other candidates) and successful implementation in other systems (e.g.Mg and Si; Wilson et al., 2022).Our calculations show small K D K for K 2 O at all conditions studied (see Table 2), consistent with previous theoretical and experimental studies (Xiong et al., 2018;Gessmann and Wood, 2002;Hirao et al.,

Table 2
Ab initio excess chemical potential differences δμ X (eV) between the metal and the silicate phases for various compositions, pressures and temperatures.Also reported are the difference in the heat of reactions δH X (eV) and the partition coefficients K D .2006; Bouhifd et al., 2007;Blanchard et al., 2017;Chidester et al., 2022).Fig. 2 compares K D K from experiment with our results and illustrates that solubility is not dependent on O concentration for reasonable values.Xiong et al. (2018) find that for a liquid outer core with 23 mol% O at 4000 K only 30 ppm of K would be soluble in the core, or perhaps as low as 1 ppm, amounting to a negligible radiogenic contribution to core power sources.We note that a different formulation of K D K is used in this study and our results would also predict 1-30 ppm equilibrium concentration of K in the core given this formulation.Instead we simply consider a dissociation reaction (Eq.15).Our calculations suggest a K D K = 1 × 10 − 5 − 6 × 10 − 3 and a significant oxygen dependence above 13 mol % (Fig. 2).For a primitive mantle with a K concentration of 240 ppm (McDonough and Sun, 1995), increasing core oxygen concentration from 3 to 30 mol% promotes K solubility by ~10× from 250 to 2200 ppm.Despite the elevated solubility, this demonstrates that even for all but unrealistic compositions (30 mol% O) the maximum concentration of K in the core is small.

Discussion
To examine the effect of MgO and K 2 O on the power available to the ancient geodynamo in a moderately high k scenario (70 W m − 1 K − 1 everywhere in the core), we simulate the thermal history of the deep Earth whilst varying these contributions.We use the thermodynamic model of Mg solubility from Badro et al. (2018) to define precipitation of light elements from the liquid core.The removal of these light elements at the CMB leaves an iron-rich density anomaly at the top of the liquid core which provides additional convective power.For K, we use our results to define the temperature dependence of solubility and set initial core compositions.The decay of 40 K then heats the core, helping to offset the greater conductive power loss in a high k scenario.

MgO
Thermal histories of the Earth's core are modelled using coupled 1D parameterisations of the core (Davies, 2015) and mantle (Driscoll and Bercovici, 2014).Following Gubbins et al. (2004), the power available to convection in the core can be evaluated (whilst ignoring small terms) through where Q cmb is the heat flux across the CMB, Q s is the power from the gradual loss of primoridial heat from the core, Q L is the power released as latent heat due to inner core growth, Q g is the gravitational power generated from the preferential partitioning of O into the lowermost liquid core upon freezing and A and B are integrals of known core properties.E j and E k are the entropy due to ohmic dissipation and the entropy produced by thermal conduction in the core respectively and all other entropy terms have the same notation as the equivalent power term.By evaluating these power sources in the core, and Q cmb due to the mantle, we can evalutate the core cooling rate, which allows calculation of E j .The core and mantle models are joined at the CMB where the core defines the CMB temperature and the mantle defines Q cmb .There are several ways in which the mantle model can be implemented to define Q cmb .One option would be to use simulations of mantle convection to define the heat flux.This approach has the advantage of including all of the relevant physics of the mantle but is prohibitively expensive and would not allow us to explore parameter space as we have done here.A second option is directly opposite to the first, one would simply parameterise Q cmb in the model (for example as was done by Davies et al., 2022).This minimises the number of uncertain parameters but neglects to include any of the physics of the mantle.The obvious third option is a compromise of these two extremes where Q cmb is parameterised indirectly via a parameterisation of mantle convection.This has the benefits of being computationally efficient enough to allow exploration of parameter space and also including some of the physics of

Table 3
Excess chemical potential differences δμ MgO (eV) for MgO between liquid iron with various compositions and the solid B1 structure.Calculations have been performed with Method 1 [Method 2] for the liquid, and with Method 1 for the solid, using the inverse power and the harmonic potential as reference systems, respectively.Also reported are the difference in the heat of reaction δH MgO (eV), mantle convection.We opt for the third choice in this study, which has been the subject of a significant body of work in the literature, because of these benefits.Driscoll and Bercovici (2014) present a mantle evolution model based on the classic boundary layer theory.The model assumes a single convective mode at all times, plate tectonics.We do not include melting, because the parameterisation is uncertain; but the power from melting is found to be small when evaluated (0.7 TW).Other complex components of mantle convection including dehydration stiffening (Korenaga, 2006) and heating due to plate bending (Conrad and Hager, 1999) are not included in this model.With the default values of Driscoll and Bercovici (2014) the model can employ a present day Urey ratio of 1 3 (as inferred for Earthsee papers by Korenaga (2008) and Jaupart et al. (2007)) and a classical exponent of β = 1 3 while avoiding the mantle thermal catastrophe with a high present Q cmb and 2 TW of radiogenic heat in the core.We employ these default values and all the other values (unless otherwise stated) in the Driscoll and Bercovici (2014) parameterisation of the mantle (Table 3) except we do not include any radiogenic heating of the core.The same is true of the core model, where values of the original study (Table 1 Davies, 2015) remain unchanged unless otherwise stated.
Additional entropy is produced in the liquid core from the precipitation of Mg at the CMB where, T cmb is T at the CMB, V c is the volume of the core, ψ is the gravitational potential, ρ is density, Mg dT cmb is the precipitation rate of Mg (in wt% K − 1 ), and chemical expansivity is α ppt Mg = 1.12 (O'Rourke and Stevenson, 2016).The power associated with precipitation of light elements at the CMB is then Precipitation is quantified by removing Mg from the core until the equilibrium concentration at the CMB is achieved.This assumes that the core is thoroughly mixed on timescales far shorter than the cooling rate and that precipitation occurs at the coolest part of the core.We treat the mantle which interacts with the core as a constant composition, equal to the bulk composition, meaning that the mantle similarly sweeps precipitate from the CMB on timescales shorter than the cooling rate and that the volume of the mantle is sufficiently large that the precipitate will make little difference to the bulk composition.
To define the equilibrium concentration of Mg in the core, we use the interaction parameter model of Ma (2001) with the parameters found by Badro et al. (2018), who also studied this problem but inferred thermal evolution outcomes from their own temperature dependence of solubility and the core cooling rate of O' Rourke et al. (2017).This model is based on liquid-liquid interaction between metal and silicate and defines both the composition of the liquid core and the precipitation rate C Mg .We set initial compositions for the core and mantle for each case and tune models to satisfy present day constraints of heat flux from the convecting mantle (38 ± 3 TW, Jaupart et al., 2007), mid-mantle temperature (2320 ± 50 K, Katsura et al., 2010) and inner core radius (1221 ±10 km) and maintaining positive entropy for powering the geodynamo during the last 3.5 Gyrs.Freezing of the inner core is controlled by the melting curve of pure iron from Alfè et al. (2002), 6125 K at 360 GPa, with a linear correction to the chemical potentials of solid and liquid due to O content (according to Alfè et al., 2002), following the method of Davies (2015).We vary O from 5 to 20 mol% giving present day inner core boundary temperatures ranging from 5870 K to 4950 K.We vary the initial CMB temperature (T cmb ) and ratio of upper to lower mantle viscosity (f viscosity ) where the cooling rate of the core is tuned by f viscosity and the initial temperature of the core is set by T cmb .
We consider two precipitation scenarios and three different initial compositions for each of the mantle and the core.The precipitation scenarios explore how Mg is removed from the core.In the first scenario, Mg and O are removed to the mantle in equal proportion as Mg is precipitated, ensuring charge balance and also reducing the melting point depression associated with O content of the core.For this case, in Eq. 17 α ppt becomes an average of α ppt Mg and α ppt O (1.1), and C Mg is doubled.
In the second scenario, only Mg is removed, accounting for the  Gessmann and Wood, 2002;Hirao et al., 2006;Bouhifd et al., 2007;Blanchard et al., 2017;Chidester et al., 2022).Experimental K D K is calculated from element concentrations in recovered quenched samples (Eq.15).Black line is the best-fit of temperature dependence logK D K = A*T + B. Right: compositional dependence of K solubility in the core, where K is assumed to enter the core through dissociation.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) A.J. Wilson et al. possibility that other reactions may account for charge balancing.Our calculations do not define what is precipitating from the core, only the stable fraction of Mg in the liquid metal.These two scenarios explore the uncertainty of how Mg is transported across the CMB whilst also examining the extremes of E ppt .When removing Mg and O from the core through Mg precipitation, significantly more convective power is contributed to the core because far more mass is being transferred across the CMB.We do not consider any power arising from the enthalpy of fusion or dissolution for Mg or MgO after precipitation and release to the mantle because the mechanism of incorporation into the mantle is not understood.We expect this effect to be small due to the relatively low volumes being released (0.3 mol% over 4.5 Gyrs).Both cases see all terms in Eq. 17 evaluated according to the present day core composition, which evolves due to precipitation and inner core growth.
We consider a mantle composition of 30, 50 (pyrolitic) and 70 mol% MgO, where higher concentrations allow more Mg to be dissolved into the core (Eq.7).We also consider three initial oxygen concentrations of the core (5, 10 and 20 mol%).The model of Badro et al. (2018) predicts that Mg and O are mutually beneficial for solubility in the core, meaning a more oxidised core can host higher Mg concentrations and therefore provide more power from precipitation.The initial Mg concentration of the core is the equilibrium concentration at the CMB for the initial conditions for our cases.A maximum of 2 mol% Mg is dissolved into the core (for a 20 mol% O in the core and 70 mol% MgO in the mantle) and a minimum of 0.4 mol% (5 mol% O and 30 mol% MgO).These core and mantle compositions represent a reasonable range as lower concentrations result in limited Mg being dissolved into the core at early times, and therefore limited power from precipitation immediately before inner core nucleation.Higher O concentrations produce significant melting point depression of the iron alloy in the core, meaning that for the core to be cold enough to freeze the inner core, the power from both secular cooling and precipitation are greatly reduced.
We evolve these nine starting compositions and two precipitation scenarios both with and without the power from precipitation included (setting Q ppt ∕ = 0 and Q ppt = 0 respectively).This gives thirty six unique models for f viscosity and T cmb , with which to best fit the present day constraints.If a case is unable to maintain E j > 0 for all time after 3.5 Ga, it is possible to increase the available entropy, for example through a hotter initial T cmb , however this will result in a smaller than observed present day inner core.All cases presented here have been tuned to best match present day constraints and whilst some may fail by a small margin, this represents the most optimised case.Fig. 3 shows two example cases where an initial core O concentration of 10 mol% and pyrolitic mantle is evolved with just Mg being extracted and with Mg precipitation also removing O to maintain charge balance.These cases illustrate the additional convective power supplied by the precipitation of MgO over Mg.At the time when the inner core forms, precipitation of Mg is slow, meaning that for both cases there is little difference in the rate of change of O concentration because this is controlled predominately by inner core growth.The inner core age does not differ greatly between these two cases, highlighting that whilst Q ppt can be small, the effect of E ppt remains significant.
All models where Q ppt = 0 fail to maintain E j > 0 prior to inner core nucleation (O'Rourke et al., 2017;Driscoll and Davies, 2023).Similarly, all models where only Mg is extracted through precipitation (O concentration only changes due to inner core growth) also fail in this regard.All models where both Mg and O are removed from the liquid core with cooling are able to maintain surplus power for the geodynamo for all time after 3.5 Ga but those with an initial O concentration of 20 mol% struggle to reproduce the present day mantle temperature.Fig. 4 compares the age of the inner core and the CMB temperature at 3.5 Ga (indicative of early core conditions) for these nine successful cases out of the thirty six cases studied.The high temperatures of the early core are consistent with the temperatures our ab inito calculations were Fig. 3. Thermal histories including the precipitation of Mg from the core.Successful models are those which reproduce properties of the deep Earth (diamonds, present day constraints; surface heat flow, inner core radius and mid-mantle temperature) whilst consistently providing power (E j > 0) for the geodynamo (black dashed line, constraint for the past 3.5 Gyrs).Example cases are shown for a pyrolitic mantle and a 10 mol% initial oxygen content of the core.When O and Mg are extracted in equal proportion to the mantle via precipitation of Mg (solid lines), each contribute similar convective power to the outer core and E j > 0 for all time.When only Mg is removed (dotted lines) there is insufficient power for the geodynamo for 700 Myrs prior to inner core nucleation.Once the inner core forms precipitation rate is low in both cases and changes in oxygen concentration are primarily driven by inner core growth.
performed at and suggest that liquid-liquid interactions at the CMB were long lived.O content of the core provides a strong control on the core temperature due to melting point depression of the iron alloy; more O rich models must be cooler in order to freeze the inner core.When the mantle contains less MgO compared to pyrolite, or the O concentration of the core is lower, there is less Mg dissolved into the core and precipitation rates are lower.For E j to be greater than zero prior to inner core formation, these cases must extract slightly more heat through secular cooling meaning that the inner core is older by 20-50 Myrs.The opposite is true, to a lesser degree for mantle compositions with more MgO than pyrolite, producing elevated precipitation rates.These differences are relatively small, meaning that ancient core temperatures are similar.
We define the success of these models based on satisfying constraints of positive dynamo entropy, inner core size and mantle temperature and convective heat flux but note that additional properties of the deep Earth might be useful in assessing compatibility with the present day deep Earth.Q cmb can be used in this regard, but does not provide a direct constraint because it remains largely unknown.We find values of 6-8 TW in our successful models, which is consistent with our previous work (Wilson et al., 2022).Previous estimates of Q cmb commonly range between 5 and 15 (Jaupart et al., 2007;Lay et al., 2008;Nimmo, 2015c;Frost et al., 2022), with the upper range of estimates only being required for higher values of κ than applied here and without precipitation contributing to core convection.

Potassium
We apply the same model as detailed above for studying the effect of 40 K on the thermal evolution of the core, except without the precipitation of MgO and with the heat from 40 K decay.There are relatively few experimental metal-silicate partitioning studies which include K compared with those including Mg (due to ubiquity in silicates).This, combined with the low solubility of K in liquid Fe, makes the construction of a thermodynamic model (such as those by Badro et al., 2018;Wilson et al., 2022) challenging.Parameters defining the various possible interactions of potassium in iron rich liquids are not yet resolvable.As such we are unable, with presently available experimental data, to construct a thermodynamic partitioning model for K of any utility.Instead we define the temperature dependent solubility of K in the core completely with our own calculated K D K (Table 2).We find that decay rate exceeds the exsolution rate due to core cooling in all cases meaning that we do not include a temperature dependence of solubility and simply set an initial K concentration for the core.
In Fig. 5 we show three cases with different initial K concentrations; 0 ppm as a reference, 100 ppm as a low concentration, 250 ppm as the maximum soluble concentration allowed by our K D K with reasonable O concentrations (e.g.Badro et al., 2015;Davies et al., 2015).Power is integrated into our model assuming 2.6×10 − 15 W kg − 1 ppm − 1 (Clauser and Gupta, 2011).The temperature and E j of the core are extremely similar across all of these cases with even the highest concentration only contributing 1.3 TW at the present day.Decay of 40 K is, even in the most optimistic case, just short of providing the required power to sustain the geodynamo prior to inner core nucleation.This agrees with the findings of (Driscoll and Davies, 2023) who found that a present day radiogenic heat production of ≥ 2 TW and an f viscosity ≤ 10 (3-8 in our cases) is needed for E j > 0 prior to inner core formation and a present day inner core radius of 1221 km.

Conclusions
In this work we have calculated the solubility of MgO and K 2 O in liquid Fe mixtures representative of Earth's core composition and in equilibrium with both a solid oxide and liquid silicates representative of an early magma ocean, at times when the bottom of the mantle may have been completely molten.The methods used are similar to those we previously employed to investigate FeO solubility (see Pozzo et al., 2019), with solubility data here determined by computing the chemical potentials of MgO and K 2 O at high P, T and variable composition in liquid Fe mixtures and silicate solids and melts.Our results, in addition to previous implementations of these methods (Pozzo et al., 2019;Wilson et al., 2022), show that this approach of calculating K D from ab initio simulations is consistent with experimental studies of metalsilicate partitioning at high pressure and temperature.
By computing chemical potentials we have established that K D K for the core is small.These results are consistent with those of Xiong et al. (2018), although implementing different reactions mean a ~10×  A.J. Wilson et al.

Fig. 1 .
Fig. 1.Temperature dependence of MgO partition coefficients for dissociation (top) and dissolution (bottom) reactions.Our results explore four oxygen concentrations; 4 mol% (light blue), 13-15 mol% (light green), 30 mol% (dark green), and three pressures; 124 GPa (squares), 58 GPa (diamond), 50 GPa (triangle)).Lines show extrapolation of our results via δμ x δT = μ x − Hx T for ±500 K.The dataset and model of Badro et al. (2018) are used for experimental ideal γ i = 1 (red points) and non-ideal γ i ∕ = 1 (white points), where compositional effects are included.Fits of ideal and non-ideal K D are shown as read and black lines respectively.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Left: Potassium partition coefficient as a function of the inverse of temperature from this study (green and turquoise squares) compared with another ab initio study (pink squares; Xiong et al., 2018) and several experimental partitioning studies (circles;Gessmann and Wood, 2002;Hirao et al., 2006;Bouhifd et al., 2007;Blanchard et al., 2017;Chidester et al., 2022).Experimental K D K is calculated from element concentrations in recovered quenched samples (Eq.15).Black line is the

Fig. 4 .
Fig. 4. Inner core age and CMB temperature at 1 Ga for thermal histories including the precipitation of Mg from the liquid core.Initial oxygen content of 5, 10 and 20 mol% are shown as light, medium and dark colours respectively.Mantle compositions with 30 (greens), 50 (reds, pyrolitic) and 70 (blues) mol.% MgO are compared.Only successful cases are shown (solid colours), where Mg and O are extracted from the core in equal proportion (triangles) and produce E j > 0 prior to inner core formation.Also shown are 3 cases where only Mg is removed to the mantle (squares) which do not maintain a dynamo (transparent colours).

Fig. 5 .
Fig.5.Radiogenic power in the core from decay of 40 K for parameterised thermal evolution models.Thermal histories are shown with initial K concentrations of 0 (black), 100 (red) and 250 ppm (orange) which correspond to the equilibrium concentration at the CMB for a temperature of 5500 K, as an upper limit case.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Electronic configurations and core radii of the PAW potentials employed in this work.
Badro et al. (2018)019), log K D increases with increasing oxygen concentration in the metal.Our largest O concentration is 30 mol% while our smallest O concentration is 4 mol%, which encompass the values in theBadro et al. (2018)database.The variation in K D over this range of O fromDu et al.  (