© ESO, 2012

1. Introduction

According to our present understanding of the formation of terrestrial planets, planetesimals from the size class of 100 km form an important transition state during the formation of protoplanets and planets (e.g. Weidenschilling & Cuzzi 2006; Nagasawa et al. 2007). These bodies are sufficiently big that heat liberated by decay of short lived radioactive nuclei, e.g. ²⁶Al, does not easily flow to the surface and radiates away. Instead they heat up and the bigger ones start to melt in their core region. The pristine dust material from the accretion disk undergoes a number of metamorphic processes in this way until it evolves into planetary material. A few such bodies that represent this intermediate stage of the planetary formation process have survived in the asteroid belt and material from these bodies is available on Earth as meteorites. Meteorites preserve in their structure and composition quite detailed information on the processes that were at work during planetary formation. Recovering this information from these witnesses of the early history of our solar system requires the modelling of the structure, composition, and thermal history of meteorite parent bodies to put the analytical results of laboratory investigations of individual meteorites into a more general context.

There are fundamental differences in the composition of ordinary and carbonaceous chondritic meteorites that are related to the absence or presence of ice, respectively, in the region of the protoplanetary disk where their parent bodies formed. We aim to study the parent bodies of ordinary chondrites. They were essentially ice-free, and are in this respect more similar to planetesimals in the inner solar system, where the terrestrial planets formed. In particular we concentrate on the H and L chondrites that form two fairly homogeneous classes and seem each to descend from a single parent body in the asteroid belt.

The model calculations for the thermal evolution of asteroids made so far were reviewed by McSween et al. (2003) and Ghosh et al. (2006). Many of the calculations were based on the analytical solution of the heat-conduction equation with constant coefficients for a spherical body heated by a homogeneously distributed and exponentially decaying heat source (²⁶Al) given by Miyamoto et al. (1981). These models have the advantage of being simple and easily applicable for discussions of special problems, but they cannot be extended to more realistic cases where material properties like heat-conduction, specific heats and others are neither spatially nor temporally constant, or to include additional physics. It was possible, however, to show that H and L chondrites of petrologic classes 3 to 6 originate from bodies of about 100 km size, heated by radioactive decay of short-lived nuclei, and that the different petrologic classes correspond to layers at different depths within the same body, which experienced different thermal histories during the initial heating and subsequent cooling of the body (Göpel et al. 1994; Trieloff et al. 2003; Kleine et al. 2008; Harrison & Grimm 2010).

If more physics is to be included in a model, a numerical solution of the basic set of equations is necessary. One of the first models of this kind was the model calculation by Yomogida & Matsui (1984) who used data for the temperature and porosity dependence of material properties which they determined by laboratory measurements on H and L chondrite material (Yomogida & Matsui 1983). A central point of this model calculation was that it addressed for the first time quantitatively the process of sintering of the initially strongly porous material under the action of pressure resulting from self-gravity once the body is heated to high temperatures ( ≳ 700 K). As a result, the body develops a strong zoning of the material properties: a highly consolidated core with high thermal conductivity if temperature and pressure become sufficiently high during the course of the thermal evolution of the body, and a porous outer layer with low heat-conductivity where temperatures and pressures never climbed to the level where compaction by sintering becomes possible.

This model shows distinct basic differences to a model with constant coefficients like that of Miyamoto et al. (1981) because it has an extended inner region with only a slow outwards drop of temperature where heat-conductivity is high, and a shallow outer layer where the temperature rapidly drops to the outer boundary temperature. The results for the depths of formation of different petrologic types and the predicted spatial extent of these zones are very different in the models of Miyamoto et al. (1981) and Yomogida & Matsui (1984). This demonstrates the importance of including the consolidation of granular material by “hot pressing” and the resulting changes of thermal conductivity in the models. Unfortunately, the results of Yomogida & Matsui (1984) are somewhat impaired because the important heating source ²⁶Al was not included in the calculation. In the sequel the sintering processes seem to have been included in the modelling of thermal evolution of asteroids only in the calculations of Akridge et al. (1997) and of Senshu (2004). Other studies of the blanketing effect of porous outer layers on thermal evolution (Akridge et al. 1998; Hevey & Sanders 2006; Sahijpal et al. 2007; Gupta & Sahijpal 2010) assume the consolidation of these layers at some characteristic temperature, but do not include this process as part of the model calculation.

We study the thermal history of the parent bodies of ordinary chondrites using new data for the thermal conductivity of granular material (Krause et al. 2011a,b) and on the compaction by cold isostatic pressing before the onset of sintering (Güttler et al. 2009). The modelling of the sintering process is based on the same kind of theory (Rao & Chaklader 1972) as in Yomogida & Matsui (1984). This theory does not correspond to the more recent approaches used for modelling of hot pressing in technical processes (e.g. Arzt et al. 1983; Fischmeister & Arzt 1983; Larsson et al. 1996; Storåkers et al. 1999), but it seems more appropriate for the lower pressure regime encountered in asteroids.

The main purpose of this paper is to develop improved models for the thermal history of parent bodies of ordinary chondrites and to compare the model results with cooling histories of H and L chondrites of different petrologic types determined from thermochronological methods in cosmochemistry. This allows us to better constrain the size of the parent bodies and their formation times.

2. Thermal history of chondrite parent bodies

The heat source for early differentiation and metamorphism of meteorite parent bodies was decay heat of short-lived nuclides like ²⁶Al or ⁶⁰Fe (Göpel et al. 1994; Trieloff et al. 2003; Kleine et al. 2008; Bouvier et al. 2007). If planetesimals are larger than tens of km, the maximum degree of internal heating is given by the initial ²⁶Al abundance, i.e., mainly formation-time. Planetesimals formed at  ≈2 Ma after calcium-aluminium-rich inclusions (CAIs – commonly regarded as oldest objects in the solar system), heat up without differentiation, yielding thermally metamorphosed chondritic parent bodies (ordinary chondrites, enstatite chondrites, or more strongly heated acapulcoites and lodranites). Primitive chondrites (such as petrologic type 3) can survive in the outer cool layers of larger parent bodies (following the onion shell model, see, e.g., Göpel et al. 1994; Trieloff et al. 2003), or in bodies that never grew larger than 10 − 20 km in size (Hevey & Sanders 2006; Yomogida & Matsui 1984), or in bodies that formed relatively late.

These scenarios can be verified via their thermal history and cooling paths. These paths can be evaluated by a variety of thermochronological methods applying chronometers with different closure temperatures, which means the temperature where no parent or daughter nuclide is lost from their host minerals. The Hf-W dating method (e.g. Kleine et al. 2005, 2008) is useful for cooling below 1050 to 1150 K. The closure temperature for the U-Pb-Pb system in phosphates is  ≈720 K (e.g. Göpel et al. 1994; Bouvier et al. 2007), Ar-Ar ages reflect the cooling below  ≈550 K for oligoclase feldspar. The annealing temperature for ²⁴⁴Pu fission tracks (e.g. Pellas et al. 1997; Trieloff et al. 2003) in orthopyroxene is 550 K, for the phosphate merrillite  ≈ 390 K. The latter yield a relative cooling rate between 550 and 390 K for the respective rock (e.g. Pellas et al. 1997; Trieloff et al. 2003). Another method to determine cooling rates between 800 and 600 K are metallographic cooling rates, which use Fe-Ni diffusion profiles in metal grains consisting of the two different Fe/Ni phases kamacite and taenite (e.g. Herpfer et al. 1994; Hopfe & Goldstein 2001). These data are the basis for the simulation of the cooling of chondritic parent bodies in Sect. 8.

3. Composition of planetesimal material

3.1. Mineral composition of planetesimal material

From the pristine dust material in the accretion disk, which is formed simultaneously with the formation of a new star, finally planetesimals are formed by a complicated agglomeration process. This process is not modelled in our calculation (for a review on this part of the problem see, e.g., Nagasawa et al. 2007). We concentrate here on the evolution of the internal constitution of the planetesimals, the successors of which can still be found in the asteroid belt. In particular we concentrate on the class of planetesimals that are the parent bodies of the H and L chondrites. Their composition is well known from studies of meteorites. The matrix and chondrules that form the material of the ordinary chondrites consist of a mixture of minerals that themselves in most cases are solid solutions of a number of components. We considered an average composition for the material contained in the parent bodies of the H and of the L chondrites that concentrates on the few main constituents that represent almost 100% of the matter. The many more additional compounds with small abundances found in the material are not important for the bulk properties of those materials that determine the constitution of the planetesimals and their evolution.

The pure mineral components used in our model for the composition of planetesimals and their mass-densities are listed in Table 1. The mixture of minerals that is assumed to form the planetesimal material is listed in Table 2. The composition is given in the notation where, e.g., Fo₈₀Fa₂₀ denotes that 80% of all basic building blocks of the mineral correspond to forsterite and 20% to fayalite. The quantity X_min denotes the mass fraction of the component in the mixture. The composition of the mixture components and the fractions X_min are taken from Yomogida & Matsui (1983); the average composition of chondrites given by Jarosewich (1990) is essentially the same, only the Fe and Al mass-fractions are slightly lower than in Yomogida & Matsui (1983).

The density of the solid solutions given in the table is simply calculated as the average of the densities of the pure components, weighted with their mole fractions in the solid solution. This assumes that non-ideal mixing effects are too insignificant to be considered for our calculations. The average density ϱ_bulk of the mixture is calculated from (1)This is the density of the consolidated meteoritic material, i.e., without pores.

Table 2 presents mass-fractions X_el of the elements for which radioactive isotopes of sufficient abundance exist that their energy release during decay contributes significantly to the heating of the planetesimals. The element abundances used for the calculation are from Lodders et al. (2009).

3.2. Granular nature of the material

The main constituents present in chondritic meteorites are chondrules and a matrix of fine dust grains. The typical relative abundance of these constituents is shown in Table 3; they have typical sizes of micrometer in the case of relatively unprocessed matrix, and typical sizes of millimetres in the case of chondrules, with variations specific to each chondrite group (Scott 2007). Though chondrules may also contain fine-grained mineral assemblages, they entered the meteorite parent body as solidified melt droplets after they formed by unspecified flash heating events in the solar nebula.

The size of the fine dust grains of the matrix before compaction and sintering of the meteoritic material is not known. The observed sizes of matrix particles may be not representative because coarsening may have occurred during heating of the parent body (Ostwald ripening). Probably the very fine granular units in interplanetary dust particles (IDPs) indirectly indicate the pristine size of dust particles in protoplanetary discs before incorporation into the forming bodies of the solar system. Interplanetary dust particles seem to represent cometary dust (Bradley 2010) and this seems to be the least processed dust material handed down from the early solar nebula. The study of Rietmeijer (1993) shows that the granular units forming IDPs are sub-micrometer sized, most of them having diameters smaller than 0.5   μm and many being nano-meter sized,  ≪ 0.1   μ, with a small population ranging in size up to several μm. Typical particle radii then are roughly of the order of 0.1   μm ... 0.2   μm.

It is generally assumed that the planetesimal material initially is very loosely packed with a high degree of porosity. As porosity φ we denote the volume fraction not filled by solid dust material. The volume not filled by the solids is called the pore space or the pores. The pores may form an interconnected network of voids and channels at high porosity, φ ≲ 1, or isolated voids at low porosity, φ ≪ 1. The pores may be empty (vacuum) or may be filled with gas or something else. If there are no pores, the material is called compact.

The compact solid material has a density ϱ_bulk. This material is a complex mixture of minerals, which is described in Sect. 3.1. The porous material has a density (2)where ϱ_p is the density of the material in the pores. If the pores are filled with gas, ϱ_p is small enough that we may assume ϱ_p = 0. At sufficiently low temperatures the pores may alternatively be filled by water and ice, in which case ϱ_p may not be neglected, but here we do not consider these cases. Then (3)In addition to porosity φ we also consider the filling factor (4)Then (5)For this reason D is also called the relative density of the porous material.

The porous material may approximately be described as a packing of spheres. With respect to the chondrules this is not unrealistic, because they show nearly spherical shape. For the matrix grains this may be taken as a rough approximation to describe some basic properties of the packing. For a random packing of equal-sized spheres, Arzt (1982) found a relation between the coordination number Z, i.e., the average number of contact points of a particle with its neighbours, with the filling factor or relative density D: (6)Here, Z₀ = 7.3 and D₀ = 0.64. These values refer to the coordination number of a random dense packing of spheres (cf. Jeager & Nagel 1992). The constant C = 15.5 is the slope of the radial density function (distribution of centre distances). The approximation holds up to D > 0.9.

For a packing of equal-sized spheres two critical types of packing exist. One is the random close packing with a filling factor of D = 0.64, i.e., a porosity of φ = 0.36 (Scott 1962; Jeager & Nagel 1992), and an average coordination number Z = 7.3. Its formation requires taping or joggling of the material. The other one is the loosest close packing that is just stable under the application of external forces (in the limit of vanishing force), which has D = 0.56 or φ = 0.44 (Onoda & Liniger 1990; Jeager & Nagel 1992) and average coordination number Z ≈ 6.6. With respect to chondritic material this means that volume filling factors of more than D = 0.64 for chondrules, like that given in Table  3, require that compaction by sintering or compaction by crushing of particles through strong pressure occurred.

For volume filling factors of chondrules between the random close packing (D = 0.64) and the random loose packing (D = 0.56) some pre-compaction by applied forces occurred. The high volume filling factors of chondrules in the material of the chondrites shows therefore that this material was already significantly compressed. It is not possible to reconstruct from this the initial properties of the material in the bodies that contributed to the growth of the parent bodies. We have to estimate this.

At the basis of the planetesimal formation process there is the agglomeration of fine dust grains like that found in IDPs into dust aggregates of increasing size. This agglomerated material from very fine grains that was not subject to any pressure has a porosity as high as φ ≈ 0.8...0.9 (e.g. Yang et al. 2000; Krause et al. 2011a). The number of contacts of a particle to neighbouring particles then is on the average equal to about two, i.e., the fine dust grains form essentially chain-like structures. Such high-porosity material seems to be preserved in some comets (Blum et al. 2006).

Collisions of dust aggregates during the growth process of planetesimals leads to compaction of the material. The experiments of Weidling et al. (2009) have shown that the porosity can be reduced to φ ≈ 0.64 (or even lower, see Kothe et al. 2010) by repeated impacts. This porosity is still higher than the random loose packing and is therefore not guaranteed to be mechanically stable. The average coordination number is Z ≈ 5. Lower porosities of φ ≲ 0.4 were obtained by applying static pressures of more than 10 bar (Güttler et al. 2009). Because the planetesimals form by repeated growth process by collisions with other growing planetesimals and impact velocities can be quite high, we assumed in the model calculations that the initial porosity of the material from which the parent bodies of the asteroids formed is already compacted to some extent, and assumed a porosity of the order of φ = 0.6...0.5.

4. Internal pressure within planetesimals

4.1. Hydrostatic equilibrium

The planetesimals are subject to the mutual gravitational interaction between their different parts, which results in an inward-directed gravitational force at each location. As long as the material is highly porous and has not yet been subject to sintering processes, the granular material may start to flow (by rolling and gliding of the powder particles) and evolve into a state with the densest packing of the granular material. If the planetesimal material is approximated by a random packing of spheres of equal size, its porosity in this state would be about φ = 0.37 (Yang et al. 2000). In this state the forces of the mutual gravitational attraction of all other particles are compensated for by the reactive forces from internal stresses that are built up within the particles as a result of the forces acting between contact points to neighbouring particles. The material can then exist in a state of hydrostatic equilibrium with no relative motion between grains. This equilibrium state, however, is of a labile character because at contact points between grains stresses may be built up that are strong enough that particles may break apart into smaller pieces and some further compaction of the material is possible by application of high pressures. For planetesimal material this kind of compaction is probably only relevant for impact problems. During the gradual build-up of planetesimals, temperatures become already sufficiently high by internal heating for sintering by internal creep to occur before really high pressures are built up by which the material may be crushed.

We assumed in our considerations that during the very initial build-up of planetesimals the initially very loosely packed granular material, under the action of its own gravitational attraction and/or during collisions with other planetesimals, has already evolved by granular flow to a state where all particles have a sufficient number of contacts to neighbours such that additional relative motions are effectively blocked (φ ≲ 0.44). Further compaction then requires considerable pressures because particles have to be broken for this. Additionally, we neglected that the planetesimals may rotate and therefore neglected any deformation of the body resulting from this. Then we may assume that the planetesimals have a spherically symmetric structure and the distribution of internal pressure in the porous dust ball is described by the well-known hydrostatic pressure equation (7)where (8)The density ϱ is the mass-density of the material.

As long as a planetesimal is embedded in an accretion disk, it is subject to the external pressure in the disk. The equation of the hydrostatic pressure stratification in the planetesimal therefore has to be solved with the outer boundary condition at the outer radius R(9)with p_ext being the pressure in the accretion disk. The pressure in the midplane of the solar nebula at a number of instants as calculated from an evolution model of the disk (Wehrstedt & Gail 2002, 2008) is shown in Fig. 1b. This defines the outer boundary condition p_ext. After dissipation of the accretion disk one has to use (10)instead.

In order to estimate the magnitude of pressure let us consider a body with constant density ϱ. We have in this case (11)For a body with radius R and external pressure p_ext integration of the pressure equation yields (12)and then (13)For the central pressure p₀ at r = 0 we have numerically (14)The external pressure in the accretion disk is shown in Fig. 1b and we see that already for bodies as small as 1 km radius the central pressure in the body is significantly higher than the pressure in the accretion disk. In bodies with radii of 10 km and bigger one encounters central pressures of the order of 1 bar and higher.

4.2. Isostatic pressing of the granular material

The behaviour of granular material under pressure may be quite complex. This has been discussed by Güttler et al. (2009) with particular emphasis on impact processes during planetesimal growth. These authors also performed laboratory experiments for the behaviour of fine grained material under the action of static pressure and investigated how the porosity is reduced if the material is compacted by pressing. They derived a relation between the applied pressure p and the relative density (filling factor) D that is observed after the material has come to rest¹(15)This relation only holds for D > 0.15 because D ≈ 0.15 was the initial value of the relative density of the material in the experiments before pressing. Its validity is also limited to pressures where the granular material is not yet compacted to the relative density D ≈ 0.64 of a random densest packing of equal-sized spheres.

The meaning of Eq. (15) is that it describes the porosity φ = 1 − D of a fine powder that experienced the maximum pressure p. With respect to the material in a planetesimal, it gives the porosity distribution within the planetesimal as function of depth for those parts where Eq. (15) predicts a lower porosity than the porosity of the surface material. The latter is determined by the continued impact processes during particle growth that also result in a compaction of the powder material (Weidling et al. 2009), a process, however, that is not described by Eq. (15). The porosity φ_srf of the surface layer material has to be described in a different way. Hence we have to determine the porosity in the planetesimal from the following equation (16)This holds before sintering (see Sect. 6) becomes active. Figure 2 shows in the upper part the variation of relative density D with maximum experienced pressure according to Eq. (16). The lower part shows as an example the distribution of relative density within planetesimals with 3 km to 50 km radius resulting from cold isostatic pressing, if no upper limit for the surface porosity is prescribed. The solutions are obtained by integrating Eq. (7) with Eqs. (5) and (15) as equation of state. This shows that up to planetesimal sizes of about 7 km the material is not substantially compacted by the self-gravity of the body. For planetesimals bigger than 10 km most part of the body is already compacted by cold isostatic pressing caused by self-gravity to about a density close to that of the densest random packing of equal-sized spheres. Only a shallow surface layer of highly porous material remains in these bodies. This is evident, but here a quantitative prediction can be made.

4.3. Exosphere

During the first few million years all bodies in a nascent planetary system are immersed in the protoplanetary accretion disk. Sufficiently massive bodies are able to gravitationally bind some of the gas. The minimum requirement for this is that the escape velocity of a gas particle from the surface of a body exceeds its average kinetic energy of thermal motion corresponding to the disk temperature (17)where R is the planetesimal radius and m_g the average mass of the gas particles. From this one has for the lower limit of the planetesimal radius from which on gravitational bonding of an atmosphere starts to become possible (18)For a body with constant density ϱ this means that the radius has to exceed a radius of (19)in order to gravitationally bind gas from the accretion disk and increase the gas pressure at their surface over the local pressure in the disk. This is only relevant for protoplanets with radii R ≳ 1000 km.

5. Temperature structure of planetesimals

5.1. Heat-conduction equation

For the kind of bodies that we considered, the material has the structure of granular matter. The internal structure of this material is not isotropic and its properties are subject to strong local variations on the scale of particle sizes. As a result the temperature will also show local variations on the same length scales. The particles that form the granular material, however, are very small and they are in particular extremely small compared to typical length scales over which macroscopic properties of the bodies may vary. In this case we may average the microscopic properties of the material and also the temperature over volumes that contain a big number of particles and at the same time have dimensions small to the characteristic scale lengths for changes of the values of variables like average temperature T or average density ϱ. Accordingly we only considered these average quantities, for which we could assume that they are isotropic after averaging over all possible particle orientations.

With this approximation the equation of energy conservation for a spherically symmetric body, expressed as an equation for the temperature T of the matter, averaged over the microscopic fluctuations, is (20)where c_v is the average specific heat of the granular material per unit mass, q_r is the radial component of the average heat flow vector, h is the average source term for heat production or consumption per unit mass, the pdV-term is the compressional work done per unit mass, and the last term is the work done by external forces. It is assumed that the external forces, F_r, are species-independent and that no differential motion between the components of the material occurs (flow of gas or water through the porous matrix is presently not considered). The last two terms have not been considered so far in model calculations for thermal evolution of asteroids because they are negligible if the material is practically incompressible. However, if sintering is considered, the material becomes strongly compressed during the course of evolution and these terms have to be included.

The quantity v_r is the uniform radial velocity of the mixture components. If shrinking of the body by sintering is the only kind of motion of the otherwise stationary structure of the body, the velocity v_r is obtained by differentiation of Eq. (8) for fixed M_r with respect to time. There follows (21)We considered models of planetesimals that are in hydrostatic equilibrium without internal flows. There may be some very slow radial motion of the matter if the material starts to shrink by sintering at high temperatures. This kind of extremely slow motion is completely negligible in the substantial derivatives d/dt = ∂/∂t + v_r∂/∂r in Eq. (20). However, one consequence of shrinking is not negligible: If the body shrinks, the gravitational potential energy decreases and the corresponding amount of energy is transferred to the matter as heat. This is described by the term ϱv_rF_r where F_r is the local gravitational acceleration (22)The term corresponding to the work done by the forces has to be retained. With this approximation we have the following equation for the temperature (23)The variation of ϱ with time is discussed in Sect. 6.

The heat flow vector has contributions from a number of processes. For the solid component of the material there is a contribution from heat-conduction by phonons or, in the case of electric conductors (e.g. iron), from conduction electrons. For a porous material there is also a contribution from the heat-conduction by the gas in the pores. If the material is translucent, one has also to consider a contribution to heat-conduction by radiative transfer. All these processes have the property that their contribution to the total heat flow is proportional to the gradient of the temperature. Generally. the coefficient of proportionality is a second rank tensor, unless if the properties of the material are isotropic, in which case it degenerates to a simple scalar factor. For granular material the local transport properties for heat are by no means isotropic. We assumed, however, that after averaging the average heat flow vector is proportional to the gradient of the average temperature. The radial component of the heat flow vector then takes the specific form (24)with some average heat-conduction coefficient K that is different for each of the different transport processes.

The essential material properties that enter into Eq. (23) are the specific heat-per-unit mass, c_v, the heat-conductivity K, and the heat production. In the following subsections we describe how we determined these quantities for the material of the parent bodies of ordinary chondrites.

The body is also subject to heat and matter exchange with its environment. This is treated by defining appropriate boundary conditions for Eq. (23) that are discussed after our discussion of the material properties.

5.2. Heat capacity of material

The heat capacity c_v of a mineral mixture is simply the weighted sum of the heat capacities of its components. It is calculated in our model calculations from (25)where X_min,i is the mass-fraction of the ith component in the mixture of solid components and c_v,i the heat-capacity-per-unit mass of component i. The quantities X_min,i are given in Table 2. Because we intend to consider bodies of a size of no more than a few 100 km, pressures remain below the kbar-range (see Eq. (14)) and compression of solid material under pressure is negligible because of the low compressibility of minerals. Under these conditions there is no significant difference between c_p and c_v and we may use c_p instead of c_v. For our model calculations data for c_p,i(T) were taken from the compilation of Barin (1995), who gives the heat capacities c_p per mole. These quantities are converted to heat-capacity-per-unit mass by dividing by the mole mass M_i of species i: (26)Values of c_v for the required temperatures are determined by interpolation in the tables for c_p,i. The heat capacity for solid solutions is calculated as the weighted mean of the heat capacities of the pure components taking their mole fractions as weights.

The variation of the specific heat of the mixture is shown in Fig. 3. Because some of the minerals suffer structural transitions at certain temperatures and because c_p of iron has a cusp at the Curie-temperature of 1042 K, the temperature variation of c_p shows some kinks and jumps (cf. also Ghosh & McSween 1999). They might be sources of numerical problems. For comparison the figure also shows the analytical approximation for c_p for the bulk material given by Yomogida & Matsui (1984). This is an alternative if the jumps in the temperature variation of c_p result in numerical problems, but in our calculations it was not necessary to take recourse to this approximation.

5.3. Heating by radioactive nuclei

Next we consider the source term h in Eq. (23). There are essentially two different kinds of sources and sinks of heat within the planetesimal bodies. One source is the energy liberated during the decay of radioactive isotopes of a number of elements. The other one is the consumption of energy during the melting of planetesimal material or the liberation of energy during the solidification of the melt. Melting is not considered because we aim to study parent bodies of undifferentiated meteorites.

The main sources of heat input by radioactives during the early heating of planetesimals and the subsequent cooling phase are decay of ²⁶Al and possibly ⁶⁰Fe (cf. the discussion in Ghosh et al. 2006). More long-lived radioactive nuclei, essentially ²³²Th, ^235,238U, and ⁴⁰K, are responsible for the later long-term evolution of the temperature. For the nuclei decaying by β-decay (²⁶Al, ⁶⁰Fe, ⁴⁰K) the energy of the fast electrons and of the emitted γ-photons is absorbed within the planetesimal material and is converted to heat, while the neutrinos leave the bodies and carry away their part of the energy. For the nuclei decaying by α-decay the whole decay energy is absorbed by the planetesimal material. We assume that no chemical differentiation occurs in the bodies that we consider. Hence, after averaging over the inhomogeneous microstructure of the material, the heat-producing nuclei are homogeneously distributed over the bodies. The heat production rate by these nuclear processes is (27)The sum is over all nuclei that contribute to heating, X_el,i denotes the mass-fraction of the corresponding element in the material of the planetesimals (see Table 2), m_el,i the atomic mass of the element for the isotopic mixture at the time of formation of the solar system, f_i the fraction of the isotope of interest at the time of formation of the solar system, τ_i the decay time-scale for e-fold decrease of abundance of the isotope, and t is the time elapsed since formation of the solar system.

The constants for calculating h_nuc are given in Table 4. The element abundances used for calculating X_el,i for the mineral mixture given in Table 2 are taken from Lodders et al. (2009). Isotopic abundances of K and U at the time of the solar system formation are taken from Anders & Grevesse (1989). The abundance of ²⁶Al is that given by Nyquist et al. (2009). The abundance of ⁶⁰Fe is disputed in the literature. Table 4 gives the probably highest value from Birck & Lugmair (1988) and the lower limit according to Quitté et al. (2007). There are indications that ⁶⁰Fe was not homogeneously distributed in the early solar system (Quitté et al. 2010), which means that the initial ⁶⁰Fe abundance in the parent bodies is not known a priori and is an additional free parameter for the modelling. For the decay time of ⁶⁰Fe the recent revised value for the half-life τ_1/2 = 2.62 ± 0.04 Ma of Rugel et al. (2009) is used.

Figure 4 shows the variation of the heating rate per unit mass with time elapsed since formation of CAI, calculated with the high ⁶⁰Fe abundance (see Table 4 for this). The dominant heating source is ²⁶Al at the time of formation of planetesimals (t ≲ 5 Ma), but ⁶⁰Fe dominates as heat source for an extended period from  ≈ 5 Ma to  ≈ 20 Ma because of the revised ⁶⁰Fe half-life.

For comparison Fig. 4 also shows the contribution of the release of gravitational energy to heating during shrinking of the body, resulting for the model of Sect. 8. For the small bodies that are considered in this paper this heating source is not important (but included in the model).

5.4. Heat-conduction by the porous solid material

For the heat-conductivity K of the chondritic material we used two different types of experimental data. For low porosities from the range 0 < φ ≤ 0.2 we used data measured for a number of ordinary H and L chondrites by Yomogida & Matsui (1983). For high porosities φ > 0.4 we used the data for silica powder derived from laboratory measurements by Krause et al. (2011a). All these measurements were conducted under vacuum conditions to exclude any contribution from heat transport from gas-fillings in the pores.

Figure 5 shows conductivity K plotted versus porosity for chondritic material at T = 300 K. The data for H and L chondrites scatter significantly, and because of the small number of available data points, no obvious systematic difference between the two types of material can be recognized. Therefore we fitted both sets of data with a single analytical approximation (28)with two constants K_b and φ₁. This exponential form enables a reasonable fit of the data. The constant K_b may be interpreted as the extrapolated average thermal conductivity of the bulk material at vanishing porosity, for which we obtain K_b = 4.3 W m^-1 K^-1. For the second constant we chose φ₁ = 0.08 (see also Krause et al. 2011b). In Fig. 5 the data are normalised with the value of K_b. The dashed line running through the data points for the chondrites gives our approximation K₁(φ).

At high porosities Fig. 5 shows the conductivity K for a silica powder that consists of equal-sized spheres with 1.5   μm diameter. By the nature of the experimental method the data of Krause et al. (2011a) do not refer to a well-defined temperature but the heat-conductivity was derived by analysing the cooling behaviour of their sample. Therefore the value of K is some average value over the raise and fall of the temperature in the experiments, which is somewhat above 300 K temperature. The data are fitted with an analytical approximation of the form (29)with two constants a and φ₁ and the same value of K_b as before. This type of exponential dependence on φ allows a reasonable fit of the data points in this case as well. The constants are found to be a = 1.2 and φ₂ = 0.167 (in an earlier version, Krause et al. 2011b, slightly different values of the constants were given). This fit is shown as the second dashed line in Fig. 5.

The experiments of Krause et al. (2011a) are conducted with very fine grained silica powder. This is not the same kind of material as is found in chondrites, but for two reasons it may be considered as a reasonable proxy for chondrite material before strong compaction. First, the mineral mix in chondrites is dominated by silicates and all silicates have similar heat-conduction coefficients. Second, the heat-conduction of very loosely packed material is via the tiny contact regions between the particles. In chondrites part of the material are the big chondrules (size  ≈ 1 mm), but as long as the material is not strongly compacted and the chondrules are well separated by the very small grained matrix (particle sizes  ≲ 0.5   μm), the heat-conduction is obviously governed by the heat flow through the contact points between the tiny matrix particles. In this respect the basic physics of the heat transport in the experiment of Krause et al. (2011a) should be very similar to that in chondrite material before strong compaction.

The powder particles used in the experiments are about a factor of five to ten times bigger than the matrix particles in chondrites (Rietmeijer 1993). Some indications on the influence of particle sizes may be obtained by comparing the results of Krause et al. (2011a) with the results of heat-conduction measurements of Fountain & West (1970) for powders of basaltic particulates that are much coarser grained. Figure 5 shows results for their size separate with average particle size of 50   μm. The variation of K with porosity for the Fountain & West (1970) granular material is very similar to the silica powder used by Krause et al. (2011a), except that the conductivity is lower by a factor of somewhat less then a factor of ten. The particle sizes, on the other hand, are bigger on average by more than a factor of thirty. This suggests that the heat-conductivity measured by Krause et al. (2011a) for the highly porous silica powder underestimates the conductivity of the matrix material of chondrites at the same porosity, but probably not by a big factor. For a power-law dependence of K on particle sizes one may speculate that a conductivity higher by a factor of about three of meteoritic material than for the silica powder may be an appropriate estimation, but without more definite information we worked with the values of Krause et al. (2011a) for our model calculations.

The two fits, Eqs. (28) and (29), for φ < 0.2 and φ > 0.4, respectively, are combined into a single analytical approximation for K(φ) by (30)to smoothly interpolate between the two limit cases, in particular in the intermediate range of porosities 0.2 ≤ φ ≤ 0.4. This approximation is shown as the full line in Fig. 5.

The bulk conductivity K_b in Eq. (28) is temperature-dependent. Figure 6 shows data for K for H and L chondrites as given by Yomogida & Matsui (1983), divided by K₁(φ), given by Eq. (28). The data for K(T)/K₁(φ) show for each of the meteorites a clear systematic variation with temperature. The extent of these variations, however, is much less than the scattering between the different meteorites, which amounts to variations by a factor of about two. The origin of these meteorite-to-meteorite variations is not known, but most likely they have their origin in variations in the composition and structure of the meteoritic material. These variations can presently not accurately be accounted for and will not be considered in our model calculations. Therefore we will also neglect the small variation of K with T and take K_b in our calculations to be temperature-independent.

5.5. Heat-conduction by other processes

The material of the planetesimals is dominated by minerals that are transparent in the far-infrared region. Therefore energy transport by radiation is possible. The energy flux by radiation increases strongly with increasing temperature, and for temperatures of the order of about 800 K and higher the heat flux by radiation becomes non-negligible. The complicated structure of the material (chondrules densely packed in a porous dust matrix) makes it difficult to calculate this from first principles. We therefore follow the proposal of Yomogida & Matsui (1984) to take the results of the laboratory measurements of Fountain & West (1970) of radiative heat-conduction in basaltic powders as an approximation for meteoritic material. The model calculations show that for the bodies of interest (temperature stays below melting temperature) radiative heat-conduction never becomes an important energy transport mode. This is because if temperature becomes high enough for radiative heat-conduction to contribute somewhat to the net heat flux, the onset of sintering above  ≈ 700 K results in a strong increase of heat-conductivity by phonons that outnumbers radiative contributions again.

A possible contribution to heat-conduction by pore gas is small and is neglected in this paper; details of this will be discussed elsewhere.

5.6. Boundary conditions

For solving the heat-conduction Eq. (23) one has to specify an initial condition T(r,t₀) at some initial instant t₀ and boundary conditions, in our case at r = 0 and r = R.

The inner boundary condition is that there is no point source for heat at the centre, which translates into the Neumann boundary condition (31)The temperature T_s of the surface layer of the body is determined by the equilibrium between (i) the energy fluxes towards the surface from the interior and from the outside, and (ii) the energy fluxes away from the surface to exterior space (cf., e.g., Grimm & McSween 1989; Ghosh et al. 2003). This results in a mixed boundary condition at the surface (32)The left-hand side describes the energy flow from the interior towards the surface by heat-conduction. The first term on the r.h.s. is the rate of energy loss by radiation from the surface to exterior space, the second term, F_ext, is the rate of energy gain by outer sources. During the first few million years, F_ext is given by , with T_c being the temperature at the midplane of the disk (see Fig. 1a), if the planetesimals are still embedded in an optically thick accretion disk. After disk dispersal the planetesimal is irradiated by the proto-sun, and the rate of energy input, F_ext, is given by . Here a is the (average) distance to the star and A_surf is the albedo of the planetesimals surface.

Alternatively, one can consider the Dirichlet condition at the surface (33)with some prescribed value T_b. This has been done in a number of published model calculations, where some fixed value for T_b was assumed.

5.7. Initial condition

Large planetesimals with radii of the order of 100 km occur as transition states during the growth from km-sized bodies of the initial planetesimal swarm to protoplanets with sizes of the order of 1000 km. The growth initially proceeds slowly on timescales of a few 10⁵ years until run-away growth commences after the biggest planetesimals reach sizes of about 10 − 20 km (e.g. Weidenschilling & Cuzzi 2006; Nagasawa et al. 2007). During run-away growth the mass rapidly increases within less than 10⁵ years to the size of a protoplanet. This means that the bodies are formed on timescales shorter than the decay time of  ≈ 1 Ma for ²⁶Al and also that they collect most of their mass within a period even much shorter than this. There is not sufficient time available for strong heating by ²⁶Al decay of those planetesimals that contributed to the growth of a 100 km-sized body; most of the heat released by radioactive decay is released after its formation.

Therefore we will base our calculations in this paper on the “instantaneous formation” approximation, where it is assumed that the body is formed within such a short period that all heating occurred after its formation. Within the framework of the instantaneous formation approximation it is appropriate to prescribe the temperature T_c of the disk material at the formation-time of the planetesimals as initial value of the temperature.

Modifications resulting from a finite duration of the growth period will be considered elsewhere.

6. Sintering

The compaction of material in planetesimals is a two-step process. The initially very loosely packed dust material in the planetesimals comes under increasing pressure by the growing self-gravity of the bodies. The granular material can adjust by mutual gliding and rolling of the granular components to the exerted force and evolves into configurations with closer packing. The ongoing collisions with other bodies during the growth process enhances this kind of compaction of the material. This mode of compaction, “cold pressing”, by its very nature does not depend on temperature and operates already at low temperature; it was considered in Sect. 4.2. A second mode of compaction commences if radioactive decays heat the planetesimal material to such an extent that creep processes are thermally activated in the lattice of the solid material. The granular components then are plastically deformed under pressure and voids are gradually closed. This kind of compaction by “hot pressing” or “sintering” is what obviously operated in ordinary chondrite material, and the different petrologic types 4 to 6 of chondrites are obviously different stages of compaction by hot pressing.

Yomogida & Matsui (1984) were the first to perform a quantitative study of this process by applying early theories of sintering developed in material science. We follow here essentially the same approach because more advanced modern theories of hot pressing are developed to model metallurgical processes that apply generally much higher pressures ( ≫ 1 kbar) than what is typically encountered in compaction of planetesimal material ( ≲ 10² bar) and are mainly concerned with the final stages of the process. Because the rate of increase of temperature in planetesimals is very low (of the order of 10^-3 K yr^-1), the creep processes result in finite deformations of the material already at quite low temperatures or pressures, where under laboratory conditions no effect is seen. The more simple early theories of hot pressing seem to fit to the situation better.

6.1. Equations for hot pressing

For describing the sintering process we initially assume a dense packing of equal-sized spheres with initial radius R₀. The packing is sufficiently dense that no higher compaction can be achieved by pressing without crushing the spheres. On average, the individual granular units will touch each other at Z contact points. At sufficiently high pressure and temperature the individual spheres will plastically deform at the contact points by creep processes and contact faces will form between adjacent particles, while the volume of the particle will remain constant. As sintering proceeds, the voids between the spheres become smaller and the sphere centres get closer.

There are two stages for this process. In the first stage the voids form an interconnected network between the granular units. This closes at some stage of the sintering process, and isolated pores remain, that have to close by additional sintering in a second stage. The relative density D at the transition between stages one and two depends on the type of packing. The following approximations are for the first stage.

The first basic assumption of the deformation theory of hot-pressing by Kakar & Chaklader (1967), on which the work of Yomogida & Matsui (1984) is based, is a purely geometrical one. It is assumed that the formation of the contact faces can be conceived as if at each contact point a cap would have been cut off from each of the two contacting grains. Then, for grains of equal radius, the contact areas are circular areas with radius a. It is assumed that all cut-away caps have the same height h and radius a at their base. The volume of one such cap is (34)and its height (35)The granular units then are (by assumption) spheres with Z caps cut off. To conserve the original volume of the sphere, the radius R of such a truncated sphere has to be bigger than the pristine radius R₀. Conservation of volume requires (36)This holds as long as the contact areas do not come into contact with each other. The relation to the relative density D is (Kakar & Chaklader 1967) (37)Here D₀ is the relative density of the initial packing of spheres with radius R₀. For a given number of contact points Z and given D₀, R₀, Eqs. (34) to (37) define R and a in terms of the relative density D.

In the theory of Rao & Chaklader (1972) a number of regular packings of spheres was considered for which the number of contact points Z is fixed (Z = 6,8,12). In particular the authors favoured the “hexagonal prismatic” packing with Z = 8 and gave their formula for this case. This is the model that has been used by Yomogida & Matsui (1984) in their study of sintering of planetesimals. These authors argued that many experiments on the packing of small spherical particles of constant size show that the porosity achieved after sufficient tapping would be near 40%, with an average of about eight contact points per particle. Because a regular hexagonal prismatic packing of spheres also has a coordination number of eight and a porosity of 39.5%, and a random close packing has a porosity 36% and on average Z = 7.3 (see Sect. 3.2), they used that packing model in their sintering models. For a discussion of the more general case of random packings see, e.g., Arzt (1982); Arzt et al. (1983); Fischmeister & Arzt (1983); the equations obtained in that case are more involved, but there are no basic differences.

The second basic assumption in the theory of Rao & Chaklader (1972) of hot pressing is that the strain rate is related to the applied stress by the relation for power-law creep (38)and that is given in terms of the rate of change of relative density as (39)The stress σ₁ is the pressure acting at the contact faces of the granular units. The quantities A and n have to be determined experimentally for each material. The quantity A depends on temperature.

The stress σ₁ is given by the pressure acting at the contact areas between the granular units. It is assumed that this is given in terms of the applied pressure p and the areas of contact faces, πa², and average cross-section of the cell occupied by one granular unit, C_av, as (40)In the hexagonal prismatic packing model favoured by Yomogida & Matsui (1984), the cross-section C_av is given by (41)Via the dependence of R and a on D this is a function of D. Values for other packing models (e.g. Kakar & Chaklader 1967) are within O(1) of this. In the initial stages of sintering (small a) the stress σ₁ is much higher than p, in the final stages it approaches p.

Equations (38)–(40) result in the following differential equation for the relative density (42)With the above relations between R, a and D this is a (closed) set of equations for calculating the time evolution of D, which has to be solved together with the other equations for the structure and evolution of the planetesimal which define the pressure and temperature.

The transition to stage 2 by closure of voids is assumed in models of hot pressing to occur at D ≳ 0.9 (e.g. Arzt et al. 1983). The equation for Ḋ becomes Ḋ ∝ 1 − D for this case (cf. Wilkinson & Ashby 1975; Arzt et al. 1983). Because the corresponding full equations are similar in structure to Eq. (42) except for the factor 1 − D, we included the final pore-closing stage simply by multiplying for D > 0.9 the r.h.s. of Eq. (42) by a factor (43)to obtain a continuous transition between both cases.

6.2. Data for olivine

The pre-factor A and the power n in Eq. (38) have to be determined by laboratory experiments. Yomogida & Matsui (1984) used data from Schwenn & Goetze (1978) for olivine. No other data for materials of interest for planetesimals seem to have been determined since then. Schwenn & Goetze (1978) gave the following fit to their experimental data for small spheres of olivine (R < 53,μm): (44)with σ stress on contact faces (in bar), E_act the activation energy for creep, T the temperature, R_gas the gas constant, and R the radius of granular units (in units cm).

For the activation energy a value of E_act = 85 ± 29   kcal   mol^-1 was given by Schwenn & Goetze (1978). In the model calculations we used the value E_act = 85   kcal   mol^-1. For the pre-factor A a range of values from 1.6 × 10^-5 to 5.4 × 10^-5 was given by Schwenn & Goetze (1978); as a compromise we used a value of A = 4 × 10^-5 in our model calculations.

Note that Yomogida & Matsui (1984) chose to use for the approximation with n = 1 given in Eq. (7) in Schwenn & Goetze (1978), while we prefered to use the approximation given by Eq. (8) of Schwenn & Goetze (1978), because these authors explicitly state that this describes their measured σ-dependence.

7. Results for thermal evolution of planetesimals

7.1. Model calculation

The calculation of a model requires solving the differential equations for heat-conduction, Eq. (23), hydrostatic equilibrium, Eq. (7), for the evolution of porosity, Eq. (44), together with equations for the material properties, the equation of state, Eq. (3), the equations for heat-conductivity, Eq. (30), and heat capacity, Eq. (26), and together with appropriate initial and boundary conditions.

The heat-conduction equation and the pressure equation are re-written in terms of M_r, defined by Eq. (8), as independent variable and are discretised for a set of fixed mass-shells with masses ΔM_i (i = 1,I) and shell boundaries r_i (i = 0,I). The M_r-coordinate corresponds to a Lagrangean coordinate that is fixed to the matter. For this choice of coordinates there is no flow of matter across cell boundaries. This enables a simple treatment of growth of the body, if this is considered, and it avoids problems with numerical diffusion in case of inhomogeneous composition (e.g., radial variation of porosity).

The heat-conduction equation is solved by a fully implicit finite difference method with Neumann boundary condition, Eq. (31), at the centre and a Dirichlet boundary condition, Eq. (33), at the surface. The first-order ordinary differential equation for φ, Eq. (44), is solved by the fully implicit Euler method.

To account for the non-linear coupling between the different equations, we performed a fixed-point iteration where we solved the equations in turn as follows:

• 1.

  Given are values of φ_i, T_i, ΔM_i for each mass shell i at some instant t_k − 1. We have to calculate new values at the next instant t_k = t_k − 1 + Δt. The values of φ_i, T_i at t_k − 1 were used as starting values for the iterative calculation of φ_i, T_i at t_k.

• 2.

  The heat production by radioactive decays over the period Δt was calculated for each shell.

• 3.

  For a given porosity φ_i one finds ϱ_i from Eq. (3), and with a given mass ΔM_i we calculated in each mass shell i the shell boundaries r_i at t_k, starting from the centre.

• 4.

  From the change of r_i over time Δt we found the grid velocity and the heat production by release of gravitational energy for each shell (last term in Eq. (23)).

• 5.

  We calculated the pressures p_i at the shell boundaries from the discretised pressure equation, starting with the given external pressure at the surface.

• 6.

  We solved for given temperatures T_i and pressure p_i at each grid-point Eq. (42) for the porosity over time interval Δt to determine an updated value of φ_i at t_k. The corresponding non-linear equations were solved iteratively with an accuracy of better than 10^-12.

• 7.

  We calculated the heat-conductivity from the updated porosity and pressure.

• 8.

  We calculated the heat capacity for given T_i.

• 9.

  The surface temperature T_s was determined from Eq. (32). This equation was solved for T_s, using on the l.h.s. the values for T_i and K from the last iteration step.

• 10.

  Updated values of temperatures T_i at t_k were calculated for all i from the difference equations resulting from the heat-conduction equation.

• 11.

  We checked, if deviation of new values of T_i and φ_i at t_k from current values was sufficiently small.

• 12.

  If this was not the case, we repeated the calculation from step 3 on.

• 13.

  If it was the case, we determined a new stepwidth Δt (see below), advanced k by one, and repeated from step 1 on for the next time step.

This kind of iteration converges usually within about ten to twenty steps. The accuracy requirement was that the relative change of T_i, φ_i between subsequent iteration steps was less than 10^-8. This simple scheme works efficiently because the coefficients in the heat-conduction equation do not strongly depend on temperature. Then this equation can be solved by the simple fixed-point iteration described. Test calculations made with a complete linearisation of the non-linear equation showed that this did not significantly improve the efficiency of the solution method in our particular case. The advantage of our method is that it poses less stringent requirements on the existence and continuity of derivatives than the Newton-Raphson method for convergence of the iteration method.

The boundary condition given by Eq. (32) should in principle be integrated in the difference equations for the heat-conduction equation. Numerical experience showed that this occasionally resulted in stability problems. Our present method is to solve Eq. (32) as a separate equation, using at the current iteration step a value for the conductive heat flux at the surface that was calculated from the quantities of the last iteration step. The resulting value of T_s is prescribed as Dirichlet boundary condition, Eq. (33), for Eq. (23). The temperatures T_s calculated this way at each iteration step converge to the solution of Eq. (23) subject to Neumann boundary condition Eq. (32). This method worked without problems.

The time steps Δt were chosen such that the relative change of the variables over Δt was about 3%. This is sufficiently small that an additional reduction of the stepwidth did not significantly improve the accuracy of the numerical solution; a reduction of the admitted stepwidth by a factor of two resulted in our case in relative changes of the numerical values of the variables by a few times of 10^-4. If the number of iteration steps becomes too big (e.g.  > 20) with this choice, the stepwidth Δt is reduced by a factor of two until the number of iteration steps does no longer exceed the limit. Because we used an implicit solution method, there is no limitation for the stepwidth from stability requirements.

The initial model for the start-up of the solution method assumes a fixed temperature ( = T_s at initial time) within the body. An appropriate set of masses ΔM_i was chosen which resulted (i) in a sufficiently fine grid at the surface to resolve the rapid temperature variations at the surface, and which (ii) is sufficiently fine for calculating the derivative ∂T/∂M_r at the centre with sufficient accuracy. For the initial model the porosities φ_i and radii r_i for the set of mass-shells were calculated from hydrostatic equilibrium and the equation of state for cold pressing, as described in Sect. 4.2.

If a fixed temperature T_b is to be prescribed at the outer boundary, this is technically achieved within the frame of our solution algorithm by letting in Eq. (32).

The solution method also allows us to consider growing bodies by increasing the mass of the outermost shell according to some prescribed growth-law and splitting this shell into two shells at each instant where its mass exceeds some threshold value. This option in our code was not used in the model calculations presented in this paper.

7.2. Some sample calculations

7.2.1. The model of Miyamoto et al.

As a first test we calculated a model with the code using, the same model parameters as in Miyamoto et al. (1981). The basic parameters of the model are given in Table 5 in the column marked with MFT81. The model of Miyamoto et al. (1981) assumes a completely homogeneous body and does not consider the effects of porosity and the possibility of sintering. The data assumed for K and ϱ correspond to average properties of L chondrites that, in fact, have low but non-vanishing porosities, scattering around about φ = 10% (Yomogida & Matsui 1983). The true bulk density and heat-conductivity of completely consolidated chondritic material is higher (ϱ_bulk = 3.6 g cm^-3, K = 4.27 W kg^-1 K^-1, see Yomogida & Matsui 1983, their Table 5). The model is calculated by choosing as initial value φ = 0, which guaranties that during the course of the calculation the porosity remains zero. Heating is only by decay of ²⁶Al. The result for the temperature evolution in the centre of the body and a number of selected radii is shown in Fig. 7a. This is almost identical to the result obtained by Miyamoto et al. (1981) from the analytical solution of the heat-conduction equation, i.e., our code reproduces this exact analytical result.

7.2.2. Model of a porous body

In Fig. 7b we show the results for the temperature evolution of a body having the same size and using a similar set of parameters, but now considering that the heat-conductivity of the porous material, Eq. (30), is different from the value of heat-conductivity used by Miyamoto et al. (1981), and that the material from which the body forms is initially porous and consolidates by sintering. The parameters of the model are given in Table 5 in the column marked with PL0.

We assumed that the porosity of the surface layers at low pressures is φ_srf = 0.6, corresponding to the degree of compaction found in Weidling et al. (2009) for powder material that was subject to numerous impacts. This is what one expects for the early formation-time of asteroids where they grow by repeated slow impacts of much smaller bodies. In deeper layers of the body, where pressures are high because of self-gravity, the material is compressed by isostatic pressing to higher densities up to a limiting value of φ ≈ 0.4, see Sect. 4.2. The corresponding initial distribution of porosities in the interior was calculated as described in Sect. 4.2. For typical results see Fig. 2. This kind of compaction is a purely mechanical effect caused by mutual rolling and gliding of the powder particles driven by an applied pressure, which requires no elevated temperatures and acts therefore already in cold bodies (cold pressing). Starting from this initial configuration, we calculated the evolution of the model. The porosity dependence of the heat-conduction was taken into account by means of Eq. (30). The surface temperature was taken to be constant over time and equal to T_s = 150 K.

Figure 7b shows that the peak temperature achieved in the centre of the body is lower in this model than in the model of Miyamoto et al. (1981). This is because of the higher heat-conductivity in our model after complete sintering (K_b = 4 vs. K_b = 1). In contrast to this, the temperature in the outer layers of the model increases more rapidly than in the model of Miyamoto et al. (1981) because the high porosity outer layer acts as an insulating shield that prevents an efficient loss of heat to outer space. At a temperature of about 700−750 K (depending on pressure), sintering by hot pressing becomes active and the porosity rapidly approaches φ = 0 at higher temperatures. The temperature structure then becomes nearly isothermal in the interior of the body and the temperature drops to the outer boundary value within a layer of only a few km thickness. This is shown in detail in Fig. 8.

Figure 9 shows the evolution of the porosity during the earliest evolutionary phase. An outer shielding powder layer persists during the whole evolution of the body because cooling of the outer layers prevents the material to warm up to the threshold temperature at about 750 K for compaction by hot pressing at low pressures. This behaviour is completely analogous to what was found in Yomogida & Matsui (1984) and can be compared with results by Sahijpal et al. (2007) and Gupta & Sahijpal (2010). They considered gradual sintering in the temperature range of 670–700 K by a smooth interpolation recipe, reducing the porosity from 55% to 0% by compaction of the body, and took into account the respective changes in thermal diffusivity.

Because the porosity approaches zero everywhere in the body where the temperature exceeds the threshold for sintering by hot pressing, the radius of the body shrinks significantly, typically by 20% of its initial value. This occurs after about 0.6 decay timescales of the main heat source, ²⁶Al, in this model, cf. Figs. 8 and 9. The size of the model, marked as PL0 in Table 5, of 85 km corresponds to the final radius that the body would have after compaction to zero porosity (the initial radius before sintering is  ≈ 105 km). The final radius of the body almost exactly equals the hypothetical final radius at zero porosity, because the powder layer remaining at the surface is thin. Moreover, the temperatures shown in Fig. 7b for a number of depth points below the surface correspond to that Lagrangean M_r-coordinate, which after compression to zero porosity would have the given value of the radius coordinate. Before the body shrinks, these points have somewhat bigger depths below the surface.

In Fig. 7c we show the temperature evolution in a model with the same set of parameters as the previous model (see Table 5, model PL), which considers heating by decay of ⁶⁰Fe and long-lived radioactive nuclei as additional heat sources, using an ⁶⁰Fe/⁵⁶Fe ratio at the upper limit of observed values (see Sect. 5.3). The peak central temperature is about 30% higher than without this heat source.

7.2.3. Mass-shells and time steps

Since the mass of the body is constant in the present models, a fixed grid of mass-shells was used in the calculations. This grid was constructed as follows: Starting from an outer layer with some (small) thickness, the thickness of the layers from the outside to the interior increases by a constant factor from one shell to next such that for some given number K of mass-shells the radius of the innermost sphere is 100 m; this fixes the width of the outermost layer. The models of this paper were calculated with K = 300, in which case the outermost layer has a thickness of  ≈ 3 m. This choice of grid guarantees a sufficiently high resolution of the thin outer region of the rapid drop of temperature. An increase of K to 600 did not result in significant changes of the model characteristics; the relative change of central temperature, e.g., is  ≈ 10^-4.

In the model calculations used for fitting models to empirical cooling histories of meteorites described in Sect. 8, the total number K of shells is increased to K = 1200. This was necessary if one requires that even in the region of the steepest temperature decrease towards the surface the relative changes of temperature at some fixed mass-coordinate are at most of the order of a few times 10^-4 if the number of grid points is doubled.

The timesteps Δt were chosen according to the strategy described in Sect. 7.1. The timesteps during the initial heat-up phase were typically a few thousand years. Once sintering commenced, the step size decreased to about 100 years and varied around this value until complete sintering of the body (except for the outermost layers). Then Δt increased continuously and during the final phase was limited to a maximum value of 1 Ma to obtain steps not too widely spaced for plotting purposes. The total number of timesteps required for a complete model run for an evolution period of 100 Ma is between 3500 and 4000.

7.2.4. Dependence on parameters

To assess how the model depends on some important parameters, we show in Fig. 10 results for the same kind of model as model PL, but calculated with two different values of the surface temperature T_s and bulk heat-conductivity coefficient K_b. The models in the left panel were calculated with a fixed surface temperature of 150 K, the models in the right panel with a fixed surface temperature of 250 K. These two values encompass the possible values of the surface temperature for bodies in the region of the asteroid belt for both cases, if either the body is embedded in the accretion disk (cf. Fig. 2) or the accretion disk disappeared and the body is illuminated by the radiation of the young sun. There are only small differences between the run of the temperature at different depths, i.e., the temperature evolution below the immediate surface layer does not critically depend on the surface temperature, at least not for bodies with radii of the order of 100 km, which are our main topic. Therefore we did not attempt to calculate T_s as precise as possible from Eq. (32) and simply assumed a reasonable but fixed value for all times.

The upper panel in Fig. 10 was calculated with a value for the heat-conduction coefficient of K_b = 4 W m^-1 K^-1, the lower panel with K_b = 2 W m^-1 K^-1. The first value corresponds to what has been found as the average value for H and L chondrites if the measured values are extrapolated to zero porosity, see Sect. 5.4. As one can see from Fig. 6, there is significant scatter in the heat-conduction coefficient (of presently unknown origin) and it is not clear whether the investigated sample of H and L chondrites are representative for the whole material of the parent body of the H or L chondrites. The value of K_b = 4 W m^-1 K^-1 corresponds to typical values of K for pure silicate minerals (cf. Yomogida & Matsui 1983) and therefore probably represents the upper limit for the possible values of K_b. Lower values may therefore also be of interest for real planetesimals. Figure 10 shows that the results of the model calculation depend significantly on the value of K_b. Because it is presently not possible to determine K_b from first principles for chondritic material, we considered K_b in our later model calculations to be a free parameter (but, of course, restricted to the range of values found for chondrites).

7.3. Maximum central temperatures

The maximum central temperature that is reached during the course of the evolution of a planetesimal is an indicator of what kind of changes the material may undergo. If the central temperatures exceed the solidus temperature of chondritic material of about 1400 K (Agee et al. 1995), partial melting occurs and the body starts to differentiate. If the temperature stays below the threshold temperature of about 700 K (at  ≈ 0.1 bar) for sintering, the whole body retains its porous structure. The maximum central temperature T_c,max depends mainly

• 1.

  on the radius, R, of the body, because this determines how efficiently heat can be removed from the central region by heat-conduction; and

• 2.

  on the formation-time, t_form, because this determines how much of the initial inventory of radioactive material already decayed before the body grew to its final size.

Figure 11 shows the dependence of T_c,max on R and t_form for models of initially porous bodies and, for comparison, of bodies with pore-free material. Obviously, the thermal evolution history of initially porous bodies is very different from history of equal-sized compact bodies. Models are shown for three different assumptions on the abundance of ⁶⁰Fe as additional heat source besides ²⁶Al.

For porous bodies smaller than  ≈ 5 km radius the initial porosity is very high because they are not even compacted by cold pressing (cf. Fig. 2). Because of their low initial heat-conductivity even small bodies (R ≳ 0.5 km) heat up at least to the threshold temperature for sintering and become compacted in their interior, because the strongly insulating powder layer on their surface prevents their cooling. Completely compact bodies would reach central temperatures higher than 700 K only for radii exceeding  ≈ 5 km because of the much more efficient heat-conduction.

For initially porous bodies bigger than  ≈50 km radius the initial porosity is already low throughout almost all of the body because these bodies are already strongly compacted by cold pressing (cf. Fig. 2) and the remaining porosity rapidly disappears by sintering. Their thermal evolution is essentially identical to that of completely compact bodies, except, of course, in the layers close to the surface that retain part of their initial porosity.

Porous bodies with radii between  ≈5 km and  ≈20 km are already significantly compacted in their central part by cold pressing (cf. Fig. 2) but have initially an extended low-porosity outer mantle. Porous bodies with radii between  ≈20 km and  ≈50 km also are already compacted throughout the body by cold pressing (cf. Fig. 2), except for the outer  ≈10% of their radius. They show the most complex dependence of T_c,max on R and formation-time.

Temperatures above T_c,max1500 K were not considered because here we considered parent bodies of undifferentiated meteorites. At a temperature of T_c,max ≳ 1400 K the silicates start melting (Agee et al. 1995) and differentiation is unavoidable.

8. Thermal history of the H chondrite parent body

8.1. Empirical data for cooling history

Most meteorites reach the Earth as cm- or m-sized rocks, because they are the result of repeated impact fragmentation of the initially much larger parent bodies. Because these events can be highly energetic, they change the chemistry and the structure, and also disturb the isotopic memory (i.e., the age information). Hence, information of cooling histories extending back to the origin of the solar system must be restricted to meteorites that

• 1.

  show extraordinarily low mineralogical and structuralcharacteristics of shock metamorphism induced by asteroidcollisions;

• 2.

  were dated with high precision; and

• 3.

  were dated simultaneously by a set of various high- and low-temperature chronometers tracing the cooling history from high temperatures (>1000 K, e.g., Hf-W) down to intermediate (e.g. U-Pb-Pb or K-Ar ) and very low temperatures (if possible  <400 K, e.g. ²⁴⁴Pu fission tracks).

Such high quality data on un-shocked chondrites are quite limited, in spite of the high number of meteorites available in terrestrial collections. While in the case of L chondrites, a major impact 470 Ma ago (Korochantseva et al. 2007) seems to have deeply erased the primordial low temperature cooling history, H chondrites seem more promising. For a number of seven – noticeably un-shocked – H chondrites, complete high precision Hf-W, U-Pb-Pb, Ar-Ar and ²⁴⁴Pu fission track data along with metallographic cooling rate data exist. Table 6 shows such data for the H6 chondrite Kernouvè and the H5 chondrite Richardton, which we can use for a preliminary sample calculation.

8.2. Fit of selected H chondrites

A fit to the data in Table 6 is shown in Fig. 12. The chronological data for these two meteorites best fit to the cooling curves in an asteroid at 36 and 8.9 km depth. The properties of the parent body in Table 7 were obtained according to the following fit-procedure: the initial abundance of ²⁶Al is roughly constrained by the ²⁶Al/²⁷Al ratio of ordinary chondrite chondrules, which is an upper limit shortly before parent body formation. Chondrule measurements indicate an ²⁶Al abundance corresponding roughly to 2−3 Ma formation-time after CAIs. Similarly, the ⁶⁰Fe abundance is constrained by primitive type 3 ordinary chondrites (sulphides, see Tachibana & Huss 2003). Furthermore, abundances of ²⁶Al and ⁶⁰Fe (or, in other words, the formation-time of the parent body) must be such that sufficiently high maximum metamorphic temperatures in the asteroid are achieved to allow strong type 6 metamorphism, but also not too high to prevent partial melting, metal silicate separation, and differentiation. Heat-conductivity and radius mainly influence the total duration of the cooling time of the parent body. We arbitrarily chose 100 km, although a smaller asteroid would also allow extended cooling as observed for Kernouvè, which in this model needs only a burial depth of about 36 km. The boundaries separating H6 from H5, H4 and H3 material are relatively shallow in the model, because of the insulating porous thin outer layer.

8.3. Discussion

The most prominent feature in our new model is the possibility of relatively small parent planetesimals with significant heat retention. In the H chondrite parent body model, this shows up in relatively thin layers of less heated or metamorphosed material. Moreover, the relatively fast cooling required to achieve temperatures below Hf-W and U-Pb-Pb closure for the H5 chondrite Richardton (within 3 and 13 Ma, respectively) sets an upper limit to the contribution of decay heat of ⁶⁰Fe (roughly about 20−30%). For example, if ⁶⁰Fe contributed all decay heat of ordinary chondrite parent body metamorphism, sufficiently fast cooling of H5 metamorphosed material would not have been possible, because the half-live of ⁶⁰Fe is 2.6 Ma (a new revised value) versus 0.72 Ma for ²⁶Al, implying significantly longer heat production and implicit slower cooling in the first few Ma. This result is in line with the initial ⁶⁰Fe concentration found in primitive type 3 ordinary chondrites (Tachibana & Huss 2003), lower than initial values previously obtained for CAIs (Birck & Lugmair 1988), and supports the view that ⁶⁰Fe was likely not distributed homogeneously in the early solar system. A more detailed H chondrite parent body modelling will be presented elsewhere.

9. Concluding remarks

We constructed a model for the thermal evolution of the parent body of H chondrites. The model considers compaction by cold pressing and sintering by “hot pressing”. The heat-conductivity of the porous material was determined by combining new data obtained by Krause et al. (2011a) for high-porosity material with data for porous chondrites (Yomogida & Matsui 1983). A model was fitted to data on the cooling history for two H chondrites, Kernouvé (H6) and Richardton (H5), for which particularly accurate data are available. We showed that it is possible to find a consistent fit for the parent body size and formation-time that reproduces the empirically determined cooling history of H5 and H6 chondrites with good accuracy.

For obtaining our consistent fit, it was necessary to include (besides radius of the body and formation-time) also the abundance of ⁶⁰Fe into the fitting procedure. A value of ⁶⁰Fe/⁵⁶Fe was determined, which is within the range of values reported in the literature for different meteorites. No other arbitrary fit parameters are required; all other properties of the model are fixed by the physics of the problem.

The new model predicts shallow outer layers for petrologic types 3 to 5 because of the strong blanketing effect of an outer powder layer, which escapes sintering. In this respect the model deviates considerably from previous models that are based on the analytical model of Miyamoto et al. (1981). Other properties of the model are similar to older models; in particular radius and formation-time are not substantially different.

The present model, though relaxing some earlier shortcomings, still has a number of shortcomings. The most stringent is the instantaneous formation hypothesis, which needs to be relaxed because the formation-time of the body by run-away growth (of the order of 10⁵a) is shorter than the decay time of ²⁶Al, the main heating source, but probably not short enough to be completely negligible. The second severe shortcoming is that heat-conductivity of porous media cannot yet be treated from first principles on. This is not possible with the presently available computing resources, but this may change in the near future. Other shortcomings are a simplistic treatment of sintering and of the outer boundary temperature. We modelled sintering with the same kind of theoretical description as Yomogida & Matsui (1984) to facilitate a comparison of our results with that model. This treatment of sintering should be replaced by more elaborated modern theories of hot pressing.

---

Acknowledgments

We thank the referee S. Sahijpal for a constructive and helpful referee report. This work was supported in part by “Forschergruppe 759” and in part by “Schwerpunktprogramm 1385”. Both are supported by the “Deutsche Forschungsgemeinschaft (DFG)”.

References

All Tables

All Figures