1 Introduction

How our Earth formed is a question of major importance in contemporary astrophysics, cosmochemistry, and geochemistry, not only to explain the origin of our own planet but also to understand the prevalence of potentially habitable planets orbiting ot...

2 The ADAP code

We present here the ADAP code (Accretion and Differentiation of Asteroids and Planets). The code solves the 1D thermal conduction equation on a spherically symmetric grid. Importantly, the code includes the time-dependent melting of the layers and the sinking of denser material through lighter material, while conserving the total energy. The code considers three heat sources: (1) decay of short-lived radionuclides, (2) accretion and (3) differentiation.

2.1 Basic materials

ADAP bundles, for simplicity, the interior composition of asteroids and planets into five basic materials: Here, Met represents metal (consisting of iron Fe, iron-sulfide FeS and nickel Ni) that will enter the core upon planetary melting due to their ...

2.2 Material transitions

We show an overview of how ADAP treats materials and their transitions in Fig. 1. The three fundamental species (Met, Sil and Wat) combine into the primitive starting material MetSilWat. We assume that this material melts in several distinct stages. At...

2.3 Material densities

For the five basic material densities we choose: Here the low density of Met reflects its high contents of sulfur (Morard et al. 2018; Johnson et al. 2020) and the density of Sil is based on the uncompressed density of iron-depleted mantle material (Elkins-Tanton & Seager 2008). We assume for simplicity that the densities are constant irrespective of temperature and pressure. The densities of composite layers made out of these five base materials are constructed using volume partitioning, ρ=1∑ fi/ρi,(3)where fi is the normalized mass fraction, with ∑i fi ≡ 1 for each layer. Clay contains 15% by mass; we include here the full mass of the H2 molecule and ignore any outgassing of hydrogen in the clay formation process. The formation of oxidized metal MeO, on the other hand, is associated with hydrogen loss. We do not keep track of this lost hydrogen since clay formation happens in the earliest stages of planetary growth where the protoplanet cannot yet hold onto an outgassed atmosphere.

2.4 Differentiation

We require that the material density decreases outwards from the centre for stability. If a denser layer lies on top of a lighter layer, then we swap the two layers. We calculate the new gravitational binding energy after the swap and add the binding energy difference to the thermal energy of the two layers. We require the two layers to maintain a shared temperature after the addition of the thermal energy. This approach conserves total energy (gravitational plus thermal) in the differentiation process.

2.5 Heat transport

We solve for the total energy e present in each spherical shell. We choose to have e as the evolved variable because this eases the division of the total energy between latent heat and temperature (see Sect. 2.6). We numerically solve the heat transfer equation ℱ=−K∂T∂r.(4)Here ℱ is the heat flux and K is the heat conductivity. We follow Desch et al. (2009) and discretize the heat transport across the interface between shell i and shell i + 1 in an upwinded fashion, ℱi=−0.5(Ki+Ki+1)Ti+1−Tiδr.(5)Here δr is the radial grid size. We consider the mean of the heat conductivities across the interface to avoid building up spuriously large temperature gradients at interfaces between liquid and solid layers. The evolution equation for the energy of cell i, ei, becomes e˙i=4πri−12ℱi−1−4πri2ℱi.(6)

2.6 Heat capacities and latent heat

We translate the energy in a cell ei to temperature Ti taking into account both the heat capacity and the latent heats. Heat capacities and latent heats are assumed to be constant, independent of pressure and temperature. This way the translation from energy to temperature can be calculated quickly and precisely using analytical formulae. We use the NIST Chemistry Handbook (Linstrom & Mallard 2021) and Lange & Ahrens (1982) to identify approximate values for the heat capacity relevant at high temperatures: We ignore the increase in cp as the temperature reaches the critical temperature for the material. The combined heat capacities of composite materials are calculated from cp=∑ifici(7)where ƒi are the mass fractions and ƒi = 1.The latent heats of melting of the basic materials are assumed to take the constant values (Sahijpal et al. 2007) We ignore heat consumption in clay desiccation and set HCla = 0. The latent heats of the composite materials are calculated from L=fjHj(8)where j = 1, 2, 3 represents Met, Sil and Wat and ƒj their total mass fractions in a composite material. In this approach, the O in MeO is inherited from oxidation of Fe and Ni by water and hence the translation of the energy of a liquid MeO layer to temperature must, for consistency, take into account the latent heat of both iron and water. Similar considerations apply to the Cla material component.

2.7 Thermal expansion coefficient

The coefficient of thermal expansion is defined as α=1V(∂V∂T)P=1V∑i(∂Vi∂T)P=∑iViVαi,(9)where V denotes the volume. The coefficient enters calculations of the adiabat and the convective heat conduction. For the basic materials we take pressure and temperature-dependent parametrizations for αMet from Chen et al. (2007), αSil from Abe (1997) and αWat for ice from Desch et al. (2009) and for water from the IAPWS database (see Johansen et al. 2023a, hereafter Paper II).

2.8 Conduction and convection

Equation (4) describes heat transport by conduction with conductivity K. We take approximate heat conductivity values from Desch et al. (2009) and Neumann et al. (2012): The heat conductivity of mixed layers is calculated using the geometric mean model (Neum...

2.9 Radioactive heating

We add to the energy equation (Eq. (6)) also the heating by decay of 26Al, e˙26,i=MSillE26A26exp[ −t/τ26 ].(18)Here E26 = 3.12 MeV is the energy released by each decay minus the energy carried by neutrinos (Castillo-Rogez et al. 2009), A26 = fAlf26/(m26τ26) is the decay rate per kg of Sil material at t = 0, m26 is the mass of the 26Al atom and τ26 = 1.05 × 106 yr is the decay constant. We have normalized A26 to the Sil component only. Thus fAl is the fraction of aluminium in the silicate part the solar composition; we calculate fAl = 0.022 by normalizing first relative to hydrogen and helium and then relative to the total silicate contents (Lodders 2003). We take a base value of f26 = 5 × 10−5, but note that the initial amount of 26Al may have been lower in the terrestrial planet forming region than in the CAI forming regions closer to the star (Larsen et al. 2011; Connelly et al. 2012; Schiller et al. 2015). This gives an initial heating rate of 5 × 10−7 Wkg−1 for the silicates. By construction, the Met and Wat components have no radiogenic heating. Normalized to MetSil the heating rate is lower, 3 × 10−7 Wkg−1.This is about 1.5 times the value considered by Hevey & Sanders (2006) for CI chondrites. However, the CI chondrites have a density of only 2.2 × 103 kg m−3 and contain a significant amount of oxidized metal and clay.

2.10 Accretion of mass and energy

ADAP allows the planetesimal to grow in mass by accretion and includes the accretional heating of the surface. The mass accreted in each time-step is accumulated until it reaches a high enough value to create an additional grid shell on top of the existing planetesimal. This new cell inherits the temperature of the underlying grid cell. We add in each time-step accretion heat to the outermost cell, e˙acc=GMM˙R.(19)Here G is the gravitational constant, M is the mass of the protoplanet, M is the growth rate and R is the radius of the protoplanet. We nevertheless ignore the accretion phase here in Paper I and we checked that indeed the accretion energy is negligible for the small planetesimals considered here. The accretion heat does become important for massive, accreting protoplanets when released under the thermal blanketing by a dense outgassed atmosphere (Matsui & Abe 1986a,b). We discuss this further in Paper II.

2.11 Removal of outgassed atmosphere

Planetesimals and terrestrial planets forming by pebble accretion grow and evolve within the protoplanetary gas disc and hence the outgassed volatiles can escape by gas drag for low protoplanet masses. The gas flow around the protoplanet is dominated by...

3 Heating and differentiation of planetesimals

We consider Vesta as a prototype of the population of protoplanets (Russell et al. 2012) that formed in the terrestrial planet region and competed for pebbles to grow to planetary sizes. Vesta is differentiated and has an iron core of approximately 110 km in radius, giving a core mass fraction of approximately 18% (Russell et al. 2012). Its thermal evolution may be connected to that of Ceres, which remained volatile-rich and avoided core formation, with the main difference being the formation time and hence the amount of radiogenic heating available (McCord & Sotin 2005). The mantle of Vesta contains a large fraction of oxidized iron FeO, approximately 24% inferred from HED meteorites, observations with the Dawn spacecraft and comparisons with the H-chondrites (Toplis et al. 2013; Trønnes et al. 2019). We take the high oxidation as evidence that Vesta formed exterior of the ice line with a substantial water fraction. This ice then formed clays and oxidized metallic iron when it melted, leaving the core to form mainly from FeS, which is siderophile and provides the main sulfur carrier at high temperatures (Scott et al. 2002). The clays would have subsequently dried out after heating above 900 K, leaving water only as structurally bound remnants in apatite minerals (Sarafian et al. 2014).

3.1 Planetesimal model setup

To facilitate comparisons with Vesta data and to differences between Vesta and Ceres, we fix the planetesimal radius to R = 250 km and the background temperature of the protoplanetary disc is set to T0 = 135 K. The radial resolution element is chosen to be ...

3.2 Interior temperature and composition of planetesimals

We show space-time plots of the results of our heating and differentiation model in Fig. 3 for two different starting times. Our planetesimals start with a significant fraction of water mixed with metal and silicates (MetSilWat). As the water melts at...

3.3 Varying the 26AI abundance

In Fig. 4 we show the maximum temperature of the planetesi-mal as a function of the starting ratio of 26Al/27Al. The starting amount of 26Al plays a decisive role in the final composition of the planetesimal. Planetesimals that form within two half-lives of 26Al undergo the sequence of melting and differentiation that happened to Vesta; the water is then only present in the form of FeO in the final body. Forming in the interval between two and three half-lives instead leads to the formation of a desiccated body that lost its water but did not heat enough to differentiate into a silicate mantle and a metal core. Formation in the interval between three and four lifetimes gives maximum temperatures below the desiccation temperature of clay. These Ceres analogues are dominated in their interiors by clay and the mantle consists of the excess water that could not enter the phyllosili-cates. Finally, primitive bodies akin to comets and Kuiper belt objects form after four half-lives of26Al.

3.4 Mapping the thermal evolution of planetesimals

The maximum interior temperature and the initial composition can be mapped onto the thermal evolution of the inner regions of the solar protoplanetary disc. We adopt the protoplanetary model of Johansen et al. (2021), which considers the irradiation from the forming star as well as viscous heating close to the star where the ionization degree is high enough to trigger the magnetorota-tional instability (Desch & Turner 2015). We take into account here that the initial 26Al/27A1 ratio in the formation region of asteroids in the Solar System appears to have been lower outside of the CAI-forming region (Larsen et al. 2011; Connelly et al. 2012); we thus use 1/4 of the CAI value to be consistent with these measurements (Schiller et al. 2015). This yields the planetesimal map displayed in Fig. 5. Exterior of the water ice line, FeO-rich planetesimals form by melting of the ice, unless the planetesimal forms after 1.8 Myr or so. Between the silicate sublimation front and the water ice line, we have FeO-poor bodies akin to enstatite chondrites or aubrites (Keil 2010).

4 Core mass fraction and FeO fraction in the mantle

We will present full models of heating and differentiation of terrestrial planets forming by pebble accretion in Paper II and Paper III. Here we show an important intermediate result, namely that the core mass fraction and FeO mantle fraction of planets can be generalized directly from the planetesimal heating model presented in the previous section.

4.1 Core and mantle iron in Solar System bodies

The distribution of iron between metal core and mantle FeO within a differentiated body carries a history of the amount of water that was accreted. The fraction of FeO in the mantle is known for Venus, Earth, Mars and Vesta. The FeO fraction increases wi...

4.2 Light elements in the core

When discussing core mass fractions, the possible presence of lighter elements in the core becomes important to take into consideration. The density of the Earth’s core is consistent with incorporation of both Si and O into the metallic melt during differentiation (Badro et al. 2015). The cores of Vesta and Mars could both contain approximately 15% weight sulfur (Rivoldini et al. 2011). The sulfur weight fraction can be maximally 36% if the core consists entirely of FeS. Using the solar composition (Lodders 2003) we have S/Fe = 0.53 in number, which would yield a core mass fraction of at least 27% if all S is transported to the core in the form of siderophile FeS. This is clearly above the estimate for Vesta and Mars. Sulfur could nevertheless have entered the core and expelled later due to its incompatibility with solid iron (Johnson et al. 2020). Steenstra et al. (2019) infer 15% sulfur in the core of Vesta from chalcophile elements. This is consistent with measurements of the original S contents of iron meteorites (Chabot 2004), later expelled due to incompatibility with the solid iron phase. We therefore consider 53% of the iron to have been originally in the form of FeS. This attachment is important, since the iron is protected from oxidation by the sulfur. We additionally consider 5% Ni in the cores and furthermore 10% light elements by mass for Venus, Earth and Vesta. For Mars, based on seismic measurements by the Insight mission, we take 5% Ni and 15% light elements in the core (Stähler et al. 2021).

4.3 The sulfur content of chondrites

The CI chondrites contain an approximately solar composition of S, with S/Fe = 0.5 - all other meteorite classes are depleted relative to this value (Ebel 2011). FeS sublimates at approximately 700 K under nominal protoplanetary disc abundances and pressures (Lodders 2003). This temperature makes S volatile in the envelopes of the more massive planets. S should therefore strongly record the temperature of the envelope of the parent body, with asteroids accreting nominal values and planets successively missing more S with increasing mass. Earth’s mantle contains only approximately 200 ppm of sulfur (Jackson et al. 2021), while the core of Earth may hold significant S due to metal-silicate partition coefficients in the range 100-1000 for the relevant pressures at the core-mantle boundary (Boujibar et al. 2014; Jackson et al. 2021). Chondrules were already depleted in sulfur (Marrocchi & Libourel 2013) and hence accretion of thermally processed chondrules constitutes an additional pathway to S depletion before accretion.

4.4 Planet formation model

We use the planet formation model of Johansen et al. (2021) to test if the pattern evident in Fig. 6 is consistent with oxidation of iron by interaction with water outside of the water ice line. This model considers the formation and migration of terrestrial planets that accrete both pebbles and planetesimals. The plan-etesimals are assumed to form within a ring of width 0.1 AU at a distance of 2.3 AU from the host star, hence the protoplanet can only accrete planetesimals until it migrates out of the ring. The pebbles are assumed to drift through the protoplanetary disc. The planetesimal contribution is high (~50%) as a protoplanet grows up to the mass of Mars; for higher masses pebble accretion strongly dominates.We calculate for each time-step in the planet formation model of Johansen et al. (2021) the accreted iron mass, with the iron fraction based on the total amount calculated per body in Fig. 6 to take into account a varying iron content in the accreted material for the different bodies. We assume that the iron in the accreted pebbles has a low oxide fraction (FeO, 5%) and high sulfide fraction (FeS, 45%), broadly consistent with spectral analysis of interstellar dust (Westphal et al. 2019)1. We assume that the iron accreted together with water becomes 100% oxidized. We furthermore assume that iron in the accreted planetesimals is 100% oxidized, except for the iron bound in FeS. All the protoplanets are started at r0 = 1.6 AU.

4.5 Oxidation state of Vesta

The oxidation of iron by dissolution in hot water through the Schikorr reaction (Eq. (1)) should yield an average iron oxidation similar to magnetite (Fe3O4), containing two Fe3+ and one Fe2+ to balance the four O2−. This would imply a strongly oxidizing ma...

4.6 Calculated FeO fractions

In Fig. 7 we show the calculated FeO mantle fraction and core mass fraction from the planet formation model of Johansen et al. (2021). We show in the figure both the total FeO mantle fraction in the model and that of planetesimals separately, to illustrate th...

4.7 Varying the pebble composition

The final FeO mass fraction depends sensitively on our assumption that the accreted pebbles were very reduced. This assumption is partially balanced by the oxidation of this iron by dissolution in hot water; particularly the early-formed planetesimals de...

5 Discussion and implications

We demonstrate in this paper how a high primordial ice content of planetesimals is remembered even after the bulk of the water has been expelled, through the high fraction of oxidized iron (FeO) in the mantle. Lichtenberg et al. (2019) proposed that th...

Acknowledgements

We thank an anonymous referee for carefully reading the three papers in this series and for giving us many comments and questions that helped improve the original manuscripts. We are also grateful to the second referee for reading the revised manuscript carefully. A.J. acknowledges funding from the European Research Foundation (ERC Consolidator Grant 724687-PLANETESYS), the Knut and Alice Wallenberg Foundation (Wallenberg Scholar Grant 2019.0442), the Swedish Research Council (Project Grant 2018-04867), the Danish National Research Foundation (DNRF Chair Grant DNRF159) and the Göran Gustafsson Foundation. M.B. acknowledges funding from the Carlsberg Foundation (CF18_1105) and the European Research Council (ERC Advanced Grant 833275-DEEPTIME). M.S. acknowledges funding from Villum Fonden (grant number #00025333) and the Carlsberg Foundation (grant number CF20-0209). The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC), partially funded by the Swedish Research Council through grant agreement no. 2020/5-387.

References

All Figures