1. Introduction

Most meteorites are considered to be fragments of asteroids (McSween, 1999). A solar system body from which a meteorite originates is called the parent body of the meteorite. There have been a number of studies on the thermal evolution of the parent bodies of meteorites (e.g., Miyamoto et al., 1981; Grimm and McSween, 1989; Ghosh and McSween, 1998; LaTourrette and Wasserburg, 1998; Young et al., 1999; Ito and Ganguly, 2004). The decay energy of 26Al is the main heat source for such bodies (Urey, 1955). An analytic solution of the thermal history of a body heated by radioactive nuclei was firstly obtained and applied to the parent bodies of ordinary chondrites by Miyamoto et al. (1981), and some authors have applied this solution to determine the thermal histories of homogeneous bodies (LaTourrette and Wasserburg, 1998; Ito and Ganguly, 2004). On the other hand, thermal evolutions, which include ice melting and a hydrothermal reaction (Grimm and McSween, 1989) or rock melting and differentiation (Ghosh and McSween, 1998), are very complicated and are determined numerically. Grimm and McSween (1989) considered the thermal evolution of parent bodies for carbonaceous chondrites, containing ice (e.g. H2O). Key processes included the effects of ice melting and hydration reactions, and the thermal and chemical evolutions of parent bodies were numerically determined (Grimm and McSween, 1989; Young et al., 1999; Cohen and Coker, 2000; Young, 2001).

It is important to consider hydration reactions in bodies which initially contain ice. The temperature increases by the hydration reaction, as a result of which serpentine forms from forsterite and H2O. The hydration reaction causes the runaway melting of ice (Cohen and Coker, 2000). On the other hand, the increase in temperature is suppressed due to ice melting and hydrothermal convection. Previous studies on the thermal evolution of parent bodies of CM and CI chondrites performed numerical simulations including hydrothermal convection and gas diffusion (Grimm and McSween, 1989; Cohen and Coker, 2000); the hydrother-mal convection through the pores of rocks occurs when the Rayleigh number Ra exceeds (Turcotte and Schubert, 2002), and gas diffusion is controlled by Darcy’s law. Several authors have considered hydrothermal convections using three-dimensional models (Travis and Schubert, 2005; Palguta et al., 2010), but one study has reported an unlikeliness of large scale fluid flow (Bland et al., 2009).

The planetesimals in the early solar system would have had a wide range of initial size. Recent studies into the formation of planetesimals show that streaming instability accelerates dust concentration (Youdin and Goodman, 2005), and 1000-km- and 100-km-size planetesimals are formed by gravitational instability if the magneto-rotational instability is active (Johansen et al., 2007), and passive (Johansen et al., 2009), respectively. Other work has reported that 10–100-km-size planetesimals are formed by dust concentrations triggered by turbulence (Cuzzi et al., 2008; Chambers, 2010).

The thermal evolution of planetesimals of a wide range of sizes, including the effects of the hydration reaction, must therefore be considered. Previous studies, which include the hydration reaction, mainly performed numerical simulations of parent bodies having a 50-km radius, because this radius was considered to be the most appropriate for the parent bodies of CI and CM chondrites (e.g. Grimm and McSween, 1989; Cohen and Coker, 2000). The thermal evolutions of larger icy bodies (e.g. large satellites of outer planets and Kuiper belt objects) have been investigated by considering the rocky-core formation due to radiogenic heating (Schubert et al, 2007; Desch et al, 2009). However, these studies did not include the thermal effect of hydration reactions.

In this paper, a thermal evolution model of icy planetesi-mals is constructed, and the temperature evolutions are numerically solved with various values of parameters (size, initial abundance of 26Al and initial temperature). Many physical processes are involved, such as convection and rocky-core formation, as well as chemical processes, such as aqueous alteration and dehydration reaction. We take the effects of both chemical and physical processes into account to develop a thermal evolution model of icy planetesimals. We treat an icy planetesimal in the early solar nebula generally, located in a region where H2O almost completely condenses.

2. Methods

2.1 Compositions and parameters

A planetesimal in the early Solar System would have been made of elements which could condense at a distance from the Sun among the Solar System abundances of elements (Anders and Grevesse, 1989). The major elements and molecules of the solar nebula are shown in Table 1, in which the molecular abundances are calculated following Yamoto (2003). We assume all metals in rocks are oxidized. The residual oxygen is shared with a ratio of H2O/CO = 5, which is derived from the data of comet Hale-Bopp (Crovisier, 1998). We assume carbon atoms, except those in CO, are in hydrocarbons (CH2). The carbon content in the Tagish Lake meteorite is only 5.81 wt% and much lower than the solar abundance of elements (Grady et al , 2002). Although it is not clear whether the Tagish Lake meteorite represents icy planetesimals, hydrocarbons and carbon oxides (see below) are assumed to be in the gas phase, for simplicity. We assume that planetesimals do not contain any carbons. As seen in Table 1, the mass ratio of H2O/rocks (oxidized metal) is 1.72. We assume the initial rocks contain major refractory elements with the ratio Si: Fe : Mg = 1:1:1, which is given by the approximate ratios of the solar abundance, and a mineral assumed to be olivine ((Mg0.5,Fe0.5)2SiO4).

Table 1 Abundances of elements and molecules in the solar nebula. The number of atoms are normalized to Si (= 1.00 × 106).
Full size table

The values of the initial temperatures of icy planetesimals are determined so that almost all H2O molecules condense, but other molecules which may form ices (e.g. carbon oxides and nitrogen) do not. The minimum-mass solar nebula model (Hayashi, 1981) is used here. The temperature of the minimum-mass solar nebula (Hayashi, 1981) and the partial and saturated pressures of H2O (Murphy and Koop, 2005) are given by

(1)
(2)
(3)

respectively, where RAU is the distance from the Sun, (see Table 1), is the gas density in the minimum-mass solar nebula, is the Boltzmann constant, is the mean molecular weight and is the mass of the H atom. If the ratio of exceeds, say, one hundred, most H2O molecules would condense; such a temperature in the solar nebula is 130 K. Most of the CO molecules would be in the gas phase for (Aikawa et al., 1999). The minimum temperature for carbon oxides to avoid condensation is assumed to be 70 K. We consider several values of the initial temperature of icy planetesimals between 130 K and 70 K. These values correspond to 4 AU and 16 AU, respectively, in the solar nebula (Eq. (1)). We also assume that all condensed H2O is crystalline ice.

We take several values between 10 km and 1000 km for the planetesimals’ radius (Johansen et al., 2007; Cuzzi et al, 2008; Chambers, 2010). The maximum value, 1000 km, is obtained by Johansen et al (2007) by the numerical simulation of planetesimal formation due to gravitational instability, assuming that the magneto-rotational instability is active. Although the formation time of planetesimals is still unknown, we assume planetesimals are formed instantaneously when the condition for the formation of planetes-imals is satisfied at an orbital radius in the solar nebula. We assume a planetesimal has a constant radius and we do not take into account the growth of a planetesimal. If we consider the growth of a planetesimal, the amount of radionu-clide on the additional surface due to the growth would be lower than those in early-formed internal parts. This might lead to a lower temperature evolution on the surface, but this would not affect the center of the planetesimal as long as it is large enough. Compression due to gravity might be significant with a radius larger than 500 km (see p. 247 of Lewis, 1997). We omit the effects of gravitational compression for simplicity.

The main heat source is the decay energy of the shortlived nuclide 26Al, the half-life of which is 0.72 Myr. Another short-lived nuclide, 60Fe could also be a heat source for the planetesimals (e.g. Ito and Ganguly, 2004). We will discuss the effect of 60Fe in Section 4.4. The highest initial ratio of 26Al/27 Al in the solar system would be 5.0 × 10-5, the canonical value of Ca-Al-rich inclusions (CAIs) (e.g. MacPherson et al., 1995). CAIs are thought to be the earliest solids in the Solar System (4567 Ma Amelin et al., 2010). We take several values between 5.0 × 10-5 and 8.0 × 10-7 for the initial ratio of 26Al/27Al. We take the upper limit to be equal to the canonical value of CAI, which means a planetesimal formed just after the formation of the solar nebula. We choose the lower limit because the temperature does not reach the melting point of ice if the initial ratio is smaller than that. These values correspond to the time of the CAI formation and 4.3 Myr after the CAI formation, respectively.

2.2 Basic equations

We use the heat-conduction equation for a spherically-symmetric body to calculate the temperature evolutions. The equation is written as follows,

(4)

where is the density, c is the specific heat, T is the temperature, t is time, r is the radius from the center, K is the thermal conductivity, A is the initial radiogenic-heat-generation rate per unit volume, and is the decay constant of a radionuclide. The initial condition (at ) is where is the initial temperature, for which we take two values: 130 K and 70 K (see Section 2.1). The boundary condition at the center is and the boundary condition at the surface is assumed to be , i.e., the planetesimal has a uniform temperature identical to the surrounding nebula. When the melting of ice and the chemical reactions (described in Section 2.4) take place, the latent heat of ice melting and the chemical-reaction heats are taken into account. We solved Eq. (4) using a finite difference method, and an explicit method of time integration is used (see Section 2.5, where we discuss the time steps and the stability). The density, specific heat and thermal conductivity of each material (ice, water and rocks) are shown in Table 2, and all of these values are assumed to be constant. We assume that water and ice have the same density, for simplicity, and the values for rocks are based on Yomogida and Matsui (1983) for Leoville (CV3), and the values for ice and water are based on Murphy and Koop (2005) and Chronological Scientific Tables (2010). We use the mass-weighted mean specific heat, the volume mean density and thermal conductivity for a region where several materials are mixed. For example, an initial icy planetesimal contains ice and rocks (see Section 2.3), and the mean values are given by,

(5)
(6)

and

(7)

where represents the mass ratio of rocks to the total, represents the volume ratio of rocks to the total, the symbol “ -” denotes the mean value, and the subscripts i and r denote ice and rocks, respectively. The values obtained are and for the values of the parameters in Table 2.

Table 2 Density ρ, specific heat c, and thermal conductivity K of ice, water and rocks.
Full size table

The only radionuclide considered in this work is 26Al, the decay constant of which is 9.63 x 10-7 year-1.We take the thermal energy per26 Al decay (Edecay) as 3.16 MeV (Schramm et al., 1970), and the abundance of Al () from the CI carbonaceous chondrite (Lodders and Fegley, 1998). The initial heat-generation rate per unit volume is given by, where m26 is the mass of an26 Al nucleus. Note that the initial heat generation rate per unit volume depends on the initial ratio of 26Al/27Al, which we choose to be in the range 5.0 x 10-5 to 8.0 x 10-7 (see Section 2.1).

2.3 Core formation model of icy planetesimals

An initial icy planetesimal is made of ice and rocks; the latter would be small grains (≥ 1mm), chondrules and CAIs embedded among the ice. The melting temperature, Tmelt, of the ice depends on the pressure, but we assume that the melting temperature of the ice is 273 K, for simplicity. The latent heat of the ice melting is 3.34 x 105 J kg-1 (Chronological Scientific Tables, 2010). The internal ice of a planetesimal melts due to the radiogenic heat of26 Al, and rocks would sink toward the center. The planetesimal would then have a structure with a rocky core, a water mantle, an ice melting region and a frozen crust (see Fig. 1).

Fig. 1
figure 1

Illustrations of the evolution of icy planetesimals. First stage (top left) shows an icy planetesimal made of ice and rocks. The second stage (top right) of an icy planetesimal appears when the temperature reaches Tmelt of ice (273 K). At this stage, an icy planetesimal has a frozen crust, which retains the initial state, and an ice melting region, in which ice, water and rocks coexist at Tmelt of ice. The third stage (bottom) shows an icy planetesimal after the complete melting of ice and the settling of rocks have occurred, except at the surface regions. An icy planetesimal has 4 regions from the center to the surface: the rocky core, the water mantle (273 K), the ice melting region (273 K), and the frozen crust.

Full size image

It is a complicated problem to solve the convective motion of the mixture of rocks, water and ice and settling of rocks. We make an assumption that the rocks begin to settle after the ice is completely melted, and that an icy planetesimal has an initial mean density during the settling, i.e. . The mass of a rock m is , where is mean material density of the rock and a is the radius of the rock. If the drag force is given by the Stokes formula, the equation of motion for a settling rock in water is given as:

(8)

where is the viscosity of water, and G is the gravitational constant. Substituting the above expressions into Eq. (8), we obtain

(9)

The solution is given by , where r0 is the initial radius at which the rock begins to settle and is the settling time: . From the above, , for 1-mm-size rocks, is 7 × 106 s. Note this time scale decreases if we take into account the increase in density of the icy planetesimals due to the settling, instead of using in the expression of M(r) . In the following, we assume that a rocky core, and a water mantle, form instantaneously after the complete melting of ice, because is much shorter than the thermal evolution time scale of the order of 1 Myr. The heat flux from the rocky core only raises the temperature of the water mantle by, at most, several degrees above 273 K. Because the increase of the water temperature is small, we simply assume that the temperature of the water mantle is kept equal to 273 K. The region where the temperature could not reach Tmelt remains frozen with the initial composition, which we call the frozen crust.

2.4 Chemical reactions

As for the chemical reactions between the rocks and the water, the aqueous alteration and the dehydration reaction are considered. We assume the initial rocks are made of olivine (Mg0.5,Fe0.5)2SiO4 only, for simplicity (see Section 2.1). When ice melts due to 26Al heating, olivine and water react, and aqueous alteration takes place. In our model, we assume the rocky core and the water mantle form after complete ice melting, and that rocks and water could not coexist after the core formation. The ice melting region is the only region where the rocks and water could coexist, and the temperature of the ice melting regions is 273 K. The aqueous alteration occurs during a period of 15 days at 200°C (473 K) and 60 days at 150°C (423 K) (Jones and Brearley, 2006). We obtained the duration of the aqueous alteration at 273 K, calculating the activation temperature Tact based on their results. The reaction-rate coefficient is assumed to be given by . The ratio of product to the original reactant during the reaction time is written . We found that . By using these relations, we find that aqueous alteration occurs during a period of 210 yr at 273 K; this period is very short compared to the thermal evolution time scale (of the order of 1 Myr). We assume that aqueous alteration takes place at 273 K in the ice melting regions.

Serpentine , saponite and magnetite (Fe3O4) in CI chondrite are the main products due to aqueous alteration (Tomeoka and Buseck, 1988). In this study, we take into account only one chemical-reaction formula for the aqueous alteration:

The values of the iron ratio, Fe/(Mg+Fe) of serpentine and talc are determined based on Tomeoka and Buseck (1988).

We assume the dehydration reaction takes place at 873 K, based on experimental results (Nozaki et al , 2006). When the temperature reaches 873 K in the rocky core, serpentine and talc would be dehydrated, and dehydrated olivine, silica and water would be produced; on the other hand, magnetite could remain at 873 K. We make the assumption that the value of the iron ratio of dehydrated olivine is the same as serpentine and talc, for simplicity. The chemical reaction formula could be

H2O is produced in the dehydration reaction; however, the mass ratio of H2O/rocks is small (about 0.07). The effect of water after the dehydration reaction is not considered in this work. The products and reactants of rocks are summarized in Fig. 2. Due to the aqueous alteration, the initial olivine alters into serpentine, talc and magnetite. The dehydrated olivine and silica are produced from serpentine and talc by the dehydration reaction. If the temperature decreases after the dehydration reaction, which releases water around the dehydrated rocks, the aqueous alteration may occur again. However, for simplicity, we do not take this reaction into account.

Fig. 2
figure 2

Ternary diagram of rocks. Initial olivine (diamond in black), serpentine (down triangle in black), talc (up triangle in black), magnetite Fe3O4(square in white), dehydrated olivine (diamond in white) and silica SiO2(circle in white) are plotted.

Full size image

The chemical-reaction heats are included in the calculation and those are calculated from the enthalpy of each mineral. We take the enthalpies of minerals from Holland and Powell (1998). Calculated results show that the aqueous alteration is an exothermic reaction and its heat is 2. 77 × 105 Jkg-1, and the dehydration reaction is an endothermic reaction and its heat is 4.17 × 105 J kg-1. The dehydration reaction needs more heat than aqueous alteration produces.

The thermal effect of the aqueous alteration is described as follows. From the formula of the aqueous alteration, we found that water and an anhydrous mineral (olivine) are consumed with a mass ratio of 0.1. As a result, the mass of the hydrous minerals (serpentine, saponite and magnetite) is 1.1 times as large as the mass of the original anhydrous mineral. The heat production due to the aqueous alteration is nine times larger than the latent heat of the ice melting required to supply the water consumed by the aqueous alteration. Once the ice begins to melt and a small amount of water is generated, the aqueous alteration which takes place produces enough heat to melt more ice which, in turn, causes further aqueous alteration. Hence, once the ice begins to melt, all the anhydrous minerals are converted to hydrous minerals instantaneously. The heat produced by the aqueous-alteration reaction melts ice of mass nine times as large as the mass of water consumed by the aqueous alteration, as mentioned above. Therefore, just after the instantaneous reaction, water appears. The mass of the water appears after the aqueous alteration is eight times as large as the mass of water consumed in the aqueous-alteration reaction.

2.5 Numerical methods

We solved the basic equation (Eq. (4)) as follows. Hereafter we write , which represents the time at the n-th step and position of j-th mesh point. The first term of the RHS of Eq. (4) is approximated by using a finite-difference method with . We have

(10)

for , and

(11)

for j = 0. Equation (11) is derived as follows. The Taylor expansion of T r around r = 0is

(12)

and the limit of leads

(13)

On the other hand, we have

(14)

From Eqs. (13) and (14), we find Eq. (11).

In order to perform the time integration, we use the explicit method. From Eq. (4), for , we can write

(15)

From Eqs. (4) and (11), for j = 0, we can write

(16)

The analytic solution of Eq. (4) is

(17)

Equation (17) is the temperature profile of bodies with the Dirichlet boundary condition (Carslaw and Jaeger, 1959).

We compare our numerical solution to the analytic solution using the error defined by

(18)

For the stability of numerical simulations, the Courant-Friedrichs-Lewy conditions should be fulfilled. The CFL condition is given by . We examined our numerical code, which satisfies the CFL condition, for a 500-km-size rocky planetesimal. As seen from Figs. 3 and 4, the solution obtained by our numerical code is in good agreement with the analytic solution. We used for all our numerical simulations shown in Section 3 (e.g., for planetesimals with a radius of 500 km, yr and km).

Fig. 3
figure 3

Time evolutions of (defined by Eq. (18)) for a constant time step yr and various values of km (dashed line), km (dot-dashed line), and km (solid line). The body investigated is a rocky planetesimal with a radius of 500 km.

Full size image
Fig. 4
figure 4

Time evolutions of (defined by Eq. (18)) for the constant value of km and various values of yr (dotted line), yr (solid line), and yr (dashed line). The investigated body is a rocky planetesimal with a radius of 500 km.

Full size image

When the temperature T becomes slightly above the melting temperature of ice, , as a result of numerical integration, the ice partially melts and is converted to water, keeping the temperature equal to Tmelt. Then, the aqueous-alteration reaction begins to occur, and the reaction proceeds as a sort of chain reaction (see Section 2.4). Numerically, the following changes take place at the time step at which T exceeds Tmelt: (1) All the anhydrous rocks are converted to hydrous rocks. The mass of the hydrous rocks is 1.1 times as large as the original anhydrous rocks, i.e., where is the mass of the hydrous rocks among a unit mass of the mixture of hydrous rocks, ice and water. (2) The mass of ice is reduced by 0.9 times the mass of the original anhydrous rocks, i.e., where is the mass of ice among a unit mass of the mixture (ice, rock and water). (3) The mass of water is 0.8 times the mass of the original anhydrous rocks, i.e., , where is the mass of water among a unit mass of the mixture. Note that before aqueous alteration, , and after aqueous alteration, , and the law of mass conservation is satisfied.

After the completion of the aqueous alteration, the remaining ice begins to melt to form water without being consumed by the aqueous-alteration reaction, keeping the temperature equal to Tmelt. During this stage, the changes of and in one time step are given by

(19)

When reaches 0, the ice has totally melted. Then, the hydrous rocks are assumed to fall onto the central rocky core instantaneously. In the rocky core, we solve Eq. (4) by using parameters for rocks instead of mean values of the ice and rock mixture.

When the rocky core reaches the temperature for the dehydration reaction, , we calculate the change of the mass of hydrated rock in one time step by , where H is the en-dothermic heat for the dehydration reaction, keeping (see Section 2.4). When reaches 0.3, the serpentine and talc are totally dehydrated. On the other hand, magnetite, the mass ratio of which is 0.3, remains unchanged at 873 K (see Section 2.4). After that, we solve Eq. (4) by using parameters for rocks again.

At the lower boundary of the frozen crust, which contacts with the water mantle, the heat flux in the frozen crust is supplied by the latent heat of water condensation in the water mantle. The thickness of water which condenses during a time step just below the boundary is given by . Because in our calculation, the fraction of the partial freezing is memorized and used in the calculation in the next time step. In the case where the water mantle contacts the ice melting region, a calculation similar to the freezing of the water mantle is carried out including the heating by 26Al simultaneously.

Although there is the heat flux from the inside of the rocky core at the outer boundary of the rocky core, the temperature change in the water mantle due to the flux is estimated to be, at most, several degrees of kelvin, and we just assume that the temperature of the water mantle is kept equal to 273 K.

3. Results

We solved Eq. (4) by numerical simulation; the results are shown in the following. The used model parameters of numerical simulations, and brief results, are shown in Tables 3 and 4, where is the time after CAI formation in Myr, tice is the duration of ice melting in Myr, tdehydration is the duration of the dehydration reaction in Myr, and Tpeak is the peak temperature at the center in kelvin. As seen from Tables 3 and 4, we find that icy planetesimals formed 2.4 Myr after CAI formation could not reach the melting temperature of ice. In these models, it is required that icy planetesimals are formed within 2.4 Myr after CAI formation to melt the ice.

Table 3 Model parameters of numerical simulation with Tinit of 130 K. Note that “—” means the corresponding event does not occur.
Full size table

3.1 Radius

Figure 5 shows the temperature evolutions at the centers of different-sized icy planetesimals with radii 10 km, 50 km, 100 km, 500 km and 1000 km which formed 1.0 Myr after CAI formation with 130 K for the initial temperature (models 10–10A, 50–10A, 100–10A, 500–10A and 1000–10A, respectively). Icy planetesimals over 50 km in radius experience ice melting, aqueous alteration and the dehydration reaction. The duration of ice melting beginning at 1.3 Myr and ending at 1.5 Myr after CAI formation, and the dehy-dration reaction beginning at 2.3 Myr and ending at 3.3 Myr after CAI formation are the same for icy planetesimals having radii over 50 km. This means that thermal conduction is not of effect for these sizes. By way of contrast, the 10-km-size icy planetesimal (model 10–10A) could not attain the melting temperature of ice. Thus, thermal conduction is needed for a 10-km-size icy planetesimal formed within 0.5 Myr after CAI formation to reach the melting temperature of ice (see Tables 3 and 4).

Fig. 5
figure 5

Temperature evolutions at the centers of planetesimals for models 10–10A (dot-dot-dashed line), 50-10A (dashed line), 100–10A (dot-dashedline), 500–10A (solidline) and 1000–10B (thick dotted line). The abscissa shows the time after CAI formation in units of a million years and the ordinate shows the temperature at the center in kelvin. The horizontal dotted lines represent the temperature at which the chemical reactions take place: 273 K (ice melting and aqueous alteration) and 873 K (the dehydration reaction). The results of 500–10A and 1000–10A are almost the same and those lines overlap.

Full size image
Table 4 Model parameters of numerical simulation with Tinit of 70 K. Note that“—” means the corresponding event does not occur.
Full size table

3.2 Initial temperature

Temperature evolutions at the centers of different-sized icy planetesimals, with radii 10 km, 50 km, 100 km, 500 km and 1000 km, which formed 1.0 Myr after CAI formation with 70 K for the initial temperature (models 10–10B, 50-10B, 100–10B, 500–10B and 1000–10B, respectively), are shown in Fig. 6. Compared to Fig. 5, only the initial temperature is different. Icy planetesimals over 50 km show a dehydration reaction, but could not exceed the dehydration temperature (873 K). The timing and duration of ice melting are also different, because the energy, which is required to reach the ice melting temperature and dehydration temperature, for low initial temperature models (70 K) is larger than that for high temperature models (130 K). Comparing the results between Tables 3 and 4, we note that the initial temperature also affects the starting point and the duration of melting and dehydration reaction.

Fig. 6
figure 6

Temperature evolutions at the centers of planetesimals for models 10–10B (dot-dot-dashed line), 50–10B (dashed line), 100–10B (dot-dashed line), 500–10B (solid line) and 1000–10B (thick dotted line). The abscissa, the ordinate and the horizontal dotted lines are the same as in Fig. 5. The results of 500–10B and 1000–10B are almost same and those lines overlap.

Full size image

3.3 Formation time

The results for 500-km-size icy planetesimals, which formed 1.0 Myr, 1.5 Myr, 2.0 Myr and 2.4 Myr after CAI formation with 130 K for the initial temperature (models 500–10A, 500–15A, 500–20A and 500–24A, respectively), are shown in Fig. 7. The formation time (i.e. the initial ratio of 26Al/27Al) largely affects the temperature evolution, i.e. the occurrence of the ice melting, the aqueous alteration and the dehydration reaction, as well as the beginning and the duration of chemical reactions.

Fig. 7
figure 7

Temperature evolutions at the centers of planetesimals for models 500–10A (solid line), 500–15A (dotted line), 500–20A (dot-dashed line), and 500–24A (dashed line). The abscissa, the ordinate and the horizontal dotted lines are the same as in Fig. 5.

Full size image

The details of an icy planetesimal having a 500-km radius, and 130 K for the initial temperature, formed at 1.0 Myr after CAI formation (model 500–10A), are shown in Figs. 8 to 11. Before the ice melting begins (e.g., 1.2 Myr in Fig. 9), the temperature is almost constant with radius, with the exception of a thin surface region. After the ice melting, there are two regions where the temperature is almost constant, i.e. the rocky core except for a thin surface region of the rocky core, and the water mantle; on the other hand, the frozen crust has a larger temperature gradient (see Fig. 9). A rocky core and a water mantle are formed at 1.5 Myr after CAI formation. The radii of the rocky core (270 km) and the thickness of the water mantle (220 km) are almost the same over 10 Myr (Fig. 8).

Fig. 8
figure 8

The temperature evolution of model 500–10A. The abscissa shows the time after CAI formation in units of a million years and the ordinate shows the radius from the center in km. The boundaries of a rocky core, a water mantle and the melting region (see Fig. 1) are shown (thick dashed line). The isothermal lines (thin solid line) are drawn every 100 K.

Full size image
Fig. 9
figure 9

The temperature evolution of model 500–10A. Line graphs represent 1.2 Myr (dashed line), 2.0 Myr (dot-dashed line), 2.8 Myr (dot-dot-dashed line) and 6.0 Myr after CAI formation (solid line), respectively. The abscissa shows the radius from the center in km and theordinate shows the temperature in kelvin. The horizontal dotted lines are the same as in Fig. 5.

Full size image

Figure 10 shows the mass ratio of ice, water and rocks at 6.0 Myr after CAI formation. Due to the dehydration reaction, the formation of dehydrated rocks (olivine and silica) occurs in the rocky core (260 km) except for a thin outer shell (10 km). Cumulative energies of heat sources (radiogenic heat and heat of aqueous alteration) and heat consumptions (change in internal energy, latent heat of ice, heat of dehydration reaction and radiation at the surface) are shown in Fig. 11. The heat of the aqueous alteration is released in a short period of 0.007 Myr around 1.3 Myr after CAI formation. On the other hand, the consumption of heat by the dehydration reaction takes longer than the aqueous alteration (2.3 Myr-3.3 Myr after CAI formation). The total energy (difference between heat source and heat consumption) is kept almost equal to zero, which means that the energy conservation is well satisfied in the numerical simulation.

Fig. 10
figure 10

The mass ratios of materials. Mass ratios of ice (thin solid line), water (thin dashed line), initial rocks (thick dot-dashed line), altered rocks (thick dashed right) and dehydrated rocks (thick dot-dot-dashed line) at 6.0 Myr after CAI formation of model 500–10A (solid line in Fig. 9) are shown. The abscissa shows the time after CAI formation in units of a million years and the ordinate shows the mass ratio. The right panel shows a magnified figure between 490 km and 500 km.

Full size image
Fig. 11
figure 11

Cumulative energy of an icy planetesimal for model 500–10A. The lines represent radiogenic heat (solid line), heat of the aqueous alteration (long-dashed line), change in internal energy (dot-dot-dashed line), latent heat of ice (dotted line), radiation at the surface (short-dashed line), heatofthedehydrationreaction(dot-dashed line) and “heat source-heat consumption”, i.e. the difference between the sum of the heat sources and the sum of the heat consumptions (thin solid line), respectively. The abscissa shows the time after CAI formation in units of a million years and the ordinate shows the cumulative energy in joules. Note that the heat sources (radiogenic heat and heat of the aqueous alteration) are positive in the ordinate, and that the heat consumptions (change in internal energy, latent heat of ice, radiation at surface, heat of the dehydration reaction) are negative in the ordinate.

Full size image

4. Discussions

4.1 Rocky core

The temperature evolutions at the centers of icy planetes-imals, having radii over 50 km, after the dehydration reaction are different for different radii (Fig. 5). We examine the radii of rocky cores to see the cause of this difference. The radii of rocky cores of icy planetesimals for 50 km, 100 km, 500 km and 1000 km (models 50–10A, 100–10A, 500–10A and 1000–10A, respectively) are 25 km, 53 km, 270 km and 550 km, respectively. It is expected that the temperature evolution in a rocky core would evolve in the same way as a rocky planetesimal, which has a radius of rocky core, formed 1.5 Myr after CAI formation (i.e. the time when the rocky core formed) with 273 K for the initial temperature. The temperature evolutions at the center of rocky planetesimals are shown in Fig. 12. The temperature evolutions at the centers of an icy planetesimal, and that of a rocky planetesimal, are nearly the same provided that the radius of the rocky core of the icy planetesimal is the same as the radius of the rocky planetesimal. The difference in size determines the rocky core’s radius as long as the time of formation is the same, and the radius of the rocky core, in turn, determines the temperature evolution of the rocky core.

Fig. 12
figure 12

Temperature evolutions at the centers of different-sized icy plan-etesimals and rocky planetesimals. Results for model 50–10A (thin dashed line), 100–10A (thin dot-dashed line) and 500–10A (thin solid line), and for rocky planetesimals with radii of 23 km (thick dashed line), 53 km (thick dot-dashed line), and 270 km (thick solid line) formed 1.5 Myr after CAI formation with 273 K for the initial temperature are shown. The abscissa, the ordinate and the horizontal dotted lines are the same as in Fig. 5.

Full size image

4.2 The effect of chemical reactions

In order to elucidate the effects of chemical reactions (the aqueous alteration and the dehydration reaction), we carried out a numerical simulation of temperature evolution without chemical reactions. We used the same parameters of model 500–10A (500-km-size icy planetesimal formed 1.0 Myr after CAI formation with 130 K for the initial temperature). We compare these results with, and without, chemical reactions in Fig. 13. As seen from Fig. 13, the ice melting begins at 1.3 Myr and ends at 1.8 Myr after CAI formation; this period is longer than model 500–10A with the aqueous alteration. Aqueous alteration plays an important role in ice melting (Cohen and Coker, 2000).

Fig. 13
figure 13

Temperature evolutions at the centers of planetesimals for model 500–10A, without chemical reactions (dot-dashed line), and with chemical reactions (solid line). The abscissa, the ordinate and the horizontal dotted lines are the same as in Fig. 5.

Full size image

4.3 Effect of porosity

Meteorites are composed of porous material (Yomogida and Matsui, 1983). The parent body of meteorites (planetesimals or asteroids) would also be porous. It is, therefore, worthwhile considering the effect of porosity. Thermal conductivity is reduced in porous material because of the small contact areas of the grains (Haruyama et al , 1993). However, a porous material enables vapor flow before a melting temperature is achieved, and this contributes to heat transfer. It is difficult to determine whether the effective thermal conductivity decreases or increases.

Here, we simply calculate two cases, in which the thermal conductivity is twice and a half, respectively, our standard value (see Section 2.2), in order to gauge the effect of the change of thermal conductivity in the frozen crust. We assume the pores are filled with water after the melting of ice. The effect of porosity will appear before ice melting and does not show after the melting of ice starts. Figure 14 shows the result of icy planetesimals having a 50-km radius and 130 K for the initial temperature formed 1.0 Myr after CAI formation (model 50–10A) with twice, half of, and just equal to, the standard value for the thermal conductivity in the frozen crust. The ice melting region changes as a result of the effect of porosity. The radius of the rocky core and the thickness of the water mantle in an icy planetesimal are greater when we halve the thermal conductivity, but smaller when we double the thermal conductivity. However, the beginning, and the duration, of ice melting does not change. This means that the effect of porosity in the frozen crust does not affect the timing and duration of these events. If a rocky core is porous, the core would contain water in the pores; in this case, the hydrothermal reaction might continue after the core formation. We leave this issue for future work.

Fig. 14
figure 14

Comparison of the radii of the inner and outer boundaries of water mantles (see Fig. 1) for three values of thermal conductivity: twice (dot lines), half (dot-dashed lines), and just equal to (solid lines), the standard value. Each icy planetesimal has the same parameters as model 50–10A. The abscissa and the ordinate are the same as in Fig. 8.

Full size image

4.4 60 Fe as a heat source

Some authors take into account both26 Al and 60Fe as heat sources for the parent bodies of meteorites. Ito and Ganguly (2004) used a Solar System initial ratio of 60Fe/56Fe ∼ 3 × 10−7 (Tachibana and Huss, 2003). Recent studies show a Solar System initial ratio of 60Fe/56Fe < 1 × 10-8 (Spivak-Brindorf et al., 2011; Tang and Dauphas, 2011) and < 2 × 10−7 (Telus et al., 2011). It is worthwhile considering the effects of 60Fe using recent values of the initial ratio of 60Fe/56Fe.

The heat production rate of 26Al and 60Fe is shown in Fig. 15. At first, the heat production rate of26Al is larger than that of 60Fe. But the heat production rate of 60Fe exceeds that of26 Al at 7 Myr for 60Fe/56Fe = 2.0 × 10−7, and at 11 Myr for 60Fe/56Fe = 1.0 × 10−8. This is because 60Fe has a longer half-life than26Al (0.72 Myr; see Section 2.1). The half-life of 60Fe is 2.6 Myr (Rugel et al., 2009).

Fig. 15
figure 15

The heat production rate of 26Al and 60Fe. Each line represents the heat production rate by 26Al (without 60Fe, solid line) and 60Fe (60Fe/56Fe = 2.0 × 10−7 (dot-dashed line) and 60Fe/56Fe = 1.0 × 10−8(dotted line)). The abscissa shows the time after CAI formation in units of a million years and the ordinate shows the heat production rate in watt per kilogram.

Full size image

We have calculated the thermal evolution of icy plan-etesimals with 26Al and 60Fe as heat sources. We take the thermal energy released by the decay of one 60Fe nucleus to be 2.71 MeV (Tuli, 2003), and the abundance of Fe from CI carbonaceous chondrite (Lodders and Fegley, 1998). Figure 16 shows the result for icy planetesimals having radii of 100 km and 50 km and 130 K for the initial temperature formed 1.0 Myr after CAI formation (models 100–10A and 50–10A) for three values of 60Fe/56Fe: 1.5 × 10−7,0.8 × 10−8 and 0 (without 60Fe). The values of 60Fe/56Fe at 1.0 Myr after CAI formation are calculated from the Solar System initial ratio of 60Fe/56Fe (2.0 × 10−7 and 1.0 × 10−8). In these calculations, we include aluminum heating with 26Al/27Al = 2× 10−5 at 1.0 Myr after CAI formation. The early evolution before the beginning of the dehydration reaction is not affected by the abundance of 60Fe, because 26Al heating prevails, as shown in Fig. 15. In the results for icy planetesimals with 60Fe/56Fe = 1.5 × 10−7, the peak temperature is higher and the duration of the dehydration reaction is shorter than in the results for icy plan-etesimals without 60Fe. However, the temperature profile does not change in the results for icy planetesimals with 60Fe/56Fe = 0.8 × 10−8. A large Solar System initial ratio of 60Fe/56Fe is needed to affect the temperature evolution of icy planetesimals. Further studies on the initial ratio of 60Fe/56Fe are needed in order to consider 60Fe as a heat source.

Fig. 16
figure 16

Temperature evolutions at the centers of icy planetesimals with the additional heat source 60Fe (60Fe/56Fe = 1.5 × 10−7(dot-dashed line) and 60Fe/56Fe = 0.8 × 10−8(dotted line)) and without 60Fe (solid line). Results for models of 100-km radius (thick line) and 50 km radius (thin line) formed at 1.0 Myr after CAI formation with 150K for the initial temperature are shown. The abscissa, the ordinate and the horizontal dotted lines are the same as in Fig. 5. Results with 60Fe/56Fe = 0.8 × 10−8, and without 60Fe, are almost the same and those lines overlap.

Full size image

4.5 Comparison of the results to the chronological data of meteorites

The planetesimals with H2O ice, but without CO ice, are located at 4 AU-16 AU in the early solar nebula, where most of the water molecules condense, and have an ice/rock mass ratio, of 1.72 (see Section 2.1). It is not well known if there exist meteorites which originate from such icy plan-etesimals.

We have compared the timing of aqueous alteration in carbonaceous chondrite with our results. The timing of aqueous alteration is obtained by 53Mn-53Cr age studies. The half-life of 53Mn is 3.7 Myr. From recent isotopic studies for aqueously-altered carbonaceous chondrites, the initial ratios of 53Mn/55Mn are given as (2.3 ± 0.1) × 10−6 for Mokoia (CV3) (Hutcheon et al, 1998), (2.7 ± 1.0) × 10−6 for Allende (CV3) (Trinquier et al, 2008), (2.3 ± 0.5) × 10−6 for Vigarano (CV3) (Jogo et al, 2009), (2.7 ± 0.8) × 10−6 for Murchison (CM2.5), and (2.8 ± 0.4) × 10−6 for ALH 83100 (CM2.1) (Fujiya et al, 2010). The initial ratio of 53Mn/55Mn of CAIs is (9.1 ± 1.7) × 10−6 and the estimated absolute age of these CAIs is 4568 Ma (Nyquist et al , 2009). If we use these CAIs as a time marker, aqueous alteration occurred at 6 Myr-7 Myr after CAIs. Our calculation shows that icy planetesimals formed 2.0 Myr after CAI formation at 4 AU having radii over 50 km (models 50– 20A, 100–20A, 500–20A and 1000–20A) satisfy those times of aqueous alteration. However, our model assumes that the initial water/rock mass ratio is equal to 1.72, which corresponds to the solar nebula abundance, and we also assume that all the rocks have once suffered aqueous alteration at 273 K. We should take into account different ice/rock ratios, as well as hydrothermal reactions in porous rock cores, in order to simulate the thermal evolutions of the parent bodies of carbonaceous chondrite.

5. Conclusions

The thermal history and time evolution of planetesimals in the early solar nebula is one of the important issues for the elucidation of planet formation. We have performed numerical simulations for the thermal and time evolutions of icy planetesimals to constrain on its size and time of its formation under certain conditions.

Icy planetesimals, which are larger than 50 km in radius and which formed 1.0 Myr after CAI formation with 130 K for an initial temperature, could experience ice melting, aqueous alteration and a dehydration reaction. The duration of the ice melting (1.3–1.5 Myr after CAI formation) and the dehydration reaction (2.3–3.3 Myr after CAI formation) are the same for icy planetesimals having radii greater than 50 km. Thermal conduction is not of effect for these sizes. The radius of the rocky core depends on the icy plan-etesimal radius, which, in turn, affects the temperature evolution of the rocky core. On the other hand, an icy planetesi-mal 10 km in radius, formed at 1.0 Myr after CAI formation with 130 K for the initial temperature, could not attain the melting temperature of the ice. Icy planetesimals formed at 2.4 Myr after CAI formation, or later, could not attain the melting temperature of the ice either, regardless of their size and initial temperature.

The beginning, and the duration, of the ice melting and the chemical reactions are the same for different sizes of icy planetesimals formed at the same time, as long as their radius is greater than 50 km, but are different for the same-size icy planetesimals formed at different times. The difference in the formation time of an icy planetesimal is the most influential parameter for the occurrence of ice melting, aqueous alteration and the dehydration reaction, as well as the beginning, and duration, of chemical reactions. The effects of chemical reactions are also considered. Aqueous alteration has the effect to melt ice quickly.

The growth of icy planetesimals, the differentiation of the rocky core, and the hydrothermal reaction due to water contained in the pores of the rocky core, are aspects which should be considered in future work.