© ESO, 2012

1. Introduction

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 utilize 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 heat-producing short-lived radionuclides (²⁶Al and ⁶⁰Fe), which are determined by the accretion timescale, and the terminal size, chemical composition and physical properties of the chondritic planetesimals. Here we use our newly developed code for modelling the structure and thermal evolution of planetesimals of  ≈ 100 km size (Henke et al. 2012, henceforth called Paper I) to constrain the properties and thermal evolution of the H chondrite parent body, and to evaluate the most important parameters that determined its thermal evolution. This is done by finding the optimal fit between a theoretical thermal model and empirical data for the cooling histories of H chondrites. The H chondrite parent body was chosen here as a test case, because it is the only body for which a comprehensive set of thermochronological data is available. For eight H chondrites in particular radiometric ages exist that were obtained using at least three distinct radioisotope chronometers characterised by different closure temperatures. The cooling history of these samples, therefore, is well constrained.

The well preserved and smooth cooling histories of many H chondrites were used to argue that the H chondrite parent body has not been subjected to catastrophic collisions during the first  ≈100 Ma of its history. In that case an onion shell structure develops in which the degree of thermal metamorphism is a function of depth within the parent body; material of petrologic type H6 is found in the central part, while type H3 material is located at the surface. Based on the assumption of an onion shell structure, a number of thermal evolution models were constructed that used increasingly more complex physical input (Minster & Allegre 1979; Miyamoto et al. 1981; Bennett & McSween 1996; Akridge et al. 1998; Trieloff et al. 2003; Hevey & Sanders 2006; Sahijpal et al. 2007; Kleine et al. 2008; Harrison & Grimm 2010). All the previous thermal models of the H chondrite parent body found a radius of about 100 km and an accretion time of about 2 Ma after Ca-Al-rich inclusion (CAI) formation. The success of these models in reconstructing the empirical cooling histories of the H chondrites gave much support to the hypothesis that catastrophic collisions did not significantly disturb the initial thermal structure of the H chondrite parent body. It should be noted, however, that some hints exist for the occurrence of impacts (see Wittmann et al. 2010, and references therein), but these do not appear to have resulted in catastrophic disruption.

These previous thermal models, despite their success in reconstructing the observed empirical cooling rates for several H chondrites, fall short in providing insights into the physical processes occurring during the thermal history of a  ≈ 100 km-sized planetesimals. This motivated us to use our new thermal evolution model (described in Paper I) to find the optimal fit between the empirical cooling history of the H chondrites and the thermal model. Our new thermal model applies recent data on 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), but with some improvements. This theory does not consider more recent approaches for modelling hot pressing for technical processes (e.g. Arzt et al. 1983; Fischmeister & Arzt 1983; Larsson et al. 1996; Storåkers et al. 1999) but appears more appropriate for the lower pressure regime relevant for asteroids. The treatment of heat conduction of the chondritic material and its sintering under pressure and at high temperature still has to be considered as rather crude, but at least is based on physical concepts that are successfully applied in other fields and in some laboratory measurements.

We do not consider a possible melting of the body, because the H chondrite parent body is thought to have remained undifferentiated. Our results will show, however, that this assumption may not be entirely valid, because the calculated central temperature of the H chondrite parent body is only slightly below the solidus temperature of chondritic material (1200 K, Fei et al. 1997).

Our model considers more properties of the asteroidal body than traditional approaches, and so fitting of the model to the observed cooling histories is more difficult because of the much larger parameter space. Therefore we develop an automated method to perform such a fit that is based on the evolution algorithm that mimics the concepts that are thought to rule biological evolution in nature. This is some kind of intelligent trial-and-error method by which the optimum of some quality function can be found. Its advantage is that it does not require any kind of good behaviour of the quality function. It will be shown that this method (in the version described in Charbonneau 1995) can successfully be applied to our problem.

For fitting the model we exclusively use meteorites for which at least three age determinations at different closure temperatures are available. These meteorites do not only constrain the slope but also the curvature of an individual cooling path at a certain location within the parent body. Within the framework of the analytic model of Miyamoto et al. (1981), such a data set for a single meteorite would already completely determine the radius and the formation time of the body. However, in view of the presumably complex thermal evolution of meteorite parent bodies, a larger dataset for several meteorites is necessary for a reliable fit. For the H chondrites, 31 data from 8 meteorites (2 meteorites with 5 data, 3 with 4 data, 3 with 3 data) are available for this purpose. This allows determining all the important parameters of the H chondrite parent body by an optimization procedure.

The plan of this paper is as follows: in Sect. 2 we describe the available empirical data. In Sect. 3 we give a brief overview how we calculate thermal evolution models. Section 4 describes the optimisation method. Section 5 gives our results and conclusions. The paper closes with some final remarks in Sect. 6.

2. Data for cooling histories of H chondrites

Most ordinary chondrites were affected by thermal metamorphism in their parent bodies and exhibit a range of metamorphic conditions from type 3 (unequilibrated) to type 6 (equilibrated) (Dodd 1969). This range in metamorphic conditions probably reflects different heating and cooling histories, presumably caused by varying burial depths within the parent body (Miyamoto et al. 1981; McSween & Patchen 1989; Trieloff et al. 2003; Kleine et al. 2008). The heat source for the thermal metamorphism of ordinary chondrites most probably was heating by decay of ²⁶Al and to a lesser extent also ⁶⁰Fe.

Ordinary chondrites and their components have been dated by a variety of methods (Nyquist et al. 2009). Al-Mg ages for chondrules from L and LL chondrites indicate that chondrule formation occurred approximately 2 Myr after formation of CAIs, commonly regarded as the oldest objects formed in the solar system (Kita et al. 2000; Rudraswami & Goswami 2007; Rudraswami et al. 2008). The formation age of chondrules from H chondrites is less well known because no primitive H chondrite has yet been investigated with the Al-Mg system. However, a Hf-W age for the H4 chondrite Ste. Marguerite appears to date chondrule formation at 1.7  ±  0.7 Myr after CAI formation (Kleine et al. 2008, 2009), about contemporaneously to chondrule formation in L and LL chondrites.

The ages for chondrules from ordinary chondrites mainly constrain the timescale of parent body accretion (Alexander 2005; Kleine et al. 2008). In contrast, ages for equilibrated chondrites can be used to reconstruct the cooling history of ordinary chondrite parent bodies, by applying chronometers having different closure temperatures to a set of equilibrated chondrites of different metamorphic grade. By far the most comprehensive set of such age data exists for H chondrites (see Table 1), providing temporal information for cooling over almost the entire range of temperatures prevailing from the metamorphic peak down to ambient temperature. The closure temperatures for diffusive exchange of parent and daughter elements among the different minerals in a rock (Dodson 1973; Ganguly & Tirone 2001) for radiometric methods applied to H chondrites range from  ≈ 1150 K (Hf-W) to  ≈550 K (Ar-Ar). Additional information on the cooling history is provided by ²⁴⁴Pu fission tracks in orthopyroxene (corresponding to cooling below  ≈550 K) and merrillite (corresponding to cooling below  ≈390 K) (Pellas et al. 1997; Trieloff et al. 2003), and from metallographic cooling rates (corresponding to cooling between  ≈800 and  ≈600 K), which employ Fe-Ni diffusion profiles in metal grains consisting of kamacite and taenite (e.g., Herpfer et al. 1994; Hopfe & Goldstein 2001).

The age and thermochronological data available for H chondrites are summarised in Table 1, along with information on the closure temperatures of the different methods. These age data permit a detailed reconstruction of the cooling history of H chondrites, and form the basis for numerically simulating the accretion and cooling history of the H chondrite parent body (Sect. 5).

For calculating heating rates the difference between the formation time of CAIs and the ages from Table 1 has to be formed. This sets the zero-point of our time-scale and the instant for which the initial abundance of ²⁶Al is prescribed. As reference time we take the value of 4568.5 ± 0.5 Ma from Bouvier et al. (2007).

3. Thermal evolution models

The temperature history of an undifferentiated meteoritic parent body is assumed to be determined by essentially three processes: (i) transient heating by decay of short-lived radioactive nuclei ²⁶Al and ⁶⁰Fe and long lasting heating by long-lived radioactives (⁴⁰K, ²³²Th, ²³⁵U, ²³⁸U), (ii) transport of heat to the surface by heat conduction, and (iii) energy exchange with the environment by emission and absorption of radiation energy. The heat transport in the body is governed by the fact that the material is porous, that heat conductivity of chondritic material strongly depends on porosity (e.g. Krause et al. 2011b), and that the porosity changes by sintering as the body is heated up.

For the formation of the bodies we apply the “instantaneous formation” approximation, i.e., it is assumed that the bodies acquired most of their mass during a period short compared to the half-life τ_1/2 = 7.2 × 10⁵ a of the main heating source, ²⁶Al. Present planetary formation scenarios suggest such rapid mass acquisition of 100 km sized and bigger bodies (e.g. Weidenschilling & Cuzzi 2006; Nagasawa et al. 2007). The rapid growth is approximated by the assumption that the bodies came into life at some instant t_form and have constant mass over their subsequent evolution.

For reasons of computational economics the shape and internal structure of the body is assumed to be spherically symmetric, though this is certainly a strong simplification. In detail, the following equations are solved:

• 1.

  The equation of hydrostatic equilibrium(1)for the pressure p, where (2)and (3)The density ϱ_b is the mass-density of the consolidated material and D the fraction of volume filled with matter. The void fraction or porosity of the material is φ = 1 − D.

• The heat conduction equation for the evolution of the temperature (4)where c_v is the specific heat of material per unit mass, h the specific heating rate by decay of short living radioactives, F_r the local gravitational acceleration and v_r the shrinking velocity during sintering. The heat current is (5)with the heat-conduction coefficient K. The heat conductivity used in our model depends on the porosity φ = 1 − D of the material. We use the analytic fit (6)to experimental data discussed in Krause et al. (2011a) and in Paper I, where φ₁ = 0.08, φ₂ = 0.167, and a =  −1.1 are constants given in Paper I¹. The pre-factor K_b corresponds to the heat conductivity of the bulk material (i.e., at D = 1). The average value of K_b is determined by extrapolating data obtained from a couple of ordinary chondrites (Yomogida & Matsui 1983) by means of Eq. (6) to zero porosity. Its value is K = 4 W m^-1 s^-1, but there is significant scatter.

• 3.

  An equation for the sintering of the initially porous material under the influence of pressure (7)where the right hand side is determined by solving a set of equations for the specific sintering model of Kakar & Chaklader (1967) and Rao & Chaklader (1972).

The internal structure and thermal evolution of parent bodies of planetesimals is calculated by the numerical method described in Paper I.

The radius of the body is not fixed but changes during the evolution by sintering of the initially porous material. If the body would be completely sintered the final radius would be R. This radius is related to the mass by M = (4π/3)ρ_b   R³ with ρ_b being the bulk density of the completely consolidated material. Initially the material is porous and the radius is bigger because the density of the porous material ρ = Dρ_b is lower. The initial porosity is φ_srf. This equals also the porosity of the outermost surface material which is too cool and under too low pressure for sintering.

The problem depends on a number of parameters, the most important of which are the bulk density ρ_b of the chondritic material, the mass M of the body or its radius R, the formation time t_form (measured relative to the formation time of CAIs), the heat conduction properties of the material (K_b), the mass fractions X_Al and X_Fe of Al and Fe, respectively, in the chondritic material, and the abundance ratios of short-lived radioactives, ²⁶Al/ ²⁷Al and ⁶⁰Fe/ ⁵⁶Fe, the porosity φ_srf of the surface material, and the surface temperature T_srf.

The surface temperature T_srf depends during the initial phase of evolution on the properties of the ambient dusty accretion disk, after disk dispersal on the distance to the proto-sun. Since model properties depend only weakly on T_srf, this quantity is held fixed over the whole evolution period in this study.

Some of the parameters are fixed by the properties of the H chondrite material and general properties of the solar system. They are given in Table 2 and their values are chosen as in Paper I. The remaining parameters M or R, t_form, T_srf, φ_srf, and K_b are special properties of individual bodies and are varied in our model calculations.

Since there seem to be doubts that the ⁶⁰Fe/ ⁵⁶Fe ratio (at time of CAI formation) is the same for all bodies (see Quitté et al. 2010) we prefer to consider also this quantity as a free parameter (details are given in Sect. 5.1).

Fig. 1 Schematic representation of the temperature evolution at two different depths (full lines) and at the centre (dotted line) of a planetesimal of the 100 km size class. The horizontal dashed lines correspond to fixed temperatures equal to the closure temperatures in Table 1. The circles indicate the instants at which at a given depth the temperature drops below a given temperature. These instants correspond to the closure ages of the different thermochronometers used for meteorites.

4. The fitting procedure

4.1. Comparing thermochronological data with asteroid models

The data on the cooling history of meteorites described in Sect. 2 associate with a meteorite a number of temperature-time data pairs, corresponding to the instants where the temperature at the pristine location of the meteorite inside of its parent body fell below that temperature. The set of data pairs obtained for a meteorite corresponds to a set of intersection points between the temperature evolution curve of some mass-element with the lines of constant temperature on the descending part of the evolution curve. This requires, of course, that the temperature evolution at the location of the later meteorite proceeds without substantial disturbance by external influences (impacts, for instance).

Figure 1 demonstrates this for the four closure temperatures (1150 K, 720 K, 550 K, 390 K) of Table 1. The cooling curves are shown as full lines, the lines of constant temperature as dashed horizontal lines. The circles on a cooling curve correspond to the data for a hypothetical meteorite. Since the average slope of the curve connecting two neighbouring data points significantly varies for curves corresponding to different depth’s, it is obviously possible to reconstruct for a given model from only two data points the appropriate depth for which the temperature evolution curve runs through just these two data points.

If the properties of the body are not known and only two data points are available for a meteorite, such a fit would be possible for a whole range of model parameters. This does not allow any conclusion on the properties of the body. If, however, more than two data points are available for a meteorite this restricts the possible models since the overall shape of the temperature evolution curves depends on the model properties. The average slope of that part of the curve that connects neighbouring data points on the curve varies for different such pairs with model properties in a different way. This allows to reconstruct the properties of the body from cooling data on meteorites by looking for a set of model parameters for which for each of the meteorites an appropriate depth can be found such that the cooling curve runs through all data points.

Since for real meteorites, experimentally determined closure temperatures and ages are subject to considerable errors (cf. Table 1), data sets for an as big as possible number of meteorites (with at least three data points) are required and the task is to find a set of parameters that describe the properties of the parent body such that for each of the meteorites a depth below surface can be found such that the (somehow defined) minimum distance of all of the data points of a meteorite from the corresponding temperature evolution curve becomes as small as possible; and this should hold for all available meteorites. Hence one has to solve a minimisation problem, where some quality function has to be minimised with respect to a set of data and a set of parameters.

In the past, this task has been solved mainly by trial and error. Formation times t_form and parent body radii R for some bodies have been determined this way (Miyamoto et al. 1981; Bennett & McSween 1996; Akridge et al. 1998; Trieloff et al. 2003; Hevey & Sanders 2006; Sahijpal et al. 2007; Kleine et al. 2008; Harrison & Grimm 2010). However, the essential properties of an asteroidal body that determine its thermal evolution, are not merely R and t_form but also some additional parameters, see Table 3. To these, one has to add the unknown depth’s for each meteorite as additional parameters of the optimisation problem. Here we aim to solve the problem by a general numerical algorithm that allows to treat the more extended parameter space (and possibly additional parameters in future calculations) and large data sets.

4.2. The data set

For reconstructing the parent body of a group of meteorites that are thought to derive from the same asteroidal body, we use all specimens for which at least three thermochronological data points are available. For the parent body of, e.g., the H chondrites, the eight meteorites given in Table 1 have sufficiently accurate determined data and the required minimum number of data points. Three of them have even four data points, the other five have three of them.

The burial depth of each of the J meteorites, enumerated by an index j (=1,...,J), in the parent body is characterised by the J mass-coordinates M_j in the theoretical model of the parent body, because the model calculation uses the mass-coordinate as independent variable (see Paper I). Then we have a list of I data quintuples (8)with T_i being the closure temperature, t_i the corresponding age, σ_T,i and σ_t,i the errors of T_i and t_i, respectively, and j_i a pointer to the appropriate member in the list (M_j,j ∈  { 1,...,J } ) of mass-coordinates corresponding to the different meteorites.

The mass-coordinates M_j are not known a priori. They have to be determined simultaneously with the parameter set given in Table 3 by the optimisation process.

4.3. The quality function

The quality function to be optimised can be defined in different ways. We construct in this paper one in the following way: For each of the mass-coordinates M_j it is known which of the data sets ϖ_i are associated with this M_j. During a model run, the instants when the temperature drops below the closing temperatures belonging to the data sets ϖ_i associated with j, are recorded. This defines for each i an age t_mod,i on the temperature evolution curve j where it passes through the prescribed temperature T_i. Analogously, the temperature T_mod,i is recorded at instant t_i for each of the data points.

As quality function for the optimisation process we choose the least squares deviation of the model age t_mod,i from the empirically determined age t_i, and of the model temperature T_mod,i from the closure temperature T_i: (9)this χ² is to be minimised. The parameters (R, t_form, ⁶⁰Fe/⁵⁶Fe, T_srf, φ_srf, K_b, M_j (j = 1,...,J)) are then determined such that each cooling curve corresponding to the burial depth of a certain meteorite in the model passes as close as possible to the empirically determined age-closure temperature pairs for all data points of that meteorite.

4.4. Choice of the optimisation algorithm

The basic problem is to determine the absolute minimum of the quality function Φ(ϖ) and possibly all other deep local potential minima with comparable depth. As is well known, such an optimisation problem with a large number of variables cannot be solved by analytical methods and the task of solving the problem by numerical methods becomes extremely difficult and time consuming if the surfaces Φ(ϖ) = const. have a complicated geometric structure, which is expected to be the case for our problem. No general method is known which allows to definitely

• 1.

  find all local extremal values of an optimisation problem;

• 2.

  and to decide if the real optimum indeed is contained within the set of extremal values found by the applied method.

Methods have been developed which allow one to find an in a certain sense best result, which means that for test cases with a known absolute minimum the algorithm always finds this minimum or at least a local minimum that is not much worse than the absolute one.

One method, that has successfully been applied to several difficult optimisation problems like that of the travelling salesman is the “evolution algorithm” developed by Rechenberg (1973, 1989). A detailed description of how this method can be applied to problems in astronomy and astrophysics and a code are given in Charbonneau (1995). We use this method to reconstruct the parent body’s properties of a suite of meteorites.

The basic idea is to imitate the concepts that are thought to rule the evolution of species in an algorithm to find the extremal values of some quality function. Evolution theory describes the change of the shape of creatures in a population during many generations due to the variability of the population and the dependence of the amount of the creature’s offspring on their conformity to the environment. The biologic term genotype means the whole heritable attributes of a creature, while phenotype denotes its appearance. During the development of the creature the phenotype forms according to the instructions of the genotype. The position in the genotype, where a certain property (e.g. eye colour) of the individual is encoded, is called gene. The value of the property (e.g. blue) is called allele. The amount of possible offspring (the reproduction outcome) depends on, how advantageous its alleles are. The reproduction outcome determines, how numerously certain alleles in a population are inherited compared to other respective alleles. Which processes are advantageous is determined by the environment. Thus after some generations advantageous alleles will accumulate in the population, while less advantageous ones will be depleted, such that after a certain number of generations only those alleles, the creature benefits most from, remain in the population in a noteworthy amount. This leads to a degeneration of the gene pool. But additionally mutations can cause random change of the genotype in a certain allele, which can cause change of the affected gene in a significant way and can create a completely different phenotype. Thus mutation adds novel alleles to the gene pool. These can be superior to the established alleles (even if they are destructive in most cases). So the mutations keep the variability of a population, and allow the population to adapt to changing environments.

An algorithm for finding parameters that works like evolution due to inheritance and natural selection is superior in finding of optimal parameters in ugly problems with many local minima compared to methods that, e.g., calculate in some way the derivative of the quality function and then go into the direction of downward slope and repeat calculating the derivative and so on. Those methods readily find the next local minimum which, however, might not be the absolute one. In the evolution algorithm mutations can produce extensive changes of the parameters to the parents’ ones. This can cause big jumps in the parameter space. Thus the population will not remain in a local minimum, but test remote regions again and again and perhaps find a deeper local minimum that is possibly a total one.

4.5. General implementation of an evolution algorithm

An evolution algorithm can be applied to every problem in which the variables to be optimised can be encoded to form a string of numbers (chromosome). This string then is the genotype. It is assumed in the following that for every individual the model is calculated and compared to the data by calculating a quality function, like Eq. (9), e.g., and ranked by this to other individuals. After determining which individuals are to get off-springs, their genotypes are cut at the same random digit(s) and the related parts are exchanged, so that, e.g., the initial part of parent 1 is connected to the end part of parent 2. The so generated string then represents the genotype of the new offspring individual. Mutation can be performed by giving every digit of the string of the genotype a certain possibility to take a random value.

An implementation of an evolution algorithm can be as follows:

• 1.

  A random initial population is constructed with normally noreference to the problem.

• 2.

  The fitness function of each individual is calculated by a user given algorithm which determines how well the individuals match the problem.

• 3.

  A new generation is constructed by breeding of selected individuals of the old population. There the algorithm creates new parameters by exchange of information between the parameters of the two parent individuals. Also the mutation takes place in this step. In most cases the number of individuals in each generation will remain the same.

• 4.

  The old population is replaced by the new one.

• 5.

  Then back to 2.

4.6. Application to our problem

The particular implementation of a genetic algorithm that we use is that described in Charbonneau (1995). It uses a random number generator that generates numbers between 0 and 1 to derive a set of random parameter values for formation time, radius and all other parameters of the model. Each time the genetic algorithm requires a new value of Φ for some individual it provides a set of K random numbers 0 ≤ ξ_k ≤ 1 (k = 1,...,K) from which values for the corresponding physical parametres x_k of the model are calculated as (10)Here x_k,min and x_k,max are the minimum and maximum allowed values for x_k according to Table 3. For this set of parameter values x_k an evolution model is calculated and from its results the corresponding value of the quality function Φ.

The burial depth’s of the meteorites are treated as part of the parameter set of the model. Random values are assigned to the burial depths of each of the distinct H-chondrites. In this way for each of the individuals of the evolution algorithm, a certain depth is ascribed to every chondrite. Technically this is done by assigning to each meteorite a layer characterised by an index i and a position factor y that gives its relative position inside of the layer with respect to the edges of the layer, r_i and r_i − 1, so that the real position of the meteorite is (11)Because the sintering behaviour is treated as the same everywhere within a single layer the factor y will not change during the thermal evolution. The temperature at the position of the meteorite is calculated by (12)and it is tested how good the temperature evolution of that T_met matches the empirical data of the corresponding meteorite by means of Eq. (9).

According to Charbonneau (1995) a fit achieved by a genetic algorithm is moderately good if I − K ∽ χ², where I is the number of data points to be fitted and K the number of the optimised parameters. For example, for the six parameters given in Table 3, plus the eight layers of the eight meteorites (K = 6 + 8 = 14), and the 28 data points, one gets a desired value of χ² < 14 for an acceptable fit.

5. Results for the parent body of H-chondrites

5.1. Parameters

Besides the physical parameters of the problem, that are not subject of the fitting procedure, we have a number of parameters that describe properties of the body and external quantities. These are given in Table 3, including the range of variation admitted during the optimisation process. The following parameters require some comments:

Radius: There are a number of reasons (Gaffey & Gilbert 1998; Morbidelli et al. 1994) why asteroid (6) Hebe is thought to be the parent body of the H chondrites and the IIE iron meteorites. Some general characteristics of (6) Hebe are listed in Table 4. If (6) Hebe is only a fragment of the original H chondrite parent body, then the size of (6) Hebe provides a firm lower limit for the size of the H chondrite parent body. Hebe is a somewhat irregular shaped body (Torppa et al. 2003) with dimensions 205 × 185 × 170 km corresponding to an average diameter of 186 km. Its high mass-density means that it is not a rubble pile as most other asteroids. It did not suffer catastrophic collisions that disrupted the whole body and where the debris later re-assembled to the present body. Nevertheless, a surface layer of unknown thickness has certainly been eroded over time by collisions with other bodies and the present size of Hebe defines only a lower limit to the size 4.56 Ga ago. But this lower limit should be observed in the model construction.

⁶⁰Fe/ ⁵⁶Fe-ratio: The decay of ⁶⁰Fe may have contributed to the heating of the body, in particular since the recently re-determined half-life of τ_1/2 = 2.62 ± 0.04 Ma (Rugel et al. 2009) significantly exceeds the half-life of ²⁶Al of τ_1/2 = 0.72 Ma. Though there is less ⁶⁰Fe than ²⁶Al, the heating period can be significantly prolonged in models with ⁶⁰Fe compared to models with only ²⁶Al included.

It has been questioned (cf. Dauphas et al. 2008; Quitté et al. 2010) that ⁶⁰Fe was homogeneously distributed throughout the solar system. If true, this would mean that the initial ⁶⁰Fe abundance in the parent body is not known a priori. On the other hand, a recent study by Tang & Dauphas (2012) argued for a homogeneous distribution of ⁶⁰Fe on a relatively low level of 1.08 ± 0.21 × 10^-8 for the initial ⁶⁰Fe/ ⁵⁶Fe abundance ratio at time of CAI formation. This is much lower than values previously reported, although the authors do no comment on the credibility of earlier studies.

At such a low abundance heating, by ⁶⁰Fe is nearly negligible. Because of the presently unclear situation with respect to which value should be taken for the ⁶⁰Fe/ ⁵⁶Fe abundance ratio we treat this as an additional free parameter for the optimisation process. The range of ⁶⁰Fe/ ⁵⁶Fe ratios considered (see Table 3) is between no ⁶⁰Fe and the maximum reported value (Birck & Lugmair 1988).

Fig. 2 Results for the individual meteorites used in the model fit. The values of the free parameters and the burial depth’s of the meteorites found for the optimised model are given in Table 5. The full line shows the temperature evolution at the burial depth’s of the meteorites in the asteroid, as determined by the evolution algorithm. The dashed line shows the temperature evolution at the centre and the dotted line is the surface temperature. The square boxes indicate the individual data points and the error bars the uncertainty of closing temperature and age determinations.

Surface temperature: This is determined by the unknown radius of the orbit, its eccentricity, by the rotation of the body, the surface albedo, and, during the first few million years, by the properties of the ambient accretion disc. All these factors are not known precisely. Therefore we think it best to include the surface temperature T_srf as a free parameter in the optimisation process. The resulting temperature T_srf then represents some kind of average temperature of the surface layers during the evolution period where the body first heats-up and then finally cools down below the lowest of the closure temperatures.

The resulting value for T_srf has no bearing to the present surface temperature of a possible successor body, because in the meantime the body may have changed its properties by catastrophic collisions and it may have changed to some extent its orbit during the evolution of the asteroid belt. Surface temperatures in the present-day asteroid belt therefore cannot be considered as representative for the surface temperature during the phase where the essential metamorphic processes acted in the body, as it is assumed in some thermal evolution models.

The range of surface temperatures (see Table 3) is limited by the fact that T_srf should be higher than the temperature where the body would contain ice ( ≈ 150 K) and lower than the lowest of the closing temperatures ( ≈ 400 K).

Surface porosity: During the (short) main growth phase of the body the structure of the surface layer is determined by the properties of the smaller planetesimals that form the body under consideration. Later the surface layer is gardened by continued impacts that produce a regolith surface layer for which theoretical calculations indicate a porosity in the range 20% to 30% (Richardson 2009; Warren 2011). Presently, too little is known on these processes as that they can be included in our model calculation. Therefore we include the porosity of the surface layer, φ_srf, as a free parameter in the optimisation procedure. The resulting value then is some sort of average value of φ_srf over the period where the metamorphic processes were active.

5.2. Data

The general optimisation method of the evolution algorithm is applied to the data for the set of H chondrites given in Table 1. For reasons that will become clear from the model results shown in the subsequent figures and are discussed later the two highest closing temperature data of Ste. Marguerite and the plutonium fission track datum of Nadiabondi are not used. Then we have totally 28 data points that have to be fitted. The problem depends on the 5 parameters given in Table 3 and the eight unknown burial depths of the meteorites.

The range of surface porosities (see Table 3) should in principle be restricted by the requirement, that the porosity should be higher than that of the most porous chondrites (about 20%, cf. Consolmagno et al. 2008). During the model calculations it was found that for an optimum fit values not much bigger than this minimum value are required and, therefore, no a priori restriction was imposed on φ_srf during the final optimisation runs. The optimal set of parameters, indeed, satisfies the assumption.

5.3. Model with fixed heat conductivity

An optimisation was run for the 5 + 8 = 13 parameter and the 28 given data points. The model used 250 mass-shells. The optimisation was run for about 5500 generations with hundred individuals per generation, i.e., totally a number of about 550 000 complete thermal evolution models were calculated for a complete optimisation run. This used about seven days on a normal desk-top computer equipped with a dual-core pentium processor operating with 2.6 GHz.

The finally accepted optimum fit had a value of χ² = 3.1420. This value for χ² can be considered as excellent in view of the fact that a good fit would require only χ² ≲ 15 (cf. Sect. 4.6). The model parameters of the final optimum fit for the properties of the body and the burial depths of the meteoritic material are given in Table 5.

The central temperature at maximum, T_c = 1163 K, that results for the set of optimum parameters is below the Fe-FeS eutectic at 1220 K (Fei et al. 1997). For the optimum model the temperature in the centre therefore is below the border temperature where differentiation would just become possible. In so far the model is compatible with the basic assumption of the model calculation that melting can be neglected.

The cooling curves of the meteoritic materials at the resulting burial depths and the empirical age-temperature data points given by the radiometric ages are shown in Fig. 2 for the optimised model. The solution of the optimisation problem found here meets the data for the cooling history of the meteorites of Table 1 within their 1σ-error-boxes in all eight cases, though in a few cases only marginally. The fit between model and empirical data can be considered as excellent for the three meteorites with four data points: Kernouvé, Richardton, and Estacado. Also for three of the meteorites with three data points, Guareña, Allegan, and Forest Vale, the model fits the data-points rather well. Only for Nadiabondi and Ste. Marguerite the fit is of low quality. These two cases are discussed later.

The maximum temperature, T_max, reached during the thermal history at each of the different depths are given in Table 5. For comparison, the table also gives the metamorphic temperature of the different petrologic types, T_met (Slater-Reynolds & McSween 2005). The maximum temperatures of the different meteorites at their burial depth, T_max, fit reasonably well with the requirement T_max > T_met, except for Ste. Marguerite.

The radius of 132 km required for the optimum model is significantly larger than the present day radius of (6) Hebe. If this asteroid is really the parent body of the H-chondrites this would mean, that a layer with an average thickness of about 40 km has been eroded over time from the surface of (6) Hebe. This is of the same size as the depth of the layer where the petrologic type changes from H5 to H6. The change from initial to present radius amounts to 29% or corresponds to a volume change of 64%, i.e., the body would have lost most of its initial mass. Consequently, all H3 and H4 chondrite material should have been removed, which is inconsistent with meteoritic record. Probably our optimised radius is too big and the most likely reason for this is that the treatment of heat conduction in the porous medium and of the sintering process is not yet realistic enough. The other possibility is that H-chondrites do not come from (6) Hebe, although the arguments favouring (6) Hebe as parent body are strong (Gaffey & Gilbert 1998).

The optimum value for the ⁶⁰Fe abundance is compatible with the most recent abundance determination by Tang & Dauphas (2012). At this low abundance, heating by ⁶⁰Fe contributes only marginally to the total heating of the parent body.

The thickness of the remaining porous outer layer is also given in Table 5. As lower boundary of this layer we define the radius where φ = 10%. At this degree of porosity usually the connections between pores start to close and the pore structure changes from a network of interconnected pores to isolated voids.

The low value of the surface porosity φ_srf = 0.28 may either result from the fact, that planetesimals of the 100 km size-class gain most of their mass by collisions with smaller planetesimals that have already diameters of about 10 km and more and, therefore, are already largely compacted (see Fig. 10 of Paper I) if they are incorporated into the growing larger body. A second mechanism that may be responsible for the low value of φ_srf is that our optimised value represents a mean value of the actual value of φ_srf over the initial heating and the extended subsequent cooling period of 100 Ma. During that period impacts certainly continue with a decreasing rate that do not contribute to a significant further growth, but that form a regolith surface. According to model calculations this has a saturation porosity of 20% to 30% (Richardson 2009; Warren 2011) that would fit to the value of φ_srf = 0.28 found by the optimisation process.

Fig. 3 Results for a model fit where the heat conductivity at zero porosity, Kb, is included in the optimisation process. The values of the free parameters and the burial depth’s of the meteorites found for the optimised model are given in Table 6. For meaning of symbols see caption to Fig. 2.

5.4. Model with optimised heat conductivity

The pre-factor K_b in our approximation for the porosity-dependent heat conductivity, Eq. (6), corresponds to the heat conductivity of the bulk material (i.e., at φ = 0 or D = 1), determined by extrapolating data obtained from a couple of ordinary chondrites (Yomogida & Matsui 1983) to zero porosity. Figure 6 in that paper and our Fig. 6 shows that there is considerable scatter in the data which results in some uncertainty of the value of K_b. The value for K_b used in the above model calculation (see Table 5) may therefore be an unfavourable choice from the allowed range of values compatible with experimental values. Additionally, there is the problem that the small suite of chondrites used for determining K_b may not be representative for the bulk of the material of the parent body.

Therefore we performed a model calculation where K_b was considered to be a free parameter of the problem for which an optimal value is determined by the optimisation procedure. The remaining number of data-points for closure times with which to compare is still big enough to allow for an additional parameter, because we now have 28 data points and 14 parameters (6 for the body and the 8 burial depth’s).

The model, again, used 250 mass-shells. The optimisation was run for about 5500 generations with hundred individuals per generation. The finally accepted optimum fit had a value of χ² = 6.116 in generation 2945. With this value for χ² the fit can also be considered as an excellent one in view of the fact that an acceptable fit would require only χ² ≲ 14 (cf. Sect. 4.6). How the value of χ² evolved during the optimisation process for this model is shown in Fig. 5. The somewhat “noisy” appearance of the evolution history reflects the permanent mutations and crossing of genes in the evolution algorithm that not always results in children that are as fit as their parents.

The model parameters of the final optimum fit for the properties of the body and the burial depth’s of the meteoritic material are given in Table 5 and the model results for the individual meteorites are shown in Fig. 3. The already good fit to the data points shown in Fig. 2 is once more slightly improved. The two meteorites Nadiabondi and Ste. Marguerite, however, drop out again and are not well fitted.

The thickness of the remaining porous outer layer is also given in Table 6. This layer is defined as that region where φ ≥ 0.1; its lower boarder is given by the radius where φ = 0.1. Around this border the heat conductivity K(φ) drops nearly discontinuously from the bulk value K_b to the low value K(φ_srf). The blanketing effect of this layer is responsible for the shallow burial depths of the H4 and H5 chondrites found in the model, which are also found in other models (e.g Bennett & McSween 1996; Akridge et al. 1998; Harrison & Grimm 2010) that introduce an insulating outer regolith layer instead of modelling the sintering process.

In this model the radius of the body is smaller than in the model with prescribed K_b. With R = 111 km it is only 18 km bigger than the average radius of 93 km of (6) Hebe. If (6) Hebe really is the parent body, this would mean that material corresponding to 16% of the initial radius or 41% of the initial mass has been knocked-off by collisions. This model seems to be more plausible because the body then retained most of its mass and did not lose most of its mass as in the preceding model. This makes it more plausible that we still find H4 chondrites.

The optimum value for the ⁶⁰Fe abundance agrees well with the most recent abundance determination of 1.08 ± 0.21 × 10^-8 by Tang & Dauphas (2012).

The value for K_b required for an optimum fit is significantly lower than the average value for K_b obtained from the laboratory data but still inside the range of experimental values. This is shown in Fig. 6 where we plot the ratio K_exp/K(φ) for all chondrites from Yomogida & Matsui (1983). If one would have used only the experimental values for the H chondrites, the optimised value for K_b would be definitely too low. This may result from the fact that, as already mentioned, the H chondrites used by Yomogida & Matsui (1983) for measuring the heat conductivity and its dependence on φ are unfortunately not representative for the parent body as a whole. For instance could the material in the more central part of the body have for some reason different properties than the zones probed by the meteorites. Another probability could be that we neglected the temperature variation of the heat conductivity. If one looks at Fig. 6 this variation seems to be small, but in fact the heat conductivity drops with increasing temperature and may be significantly lower at 1000 K than the values seen in Fig. 6. Therefore the low value found in the optimisation may also simply mean that if we optimise K_b, the resulting value is kind of a “best average” between the lower conductivity of the hot core region of the parent body and of the higher conductivity in the cooler mantle region. The temperature dependence of the heat conductivity therefore should be considered in future model calculations.

5.5. Model with given Radius of (6) Hebe

We also tried a model fit with the present-day radius R = 93 km of (6) Hebe to check whether a somewhat sub-optimal but, nevertheless, acceptable fit could be obtained for this case. The value of K_b is included in the optimisation process, but R is now fixed. The model, again, used 250 mass-shells. The optimisation was run for about 5500 generations with hundred individuals per generation. The finally accepted optimum fit had a value of χ² = 6.058 in generation 5461. With this value for χ² the fit is not worse than the fit of the model with fixed K in view of the fact that an acceptable fit would require only χ² ≲ 15 (cf. Sect. 4.6). The resulting set of parameters is given in Table 7 and the temperature variation at the burial depths of the meteorites is shown in Fig. 4.

The main difference between the models with optimised R and fixed R is that the resulting T_srf in the first model is 178.2 and 250.4 K in the latter, all other things being very similar. The higher surface temperature is required to prevent a too rapid cooling of the outer layers of the smaller body. In practice this means that there is some ambiguity in the models with respect to the location of the body (surface temperature should decrease with increasing distance to the sun) and its radius. The optimisation process cannot really discriminate between the two cases considered here.

Fig. 4 Results for a model fit where the radius is set equal to the observed radius of (6) Hebe. The burial depth’s of the meteorites found for the optimised model are given in Table 7. For meaning of symbols see caption to Fig. 2.

Fig. 5 Evolution history of the value of χ2 during the optimisation process with progressing number of generations, n, for the model where the heat conductivity at zero porosity, Kb, is included in the optimisation process. Shown is for each generation the particular individual with the lowest value of χ2 of all individuals of that generation, the formation time, tform, of this individual, and the radius of the body, R. For two representative meteorites, Kernouvé and Allegan, the evolution of the optimum value of the burial depth, d, is also shown.

Fig. 6 Data for thermal conductivity Kexp of ordinary chondrites (from Yomogida & Matsui 1983) divided by the analytical fit K(φ) given by Eq. (6) using a value of Kb = 1.75 obtained as optimum value for fitting the cooling history. Full black circles represent data for H chondrites, open circles data for L chondrites. The dashed line corresponds to the value of K = 4 used in model with fixed K.

5.6. Nadiabondi and Ste. Marguerite

For Nadiabondi and Ste. Marguerite, some data points need to be discussed separately.

For Ste. Marguerite the temperature increase during metamorphism was insufficient for resetting some of the radioisotopic systems, making it difficult to constrain the cooling history of this meteorite. For instance, Kleine et al. (2008) argued that the closure temperature of the Hf-W system is much higher than the metamorphic peak temperature of Ste. Marguerite. The Hf-W age of 1.7 ± 0.7 Ma after CAI formation, therefore, must date a high-temperature event prior to the metamorphism. This event most likely was chondrule formation. This interpretation is consistent with our finding that the H chondrite parent body probably accreted at about 1.8 Ma after CAI formation.

Second, some data of relatively poor precision do not really constrain the fit procedure, e.g., in the case of Nadiabondi the low temperature datum of fission track retention of the phosphate merrillite. Nominally it yields an – impossible – negative cooling rate, and only within errors just relatively fast cooling can be inferred. Similarly, both fission track and Ar-Ar data of Ste. Marguerite have large uncertainties, precluding firm constraints on the cooling history, although the nominal values plot on the fitted cooling-curves.

5.7. A Miyamoto-like model

For comparison, the optimization was also done for a model where the porosity was set to a value of φ_srf = 0.08. This is a typical value as it is observed for various H6 chondrites that have not been subject to strong shock modifications (Yomogida & Matsui 1983; Consolmagno et al. 2008). With this initial porosity the body is practically already compacted from the beginning on and sintering does not play any role. The heat conductivity then is practically constant throughout the body and for all times. The prefactor K_b of the heat conductivity was set to the same value as in Sect. 5.3 and was not included in the optimisation process. Such a model corresponds essentially to the kind of model considered originally by Miyamoto et al. (1981), the mayor difference being that heating by ⁶⁰Fe is included.

The model, again, used 250 mass-shells. The optimisation was run for about 3260 generations with hundred individuals per generation. In generation 2299 an optimum value of χ² = 5.94 was obtained that was considered as sufficiently good. The model parameters of the optimised model are shown in Table 8 and Fig. 7 shows the fit for the individual meteorites.

It is obviously possible to obtain also an acceptable fit to the observed cooling history of the meteorites with such kind of model. The parameters of such a model are somewhat different from the preceding cases because of the lack of an insulating porous outer layer and a higher heat conductivity during the initial warm-up phase. In particular a significant contribution of ⁶⁰Fe decay to the heating is required to compensate for the more easy heat-loss in the present case. Though the model as a whole does not look worser than the model with optimised K_b, it is physically not acceptable since the lack of a porous outer layer is unrealistic.

Fig. 7 Results for the individual meteorites if a Miyamoto-like model is used for the model fit. The burial depth’s of the meteorites found for this model are given in Table 8. For meaning of symbols see caption to Fig. 2.

5.8. Comparison with earlier results

There are a number of previous attempts to reconstruct the properties of the H chondrite parent body using the empirically determined cooling rates of H chondrites derived from different depth within the parent body. Table 9 gives an overview on the results and the values of some key parameters of some of these models. The table lists all the models that attempted to obtain in some way a “best fit” between properties and cooling histories derived from laboratory studies and a thermal evolution model of the parent body. The results are generally the putative radius and formation time of the body, and the burial depths of the meteorites. All studies, including ours, showed, that the observed properties of the H chondrites are compatible with a parent body of about 100 km radius and an “onion shell” model. In this model the body is first heated within a period of ca. 3 Ma by decay of radioactive nuclei, mainly ²⁶Al, and then cools down over an extended period of ca. 100 Ma, without suffering catastrophic collisions that would have disrupted the body, thereby disturbing its thermal evolution record.

Despite the different assumptions made in the model calculations, from the rather simplistic model of a homogeneous body with constant material properties, that can be treated analytically (Miyamoto et al. 1981), to models with rather complex implemented physics like that of Harrison & Grimm (2010) and ours, the derived properties of the parent body are similar. These do not seem to depend very critical on the particular assumptions. The more strongly deviating results for the burial depths of the meteorites cannot be used to discriminate between models, because no independent information on this is presently available. This would require that methods are developed to determine the maximum pressure experienced by a meteorite.

The most important difference between our model and previous models is that we included sintering of the granular material in the modelling. This process is not considered in other models, except for Yomogida & Matsui (1984) and Benoit & Sears (1996). The model of Yomogida & Matsui (1984) unfortunately does not consider the important contribution of ²⁶Al to the heating. Therefore this model cannot be compared to other models. Benoit & Sears (1996) gave only very few details on their model results such that one cannot make a comparison with our and other models also in this case. The models of Akridge et al. (1998) and Harrison & Grimm (2010) do not consider sintering, but simulate the insulating effect of a porous outer layer, that escaped sintering, over a consolidated body by introducing two outer layers with low filling factor, a “regolith” and a “mega-regolith” layer. The thickness’s of both layers were assumed by Akridge et al. (1998), but they were determined by an optimisation of the model by Harrison & Grimm (2010) and a thickness of 410 m was obtained for each of the layers. This compares well with a thickness of  <1 km of a residual porous outer layer in our model. The thermal shielding of the consolidated interior of the body by a thin insulating outer layer of low thermal conductivity is obviously essential for obtaining a good fit between thermal evolution models to meteoritic data as it is obtained in the model by Harrison & Grimm (2010) and in ours.

A major difference between our model results and others is the earlier derived formation time, t_form. This difference probably results from the lower mass-density used in other model calculations and the resulting lower mass of the body. Because of the resulting higher total heat capacity of the whole body (^∫ρc_pdV) in our model more heat from decay of radioactive nuclei is required to heat the body to the temperatures as witnessed by the observed cooling history of the various meteorites. Hence, an earlier formation time is required for our optimised model since the specific heating rate at zero time (equal to time of CAI formation) is almost the same in all models. The higher average density in our model is a result of the sintering process which results in very low or vanishing porosity over most of the volume of the body.

If (6) Hebe really is the parent body, the mass-density in our model fits well to that determined for Hebe (James & Chesley 2008; Baer et al. 2011). The densities used in other model calculations then would be too low and the earlier formation time found in the present model calculation therefore is probably more realistic than the later formation times found in other models.

The surface temperature, T_srf, is determined in our model by the optimisation process, in the other model it seems to be chosen more or less arbitrary within the frame expected for bodies in the asteroid belt and in some cases adjusted to improve the model fit between closing temperatures and closing times observed for meteorites. The dependence of the thermal evolution model on the value of T_srf is generally weak (e.g. Henke et al. 2012). The value determined by our optimisation procedure hinges mainly upon the fit to the two H4 chondrites that are most strongly affected by T_srf because of the shallow burial depths required by their petrologic type.

There are a number of other model calculations (Bennett & McSween 1996; Hevey & Sanders 2006; Sahijpal et al. 2007; Warren 2011) that also studied the thermal evolution of asteroids and give important insight into their properties and the processes operating in them, but did not attempt to obtain definite results for the properties of the parent body of H chondrites from some kind of “best fit” to meteoritic data. Only the range of possible values is studied in these papers. This does not enable a direct comparison, but the range of possible values for radii and formation times given are in accord with the results shown in Table 9.

5.9. Orbit of parent body

If the number obtained for T_srf from our optimisation process is taken seriously, it measures the average surface temperature for the first  ≈20...30 Ma of evolution where the temperature for Forest Vale and Ste. Marguerite was above the closure temperature for Pu fission tracks. Since the observed lifetime of accretion disks is less than 10 Ma (Haisch et al. 2001) the surface temperature is determined by illumination by the young proto-sun and thermal emission from the surface, in particular during the most critical phase where the temperature at the burial depth evolved through the closing temperature for Pu fission tracks. Then we can estimate the average radius of the orbit of the parent body from the relation (for a rapidly rotating body) (13)where σ is the Stefan-Boltzmann constant, A the albedo of the surface, a the distance to the protosun in AU, and S(1) the solar constant at 1 AU at this early evolution period.

For the solar constant S(1) we have to observe that the sun ignites at about 10 Ma hydrogen burning and has the luminosity of a 1 M_⊙ star on the main sequence. Its luminosity during the first 100 Ma is by about 30% less than its present value (e.g. D’Antonia & Mazzitelli 1994). The typical surface albedo for S asteroids is 0.2 (Ryan & Woodward 2010). The surface albedo of (6) Hebe is particular high with 0.291 ± 0.029 (Ryan & Woodward 2010). The albedo found today is probably not typical for the material that covered the surfaces of the bodies 4.5 Ga ago, because highly processed material and the products of its grinding are presented at the surface. It must had contained a lot of more fine-grained material that only later matured by crystal growth at elevated temperatures and probably was darker.

We assume arbitrarily A = 0.1 and find a radius of the orbit of a = 1.96 AU (an albedo of A = 0.2 yields a = 1.85 AU). This is smaller than the present average radius of the orbit of 2.4 AU. If the temperature T_srf found in optimisation is really representative for the average surface temperature during the essential phase of the thermal evolution, and if (6) Hebe is the parent body of the H chondrites, then our result means that (6) Hebe formed closer to sun than its present distance. However, a more detailed treatment of the processes determining the surface temperature is necessary, before any definite conclusions can be drawn.

6. Concluding remarks

We combined the code for modelling the internal structure and evolution of planetesimals of the 100 km size class (Henke et al. 2012) with a code for optimisation of a quality function that is based on an evolution algorithm. This allows finding a consistent set of parameters that describe the properties and evolution of a planetesimal and the burial depths of the meteorites. As such the thermal evolution of the model body reproduces very closely the empirically determined cooling history of all H chondrites with sufficient available data.

The results show, that the “onion shell” model successfully reproduces the thermal evolution history of all the meteorites used in this study. For six of the eight meteorites a burial depth can be found for which the empirical cooling data can be fitted in an excellent way with the local evolution of temperature at that burial depth. For the remaining two cases it is still possible to find a burial depth that results in a good fit if the Pu fission-track datum of Nadiabondi is rejected. This confirms earlier conclusions that the thermal evolution of the H chondrite parent body and its initial thermal structure has not been disturbed by catastrophic impacts.

The reconstructed properties of the parent body of the H chondrites are comparable to what has been found in other model calculations. The different published models treated the problem in different kinds of approximations, solving the heat conduction equation being almost the only feature common to all models, and used different values for the material properties. Our estimated size of 111 km of the parent body is somewhat bigger compared to the 100 km found in other models, and requires a slightly earlier formation time of 1.8 Ma compared to the 2.2–2.5 Ma found in other models. Such minor differences are certainly due to the different approaches used.

Our study shows, that the input physics of the models has to be improved, before more reliable data on the parent body can be obtained. Such improvements are currently underway in our group. In particular the heat conductivity has to be calculated from the specific properties of chondritic material. It will also be necessary to obtain more precise radiometric ages for a wider range of samples.

The advantage of our method is that it allows for a simultaneous and consistent fit of several model parameters to a big number of observational restrictions and is very flexible if the kind of data to be included in the optimization process is to be changed. A minor shortcoming of the optimisation method used is that a considerable number of complete evolution models has to be calculated during the course of the optimisation process, but our implicit solution method is fast enough that no restriction as to the applicability of the method to our problem results from this.

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)”. S.H. is Member of the IMPRS for Astronomy & Cosmic Physics at the University of Heidelberg.

References

All Tables

All Figures

Fig. 1 Schematic representation of the temperature evolution at two different depths (full lines) and at the centre (dotted line) of a planetesimal of the 100 km size class. The horizontal dashed lines correspond to fixed temperatures equal to the closure temperatures in Table 1. The circles indicate the instants at which at a given depth the temperature drops below a given temperature. These instants correspond to the closure ages of the different thermochronometers used for meteorites.

In the text

Fig. 2 Results for the individual meteorites used in the model fit. The values of the free parameters and the burial depth’s of the meteorites found for the optimised model are given in Table 5. The full line shows the temperature evolution at the burial depth’s of the meteorites in the asteroid, as determined by the evolution algorithm. The dashed line shows the temperature evolution at the centre and the dotted line is the surface temperature. The square boxes indicate the individual data points and the error bars the uncertainty of closing temperature and age determinations.

In the text

	Fig. 3 Results for a model fit where the heat conductivity at zero porosity, Kb, is included in the optimisation process. The values of the free parameters and the burial depth’s of the meteorites found for the optimised model are given in Table 6. For meaning of symbols see caption to Fig. 2.
In the text

	Fig. 4 Results for a model fit where the radius is set equal to the observed radius of (6) Hebe. The burial depth’s of the meteorites found for the optimised model are given in Table 7. For meaning of symbols see caption to Fig. 2.
In the text

Fig. 5 Evolution history of the value of χ2 during the optimisation process with progressing number of generations, n, for the model where the heat conductivity at zero porosity, Kb, is included in the optimisation process. Shown is for each generation the particular individual with the lowest value of χ2 of all individuals of that generation, the formation time, tform, of this individual, and the radius of the body, R. For two representative meteorites, Kernouvé and Allegan, the evolution of the optimum value of the burial depth, d, is also shown.

In the text

	Fig. 6 Data for thermal conductivity Kexp of ordinary chondrites (from Yomogida & Matsui 1983) divided by the analytical fit K(φ) given by Eq. (6) using a value of Kb = 1.75 obtained as optimum value for fitting the cooling history. Full black circles represent data for H chondrites, open circles data for L chondrites. The dashed line corresponds to the value of K = 4 used in model with fixed K.
In the text

	Fig. 7 Results for the individual meteorites if a Miyamoto-like model is used for the model fit. The burial depth’s of the meteorites found for this model are given in Table 8. For meaning of symbols see caption to Fig. 2.
In the text