Astrea

A program for modelling the internal structure and evolution of planetesimals with ordinary chondrite composition and for finding optimal fits to parent body properties from cooling data of meteorites

The numerical model program was written by Stephan Henke at the Institut für Theoretische Astrophysik, Zentrum für Astronomie of the Universität Heidelberg as part of a project of the special research programme (SPP) 1385 The first ten million years of the solar system funded by the Deutsche Forschungsgemeinschaft (DFG). The program is written in FORTRAN90. This published version is a programme which was edited for public use. The basic physical assumptions on which the model calculation is based are described and the research that was done with this model program is published in the papers:

The thermal evolution of asteroids model

This programme simulates the thermal evolution of water free undifferentiated asteroids with sizes of the order of 100 km (for pressure effects in the equation of state to be negligible) in the early solar system. The model includes the following process :

Internal heating of the asteroid by decay of short and long lived radioactives.
External heating of the asteroid by solar irradiation and cooling of the body by emission of heat to the environment
Heat conduction within the body.
Processes that are active inside the asteroid due to the change of temperature such as
- sintering of the granular material, and
- melting (but no sedimentation)
of the material.

The programme is constructed to simulate the thermal evolution of ordinary chondrite parent bodies. The basic assumptions and physical processes considered are in detail:

Heating by decay of short and long living radioactives ²⁶Al, ⁶⁰Fe, ²³⁵U,²³⁸U,²³²Th, ⁴⁰K.
Sintering of granular media according to the theory of Helle.

Optional: Sintering according to the simple theory of Kakar and Rao.
Temperature dependent heat capacity for the material of ordinary chondrites.
Porosity and temperature dependence of heat conductivity based on literature data.
The pressure within the body is calculated from hydrostatic equilibrium.
The initial porosity after cold pressure compaction is calculated if the initial value is set to values higher than 40 %.
Melting of the Fe,Ni-FeS phase and of the silicate phase.
Spherical symmetry of the body is assumed.
Instantaneous Accretion is assumed.

Optional: Linear radius growth for a predefined period of time.
No melt migration is considered.
No pressure dependencies of solid state properties are considered.

The programme allows the user to easily modify the model parameters radius, formation time, initial porosity, surface temperature, bulk heat conductivity, and heat capacity to model ordinary chondrite parent bodies with different properties. Output can be produced for the temperature evolution in several different customisable depths.

Mathematical features

The spherically symmetric heat conduction equation is solved by a fully implicit finite difference method assuming a von Neumann boundary condition at the centre of the body and a fixed temperature at its surface.

The spatial steps are distributed in a logarithmic scale providing a higher resolution near the surface than at the centre since temperature changes are usually steeper near the surface. In addition the algorithm has a time step size control that checks the change in temperature, porosity (and mass if growth is considered) in each time step and adjusts the time step size according to accuracy requirements.

Optimisation

A version of the code is provided that can be used to determine optimized model parameters of a parent body for which sufficient meteoritic cooling age data exist, as is the case for the H chondrite parent body. For that purpose the code uses a genetical algorithm which is based on the principles of the evolution theory of Darwin. The special version of such an algorithm chosen for our purpose is the code PIKAIA, version 1.2, which is described here. The code was converted to Fortran 90 using the Alan Miller's tool to_f90.f90.

PIKAIA uses a random number generator. By default it uses the random number generator of Park and Miller which is a rather basic number generators that however works fine with PIKAIA, since genetical algorithms are very robust to low quality random numbers. If a better random number generator is required, one can use for example ran2 from the Numerical Recipes in Fortran 90 (Press et al. 1996).

How to use the code

This code package is provided in .zip format. After unzipping there are two source files in the main folder. The one called Asteroid.f90 calculates the thermal evolution of an asteroid. The other one called Evolution.f90 finds the best asteroid parameters which match given cooling ages and temperatures for several meteorites.

Both source files are tested to compile correctly with the GNU fortran compiler gfortran as well as the Intel Fortran Comiler ifort. But any compiler should work.

Modifying input files

In case that the datasets have to be frequently changed the input filenames and paths can be modified in the file files_used.dat found in the subfolder subroutinen. In this description the files are however called by their default name.

The Asteroid code

The program called Asteroid.f90 calculates the thermal evolution of an asteroid. The parameters of the body and some basic model settings can be changed in the file Daten.dat. These settings are intended for every day use and contain only very few parameters. Many additional parameter settings can be found in the subfolder Subroutinen in the file Daten.dat. The data in the former file will overwrite any in the latter file during program start. Furthermore for each parameter there is a default value defined in the source code to ensure that each parameter has a reasonable value in case a parameter has been forgotten to set in the data files by the user.

Asteroid.f90 produces several output files that can be found in the subfolder Ausgabe:

Asteroid.dat shows some properties at each radial grid point within the modelled body as they are found after completion of the thermal evolution. The file contains the radial position, density, cumulative mass, pressure and maximum temperature.
Furthermore a summary of the model settings is shown as well as the depths and maximum temperatures of the nine layers that are chosen in the input parameters.
Temperature.dat, Pressure.dat and Filling_density.dat list the temperature, pressure and filling density, respectively, as function of time for the centre and the nine layers that are chosen in the input parameters.
Also a file Initial_temperature.dat is provided which contains the surface temperature of the asteroid, which is only needed for plotting.

Also a gnuplot script Thermal_Evolution_Gnuplot is provided, which produces a figure that shows the evolution of the temperature at the centre and the nine layers that are chosen in the input parameters. It uses as input files Temperature.dat and Initial_temperature.dat. The gnuplot script is tested with gnuplot version 4.6. Unfortunately it is not guaranteed that the script works properly with other versions.

The Evolution code

The programme Evolution.f90 calculates thermal evolutions for many bodies with different model parameters to determine which model parameters reproduce the data stored in file Alter.f90 in the subfolder Subroutinen. These data consist of cooling ages and temperatures for nine different H chondrites that are assumed to originate in the same parent body. These data can be modified to correspond to any other meteorite parent body for which sufficient data exist. Which parameters are to be modelled and in which range is user defined (see Modifying data).

Parallel

This program can also be run in parallel mode using openmp for example by: ifort Evolution.f90 -openmp. By default the program will then use two cores. This can however be easily changed by changing the variable Kerne to the desired number of used cores. This variable can be found in Evolution.f90 in the module Parameter around line 80, where it is the first defined variable in this module.

Modifying data

Modifying the cooling age data used can be done in the file Subroutinen/Alter.dat. The first block in this file reads in the meteorite names for which data are available. These names are used throughout the program to identify the cooling ages and temperatures. The second block reads in the names of the radioactive decay systems used for the age dating. Note that the integer variables MetZahl and SchliessTempZahl found in Subroutinen/Daten.dat have to match the numbers of meteorite and decay system names respectively. Feel free to edit the numbers to your needs.

In the next MetZahl blocks the data for the meteorites is read in. The first column contains the name of the decay system, the second the corresponding age of the meteorite, the third the error of the age and the last column the temperature for the decay system. Furthermore for each meteorite an allowed range for the maximum temperature is defined (from Tmin to Tmax). First number is the temperature, second the error. The block named CoolAgeError contains the errors to the decay system temperatures. The reason for this separating is historical since during the code development the errors were assumed to be the same for each decay system for all meteorites but the temperatures themselves not. The small following block just contains the CAI-age used in the program. The last block named ParBereich contains the ranges of the model parameters in which the program is allowed to vary them.

The sequence of the blocks does not matter, however, the blocks by themselves are not allowed to be changed.

Reasonable parameter ranges

During the optimisation every parameter is varied within a certain range which is given by the user. Those ranges are also provided in Subroutinen/Alter.dat in the block titled ParBereich. The lower and upper range of the parameters have to be specified here. If both values are the same, the parameter is not varied. Ranges can be set for formation time, radius, ⁶⁰Fe/⁵⁶Fe ratio, initial porosity, surface temperature, bulk heat conductivity and length of accretion period.

(By default, surface temperature before disk dispersal is set to the temperature after dispersal. Changing this requires changes in Evolution.f90 in the function OptA at the rather top around line 295. There the assignment for Parameter_ein(8) has to be changed to that in the comment.)

For the H chondrite parent body, e.g., typical formation times found by models usually lie around 2 Ma after CAI-formation and radii are between 100 and 200 km. The upper and lower limits should be selected generous such the algorithm can check the likely range of values.
The formation time has to be given in units of Ma, the radius in units of km.
The ratio of ⁶⁰Fe/⁵⁶Fe at time of CAI-formation can be varied. It has to be given in units of 10^-6.
Since newer research showed that its value lies around the default value of 1.15×10^-8 it is recommend to use that one.
The initial porosity of the material is only varied for a one component material (see Different code modes). It has to be given in volume fractions (not in percent).
Findings in H chondrites constrain the initial porosity to be not significantly below 20%. Default value is set to 30 %. Do not chose porosities higher than 80% for sintering is not defined then.
The surface temperature is limited between 150 K, since below of this significant amounts of water ice are to be expected in the material which is not found, and 400 K, since the lowest closure temperatures lie at around 400 K.
The bulk heat conductivity has to be given in units of W/mK. Newer research by us (Henke et al. (2016))) shows that the bulk heat conductivity likely has the default value K = 4.3 W/mK. But we also often used a range between 1 and 9 W/mK.
Duration of accretion in units of Ma. In Paper 3 we showed that the instantaneous formation hypothesis (formation ca. 0.1 Ma) holds for the H chondrite parent body.

Different code modes

The code allows the user to select different modes of what physical principles are used:

There are two different possibilities of how the heat capacity is determined for the considered material which can be selected by the logical variable cp_genau. If c_p=.false. a simple arithmetic fit to meteoritic data is used which is published in Yomogida & Matsui (1983). If c_p=.true. heat capacity vs. temperature data are read in from file Subroutinen/Cp_H.dat which have been calculated for H chondrite material. This is a much more accurate treatment of the heat capacity for H chondrite material, so we recommend to use this mode, if a H chondrite parent body is modelled. Additional data files for L chondrite material and Acapulco material are provided (Subroutinen/Cp_L.dat, Subroutinen/Cp_A.dat).
Two different sintering modes are implemented and can be selected by the logical variable Helle. If .true. the sintering theory of Helle et al. (1985) is used. If .false. the sintering theory of Kakar & Chaklader (1967) is used. The second method corresponds to the method first applied in Yomogida & Matsui (1984). We consider the theory of Helle to be the more accurate one which was implemented later Henke et al. (2015) during our research activity.
The user can chose between a one component or a two component material consisting of spheres of only one size and of two different sizes respectively. The latter case would resemble a mixture of chondrules and matrix material. The mode can be selected by the logical variable Gemisch where .false. stands for the one component mixture and .true. for the two component mixture. If a two component mixture is chosen the ratio can be set by the variable f_ma which denotes the fraction of matrix material in the mixture. But only values between 0.0 - 0.265 and 0.444 - 1.0 are allowed (for the reasons given in appendix A of Henke et al. (2015)). Only use Gemisch=.true. together with Helle=.true..

Download:

Files