Geo- and Cosmochemistry

DFG SPP 1385

The First 10 Million Years of the Solar System - a Planetary Materials Approach

Downloads/Results

1. Project: Modelling the internal constitution and thermal evolution of planetesimals and applications to geochronology (Hans-Peter Gail (U Heidelberg), Mario Trieloff (U Heidelberg), Thorsten Kleine (U Münster), Stephan Henke (U Heidelberg)) (Gepris: GA 162/4-1/2/3):

   Preface

   The numerical model presented below came into being during a project funded by the DFG during the SPP 1385 The first ten million years of the solar system and was develloped by Stephan Henke of the Institut für Theoretische Astrophysik of the Universität Heidelberg. The version published here is a programme which was edited for the use by others. The research that was done with this model is published in several papers:
   • Thermal evolution and sintering of chondritic planetesimals (Paper 1)
   • Thermal history modelling of the H chondrite parent body (Paper 2)
   • Thermal evolution model for the H chondrite asteroid - instantaneous formation versus protracted accretion (Paper 3)
   • Thermal evolution and sintering of chondritic planetesimals II. Improved treatment of the compaction process (Paper 4)

   The thermal evolution of asteroids model

   This programme simulates the thermal evolution of asteroids with sizes of the order of 100 km in the early solar system. Thermal evolution of asteroids describes the following process :
   • Heating of the asteroid by decay of mainly ²⁶Al (half-live 0.7 Ma) which was abundant enough in the early solar system to cause significant heating of the order of 1000 K in the first few Ma after CAI formation. Other isotopes or impacts did not produce significant amounts of heat for asteroids of the considered size class.
   • Cooling of the body by emission of heat to the environment which basically means black body radiation into space under consideration of the surface temperature which is given by the equilibrium of emission and irradiation from the sun.The protoplanetary disk is only present during the first 10 Ma of the solar system and thus of minor importance.
   • Processes that are active in the asteroid due to the change of temperature such as sintering, chemical reactions, change of aggregate state of the material components or escaping of volatiles.
   The programme offered here is specialised to simulate the thermal evolution of ordinary chondrite parent bodies. The physical processes considered here are:
   • Heating by decay of ²⁶Al, ⁶⁰Fe, ²³⁵U,²³⁸U,²³²Th, ⁴⁰K.
   • Sintering according to the theory of Kakar and Rao or according to the theory of Helle.
   • A temperature dependent heat capacity matching the material properties of ordinary chondrites.
   • A porosity and temperature dependent heat conductivity based on literature data.
   • The pressure within the body is calculated from the hydrostatic pressure equation.
   • The pressure dependent initial porosity if the initial value is set to higher than 40 %.
   • Melting of the Fe,Ni-FeS phase and of the silicate phase.
   • Optional: Linear radius growth for a predefined period of time.
   • Since there is no evidence for significant abundace of water or other volatiles in ordinary chondrites, processes such as freezing melting or vaporisation of water as well as volatile loss by venting and other causes are not considered in this model.
   The programme allows the user to easily modify the model parameters radius, formation time, initial porosity, surface temperature and bulk heat conductivity to model ordinary chondrite parent bodies with different porosities. Output can be produced for the temperature evolution in several different customisable depths.

   Mathematical features

   This model is considered to be a spherically symmetric problem. Therefore in this programme the spherically symmetric heat conduction equation is solved by a finite difference method assuming a von Neumann boundary in the center of the body and a fixed temperature at its surface.

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

   Optimisation

   We also provide a version of the code that can be used to determine model parameters of a parent body for which meteoritic cooling age data exist, as is the case for the H chondrite parent body. For that the code uses a genetical algorithm which is based on the principles of the evolution theory of Darwin. We use the genetical algorithm PIKAIA which can also be found here. We used the PIKAIA version 1.2 which we 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 rather basic number generators that however works fine with PIKAIA, since genetical algorithms are very robust to low quality random numbers. If you want a better random number generator, 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 you will find 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 user wants to change datasets frequently 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 model settings can be changed in the file Daten.dat. These settings are intended for every day use and contain only very few parameters. A lot more 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 of each spatial step within the modelled body as they are found after completion of the thermal evolution. Those parameters are 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 outputs the surface temperature of the asteroid, which is only needed for plotting.
   Also a gnuplot script Thermal_Evolution_Gnuplot is given, which produces a figure that shows the evolution of the temperature for 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 garantied 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 corresponding age for the meteorite, the third the error for the age and the last column states 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 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 devided.

   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 te block titled ParBereich. Given there are the lower and upper range of the parameter. If both values are the same, the parameter is not varied. Ranges can be set for formation time, radius, Fe-60/Fe-56 ratio, initial porsity, surface temperature, bulk heat conductivity and accretion duration.

   (By default surface temperature before disk dispersal is set to the temperature after. 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 formation times lie usually around 2 Ma after CAI-formation and radii are between 100 and 200 km.
   • The ratio of Fe-60/Fe-56 during CAI-formation can be varied. But since newer research showed that its value lies around the default value 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). 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 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.
   • Newer research by us (Paper 5) shows that the bulk heat conductivity likely hat the default value. But we also often used a range between 1 and 9 W/mK.
   • 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 to different possibilities of how the heat capacity is determined for the considered material which can be set by the logical variable cp_genau. If .false. a simple arithmetic fit to meteoritic data is used which is published in Yomogida/Matsui. If .true. heat capacity temperature data are read in from file Subroutinen/Cp_H.dat 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 that, if an H chondrite parent body is modelled. Furthermore also data files for L chondrite material and acapulco material are provided (Subroutinen/Cp_L.dat, Subroutinen/Cp_A.dat).
   • Two different sintering modes are supported and can be changed by the logical variable Helle. If .true. the sintering theory of Helle is used. If .false. the sintering theory of Kakar and Rao is used. We consider the theory of Helle to be the more accurate one which also was implemented later during our research activity. The theory of Kakar and Rao was used earlier.
   • We allow the user to assume 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. .false. stands for the one component mixture, .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. Only use Gemisch together with Helle.

   Download:

   Files