© ESO 2019

1 Introduction

Radiometric ages for chondritic meteorites and their components provide information on the accretion timescale of chondrite parent bodies, and on the cooling history within certain regions of these bodies. However, to apply this age information to constrain the internal structure, and the accretion and cooling history of the chondrite parent bodies, the empirical cooling paths obtained by dating chondrites must be combined with theoretical models of the thermal evolution of planetesimals. Important parameters in such thermal models include the initial abundances of heat-producing, short-lived radio nuclides (²⁶Al and ⁶⁰Fe), which depend on the accretion time of the parent body, the terminal size of the parent body, and the chemical composition and physical properties of the chondritic material. Some of these parameters, like material properties, can be determined from laboratory investigations of meteorites, and others, like accretion time and radius, have to be found by comparing evolution models with empirical cooling histories of the meteorites. For L chondrites and their parent body this has been done in various degrees of approximation by Miyamoto et al. (1981), Bennett & McSween (1995, 1996), Benoit et al. (2002), Bouvier et al. (2007), Sprung et al. (2011), Mare et al. (2014), and Blackburn et al. (2017). Here we use the more detailed method for modelling the structure and thermal evolution of planetesimals of ~100 km size of Henke et al. (2012a, henceforth called paper I) and its extensions (Henke et al. 2012b, 2013; Gail et al. 2015; Gail & Trieloff 2018) to reconstruct the properties of the L chondrite parent body and to evaluate the most important parameters that determined its thermal evolution.

For comparison with theoretical models we searched the literature on L chondrites for published data on closure times for the diffusion of decay products of radioactive nuclei out of the carrier phases. The amount of useful data is small because most of the L chondrites were heavily shocked by a catastrophic impact event about 470 Ma ago, which likely disrupted the parent body and is thought to be responsible for the high concentration of fossil meteorites in mid-Ordovician marine limestone in southern Sweden (Heymann 1967; Heck et al. 2004, 2008; Korochantseva et al. 2007). For most of the L chondrites the corresponding data are reset by this event and not useful for our purposes. What one needs are meteorites with a low shock grade for which accurate closure temperature data have been determined for at least two different decay systems. One is required for fixing the burial depth of the meteorite in the parent body, a second one (and ideally further ones) is required to fix the properties of the parent body. Data that satisfy this requirement are found for only five meteorites. These are just sufficient to determine the radius and formation time of the L chondrite parent body. For fourteen further meteorites only data of insufficient accuracy or only for one decay system are found.

In view of the meager data set, we construct a somewhat simplified thermal evolution model, which only considers thermal conductivitywith temperature dependent heat conductivity and heat capacity, and sintering of the initially porous material. The only free parameters in our model are the radius and the formation time of the body; all other parameters are set to values known from other sources or to plausible values. Additionally we allow for the presence of a regolith layer at the surface with strongly reduced heat conductivity due to impact-induced micro-porosity. Because no theoretical model exists for predicting the thickness of such a layer, this is treated as a free parameter to be determined by the model fit. We fit such models to the set of high-quality thermochronological data by means of a genetic evolution algorithm (Charbonneau 1995). We compare the best models found with the data for the individual meteorites and show that it is possible to find a thermal evolution model that reproduces the empirical cooling histories of all the meteorites used here.

Different from the case of H chondrites, two different models with different radii but similar formation times are found to reproduce the empirical data reasonably well, which is not unusual for non-linear optimisation problems. For the smaller one the central temperature stays well below the eutectic melting temperature of the Ni,Fe-FeS system, which is in accord with the fact that no primitive achondrite is known to be compositionally related to L chondrites and may be derived from the same parent body. For the model with bigger radius, a small part of the volume shows temperatures above the eutectic melting, but not to such an extent that differentiation can also be expected. This opens the possibility that primitive achondrites from the parent body of the L chondrites might exist, though for some reason none are represented in our meteorite collections.

We also searched the literature on L chondrites for published data on cooling rates derived from Ni concentrations in taenite, from Pu fission tracks, and from W diffusion between silicates and metal by the new approach of Hellmann et al. (2019). Such data are found for fourteen meteorites, most of them metallographic cooling rates, but all of them are of insufficient accuracy to be useful forthe reconstruction of the parent body. They can only be used to check if they are consistent with the constructed model. Despite this limitation this is of interest because it is claimed in Taylor et al. (1987; and emphasised in Scott et al. 2014, for the case of H chondrites) that metallographic cooling rates are in contradiction to the onion-shell model according to which chondrites of petrologic types L3 to L6 originate from a sequence of layers of increasing temperature from the surface to the interior. We find here, as far as can be judged on the basis of the little available data, that within the accuracy limits of thedata no discrepancy between empirical cooling rates and those found for the parent body model exists.

The plan of our paper is as follows. Section 2 presents the available meteoritic data, Sect. 3 describes how the evolution model of the parent body is calculated, and Sect. 4 describes how the models are fitted to the empirical data. Section 5 presents some properties of our preferred model for the parent body of the L chondrites, Sect. 6 compares the model to previous published models, and Sect. 7 gives some final remarks. Appendix A describes details of how the required material properties of chondritic matter are determined.

2 Data for chondrites

2.1 Thermochronological data

Table 1 summarises isotopic age data for L chondrites that are > 4.4 Ga old and could principally be related to parent body cooling during the first 200 Ma of solar system history. To constrain a cooling path of a certain meteorite within its parent body, at least two data points are needed: two ages of isotopic systems with distinct closure temperature. These meteorites are listed in the upper part of Table 1: Bruderheim, Saratov, Elenovka, Homestead, Barwell, Ladder Creek, and Marion(Iowa). Other meteorites where only one isotopic age data point is available are listed in the lower part of Table 1 for completeness, but not used as a modelling constraint.

Table 1 excludes isotopic ages younger than 4.4 Ga for L chondrites, particularly ⁴⁰Ar-³⁹Ar ages clustering around 0.47 Ga (see e.g. Korochantseva et al. 2007, and references therein). These young ages are frequently associated with high shock stages and interpreted as major impact or even impact disruption of the L chondrite parent body(e.g. Heymann 1967). Such a late reset of radiometric clocks primarily affects isotopic systems with a low closure temperature such as U,Th-He or K-Ar chronometers, and may also disturb – though not reset – chronometers with higher closure temperature such as U-Pb. On the other hand, some rocks may have preserved their original isotopic memory by escaping the thermal effects of late impacts, although some shock effects are discernible in the rocks. Examples of such rocks could be L6 chondrites Bruderheim and Marion (Iowa) in Table 1 for which some studies ascribed shock stage S4. An indication of preservation of the earliest isotopic memory is the presence of decay products of short-lived nuclides like ¹⁸²Hf and ¹²⁹I. While Hf-Wages can be used for cooling curve reconstructions, I-Xe ages are hardly suitable, as for the meteorites in the upper part of Table 1 (Bruderheim, Saratov, Elenovka) I-Xe ages were obtained on whole rock samples, where the iodine carrier phases and closure temperature are unknown. However, the fact that I-Xe ages could be obtained at least indicates that possible late disturbances did not erase the isotopic memory of the in situ decay of ¹²⁹I during the first tens of millions of years of parent body evolution. All in all, most of the isotopic age data used for the reconstruction of the early parent body history are Hf-W and U-Pb ages with rather high closure temperature, that is, isotopic systems that are least susceptible to impact-related disturbance.

2.2 Cooling rates

Reconstructed cooling curves obtained by modelling can further be compared with cooling rates listed in Table 2. Cooling rate data could be found for twelve meteorites, but only for Bruderheim, Barwell, and Elenovka are thermochronological data also available. Only these can be used for comparison. The other meteorites with only coolingrate data are listed in Table 2 for completeness, but not used for comparison with models.

3 Thermal evolution model

We aim to calculate the temperature evolution within the parent body of the L chondrites and to compare this with the thermochronological data of the set of meteorites discussed in the preceding section. Because we have only relatively little data at hand to compare with, we construct a simple model that depends only on a small set of parameters, which we attempt to fix by comparison with the observed properties of the L chondrites.

3.1 Basic assumptions

The basic assumptions are the same as in almost all model calculations of the internal evolution of parent bodies of ordinary chondrites (cf. the review of McSween et al. 2003):

• 1.

  The thermal evolution of the body is dominated by the decay of short- and long-lived radioactive nuclei, in particular ²⁶Al and ⁴⁰K, and others. Other possible heat sources, in particular external heating, are assumed to have been not important.

• 2.

  The body is formed within a period that is short compared to the decay time of ²⁶Al from the mixture of dust and chondrules in the solar nebula, such that the details of the short initial formation period are not important with respect to the long-term evolution of the body, which is the main subject of the model calculation. In the case of the H chondrite parent body, such a rapid accretion was found in Henke et al. (2013) and Monnereau et al. (2013) to provide the best fit to the meteoritic record. Because the parent bodies of L and H chondrites likely formed in close-by zones to the solar nebula, growth conditions should be similar for both bodies. Then we can start the model calculation with an initially cold body that has already achieved its final radius (instantaneous formation approximation).

• 3.

  We assume that the parent body of the L chondrites has a spherically symmetric structure such that the temperature and all other physical variables depend only on the radial distance from the centre, r, and on time, t.

• 4.

  The initial composition of the matter is assumed to be homogeneous within the body. Also the heat production rate per unit mass is assumed to be initially homogeneously distributed within the body and continues to be so during the whole evolution (no differentiation).

• 5.

  We do not consider melting since presently there are no indications that the parent body of the L chondrites was molten and possibly differentiated in its interior. No achondrites or at least primitive achondrites are presently known that can definitely be related to L chondrites, and the L7 chondrites showing evidence of melting are interpreted as impact melts (Mittlefehldt & Lindstrom 2001). Recently, however, a primitive achondrite with oxygen isotopic ratios typical for the L group was found (Vaci et al. 2017a), which may indicate that there was incipient melting deep in the interior of the parent body of the L chondrites. It remains to be shown, however, that this meteorite is really related to the same parent body like the other L chondrites.

• 6.

  We also do not consider giant impacts by other planetesimals that may have already disrupted the body during the early period of its evolution. This is because we aim to study whether it is possible to explain the meteoritical record with respect to temperature evolution with a relatively undisturbed evolution of the parent body over the initial period of heating and cooling during the first approximately 100 Ma of its evolution. This is not unlikely since, according to dynamical calculations, 100-km-sized bodies have a high likelihood of surviving the first 100 Ma un-disrupted (Davison et al. 2013).

3.2 Equations

With our assumptions the temporal evolution of T(r, t) is governed by the heat conduction equation $${\varrho }_{\text{b}}{c}_p\frac{\partial T}{\partial t}=\frac{1}{r^2}\frac{\partial }{r}{r}^{2}K\frac{\partial T}{\partial r}+{\varrho }_{\text{b}}h,$$(1)

where ϱ_b is the mass density of the asteroid material, c_p its specific heat capacity per unit mass, T the temperature, K the heat conduction coefficient, and h the rate of heat production per unit mass. The quantities ϱ_b, c_p, K, and h depend on the properties and composition of the material from which the body is built, and in part on temperature and porosity. The dependence of the heat conductivity K on porosity Φ and T and of the heat capacity c_p on T is complex, as found by experimental determinations of K and c_p for chondrites (Yomogida & Matsui 1983; Opeil et al. 2010, 2012) and in the theoretical investigation by Henke et al. (2016) and Gail & Trieloff (2018). The details of how the properties of the material are calculated are given in Appendix A.

Since the material in ordinary chondrites is subject to heating up to temperatures above 1000 K, the composition and physical properties of the mineral mixture change significantly. Besides the temperature dependencies of K and c_p, the most important change is reduction of the initial pore space (sintering) since this results in particular in a considerable increase of the heat conductivity, K, and, to a lesser extent, to an increase of the bulk density, ϱ_b, of the mineral mixture. The bulk density, the intrinsic density ϱ_i of completely compacted material, and the volume fraction of voids, Φ, are related by $${\varrho }_{\text{b}}=\left(1-\Phi \right){\varrho }_{\text{i}}.$$(2)

The variation of the porosity depends on temperature, T, and the applied pressure, p. Theories for calculating the change of Φ were developed for technical purposes in material science where a set of equations is derived for calculating $$\frac{\partial \Phi }{\partial t}=F\left(\Phi ,T,p\right).$$(3)

Here we apply the sinter model of Helle et al. (1985). The application of this model to parent bodies of ordinary chondrites is discussed in Gail et al. (2015); it requires the solving of a differential equation for Φ simultaneously with the heat conduction equation.

The pressure p required for calculating the sintering process is determined by the hydrostatic equilibrium equation $$\begin{array}{l}\frac{\partial p}{\partial r}=-\frac{G{M}_r}{r^2}{\varrho }_{\text{b}}\\ \frac{\partial M_r}{\partial r}=4\pi r^2{\varrho }_{\text{b}}.\end{array}$$

This equation also has to be solved simultaneously with the heat conduction equation and the equation for sintering.

3.3 Initial and boundary conditions

The heat conduction equation is subject to initial and boundary conditions. Within the instantaneous formation approximation the initial condition is $$T(r,0)=T_0(0\le r\le R),$$(6)

where R is the radius of the body and T₀ the initial temperature at t = 0. At the centre of the body the solution has to satisfy the boundary condition $$\frac{\partial T(0,t)}{\partial r}=0(t>0)$$(7)

to exclude the presence of a point source for heat at the centre. At r = R we would have to determine the surface temperature from the equilibrium between heating and cooling processes at the surface (see e.g., Henke et al. 2012a). In practice, insufficient information is available for this and we only apply $$T(R,t)=T_{\text{s}}(t)(t>0)$$(8)

with given surface temperature T_s(t) as outer boundary condition.

In view of the uncertain initial conditions at the time of planetesimal formation and the uncertainties with respect to the environmental conditions experienced by the body during its later evolution, we can only prescribe a plausible value T_s for the surface temperature and a spatially constant initial temperature T₀. A choice for T₀ and T_s only needs to obey two conditions. First, the initial temperature and, at least as long as the accretion disc still exists, also the boundary temperature have to be higher than ~150 K because otherwise the parent body of the L chondrites would have been formed with abundant water ice, which clearly was not the case. Second, the boundary temperature has to be well below the annealing temperature of ~ 390 K for Pu fission tracks, because such fission tracks are observed for some L chondrites (Pellas & Störzer 1981). We chose in our model calculations T₀ = T_s = 200 K.

The initial porosity Φ₀ of the chondritic material can be estimated from the properties of the unequilibrated meteorites. This is discussed in Appendix A.1. The estimated value of Φ₀ = 0.25 for the parent body of the L chondrites is used as an initial value for the solution of Eq. (3).

3.4 Solution of the model equations

The simplest possibility to obtain a solution of the model equations would be to restrict ourselves to the kind of model considered first by Miyamoto et al. (1981). They considered only the heat conduction problem and use for the material properties ϱ, c_p, and K constant values as derived from measured properties of meteorites at room temperature. For constant coefficients, Eq. (1) can be solved analytically, but in reality the corresponding quantities are by no means constant, which significantly changes the thermal structure. Hence, for obtaining physically reliable models we solve Eqs. (1)–(5) numerically, which allows us to consider the dependency of c_p and K on temperatureand other physical properties of the material. The heat conduction equation and the sintering equation are solved by an implicit finite difference method. The non-linear character of the equations is accounted for by iteration to self-consistency by a fixed point iteration.

4 Model fit

4.1 Parameters

The thermalevolution model depends on a number of parameters. Some of these parameters can be derived from the properties and composition of meteorites, while some parameters are not known a priori since they depend on the unknown formation history and formation site of the parent body of the meteorites. The most important unknown parameters are the radius of the body, R, its birth time, t_b, and the thickness, d, of an outer impact-modified layer with reduced heat conductivity. Also the initial temperature, T₀, at formation time and the radiative equilibrium temperature, T_s, of the surface at the (unknown) distance from the sun are not known. These five unknown parameters can in principle be fixed by determining a set of parameter values for which the deviation between the empirical information on the cooling history of the meteorites and a thermal evolution model takes a minimum, provided a sufficient number of empirical data is available. This is how we proceeded for the H chondrites (Henke et al. 2013). We think it better, however, not to overburden the simple thermal evolution model described above by considering details that have merely a limited influence on the model. This holds in particular for the initial and boundary temperatures, T₀ and T_s, which are therefore not included in the optimisation but instead are fixed to typical values. Therefore, we restrict the set of model parameters to the birth time, t_b, the radius of the parent body, R, and the thickness of the insulating layer, d. This set of parameters is considered as the best information on the likely properties of the parent body of the L chondrites that can presently be obtained.

4.2 Comparison with meteoritic data

The data set for L chondrites to be compared with evolution models for their parent body is summarised in Sect. 2. We can only use the seven L chondrites Bruderheim, Saratov, Elenovka, Homestead, Barwell, Ladder Creek, and Marion (Iowa), for which useful data on closure ages for at least two different thermochronometers with different closure temperatures are available. In the fitting procedure between models and meteoritic data we omit, however, the data for Saratov and Homestead because their Hf-W ages may not have been completely reset by mild metamorphism and they may partly keep the memory of chondrule formation processes (J. Hellmann, priv. comm.). The final data set used for model fitting is shown in Table 3; data are much more limited than for H chondrites, but should in principle suffice to perform a meaningful fit.

In the following, m enumerates the meteorites used. For each of these we know for totally N(m) different thermochronologic systems closure temperatures T_m,n and closure times t_m,n and associated errors σ_T,m,n of the closure temperature, where n =1, …, N(m) enumerates these systems. The burial depths of the meteorites, d_m, are determined by first determining for a given model and a given meteorite that depth for which $${\chi }_m^2={\displaystyle \sum _{n=1}^{N(m)}\frac{{\left(T(t_{m,n})-T_{m,n}\right)}^2}{{\sigma }_{T,m,n}^2}}$$(9)

takes a minimum. Here T(t_m,n) is the modeltemperature at the burial depth of the metorite. The minimum is determined by calculating for a grid of burial depths the values of ${\chi }_m^2$ and determining the minimum from fitting a parabola through the lowest value and that of the two adjacent points of the grid. The location at which the minimum is taken corresponds to the burial depth of the meteorite. The quality function then is $${\chi }^2={\displaystyle \sum _m}{\chi }_m^2+P$$(10)

with the values of ${\chi }_m^2$ as found in the determination of the burial depths.

We can add a penalty function $\mathcal{P}$ to χ² to enforce the model to obey additional restrictions, for example to guarantee that the central temperature stays below the eutecting melting temperature of the Ni,Fe-FeS system. In this latter case we could chose $\mathcal{P}$ for example as follows $$\mathcal{P}=\max \left(T_{\text{c}}-T_{\text{m}},0\right),$$(11)

where T_c is the maximum central temperatureof the model and T_m the melting temperature.

To obtain a best fit between data and model results, and thus to minimise the quality function, χ², a fitting routine is needed for which we adopt the genetic algorithm PIKAIA of Charbonneau (1995), which, like any genetic algorithm, is in principle able to find the global minimum even for highly complicated problems. This has also been used in our preceding papers and found to be very useful for this purpose.

In our calculation we vary three model parameters to optimise a model: the radius, R, the birth time of the body, t_b, and the thickness of an impact-damaged surface layer, d. The parameters are varied each within a range of physically meaningful values, which are given in Table 4. This limitation of the range of variation serves to prevent the algorithm from trying to explore unphysical regions of the parameter space.

To achieve optimisation, we ran the genetic code for at least 40 generations with 50 individuals. We always performed at least ten optimisations with different seeds for the random number generator in order to avoid the algorithm getting stuck for a large number of steps in a local minimum and the calculation being stopped too early before a mutation pushes the algorithm to continue the search for a better result in another region of the parameter space. We used a resolution of the parameter values of four digits (see Charbonneau 1995, for the meaning of this).

4.3 Optimised model

Ten optimisations with different seeds of the random number generator were run for the range of reasonable parameter values, as defined in Table 4, with 40 generations and 50 individuals each. The optimisation resulting in the lowest value of χ² (usually there are more than one with an almost equal value of χ²) was re-run with 100 generations. The optimum model was found in step 1884 of 5000 and no improvement was found in the following steps. This suggests that we really have found the optimum model corresponding to the lowest value of χ². We found that within the range of admitted parameter values some additional well-defined local minima of the quality function χ² exist where χ² takes a localminimum value comparable to that of the optimum model. The set of parameter values corresponding to the optimum model and to two of the alternative local minima are shown in Table 5. These models are labelled LM1, LM2, and LM3, with LM3 being the optimum model. If we apply the penalty function Eq. (11), only modelsLM1 and LM2 are obtained. Models LM1 and LM2 have similar low values of the quality function χ² to the optimal model LM3 and therefore fit the empirical data nearly equally well. Table 6 presents the resulting burial depths of the meteorites and the peak temperatures achieved during the evolution at the burial depths.

Table 5 also shows the normalised fit quality $${\chi }_{\text{n}}^2={\chi }^2/(D-P),$$(12)

which is a better indicator for a proper characterisation of the minimum. Here D is the number of data points fitted (in this case D = 10) and P the number of parameters (which is P = 8, the three varied model parameters plus the five burial depths that have to be determined). For a resulting value of ${\chi }_{\text{n}}^2<1,$ the model can be considered as agreeing with the data sufficiently well such that we can assume that it reproduces the meteorite cooling curves as documented by the thermochronological data, whereas a value of ${\chi }_{\text{n}}^2\gtrsim 1$ would mean that there is a significant misfit between model and data. The value ${\chi }_{\text{n}}^2~0.0109$ for the optimum model LM3 given in Table 5 satisfies ${\chi }_{\text{n}}^2\ll 1,$ which means that for our evolution model a set of parameter values has been found for which the corresponding model reproduces the data. Models LM1 and LM2 also satisfy this condition.

To get an impression of how well it is possible to fit the set of available thermochronological data of the L chondrites using our thermal evolution model, we show in Fig. 1a the χ² (R, t_b) surface and the contour lines of constant χ² in some range around the optimum set of parameter values for varying R and t_b and for a fixed value of d equal to its optimum value. We easily recognise the well-defined minimum in the χ² surface around the optimum value of R of 160 km and a formation time of 1.84 Ma. Figure 1b shows a more close-up view of the projections of a number of contour lines onto the R–t_b plane around the position of the optimum R–t_b values at fixed d. An inspection shows that the formation time and the radius of the body are rather well determined if d is held fixed at its optimum. A deviation of χ² by a factor of two from its optimum corresponds to a deviation of the formation time from its optimum value by about ± 0.01 Ma, and a deviation of the radius from its optimum value by about ± 10 km. This indicates the accuracy with which the formation time and the radius of the parent body can be determined in our case.

Fig. 1 Variation of χ2 with radius R of the body and birth time tb for fixed thickness d of the outer impact-damaged layer. Left part: shape of the χ2 surface. Alsoshown are contour lines of constant values of χ2 on the χ2 surface and their projection onto the basis plane. Right part: contour lines of χ2 in the R-tb plane.

4.4 Alternative models

Additionally one finds minima of χ² at radius >200 km and at R ≲ 140 km, but with a somewhat bigger value of χ² at minimum. The bigger body would be of the size of Vesta and as such is unlikely to have been completely disrupted during the lifetime of the solar system. Since Vesta is not the source of the L chondrites and no other body (except for Ceres) of this size is found in theAsteroid belt, the bigger model is rejected.

For the two other local minima in the χ² surface at R ≲ 140 km, we performed an optimisation with a restricted search volume in parameter space. The resulting optimum parameter values and the corresponding values of χ² are given in Table 5 as models LM1 and LM2. In both cases ${\chi }_{\text{n}}^2$ is quite small, though not as small as for model LM3, but both models also fit the empirical data quite well. For models LM1 and LM2 there is an important difference in the model properties compared to LM3: the central temperature definitely stays below the temperaturefor eutectic melting of the Ni,Fe-FeS system of 1223 K (Mare et al. 2014), while for model LM3 it definitely exceeds this melting temperature.

If model LM3 correctly describes the properties and early evolution of the L chondrite parent body, then there could exist primitive achondrites with signatures of incipient eutectic melting of Ni,Fe-FeS that show a clear relation to L chondrites according to their composition. If, on the other hand, no such primitive achondrites related to L chondrites are found, models LM1 or LM2 probably describe the properties and early evolution of L chondrites. In fact, despite L chondrites being the most abundant known ordinary chondrites, with presently about 7100 known specimens, no primitive achondrite clearly related to the class of L chondrites has ever been found. This suggests models of this type to be appropriate. Because of its slightly better fit quality, we prefer model LM2 in the following.

Recently, however, a primitive achondrite with oxygen isotopic ratios typical for the L group was found (Vaci et al. 2017a), which may indicate that after all there was incipient melting deep in the interior of the parent body of the L chondrites. This would point to model LM3 as the correct one, but it remains to be shown that this meteorite is really related to the same parent body as the other L chondrites.

5 Model results

5.1 Fit of meteoritic data

In order to demonstrate how the models fit to the meteoritic data, Fig. 2 shows the variation of temperature at a number of selected depths during the course of the evolution for models LM2 and LM3. Also shown are the thermochronological data for the meteorites for which at least data of acceptable accuracy for a minimum of two thermochronological systems are available, including the meteorites that are not used in the optimisation. The thermal evolution of the five meteorites from Table 6 is shown as black lines. The temperature curves for each of these meteorites pass right through the centre of the error boxes for almost all of their thermochronological data (see Fig. 2). The fit therefore is quite good, which may in part result from the circumstance that we only use those data for the optimisation that can be considered as reliable and highly accurate. If one also includes data of lower accuracy and unclear reliability from Table 1, the fit quality would become much worse.

The quality of the fits to the data is practically indistinguishable for the two kinds of models with either no melting in the core region or with incipient melting in the core. The presently available data do not allow us to discriminate between the smaller or bigger sized body on the basis of the thermochronological data. It is only the lack of primitive achondrites related to L chondrites that favours model LM2.

We now consider the individual fits between empirical data for the cooling history of a meteorite and the thermal evolution model of the parent body for the five meteorites used for the optimisation. In the upper part of the pictures, Fig. 3 shows the temperature evolution of the favoured model LM2 at the burial depth of the meteorites. In the lower part of the pictures, the cooling rate after culmination of the temperature is shown. In the upper part we also show (as red dots) the data points of the closure time and closure temperature of the radioisotopic clocks used for the optimisation, with error bars for closure time and closure temperature. In the lower parts of the picture we show data for cooling rates for the meteorites; besides radioisotopic clock data, this kind of data are also available for three meteorites. The figure confirms again that it is possible to obtain a very good fit between the highly accurate Hf-W and Pb-Pb data and an evolution model for the parent body.

If such data exist, cooling rates are not used for the model fit because of their generally low accuracy. For the meteorites shown in Fig. 3, metallographic cooling rates are published for Elenovka and Bruderheim, and cooling rates determined by the new method of Hellmann et al. (2019) for the Hf-W system are available for Barwell and Bruderheim. In the case of the metallographic cooling rates we derive a closure time for a closure temperature of 773 K for Ni diffusion in taenite from the temperature run at the known burial depth, and we plot the corresponding data point in the lower part of the pictures. The accuracy of the cooling rates derived from Ni diffusion in taenite is estimated by Wood (1967) as being about a factor of 2.5. For the Hf-W system, the time to which the cooling rate refers is given directly in Hellmann et al. (2019).

In total,we have four pieces of data on the cooling rates of three meteorites that can be compared to our parent body model. For these empirical cooling rates, the cooling rate curve at the burial depth of the parent body model agrees within the error limits with the observed cooling rates. This suggests that the cooling rates of the model for the parent body also agree reasonably well with the empirical cooling rates. In summary, we conclude that it is possible to find a model that reproduces the observed cooling history of all L chondrites for which presently accurate thermochronological data are available.

Fig. 2 Temperature evolution at different depths below the surface (grey lines) and evolution of the central temperature (dashed grey line). The symbols with error bars correspond to the closure temperatures and closure times of L chondrites. The black lines correspond to the temperature evolution at the burial depths of the five meteorites from Table 6 for model LM2 (a) and model LM3 (b). The data of Homestead and Saratov are not used for the optimisation. Grey dots show chondrule ages for L chondrites from Pape et al. (2019).

5.2 Age of chondrules

As chondrules are formed as individual millimeter sized objects in the solar nebula before accretion of the parent asteroid, the parent body formation time resulting from our model serves as a lower limit for the age of chondrule formation. Figure 2 shows ²⁶Al/²⁶Mg ages for a number of individual chondrules from types L3.00 to L3.15 from Pape et al. (2019). These are probably the most accurate chondrule ages presently known. They range from an average age of ${1.80}_{-0.10}^{+0.11}$ Ma for the oldest chondrules up to ~ 3 Ma after the formation of CAIs. They are plotted in Fig. 2 with error bars corresponding to the possible range of initial temperatures for ordinary chondrite parent bodies. The age of the oldest chondrules is interpreted by Pape et al. (2019) “to record the major and relatively punctuated onset of chondrule formation”. If this interpretation is correct, our result on the birth time of the parent body of the L chondrites suggests that the formation of the body followed almost immediately the onset of chondrule formation, within ~ 0.03 Ma if model LM3 is applicable or within ~0.08 Ma if model LM2 is applicable. This means that it is within the errors bounds for the time of onset of chondrule formation.

The ages of the chondrules in Fig. 2 are for chondrules from the most primitive un-equilibrated meteorites with the lowest petrologic types L3.00 to L3.15 of ordinary chondrites. They resided in layers within a few kilometres below the surface of the body, which most likely acquire some additional later-formed material even after the main formation period of the body has terminated. The spread in chondrule ages over a period of more than 1 Ma, therefore, is not necessarily in conflict with the instantaneous formation period hypothesis used in our model calculations, but hints either at an extended period of chondrule formation over 1.2 Ma duration (Pape et al. 2019) or subsequent disturbance of the Al-Mg system.

Fig. 3 Temperature evolution at the burial depths of the meteorites and the corresponding cooling rates (in K per Ma). The dots with error bars show the empirical data. Red dots show the data used for the optimisation, blue dots show the cooling rates.

5.3 Properties of the parent body

In the following we consider some consequences of the thermal evolution of the body. The variation of temperature at fixed radius during the course of the thermal evolution of the body is shown for some selected radii in Fig. 2. There is a rapid increase in temperature within ~3 Ma immediately after formation of the body, a period of high temperature lasting several to tens of Mega years, and a final cooling period lasting tens to hundreds of Mega years.

The evolution of the temperature at the centre of the body is shown by the dashed line. The central temperature has a maximum at about 40 Ma for model LM2 and 80 Ma for model LM3. The corresponding flat hump in Fig. 2b is due to the contribution of ⁴⁰K to the heating of the body. Without this the temperature inside the body would decrease significantly earlier, which would also modify the thermal evolution at the lower depths from which the meteorites come (cf. Fig. 8).

An important aspect of the thermal history is the highest temperature, T_max, found at a given radius since this determines the degree of metamorphism from the complex initial mixture of dust and chondrules seen in petrologic type 3 chondrites to the compacted and crystallised rocky material seen in type 6 chondrites, or, maybe, even to partially molten material of primitive achondrites. The radial variation of this maximum temperaturein the model is shown in Fig. 4a for model LM2 and in Fig. 5a for model LM3.

There are three characteristic temperatures where the structure of the material and its heat conductivity significantly change. Thefirst is where the initially porous matrix material is compacted by surface diffusion (Gail et al. 2015), the second is where the ensemble of chondrules is compacted by deformation creep (Gail et al. 2015), and the third where the first components of the complex mixture of minerals and metal of the chondritic material start melting.

The compaction of matrix and chondrules extends over a temperature range that is indicated in Figs. 4a and 5a by two red strips. Each of the two processes result in a significant increase of the strength of the material. At about 560 ± 20 K the probably rather loosely bound initial dust and chondrule mixture bake together into a sandstone-like material. Around 1060 ± 30 K creep processes of the minerals in chondrules become thermally activated such that under the action of pressure inside the body the remaining pore space between chondrules disappears. At the same time, at such a high temperature crystal growth of the mineral components and diffusional re-distribution of iron into coarser inter-chondrule-space-filling iron grains becomes active, which altogether act to form a rocky material of significant strength.

The overall run of T_max(r) with distance r from the centre as seen in Figs. 4a and 5a is rather flat across most parts of the body and steeply drops to the low surface temperature within a thin surface layer. This nearly isothermal interior of the body and the steep gradient near the surface results from efficient heat conductivity in the compacted interior of the body and low heat conductivity in the porous and impact-damaged structure of the surface material.

The variation of K with radius in our model of the L chondrite parent body is shown in Figs. 4b and 5b. Since the heat conductivity depends on temperature, and thus on time, we show the heat conductivity at 30 and 65 Ma after the formation of the body for models LM2 and LM3, respectively, which is representative for the period where the temperature in the central region is highest. The transition from the porous outer material to the completely compacted material in the interior occurs at a depth between 11.0 and 17.5 km below the surface in model LM2 and between 8.9 and 11.5 km below the surface in model LM3. The second significant change in K occurs at a depth of 2.9 km below the surface in model LM2 and at 3.5 km below the surface in model LM3. This is because in our model we allowed for a surface layer where the material is heavily damaged by impacts and where numerous micro-cracks in the mineral component of the material significantly reduce the heat conductivity to a low value as it is also measured for part of the chondrites (see Fig. A.2b).

According to the model optimisation, such a layer is required to obtain a good fit to the empirical data. Since we presently cannot predict d from a theory of the impact-induced modifications of the properties of the surface material on an asteroid, it is hard to decide whether our optimum value of d obtained from the fitting process is physically justified or simply allows us to compensate for other deficits of the model. However, at least one point seems clear: a reliable fit with only a very thin or even absent insulating surface layer cannot be obtained for the L chondrite parent body. The situation would probably improve if data both for L5 and L6 chondrites and also for undisturbed L4 chondrites were available. They are more sensitive to the temperature at shallower surface layers and would allow us to better to constrain the properties of the outer part of the body.

Fig. 4 Radial variation of maximum metamorphic temperature and heat conductivity. (a) Maximum of temperature at radius r during the undisturbed thermal evolution of the body (black line) for model LM2. The two broad red strips indicate thetemperature range over which the chondrule assembly (upper strip) and the initially fine-grained matrix material (lower strip) are compacted. The hatched area is the temperature range between the upper metamorphic temperatures of L4 and of L5 chondrites (973 to 1023 K, respectively, according to McSween et al. 1988) where the inter-granular voids between chondrules disappear. The upper blue line shows the eutectic melting temperature of 1223 K for the Ni,Fe-FeS system at the Ni content of L chondrites. (b) Radial variation of the heat conductivity K at 30 Ma after formation of the body.

Fig. 5 As Fig. 4, but for model LM3.

5.4 Onion-shell model

It is generally assumed that the metamorphism of the chondritic material from the highly pristine material seen in L3.0 meteorites to the highly equilibrated material seen in L6 meteorites is essentially determined by the peak temperature experienced by the material. According to the onion-shell model hypothesis, this thermal metamorphism proceeded in different depth regions of the parent body that experienced different maximum temperatures (as shown in Fig. 4) during the course of the thermal evolution of the body. Meteorites of type L6 correspond in this model to material from the core region, type L3 to material from a layer just below the surface, and types L4 and L5 from layers of increasing depth and temperature below the surface. The basic assumption then is that the body experienced no disruptive collision within at least the first 100 to 200 Ma of its internal thermal evolution, such that all layers cooled below the closure temperatures of the thermochronological systems, and at least some meteorites can be found where their thermochronological clocks remained essentially undisturbed if the body was shattered into a rubble pile.

Since the seminal paper of Yomogida & Matsui (1983), the boundaries between regions where the material shows characteristics corresponding to the different petrologic types are calculated from maximum metamorphic temperatures of a given petrologic type. We take data for these temperatures for type L3 and L4 (873 and 973 K, respectively) from McSween et al. (1988) and for L6 (1173 K)from Slater-Reynolds & McSween (2005). The maximum peak temperatures are generally highly uncertain (about ± 50 K, e.g. Kessel et al. 2007), in particular for petrologic types 4 and 5 for which they are only estimates. Especially the temperature of 1023 K for the border between L5 and L6 as given by McSween et al. (1988) seems problematic, since with this value for the border temperature the width of the temperature range for type L5 between the lower temperature border separating L5 from L4 and the higher temperature border separating L5 from L6 is of the same order of magnitude as the error of the two border temperatures. Therefore, we proceed in this case as follows. According to the classification of van Schmus & Wood (1967), one major difference between petrologic types 5 and 6 is that for type 6, chondrules are poorly defined, while for type 5, chondrules are readily delineated. We then identify the border between types 5 and 6 by the disappearance of chondrules by sintering as easy-to-recognise individuals and take (somewhat arbitrarily) a residual pore space of the granular material of 5% as limit, which translates according to our model into a temperature at the border between L5 and L6 of 1 100 ± 25 K (the error limits corresponds to almost none to almost complete sintering).

For our preferred model LM2 the depths at transitions between the different petrologic types with these values of temperature are given in Table 7. A schematic sketch of the distribution of petrologic types in the body is shown in Fig. 6. The volume fractions of types L3–L6 are considered in Sect. 6.4.

Fig. 6 Schematic sketch of the onion-shell structure of the L chondrite parent body.

5.5 Melting and differentiation

The first melting process in chondritic material would be the eutectic melting in the Ni,Fe-FeS system. It is usually assumed that central temperatures in the ordinary chondrite (H, L, LL, EH, EL) parent bodies always stayed below the temperature required for the onset of melting. The reason for this assumption is the lack of ordinary chondrites with indications of melting or at least partial melting. Recently, primitive achondrites have been detected whose oxygen isotopic ratios plot within the corresponding fields of L and LL chondrites (Vaci et al. 2017a,b, 2018). In particular NWA 11 042 was interpreted as being derived from the L chondrite parent body. Its texture is achondritic and its elemental compositionresembles that of L chondrites, but it is depleted of siderophile elements compared to L chondrites and the Ni,Fe-metal content is very low. This is interpreted as meaning that the material of NWA 11 042 is derived from L chondrite material by melting and differentiation. Vaci et al. (2017a) speculated that this fits within the frame of the thermal evolution of the parent body according to the “onion-shell” model, but requires higher central temperatures and a bigger parent body radius than existing models. On the other hand, Wu & Hsu (2017) alternatively propose that NWA 11 042 derives from a massive impact-produced magma chamber on the L chondrite parent body.

Here wecheck whether it is possible that the core region in our model of the parent body of L chondrites may be the source of primitive achondrites. This requires that the central temperature well exceeds the melting temperatureof the Fe,Ni-FeS eutectic, but does not result in extended silicate melting, since only partial melting of the silicates is observed. The eutectic melting temperature of Ni,Fe-FeS at contact between Fe,Ni-metal and FeS grains was studied in Tomkins (2009) and Mare et al. (2014). According to them the addition of Ni to the Fe-FeS system lowers the eutectic melting temperature from 1 261 K to 1 223 K. In Figs. 4 and 5 we show the maximum temperature achieved during the thermal evolution of the body for each radius r for models LM2 and LM3, respectively. In our optimum model LM3 the melting temperature is indeed exceeded in a wide inner core region extending from the centre to about one half of the radius¹. Hence, the parent body of the L chondrites could have experienced incipient eutectic melting of FeS in its core region. This region would extend from the centre up to r = 82.0 km, which means that an inner core region of 13% of the total volume would contain material that would be classified as a primitive achondrite.

The absence of primitive achondrites related to L chondrites contradicts such a model. Also the empirically determined upper metamorphic temperature of L6 chondrites of about 1173 ± 30 K given by Slater-Reynolds & McSween (2005) is in conflict with such a model, because the parent body of the L chondrites was most likely completely disrupted, in which case one should have meteoritic samples also from the deepest layers with temperature ≳ 1223 K. The maximum metamorphic temperatures, however, fit reasonably well to the maximum temperature at the centre of 1186 K for our model LM2.

6 Comparison with previous model calculations

6.1 Method of modelling

The numberof published model calculations for the thermal evolution of the L chondrites parent body (Miyamoto et al. 1981; Bennett & McSween 1995, 1996; Benoit et al. 2002; Bouvier et al. 2007; Sprung et al. 2011; Mare et al. 2014; Blackburn et al. 2017) is small, probably because of the little available data on L chondrites suited for a comparison of models with laboratory data. The basic characteristics of the published models are summarised in Table 8, except for the model of Benoit et al. (2002) for which no details are published. Most of these models assume a short duration of the formation period of the body compared to the characteristic heating time by ²⁶Al such that the instantaneous formation hypothesis can be applied. Modern theories of planetesimal formation support this assumption (e.g. Johansen et al. 2015). Only the model calculation of Bouvier et al. (2007) includes growth of the L chondrite parent bodyover an extended period, which results in a somewhat different evolution from that in the other models. In the following weconsider these models in detail.

The model of Miyamoto et al. (1981; and the preliminary communication by Miyamoto et al. 1980) was the first complete model calculation for L chondrites parent bodies, based on internal heating by the short-lived radio nuclide ²⁶Al (it also discusses the evolution of H chondrites). It follows essentially the method proposed by Minster & Allegre (1979) to calculate the thermal history of the H chondrite parent body.

The model calculation is based on the analytic solution of the heat conduction Eq. (1), as given in Carslaw & Jaeger (1959), for the case of spatially and temporal constant values of heat conductivity K, specific heat capacity c_p, density ϱ, and initial temperature T₀, and a spatially constant and exponentially decaying heat source. This analytic model also forms the basis of almost all subsequent models for the thermal evolution of the L chondrite parent body considered so far (Bennett & McSween 1995, 1996; Sprung et al. 2011; Blackburn et al. 2017). The different models merely differ with respect to the particular choice of the constant input parameters.

In the model calculations of Benoit et al. (2002), Mare et al. (2014), and the present work the heat conduction equation is solved numerically. This relaxes all restrictions with respect to the variation of the coefficients and of the initial and boundary conditions and allows us to couple the evolution of temperature with other processes like melting (in Mare et al. 2014) and sintering of the initially porous material (this work).

6.2 Material properties

The requirement that spatially and temporally constant values for the temperature- and porosity-dependent parameters that are employed in the heat conduction equation have to be used in the analytic solution makes it necessary to choose a fixed value for these quantities. These quantities are then thought to be representative for all zones of the body and during the complete heating and cooling period. This casts some doubts on the reliability of the models obtained by using the analytic solution because in the interior of the body the initially porous material will be compacted by sintering at temperatures of 700 to > 900 K (depending on the chondrule diameters), while in the cold outer layers the material remains loosely compacted. Therefore, the heat conductivity will be quite different in the core and in the outer layers of the body. For this reason Bennett & McSween (1996) calculated two different models with constant values of the material properties throughout the body as limit cases, one for a strongly compacted material (tagged “comp.” in Table 8) and another one for a porous material (tagged “por.” in Table 8), assuming that a solution accounting for the variability of material properties would be found somewhere in between. For H chondrites, two-zone models with a conductive core and an insulating mantle to handle this problem have been studied by Akridge et al. (1998) and Harrison & Grimm (2010). For L chondrites such a type of study has only been conducted by Benoit et al. (2002), but no details are given for that model.

The density is chosen in all models according to observations on L chondrites. Yomogida & Matsui (1983) give the average chemical composition of L chondrites, which is also adopted in our model calculation. From this they calculate an intrinsic density, ϱ_i, of 3610 kg m⁻³. This value is used in Miyamoto et al. (1981), Bennett & McSween (1995), the compacted model of Bennett & McSween (1996), and in our model. For the bulk density, ϱ_b, the value of 3200 kg m⁻³ used in Miyamoto et al. (1981) and Sprung et al. (2011) is the average of the H chondrites from Antartica listed in Yomogida & Matsui (1983). This corresponds to an average porosity of about 11%, a value that deviates from values given in Britt & Consolmagno (2003) for a much bigger sample, but since the solution depends only weakly on density its assumed value is not so important for the results. The bulk densities and intrinsic densities in Bennett & McSween (1996) for their two models correspond to the data of the two meteorites Y750977 (porous material) and Arapahoe (compacted material) considered as representatives for loosely and strongly compacted material. The bulk densityused in Mare et al. (2014) is taken from the compilation in Wilkison et al. (2003), no value for ϱ_i or porosity is given. In Blackburn et al. (2017) values for the average bulk density and the intrinsic density are taken from Britt & Consolmagno (2003). The density used in Benoit et al. (2002) is not specified.

A significant variation with porosity is observed for the heat conductivity K, which has a strong impact on the characteristic cooling timescale. Measured values for L chondrites show a wide and unsystematic variation with porosity (between ~0.4 and ~3 W m⁻¹K⁻¹ at 300 K, see Fig. 7, which is about a factor of eight between the lowest and highest value) and quite different types of temperature variation (Yomogida & Matsui 1983; Opeil et al. 2010, 2012). The origin of these differences could be explained in Henke et al. (2016) and Gail & Trieloff (2018). The published model calculations for thermal evolution assume quite different values for the heat conductivity K, ranging between the lower end of observed values for L chondrites in the models Mi81, Be96(por), and Sp11, up to the value of the fully compacted material in Bl18 (4.17 W m⁻¹K⁻¹, see Henke et al. 2016), and a more intermediate value in Be95 and Be96(comp). Using the low values for the heat conductivity means that more stress is laid on the heat transport in the surface-near layers, which is responsible for the global evolution of the heat contentof the body, while high values better account for the temperature structure in the interior of the body where the material is compacted by sintering.

An equation for the global heat content of the body is obtained by integrating the heat conduction Eq. (1) over the whole volume of the body, $$\frac{\text{d}U}{\text{d}t}=-4\pi R^{2}K{\frac{\partial T}{\partial r}|}_{r=R}+M\cdot h,$$(13)

where M is the total mass of the body, U its total heat content, and h the specific heating rate per unit mass. This shows that the long-term thermal evolution is crucially determined by the heat conductivity of the surface layers and that there is some justification in using a representative value for this in the model calculations. A value corresponding to the conductivity of the compacted material as used in model Bl18 in any case overestimates the rate of heat loss. The lower values chosen in the models Mi81, Be95, Be96, and Sp11 seem more appropriate, but the precise reasons that led the authors to their preferred choice are not communicated in the papers. A good deal of the significant discrepancies between the radii found for the parent body of the L chondrites in the different models rests on the different choices for K. In particular the big radius found in model Bl18 is required to compensate for the high rate of heat loss.

Much less disagreement exists between the different choices for a representative value for the heat capacity in models Mi81, Sp11, Be95, Be96, and Bl18 (see Table 8). Values given in italics are not explicitly stated in the papers but calculated from the values given for heat diffusivity, κ = K∕ϱ_bc_p, and heat conductivity K. Chosen values are in the middle range between a lower value of ~200 J kg⁻¹K⁻¹ at a temperature of 200 K at the surface and the highest value of ~ 1200 J kg⁻¹K⁻¹ at a peak temperature of ~1200 K achieved at the centre. Also, in the case of heat capacities, the reasons which led the different authors to their preferred choice are not communicated in the papers. A variation of the cooling timescale resulting from the temperature variation of the heat capacity is moderate and, thus, much less influential for the model than the variation resulting from variations in the heat conductivity.

Most of the models considered assume that the heating of the body is due to the decay of radioactive nuclei and that the heating is dominated by the decayof ²⁶Al, such that itsuffices to consider only this contribution. A possible contribution of ⁶⁰Fe to the heat budget of the body was considered in Bouvier et al. (2007) and Mare et al. (2014). Its importance depends on the assumed initial abundance of ⁶⁰Fe, which is disputed, and on the half-life of ⁶⁰Fe, which was corrected in Rugel et al. (2009). Older determinations of the ⁶⁰Fe/⁵⁶Fe ratio found a value of 9.2 × 10⁻⁷ (e.g. Mostefaoui et al. 2005), which was used in the model of Bouvier et al. (2007, who also used an outdatedvalue for the half-life), while more recent determinations gave a much lower value of 1.15× 10⁻⁸ (Tang & Dauphas 2015), which is used in Mare et al. (2014) and in our calculation. With the low value of the ⁶⁰Fe/⁵⁶Fe ratio the contribution of ⁶⁰Fe to the heating is not important. This is in accord with our finding for the H chondrite parent body that if the ⁶⁰Fe abundance is considered as a free parameter and determined from a best fit of models to meteoritic data, then it is found that the observed thermal evolution is best reproduced without a contribution of ⁶⁰Fe to heating (Henke et al. 2013).

There is, however, some contribution to heating of the body by long-lived radioactives, in particular by ⁴⁰K. For asteroid-sized bodies the main effect of these heat sources is a certain delay in the onset of cooling in the central region of the body. This manifests in the slight increase of the central temperature during the first 100 Ma recognisable in Fig. 2. These heatsources are neglected in all previous model calculations because they are not important for bodies of the 100 km size class, but are included in our model to enable also the modelling of bigger-than-usual bodies for which they may become important.

Fig. 7 Measured heat conductivities and porosities of L chondrites at 300 K (Yomogida & Matsui 1983; Opeil et al. 2010, 2012) for which both quantities have been measured. The dashed lines correspond to the heat conductivity according to the conduction model (Appendix A.4) for a porous material with abundant micro-cracks. The indicated values correspond to the reduction factor, s, for phonon scattering length. The blue line correspond to the heat conductivity of a sand-stone like porous material according to Eq. (A.4).

Fig. 8 Importance of 40K as heat source. Upper panel: temperature evolution of the parent body model at a set of selected depths (km) considering the contributions of the full set of heat sources (red lines) for our optimised model, and a model with the same formation time and radius but with the sole contribution of 26Al to heating (grey lines). Also shown for comparison are Hf-W and Pb-Pb ages and the corresponding closure temperatures for the two indicated meteorites. Lower panel: cumulated heat released, Q, up to instantt for the different contributing radioactive species.

6.3 Contribution of different heat sources

Figure 8 demonstrates the effect of including other heat sources other than ²⁶Al by comparing the temperature evolution at some selected depths below the surface for our preferred model LM2 with that in a model where only ²⁶Al is considered. During the cooling phase of the body, there are moderate but significant differences of the instant at which the temperature at a depth drops below a given temperature, for example, the closure temperature of a thermochronological system. The difference between the two models is essentially due to the contribution of ⁴⁰K; the contribution of U and Th becomes only important at times when temperatures dropped already well below the closure temperatures of the thermochronological systems of interest, and ⁶⁰Fe is never important at the presently favoured low abundance in the solar nebula. For accurate thermochronological data (like those we used) the accuracy with which closure times for a specific radiometric clock are determined are much better than the difference between the instances where the cooling curve for a meteorite passes through the corresponding closure temperature for different parent body models that either include heating by ⁴⁰K or not. If detailed models for thermal evolution of asteroids are compared to such accurate data, the contribution of ⁴⁰K has to be considered. If, on the other hand, more qualitative models based on a solution of the heat conduction equation with constant coeffcients are considered, it suffices to consider only ²⁶Al as heat source.

6.4 Abundance of petrologic types

It is generally assumed that the metamorphism of the chondritic material from the highly pristine material seen in L3 meteorites to the highly equilibrated material seen in L6 meteorites is essentially determined by the peak temperature achieved by the material during the thermal evolution of the parent body. For the given maximum metamorphic temperature of a particular petrologic type and a temperature evolution model of the parent body, one readily derives the depth of the borders between the zones where the structure of the material corresponds to a given petrologic type and the next higher one, and then the volume percentages of the parent body occupied by the different petrologic types.

In this way Miyamoto et al. (1981) derived the volume fractions of L3 to L6 material in the parent body of the L chondrites and compared this to the observed fall statistics of L chondrites. They used a value of 500 K for the maximum metamorphic temperatureof L3 chondrites, a maximum temperature of 800–850 K for type L4, and a maximum temperature of 1150 K for type L6. The border between type 5 and type 6 was somewhat arbitrarily set to the point of steepest descent of the T_max curve.

The calculated volume fraction of L3 was used by Miyamoto et al. (1981) as additional information to fix the boundary temperature T_b of the model, because the volume fraction of L3 material in the body strongly depends on the surface temperature. The observed volume fraction of 5–10% L3 chondrites from all antarctic L chondrites motivated the choice of a surface temperature of 200 K for their models.

To calculate volume fractions, Bennett & McSween (1996) used maximum metamorphic temperatures of 873, 973, and 1023 K for L3, L4, and L5 chondrites, respectively, from McSween et al. (1988), and of 1273 K for L6 chondrites from Harvey et al. (1993). Mare et al. (2014) used the values for maximum metamorphic temperature as proposed by Harrison & Grimm (2010): 948 K for type 3 and 1273 K for type 6. They do not discriminate between types 4 and 5; the boundary between type 5 and type 6 is assumed to be at 1138 K. The data we adopt for the maximum metamorphic temperatures to calculate volume fractions for our model LM2 are shown in Table 7. The resulting volume fractions for this and the previous models are given in Table 9 where they are compared to the fractional abundance of L3 to L6 chondrites as derived from the numbers of known meteorites of type L3 to L6. For comparison we also show the corresponding numbers forthe case of the H chondrites, the only other meteorite class for which such information is presently available.

In all cases, both for L and H chondrites, no good agreement was found between calculated and observed volume fractions. This is no surprise because the observed relative proportions of the different petrologic types depend on the eventualities of the collision history of large to small bodies and of collision fragments in the planetary system during the last 100 million years or so, and on the diverse biases for recovering meteorites of different types. It is presently not possible to reconstruct the volume percentages of petrologic types 3 to 6 in the parent body from that observed for meteorites.

7 Concluding remarks

We attempted to reconstruct the parent body of the L chondrites from empirical data on the cooling history of meteorites. From the literature we collected the available data on closure times for different thermochronological systems and for cooling rates. The number of useful data found for L chondrites is far less than for the much better studied H chondrites, despite the fact that L chondrites are the most frequent ordinary chondrites found. This certainly results from the fact that the parent body of the L chondrites was catastrophically disrupted 470 Ma ago by a collision that is responsible for the high concentration of fossil meteorites in mid-Ordovician marine limestone in southern Sweden (Heck et al. 2004, 2008) and that only a few L chondrites have been found for which the radioactive decay clocks were not reset by this event. In total, high-quality thermochronological data are available for five meteorites for two thermochronological systems, the Hf-W and U-Pb-Pb systems. This is just sufficient observational material to pin down the two most important parameters of the parent body, its radius and its formation time, by fitting models for the internal constitution and evolution of small bodies from the Asteroid belt to the meteoritic record. Some additional data material on the cooling history of L chondrites is available that can be compared with parent body models, but it is of insufficient accuracy to be used to narrow down the parent body properties.

Based on the assumption of a spherically symmetric body that formed instantaneously and is heated by the decay of long- and short-lived radioactives, in particular by ²⁶Al, we solved the heat conduction equation coupled with equations for the sintering of the initial granular chondrule-matrix mixture and for the porosity and temperature dependence of the heat conductivity. Additionally we allowed for a surface layer of reduced heat conductivity due to impact-generated cracks. By varying the radius and formation time (after CAI formation) we determined the parameter combination that fitted the set of high-quality thermochronological data as well as possible.

It turned out that two different types of models are compatible with the available data. One solution is a body with a radius of 115 km and formation time of 1.89 Ma after CAI formation, another solution is a body with a 160 km radius and 1.835 Ma formation time. The basic difference between the alternative models is that for the bigger model the core region of the body shows incipient eutectic melting of Ni,Fe-FeS, while for the smaller model the central temperature remains well below the melting temperature. The lack of primitive achondrites that are related on compositional grounds to L chondrites presently favours the smaller of the two models.

Acknowledgements

We acknowledge J. L. Hellmann, T. S. Kruijer, J. A. Van Orman, K. Metzler, and T. Kleine for making available to us their data on Hf-W ages and cooling of L chondrites prior to publication. This work was performed as part of a project of “Schwerpunktprogramm 1385”, supported by the Deutsche Forschungsgemeinschaft (DFG)”. M.T. acknowledges support by the Klaus Tschira Stiftung gGmbH. This research has made use of NASA’s Astrophysics Data System.

Appendix A Properties of chondritic material

A.1 Porosity

The composition and structure of the material of the parent body is inferred from the properties of the meteorites. The material of the L chondrites is a mixture of millimetre-sized chondrules embedded in a matrix of micrometre-sized dust grains. This material is metamorphosed to varying degrees by the action of high temperatures and pressures due to internal heating of the body and the action of impacts on surface-near material. Its properties form a sequence between a weakly consolidated and highly porous material composed of clearly discernable chondrules and matrix (petrologic type 3) and a highly compacted and re-crystallised material with hardly discernable chondrules (petrologic type 6), with all possible intermediate states between the two extremecases.

It is assumed that the most porous material found in L3 chondrites represents the composition and structure of the material found in surface-near layers of the parent body and essentially represents the initial granular state of the parent body material. The consolidated and crystallised material found in L5 and L6 chondrites is assumed to be the product of sintering and crystal growth at high temperature and moderate pressure during transient heating in deeper layers of the parent body, resulting from the same kind of material as found in L3 chondrites.

The material of chondrites originating from surface-near layers of the parent body, however, does not simply represent the pristine structure of the parent body material. During the ~ 4.6 Ga period between formation of the body and the excavation of a meteorite from the surface of the body or of a major fragment of it, material from the first few kilometres beneath the surface has been subject to many shock waves emanating from impacts of small to large bodies from the Asteroid belt onto the parent body. If the peak pressure in the shock wave exceeds about 5 GPa (Stöffler et al. 1991), existing pore space in the material is reduced or closed (by rolling, gliding, crushing, transient melting), while at the same time numerous micro-cracks are freshly generated in the mineral crystals.

In orderto determine the initial porosity of the parent body material, we have to look for the highest porosities observed in meteorites that show the least indications of shock metamorphosis (shock stage S1, cf. Stöffler et al. 1991)because these also show the highest porosities (Consolmagno et al. 2008). These high porosities are not necessarily found, however, amongst the chondrites of the lowest petrologic type, because they are just those most affected by impact compaction. Intermediate types (L4,L5) are more promising for this, while type L6 meteorites are already heavily compacted by sintering. The investigations of the porosity of meteorites by Consolmagno et al. (1998, 2008) found a broad distribution of porosities extending from nearly zero values up to ~ 0.2–0.25, with some individuals showing porosities up to 0.3. On average, the porosities are highest for the lowest shock stage, as expected. We interpret the findings in Consolmagno et al. (1998) as showing the important role played by impacts and sintering on a material with an initial porosity corresponding to the highest observed values. We take a porosity of Φ =0.25 as the initial value of the average porosity of the pristine material of the L chondrite parent body.

This value of Φ fits well to the porosity of Φ₀ = 0.248 expected for a binary granular mixture consisting of two components with significantly distinct diameters of the granular units of the two components (cf. Gail et al. 2015). For the L chondrites the chondrules have a typical diameter of 0.5 mm with a rather narrow distribution around the mean (Friedrich et al. 2015), while the matrix particles have diameters on the order of a few micro-metres (we chose a value 2 μm). These values are used as initial values for the calculation of the porosity.

A.2 Composition

With respect to the average composition of the main components of the mineral mixture of L chondrites, we assume modal mineral abundances as given in Yomogida & Matsui (1983), which are almost the same as those given in the more recent compilation in Schaefer & Fegley (2007). This composition is used to calculate the heat capacity and heat conductivity of the material.

To calculate the heat production by decay of radioactives we need to know the average mass fraction of the corresponding elements. For the elements Al, K, and Fe we use data for the composition of L chondrites from Table 1 in Schaefer & Fegley (2007) from which we calculate the mass fractions X_i of the elements in L chondrites given in Table A.1. These are used by us to calculate the heat production rate.

The mass fractions of elements calculated for the composition in Yomogida & Matsui (1983) are slightly different from those given in Wasson & Kallemeyn (1988) and from the values calculated from Schaefer & Fegley (2007), but they are usually within 20%, which is less than the uncertainty of the averages. We prefer the values derived from Schaefer & Fegley (2007) because they are based on a broader data set.

A.3 Heat production

For the specific heat production rate, h, one has to consider the contributions of short- and long-lived radioactives, in particular that of ²⁶Al, ⁶⁰Fe, ⁴⁰K, which determine the temperature evolution during the first 200 Ma of evolution. This is the main period of interest for us. For the evolution on longer terms one also has to consider Th and U. The heating rate, h, in the heat conduction equation has the form $$\begin{array}{ll}h\hfill & ={\displaystyle \sum _i}{h}_i{\text{e}}^{-{\lambda }_i(t-t_0)}.\hfill \end{array}$$(A.1)

The time t₀ is the time of solar system formation, which is identified with the time of CAI formation. The quantities λ_i are the decay constants of the nuclei and $$\begin{array}{ll}{h}_i\hfill & =\frac{X_i}{m_i}{f}_i{\lambda }_i{E}_i\hfill \end{array}$$(A.2)

is the specific heating rates as calculated from the composition of the chondritic material and the initial abundances of the radioactive nuclei at time t₀. The decay constants are λ_i = ln 2∕τ_1∕2,i, where τ_1∕2,i is the half-life of the isotopes of interest, X_i is the mass fraction of the corresponding element and m_i its atomic mass, f_i the fraction of the isotope at the time of solar system formation, and E_i the energy released (without neutrino energies) by the decay into the final stable state. The sum runs over the different radioactive nuclear species that are important for heating. The heating rates h_i do not depend on time and temperature. They also do not depend on r as long as no differentiation occurs. The values of the coefficients used in the calculation are shown in Table A.1. The heat production rate may also include other heat sources, for example the latent heat during melting and crystallisation, but apparently this is not of interest for the L chondrite parent body and we neglect them.

A.4 Heat capacity and conductivity

For the kind of investigations that we have in mind, it is important to base a calculation of the temperature evolution on heat capacities, c_p, and heat conductivities, K, that take account their temperature dependencies. The temperature structures of models calculated with constant values of c_p and K deviate significantly from models with realistic temperature variation and porosity dependence, as was already found by Yomogida & Matsui (1984) and demonstrated in detail by Ghosh & McSween (1999) and Akridge et al. (1998), and was also found in our previous model calculations. Even in a simplified model calculation one has to observe this.

A.4.1 Heat capacity

The specific heat capacity c_p(T) of the chondritic mixture of minerals and metal is calculated from the specific heat capacity c_p,i (T) of the components i and their mass fractions X_i in the mixture found in L chondrites as $$c_p(T)={\displaystyle \sum _i}{X}_i{c}_{p,i}(T).$$(A.3)

The heat capacity is calculated for the standard pressure of one bar since for the pressure range of concern for bodies not bigger than a few hundred kilometre radius c_p is practically pressure independent. Further details are given in Henke et al. (2012a). The result is shown in Fig. A.1, which shows for comparison also the result of a laboratory measurement of c_p for T ≥300 K of the L6 chondrite Sołtmany by Wach et al. (2013) and for some L chondrites for T ≤300 K by Macke et al. (2016). We see reasonable agreement between c_p calculated according to Eq. (A.3) and measured heat capacities.

Fig. A.1 Specific heat of L chondrite material, as calculated from average composition (black line), average values of experimental data (red dashed line) according to Macke et al. (2016), and data measured for an L6 chondrite (Wach et al. 2013; blue squares). The discontinuities and the cusp at certain temperatures are due to phase transitions in some components of the mixture.

In the temperature range T > 300 K some discontinuities and a cusp are found in the temperature variation of c_p. These are due to phase transitions (mainly in the FeS component) and the Curie point of iron. The heat capacity shown in Fig. A.1 does not include a contribution of molten components, because we assume that the temperature at the centre of the L chondrite parent body stays below the Fe,Ni-FeS eutectic melt temperature of 1 223 K at the Ni content of L chondrites (Tomkins 2009). In the model calculation, first a table for c_p (T) is generated and during the model calculation the value of c_p(T) at the required temperature is determined by interpolation from the table.

A.4.2 Heat conductivity

The temperature dependence of the heat conductivity of chondrites was studied in Gail & Trieloff (2018). We found that the experimental findings for the heat conductivity of chondrites (Yomogida & Matsui 1983; Opeil et al. 2010, 2012) can be explained rather well by the heat conductivities of the individual components of the mixture if one appropriately considers the influence of isolated pores and in particular the network of abundant micro-cracks in the silicate mineral components. The model accounts for the observation that the observed heat conductivity of meteorites is far below that of the bulk material and it explains the large differences in the observed temperature variation of heat conductivity in meteorites of otherwise similar properties.

Figure A.2a shows the calculated effective heat conduction coefficient of chondritic material according to the model of Gail & Trieloff (2018), with a composition corresponding to L chondrites for the compact material (solid black line), for a material with 25% porosity by macro-pores (solid red line), and for a material where in addition of the scattering of phonons in the silicate minerals is enhanced by numerous micro-cracks (solid blue line). Here we used a modified fit to the heat conductivity data of enstatite. The corresponding coefficients for the Callaway model (see Gail & Trieloff 2018, for their meaning) of phonon heat conductivity are given in Table A.2. This reduces the deviation between the model and the laboratory data that remained in Gail & Trieloff (2018) for the case of the two enstatite chondrites from Opeil et al. (2010, 2012). With this modificationit is possible to obtain a nearly perfect fit of the heat conductivities of L chondrites. Figure A.2b shows measured heat conductivities in the temperature range from 10 to 300 K for three L chondrites from Opeil et al. (2010, 2012, data made kindly available by G. Consolmagno) and the temperature variation according to the model of Gail & Trieloff (2018) if the two basic parameters of the model, the porosity due to isolated pores, denoted by Φ, and that due to micro-cracks, are determined by a fit to the empirical data. The close fit between the model and the measured data suggests that the theoretical model is suited to describe the temperature variation of the heat conductivity over the temperature range between about 20 and 1200 K, which is relevant for the modelling of the thermal structure of ordinary chondrite parent bodies if the porosity, Φ, due to macro-pores, the enhancement factor, s, for phonon scattering (defined in Gail & Trieloff 2018), and the aspect ratio, α, of the pores are chosen appropriately.

Fig. A.2 Heat conductivity of L chondrite material. (a) Heat conductivity by phonon-transport for the bulk material as calculated from its composition (black line), conductivity for a porous material with a 25% porosity (red line), and conductivity for material with increased phonon scattering by micro-cracks (blue line) with phonon scattering length reduced by a factor s = 30. The dashedline shows the contribution of radiative transfer to heat conductivity. (b) Comparison of the theoretical model (Gail & Trieloff 2018) for L chondrites (black lines) to measured heat conductivities of three L6 chondrites from Opeil et al. (2010, 2012; blue circles) by fitting the model parameters’ porosity (due to isolated pores) and scattering enhancement factor (accounting for micro-cracks).

The most important factor determining the strong scatter of observed heat conductivities was found to be the reduction of phonon-scattering length by micro-cracks. Such micro-cracks are abundant in the present-day meteorites because of the many Gigayears long history of weak to strong impacts on the parent body. It is unlikely that they were already abundant during the early evolution period of asteroids, because during planetesimal formation mass is either acquired by low-velocity (< 1 km s⁻¹ relative velocity) planetesimal–planetesimal collisions, if this is the dominating formation route, or by low-velocity acquisition of pebbles, if this is the dominating formation route. In our model calculation we assume that initially only the macro-porosity due to the formation of the parent body from a granular material of chondrules and matrix is responsible for a reduction of the heat conductivity compared to that of the completely compacted material. In Henke et al. (2016) it was found that the heat conductivity in porous material with macro-pores and a porosity not bigger than ~ 0.3 can be approximated by the expression $$K(T,\Phi )=K_{\text{b}}(T)\left(1-2.216\Phi \right),$$(A.4)

where K_b is the heat conductivity of the bulk material (Fig. A.2a, the black line). This relation is used to calculate the heat conductivity in Eq. (1). For the initial porosity Φ₀ = 0.25, this means for instance that the heat conductivity is reduced to 45% (the red curve in Fig. A.2a) of its value for the compacted material (the black curve in Fig. A.2a). The gradual vanishing of this macro-porosity at elevated temperatures due to sintering is calculated as part of the model calculation. In the model calculation first a table for K_b (T) is generated and during the model calculation the value of K_b at the required temperature is determined by interpolation from the table.

A.4.3 Regolith layer

After formation of Jupiter, relative velocities between planetesimals are pumped up to several km s⁻¹ (e.g. Davison et al. 2013) and such collisions would damage the mineral structure and generate abundant micro-cracks at the impact locations. This requires peak pressures ≳5 GPa at impact (Stöffler et al. 1991). By this a surface layer of some thickness d builds up where the material structure is strongly modified by impacts. In this layer, the heat conductivity probably approaches the low values observed for present day asteroids (the blue curve in Fig. A.2a).

However, even then it would be a considerable time until most of the surface of the body is affected by this and the temperature evolution of the body is modified by the presence of a thermally insulating layer. Unfortunately, there seems to exist no model for the evolution of surface cratering for asteroids over time, in particular over the asteroid’s early phases, except for the outdated model of Housen & Wilkening (1979), from which one could derive estimates for the evolution of the density of micro-cracks. The model of Housen & Wilkening (1979) predicts > 50 Ma for evolution of a significant surface cratering, which seems to be generally in line with the model results in Davison et al. (2013). Therefore, it seems questionable whether the micro-structure of present-day meteorites is representative of the micro-structure of the surface material of planetesimals during the crucial thermal evolutionary phase of the first ~ 100–200 Ma, where temperatures drop below the closure temperatures of the thermochronological systems.

In view of this uncertainty we consider the effect that a surface layer with a reduced heat conductivity due to micro-cracks has on the thermal evolution model by using in Eq. (A.4) for K_b (T) the heat conductivity for impact-damaged material for a radius range between the surface and a depth d below the surface. To calculate this conductivity we assumed a reduction factor of s = 30 and a value of the aspect ratio of the micro-cracks of α = 25 (see Gail & Trieloff 2018). The heat conductivity for this case is shown by the blue line in Fig. A.2a. The value of d is determined such that we obtain as good a fit as possible of the model results as compared to the meteoritical record.

A.5 Sintering

The porous granular mixture of minerals and the metal from which the parent body of the meteorites of an ordinary chondrite group initially is composed is compacted to a state with vanishing void space by creep processes and surface diffusion at elevated temperatures and/or high pressures (sintering). At the same time its mineral composition is equilibrated. We include the process of compaction in our thermal evolution model of theL chondrite parent body in the same way as described in Gail et al. (2015) for the case of the evolution of the H chondrite parent body. The calculation is based on a theory of isostatic hot pressing used for the modelling of sintering in technical processes. All necessary equations and input data are described in Gail et al. (2015). We have only to observe that the average chondrule diameter in L chondrites is 0.5 mm, that is, bigger than the average diameter of 0.3 mm for chondrules in H chondrites. Because the time required for sintering a granular material strongly increases with the size of its granular units and with decreasing temperature, this means that the compaction of chondrules in L chondrite material requires a higher temperature than for chondrules in H chondrite material, to compensate for the bigger size.

The sintering operates in the temperature range above 600 to ≳1000 K, depending on the kind of material and particle sizes. Since temperatures close to the surface of a body always stay below this threshold temperature, there remains an outer layer of considerable thickness with barely compacted porous material. It is important to include the compaction process into the modelling of the thermal evolution of the body because from the compaction at high temperatures there results a structure of the body where an extended compacted core with high heat conductivity is mantled by a layer of residual in-compacted material with low heat conductivity. Since the study of Yomogida & Matsui (1984) on compaction of parent bodies of chondrites and even more since the paper of Akridge et al. (1998), it is clear that this core-mantle structure significantly modifies the thermal structure of such bodies. In such a case one has a nearly isothermal core and a surface layer with a rather strong temperature gradient, as opposed to the case where heat conductivityis constant or varies only gradually across a body. In this case the temperature varies uniformly between surface and centre.

The variation of porosity Φ of a granular material is described by Eq. (3), where the right-hand side contains contributions from different mechanisms resulting in a compaction of the material (see Gail et al. 2015, for details). For simplicity we consider in the model calculation only the contribution of the creep process, which is responsible for the deformation and compaction of the chondrule aggregate.

We also consider in part of the calculations an approximation where we do not model the sintering process explicitly, as is done for part of our model calculation, but use the simple prescription that $$\Phi =\{\begin{array}{ll}{\Phi }_0\hfill & \text{if}{T}_{\text{max}}\le T_{\text{sl}},\hfill \\ \frac{T_{\text{sh}}-T}{T_{\text{sh}}-T_{\text{sl}}}{\Phi }_0\hfill & \text{if}{T}_{\text{sl}}<T_{\text{max}}\le T_{\text{sh}},\hfill \\ 0\hfill & \text{if}{T}_{\text{sh}}\le T_{\text{max}},\hfill \end{array}$$(A.5)

where T_max is the maximum temperature experienced by a mass element during its past history. This saves much computing time. The prescription is guided by experience with complete numerical simulations of the sintering process (Gail et al. 2015) that in models for chondrite parent bodies the loss of porosity by sintering occurs over a rather narrow temperature interval. The prescription mimics this by switching continuously between the limit cases of low conductivity for porous material and high conductivity for compacted material within the temperature interval from T_sl to T_sh. A similar approach has already been used by Hevey & Sanders (2006).

Fig. A.3 Lines of constant time required for compaction of chondritic material of H chondrites at different temperature and pressure conditions. The upper group of lines corresponds to the sintering of chondrules with 0.5 mm diameter. The middle group corresponds to the sintering of matrix material with particle diameters of 2 μm. The numbers at the lines show the time required for compaction (in Ma). The lower group corresponds to the sintering of fine dust material with particle diameters of 0.2 μm. The blue dashed line corresponds to the highest temperatures achieved at a radius during the evolution of model LM3 for the parent body of the L chondrites and the pressure at that depth. The crosses correspond to radial positions with a 1.5 km distance from each other. The horizontal dash-dotted lines correspond to the incipient melting of a Fe,Ni-FeS eutectic at 1223 K for L chondrites (Tomkins 2009) and that of silicates at 1440 K.

To justify this approximation, Fig. A.3 shows lines of constant time required for sintering of chondritic material from an initial state of random close packing with Φ ~ 0.34 to a state with Φ = 0.05 where the pore space loses its connectivity. The lines are constructed by integrating the equations for porosity evolution during hot pressing given in Gail et al. (2015) from the initial to the final state at constant T and P for chondrules, matrix material, and for fine dust material. Most important for us is the compaction of the chondrules that occurs by creep processes. The figure also shows the pressure-temperature stratification for the optimised model of the parent body of the L chondrites. This line intersects the line for a duration of the sintering process of 10⁶yr at a temperature of ~1060 K. This duration corresponds to the width of the temperature peak during the time evolution of temperature in the surface-near layers, which is essentially ruled by the decay-time of ²⁶Al.

The matrix already re-crystallises and sinters at lower pressures and temperatures by surface diffusion (Gail et al. 2015). According to our sinter model, the initially fine-grained matrix starts to coarsen and sinter at temperatures between 520 and 550 K (see Fig. A.3). This does not result in a significant compaction of the material, but results in a baking together of the chondrule-matrix mixture into a sandstone-like material that shows significant cohesion, different to the sand-like character of the material if no sintering of the matrix material occurred. The effective porosity of the mixture would change if matrix material in the contact region between chondrules sinters, but the effect must be small for the parent bodies of ordinary chondrites because of the low matrix abundance (~ 10 vol%). We neglect this. In our model calculation we chose the following constants in Eq. (A.5): T_sl = 1030 K and T_sh = 1090 K for chondrule sintering, and T_sl = 540 K and T_sh = 560 K for matrix sintering.

References

All Tables

All Figures

	Fig. 1 Variation of χ2 with radius R of the body and birth time tb for fixed thickness d of the outer impact-damaged layer. Left part: shape of the χ2 surface. Alsoshown are contour lines of constant values of χ2 on the χ2 surface and their projection onto the basis plane. Right part: contour lines of χ2 in the R-tb plane.
In the text

Fig. 2 Temperature evolution at different depths below the surface (grey lines) and evolution of the central temperature (dashed grey line). The symbols with error bars correspond to the closure temperatures and closure times of L chondrites. The black lines correspond to the temperature evolution at the burial depths of the five meteorites from Table 6 for model LM2 (a) and model LM3 (b). The data of Homestead and Saratov are not used for the optimisation. Grey dots show chondrule ages for L chondrites from Pape et al. (2019).

In the text

	Fig. 3 Temperature evolution at the burial depths of the meteorites and the corresponding cooling rates (in K per Ma). The dots with error bars show the empirical data. Red dots show the data used for the optimisation, blue dots show the cooling rates.
In the text

Fig. 4 Radial variation of maximum metamorphic temperature and heat conductivity. (a) Maximum of temperature at radius r during the undisturbed thermal evolution of the body (black line) for model LM2. The two broad red strips indicate thetemperature range over which the chondrule assembly (upper strip) and the initially fine-grained matrix material (lower strip) are compacted. The hatched area is the temperature range between the upper metamorphic temperatures of L4 and of L5 chondrites (973 to 1023 K, respectively, according to McSween et al. 1988) where the inter-granular voids between chondrules disappear. The upper blue line shows the eutectic melting temperature of 1223 K for the Ni,Fe-FeS system at the Ni content of L chondrites. (b) Radial variation of the heat conductivity K at 30 Ma after formation of the body.

In the text

	Fig. 5 As Fig. 4, but for model LM3.
In the text

	Fig. 6 Schematic sketch of the onion-shell structure of the L chondrite parent body.
In the text

Fig. 7 Measured heat conductivities and porosities of L chondrites at 300 K (Yomogida & Matsui 1983; Opeil et al. 2010, 2012) for which both quantities have been measured. The dashed lines correspond to the heat conductivity according to the conduction model (Appendix A.4) for a porous material with abundant micro-cracks. The indicated values correspond to the reduction factor, s, for phonon scattering length. The blue line correspond to the heat conductivity of a sand-stone like porous material according to Eq. (A.4).

In the text

Fig. 8 Importance of 40K as heat source. Upper panel: temperature evolution of the parent body model at a set of selected depths (km) considering the contributions of the full set of heat sources (red lines) for our optimised model, and a model with the same formation time and radius but with the sole contribution of 26Al to heating (grey lines). Also shown for comparison are Hf-W and Pb-Pb ages and the corresponding closure temperatures for the two indicated meteorites. Lower panel: cumulated heat released, Q, up to instantt for the different contributing radioactive species.

In the text

	Fig. A.1 Specific heat of L chondrite material, as calculated from average composition (black line), average values of experimental data (red dashed line) according to Macke et al. (2016), and data measured for an L6 chondrite (Wach et al. 2013; blue squares). The discontinuities and the cusp at certain temperatures are due to phase transitions in some components of the mixture.
In the text

Fig. A.2 Heat conductivity of L chondrite material. (a) Heat conductivity by phonon-transport for the bulk material as calculated from its composition (black line), conductivity for a porous material with a 25% porosity (red line), and conductivity for material with increased phonon scattering by micro-cracks (blue line) with phonon scattering length reduced by a factor s = 30. The dashedline shows the contribution of radiative transfer to heat conductivity. (b) Comparison of the theoretical model (Gail & Trieloff 2018) for L chondrites (black lines) to measured heat conductivities of three L6 chondrites from Opeil et al. (2010, 2012; blue circles) by fitting the model parameters’ porosity (due to isolated pores) and scattering enhancement factor (accounting for micro-cracks).

In the text

Fig. A.3 Lines of constant time required for compaction of chondritic material of H chondrites at different temperature and pressure conditions. The upper group of lines corresponds to the sintering of chondrules with 0.5 mm diameter. The middle group corresponds to the sintering of matrix material with particle diameters of 2 μm. The numbers at the lines show the time required for compaction (in Ma). The lower group corresponds to the sintering of fine dust material with particle diameters of 0.2 μm. The blue dashed line corresponds to the highest temperatures achieved at a radius during the evolution of model LM3 for the parent body of the L chondrites and the pressure at that depth. The crosses correspond to radial positions with a 1.5 km distance from each other. The horizontal dash-dotted lines correspond to the incipient melting of a Fe,Ni-FeS eutectic at 1223 K for L chondrites (Tomkins 2009) and that of silicates at 1440 K.

In the text