Stochastic accretion of the Earth.

Chondritic meteorites are thought to be representative of the material that formed the Earth. However, the Earth is depleted in volatile elements in a manner unlike that in any chondrite, and yet these elements retain chondritic isotope ratios. Here we use N-body simulations to show that the Earth did not form from chondrites, but rather by stochastic accretion of many precursor bodies whose variable compositions reflect the temperatures at which they formed. Earth's composition is reproduced when initial temperatures of planetesimal- to embryo-sized bodies are set by disk accretion rates of (1.08±0.17)×10-7 solar masses/yr, although they may be perturbed by 26Al heating on bodies formed at different times. Our model implies that a heliocentric gradient in composition was present in the protoplanetary disc and that planetesimals formed rapidly within ~1 Myr, in accord with radiometric volatile depletion ages of the Earth.

Introduction

The terrestrial planets of our Solar System, including Earth, are thought to represent the sum of progressive collisional growth of smaller bodies, whose size distribution evolved through time and space. In the classical model, collisions involve mostly nearby objects and lead to the formation of tens of Moon- to Mars-sized embryos within the lifetime of the protoplanetary disk (a few Myr)[1]. Only after the disappearance of gas from the disk, which causes embryos to become dynamically unstable, is mixing between objects from different heliocentric distances expected over 107 to 108 yr timescales[2,3]. Alternatively, recent models suggest that the accretion of most of the Earth’s mass occurred within the lifetime of the disk via the capture of sunward-drifting cm-sized pebbles that initially formed at the ‘snow-line’[4,5].

Such models of Earth’s accretion may be tested by examining the chemical and isotopic composition of the Bulk Silicate Earth (BSE), which has been estimated from geological evidence[6]. The BSE is the Earth excluding its metallic core, which is inaccessible to detailed chemical scrutiny. In the BSE, refractory lithophile elements (RLEs[7]) occur in solar proportions relative to one another, whereas the abundances of moderately volatile lithophile elements decrease according to their calculated 50 % nebular condensation temperatures (T) (Fig. 1). These elements have T < 1300 K and include the lithophile elements K, Rb, Zn and Cl, whose isotopic constitutions are similar to those measured in chondrites[8-11]. This observation precludes their depletion as being the result of partial evaporative loss from the Earth or its precursors, which would have caused detectable isotopic fractionation[12,13].

Figure 1

Fits to the bulk compositions of the mantles of Earth and Vesta.

CI-, Al-normalised abundances (f) and 1σ (1 standard deviation) uncertainties of lithophile elements in Bulk Silicate Earth[6] (green) and Vesta[18] (yellow) plotted as a function of their 50 % nebular condensation temperatures[19], with Si, Na, K, Zn and Cl updated from ref. 20. Both bodies have abundances that approximate a cumulative normal distribution (eq. 1). The corresponding probability distribution function (PDF) is shown in the top panel. The curves are determined by a least-squares fit to error-weighted elemental abundances, where is the modelled value, f is the abundance, and s its 1σ uncertainty (Supplementary Table 2). Best-fit values of T and σ are listed on the plot and yield Root-Mean Squared deviations of 0.071 (n = 35) for Earth and 0.028 (n = 25) for Vesta. Manganese falls significantly below the trend for Earth’s mantle, likely owing to partial incorporation into Earth’s core[21]. The linear scale highlights gradual elemental depletion with T in the BSE, arguing against scenarios in which the Earth formed by mixing of two components (85% volatile-poor; 15 % volatile-rich)[22,23]. A log-scale version of Fig. 1 is shown as Supplementary Figure 3.

Instead of partial evaporation, element depletion could occur without isotopic fractionation if the Earth accreted from mixtures of components in which an element is either present in solar abundances relative to RLEs (and therefore with chondritic isotope ratios) or almost entirely absent, with the likelihood of its presence/absence being a function of its condensation temperature. As most elements are calculated to condense from the solar nebula gas over a narrow temperature window (a few tens of Kelvin), the condensation of one element should be nearly complete before that of the next element has begun[13,14]. If condensation ceases at some threshold temperature, T, a step-function pattern in the planetary body is produced when the CI- (Ivuna-type carbonaceous chondrites), Al-normalised abundance of element i, f, is plotted against T [15]. Because such step-like depletion patterns are observed in small, rocky bodies (e.g., Vesta; Fig. 1), this feature was likely commonplace in the pebble- and planetesimal-sized building blocks of Earth. The Earth instead exhibits a smoother decline in elemental abundance over a wider range of T (Fig. 1) that differs from any chondrite. Indeed, no mixture of the extant suite of chondritic meteorites can simultaneously reproduce the isotopic and chemical composition of the Earth[16,17].

Abundances (x) of element i normalised to CI and Al in both the Earth and Vesta increase with increasing T in a manner approximated by a cumulative normal distribution;

The inflection point reflects the mean temperature, T, of the distribution. Weighted fits yield 1144±12 K and 1081±29 K, for Earth and Vesta, respectively. The steepness of the sigmoid is described by σ, the standard deviation, approximating a step-function at low σ (57±17 K, Vesta) or a gradual curve at high σ (226±16 K, Earth).

The mean temperatures experienced by the building blocks of Earth (1 AU) exceeded those of Vesta (2.36 AU), consistent with, but not indicative of, a heliocentric gradient in composition. Any such gradient among precursor bodies may have been subsequently perturbed by the [26]Al heating they experienced and/or their specific orbital and collisional histories. The higher standard deviation in T of the Earth indicates that it collected material with a wider range of temperatures than did Vesta. This result is in accord with dynamical models that predict larger bodies accumulate material from a broader range of heliocentric distances than do smaller ones[2,3]. On this basis, we propose that Earth did not accrete from chondrites, but rather stochastically from a summation of smaller precursor bodies, each with Vesta-like step-function element abundance patterns (σ tends to 0) and variable threshold temperatures (T) that reflect the temperatures at which they formed.

In order to quantitatively test this hypothesis, we examine four different N-body simulations[24] of the ‘Grand Tack’ scenario that reproduce the approximate number, masses and semi-major axes of the terrestrial planets. In this class of models, the migration of Jupiter, first inwards, then outwards by its resonant interaction with Saturn clears the region between 3 to 6 Astronomical Units (AU) of planetesimals and embryos, and scatters them into the inner disk between 0.7 and 3 AU, leaving a remnant population of planetesimals in the asteroid belt and beyond 6 AU[24,25] (Table 1). Any other dynamical model resulting in a comparable radial distribution of accreted material (Extended Data Fig. 1) would give equivalent results.

Table 1

Input parameters and final masses and semi-major axes of the analogue planet of the Earth resulting from the N-body simulations.

Simulation	Embryos(0.7 - 3 AU)		Planetesimals(0.7 – 3 AU)		Planetesimals(6 – 8 AU)		Earth
	Mass (ME)	n	Mass (ME)	n	Mass (ME)	n	Mass (ME)	Semi-Major Axis (AU)
GT/8:1/0.025/R2	0.025	212	3.81×10-4	1750	3.81×10-4	500	1.107	0.936
GT/4:1/0.025/R7	0.025	170	3.81×10-4	2783	4.67×10-5	1563	0.957	0.907
GT/4:1/0.05/R8	0.05	87	3.81×10-4	2836	5.15×10-5	1563	0.937	0.969
GT/8:1/0.08/R18	0.08	67	3.81×10-4	1750	3.45×10-4	500	1.074	1.179

The naming convention of N-body simulations follows i) the Grand Tack (GT) model, ii) the mass ratio of embryos:planetesimals, iii) the mass of embryos in ME and iv) the run number with those parameters[25].

Extended Data Figure 1

The sensitivity of the chemical and isotopic make-up of the Earth to the compositional model, described below, was evaluated by performing 1000 Monte-Carlo (MC) simulations with three free parameters for each N-body run (see Methods). It is assumed that all precursor bodies formed at the same time, meaning their initial temperature (T) reflects their heliocentric distance, d, and is calculated according to the radial temperature profile of a high-opacity nebular midplane, T, set by the disk in steady-state mass accretion[26,27] (Extended Data Fig. 2): where τ is the optical depth at the midplane (see Methods), M⨀ is the mass of the Sun and σ the Stefan-Boltzmann constant. The mass accretion rate to the Sun, (Ṁ), the first free parameter, is varied randomly within a normal distribution of 1(±0.25, 1σ)×10-7 M⨀/yr. Observations of T-Tauri stars show that Ṁ (M⨀/yr) decays as a function of time (Myr)[28]: indicating the time of planetesimal- or embryo formation. Given T, planetesimals contain CI-normalised abundances (f) of 30 fictive elements, defined with T (at 10-4 bar) ranging from 350 to 1800 K at 50 K intervals (see Methods):

Extended Data Figure 2

Fictive elements render the models independent of uncertainties or revisions in T of real elements[19,29], and permit quantification of systematic trends in their behaviour. If T > T of element i, f = 0, otherwise, f = 1, defining a ‘step-function’ abundance pattern in planetesimals. This distribution reflects i) their thermodynamic properties resulting in their perfect retention or evaporation[15] and ii) the inability of small (< 1000 km) bodies to retain an atmosphere[30]. The embryos comprise numerous planetesimals, each with a T sampled normally about the T of the embryo with a standard deviation, σembryo (cf. eq. 1). This, the second free parameter, has a mean of 87 K with standard deviation of 24 K (1σ) changes the smoothness of embryo compositions. This temperature range equates to a feeding zone width[31] of roughly ±0.1 AU for an embryo forming at 1 AU.

To assess the extent to which partial evaporation, and hence isotopic fractionation occurs during planetary accretion, the model computes a surface temperature increase in each body experiencing an impact that is proportional to the kinetic energy imparted by that impact (see Methods). Each fictive element vaporises with a partial pressure determined by its thermodynamic properties, with more volatile elements evolving higher partial pressures (see Methods). Given the surface atmospheric temperature and pressure, the temperature-pressure profile of a convective atmosphere with height following the dry adiabat is determined, before transitioning to an isothermal region above the homopause[32]. Elements in the atmosphere may be lost to space via Jeans and hydrodynamic escape. Which regime prevails depends on the escape parameter (Extended Data Fig. 3): where 𝑣 is the escape velocity of the body (= √2GMbody/rbody), k is the Boltzmann constant, T is the temperature at the escaping surface and m is the mean atmospheric mass. The value of, the third free parameter, is sampled randomly from a normal distribution with mean 33.5 g/mol and standard deviation = 3.5 g/mol, with limits corresponding to 23 g/mol (Na) and 44 g/mol (SiO). Jeans escape occurs for λ > 3, whereas hydrodynamic escape prevails for λ < 3 (see Methods)[33], while escape by other, non-thermal processes such as photoevaporation or ionisation are neglected. Escape fluxes are calculated at 10 K steps, assuming grey-body cooling of the atmosphere (see Methods). The composition of the residue is re-computed at each cooling step, and assumes element abundances equilibrate in the entire body.

Extended Data Figure 3

Results

Temperatures of precursor bodies set by a heliocentric gradient

In all simulations, refractory elements accrete in concert with the total mass of Earth (Fig. 2; Extended Data Fig. 4), because they are present in all planetesimals and are unaffected by atmospheric losses. By contrast, more volatile elements lag behind their refractory counterparts. Simulations with high atmospheric loss rates (low m) are indistinguishable from cases in which atmospheric escape is suppressed (high m), meaning the differential timing of volatile accretion depends only on precursor body compositions set by the value of Ṁ and its original heliocentric distance, d.

Figure 2

Timing of element accretion to the Earth as a function of volatility.

Points represent the time in the N-body simulation (ordinate) at which a given element (abscissa) reaches a threshold value (10 %, red; 50 %, green; or 90 %, purple) of its final abundance. These points are % accreted = 100×(x)(x), where the subscripts t and f denote the threshold time and the final time, respectively, x is the mole fraction of element i and M is the mass of the growing Earth. Error bars (1σ) denote the range given by Monte Carlo models for different Ṁ, σembryo and m across 1000 simulations. No error bar indicates that the element reaches the given % accreted value at the same time in each MC simulation. The mass of the Earth, M, through time is delineated by horizontal lines at 10 %, 50 % and 90 %. Offset between these lines and the coloured points is the differential accretion time as denoted by the grey shaded fields. The large error bars on refractory elements arise from their dilution/enrichment in planetesimals and embryos with/without the more abundant elements with T near 1350 K (Supplementary Table 1).

Extended Data Figure 4

Our results indicate that, while the proto-Earth accreted ~80-90% of its final total mass by 25 – 100 Myr, only ~50 % of the most volatile elements (Tc50 < 500 K) accreted within this timeframe. This reflects the preferential accretion of volatile-undepleted (i.e., CI-like) planetesimals from beyond 6 AU later in Earth’s growth. Late accretion of this material (e.g., during a giant impact) provides a physical explanation for the requirement that volatile-rich material constituted the final ~10 % of Earth’s mass and accreted after 30 Myr in order to explain the Earth’s carbonaceous chondrite-like 107Ag/109Ag ratio, and to reconcile Hf-W and Pd-Ag ages[34].

Once fully accreted, modelled elemental abundances in the BSE adhere to a gradual decline with decreasing T that parallel observations (Fig. 3a, b), independent of the specific N-body simulation (Extended Data Fig. 5). Although a constant-abundance plateau is observed in some simulations, it extends only to elements with T < 450 K, and never up to 750 K or 1000 K as suggested in recent models[17,23,35]. Monte-Carlo simulations indicate that the composition of the BSE depends strongly on Ṁ but is nearly independent of σembryo and m (Extended Data Figs. 6). Therefore, its final composition is inherited from its precursors, rather than modified subsequently by partial evaporation during collisions.

Figure 3

Elemental abundances in the fully accreted Earth.

CI-, Mg-normalised elemental abundances in the BSE[6] as a function of 50% nebula condensation temperatures according to Wood et al., red[29]; and Lodders with updates from Fegley and Schaefer, purple[19,20], compared with abundances predicted by 1000 Monte Carlo simulations per N-body simulation plotted on a) a linear scale and b) a logarithmic scale. The average composition over 1000 runs is given by the solid blue line with the percentile ranges of the simulations (25th to 75th, light blue and 10th to 90th, dark blue). The dashed lines show the simulations that minimise the Root-Mean-Square (RMS) deviation from the data. The RMS values are listed in the top left-hand corner, along with those of the three random variables to which they correspond. Only lithophile elements were plotted, as siderophile elements are further depleted by core formation, which is not considered in our model.

Extended Data Figure 5

Extended Data Figure 6

The elongated (high-σ) shape in the BSE arises from the summation of precursor bodies with a wide range of T, which, for a given Ṁ, relates to heliocentric distance (Extended Data Fig. 2). This follows from the condition applied to planetesimals that f = 0 when T > Tc50, which becomes increasingly probable as T decreases. Even though planetesimals initially have step-function patterns, this process produces a smooth pattern in the BSE by the Central Limit Theorem. The corollary is that step-function patterns (low-σ) should be more evident in smaller bodies that form in local feeding zones, in which precursors have a restricted range of T, as observed for Vesta[15] (Fig. 1).

Best fits to the lithophile elemental abundances in the BSE are determined by minima of the Root-Mean-Squared (RMS) deviation of models from the data (Fig. 3; Extended Data Fig. 5). The RMS is minimised for Ṁ = 1.08(±0.17, 1σ) × 10-7 M⨀/yr (Extended Data Fig. 6), and is insensitive to the specific accretion history in the subset of the four N-body simulations investigated here. Assuming gradual, rather than step-like, patterns for planetesimals (as per those of embryos) does not change the best fit value of Ṁ. These mass accretion rates yield T = 1410±60 K at 0.7 AU, coinciding with the T of the major rock-forming elements, Mg, Si and Fe[19].

Observations of T-Tauri stars show that younger disks have faster mass accretion rates than do older ones[36] (eq. 3). By analogy, accretion rates of ~10-7 M⨀/yr indicate rapid establishment, ~0.2±0.1 Myr, of the compositional architecture of the inner disk; within the lifetime of the nebular gas. Detailed numerical models also support early formation of planetesimals in the inner disk at ~0.5 Myr[37] or ~1 Myr[38]. Moreover, this timescale matches volatile depletion ages based on Rb-Sr and Mn-Cr systematics of small differentiated bodies (<1.5 Myr[39,40]) and of the Earth (1-2 Myr[41]). Subsequent, post-nebular orbital evolution produces its smooth (high-σ) pattern over longer timescales, nearing ~10[8] yr.

The stable isotope composition of the Earth represents the sum of contributions from partial evaporative losses and accretion of unevaporated materials. The isotope fractionation evolved is proportional to the amount of the element lost and the fractionation factor. In the N-body simulations, isotope ratios for each element in planetesimals and embryos are initially CI chondritic, and atmospheric escape is calculated for all bodies that experience a collision. A peak is observed in elements with 450 < T (K) < 800 across all simulations (Fig. 4; Extended Data Fig. 7). These elements are not so volatile so as to have been lost completely (and thus prone to overprinting by late-accreted material) but still sufficiently volatile to record partial evaporation over Earth's accretion. The magnitude of isotope fractionation in the fully accreted Earth is proportional to 1/m (Fig. 4, Extended Data Fig. 8).

Figure 4

Stable isotope fractionation modelled in the fully accreted Earth as a function of volatility relative to observations.

Isotopic fractionation for elements for which the difference between Earth and CI chondrites is well-constrained, Cl10, Zn8, Rb11, K9, Si42, Fe43 and Mg44, is shown in the top panel together with their 1σ uncertainties. The model results are displayed in the bottom panel, where the red line indicates the mean isotopic fractionation (in per mille relative to the initial value, i.e., CI chondrites) across 1000 Monte Carlo runs. The isotope ratios in g/mol (m = numerator isotope, Δm = 2 atomic mass units, amu) are listed in the legend. The vertical dashed lines delimit coloured zones pertaining to Late Accretion (blue, volatile elements whose isotopes are present in chondritic ratios), Partial Evaporation (green, elements that bear an isotopic signature of partial evaporation during atmospheric loss) and Non-Volatile (red, elements that were not lost and remain chondritic). The scale of isotopic variation observed exceeds that produced in the Monte-Carlo runs, except for low m runs in N-body simulation GT/8:1/0.025/R2 (Extended Data Fig. 7). Moreover, the observations do not follow the behaviour expected for partial evaporation, which most readily generates heavy isotope ratios for elements with 450 < T (K) < 800.

Extended Data Figure 7

Extended Data Figure 8

Such isotope fractionation patterns are explicable through the dynamics of planetary impacts. The tendency for atmospheric loss (and hence isotope fractionation) to occur is quantified by the escape parameter, λ (eq. 5). As the Earth grows, λ gradually increases, upon which are superimposed embryo collisions (ME > 0.05) that produce transient decreases in λ via impact heating (Extended Data Fig. 3), during which atmospheric loss is most pronounced. For any given N-body simulation, there are typically one- to two collisions that induce hydrodynamic escape. These events are restricted to the nascent stages of Earth’s formation, during which its mass is ~0.05 to 0.2 ME (Extended Data Fig. 3). Even though temperatures in later, giant impacts may be higher, the proclivity for atmospheric loss is offset by the larger escape velocity, as λ depends on 𝑣2, but only on 1/T (eq. 5).

As such, chemical and isotopic fractionation generated by early, high-energy events is readily overprinted by accretion thereafter. The most volatile elements are delivered late by material that has seen no evaporation, and hence has chondritic isotope ratios (Fig. 4). The isotope composition of elements with 450 < T (K) < 800 depends largely on the Earth’s accretion path and hence the N-body simulation. A small fraction (<5 %) of runs with the lowest m in the most energetic simulation, GT/8:1/0.025/R2, reach 3 ‰ (Extended Data Fig. 7; 8). For all other runs, the limited degree of isotopic fractionation is consistent with the observed <0.4%[8-11] deviation of the BSE from CI chondrites in K (~1000 K), Rb (~800 K), Zn (~700 K) and Cl (~400 K) isotope compositions (Fig. 4). Therefore, partial atmospheric escape from large bodies (≥0.2 ME) had minimal influence on the final compositional and isotopic make-up of Earth. By extension, isotopic variations of less volatile elements, including Mg, Fe and Si[42-44] (Fig. 4) are unlikely to originate from partial evaporation of embryos. Instead, they could have been inherited earlier from nebular condensation[42] or planetesimal evaporation[45]. Nevertheless, consideration of the specific thermodynamic properties of these elements over a wider range of N-body simulations is required. Conversely, the isotopic imprint of atmospheric loss is more readily preserved on small bodies for which escape is both more pronounced and less overprinted by continued accretion, consistent with enrichment in heavy isotopes of K, Cl and Zn of the Moon and Vesta[46-48].

Temperatures of precursor bodies set by radiogenic heating

Aluminium-26, which decays to [26]Mg (t1/2 = 0.7 Myr), was abundant enough to have caused melting and vaporisation on small rocky bodies, provided they formed sufficiently early (≤2 Myr)[49]. Therefore, a variation on the stochastic accretion hypothesis arises if the temperatures of precursor bodies reflect [26]Al decay in addition to, or, in the extreme case, instead of the smooth heliocentric decrease in disk temperature modelled above. The endmember case is investigated by randomly assigning each body an initial temperature selected from a normal distribution with mean, T = 1000(±50, 1σ) K and standard deviation, σ = 175 K (Extended Data Fig. 5). The T value can reach ~1600 K, consistent with numerical models[50], and replaces Ṁ as the first random variable in the MC simulations. The remainder of the simulation is identical to that described above.

This temperature distribution also explains the abundance pattern of the BSE (Extended Data Fig. 9). However, because the compositions of planetesimals are independent of heliocentric distance, there is no systematic distinction in the accretion times of volatile and non-volatile elements, a result that contravenes evidence from the Pd-Ag and Hf-W systems[34]. This inconsistency can be resolved by numerical models that predict slower accretion of planetesimals with increasing heliocentric distance, causing the effects of [26]Al heating to diminish[38,49]. In this scenario, a planetesimal population in the inner disk is predicted to have formed at 0.6 – 1.3 Myr[38]. Therefore, although differential [26]Al heating perturbs the thermal structure imposed by the disk, it also produces the chemical and isotope composition of Earth by stochastic accretion, provided a heliocentric gradient is present.

Extended Data Figure 9

Growth of Earth by pebble accretion

The formation of Earth via pebble accretion[51], a popular alternative mechanism for planetary growth[52], is simulated by prescribing a mass accretion rate that scales with ME2/3, typical of 2D-pebble accretion in the Hill regime[53]. Pebble fluxes are assumed to mirror mass accretion rates to the Sun (eq. 3) and result in the growth of Earth to 0.9 ME in 3.85 Myr, consistent with disk accretion models (~5 Myr[5]). Pebbles initially form at the snow line[4] and subsequently drift inwards by interactions with the surrounding nebular gas. Because pebbles are small, they thermally equilibrate with the gas at 1 AU, whose temperature is calculated over time with eqs. (2, 3). Thus, elemental abundances in pebbles at the time of their accretion by the Earth are determined as per planetesimals, but using the temperature at 1 AU as T.

Because pebbles accrete continuously in time, Earth’s surface reaches only modest temperatures compared to planetesimal accretion, ~1750 K, sufficient to melt and vaporise the pebbles. However, the escape of vaporised material occurs in the upper atmosphere, near the Hill radius, or deeper, at the Bondi radius, if Earth is still embedded in the disk[54]. At these locations, temperatures drop to ~400- to 100 K[5], meaning elements with T >400 K are solids. Consequently, atmospheric loss and isotopic fractionation is muted (Extended Data Fig. 10). Pebble accretion models produce concave-down patterns that decline at T lower than observed in the BSE (Extended Data Fig. 10). This results from the dependence of accretion rate on ME2/3, meaning the majority (85 %) of Earth’s mass accretes between 1 and 3.85 Myr, when temperatures at 1 AU are ≤500 K (Extended Data Fig. 10). As such, depletion of elements with T above this temperature is limited. Therefore, unless pebbles themselves are already impoverished in volatile elements, pebble accretion alone does not produce the volatile depletion pattern of the Earth.

Extended Data Figure 10

Discussion

The distribution of initial temperatures, and hence compositions of accreted planetesimals and embryos, controlled the composition of the Earth. In order to achieve sufficient temperatures in these precursor bodies, either by heating in the disk and/or [26]Al decay, they must have formed by ~1 Myr, a timescale that is consistent with radiometric ages for volatile depletion[39,41] and protoplanetary disk models[37,38]. Whether precursor body compositions were set by nebular condensation or planetesimal evaporation could be discriminated by holistic models of the growth of planets from the dust- to embryo stages. In the context of the Grand Tack, and other models in which the Earth accretes from a radial distribution of material (Extended Data Fig. 1), some form of heliocentric temperature gradient is necessary to explain the late accretion of volatile-rich bodies (Fig. 2) that satisfy the near-chondritic radiogenic[34] and stable[8-11] (Fig. 4) isotope composition of volatile elements in the Earth.

The mean volatile content of a planet, quantified in its T, is expected to increase with heliocentric distance. Should a heliocentric gradient have been established by disk accretion, then a continuous decrease in T of the terrestrial planets would be observed, whereas [26]Al heating may produce a dichotomy in the compositions of inner- and outer Solar System planetesimals due to their differential formation times. As they grow, planets accrete material from an expanding range of heliocentric distances[2,3]. This translates into a greater distribution of T in precursor bodies, implying that larger planets should have more gradual depletion patterns (high-σ) than smaller ones (low-σ) (Fig. 1). Moreover, smaller planets are more likely to preserve isotopic imprints of atmospheric loss, resulting in heavy stable isotope compositions of volatile elements relative to their building blocks. Potassium isotopes in the Earth, Mars, the Moon and Vesta bear out this expectation, as they become increasingly non-chondritic with decreasing planetary mass[55]. The compositions of the terrestrial planets, together with their masses and heliocentric distances, offer the means by which to test these predictions.

Methods

N-body simulations

N-body simulations proceed by perfect-merger impacts that result in a net increase in mass of larger bodies at the expense of smaller ones (oligarchic growth) to produce three- to four terrestrial planets with masses and semi-major axes akin to those observed at the present-day[24].

We treat four simulations, labelled GT/4:1/0.05/R8, GT/4:1/0.025/R7, GT/8:1/0.08/R18 and GT/8:1/0.025/R2. See Table 1 for the masses and distribution of embryos and planetesimals, and refs. 24,25 for a detailed description of the simulations.

Disk structure and initial composition

All bodies (embryos and planetesimals) are assigned an average density of 3000 kg/m3. They are prescribed to initially contain some fraction of the CI-normalised abundances (f) of 30 fictive elements (i) defined with T (at 10-4 bar total pressure) ranging from 350 to 1800 K at 50 K intervals, where: and x corresponds to the mole fraction. Where x, the body has its full complement of element i, and is equivalent to a ‘CI-chondritic’ composition. The initial temperature of the planetesimal or embryo is determined by its heliocentric distance. The temperature profile of a high-opacity nebular midplane in steady-state mass accretion is[26,27]: where i is the optical depth, and is given by τ = κΣ/2 (ref. 56), where κ is the Rosseland mean opacity (~0.3 m[2]/kg in the range 100 – 1500 K[57]) and the surface density of the gas, Σ = 3500(Ṁ/10-8)(d/AU)-3/4 in kg/m[2], suitable for the optically-thick accretion phase of the star, Ṁ > 10-9 solar mass/yr, in the planet-forming region[27,58]), M⨀ the mass of the Sun, Ṁ the mass accretion rate, σ the Stefan-Boltzmann constant, and d the heliocentric distance of the body given in the N-body simulation. Ṁ is the first random variable for the Monte-Carlo simulations (see Monte-Carlo simulations).

We impose the constraint that, for planetesimals, if T > Tc50 of element i, f = 0, otherwise, f = 1, defining a ‘step-function’ pattern. As embryos are larger bodies, they themselves are expected to have been comprised of planetesimals with step-functions. Hence, abundances in the embryos reflect the summation of a number, N = (Membryo/Mplanetesimals) of planetesimals, each with a T sampled normally about the T of the embryo, with standard deviation σembryo. If the σembryo is low, the depletion pattern in the embryo is closer to a step-function, while high σembryo leads to a smoother curve. σembryo is the second random variable for the Monte-Carlo simulations (see Monte-Carlo simulations).

Energy budget during impact heating and cooling

Collisions result in a temperature increase according to the conversion of kinetic energy to heat. The energy released upon impact of two bodies of masses M and M travelling at velocities ν and ν, respectively, is the sum of their kinetic energies:

However, only a fraction of the total kinetic energy is converted into heat, depending on the characteristics of the impacting material and obliquity and mass of the impact[59,60], while the heated region depends on the amount of energy transmitted beyond the isobaric core of the body. These poorly-known terms are consolidated into a single empirical constant, h, which ≈0.5[61]. The temperature increase is the total heat distributed throughout the affected mass, divided by the heat capacity of the material C :

Where C is that for molten peridotite near its liquidus (1800 J/kg.K[62]), M is the final mass of the planet, and ν is the impact velocity given in the N-body simulation. Temperature increases if the energy delivered is in excess of the latent heat of fusion (ΔH = 4 × 105 J/kg at 1400 K[63]), then vaporisation (ΔH = 5 × 106 J/kg at 2000 K[61]) and finally an upper limit at 6000 K, whereupon behaviour becomes supercritical[64,65]. The mass of material that experiences peak temperatures (the isobaric core) scales with impactor radius[44,61]. However, because N-body simulations represent a swarm of planetesimal impacts, the heated area is likely to have been uniformly distributed in the body.

Following heating to maximum temperatures, the planetary body cools by a combination of black body radiation, latent heat of vaporisation and gravitational potential energy, with the energy balance given by[66]:

The contribution of the evaporative cooling and escape term, is small and is neglected. Here, r is the planetary radius (assuming a radiating atmosphere from the surface), σ the Stefan-Boltzmann constant, and the surface temperature, T. The emissivity, ϵ, takes into account that planets are not perfect black bodies, and is set to 0.05 to replicate cooling timescales for terrestrial planets with thick atmospheres[67-69].

Temperature is calculated at 10 K cooling intervals throughout the simulation. This temperature is then fed in as an input to calculate the composition of the atmosphere (Vapour pressure), its subsequent loss at each temperature step (Atmospheric loss), and thus to re-calculate a new bulk composition (Compositional evolution) in a loop.

Vapour pressures

For silicate materials, vaporisation reactions can be generalised as: where is the metal oxide species in the silicate melt, is the stable gas species with metal oxidation state of x, balanced by moles of O2 gas for a congruent reaction (i.e., in which the gas(es) produced have the same composition as the reactant(s)).

At equilibrium, the partial pressure (assuming ideal gas behaviour such that f = p), is given by: where, for which ΔG( is the Gibbs Free Energy of the vaporisation reaction, R the gas constant and T the equilibrium temperature set by the magma ocean at its surface in Kelvin. At low pressures in which volume changes are insignificant relative to the standard state (here taken at 1 bar),

In which ΔH is the enthalpy change and ΔS the entropy change of reaction. Substituting this equation into eq. (12), one obtains:

From experimental assessment of the equilibrium vapour pressures of metal-bearing gas species above metal oxides or metal oxide components in silicate melts, ΔS is near-constant[15,70], which results in slopes in logp space that are sub-parallel for almost all elements (Supplementary Figure 1). Element volatilities are defined by their 50% condensation (or evaporation) temperature in the solar nebular gas[15]:

Here, f = the fraction in the vapour (=0.5 for 50% condensation). The fO2 is that of a cooling solar nebula[71], while the total pressure is 10-4 bar in keeping with convention[19,29]. The T of a given fictive element is defined by holding ΔS constant and changing the value of ΔH in eq. (16). Their properties are shown in Supplementary Table 1.

The partial pressure is determined as a function of temperature, element activity, fO2 of the atmosphere and the stoichiometry of the vaporisation reaction, n (eq. 11). This links the T of an element with its general volatility behaviour under planet-forming conditions. Because most elements evaporate according to n = 2 stoichiometries[15], including all five major cation-forming elements in the BSE composition (Fe, Mg, Si, Ca, and Al), we adopt n = 2 for all fictive elements (Supplementary Table 1). Hence, relative fractionation by volatility is insensitive to fO2, which is not the case for the alkalis (n = 1), or for the Group VI metals (Cr, Mo, and W) for which n < 0[15]. However, this simplification is apt to establish the overall trend of element depletion, and not to account for anomalies in the abundances of certain elements. Moreover, it is independent of any eventual revisions to T values.

The sum of the partial pressures dictate the total pressure at any given T and fO2: since depends on (fO2)-n/4, congruent evaporation (eq. 11) demands that oxygen fugacity be equal to n/4 times the sum of the partial pressures of the metal oxide species. Because n = 2 in our treatment (Supplementary Table 1), then:

Thus, the value of fO2 is solved iteratively until the term 0.5P is equal to ½ the sum of the individual partial pressures of each component at a given temperature. Here, we constrain element activities by fitting our value of P to that calculated for evaporation of silicate mantles[72]. To do so, we assign an activity to each fictive element as a function of its T value. As per planetary materials, the three most abundant elements (Fe, Mg and Si) have T ~ 1300 – 1400 K. We use a Gaussian distribution to model a peak element activity about 1350 K. The rate at which element activity declines is found by minimisation to the P determined in ref. 72 over the range 1800 < T (K) < 3000, using the objective function:

By iterating the value of c:

We obtain a best-fit value to the total pressure above the bulk silicate Earth[72] (Supplementary Figure 2) at c = 100, yielding the initial activities, (xγ)initial, shown in Supplementary Table 1.

Atmospheric loss

Escape of the atmosphere is calculated self-consistently based on its pressure-temperature structure. The total pressure at the surface, P, (eq. 17) is a direct result of the planetary surface temperature, T, given that any crustal boundary layer is likely to be ~ a few cm thick[44,73] and prone to foundering. Even for small bodies ½ the mass of Pluto, the evaporation rate quickly attains a steady state with loss rate, such that the surface pressure is that at equilibrium[45]. A rapidly convecting terrestrial magma ocean (viscosities ≈ 0.1 Pa.s[74]) experiences turnover in days, meaning evaporation is not diffusion-limited by transport of material to the surface.

Throughout Earth’s accretion at any given temperature step, we model two thermal escape regimes, Jeans and hydrodynamic escape. Energy-limited or ionisation-driven escape are not considered as they require knowledge of solar energy flux. The prevailing mode of escape is determined by the competing influence of gravity and thermal energy, encapsulated in the escape parameter, λesc[75]: where υ is the escape velocity of the planet at its exobase, m is the average molar mass of the atmosphere, k is the Boltzmann constant, and T is the temperature at the exobase. The quantity m is the third random variable for the Monte-Carlo simulations (see Monte-Carlo simulations). Atmospheric structure is described by its scale height, H, the height over which its number density declines by a factor of 1/e. For a well-mixed atmosphere, H is evaluated at the tropopause by: where, where r is the radius at the surface, and r is the radius at the tropopause (rtrop = rs+z), M is the mass of the planet. Thus, the temperature, T and pressure, P set by evaporation at the surface need to be extrapolated upwards to the tropopause. This is achieved using the expression for a convective, (dry adiabat) optically-thin troposphere[32]:

Here, R = the gas constant and C is the heat capacity of an ideal, complex gas at constant pressure (=4R). For an optically-thin atmosphere in the absence of external heat sources, the temperature tends to a constant value given by the ‘skin temperature’[32] (=T),

Combining eq. 25 with eq. 24 leads to P = 0.5 P. This limit, the tropopause, defines the transition from the convective troposphere to the diffusive stratosphere that overlies it. The adiabatic expansion of an ideal gas in hydrostatic equilibrium as a function of height above the surface (z) is: where z at the tropopause is z, and is therefore given by:

As H appears in both eq. 27 and eq. 22, these equations are resolved iteratively until they converge upon a unique value of z and H. Above this point, in the stratosphere, temperature as a function of height (z) is constant, reflecting the ‘skin temperature’. The number density at the exobase, which is defined as the point at which the mean free path (ℓ = 1/nσ) equals the scale height (H), is: where σ is the collision integral:

The escape velocity at the exobase is then calculated:

Where the radius of the exobase is given by adiabatic expansion:

The number density at the surface is calculated by the ideal gas law:

The mean molar mass of the atmosphere, m, is held constant throughout the simulation. The escape parameter, λ, at the exobase of the body can then be evaluated. If λ > 3, the position of the exobase is further out than the radius at which the gas exceeds the sound speed, and Jeans escape prevails. In this regime, the mass loss rate is: or, in kg/s ;

If λ < 3, the gas exceeds the speed of sound and a transonic outflow of the atmosphere results. The gas velocity is determined by the temperature (T), total pressure (P) and density (ρ) of its constituents at any given point, constrained by mass and momentum conservation[76]. At the ‘sonic point’, the velocity of the gas reaches that of the speed of sound in the medium (Cs). For an ideal gas, Cs is given by:

Where the factor 7/5 is the value of the adiabatic index for a diatomic gas (5 degrees of freedom), and accounts for the heat of compression. The sonic point, or Bondi radius is:

The hydrodynamic mass loss rate depends on the density of the gas phase at the Bondi radius. In order to estimate this value, the pressure of the gas at the surface is integrated outwards to r. For an adiabatically expanding gas, the relation is:

The number density (units of atoms/m3) at the sonic point is:

The molecular loss rate (molecules/s) is then: or, in terms of kg/s ;

In order to estimate isotopic fractionation during atmospheric loss, two simulations are run concurrently by varying the term m (which is a bulk property of the atmosphere in the model framework and cannot be independently varied for each element) according to:

Here, the quantity (i – j) is set to 2, reflecting the isotopic mass difference between most common systems (e.g.,[66]Zn/[64]Zn; [30]Si/[28]Si, [41]K/[39]K). The scale height of both isotopes, i and j are calculated independently, such that eq. 22 becomes: simulating the separation of two isotopes by their scale heights in the stratosphere. Therefore, both Hydrodynamic and Jeans escape may cause isotope fractionation. Jeans escape rate is further modified due to the dependence of the relative mass loss rate dm on (m)0.5. As such, eq. 33 becomes:

The mean molar mass factor used to convert the loss rate from molecules/s to kg/s in both Jeans and Hydrodynamic escape is left constant (m), as m is assumed to be an infinitely dilute component of the atmosphere only (i.e., x >>x). Thus, the theoretical isotope fractionation is calculated for each element, presuming they all have the same value of m and m.

Compositional evolution

The activity, a, of each element, i, in any planetary body is given by: where (x) is given for each i by a Gaussian distribution centred at 1350 K, whose coefficients are adjusted such that the total pressure evolved as a function of temperature fits that for thermodynamic models for vaporisation of a bulk silicate Earth composition[72] (Supplementary Figure 2, eq. 20). For simplicity, here we assume ideal solution of each element, γ = 1. In the entire planet, therefore;

The composition and total pressure of the gas phase in equilibrium with the planet at each temperature step is calculated according to the section Vapour pressure calculations, a proportion of which is lost by Jeans or hydrodynamic escape, for which a mass flux is determined (i.e., a loss of mass over a time interval, see Atmospheric loss). As the mass of the planet is known, this can be expressed in terms of the total planetary mass, which gives the normalised loss rate in seconds:

In the model, the temperature is evaluated for a fixed time interval, dt, dictated by the time taken for the planet to cool by 10 K. Hence, multiplying the normalised loss rate by the seconds elapsed during the cooling interval gives the integrated mass loss (expressed as fraction of planetary mass) over that time interval:

The new activity of a given element is:

Where is the mole fraction of the element i in the gas phase being lost. Here, we assume that there is no fractionation of elemental mole fractions in the gas at the loss surface with respect to that set by equilibrium with the magma ocean at the planetary surface. That is, we neglect processes in the stratosphere that may fractionate elements according to their scale heights because the molar mass of the atmosphere is assigned as a bulk property, not a sum across each fictive element.

The output composition is then normalised to the new total mass of the body i.e.,

This process is repeated, manifest as a compositional evolution of the planet dictated by the element’s volatility and the rate of atmospheric loss. Compositional evolution with time is shown in Supplementary Figure 4. For isotope fractionation, this process is performed for elements i and j simultaneously throughout, such that the isotope composition (in per mille) of the body can be calculated at any time according to:

Monte-Carlo simulations

For each simulation, there are three random variables:

Ṁ. The stellar mass accretion rate. It sets the initial abundances of all planetesimals and embryos. Ṁ is randomly varied according to a normal distribution with a mean of 1 × 10-7 M⨀ and a standard deviation (σ) of 0.25 × 10-7 M⨀ in accord with observations of T-Tauri stars[77]. These accretion rates yield maximum temperatures (at 0.7 AU) of K (1 σ), bracketing the T of the major rock-forming elements. This range of temperatures is a robust constraint because refractory lithophile elements (T ≥ 1400 K) are unfractionated from one another in planetary materials[78], providing empirical evidence that there was no significant accretion of components with fractionated abundances of elements with T ≥ 1400 K.

σembryo. The standard deviation for the temperature distribution of planetesimals that comprise a given embryo. The standard deviation, σembryo, is randomly varied with a normal distribution with mean 87 K and σ = 24 K about the mean Tmidplane value of the embryo in each simulation. This physically relates to embryo accretion from planetesimals sourced from variably wide heliocentric feeding zones[31,79]. Low values of σembryo result in a ‘step-function’ distribution of elements, similar to planetesimals, whereas high values of σembryo pertain to accretion of planetesimals from a wider temperature range, leading to a smoother element distribution pattern in embryos as a function of T.

m The mean atmospheric molar mass. A bulk property, it is sampled randomly from a normal distribution with mean 33.5 g/mol and σ = 3.5 g/mol, with minima and maxima corresponding to the masses of Na(g) (23 g/mol) and SiO(g) / CO2(g) (both 44 g/mol), respectively. These are likely to represent the limits of the most abundant elements in a high-temperature vapour in equilibrium with a silicate mantle[72].

Monte Carlo simulations are performed 1000 times for each N-body simulation. The three variables are fixed for each simulation, selected according to the range and standard deviation for each normal distribution. Sensitivity tests as to how each parameter influences the final elemental abundance of the Earth is shown in Supplementary Figure 5. The entire process described above has been coded together in a Python script, and was run on the Ludwig-Maximilians-Universität (LMU) supercomputing cluster. Each simulation takes roughly 1 minute to complete (1000 simulations = 1000 minutes).

Pebble Accretion Simulation

To simulate the accretion of the Earth predominantly by pebbles[80], a simple model is presented, in which the Earth grows to 90 % of its current mass solely by pebble accretion in 3.85 Myr. This process imparts heating due to the release of gravitational potential energy; from which a surface temperature is determined.

Because pebbles are small (~cm-sized) they thermally equilibrate with the surrounding nebular gas at the locus of accretion, that is, at 1 AU for Earth. In order to calculate the temperature at 1 AU, an estimation is required of how the mass accretion rate to the Sun decays as a function of time. To do so, eq. 4 of ref. 28, based on observations of the mass accretion rates of T-Tauri stars, is used to calculate Ṁ (eq. 3). For each given time step (t), the corresponding value of Ṁ is inserted into eq. 2 to calculate the temperature of the accreting pebbles. Pebble fluxes are assumed to mirror mass accretion rates (eq. 3) and scale with ME2/3, typical of 2D-pebble accretion in the Hill regime[53]. Given a temperature, their composition is determined in a manner analogous to that for planetesimals (see Disk structure and initial composition). The first and second terms of eq. 3 are varied within their uncertainties in a Monte Carlo simulation to evaluate the sensitivity of Earth’s composition to mass accretion rate.

Because eq. 2 predicts very high temperatures in the first few thousand years (up to 3800 K at 5000 yr), a saturation limit is placed at 1350 K, which affects only the first 1.5±0.5 % of material accreted to the Earth. The vapour pressures, atmospheric loss and compositional evolution of Earth is treated in an identical manner to that outlined above, determined by the surface temperature given in the pebble accretion simulation.