# Early accretion of planetesimals unraveled by the thermal evolution of the parent bodies of magmatic iron meteorites

Edouard Kaminski, Angela Limare, Balthasar Kenda, Marc Chaussidon

► To cite this version:

Edouard Kaminski, Angela Limare, Balthasar Kenda, Marc Chaussidon. Early accretion of planetesimals unraveled by the thermal evolution of the parent bodies of magmatic iron meteorites. Earth and Planetary Science Letters, 2020, 548, pp.116469. 10.1016/j.epsl.2020.116469. insu-02931908

HAL Id: insu-02931908

https://insu.hal.science/insu-02931908

Submitted on 12 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Early accretion of planetesimals unraveled by the thermal evolution of the parent bodies of magmatic iron meteorites

Edouard Kaminski<sup>a,*</sup>, Angela Limare<sup>a</sup>, Balthasar Kenda<sup>a</sup>, Marc Chaussidon<sup>a</sup>

<sup>a</sup>Université de Paris, Institut de Physique du Globe de Paris, CNRS, 1 rue Jussieu, F-75005 Paris, France

## Abstract

The timing of formation of 100-300 km size planetesimals in the protoplanetary disk remains largely unconstrained. Recent models show that gravitational collapse of boulders in overdense regions of a dusty accretion disk can overcome the meter-sized barrier and lead to rapid formation of planetesimals with size of several km that further grow by pebble accretion. Hf/W ages indicate that the first large planetesimals to form could be the parent bodies of magmatic iron meteorites. These ages have been so far used to constrain timing of accretion considering (i) instantaneous accretion, and (ii) purely conductive heat transfer in the planetesimal. To relax these hypotheses we model the thermal evolution of a planetesimal in course of accretion and we take into account the possibility of convection onset. Our model is further based on considering the possibility of a common thermal evolution for all the parent bodies of iron meteorites. Within that framework we show that

<footnote>*Corresponding author
Email address: kaminski@ipgp.fr (Edouard Kaminski)</footnote>
the planetesimals could have grown following a universal accretion law start-
ing at the very beginning of the history of the disk by a nearly instantaneous
formation of 60±30 km size nuclei, followed by a growth via pebble accretion
at a much slower pace to reach final sizes of 150–300 km in about 3 Myr.
In this universal scenario, complete melting and total differentiation are not
bound to happen in the parent body due to the continous accretion of cold
pebbles. The model, though calibrated here on iron meteorites, is general
and can in principle be applied to other types of planetesimals such as for
instance the parent bodies of CV chondrites.

***Keywords:*** planetesimals, accretion, convection, analog experiments, magmatic iron meteorites

**1. Introduction**

The formation of planetesimals is a crucial stage in the evolution of the
protoplanetary disk that leads from dust to planets. Some of its fundamental
aspects remain the object of debate, such as the range in initial and final
sizes of the planetesimals, and more generally their mechanism and rate of
accretion (e.g., Xie et al., 2010; Morbidelli and Raymond, 2016; Birnstiel
et al., 2016).

A key and long-standing question has been the process of formation of
the initial seeds of the planetesimals. Models of collision between aggregates
of dust held by electrostatic forces fail to produce planetesimals, because
the particles start to bounce upon collision when they reach few centimeters
and migrate rapidly toward the star as they grow (e.g., Morbidelli and Ray-
mond, 2016). Meter-size boulders migrate so fast that they are eventually
lost to the sun. One way to overcome this meter-size barrier (Weidenschilling, 1977) is to consider a weakly turbulent nebula where large and dense clumps of chondrules can form on the smallest scale of the flow, and be preserved by self-gravity to later produce platenesimals by sedimentation towards their center (Cuzzi et al., 2008). Numerical simulations (e.g., Johansen et al., 2015; Simon et al., 2017) show that other mechanisms such as quasi-stable pressure bumps and/or streaming instabilities can produce gravitationally-bound pebble-swarms collapse under their own gravity, yielding initial planetesimal sizes up to many tens of kilometers across.

Once the seeds of planetesimals are formed, they are supposed to continue to grow by pebble accretion (Johansen and Lambrechts, 2017). This mechanism is controlled by the interplay between gravity, aerodynamics deflection and gas drag (e.g., Ormel and Klahr, 2010). Models predictions give the minimum planetesimal radius for the onset of pebble accretion as well as the growth timescales for a given radius at a given distance from the sun (Visser and Ormel, 2016). But these models have not been used yet to predict complete accretion scenarios for the planetesimals and to confront the output with independent constraints.

Constraints on the accretion of planetesimals in the disk can be obtained by astronomical observations (e.g., van der Marel and et al., 2013; Delbo' et al., 2017) as well as by the study of meteoritic fragments that sampled the early stages of the disk evolution during which planetesimals formed. The very early ages demonstrated by Hf/W studies of iron meteorites presumably originating from planetesimals accreted inside and outside the orbit of the growing Jupiter, indicate that these planetesimals and a Jupiter's core of ≈20
Earth mass may have formed within the first million year (Myr) of the solar system (Birstel et al., 2016; Kruijer et al., 2014, 2017). At variance, the presence in chondritic meteorites of Ca-, Al-rich inclusions (CAIs) and chondrules spanning in age the first 4 Myr of the solar system (e.g., Villeneuve et al., 2009) shows that their parent planetesimals completed their accretion significantly later. Because no phenomenological law for planetesimal accretion has ever been derived from data on meteorites, the first principles behind these different accretion times still remain inaccessible.

One way to unravel the accretion history of planetesimals is to study their thermal evolution, which is strongly controlled by their size (e.g., Merk et al., 2002) that sets the ratio between internal production of energy (proportional to the cube of the size) and the surface heat loss (proportional to the square of the size). Previous models have essentially considered that the surface heat loss could be estimated based on a conductive profile in the planetesimal (e.g., Ghosh et al., 2003). It has been shown however that melting-induced convection could occur ≈1 Myr after the accretion of a planetesimal (Hevey and Sanders, 2006). The conductive hypothesis is thus not valid over the whole thermal history of a planetesimal and more efficient heat transfer by convection must be taken into account. To our knowledge, the recent study of Lichtenberg et al. (2016) is the only one to explicitly model convection in the planetesimal, using 2D and 3D numerical simulations. However this study considers instantaneous accretion and does not tackle the question of the thermal evolution of the planetesimal in the route of accretion. We propose to do so in the present study by using scaling laws for the transient thermal evolution of an internally heated convective body validated by lab-
scale experiments.

## 2. Thermal evolution of a planetesimal in the course of accretion: principles and constraints

### 2.1. Conservation of energy

The average internal temperature $T(t)$ of a spherical body of radius $R(t)$ can be defined based on its total thermal energy $E(t)$ as

$$T(t) = \frac{E(t)}{4/3\pi\rho C_p R^3(t)}, \qquad (1)$$

with $\rho$ the density and $C_p$ the specific heat (Table 1). The time evolution of this average internal temperature is given by the balance between internal radioactive heating due to decay of short-lived $^{26}$Al, $H(t)$, and the surface heat loss, $Q_s(t)$,

$$\rho C_p \frac{4}{3} \pi \frac{dR^3(t)T(t)}{dt} = \frac{4}{3}\pi R^3(t)H(t) - 4\pi R^2(t)Q_s(t) + E_{int}, \qquad (2)$$

where $E_{int}=E_{lat}+E_a$ is a term corresponding to the heat absorbed by partial melting ($E_{lat}$) or brought by surface impacts ($E_a$), as discussed later (Section 3.1). For a given accretion law $R(t)$, $T(t)$ can be calculated using the energy balance equation (2) if the surface heat flux $Q_s(t)$ is known. Equivalently, but the other way around, if a thermal record $T(t)$ is given, equation 2 can be inverted to obtain the accretion curve $R(t)$.

In the purely conductive regime $Q_s(t)$ is analytically known (Carslaw and Jaeger, 1959, p. 245). In a convective regime, the determination of $Q_s(t)$ requires either the numerical resolution of the full set of conservative equations in the convective fluid (Lichtenberg et al., 2016), or the use of a
scaling law relating the surface heat flux to the fluid properties, the internal
heating rate, and the geometry of the system; an approach we will detail and
follow in section 3.

## 2.2. Constraints brought by the magmatic iron meteorites

### 2.2.1. Temperature and age of formation

The ¹⁸²Hf-¹⁸²W age of magmatic iron meteorites can be taken as the formation time of their core, or, more generally, the time of metal-silicate differentiation and metal segregation (Lee and Halliday, 1997; Kleine et al., 2002; Qin et al., 2008; Kruijer et al., 2014). Data for five groups of magmatic iron meteorites (IIAB, IID, IIIAB, IVA, IVB) provide ages spanning from 0.7±0.3 Myr after CAIs for IIAB to 3.1±0.8 Myr after CAIs for IIB (Kruijer et al., 2014). Their temperatures of formation, taken as the liquidus temperature of the iron-sulfur alloy, range from 1325 °C to 1615 °C as a function of their sulfur content (see Table 2 from Supplementary material). We note in Figure 1 that a single exponential function fits within errors bars the data on the five groups of magmatic iron meteorites $T(t) = T_\infty[1 - \exp(-t/\tau)]$. We used an Orthogonal Distance Regression (ODR) method based on a Levenberg-Marquardt type algorithm that allows non-linear fitting with errors on both variables, and we obtained the following values: $\tau = 0.424 \pm 0.108$ Myr and $T_\infty = 1524 \pm 46$ °C. Rather than considering a specific thermal evolution hence a somewhat open ended accretion history for each magmatic iron meteorite, such a unique fit opens the possibility to look for a tightly constrained universal accretion scenario that can explain all the data. To that aim, we will use in the following (section 4) the best fitting exponential evolution of the temperature to obtain accretion curves $R(t)$ - if they exist - consistent
109 with this thermal evolution.

110 **2.2.2. Metallographic cooling rates**

111 To further constrain the accretion curves, we will rely on the metallo-
112 graphic cooling rates to infer the size of the parent bodies of magmatic iron
113 meteorites. These correspond to the latest (conductive) cooling phase of a
114 metallic pocket (or core) once it has solidified. They are thus independent
115 from the thermal models considered in this study which will only describe
116 the heating (and convective) stage of the accreting planetesimals. The size
117 of the planetesimals are poorly constrained and have been estimated based
118 on two assumptions (see Table 3 from Supplementary material).
119

120 The IVA meteorites have the highest cooling rates, interpreted as result-
121 ing from the conductive cooling of a 150±50 km metallic core with its mantle
122 ripped off, yielding a parent body size of ≈300 km (Yang et al., 2008). Cool-
123 ing rates for IVB meteorites give a parent body size of 140±30 km (Yang
124 et al., 2010).

124 Other magmatic iron meteorites display smaller cooling rates (Yang and
125 Goldstein, 2006; Goldstein et al., 2009) which reflects the insulating effect
126 of a mantle and possibly of a regolith. Their size depends on the thickness
127 and thermal properties of the rocky envelopes around the core (Haack et al.,
128 1990). Based on the recorded cooling rates and on Figure 7.1 of McSween
129 (1999) we estimate the radii as 125±50 km for IIIAB, 235±65 km for IID,
130 and 165±65 km for IIAB.
### 3. Transient thermal evolution in an internally heated convective system

The resolution of the energy balance equation 2 requires a knowledge of the surface heat flux $Q_s$ which is different in a conductive and in a convective thermal regime. In the case of Rayleigh-Bénard convection, where convection results from cooling at the top of the system and heating at its base, scaling laws are expressed using the Nusselt number, $Nu$, defined as the ratio between the surface heat flux measured in the convective regime and the surface heat flux in the conductive regime. At steady state, $Nu$ scales as the 1/3rd power of the Rayleigh number of the well-mixed convective layer. In internally heated fluids, convection depends on the so called Rayleigh-Roberts number, defined as

$$Ra_H = \frac{\rho g \alpha H R^5}{k \kappa \nu}, \qquad (3)$$

with $g$ the gravity, $\alpha$ the coefficient of thermal expansion, $k$ the thermal conductivity, $\kappa$ the thermal diffusivity and $\nu$ the kinematic viscosity (Table 1). In this setting, at steady state the Nusselt number is unity and the average internal temperature at steady state in an internally heated convective fluid (Parmentier and Sotin, 2000) follows a power scaling law

$$T/T_H = cRa_H^{-\beta}, \qquad (4)$$

with $c$ and $\beta$ dimensionless constants, and $T_H = \rho H R^2/k$ the internal heating temperature scale with $R$ now the fluid thickness in the experiments. Experiments yield $\beta=1/4$ at steady state (Limare et al., 2015), in agreement with theory (Vilella and Kaminski, 2017). In principle, this law holds for the temperature difference across the unstable boundary layer. It turns out,
however, that the temperature at the edge of the boundary layer and the average internal temperature can both be scaled to $Ra_H^{-1/4}$ and hence are proportional to one another (Vilella et al., 2018; Limare et al., 2015). We confirmed that this scaling is valid in our experiments over several orders of magnitude in $Ra_H$ as can be seen in the accompanying MethodsX paper (Limare et al., 2020) with $c=3.45\pm0.15$.

As we consider the transient regime in the experiments, i.e., before the steady state is reached, a secular term $dT/dt$ has to be introduced in the scaling law. Assuming that the secular term can be treated as an additional contribution to internal heating, equation 4 becomes

$$T = cRa_H^{-\beta}\rho\left(H - C_p\frac{dT}{dt}\right)R^2/k, \quad (5)$$

where $Ra_H^*$ is a modified Rayleigh-Roberts number defined as

$$Ra_H^* = \frac{\rho g \alpha (H - C_p \frac{dT}{dt}) R^5}{k \kappa \nu}. \quad (6)$$

The two previous equations can be combined to yield an expression for the time derivative of the average internal temperature,

$$\rho C_p \frac{dT}{dt} = \rho H - \left(\frac{kT}{cR^2}\right)^{\frac{1}{1-\beta}} \left(\frac{g\alpha R^5}{k\kappa\nu}\right)^{\frac{\beta}{1-\beta}}. \quad (7)$$

Further noting that the conservation of energy per unit surface writes

$$\rho C_p R \frac{dT}{dt} = \rho H R - Q_s, \quad (8)$$

we obtain an expression for the surface heat flux in the convective regime,

$$Q_s = \left(\frac{kT}{cR}\right)^{\frac{1}{1-\beta}} \left(\frac{g\alpha R^4}{k\kappa\nu}\right)^{\frac{\beta}{1-\beta}}, \quad (9)$$
that is found in very good agreement with the experimental results over
several orders of magnitude (Limare et al., 2020). We will now apply these
experimentally validated theoretical results to systems where the secular evo-
lution of the average internal temperature is also related to the time evolution
of internal heating rate as in the case of a planetesimal.

3.1. *Adaptation of the scaling-laws to the thermal evolution of a planetesimal*

To apply the experimentally determined scaling law to the modeling of the
thermal evolution of a planetesimal, we need to consider several differences
between the experimental setting and the natural case. First, we multiply
the literal expression we have obtained in equation (9) by a factor 1/3 corre-
sponding to a spherical geometry, as established by Deschamps et al. (2012).
Second, we take into account the decreasing heating rate induced by the pro-
gressive disintegration of <sup>26</sup>Al and the implications on the material properties
of a planetesimal given in Table 1.

In the parent body, viscosity is not constant but changes as a function of
the average internal temperature and – when relevant – the fraction of melt,
and is given by

$$
\mu = \mu_0 f(\phi^{\text{melt}}) \exp \left( \frac{E_{\text{act}}}{R_g T} - \frac{E_{\text{act}}}{R_g T_0} \right), \quad (10)
$$

where $\mu_0$ is a reference viscosity at $T_0=1000$ °C (that will be set between
$10^{18}$ and $10^{21}$ Pa s, see section 4), $R_g$ is the ideal gas constant, $E_{\text{act}}$ is the
activation energy (Table 1), and $f(\phi^{\text{melt}})$ is defined as

$$
f(\phi^{\text{melt}}) = \exp(-\sigma_{\text{sil}}\phi_{\text{sil}}^{\text{melt}}) \exp(-\sigma_{\text{Fe}}\phi_{\text{Fe}}^{\text{melt}}), \quad (11)
$$
where the fractions of melt are taken as a linear function of temperature,

$$ \phi_{sil,Fe}^{melt} = \frac{T - T_{s}^{sil,Fe}}{T_{l}^{sil,Fe} - T_{s}^{sil,Fe}}, \qquad (12) $$

if $T_s^{sil,Fe} < T < T_l^{sil,Fe}$ with $T_s^{sil,Fe}$ and $T_l^{sil,Fe}$ the solidus and liquidus temperature of silicate and metal, respectively, and $\sigma_{Fe}=4$ and $\sigma_{sil}=21$ (Scott and Kohlstedt, 2006; Hustoft et al., 2007).

Because in silicate bodies the viscosity is strongly temperature dependent, we must also account for the presence of a conductive viscous lid at the top of the convective system (e.g. Davaille and Jaupart, 1993). In that case, four additional parameters have to be introduced: the thickness of the lid and the temperature drop across the lid, the thickness of the thermal boundary layer at the top of the convective fluid and the temperature drop across the thermal boundary layer. The thickness and temperature drop across the thermal boundary layer at the top of the convective fluid are given by well established scaling laws for thermal convection with large viscosity variations (Davaille and Jaupart, 1993; Choblet and Sotin, 2000). These laws introduce a “viscous scale” $\Delta T_{visc}$ for the temperature,

$$ \Delta T_{visc}(T) = -2.21 \frac{\mu(T)}{d\mu/dT(T)}. \qquad (13) $$

The thickness $\delta_{lid}$ of the viscous lid and the temperature $T_{lid}$ at its base are then obtained by imposing the continuity of temperature and heat flux at the base of the lid. To that aim we use $T = T_{lid} + \Delta T_{visc}$ as the average internal temperature of the convective fluid and $\delta_{lid} = kT_{lid}/Q_s$, while replacing $R$ by $R - \delta_{lid}$ in equation (9).

The energy absorbed by partial melting, $E_{lat}$, and the energy brought by surface impacts during accretion, $E_a$, are the last ingredients to take into
account in the conservation of energy in the case of a planetesimal. The rate of heat absorption by partial melting is given by

$$E_{lat} = L_{Sil} \frac{d\xi_{Sil}}{dt} + L_{Fe} \frac{d\xi_{Fe}}{dt}, \qquad (14)$$

with $L_{Sil}$ and $L_{Fe}$ the latent heat of the silicate and metal (Table 1), respectively, and $\xi_{Sil}$ and $\xi_{Fe}$ their melt fraction. The rate of variation of the melt fractions for $T_s^{sil,Fe} < T < T_l^{sil,Fe}$ is derived from equation (12) and writes

$$\frac{d\xi^{sil,Fe}}{dt} = \frac{dT/dt}{T_l^{sil,Fe} - T_s^{sil,Fe}}, \qquad (15)$$

with $T_s^{sil}=1150$ °C and $T_s^{Fe}=990$ °C the silicate solidus and the Fe-FeS eutectic, respectively, and $T_l^{sil}=1800$ °C and $T_l^{Fe}=1615$ °C for the silicate liquidus and the pure Fe melting temperature, respectively. To take into account surface impacts, we further consider that an impactor of mass $M_i$ and temperature $T_i$ brings a gravitational energy $E_G = GMM_i/R$, with $M$ and $R$ the mass and radius of the target, respectively, and a thermal energy $E_T = M_iC_pT_i$. The fraction of energy that is transferred to the target rather than radiated into space, $E_a$, is then set to a nominal fraction of 30% of ($E_G + E_T$) similar to the value used by Sramek et al. (2012). The effect of this fraction is discussed in the Supplementary Material Figure 11.

# 4. Accretion curve of the parent bodies of the magmatic iron meteorites from their thermal evolution

## 4.1. Instantaneous accretion

As a first step in the study of the thermal evolution of planetesimals we consider, as in previous models, instantaneous accretion and a purely conductive model. To be consistent with the sizes inferred from the metallographic
cooling rates we calculate the thermal evolution of a 300 km body as a func-
tion of its time of accretion. We further consider a constant temperature of
0 °C at the surface of the planetesimal. The results shown in the inset of
Figure 2 imply that the parent bodies of non-carbonaceous magmatic iron
meteorites accreted earlier (≈0.5 Myr after CAIs formation) than the ones
of carbonaceous magmatic iron meteorites (≈1 Myr after CAIs formation).
These results are consistent with the conclusions obtained for a 40 km body
by Kruijer et al. (2014) which shows that within the framework of instanta-
neous accretion and a purely conductive model the size of the body is not a
key parameter for the thermal evolution.

Conductive models predict extreme temperatures in the planetesimals,
larger than the liquidus of the metal and silicate phases, hence total melt-
ing and differentiation of the non-carbonaceous magmatic iron meteorites
(inset in Figure 2). Hence at least melting-induced convection should be
taken into account in the thermal modeling. Both the increase of tempera-
ture and the occurrence of melting actually decrease the viscosity (Supple-
mentary Material Figure 8) and increase the value of the Rayleigh-Roberts
number in the planetesimal. We obtain that the critical value for convection,
RaH,cr = 5758 for a rigid boundary condition at the surface (Schubert et al.,
2001), is reached at ≈0.1 Myr after CAIs (if accretion also starts at CAIs).
Once RaH ≥ RaH,cr convection should start. However, in a fluid with a
strongly temperature dependent viscosity the onset of convection is signifi-
cantly delayed (Davaille and Jaupart, 1993; Choblet and Sotin, 2000). We
calculate the onset time of convection using the scaling given by experiments
and numerical simulations (Davaille and Jaupart, 1993; Choblet and Sotin,
256 2000),

$$t_0 = \frac{1}{\pi} \frac{R^2}{\kappa} \left( \frac{Ra_{H,cr}}{Ra_H} \right)^{-2/3} \left( \frac{T}{\Delta T_{visc}} \right)^{8/3}. \quad (16)$$

The onset time evolves as a function of time following the increase of temperature in the planetesimal. We consider that convection starts when the age of the body becomes larger than $t_0$ and we obtained an onset time of ≈ 0.2 Myr. This can be taken as an upper bound as convection may start as soon as some melt segregation of iron triggers a density contrast driven Rayleigh-Taylor instability.

Thermal evolution in the case of instantaneous accretion, when convection is taken into account, is shown in Figure 2. Convection increases heat transfer in the body and results in significantly cooler average internal temperatures than in the purely conductive case. But the main conclusions reached in the conductive case still hold: the parent body of non-carbonaceous magmatic iron meteorites accreted earlier and underwent a hotter thermal evolution leading to total melting and complete differentiation. Convection does not rule out the capability of instantaneous accretion models to explain the magmatic iron meteorites recorded ages and temperatures. Instantaneous accretion cannot however yield the common thermal evolution we have proposed Figure 1 and we will look now for accretion scenarios that can produce it.

## 4.2. Continuous pebbles accretion

Using the time evolution of the average internal temperature proposed in Figure 1, we invert the energy balance equation 2 as a function of the thermal regime (conduction or convection) to obtain the accretion curves $R(t)$. In this case $E_{int}=E_{lat}$ since $E_a$ can be considered 0 for pebble accretion. Within
our framework, the free parameter in the inversion routine is the initial size
of the planetesimal, $R_0$, that we will refer to as the nucleus, produced by
gravitational collapse at the early beginning of the accretion. Because this
phenomenon is very rapid compared to the time scale of the thermal evo-
lution recorded by the magmatic iron meteorites, it can be considered as
instantaneous in the present model.

In the first stage of accretion, the radius of the body is small and heat
transfer occurs by conduction. For the sake of the argument we first con-
sider a purely conductive regime to obtain $R(t)$ from the required thermal
evolution (figure 1). The results of the inversion are shown in figure 3 for a
nucleus radius of 35 and 90 km. Based on these results we can calculate the
evolution of $Ra_H$ during the accretion of the body (Supplementary Material
Figure 9) and further consider the influence of convection on the model pre-
dictions. We calculate the onset time of convection along the accretion curve
using the scaling of equation 16. We consider that convection starts (hence
that Equation 9 applies) when the age of the body becomes larger than the
onset time.

Figure 3 shows two examples of accretion curves resulting from the inver-
sion of the time evolution of the average internal temperature once convec-
tion has been introduced. The shape of the accretion curve is the same as in
the case of a purely conductive regime (notably because the initial thermal
regime is always conductive). However, because convection corresponds to
a more efficient heat transfer, it is not possible in this regime to reach the
thermal plateau observed at large times in Figure 1 if cold pebbles continue
to accrete at the surface of the planetesimal.
In some cases it is even necessary that the radius of the body decreases
in order to counterbalance convection-induced cooling. This is due to the
fact that the internal energy varies with volume (which depends on the third
power of radius) and the cooling rate is proportional to the surface (second
power of radius). Due to convection the body might cool too fast, therefore
the model calculates a radius that start to diminish in order to restore the
almost constant temperature. This effect is more important for low values of
initial radius (cooling/internal heat ~ 1/R).

This decrease of radius can correspond either to some erosion, or to sintering at the end of accretion due to compaction of the planetesimal following partial melting. In the following we shall reject the accretion curves yielding to more than 10% of sintering as this is the bulk porosity value we consider in the present model. Based on all the constraints we conclude that the nucleus size of the parent body was between 35 and 90 km (figure 3).

As already mentioned, accretion curves depend only on the initial (nucleus) size $R(t=0) = R_0$. When made dimensionless by dividing the radius of the body by its final value, all acceptable solutions (i.e., fitting both the thermal evolution recorded by the magmatic iron meteorites and their sizes) appear to share the same evolution (figure 4). This implies that both the shape of the accretion curve, and the associated thermal evolution, do not depend on the final size of the body (which itself depends on the initial nucleus radius), and the growth curves in figures 3 follow a single parameterized accretion law (figure 4). Three stages of evolution can be introduced to describe the universal dimensionless accretion curve: (i) an initial quasi-instantaneous accretion of a nucleus with $35 \le R_0 \le 90$ km, too fast to be resolved at the
scale of our model, (ii) a following slower accretion that reaches a plateau
or even ends at ≈3 Myr, and (iii) a sintering due to the compaction of the
planetesimal following partial melting. Compared to the more simple case
of a purely conductive regime, the end of accretion at ≈3 Myr appears well
resolved in the complete model including convection.

334 **5. Discussion**

335 **5.1. Alternative accretion scenarios**

To evaluate further the robustness of our model predictions, we study two
alternative accretion scenarios: (i) successive collisions with impactors with
sizes of tens of kilometers, and (ii) the fragmentation of a parent body.

In scenarios of discontinuous accretion, the planetesimal grow through
successive collisions with impactors of a given size Rᵢ. Here we consider that
the impactors are formed with their final size Rᵢ at t=0 and we follow their
average internal temperature using the conservation of energy as given by
equation 2. We model two different types of accretion: orderly or runaway
(figure 5). The resulting thermal evolution is illustrated figure 5 for a final
parent body radius of 300 km and Rᵢ=40 km or Rᵢ=10 km. None of
these scenarios yield a thermal evolution consistent with the magmatic iron
meteorites record.

The results of scenario of discontinuous accretion – in the case of large
impactors – are a function of the fraction of energy that is transferred to the
parent body during accretion. We study the influence of this parameter for
the thermal evolution corresponding to orderly accretion of impactors of 40
km formed at time 0 producing a final body of 300 km in 4 My (Supple-
mentary Material Figure 11). We obtain that in the case of a perfect energy transfer to the parent body, or in the case of total loss of the impactor energy to space, the thermal evolution predicted is not consistent with the constraints provided by the magmatic iron meteorites. We can then conclude that our preferred scenario of continuous accretion is both robust and the only one able to yield a thermal evolution consistent with the constraints provided by the magmatic iron meteorites.

In the alternative scenario of fragmentation, the planetesimal has a large size at time 0 and is later fragmented into smaller bodies. For the sake of the argument we consider here an initial body size of 300 km, but the conclusions we reach are general. Because of the large initial size, the early thermal evolution is very fast, the surface heat loss ($\propto R^2$) being much smaller than the internal heat sources ($\propto R^3$). As a consequence total melting occurs early, even before the formation of IIAB (Supplementary Material Figure 12). The most favorable end member fragmentation scenario would be a very early fragmentation into small bodies, but this would imply in turn a subsequent episode of continuous accretion and would be equivalent to the formation of a nucleus followed by continuous pebbles accretion.

## 5.2. Sensibility of the model of pebbles accretion with respect to model parameters

Figure 3 shows that taking into account convection in a planetesimal under continuous accretion has a small impact on model prediction, limited to an end of accretion delayed from ≈3 Myr to ≈4 Myr. This reflects the strong control of the model outputs by the thermal evolution of Figure 1. As discussed above, this rather “cold” thermal evolution results mainly from the
378 accretion of cold pebbles, and the accretion curve is only tuned as a function of the thermal regime. This implies that the model predictions are robust 379 and do not change significantly when we change the parameters used in the 380 scaling laws. For example, a change of the reference viscosity by an order 381 of magnitude has an almost negligible effect (Figure 3). We obtain also a 382 negligible effect of a 50 % change in the values of activation energy and of 383 melt sensitivities σ<sub>Fe</sub> and σ<sub>sil</sub>, or for a surface temperature of -150 °C instead 384 of 0 °C (Supplementary Material Figure 10). 385

386 We also allowed the system to differentiate, considering iron segregation 387 and migration towards the center of the parent body to form a metallic core. 388 There are three effects of differentiation for the convecting silicate mantle: 389 (i) a prefactor in the Rayleigh number (Equation 3) equal to (1 + *f* + *f*²)/3, 390 where *f* = *R*<sub>core</sub>/R (Deschamps et al., 2012) with *R*<sub>core</sub> the core radius given 391 by the mass fraction of metal, (ii) a concentration of the heat source, hence 392 an increase of *H*(*t*), in the mantle, as no Alumium is partitioned into the 393 core, and (iii) a decrease of its thickness from *R* to *D* = *R* − *R*<sub>core</sub>. The two 394 last effects tend to counterbalance each other and the resulting bulk impact 395 of the differentiation on the accretion curves is not significant, as can be seen 396 in Figure 6. 397

5.3. *Comparison with astrophysical models*

398 The first two stages of our accretion law are consistent with astrophysical 399 models that predict (i) a fast initial accretion due to gravitational collapse of 400 meter to kilometer size boulders (Johansen et al., 2015; Simon et al., 2017) 401 and (ii) a rather cold and continuous snowball sweeping of dust and pebbles 402 (Xie et al., 2010). The second stage, corresponding to pebble accretion, can
further be interpreted in terms of the evolution of the fraction of pebbles in the disk.

Visser and Ormel (2016) performed numerical simulations of pebble accretion as a function of the sizes of the accreted particles and of the accreting planetesimal as well as of the distance from the star. They define and compute a pebble accretion time scale,

$$t_{growth} = \frac{M}{\dot{M}} = \frac{4\rho_{\bullet}R}{3v_{hw}\rho_{p}f_{coll}}, \qquad (17)$$

with M the mass of the accreting planetesimal, $\dot{M}$ its time derivative, $\rho_{\bullet}=1$ g cm$^{-3}$ the internal density of the pebbles, $v_{hw}=5780$ cm s$^{-1}$ the magnitude of the headwind faced by the planetesimal, $f_{coll}$ the collisional efficiency, and $\rho_{p}$ the density of pebbles in the disk (expressed as a relative density compared to the gas density $\rho_{g}$). The value of $t_{growth}$ depends on the size of the planetesimal, on the solid-to-gas ratio $\rho_{p}/\rho_{g}$, and on the distance from the star. For the sake of the argument, we shall use the results obtained by Visser and Ormel (2016) at 3 AU and for a typical planetesimal size of 100 km which is relevant for our study and corresponds to a minimal dependence of the result on the size of the pebbles. For this choice of parameters, and for a typical solid-to-gas $\rho_{p}^{ref}/\rho_{g} = 1/100$, Visser and Ormel (2016) show in their figure 4 that the growth timescale $t_{growth}^{ref}$ ranges between 40 Myr and 100 Myr for pebbles radii between 0.01 cm and 30 cm. These time scales are one order of magnitude larger than the one estimated by our model. This implies that the solid-to-gas ratio during the accretion of iron meteorites was also one order of magnitude larger.

For a given accretion curve $R(t)$ obtained by the inversion of the thermal evolution recorded by the magmatic iron meteorites it is possible to
calculate the evolution of the accretion time scale $t_{growth}(t) = \frac{1}{3}R(t)/\dot{R}(t)$.
In the conditions described above, the ratio $t_{growth}(t)/t_{growth}^{ref}$ is directly the ratio $\rho_p^{ref}/\rho_p(t)$ with $\rho_p(t)$ the time-varying fraction of pebbles seen by the planetesimal during its accretion. As an example, we show in figure 7 the evolution of $\rho_p(t)/\rho_g$ for a nucleus of 65 km. As expected, the density of pebbles starts with a very high value, consistent with the production of a nucleus by gravitational collapse in overdense zones in the disk, and sharply decreases around 3 Myr at the end of the accretion. One should note however that this marks also the appearance of a new mechanism (sintering and/or erosion) that would artificially decrease the effective density of pebbles “seen” by the planetesimal.

## 5.4. Consequences for the differentiation of the parent body of magmatic iron meteorites

Under canonic nebular conditions, iron can be anticipated to be present within these pebbles and dust as metallic Fe (dominant form), FeS and oxidized Fe in silicates (Ebel, 2006). Upon heating within the planetesimal, it is expected that FeS melts first at the corresponding Fe-FeS eutectic ≈1000°C (Brett and Bell, 1969) (see Supplementary Material, Figure 13). As the temperature rises, the melt becomes more iron rich. As Kruijer et al. (2014) we further assume that the metallic melt remains in contact with the silicate matrix until the total melting of FeS alloy. Under this assumption, the recorded temperature corresponds to the liquidus temperature set by the Sulfur content and the Hf-W age to the time at which this liquidus temperature is reached. This implies that the systematic difference in age between the old “non-carbonaceous” (NC) magmatic iron meteorites and the younger
452

"carbonaceous" (CC) magmatic iron meteorites, is not due to a difference in accretion age, but simply related to the fact that NC are richer in S, and thus reached their FeS liquidus temperature earlier than CC.

453

454

455

It is worth noting that metallic melt separation that occurs at liquidus temperature does not necessarily correspond to total differentiation, thus the Hf-W age does not correspond to the age of core formation, but rather to the formation of metal pockets large enough to prevent metal-silicate W re-equilibration. Eutectic melting of Fe-FeS can take place much below the temperature of melting of the silicate framework which may indeed hamper the percolation of the Fe-FeS liquids to the core. The problem of the differentiation of an accreting planetesimal has been previously addressed using a purely conductive model, pure Fe metal only and no surface tension (Sramek et al., 2012) for a 500 km planetesimal that accreted within 3 Myr. Their results can be used to discuss the process of compaction and core formation in our scenario. At 3 My, about half of the metal is molten and has segregated to form a central core and a non negligible part of the molten metal is furthermore trapped in the solid silicate matrix as lenses at intermediate depths. It is very likely that metal extraction will be even less efficient in our model for two reasons. First, all things being equal, the average internal temperature in a convective planetesimal (our case) is smaller than in a conductive one (the case studied in Sramek et al. (2012)). Hence the rate of production of molten metal, and the subordinated rate of metal extraction, will be smaller in the convective than in the conductive regime. As a consequence differentiation will take a longer time in a convective planetesimal and will have reached a less advanced stage when cooling and solidification

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476
will occur. Second, it has been theoretically established that high surface tension (such as the one between solid silicate and Fe-FeS metal) enhances local melt concentration rather than global melt extraction (Ricard et al., 2001). Experiments have further shown that a threshold fraction of about 15 vol% of Fe-S melt was required to overcome surface tension effects and to extract the melt from the solid silicate matrix (Bagdassarov et al., 2009). We can thus conclude that in our model complete melting of the Fe-FeS mixture is a necessary condition for an efficient compaction leading to the formation of a core. However, melting of more than 50% of the silicate matrix will be further necessary to yield complete metal extraction and full differentiation (Sramek et al., 2012), as in that case percolation is not limited by the rheology of the matrix anymore.

Based on the previous arguments, we conclude that the parent bodies of magmatic iron meteorites might be only partially differentiated, some iron meteorites being possibly fragments of local ponds of molten metal rather than pieces of the central core (figure 4). This conclusion is in fact not specific to the present accretion model. In the models considered up to now for the formation of the parent bodies of iron meteorites, i. e., instantaneous accretion at given times (Kruijer et al., 2014), the thermal evolution calculated taking into account convection predicts much lower average internal temperatures than previously anticipated in purely conductive models (Figure 2). In a scenario of instantaneous accretion, total melting and full differentiation occurred only for the “non carbonaceous” magmatic iron meteorites.

The present accretion model, while developed to reproduce the thermal histories of the parent bodies of magmatic iron meteorites, seems to explain
also correctly two major characteristics of the parent bodies of CV carbonaceous chondrites, namely (i) the age of the end of their accretion and (ii) the coexistence of an external undifferentiated shell with a central convective metallic core potentially powering a dynamo (Gattacceca et al., 2016). In our scenario, convection develops in the planetesimal underneath a stable conductive viscous lid of a few kilometer thickness (Supplementary Material Figure 14). The interior of this viscous lid is by definition much colder than the interior of the body hence cannot reach the solidus and remains undifferentiated. The distribution of $^{26}$Al ages in chondrules from CV chondrites implies that accretion of their parent body was still ongoing ≈4 Myr after CAIs (Villeneuve et al., 2009). In our model this would correspond to the very last stage of limited accretion occurring between ≈ 3 and 4 Myr, once sintering has become the dominant feature. At 4 Myr, the 300°C isotherm in the viscous lid is located ≈2-3 km below the surface, allowing the preservation of primitive accreted materials including chondrules and of undisturbed $^{26}$Al isochrons in these chondrules (Supplementary Material Figure 14). Such a planetesimal would have developed a metallic core which could remain partially melted long enough for a dynamo to develop in agreement with recent paleomagnetic observations in chondrules from CV chondrites (Gattacceca et al., 2016). It is worthwhile noting that the presence of a stagnant lid and of internal heating reduce the amount of heat that can be extracted from the core and transported by convection to the surface. In some cases, if the internal heating is too high, the heat flux at the boundary between the core and the mantle may even be negative, i.e. the core is heating up, which would then help to keep it molten for longer periods of time.
## 527 **6. Conclusion**

We have developed a model for the transient thermal evolution of a planetary in the course of accretion. Compared to previous approaches, the model takes into account both a conductive regime (at the beginning of accretion when temperature increase and radius are small) and a convective regime. Although convection is not necessarily required by our model, for nominal values of the parameters that enter the experimentally derived scaling laws, convection would start around 1 Myr during the accretion.

We show that the thermal evolution recorded by all the magmatic iron meteorites can be described by a unique exponential law that we latter employed to obtain accretion curves using the conservation of energy in the parent bodies of the magmatic iron meteorites. Once made dimensionless, all the accretion curves define a universal accretion law that introduces three stages: (i) instantaneous nucleus formation at $t=0$ by gravitational collapse, (ii) continuous pebble accretion at a rate decreasing with time and ending at $\approx 3$ Myr, (iii) negligible accretion accompanied by a few percent sintering by compaction.

A main difference of our model compared to previous predictions of models that considered instantaneous accretion and a purely conductive thermal regime is a colder thermal evolution controlled by the continuous accretion of cold “pebbles”, and that may yield partial differentiation only. Some of the magmatic iron meteorites may thus be local metal ponds rather than the remnant of the metallic core of their parent body.

Furthermore, regardless the cooling mode (conduction or convection) in the parent body, there will always exist a cold viscous stagnant lid of a few
552 kilometers that should have remained completely undifferentiated.

553 On another hand the presence of a stagnant lid may also delay the cooling
554 of the planetesimal even if cooling proceeds by convection. If such a plan-
555 etesimal is differentiated, then the metallic core can develop a dynamo for a
556 non-negligible amount of time.

557 Partial melting of the silicate phase will produce a certain degree of sil-
558 icate differentiation in addition to the formation of an iron core. Hence
559 multi-phase flow modeling will be required to further describe the thermal
600 evolution of the planetesimal during that stage. Such a complex modeling
601 remains however beyond the scope of this study.

562 More generally, our formalism which combines convective thermal trans-
563 fer and the evolution of the size of the body during its formation, could be
564 used to compute the thermal evolution produced by a given accretion law,
565 hence to provide additional constraints for theoretical models of planetesi-
566 mals formation and evolution.

567 Acknowledgment

568 We would like to thank the reviewers and the editor William B. McK-
569 innon for their fruitful comments that improved the manuscript. This work
570 was funded by the ANR-11-IS04-0004 project and supported by the Tel-
571 lus/PNP program of INSU-CNRS. Many thanks to Medi@Terre for improv-
572 ing Figure 4. This project contributes to the IdEx project Université de Paris
573 ANR-18-IDEX-0001.
574 References

575 Bagdassarov, N., Golabek, G. J., Solferino, G., Schmidt, M. W., 2009. Constraints on the FeS melt connectivity in mantle silicates from electrical 576 impedance measurements. Phys. Earth Planet. Int. 177, 139-146. 577

578 Birnstiel, T., Fang, M., Johansen, A., 2016. Dust evolution and the formation 579 of planetesimals. Space Sci. Rev. 205, 41-75.

580 Brett, R., Bell, P. M., 1969. Melting relations in the Fe-rich portion of the 581 system Fe-FeS at 30 kb pressure. Earth Planet. Sci. Lett. 6, 479-482.

582 Budde, G., Burkhardt, C., Brennecka, G. A., Fischer-Gödde, M., Kruijer, 583 T. S., Kleine, T., 2016. Molybdenum isotopic evidence for the origin 584 of chondrules and a distinct genetic heritage of carbonaceous and non- 585 carbonaceous meteorites. Earth and Planetary Science Letters 454, 293- 586 303.

587 Carslaw, H. S., Jaeger, J. C., 1959. Conduction of heat in solids. Oxford 588 University Press, New York.

589 Choblet, G., Sotin, C., 2000. 3D thermal convection with variable viscosity: 590 can transient cooling be described by a quasi-static scaling law? Phys. 591 Earth Planet. Int. 119, 321-336.

592 Cuzzi, J. N., Hogan, R. C., Shariff, K., 2008. Toward planetesimals: dense 593 chondrules clumps in the protoplanetary nebula. Astophys. J. 687, 1432- 594 1447.
595

Davaille, A., Jaupart, C., 1993. Transient high-Rayleigh-number thermal convection with large viscosity variations. J. Fluid Mech. 253, 141–166.

596

597

Delbo', M., Walsh, K., Avdellidou, C., Morbidelli, A., 2017. A major asymmetric dust trap in a transition disk. Science 357, 1026–1029.

598

599

Deschamps, F., Yao, C., Tackley, P. J., Sanchez-Valle, C., 2012. High Rayleigh number thermal convection in volumetrically heated spherical shells. J. Geophys. Res. 117, E09006.

600

601

602

Ebel, D. S., 2006. Condensation of rocky materials in astrophysical environments. In: Lauretta, D. S., McSween, H. Y. J. (Eds.), Meteorites and the early Solar System II. University of Arizona Press, Tucson, pp. 253–277.

603

604

605

Gattacceca, J., Weiss, B. P., Gounelle, M., 2016. New constraints on the magnetic history of the cv parent body and the solar nebula from the Kaba meteorite. Earth Planet. Sci. Lett. 455, 166–175.

606

607

608

Ghosh, A., Weidenschilling, S. J., McSween Jr, H. Y., 2003. Importance of the accretion process in asteroid thermal evolution: 6 Hebe as an example. Meteorit. Planet. Sci. 38, 711–724.

609

610

611

Goldstein, J. I., Scott, E. R. D., Chabot, N. L., 2009. Iron meteorites: Crystallization, thermal history, parent bodies, and origin. Chem Erde - Geochem. 69, 293–325.

612

613

614

Haack, H., Rasmussen, K. L., Warren, P. H., 1990. Effects of regolith megaregolith insulation on the cooling histories of differentiated asteroids. J. Geophys. Res. 95, 5111–5124.

615

616
617

Hevey, P. J., Sanders, I. S., 2006. A model for planetesimal meltdown by 26Al and its implications for meteorite parent bodies. Meteorit. Planet. Sci. 41, 95–106.

618

619

620

Hustoft, J., Scott, T., Kohlstedt, D. L., 2007. Effect of metallic melt on the viscosity of peridotite. Earth Planet Sci. Lett. 260, 355–360.

621

622

Johansen, A., Lambrechts, M., 2017. Forming planets vis pebble accretion. Annu. Rev. Earth Planet. Sci. 12451, 359–387.

623

624

Johansen, A., Mac Low, M. M., Lacerda, P., Bizzarro, M., 2015. Growth of asteroids, planetary embryos, and Kuiper belt objects by chondrule accretion. Sci. Adv. 1:e1500109.

625

626

627

Kleine, T., Münker, C., Mezger, K., Palme, H., 2002. Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry. Nature 418, 952–955.

628

629

630

Kruijer, T. S., Burkhardt, C., Budde, G., Kleine, T., 2017. Age of jupiter inferred from the distinct genetics and formation times of meteorites. Proc. Natl. Acad. Sci. USA 114, 6712–6716.

631

632

633

Kruijer, T. S., Touboul, M., Fischer-Gödde, M., Bermingham, K. R., Walker, R. J., Kleine, T., 2014. Protracted core formation and rapid accretion of protoplanets. Science 344, 1150–1154.

634

635

636

Lee, D. C., Halliday, A. N., 1997. Core formation on Mars and differentiated asteroids. Nature 388, 854–857.

637
638

Lichtenberg, T., Golabek, G. J., Gerya, T. V., Meyer, M. R., 2016. The effects of short-lived radionuclides and porosity on the early thermo-mechanical evolution of planetesimals. Icarus 274, 350–365.

639

640

641

Limare, A., Kenda, B., Kaminski, E., Surducan, E., Surducan, V., Neamtu, C., 2020. Transient lab-scale experiments of internally heated convection. submitted to MethodsX.

642

643

644

Limare, A., Vilella, K., Di Giuseppe, E., Farnetani, C. G., Kaminski, E., Surducan, E., Surducan, V., Neamtu, C., Fourel, L., Jaupart, C., 2015. Microwave-heating laboratory experiments for planetary mantle convection. J. Fluid Mech. 777, 50–67.

645

646

647

648

McSween, H. Y. J., 1999. Meteorites and their parent planets. Cambridge University Press, Cambridge.

649

650

Merk, R., Breuer, D., Spohn, T., 2002. Numerical modeling of 26Al-induced melting of asteroids considering accretion. Icarus 159, 183–191.

651

652

Morbidelli, A., Raymond, S. N., 2016. Challenges in planet formation. J. Geophys. Res. Planets 121, 1962–1980.

653

654

Ormel, C. W., Klahr, H. H., 2010. The effect of gas drag on the growth of protoplanets. analytical expressions for the accretion of small bodies in laminar disks. Astron. Astrophys. 520, A43.

655

656

657

Parmentier, E. M., Sotin, C., 2000. Three-dimensional numerical experiments on thermal convection in a very viscous fluid: Implications for the dynamics of a thermal boundary layer at high Rayleigh number. Phys. Fluids 12 (3), 609–617.

658

659

660
Qin, L., Dauphas, N., Wadhwa, M., Masarik, J., Janney, P. E., 2008. Rapid accretion and differentiation of iron meteorite parent bodies inferred from $^{182}$Hf-$^{182}$W chronometry and thermal modeling. Earth Planet. Sci. Lett. 273, 94–104.

Ricard, Y., Bercovici, D., Schubert, G., 2001. A two-phase model for compaction and damage 2. Applications to compaction, deformation, and the role of interfacial surface tension. J. Geophys. Res. 106, 8907–8924.

Schubert, G., Turcotte, D. L., Olson, P., 2001. Mantle convection in the Earth and planets. Cambridge University Press, Cambridge.

Scott, T., Kohlstedt, D. L., 2006. The effect of large melt fraction on the deformation behavior of peridotite. Earth Planet Sci. Lett. 216, 177–187.

Simon, J. P., Armitage, P. J., Youdin, A. N., Rixin, L., 2017. Evidence for universality in the initial planetesimal mass function. Astrophys. J. Lett. 847, L12.

Sramek, O., Milelli, L., Richard, Y., Labrosse, L., 2012. Thermal evolution and differentiation of planetesimals and planetary embryos. Icarus 217, 339–354.

Suzuki, A., Ohtani, E., Kato, T., 1998. Density and thermal expansion of a peridotite melt at high pressure. Phys. Earth Planet. Int. 107, 53–61.

van der Marel, N., et al., 2013. A major asymmetric dust trap in a transition disk. Science 340, 1199–1202.
682 Vilella, K., Kaminski, E., 2017. Fully determined scaling laws for volumetrically heated convective systems, a tool for assessing habitability of exoplanets. Phys. Earth Planet. Int. 266, 18–28.

685 Vilella, K., Limare, A., Jaupart, C., Farnetani, C., Fourel, L., Kaminski, E., 2018. Fundamentals of laminar free convection in internally heated fluids at values of the Rayleigh–Roberts number up to $10^9$. J. Fluid Mech. 846, 966–998.

689 Villeneuve, J., Chaussidon, M., Libourel, G., 2009. Homogeneous distribution of Al-26 in the solar system from the Mg isotopic composition of chondrules. Science 325, 985–988.

692 Visser, R. G., Ormel, C. W., 2016. On the growth of pebble-accreting planets. Astron. Astrophys. 586, A66, 1–11.

694 Weidenschilling, S. J., 1977. Aerodynamics of solid bodies in the solar nebula. Mon. Not. R. Astron. Soc. 180, 57–70.

696 Xie, J. W., Payne, M. J., Thébault, P., Zhou, J. L., Ge, J., 2010. From dust to planetesimals: the snowball phase? Astrophys. J. 724, 1153–1164.

698 Yang, J., Goldstein, J. I., 2006. Metallographic cooling rates of the IIIAB iron meteorites. Geochim. Cosmochim. Acta 70, 3197–3215.

700 Yang, J., Goldstein, J. I., Michael, J. R., Kotula, P. G., Scott, E. R. D., 2010. Thermal history and origin of the IVB iron meteorites and their parent body. Geochim. Cosmochi. Acta 74, 4493–4506.
703 Yang, J., Goldstein, J. I., Scott, E. R. D., 2008. Metallographic cooling rates

704 and origin of IVA iron meteorites. Geochim. Cosmochim. Acta 72, 3043–3061.
705
<caption>Table 1: Definition and values of the model parameters or formulae used to compute them.</caption>

References: + Sramek et al. (2012), ° Suzuki et al. (1998), * www.nucleide.org

<table><thead><tr><th>Quantity</th><th>Symbol</th><th>Units</th><th>Value / Formula</th></tr></thead><tbody><tr><td>Metal mass<sup>+</sup>, volume fraction</td><td>x<sup>met</sup>, φ<sup>met</sup></td><td>-</td><td>x<sup>met</sup>=0.18</td></tr><tr><td>Silicate mass<sup>+</sup>, volume fraction</td><td>x<sup>sil</sup>, φ<sup>sil</sup></td><td>-</td><td>x<sup>met</sup>=0.82</td></tr><tr><td>Average porosity</td><td>φ</td><td>-</td><td>0.1</td></tr><tr><td>Density of metal<sup>+</sup></td><td>ρ<sup>met</sup></td><td>kg m<sup>-3</sup></td><td>7800</td></tr><tr><td>Density of silicate<sup>+</sup></td><td>ρ<sup>sil</sup></td><td>kg m<sup>-3</sup></td><td>3200</td></tr><tr><td>Average density</td><td>ρ</td><td>kg m<sup>-3</sup></td><td>1/ρ = x<sup>met</sup>/ρ<sup>met</sup> + x<sup>sil</sup>/ρ<sup>sil</sup></td></tr><tr><td>Specific heat of metal<sup>+</sup></td><td>C<sub>p</sub><sup>met</sup></td><td>J K<sup>-1</sup> kg<sup>-1</sup></td><td>450</td></tr><tr><td>Specific heat of silicate<sup>+</sup></td><td>C<sub>p</sub><sup>sil</sup></td><td>J K<sup>-1</sup> kg<sup>-1</sup></td><td>1200</td></tr><tr><td>Average specific heat</td><td>C<sub>p</sub></td><td>J K<sup>-1</sup> kg<sup>-1</sup></td><td>C<sub>p</sub> = x<sup>met</sup>C<sub>p</sub><sup>met</sup> + x<sup>sil</sup>C<sub>p</sub><sup>sil</sup></td></tr><tr><td>Latent heat of metal<sup>+</sup></td><td>L<sup>met</sup></td><td>kJ kg<sup>-1</sup></td><td>250</td></tr><tr><td>Latent heat of silicate<sup>+</sup></td><td>C<sub>p</sub><sup>sil</sup></td><td>kJ kg<sup>-1</sup></td><td>500</td></tr><tr><td>Activation energy metal</td><td>E<sub>act</sub><sup>met</sup></td><td>kJ mol<sup>-1</sup></td><td>13</td></tr><tr><td>Activation energy silicate</td><td>E<sub>act</sub><sup>sil</sup></td><td>kJ mol<sup>-1</sup></td><td>250</td></tr><tr><td>Average activation energy</td><td>E<sub>act</sub></td><td>kJ mol<sup>-1</sup></td><td>E<sub>act</sub> = φ<sup>met</sup>E<sub>act</sub><sup>met</sup> + φ<sup>sil</sup>E<sub>act</sub><sup>sil</sup></td></tr><tr><td>Thermal conductivity of metal<sup>+</sup></td><td>k<sup>met</sup></td><td>J s<sup>-1</sup> m<sup>-1</sup> K<sup>-1</sup></td><td>50</td></tr><tr><td>Thermal conductivity of silicate<sup>+</sup></td><td>k<sup>sil</sup></td><td>J s<sup>-1</sup> m<sup>-1</sup> K<sup>-1</sup></td><td>3</td></tr><tr><td>Average thermal conductivity</td><td>k</td><td>J s<sup>-1</sup> m<sup>-1</sup> K<sup>-1</sup></td><td>k = φ<sup>met</sup>k<sup>met</sup> + φ<sup>sil</sup>k<sup>sil</sup></td></tr><tr><td>Thermal diffusivity</td><td>κ</td><td>m<sup>2</sup> s<sup>-1</sup></td><td>κ = k/(ρC<sub>p</sub>)</td></tr><tr><td>Dynamic viscosity</td><td>μ</td><td>Pa s</td><td></td></tr><tr><td>Kinematic viscosity</td><td>ν</td><td>m<sup>2</sup> s<sup>-1</sup></td><td>ν = μ/ρ</td></tr><tr><td>Coefficient of thermal expansion<sup>°</sup></td><td>α</td><td>K<sup>-1</sup></td><td>3 10<sup>-5</sup></td></tr><tr><td>Acceleration of gravity</td><td>g</td><td>m s<sup>-2</sup></td><td>g = 4/3πRGρ</td></tr><tr><td>Universal gravity constant</td><td>G</td><td>m<sup>3</sup> kg<sup>-1</sup> s<sup>-2</sup></td><td>6.674 10<sup>-11</sup></td></tr><tr><td>Initial radiogenic heat production</td><td>H<sub>0</sub></td><td>J s<sup>-1</sup> kg<sup>-1</sup></td><td>1.5 10<sup>-7</sup></td></tr><tr><td><sup>26</sup>Al decay constant<sup>*</sup></td><td>λ</td><td>s<sup>-1</sup></td><td>3 10<sup>-14</sup></td></tr><tr><td>Radiogenic heat production</td><td>H</td><td>J s<sup>-1</sup> kg<sup>-1</sup></td><td>H(t) = H<sub>0</sub> exp(-λt)</td></tr><tr><td>Time since CAI formation</td><td>t</td><td>s (Myr)</td><td></td></tr><tr><td>Parent body average internal temperature</td><td>T</td><td>°C</td><td></td></tr><tr><td>Parent body surface temperature</td><td>T<sub>S</sub></td><td>°C</td><td>0</td></tr><tr><td>Parent body surface heat flux</td><td>Q<sub>S</sub></td><td>J s<sup>-1</sup> m<sup>-2</sup></td><td></td></tr><tr><td>Parent body radius</td><td>R</td><td>m</td><td></td></tr><tr><td>Initial parent body radius</td><td>R<sub>0</sub></td><td>m</td><td>R<sub>0</sub> = R(t = 0)</td></tr><tr><td>Rayleigh-Roberts number</td><td>Ra<sub>H</sub></td><td>-</td><td>Ra<sub>H</sub> = ρgαHR<sup>5</sup>/kκν</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1: Thermal evolution of the parent body of the iron magmatic meteorites, $T(t) = T_\infty[1 - \exp(-t/\tau)]$, with $T_\infty=1524$°C and $\tau=0.424$ Myr. The red circles give the position of the “non-carbonaceous” magmatic iron meteorites and the blue circles of the “carbonaceous” magmatic iron meteorites (Kruijer et al., 2014; Budde et al., 2016).</figcaption></figure>
<figure><img src="image_3.png" /></figure>

Figure 2: Evolution of the average internal temperature of a 300 km radius planetesimal as a function of its time of instantaneous accretion, and for two values of the reference viscosity $\mu_0$. The red circles give the position of the “non-carbonaceous” magmatic iron meteorites and the blue circles of the “carbonaceous” magmatic iron meteorites. Compared to a purely conductive model (dashed lines in the inset, with colors corresponding to the same age of formation as in the main figure) the average internal temperature is much lower when convection is taken into account. This is due to more effective heat transfer to the surface in the convective regime than in the conductive one.
<figure><img src="image_4.png" /></figure>

Figure 3: Model results for continuous pebble accretion. The red circles give the position of the “non-carbonaceous” magmatic iron meteorites and the blue circles of the “carbonaceous” magmatic iron meteorites (Kruijer et al., 2014; Budde et al., 2016; Yang et al., 2008, 2010); time is given by the Hf-W age. The accretion curves are a function of the initial (nucleus) radius of the parent body. The green and magenta correspond to a nucleus size of 35 km and 90 km, respectively. The dash-dotted lines give the purely conductive solution, while the dashed and solid lines give the complete solutions with a convective stage for a reference viscosity $\mu_0$ of $10^{19}$ and $10^{18}$ Pa s, respectively. The stars indicate the onset on convection for a reference viscosity of $10^{18}$ Pa s. For nucleus radii smaller than ≈ 35 km the thermal evolution implies a too large sintering of the parent body and is not consistent with our model framework.
<figure><img src="image_5.png" /><figcaption>Figure 4: Universal dimensionless accretion curve predicted by the model (colored thick curve) compared to the parameterized accretion law (dashed white line) of the form $(R(t)/R_\infty)^3 = 1 - (1 - [R_0/R_\infty]^3) \exp(-[t/t_a]^b)$, where $t_a = 1.39\pm0.15$ Myr is the timescale of accretion, $b=1.45\pm0.15$ is a shape parameter, and $R_0$ and $R_\infty$ the initial and final radius, respectively, with $R_0/R_\infty=0.37\pm0.04$. The color of the curves gives the average internal temperature reached in the parent body at a given time (or equivalently for a given radius). The thickness of the colored curve corresponds to the second order variation of the accretion history for different sizes of the initial nucleus.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Alternative accretion scenarios. Upper panel: evolution of the radius of the planetesimal as a function of time for the two types of accretion considered: orderly (blue curve) and runaway (red). Lower panel: corresponding thermal evolution. The thin black line gives the evolution of the temperature of the impactor in the case of orderly accretion. (a,c) Results for 40 km impactors: the initial thermal evolution is too fast and leads to massive melting not consistent with the geochemical constraints. (b,d) Results for 10 km impactors: orderly accretion generates a scenario consistent with the parent bodies of NC iron magnetic meteorites only.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of the accretion curves when differentiation is taken into account (dotted lines) compared to the reference evolution without differentiation (solid lines). In these calculation differentiation occurs instantaneously when the solidus of the silicate is reached (T=1150°C). Differentiation produces a larger sintering at the end of the accretion in small bodies.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

Figure 7: Evolution of the density of pebbles in the disk at 3 AU for a nucleus of 65 km that grows by pebble accretion following the accretion curve based on the inversion of the thermal record of magmatic iron meteorites for a reference viscosity $\mu_0=10^{19}$ Pa s. Two values of the pebbles radius $r_p$ are considered, 0.01 cm and 30 cm. The pebbles density is very high at the beginning of the accretion – in agreement with the hypothesis of boulders accretion in dense parts of the disk – and sharply decreases to values $<10^{-2}$ by the end of accretion at $\approx$ 3 Myr after CAIs. The two steps in the solid line are due to the sharp decrease of viscosity in the planetesimal when iron and then silicate start to melt.
706 **Supplementary Material**

Table 2 contains the Hf-W age of formation, Sulfur content and temper-
ature of formation based on the FeS liquidus for the iron meteorites (Krui-
jer et al., 2014). As in Supplementary Material of Kruijer et al. (2014)
(www.sciencemag.org/content/344/6188/1150/suppl/DC1) we assume that
the metallic melt remains in contact with the silicate matrix until the total
melting of FeS alloy. Under this assumption, the recorded temperature cor-
responds to the liquidus temperature set by the Sulfur content and the Hf-W
age to the time at which this liquidus temperature is reached.
<table><caption>Table 2: Temperature and age of formation of iron meteorites. Data from Supplementary Material of Kruijer et al. (2014).</caption><thead><tr><td>Type</td><td>Age of formation(Myr)</td><td>Sulfur content (wt%)</td><td>T (°C)</td></tr></thead><tbody><tr><td>IIAB</td><td>0.7 ± 0.3</td><td>17</td><td>1325 ± 80</td></tr><tr><td>IIIAB</td><td>1.2 ± 0.3</td><td>12</td><td>1400 ± 80</td></tr><tr><td>IVA</td><td>1.4 ± 0.5</td><td>6</td><td>1445 ± 85</td></tr><tr><td>IVB</td><td>2.9 ± 0.5</td><td>0</td><td>1615 ± 85</td></tr><tr><td>IID</td><td>3.1 ± 0.7</td><td>7</td><td>1445 ± 85</td></tr></tbody></table>
<table><caption>Table 3: Typical size of meteorites estimated from metallographic cooling rates</caption><thead><tr><td>Type</td><td>Radius (km)</td><td>Method</td><td>Reference</td></tr></thead><tbody><tr><td>IIAB</td><td>165 ± 65</td><td>core+mantle</td><td>McSween (1999)</td></tr><tr><td>IIIAB</td><td>125 ± 50</td><td>core+mantle</td><td>Yang and Goldstein (2006)</td></tr><tr><td>IVA</td><td>300 ± 50</td><td>metallic core ripped of its mantle</td><td>Yang et al. (2008)</td></tr><tr><td>IVB</td><td>140 ± 30</td><td>metallic core ripped of its mantle</td><td>Yang et al. (2010)</td></tr><tr><td>IID</td><td>235 ± 65</td><td>core+mantle</td><td>McSween (1999)</td></tr></tbody></table>

Metallographic cooling rates were used to infer the size of the parent bodies of magmatic iron meteorites. These correspond to the latest (conductive) cooling phase of a metallic pocket (or core) once it has solidified. The size of the planetesimals have been estimated based on two assumptions (Table 3): either considering the core with its mantle ripped off, or the cooling of the core still surrounded by its mantle.
<figure><img src="image_9.png" /></figure>

Figure 8: Evolution of the viscosity in the parent body as a function of time for the thermal evolution of figure 1, calculated using equation 11 and for a reference viscosity $\mu_0 = 10^{18}$ Pa s. The progressive decrease of viscosity in due to (i) the progressive increase of average internal temperature in the body and (ii) the occurrence of Fe-FeS and silicate melting and following increase of melt fraction.

Based on the temperature evolution defined in Figure 1, the temperature variation of the viscosity translates into the evolution of viscosity with time shown in Figure 8.
<figure><img src="image_10.png" /></figure>

Figure 9: Evolution of the onset time of convection (left) and of the Rayleigh-Roberts number (right) in the parent body as a function of time for the accretion scenarios presented in figure 3 and for a reference viscosity $\mu_0=10^{18}$ Pa s (solid line) and $10^{19}$ Pa s (dashed line) and two nucleus radii, 35 km (green curves) and 90 km (magenta curves). Stars indicate onset of convection for $\mu_0=10^{18}$ Pa s.

The onset time of convection was calculated using equation 16, valid for a fluid with strongly temperature depending viscosity. We considered that convection starts when the age of the body becomes larger than the onset time. The corresponding Rayleigh-Roberts numbers are much larger than the critical value for rigid boundary conditions at the surface ($Ra_{H,cr}=5758$).
<figure><img src="image_11.png" /></figure>

Figure 10: Accretion scenarios obtained from the inversion of the magmatic iron meteorites record considering a surface temperature of -150 °C (dashed lines) rather than 0 °C (solid lines). Thermal evolution of the parent bodies (left panel). The corresponding accretion curves obtained for a model that included convection ($\mu_0 = 10^{18}$ Pa s) (right panel).

We studied the influence of the surface temperature on the thermal evolution of the parent bodies of meteorites (Figure 10, left panel). Only negligible differences appear in the accretion curves (Figure 10, right panel).
<figure><img src="image_12.png" /></figure>

Figure 11: Illustration of the relative contribution of the heat consumed by partial melting, $E_{lat}$, and of the heat brought by surface impacts, a fraction $F$ of this heat being transferred to the target; violet curve: $F=1$, blue curve: $F=0.3$, green curve: $F=0$. The black curve represents the accretion curve where $E_{int}=0$.

Two sources of energy have to be added to the conservation energy equation 2: $E_{int}=E_{lat}+E_a$, a term corresponding to the heat absorbed by partial melting ($E_{lat}$) or brought by surface impacts ($E_a$). We studied the influence of these energy terms, for an alternative scenario of discontinuous accretion: orderly accretion of 30km-size impactors that produces a parent body of 300 km in 4 Myr. The results are shown in Figure 11. This scenario produces a peak temperature more marked than in the scenario of continuous pebble/dust accretion. A fraction of heat transferred larger than ≈40% is required to fit the NC iron magnetic meteorites constraints, but this implies
a too high and too early peak temperature to fit the IVB and IID constraints.

In that scenario the IVB and IID would have been indeed fully molten and differentiated at ≈2 Myr and this age would have been recorded by the Hf/W chronometer. The contribution of the latent heat ($E_{lat}$) is shown by the difference between the green curve and the black curve ($E_{lat}=0$ for the black curve) and is of second order only. Note that the decrease of temperature at large time is due to efficient heat transfer both in the growing planetesimal and in the impactors whose temperature becomes smaller and smaller as a function of time.
<figure><img src="image_13.png" /></figure>

Figure 12: Thermal evolution of a parent body with an initial radius of 300 km fragmented into 27 bodies of 100 km radius at $t=0.4$ Ma (and using $\mu_0 = 10^{18}$ Pa s). The magenta line corresponds to the thermal evolution of the 100 km fragments and the green line to that of the residual body.

The alternative scenario of fragmentation of a parent body with an initial radius of 300 km is shown in Figure 12. Because of the initial large size of the body an early total melting occurs which is inconsistent with the ages of the magmatic iron meteorites.
<figure><img src="image_14.png" /></figure>

Figure 13: Evolution of the temperature of melting of the metallic phase at 30 kb in the parent body as a function of its sulfur content (Brett and Bell, 1969). The first melt produced is the eutectic ≈74 wt% of S, whereas the temperature of the last melt produced will be set by the total amount of S in the metallic phase (e.g., 1450 °C for 5 wt% of S).

We reproduce here (Figure 13) the phase diagram of the Fe-FeS system at a pressure of 30kb obtained by Brett and Bell (1969) and representing the data used by Kruijer et al. (2014) for the determination of the temperature recorded by the iron meteorites.
<figure><img src="image_15.png" /></figure>

Figure 14: Evolution of the viscous lid thickness and basal temperature once convection has started for the accretion scenario presented Figure 3 and for a reference viscosity $\mu_0 = 10^{18}$ Pa s (solid line) and $10^{19}$ Pa s (dashed line) and the two nucleus radius, 35 km (green curves) and 90 km (magenta curves).

For a fluid with strongly temperature depending viscosity as the case of parent bodies of meteorites, convection takes place under a viscous, stagnant lid. The largest part of the viscosity contrast is encompass by the stagnant lid; convection develops preferentially in the least viscous part of the fluid layer. The thickness of the lid and temperature at its base are calculated according to equation 13 and shown in Figure 14.