Issue 
A&A
Volume 537, January 2012



Article Number  A45  
Number of page(s)  19  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201117177  
Published online  06 January 2012 
Thermal evolution and sintering of chondritic planetesimals
^{1} Institut für Theoretische Astrophysik, Zentrum für Astronomie, Universität Heidelberg, AlbertÜberleStr. 2, 69120 Heidelberg, Germany
email: gail@ita.uniheidelberg.de
^{2} Institut für Geowissenschaften, Universität Heidelberg, Im Neuenheimer Feld 236, 69120 Heidelberg, Germany
^{3} Institut für Planetologie, Universität Münster, WilhelmKlemmStr. 10, 48149 Münster, Germany
Received: 2 May 2011
Accepted: 21 October 2011
Aims. Radiometric ages for chondritic meteorites and their components provide information on the accretion timescale of chondrite parent bodies, and on cooling paths within certain areas of these bodies. However, to use this age information for constraining 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 heatproducing shortlived radionuclides (^{26}Al and ^{60}Fe), which are determined by the accretion timescale and the terminal size, chemical composition and physical properties of the chondritic planetesimals. The major aim of this study is to assess the effects of sintering of initially porous material on the thermal evolution of planetesimals, and to constrain the values of basic parameters that determined the structure and evolution of the H chondrite parent body.
Methods. We present a new code for modelling the thermal evolution of ordinary chondrite parent bodies that initially are highly porous and undergo sintering by hot pressing as they are heated by decay of radioactive nuclei. The pressure and temperature stratification in the interior of the bodies was calculated by solving the equations of hydrostatic equilibrium and energy transport. The decrease of porosity of the granular material by hot pressing due to selfgravity was followed by solving a set of equations for the sintering of powder materials. For the heatconductivity of granular material we combined recently measured data for highly porous powder materials, relevant for the surface layers of planetesimals, with data for heatconductivity of chondrite material, relevant for the strongly sintered material in deeper layers.
Results. Our new model demonstrates that in initially porous planetesimals heating to central temperatures sufficient for melting can occur for bodies a few km in size, that is, a factor of ≈10 smaller than for compact bodies. Furthermore, for high initial ^{60}Fe abundances small bodies may differentiate even when they had formed as late as 3−4 Ma after CAI formation. To demonstrate the capability of our new model, the thermal evolution of the H chondrite parent body was reconstructed. The model starts with a porous body that is later compacted first by “cold pressing” at low temperatures and then by “hot pressing” for temperatures above ≈700 K, i.e., the threshold temperature for sintering of silicates. The thermal model was fitted to the wellconstrained cooling histories of the two H chondrites Kernouvé (H6) and Richardton (H5). The best fit was obtained for a parent body with a radius of 100 km that accreted at t = 2.3 Ma after CAI formation, and had an initial ^{60}Fe/^{56}Fe = 4.1 × 10^{7}. Burial depths of 8.3 km and 36 km for Richardton and Kernouvé were able to reproduce their empirically determined cooling history. These burial depths are shallower than those derived in previous models. This reflects the strong insulating effect of the residual powder surface layer, which is characterised by a low thermal conductivity.
Key words: minor planets, asteroids: general / meteorites, meteors, meteoroids / planets and satellites: physical evolution
© 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. ^{26}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 icefree, 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 heatconduction equation with constant coefficients for a spherical body heated by a homogeneously distributed and exponentially decaying heat source (^{26}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 heatconduction, 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 shortlived 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 selfgravity 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 heatconductivity 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 heatconductivity 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 ^{26}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 shortlived nuclides like ^{26}Al or ^{60}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 ^{26}Al abundance, i.e., mainly formationtime. Planetesimals formed at ≈2 Ma after calciumaluminiumrich 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 HfW dating method (e.g. Kleine et al. 2005, 2008) is useful for cooling below 1050 to 1150 K. The closure temperature for the UPbPb system in phosphates is ≈720 K (e.g. Göpel et al. 1994; Bouvier et al. 2007), ArAr ages reflect the cooling below ≈550 K for oligoclase feldspar. The annealing temperature for ^{244}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 FeNi 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.
Basic mineral species considered for the calculation of properties of chondrite material, their atomic weight A, massdensity ϱ, and abbreviation of mineral name.
The pure mineral components used in our model for the composition of planetesimals and their massdensities 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_{80}Fa_{20} 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 massfractions 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 nonideal 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 massfractions 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 finegrained 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 submicrometer sized, most of them having diameters smaller than 0.5 μm and many being nanometer 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.
Typical mineral composition of chondrites, massdensities ϱ of components, massfractions X_{min} of minerals, and massfractions X_{el} of the elements that release heat by radioactive decays (data for X_{min} from Yomogida & Matsui 1983).
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 equalsized 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_{0} = 7.3 and D_{0} = 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 equalsized 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 precompaction 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.
Some properties of chondrules and matrix in H and L chondrites (data from Scott 2007).
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 chainlike structures. Such highporosity 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.
Fig. 1
Variation of a) midplane temperature and b) pressure with distance from the star at different evolutionary epochs (0.2 Ma and from 0.5 Ma to 3 Ma in steps of 0.5 Ma) in a model for the evolution of the solar nebula (onezone αmodel). They define the outer boundary conditions for a planetesimal embedded in an accretion disk. In the left part the regions are indicated where the disappearance of an important absorber results in an extended region of nearly constant midplane temperature. 
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 inwarddirected 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 buildup 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 buildup 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 wellknown hydrostatic pressure equation (7)where (8)The density ϱ is the massdensity 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_{0} 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.
Fig. 2
Relation between relative density (filling factor) D and maximum experienced pressure for isostatic pressing (top) and run of relative density within planetesimals of the indicated size (bottom) that results from cold pressing caused by selfgravity. 
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^{1}(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 equalsized 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 selfgravity of the body. For planetesimals bigger than 10 km most part of the body is already compacted by cold isostatic pressing caused by selfgravity to about a density close to that of the densest random packing of equalsized 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. Heatconduction 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 pdVterm 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 speciesindependent 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_{r}F_{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 heatconduction 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 heatconduction by the gas in the pores. If the material is translucent, one has also to consider a contribution to heatconduction 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 heatconduction coefficient K that is different for each of the different transport processes.
The essential material properties that enter into Eq. (23) are the specific heatperunit mass, c_{v}, the heatconductivity 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 massfraction of the ith component in the mixture of solid components and c_{v,i} the heatcapacityperunit 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 kbarrange (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 heatcapacityperunit 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 Curietemperature 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.
Fig. 3
Specific heat of planetesimal material for mineral mixtures of H chondrites and L chondrites. The dotted line is the analytical approximation of Yomogida & Matsui (1984). 
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 ^{26}Al and possibly ^{60}Fe (cf. the discussion in Ghosh et al. 2006). More longlived radioactive nuclei, essentially ^{232}Th, ^{235,238}U, and ^{40}K, are responsible for the later longterm evolution of the temperature. For the nuclei decaying by βdecay (^{26}Al, ^{60}Fe, ^{40}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 heatproducing 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 massfraction 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 timescale for efold decrease of abundance of the isotope, and t is the time elapsed since formation of the solar system.
Data for calculating heating rates by decay of radioactive nuclei.
Fig. 4
Contributions of different heating mechanisms to the total heating rate. Time t is after formation of CAIs. The dashed line indicates the formationtime of the H chondrite parent body (2.3 Ma, see Sect. 8). The release of gravitational energy by contraction is shown for the final model of the H chondrite parent body. 
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 ^{26}Al is that given by Nyquist et al. (2009). The abundance of ^{60}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 ^{60}Fe was not homogeneously distributed in the early solar system (Quitté et al. 2010), which means that the initial ^{60}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 ^{60}Fe the recent revised value for the halflife τ_{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 ^{60}Fe abundance (see Table 4 for this). The dominant heating source is ^{26}Al at the time of formation of planetesimals (t ≲ 5 Ma), but ^{60}Fe dominates as heat source for an extended period from ≈ 5 Ma to ≈ 20 Ma because of the revised ^{60}Fe halflife.
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. Heatconduction by the porous solid material
For the heatconductivity 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 gasfillings in the pores.
Fig. 5
Variation of heat conductivity K with porosity φ. The results for finegrained silica powder (filled circles) from experiments of Krause et al. (2011a,b), and for particulate basalt (crosses) from Fountain & West (1970). The typical grain size is indicated for both cases. Open squares and open circles are experimental results for heatconductivity and porosity for ordinary H and L chondrites, respectively, from (Yomogida & Matsui 1983). The solid line is an analytical fit to the data, according to Eq. (30). 
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 φ_{1}. 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 φ_{1} = 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_{1}(φ).
At high porosities Fig. 5 shows the conductivity K for a silica powder that consists of equalsized 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 welldefined temperature but the heatconductivity 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 φ_{1} 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 φ_{2} = 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 heatconduction coefficients. Second, the heatconduction 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 heatconduction 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 heatconduction 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 heatconductivity 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 powerlaw 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 temperaturedependent. Figure 6 shows data for K for H and L chondrites as given by Yomogida & Matsui (1983), divided by K_{1}(φ), given by Eq. (28). The data for K(T)/K_{1}(φ) 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 meteoritetometeorite 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 temperatureindependent.
5.5. Heatconduction by other processes
Fig. 6
Temperature variation of the thermal conductivity K. Plotted are data for H chondrites (open rectangles) and L chondrites (open circles), normalised with the φvariation according to approximation (28). The lines connect data for each of the meteorites. 
The material of the planetesimals is dominated by minerals that are transparent in the farinfrared 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 nonnegligible. 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 heatconduction 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 heatconduction never becomes an important energy transport mode. This is because if temperature becomes high enough for radiative heatconduction to contribute somewhat to the net heat flux, the onset of sintering above ≈ 700 K results in a strong increase of heatconductivity by phonons that outnumbers radiative contributions again.
A possible contribution to heatconduction 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 heatconduction Eq. (23) one has to specify an initial condition T(r,t_{0}) at some initial instant t_{0} 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 lefthand side describes the energy flow from the interior towards the surface by heatconduction. 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 protosun, 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 kmsized 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^{5} years until runaway growth commences after the biggest planetesimals reach sizes of about 10 − 20 km (e.g. Weidenschilling & Cuzzi 2006; Nagasawa et al. 2007). During runaway growth the mass rapidly increases within less than 10^{5} 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 ^{26}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 ^{26}Al decay of those planetesimals that contributed to the growth of a 100 kmsized 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 formationtime 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 twostep process. The initially very loosely packed dust material in the planetesimals comes under increasing pressure by the growing selfgravity 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^{2} 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 equalsized spheres with initial radius R_{0}. 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 hotpressing 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 cutaway 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_{0}. 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_{0} is the relative density of the initial packing of spheres with radius R_{0}. For a given number of contact points Z and given D_{0}, R_{0}, 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 powerlaw creep (38)and that is given in terms of the rate of change of relative density as (39)The stress σ_{1} 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 σ_{1} 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^{2}, and average crosssection of the cell occupied by one granular unit, C_{av}, as (40)In the hexagonal prismatic packing model favoured by Yomogida & Matsui (1984), the crosssection 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 σ_{1} 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 poreclosing 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 prefactor 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 prefactor 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 heatconduction, 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 heatconductivity, Eq. (30), and heat capacity, Eq. (26), and together with appropriate initial and boundary conditions.
The heatconduction equation and the pressure equation are rewritten in terms of M_{r}, defined by Eq. (8), as independent variable and are discretised for a set of fixed massshells 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 heatconduction 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 firstorder ordinary differential equation for φ, Eq. (44), is solved by the fully implicit Euler method.
To account for the nonlinear coupling between the different equations, we performed a fixedpoint 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 gridpoint Eq. (42) for the porosity over time interval Δt to determine an updated value of φ_{i} at t_{k}. The corresponding nonlinear equations were solved iteratively with an accuracy of better than 10^{12}.

7.
We calculated the heatconductivity 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 heatconduction 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.
The boundary condition given by Eq. (32) should in principle be integrated in the difference equations for the heatconduction 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 startup 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 massshells 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 growthlaw 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.
Fig. 7
Temperature evolution of a body of 85 km radius at different depths from the surface and at centre: dotted line at depth 4.25 km, short dashed line at 17 km depth, long dashed line at 23 km depth. The full line shows the temperature evolution at the centre. a) The model is calculated with the same physical input as in the analytical model of Miyamoto et al. (1981). We provide the model parameters of MFT81 in Table 5. b) Similar model, but now calculated for a porous body, considering thermal conductivity of porous material according to Eq. (30), and sintering and cold pressing as described in this paper. We provide the model parameters of PL0 in Table 5. c) Same kind of model as b), but now additional heating by ^{60}Fe and longlived nuclei considered. We provide the model parameters of PL in Table 5. 
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 nonvanishing porosities, scattering around about φ = 10% (Yomogida & Matsui 1983). The true bulk density and heatconductivity 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 ^{26}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 heatconduction 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 heatconductivity of the porous material, Eq. (30), is different from the value of heatconductivity 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 formationtime 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 selfgravity, 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 heatconduction 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.
Model parameters for the model of Miyamoto et al. (1981), (MFT81) for a consolidated body (average L chondrite), and for a similar model of an initially porous body without (PL0) and with (PL) additional heating by decay of ^{60}Fe and longlived nuclei.
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 heatconductivity 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.
Fig. 8
Evolution of the radial distribution of the temperature for model PL (see Table 5 for its definition). Note the radial shrinking of the body by compaction of the initially porous material at about 0.6 decay timescales of ^{26}Al. 
Fig. 9
Evolution of the radial distribution of the filling factor D (relative density) for model PL (see Table 5 for its definition). Shown is the very initial phase of the evolution where the initially porous material is compacted by sintering. The resulting shrinking of the planetesimal size occurs at about 0.6 decay timescales of ^{26}Al. 
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, ^{26}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 ^{60}Fe and longlived radioactive nuclei as additional heat sources, using an ^{60}Fe/^{56}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. Massshells and time steps
Since the mass of the body is constant in the present models, a fixed grid of massshells 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 massshells 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 masscoordinate 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 heatup 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
Fig. 10
Temperature evolution of test models for porous bodies with modified values of parameters. Left panel with boundary temperature T_{b} = 150 K, right panel with T_{b} = 250 K. Upper panel with heatconduction coefficient K_{b} = 4 W m^{1} K^{1}, lower panel with K_{b} = 2 W m^{1} K^{1}. 
Fig. 11
Variation of maximum temperature in the centre of a planetesimal with radius R and instant of formation t_{form}. Upper panel: completely consolidated bodies with porosity φ = 0. Lower panel: models of initially porous bodies with porosity φ_{srf} = 0.5 that consolidate in the interior during the course of their thermal evolution. The lines correspond to the indicated maximum central temperature. At temperatures above ≈ 1400 K partial melting of the mineral mixture is expected. Temperatures exceeding 1500 K are excluded for this reason. Shown are models for three different initial ^{60}Fe/^{56}Fe abundance ratios: Left panel: without ^{60}Fe. Middle panel: the optimumfit value of 4.1 × 10^{7} for the H chondrite parent body, see Sect. 8. Right panel: highest value of 1.6 × 10^{6} given in Table 5. 
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 heatconductivity 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 heatconduction 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 heatconduction 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 heatconduction; and

2.
on the formationtime, 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 porefree material. Obviously, the thermal evolution history of initially porous bodies is very different from history of equalsized compact bodies. Models are shown for three different assumptions on the abundance of ^{60}Fe as additional heat source besides ^{26}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 heatconductivity 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 heatconduction.
Fig. 12
Optimumfit model for the cooling history of the parent body of H chondrites. The abscissa is the time elapsed since the formation of the CAIs. We show the evolution of temperature at a number of depths below the surface. The upper contours of shaded areas correspond (from bottom to top) to depths of 0.32 km, 2.3 km, 7.8 km, 11 km, and the highest contour to the centre. The rectangular boxes and circles correspond to the empirical data for H6 and H5 chondrites, respectively, given in Table 6. Crosses are error bars. The dashed lines correspond to the temperature evolution at depths of 8.9 km (lower line) and 36 km (upper line) below the surface. They represent our best fit to the empirical data. 
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 lowporosity 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 formationtime.
Temperatures above T_{c,max}1500 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 msized 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 lowtemperature chronometers tracing the cooling history from high temperatures (>1000 K, e.g., HfW) down to intermediate (e.g. UPbPb or KAr ) and very low temperatures (if possible <400 K, e.g. ^{244}Pu fission tracks).
Such high quality data on unshocked 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 unshocked – H chondrites, complete high precision HfW, UPbPb, ArAr and ^{244}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.
Key data for cooling history of selected H chondrites.
Properties of the optimumfit model.
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 fitprocedure: the initial abundance of ^{26}Al is roughly constrained by the ^{26}Al/^{27}Al ratio of ordinary chondrite chondrules, which is an upper limit shortly before parent body formation. Chondrule measurements indicate an ^{26}Al abundance corresponding roughly to 2−3 Ma formationtime after CAIs. Similarly, the ^{60}Fe abundance is constrained by primitive type 3 ordinary chondrites (sulphides, see Tachibana & Huss 2003). Furthermore, abundances of ^{26}Al and ^{60}Fe (or, in other words, the formationtime 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. Heatconductivity 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 HfW and UPbPb closure for the H5 chondrite Richardton (within 3 and 13 Ma, respectively) sets an upper limit to the contribution of decay heat of ^{60}Fe (roughly about 20−30%). For example, if ^{60}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 halflive of ^{60}Fe is 2.6 Ma (a new revised value) versus 0.72 Ma for ^{26}Al, implying significantly longer heat production and implicit slower cooling in the first few Ma. This result is in line with the initial ^{60}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 ^{60}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 heatconductivity of the porous material was determined by combining new data obtained by Krause et al. (2011a) for highporosity 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 formationtime 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 formationtime) also the abundance of ^{60}Fe into the fitting procedure. A value of ^{60}Fe/^{56}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 formationtime 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 formationtime of the body by runaway growth (of the order of 10^{5}a) is shorter than the decay time of ^{26}Al, the main heating source, but probably not short enough to be completely negligible. The second severe shortcoming is that heatconductivity 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
 Agee, C. B., Li, J., Shannon, M. C., & Circone, S. 1995, J. Geophys. Res., 100, 17725 [NASA ADS] [CrossRef] [Google Scholar]
 Akridge, G., Benoit, P. H., & Sears, D. W. G. 1997, Lunar Plan. Sc. Conf., XXVIII, 1178 [Google Scholar]
 Akridge, G., Benoit, P. H., & Sears, D. W. G. 1998, Icarus, 132, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
 Arzt, E. 1982, Acta metall., 30, 1883 [CrossRef] [Google Scholar]
 Arzt, E., Ashby, M. F., & Easterling, K. E. 1983, Metallurgical Transact. A, 14A, 211 [Google Scholar]
 Barin, I. 1995, Thermochemical Data of Pure Substances, 3rd ed., Vol. I, II (VCH Verlagsgesellschaft Weinheim) [Google Scholar]
 Birck, J. L., & Lugmair, G. W. 1988, Earth & Plan. Sci. Lett., 90, 131 [Google Scholar]
 Blum, J., Schräpler, R., Davidsson, B. J. R., & TrigoRodríguez, J. M. 2006, ApJ, 652, 1768 [NASA ADS] [CrossRef] [Google Scholar]
 Bouvier, A., BlichertToft, J., Moynier, F., Vervoort, J., & Albarède, F. 2007, Geochim. Cosmochim. Acta, 71, 1583 [NASA ADS] [CrossRef] [Google Scholar]
 Bradley, J. 2010, in Astromineralogy, 2nd ed., ed. T. K. Henning (Berlin: Springer), Lecture Notes in Physics, 609, 259 [Google Scholar]
 Finocchi, F., & Gail, H.P. 1997, A&A, 327, 825 [NASA ADS] [Google Scholar]
 Fischmeister, H. F., & Arzt, E. 1983, Powder Metallurgy, 26, 82 [Google Scholar]
 Fountain, J. A., & West, E. A. 1970, J. Geophys. Res., 75, 4063 [NASA ADS] [CrossRef] [Google Scholar]
 Ghosh, A., & McSween, H. Y. 1999, Meteoritics & Plan. Sci., 34, 121 [Google Scholar]
 Ghosh, A., Weidenschilling, S. J., & McSween, H. Y. 2003, Meteoritics & Plan. Sci., 38, 711 [Google Scholar]
 Ghosh, A., Weidenschilling, S. J., McSween, H. Y., & Rubin, A. 2006, in Meteorites and the Early Solar System II, ed. D. S. Lauretta, & H. Y. McSween Jr. (Tucson: Univ. of Arizona Press), 555 [Google Scholar]
 Göpel, C., Manhes, G., & Allegre, C. J. 1994, Earth & Plan. Sci. Lett., 121, 153 [Google Scholar]
 Grimm, R. E., & McSween, H. Y. 1989, Icarus, 82, 244 [NASA ADS] [CrossRef] [Google Scholar]
 Gupta, G., & Sahijpal, S. 2010, J. Geophys. Res., 115, E08001 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Krause, M., Geretshauser, R. J., Speith, R., & Blum, J. 2009, ApJ, 701, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Harrison, K. P., & Grimm, R. E. 2010, Geochim. Cosmochim. Acta, 74, 5410 [NASA ADS] [CrossRef] [Google Scholar]
 Herpfer, M., Larimer, J., & Goldstein, J. 1994, Geochim. Cosmochim. Acta, 58, 1353 [NASA ADS] [CrossRef] [Google Scholar]
 Hevey, P. J., & Sanders, I. S. 2006, Meteoritics & Plan. Sci., 41, 95 [Google Scholar]
 Hopfe, W., & Goldstein, J. 2001, Meteoritics & Plan. Sci., 36, 135 [Google Scholar]
 Jarosewich, E. 1990, Meteoritics, 25, 323 [NASA ADS] [CrossRef] [Google Scholar]
 Jeager, H. M., & Nagel, S. R. 1992, Science, 255, 1523 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kakar, A. K., & Chaklader, A. C. D. 1967, J. Appl. Phys., 38, 3223 [NASA ADS] [CrossRef] [Google Scholar]
 Kleine, T., Mezger, K., Palme, H. E. S., & Münker, C. 2005, Geochim. Cosmochim. Acta, 69, 5805 [NASA ADS] [CrossRef] [Google Scholar]
 Kleine, T., Touboul, M., Van Orman, J., et al. 2008, Earth & Plan. Sci. Lett., 270, 106 [Google Scholar]
 Korochantseva, E. V., Trieloff, M., Lorenz, C. A., et al. 2007, Meteor. Planet. Sci., 42, 113 [Google Scholar]
 Kothe, S., Güttler, C., & Blum, J. 2010, ApJ, 725, 1242 [NASA ADS] [CrossRef] [Google Scholar]
 Krause, M., Blum, J., Skorov, Y., & Trieloff, M. 2011a, Icarus, 214, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Krause, M., Henke, S., Gail, H.P., et al. 2011b, Lunar Planet. Sci. Conf. Lett., 42, 2696 [NASA ADS] [Google Scholar]
 Larsson, P. L., Biwa, S., & Storåkers, B. 1996, Acta mater., 44, 3655 [Google Scholar]
 Lodders, K., Palme, H., & Gail, H. P. 2009, in LandoltBörnstein, New Series, Group IV, Vol. 4, ed. J. E. Trümper (Berlin: Springer), 560 [Google Scholar]
 McSween, H., Ghosh, A., Grimm, R. E., Wilson, L., & Young, E. D. 2003, in Asteroids III, ed. W. F. Bottke (Tucson: Univ. of Arizona Press), 559 [Google Scholar]
 Miyamoto, M., Fujii, N., & Takeda, H. 1981, Lunar Planet. Sci. Conf. Lett., 12B, 1145 [NASA ADS] [Google Scholar]
 Nagasawa, M., Thommes, E. W., Kenyon, S. J., Bromley, B. C., & Lin, D. N. C. 2007, in Protostars and Planets V, ed. B. Reipurt, D. Jewitt, & K. Keil (Tucson: University of Arizona Press), 639 [Google Scholar]
 Nyquist, L. E., Kleine, T., Shih, C.Y., & Reese, Y. D. 2009, Geochim. Cosmochim. Acta, 73, 5115 [NASA ADS] [CrossRef] [Google Scholar]
 Onoda, G. Y., & Liniger, E. G. 1990, Phys. Rev. Lett., 64, 2727 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Pellas, P., Fiéni, C., Trieloff, M., & Jessberger, E. 1997, Geochim. Cosmochim. Acta, 61, 3477 [NASA ADS] [CrossRef] [Google Scholar]
 Quitté, G., Halliday, A. N., Meyer, B. S., et al. 2007, ApJ, 655, 678 [NASA ADS] [CrossRef] [Google Scholar]
 Quitté, G., Markowski, A., Latkoczy, C., Gabriel, A., & Pack, A. 2010, ApJ, 720, 1215 [NASA ADS] [CrossRef] [Google Scholar]
 Rao, A. S., & Chaklader, A. C. D. 1972, J. Am. Ceram. Soc., 55, 596 [CrossRef] [Google Scholar]
 Rietmeijer, F. 1993, Earth & Plan. Sci. Lett., 117, 609 [Google Scholar]
 Rugel, G., Faestermann, T., Knie, K., et al. 2009, Phys. Rev. Lett., 103, 072502 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Sahijpal, S., Soni, P., & Gupta, G. 2007, Meteoritics & Plan. Sci., 42, 1529 [NASA ADS] [CrossRef] [Google Scholar]
 Schwenn, M. B., & Goetze, C. 1978, Tectonophysics, 48, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Scott, G. D. 1962, Nature, 194, 956 [NASA ADS] [CrossRef] [Google Scholar]
 Scott, E. R. D. 2007, Ann. Rev. Earth & Plan. Sci., 35, 577 [Google Scholar]
 Senshu, H. 2004, Lunar Plan. Sc. Conf., XXXV, 1557 [Google Scholar]
 Storåkers, B., Fleck, N. A., & McMeeking, R. M. 1999, J. Mech. Phys. Solids, 47, 785 [NASA ADS] [CrossRef] [Google Scholar]
 Tachibana, S., & Huss, G. R. 2003, ApJ, 588, L41 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Taylor, G. J., Maggiore, P., Scott, E. R. D., & Keil, A. E. 1987, Icarus, 69, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Trieloff, M., Jessberger, E., Herrwerth, I., et al. 2003, Nature, 422, 502 [NASA ADS] [CrossRef] [Google Scholar]
 Van Schmus, W. R. 1995, in Global Earth Physics, A Handbook on Physical Constants, AGU Reference Shelf 1 (American Geophysical Union), 283 [Google Scholar]
 Wehrstedt, M., & Gail, H. 2002, A&A, 385, 181 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wehrstedt, M., & Gail, H. 2008 [arXiv:0804.3377] [Google Scholar]
 Weidenschilling, S. J., & Cuzzi, J. N. 2006, in Meteorites and the Early Solar System II, ed. D. S. Lauretta, & H. Y. McSween (Tucson: University of Arizona Press), 473 [Google Scholar]
 Weidling, R., Güttler, C., Blum, J., & Brauer, F. 2009, ApJ, 696, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Wilkinson, D., & Ashby, M. 1975, Acta Metallurgica, 23, 1277 [CrossRef] [Google Scholar]
 Yang, R. Y., Zou, R. P., & Yu, A. B. 2000, Phys. Rev. E, 62, 3900 [NASA ADS] [CrossRef] [Google Scholar]
 Yomogida, K., & Matsui, T. 1983, J. Geophys. Res., 88, 9513 [Google Scholar]
 Yomogida, K., & Matsui, T. 1984, Earth & Plan. Sci. L., 68, 34 [Google Scholar]
All Tables
Basic mineral species considered for the calculation of properties of chondrite material, their atomic weight A, massdensity ϱ, and abbreviation of mineral name.
Typical mineral composition of chondrites, massdensities ϱ of components, massfractions X_{min} of minerals, and massfractions X_{el} of the elements that release heat by radioactive decays (data for X_{min} from Yomogida & Matsui 1983).
Some properties of chondrules and matrix in H and L chondrites (data from Scott 2007).
Model parameters for the model of Miyamoto et al. (1981), (MFT81) for a consolidated body (average L chondrite), and for a similar model of an initially porous body without (PL0) and with (PL) additional heating by decay of ^{60}Fe and longlived nuclei.
All Figures
Fig. 1
Variation of a) midplane temperature and b) pressure with distance from the star at different evolutionary epochs (0.2 Ma and from 0.5 Ma to 3 Ma in steps of 0.5 Ma) in a model for the evolution of the solar nebula (onezone αmodel). They define the outer boundary conditions for a planetesimal embedded in an accretion disk. In the left part the regions are indicated where the disappearance of an important absorber results in an extended region of nearly constant midplane temperature. 

In the text 
Fig. 2
Relation between relative density (filling factor) D and maximum experienced pressure for isostatic pressing (top) and run of relative density within planetesimals of the indicated size (bottom) that results from cold pressing caused by selfgravity. 

In the text 
Fig. 3
Specific heat of planetesimal material for mineral mixtures of H chondrites and L chondrites. The dotted line is the analytical approximation of Yomogida & Matsui (1984). 

In the text 
Fig. 4
Contributions of different heating mechanisms to the total heating rate. Time t is after formation of CAIs. The dashed line indicates the formationtime of the H chondrite parent body (2.3 Ma, see Sect. 8). The release of gravitational energy by contraction is shown for the final model of the H chondrite parent body. 

In the text 
Fig. 5
Variation of heat conductivity K with porosity φ. The results for finegrained silica powder (filled circles) from experiments of Krause et al. (2011a,b), and for particulate basalt (crosses) from Fountain & West (1970). The typical grain size is indicated for both cases. Open squares and open circles are experimental results for heatconductivity and porosity for ordinary H and L chondrites, respectively, from (Yomogida & Matsui 1983). The solid line is an analytical fit to the data, according to Eq. (30). 

In the text 
Fig. 6
Temperature variation of the thermal conductivity K. Plotted are data for H chondrites (open rectangles) and L chondrites (open circles), normalised with the φvariation according to approximation (28). The lines connect data for each of the meteorites. 

In the text 
Fig. 7
Temperature evolution of a body of 85 km radius at different depths from the surface and at centre: dotted line at depth 4.25 km, short dashed line at 17 km depth, long dashed line at 23 km depth. The full line shows the temperature evolution at the centre. a) The model is calculated with the same physical input as in the analytical model of Miyamoto et al. (1981). We provide the model parameters of MFT81 in Table 5. b) Similar model, but now calculated for a porous body, considering thermal conductivity of porous material according to Eq. (30), and sintering and cold pressing as described in this paper. We provide the model parameters of PL0 in Table 5. c) Same kind of model as b), but now additional heating by ^{60}Fe and longlived nuclei considered. We provide the model parameters of PL in Table 5. 

In the text 
Fig. 8
Evolution of the radial distribution of the temperature for model PL (see Table 5 for its definition). Note the radial shrinking of the body by compaction of the initially porous material at about 0.6 decay timescales of ^{26}Al. 

In the text 
Fig. 9
Evolution of the radial distribution of the filling factor D (relative density) for model PL (see Table 5 for its definition). Shown is the very initial phase of the evolution where the initially porous material is compacted by sintering. The resulting shrinking of the planetesimal size occurs at about 0.6 decay timescales of ^{26}Al. 

In the text 
Fig. 10
Temperature evolution of test models for porous bodies with modified values of parameters. Left panel with boundary temperature T_{b} = 150 K, right panel with T_{b} = 250 K. Upper panel with heatconduction coefficient K_{b} = 4 W m^{1} K^{1}, lower panel with K_{b} = 2 W m^{1} K^{1}. 

In the text 
Fig. 11
Variation of maximum temperature in the centre of a planetesimal with radius R and instant of formation t_{form}. Upper panel: completely consolidated bodies with porosity φ = 0. Lower panel: models of initially porous bodies with porosity φ_{srf} = 0.5 that consolidate in the interior during the course of their thermal evolution. The lines correspond to the indicated maximum central temperature. At temperatures above ≈ 1400 K partial melting of the mineral mixture is expected. Temperatures exceeding 1500 K are excluded for this reason. Shown are models for three different initial ^{60}Fe/^{56}Fe abundance ratios: Left panel: without ^{60}Fe. Middle panel: the optimumfit value of 4.1 × 10^{7} for the H chondrite parent body, see Sect. 8. Right panel: highest value of 1.6 × 10^{6} given in Table 5. 

In the text 
Fig. 12
Optimumfit model for the cooling history of the parent body of H chondrites. The abscissa is the time elapsed since the formation of the CAIs. We show the evolution of temperature at a number of depths below the surface. The upper contours of shaded areas correspond (from bottom to top) to depths of 0.32 km, 2.3 km, 7.8 km, 11 km, and the highest contour to the centre. The rectangular boxes and circles correspond to the empirical data for H6 and H5 chondrites, respectively, given in Table 6. Crosses are error bars. The dashed lines correspond to the temperature evolution at depths of 8.9 km (lower line) and 36 km (upper line) below the surface. They represent our best fit to the empirical data. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.