Differentiation of Vesta and the parent bodies of other achondrites
Abstract
[1] Numerical simulations have been performed to understand the planetary differentiation of Vesta and other differentiated asteroids. The short-lived nuclides 26Al and 60Fe were used as the heat sources for the differentiation of asteroids of sizes 20–270 km. The essential aim was to simulate the two proposed differentiation scenarios associated with the formation of the basaltic achondrites. One of these scenarios deals with their origin by partial melting of silicates. The other scenario involves the origin of basalts from the residual melt left subsequent to crystallization in a convective and cooling magma ocean. The core-mantle differentiation in both the scenarios was commenced subsequent to 40% silicate melting. The partial melt scenario seems to contradict the general trend in the chronological sequence of core-mantle and mantle-crust differentiations observed in various differentiated meteorites unless the irons meteorite parent bodies differentiated rapidly compared to prolonged accretion and differentiation of distinct parent bodies of achondrites. Even though the geochemical records of siderophile elements in eucrites favor the residual melt scenario, based on our simulations it is difficult to rapidly and globally cool the magma oceans to produce residual basaltic melts subsequent to crystallization within initial ∼10 million years. This is in contradiction with the chronological records of the basalts that suggest their early origin. In order to make this scenario viable by invoking rapid and global cooling of the magma ocean it would be essential to substantially remove the insulating crust of the planetesimal by continuous bombardment of small planetesimals in the early solar system.
1. Introduction
[2] The planetary differentiation of the differentiated meteorite parent bodies [Haack and McCoy, 2005; McCoy et al., 2006] commenced within the initial few million years of the early solar system. The chronometric records based on 26Al-26Mg, 53Mn-53Cr, 182Hf-182W, and 207Pb-206Pb isotopic systematic of the majority of the differentiated meteorites are consistent with the differentiation time-spans inferred from the numerical simulation of the differentiation of planetesimals [Hevey and Sanders, 2006; Sahijpal et al., 2007; Kleine et al., 2009]. The core-mantle differentiation in the case of iron meteorite parent bodies probably initiated along with the early accretion of the planetesimals, perhaps within the initial million years after the formation of Ca-Al-rich inclusions. [e.g., Kleine et al., 2002, 2005, 2009]. The mantle-crust differentiation in general seems to have initiated afterwards [e.g., Nyquist et al., 2003, 2009; Wadhwa et al., 2009]. A recent compilation of the published whole-rock ages and mineral isochrone ages of the basaltic meteorites [Sanders and Scott, 2007; Kleine et al., 2009; Nyquist et al., 2009; Wadhwa et al., 2009] suggest that the basaltic achondrites formed mainly within 4–10 million years (Ma) in the early solar system, with some exceptions.
[3] The major collection of the achondrites, the Howardite-Eucrite-Diogenite (HED) meteorites is considered to have been derived from a single asteroid Vesta [Binzel and Xu, 1993; Drake, 2001; Usui and McSween, 2007; Scott et al., 2009]. Recent studies, however, indicate that some of the eucrites could have also originated on other planetesimals [Yamaguchi et al., 2002, 2006; Mittlefehldt, 2005]. Based on the elemental composition, the eucrites are classified into four main groups: the main group eucrites, the Stannern group, the Nuevo Laredo group, and the high Mg eucrites [Reid et al., 1979]. Broadly, two distinct planetary differentiation scenarios have been proposed for the origin of the eucrites [Drake, 2001; Holzheid and Palme, 2007]. One of these scenarios argue that the eucrites are the result of distinct degree of equilibrium partial melting of chondritic material [Consolmagno and Drake, 1977; Stolper, 1977; Jurewicz et al., 1991, 1993, 1995; Jones et al., 1996]. Stolper [1977] showed that the main group eucrites (e.g., Sioux County) could be produced by 15–20% partial melting of the olivine, pyroxene, plagioclase, Cr-rich spinel, and metal. Consolmagno and Drake [1977] showed that the main group eucrites (e.g., Juvinas, Sioux County) could have been derived from a source with chondritic REE abundances by 10–15% equilibrium partial melting. Jurewicz et al. [1991, 1993] showed that eucritic magma could be produced by the partial melting of anhydrous Murchison (CM) under reducing conditions.
[4] The alternative source of the eucrites was proposed to be the residual melt left at the end of fractional crystallization in a cooling magma ocean [Ikeda and Takeda, 1985; Warren, 1985; Longhi and Pan, 1988]. The depletion of siderophile elements in eucrites provides evidence that the metal had segregated into the iron core before the formation of eucrites. Newsom [1985] estimated that the depletion of siderophile element Mo in eucrites is consistent with the core formation prior to the formation of eucrites [Hewins and Newsom, 1988]. Presently, the most favored magma ocean model, proposed by Righter and Drake [1997], provides several advantages over the fractional crystallization and the partial melting model. According to this model, the residual melt left subsequent to equilibrium crystallization, at 80–85% crystallization, in a magma ocean would produce the observed composition of the main group of eucrites. Subsequent fractional crystallization of the isolated main group eucrite melt could have produced the Nuevo Laredo group and cumulate eucrites. The Stannern group eucrites have been recently proposed to have been produced by partial crustal melting above the magma ocean [Barrat et al., 2007].
[5] Apart from the eucrites, the angrites constitutes another group of basaltic achondrites. Angrites are considered to be the most primitive group of achondrites that probably formed from the partial melting on a differentiated body of primitive refractory composition [Taylor et al., 1993; Mikouchi et al., 2008; Weiss et al., 2008]. Jurewicz et al. [1991, 1993] demonstrated that the composition of the Angra dos Reis can be explained by a 20% partial melting of a parent body with CV3 (Allende) chondritic composition. Besides the angrites and eucrites, two distinct groups of achondrites, namely, lodranites and ureilites [Mittlefehldt, 2007; Downes et al., 2007] are considered to be the leftovers residues of a low degree of partial melting and basaltic melt extrusion [Taylor et al., 1993]. The acapulcoites constitutes the extreme class of achondrites that have undergone the least amount of partial melting to result in any appreciable amount of melt segregation from the rock [Taylor et al., 1993].
[6] In general, in order to produce the wide-range of achondrites among the meterorite collection, the partial melting and/or the equilibrium/fractional crystallization of a large number of planetesimals in the early solar system would be required. The radiogenic decay energy of the short-lived radionuclide was proposed as a possible heat source for the heating and the melting of the parent bodies of these achondrites [Urey, 1955]. The evidence for the presence of once live 26Al in CAIs (Calcium-Aluminum rich inclusions) from Allende meteorite was found by Lee et al. [1976]. Wide-spread distribution of 26Al was experimentally confirmed in CAIs from different classes of chondritic meteorites [e.g., MacPherson et al., 1995]. A canonical value of ∼5 × 10−5 for 26Al/27Al in the early solar system was measured by the refinements in the analytic techniques [Jacobsen et al., 2008]. However, it should be mentioned that the canonical value of 26Al/27Al inferred from CAIs does not rule out the possibility of heterogeneous distribution of 26Al in the early solar system. There could be 26Al-free nebula reservoirs in the early solar system on spatial or temporal scales [Sahijpal and Soni, 2006]. The possible implications of the evolutionary history of asteroids formed in the 26Al-free nebula reservoirs is discussed by Sunshine et al. [2008]
[7] The recent evidence for the presence of high excesses of 60Ni in the grains of trolite in Semarkona meteorite [Mostefaoui et al., 2005] suggested that 60Fe was also responsible as a heat source for the differentiation of the planetesimals along with 26Al. The initial 60Fe/56Fe of the solar system is still uncertain. The estimated values of the initial 60Fe/56Fe are both high >10−6 [Moynier et al., 2005; Quitte et al., 2005] and low ∼(1–5) × 10−7 [Tachibana and Huss, 2003; Regelous et al., 2008]. However, it has been indicated that the extent of heterogeneity of 60Fe was less than 10% in the formation region of the planetary bodies [Dauphas et al., 2008]. We can anticipate an almost identical upper limit for the heterogeneity of 26Al in case the two short-lived nuclides 26Al and 60Fe were produced by a single massive star [Sahijpal and Gupta, 2009].
[8] Several thermal models of asteroids have been developed in order to explain their thermal evolution [Sahijpal, 1997, 2006; Sahijpal et al., 1995, 2007; Ghosh and McSween, 1998; Merk et al., 2002; Hevey and Sanders, 2006; Gupta and Sahijpal, 2007, 2009]. Sahijpal [1997] developed thermal models for the radiogenic heating of asteroids with the temperature dependent thermal diffusivity. Ghosh and McSween [1998] developed a model with instantaneous accretion of planetesimal. Hevey and Sanders [2006] modeled the melting of planetesimals heated by 26Al using finite difference method and incorporated thermal convection. Sahijpal [2006] and Sahijpal et al. [2007] have recently numerically simulated the planetary differentiation of linearly accreting planetesimals of radii in the range of 20–270 km, with the gradual growth of an iron core due to the flow of (Fe-Ni)metal-FeS melt toward the center of the planetesimals. This comprehensive work deals with the numerical simulation of the core-mantle differentiation in asteroids during the melting of Fe-FeS. The mantle-crust differentiation in these models was numerically initiated subsequent to an appreciable (∼20%) amount of silicate melting. As discussed elaborately by Sahijpal et al. [2007], and recently reviewed by Kleine et al. [2009], there is a general agreement among the results obtained from the thermal models developed by Merk et al. [2002], Hevey and Sanders [2006], Sahijpal et al. [2007], and Kleine et al. [2009].
[9] In the present work, we have numerically simulated two distinct differentiation scenarios that are relevant for our understanding of the origin of achondrites as elaborately discussed above. The primary concern is to numerically simulate the planetary differentiation of planetesimals with the origin of basaltic achondrites explained by two distinct scenarios. These scenarios involve (1) the partial melt origin of basalts, and (2) the origin of basalts from the residual melts of equilibrium crystallization in a cooling magma ocean. We make an attempt to understand the thermal evolution and differentiation in Vesta and other differentiated asteroids within these two distinct scenarios. The preliminary results of this work were presented in two conferences [Gupta and Sahijpal, 2007, 2009]. The present work is mostly focused on the initial H chondrite composition of the planetesimals as the melting of these chondrites are well understood both theoretically and experimentally [Taylor et al., 1993; McCoy et al., 2006]. In order to understand the influence of the initial bulk composition [Jurewicz et al., 1991, 1993, 1995; Taylor et al., 1993; Righter and Drake, 1997] on differentiation, we have also considered six different sets of composition.
2. Methodology
[10] The numerical simulations performed in the present work involve solving the heat conduction partial differential equation of a planetesimal undergoing a linear accretion growth. The growth resulted in the formation of a porous planetesimal that subsequently suffered thermal sintering and volume loss on account of compaction. The core-mantle and the mantle-crust differentiations in the body were performed gradually according to the thermal evolution of the body. The differentiation was executed by the inward and outward flows of the melt plumes of Fe-FeS and basalts.
2.1. Heat Conduction Equation
2.2. Simulation Parameters
[12] We have considered a linear growth in the radii [Merk et al., 2002] of the unconsolidated (porous) planetesimals of sizes 26 km, 65 km, 130 km, and 351 km from the nebular dust of 55% porosity. The equation that governs the linear accretion is R(t) = Ro + αt, where Ro is the initial radius of the unconsolidated planetesimal considered at the time of the initiation of the accretion of the planetesimal. It was taken to be ∼390 m. This would eventually be reduced to 300 m on account of compaction (consolidation) during sintering. The accretion growth of the planetesimals was performed by modifying the spatial grid array of the finite difference code. This was achieved by the addition of un-sintered (porous) 390 m shells to the pre-existing planetesimal, thereby resulting in the gradual growth of the planetesimals in an equal incremental size step [Sahijpal et al., 2007]. Since the accretion of the planetesimal was considered to be occurring in an unconsolidated (porous) state of matter we did not execute any thermal re-equilibration of the pre-existing planetesimal surface with the newly accreted shell as the thermal diffusivity of the unconsolidated body is assumed to be three orders of magnitude less than that of the consolidated body [Yomogida and Matsui, 1984]. The accretion of the body initiates at a specific time interval tOnset (0–3 Ma) after the formation of the Ca-Al-rich inclusions (CAIs) with the canonical value of 5 × 10−5 for 26Al/27Al [MacPherson et al., 1995]. The accretion duration (tDuration) of the planetesimals in the present work was considered in the range of 0.001 to 1 Ma for different sizes of the planetesimals [Sahijpal et al., 2007]. The lower time interval corresponds to almost instantaneous accretion. During the growth of the body the average temperature of the accreting bodies on the planetesimals were assumed to be identical to the ambient temperature of the nebula that is considered to be 250 K. A constant surface temperature of 250 K was maintained throughout the simulation [Hevey and Sanders, 2006]. A moving boundary condition during the accretion growth of planetesimals was adopted by redefining the planetesimal surface each time the spatial grid array was enlarged [Sahijpal et al., 2007].
[13] The majority of the simulations were performed with H chondrite composition. The uniformly distributed abundances of 1.22% and 27.8% were assumed for the Al and Fe, respectively, in the undifferentiated planetesimals. We have carried out most of the simulations with a value of 2 × 10−6 for (60Fe/56Fe)initial ratio at the time of the formation of CAIs with the canonical value of 26Al/27Al. Due to the uncertainty in defining the (60Fe/56Fe)initial ratio, we have carried out additional simulations with the (60Fe/56Fe)initial ratio of 5 × 10−7 and 1 × 10−6. The simulation parameters, e.g., the decay energies of 26Al and 60Fe, the temperature dependent specific heat and thermal diffusivity, the latent heat of the [(Fe-Ni)metal-FeS] and silicate melting, the densities of silicate, [(Fe-Ni)metal-FeS], the solidus and liquidus temperature range of silicate and [(Fe-Ni)metal-FeS], etc., were adopted from Sahijpal et al. [2007, Table 1].
[15] The porous planetesimals of radii 26 km, 65 km, 130 km, and 351 km acquire their final radii of 20 km, 50 km, 100 km, and 270 km, respectively, subsequent to sintering. We adopted the final radii of the planetesimals to represent the nomenclature of the simulations.
2.3. Differentiation of the Planetesimals
[16] We have considered the melting of [(Fe-Ni)metal-FeS] and silicates in the temperature range of 1213–1233 K and 1450–1850 K, respectively. A nonlinear temperature dependence of the melt fraction of the bulk silicate generated at a particular temperature [Taylor et al., 1993] was considered in the present work instead of the assumed simplified linear temperature dependence considered earlier by Sahijpal et al. [2007]. The latent heats of the melting of 2.7 × 105 J kg−1 and 4.0 × 105 J kg−1 were incorporated into the specific heat during the solidus-liquidus temperature range of [(Fe-Ni)metal-FeS] and silicate, respectively, according to the method adopted by Merk et al. [2002]. The specific heat of ∼2000 J kg−1K−1 was assumed for the silicate and [(Fe-Ni)metal-FeS] melts. The specific heat at a particular spatial grid interval was estimated by taking the weighted average of the specific heat of the melt and the solid mass fractions [Sahijpal et al., 2007]. The density of [(Fe-Ni)metal-FeS] was assumed to be 7800 Kg m−3. The schematic of the algorithm of the numerical simulation of the planetesimal accretion, followed by sintering and initiation of melting is presented in Figure 1.
[17] Numerical and experimental studies of the iron melt segregation from silicate infer complexity that is at present difficult to quantitatively incorporate in the numerical simulations of planetary differentiation [Yoshino et al., 2003; Rushmer et al., 2005; Petford et al., 2006; McCoy et al., 2006]. Sahijpal et al. [2007] performed detailed numerical simulations of the planetary differentiation scenario involving the initiation of the iron core formation prior to the bulk silicate melting. The core-mantle differentiation in the present work was triggered subsequent to 40% bulk silicate melting in accordance with the depletion of siderophile elements in eucrite parent body [Newsom, 1985; Hewins and Newsom, 1988; Drake, 2001]. The segregation of iron melt at substantial silicate melting is compatible with the numerical and experimental studies of iron melt segregation. We estimated the mass of [(Fe-Ni)metal-FeS] and silicates in both melt and un-melt forms in each 0.3 km wide spatial grid interval adopted for the finite difference method. Subsequent to 40% melting of the bulk silicates, the dense [(Fe-Ni)metal-FeS] melt plumes within a specific spatial grid were gradually moved toward the planetesimal center at an assumed rate of one spatial grid (i.e., 300 m) transfer of the melt per temporal grid interval of one year. This assumed rate of descend was numerically achieved by taking into account mass balance calculations. The schematic of the algorithm involved in the planetary differentiation is presented in Figure 2. Within a specific spatial grid interval, the iron melt plumes arriving from the upper layer were appropriately added to the pre-existing iron melt plumes in the spatial grid interval. The net momentum inflow of iron melt depends upon two factors, namely, the velocity of the inflow and the bulk mass transfer. We assumed a velocity of 300 m per year on the basis of our choice of the spatial and temporal grid intervals. This could be varied depending upon the choice. However, as mentioned in the beginning of this section we preferred the present combination of the spatial and temporal grid intervals to avoid complexities originating from the presently employed finite difference method with the classical explicit approximation. The present technique allows for a limited range of spatial-temporal grid size combinations. However, the use of Crank-Nicholson approximation [Lapidus and Pinder, 1982] in solving the finite difference method would relax the constraints on the choice of temporal grid interval for a specific spatial grid interval as the approximation leads to unconditionally stable solutions. We could develop an elementary numerical code using Crank-Nicholson approximation to decipher the thermal profile of the planetesimals but it is difficult for us to incorporate the planetary differentiation processes in this code in contrary to our code employing classical explicit approximation. Any future attempts to incorporate planetary differentiation processes in the finite difference method with the Crank-Nicholson approximation would allow the choice of any desired descend velocity of the iron melt plumes in the differentiating planetesimals. However, it should be mentioned that the present choice of spatial and temporal grid sizes are appropriate to simulate the core-mantle differentiation in a realistic manner. This aspect is discussed in details afterwards in this section.
[18] During the segregation of iron and silicate, a set of free parameters were used to transfer a fraction of the bulk accumulated mass of iron melt from a specific spatial grid interval to the next spatial grid interval, toward the planetesimal center. These parameters are essential to control the mass flow of the melt in a systematic manner by avoiding any excess of melt accumulation/depletion in any spatial grid. The procedure was executed by monitoring the net density of matter in the melt and un-melt forms for all the spatial grids.
[19] An identical approach was also followed for the ascend of silicate melt or un-melt pockets coming from the interior of the planetesimal. The specific heat, thermal diffusivity, the elemental abundances of Al and Fe were estimated within all the spatial grid intervals at each temporal step. The dense iron melt plumes, subsequent to reaching the central spatial grid, replace the silicates and push the silicate pockets in an upward direction due to buoyancy. This process was executed by using the mass balance calculations and logical statements in the numerical code. The [(Fe-Ni)metal-FeS] melt upon reaching the center initiated the formation of the dense iron-core. The differentiation was executed from inside-out, i.e., starting from the inner most spatial grid, the consecutive spatial grid intervals acquired [(Fe-Ni)metal-FeS] saturation. The gradual descend of the various Fe-FeS melt plumes from different spatial grids toward the center of the planetesimal eventually led to the formation of a dense [(Fe-Ni)metal-FeS] core. The replaced silicate that was pushed upward constituted the aggregate mantle and crust composition. In order to avoid discontinuity in any thermodynamical quantity and clumping in mass and density within any spatial grid interval the movement of the various pockets were monitored and controlled by a suitable set of free parameters as discussed above. These free parameters control the rate of descend and ascend of the various melt and un-melt pockets. An identical set of numerical values were chosen for these parameter in all the simulations spanning over the planetesimal accretion time-spans of 0.001 to 1 Ma. Along with the major elemental redistribution during the core-mantle and the mantle-crust differentiation, the redistribution of the heat sources 26Al and 60Fe was also implemented in the simulations. The silicate matrix melt in all simulations with the H chondrite composition retained ∼8.8% Fe as FeO along with a proportionate amount of 60Fe after core-mantle differentiation.
[20] It should be noted that the numerical criteria adopted for the planetary differentiation by Sahijpal et al. [2007], and in the present work is a phenomenological approach as it is extremely difficult to pursue the problem starting from the first principle. The phenomenology is the only opted scientific approach in circumstances where it is difficult either to work-out a theory starting from the first principle, or the associated mathematical complexities are so huge that a parameter based model, like the one opted in the present work, seems to be a better alternative. It should be noted that some of the most fundamental concepts in physical sciences till date are phenomenological, e.g., the Newton laws, numerous concepts in electrodynamics, high-energy physics, cosmology, and even fluid dynamics. Further, the most fundamental Darcy's law in the fluid dynamics, relevant in the present case [Taylor et al., 1993; Yoshino et al., 2003], is in fact a phenomenological approach. In the presently developed numerical code it is easier to move the fluids across an asteroid in a parameterized manner. The fluid flow rates adopted in the present work and by Sahijpal et al. [2007] are in fact not distinct from the rates deduced by the previous fluid dynamics model [Yoshino et al., 2003]. Even though we can parametrically control the fluid flow rates to an extent as discussed previously, yet it would be difficult to precisely reproduce the thermal evolution of planetesimals that involves a complex dependence of the fluid flow rates as a function of fluid viscosity, porosity, and the size of the differentiating planetesimals [Taylor et al., 1993; Yoshino et al., 2003]. However, the timescales associated with the core-mantle differentiation, quoted with the least count of 10 kilo-years in the present work, would not be significantly influenced by the differentiation process that occurs on the timescales of ∼10 kilo-years subsequent to the onset of the differentiation. These timescales are within the range of the migration timescales of metal melts based on Darcy's law [Yoshino et al., 2003]. However, the uncertainties in the melt percolation velocities of the basaltic melt subsequent to partial melting [Taylor et al., 1993] could significantly influence the deduced timescales of the mantle-crust differentiation. The uncertainties in the timescale involved in the basaltic melt migration to the planetesimal surface could be of the order of the mean-life of 26Al. We have appropriately incorporated this issue in our simulations.
[21] As mentioned earlier the temperature dependence of thermal diffusivity and specific heat [Sahijpal, 1997; Sahijpal et al., 2007] were considered in the simulations. The thermal diffusivity of the molten [(Fe-Ni)metal-FeS] was assumed to be 5 × 10−6 m2 s−1. The influence of convection as an efficient mode of heat transfer was incorporated in an analogous manner by increasing the thermal diffusivity of the [(Fe-Ni)metal-FeS] melt core and the molten silicate magma ocean by three orders of magnitude (i.e., 5 × 10−4 m2 s−1) compared to the sintered rock. The convection is a complex process [Marsh, 1989; Hort et al., 1999; Elkins-Tanton et al., 2008], and it is difficult at present to incorporate in the model [Hevey and Sanders, 2006; Sahijpal et al., 2007]. In the presently used finite difference method with the classical explicit approximation it is not possible to raise the thermal diffusivity beyond the assumed value of 5 × 10−4 m2 s−1 for the magma ocean as it causes numerical instabilities [Lapidus and Pinder, 1982]. The Crank-Nicholson approximation would be better in handling such instabilities without compromising the use of still higher thermal diffusivity [Lapidus and Pinder, 1982]. Nonetheless, it should be mentioned that the results obtained by Sahijpal et al. [2007] using classical explicit approximations matches exactly with those obtained by Hevey and Sanders [2006] using Crank-Nicholson approximation. Hence, we do not anticipate major influence of the use of an even higher thermal diffusivity on the outcome of the present work. We were able to achieve almost uniform temperature within the spatial extent of the magma ocean.
[22] The immediate concern regarding the convection in the present work is the behavior of the cooling convective region at the base of the outer crustal layer of the planetesimal. This is the region through which the planetesimal losses thermal energy. The convection shows strong dependence on the fluid viscosity that in turn depends upon the temperature. A variation in the viscosity by more than an order of magnitude can bring down the extent of convection [Marsh, 1989; Hort et al., 1999]. The resulting margin of the magma ocean that constitutes rigid crust, mush and suspension would provide an insulating blanket to the cooling magma ocean. This thermal insulation along with the insulation provided by any survived chondritic crustal layer would eventually control the cooling and the crystallization in the magma ocean. In this present work it is not possible to numerically simulate the behavior of a realistic magma ocean by considering all these aspects. Nonetheless, the survived chondritic crustal layer of thickness 2–3 km in the present setup controls the magma ocean cooling.
[23] In order to numerically implement the high thermal diffusivity of the molten magma ocean and the molten [(Fe-Ni)metal-FeS] core, the thermal diffusivity was gradually increased subsequent to 0.5 fraction of bulk silicate melting [Hevey and Sanders, 2006] by using a three parameter Sigmoidal function in the assumed temperature range of 1680–1820 K. A distinct set of the constants were adopted for equation (2). These are A = 5 × 10−4 m2 s−1, B = 1725 K, and C = 10 K, with κi as the thermal diffusivity of the sintered body. The temperature range from 1680 to 1820 K corresponds to 0.44 to 0.71 fraction of bulk silicate melting. As mentioned earlier the use of the Sigmoidal function was made to avoid singularity or discontinuity in any simulation parameter. The identical but reverse trend was followed during the cooling of the magma ocean, i.e., the thermal diffusivity was gradually decreased from 5 × 10−4 m2 s−1 to ∼5 × 10−7 m2 s−1 as the temperature dropped from 1820 to 1680 K upon cooling of the magma ocean. Another set of constants for the Sigmoidal function was adopted for equation (2) during the cooling. These are A = 5 × 10−4 m2 s−1, B = 1600 K, and C = 20 K, with κi as the temperature dependent thermal diffusivity of the sintered body. The variations in these parameters were found to have no major influence on the outcome of the present work.
[24] In the present work, we have modeled two distinct differentiation scenarios for the origin of basaltic achondrites in order to understand the differentiation of the Vesta and other differentiated asteroids. One of these scenario deals with the origin of the basalts during the partial melting of the silicates, and the other deals with the residual melt origin of the basaltic achondrites subsequent to crystallization in the cooling magma ocean on the planetesimals.
2.4. Partial Melt Origin of the Basaltic Achondrites: PM Model
[25] In this case, the mantle-crust differentiation occurs prior to the core-mantle differentiation due to partial melting of the bulk silicate that commences at 1450 K. In these simulations, the [(Fe-Ni)metal-FeS] melt moves toward the center at 0.4 fraction of bulk silicate melting and gradually forms the core. The 26Al-rich basaltic melt plumes generated in the different regions of the planetesimals were gradually moved upward toward the surface at 0.2 fraction of bulk silicate melting [Stolper, 1977; Jurewicz et al., 1991, 1993, 1995; Jones et al., 1996]. The systematic of the algorithm for this scenario is shown in Figure 3. The migration of the 26Al-rich basaltic melt plumes generated from different spatial grid intervals was executed distinctly. A set of logical statements were used to execute the systematic upward movement of the plumes. The radioactive heating due to the transit of the various basaltic melt plumes through various spatial grid intervals were appropriately considered. In order to address the issue of the uncertainties of the basaltic melt percolation velocities subsequent to partial melting, we performed two set of simulations with distinct velocities of the basaltic melt migration. Since the uncertainties in the timescales could be of the order of the mean-life of 26Al [Taylor et al., 1993], the basaltic melt plumes in the two set of simulations were moved at the rate of 300 m year−1 and 0.3 m year−1, respectively. The plumes moved outwards by one spatial interval of 300 m in a time step of 1 and 1000 year, respectively, for the two sets of simulations.
2.5. Residual Melt Origin of the Basaltic Achondrites: RM Model
[26] The [(Fe-Ni)metal-FeS] melt generated in this scenario moves toward the center to form the core at 0.4 fraction of bulk silicate melting. Further in this scenario, we have simulated the convective magma ocean beyond 0.5 fraction of bulk silicate melting [Righter and Drake, 1997; Drake, 2001]. We maintained the convection until the temperature dropped gradually below ∼1600 K, where at least appreciable extent of equilibrium/fractional crystallization is complete to produce a residual melt for basalt [Righter and Drake, 1997; Drake, 2001]. Even though we have not numerically simulated the gravitational settling of crystals in the cooling molten magma subsequently followed by basaltic volcanism, the essential features of this scenario can be well understood based on the present simulations. The schematic of the algorithm of this scenario is presented in Figure 4.
[27] Finally, in order to explore the dependence of the differentiation on the initial bulk composition of the planetesimals [Jurewicz et al., 1991, 1993, 1995; Taylor et al., 1993; Righter and Drake, 1997] apart from the presently used H chondritic composition we have considered six different sets of the initial bulk compositions. These are L, LL, CI, and CV chondritic compositions [Jarosewich, 1990], and the two binary combinations of compositions; one with 70% of LL and 30% of CI compositions and the other with 70% of L and 30% of CV compositions [Righter and Drake, 1997]. The temperature dependence of the thermal diffusivity and the specific heat in these simulations were also included [Yomogida and Matsui, 1983]. We assumed an identical set of the remaining thermodynamical functions that were used in the case of H chondrites.
3. Results
[28] The numerical simulations of the planetary differentiation of Vesta and other differentiated asteroids have been performed in the present work. The choice of the simulation parameters, specifically the time of the initiation of the planetesimal accretion (tOnset) subsequent to the formation of CAIs, the accretion duration (tDuration) and the (60Fe/56Fe)initial was made to substantially explore the parametric space in order to understand the extent of melting and differentiation of planetesimals within distinct scenarios. In general, a rapid accretion rate was chosen for small bodies of radii 20–50 km, whereas the larger bodies of radii 100 and 270 km were accreted gradually over a time-span of one million years.
[29] The temporal growth of the [(Fe-Ni)metal-FeS] core and the extrusion of the basaltic melt during the partial melting of silicate have been presented for a representative set of simulations in Table 1. A representative set of simulations corresponding to the residual melt origin of the basalts have been also considered in the table. The nomenclature of a specific simulation was chosen according to the nature of the adopted differentiation scenario and the general attributes of the planetesimal. We have used the model name PM (Partial Melt scenario) or RM (Residual Melt scenario), followed by the radius of the sintered planetesimals (20–270 km), the onset time of the accretion of the planetesimals (0–3 Ma), the accretion time of the planetesimals (0.001–1 Ma) and the initial ratio of 56Fe/60Fe in the range 5 × 10−7–2 × 10−6. In the case of the partial melt scenario two set of simulations were run with the basaltic melt percolation velocities of 300 m year−1 and 0.3 m year−1. The simulations with the slower melt percolation velocity are marked by the suffix ‘s’ with the radius of the sintered planetesimal. The details of the various simulation parameters are also provided in the footnote of the Table 1. The thermal profiles along with the growth of the iron core have been graphically presented in Figure 5 for the partial melt (PM) scenario and Figure 6 for the residual melt (RM) scenario.
| Simulationsa | Figure References | Fe-FeS Coreb | Initiation of the Extrusion of Basaltic Meltc | ||||
|---|---|---|---|---|---|---|---|
| Radius = 20 km | 0+ km | 2 km | 4 km | 6 km | 8 km | ||
| PM20-1-0.001–2(−6) | 1.73 | 1.74 | 1.79 | 2.09 | none | 1.35 | |
| PM20-2-0.001–2(−6) | 5a | none | none | none | none | none | 2.95 |
| PM20-1-0.1–2(−6) | 1.76 | 1.78 | 1.83 | 2.13 | none | 1.37 | |
| PM20-2-0.1–2(−6) | none | none | none | none | none | 2.98 | |
| Radius = 50 km | 0+ km | 5 km | 10 km | 15 km | 20 km | ||
| PM50-2-0.001–2(−6) | 3.84 | 3.84 | 3.85 | 4.07 | 5.08 | 2.95 | |
| PM50-1-0.1–2(−6) | 1.74 | 1.75 | 1.79 | 1.83 | 2.10 | 1.36 | |
| PM50-2-0.1–2(−6) | 5b | 3.88 | 3.90 | 3.97 | 4.35 | 5.37 | 2.96 |
| PM50-1-1-2(−6) | 5c | 1.89 | 2.04 | 2.38 | 2.53 | 2.69 | 1.39 |
| PM50-2-1-2(−6) | 4.30 | 4.59 | 5.88 | 7.69 | none | 3.04 | |
| Radius = 100 km | 0+ km | 10 km | 20 km | 30 km | 40 km | ||
| PM100-1-0.1–2(−6) | 1.74 | 1.77 | 1.79 | 1.83 | 1.88 | 1.35 | |
| RM100-1-0.1–2(−6) | 1.46 | 1.48 | 1.51 | 1.54 | 1.58 | - | |
| PM100-2-0.1–2(−6) | 3.85 | 3.89 | 3.96 | 4.03 | 4.94 | 2.96 | |
| RM100-2-0.1–2(−6) | 3.34 | 3.37 | 3.44 | 3.49 | 3.65 | - | |
| PM100-1-1-2(−6) | 5d | 1.80 | 2.04 | 2.38 | 2.59 | 2.76 | 1.38 |
| PM100s-1-1-2(−6) | 1.79 | 2.02 | 2.37 | 2.57 | 2.79 | 1.37(1.55) | |
| RM100-1-1-2(−6) | 1.48 | 1.63 | 1.84 | 2.05 | 2.31 | - | |
| PM100-1-1-1(−6) | 5e | 2.43 | 2.71 | 3.26 | 3.81 | 4.60 | 1.43 |
| PM100-1-1-5(−7) | 5f | 4.50 | 5.11 | none | none | none | 1.48 |
| PM100-2-1-2(−6) | 5g | 4.05 | 4.43 | 5.47 | 6.42 | 7.32 | 2.99 |
| PM100s-2-1-2(−6) | 5h | 4.03 | 4.42 | 5.45 | 6.33 | 8.17 | 3.00 (3.33) |
| RM100-2-1-2(−6) | 3.43 | 3.75 | 4.10 | 4.42 | 4.99 | - | |
| Radius = 270 km | 0+ km | 27 km | 54 km | 81 km | 108km | ||
| PM270-1-1-2(−6) | 1.76 | 2.06 | 2.43 | 2.76 | 3.21 | 1.36 | |
| PM270s-1-1-2(−6) | 1.75 | 2.02 | 2.39 | 2.77 | 3.25 | 1.36 (2.25) | |
| RM270-1-1-2(−6) | 1.47 | 1.68 | 1.93 | 2.24 | 2.62 | - | |
| PM270-2-1-2(−6) | 5i | 3.91 | 4.50 | 5.65 | 6.87 | 8.36 | 2.98 |
| PM270s-2-1-2(−6) | 5j | 3.91 | 4.47 | 5.56 | 6.87 | 8.46 | 2.98 (3.87) |
| RM270-2-1-2(−6) | 3.37 | 3.86 | 4.35 | 4.90 | 5.69 | - | |
- a The simulations are titled according the choice of the various parameters. These parameters are separated by hyphens. In order these parameters are (1) the simulations types PM and RM and the radius of the planetesimals subsequent to complete sintering. In the case of the partial melt scenario the simulations with the slower melt percolation velocity of 0.3 m year−1 are designated with a suffix ‘s’ with the radius of the planetesimal. The remaining simulations were run with the melt percolation velocity of 300 m year−1. (2) The onset time, tOnset(Ma, Million year), to initiate the accretion of a planetesimal from the time of the formation of the CAIs with the canonical value of (26Al/27Al)initial = 5 × 10−5. (3) The accretion duration, tDuration(Ma), of the planetesimal. (4) An initial value 5 × 10−7 [represented as 5(−7)], 1 × 10−6 [represented as 1(−6)] and 2 × 10−6 [represented as 2(−6)] for the 60Fe/56Fe ratio at the time of formation of the CAIs with the canonical value of (26Al/27Al)initial.
- b Time (Ma) taken for the Fe-FeS core to grow to a specific size. Five different arbitrary choices of the core sizes have been considered. With respect to the final radius of the completely sintered planetesimal, these core sizes are expressed in percentage. In order these are 0+% (initiation of core formation), 10%, 20%, 30% and 40%. All time spans are measured with respect to the formation of the CAIs with the canonical value of (26Al/27Al)initial = 5 × 10−5.
- c Time obtained for the initiation of the extrusion of the basaltic melt for PM simulations with respect to the formation of the CAIs with the canonical value of (26Al/27Al)initial. The values in the bracket represent the time for the extrusion of the basaltic melt to the surface of the planetesimals.
[30] The Table 2 presents the dependence of the initiation of the basaltic melt extrusion for different onset times of the accretion of the planetesimals on the initial bulk composition of a 100 km sized planetesimal. The basaltic melt can move to the planetesimal surface over the timescale with uncertainty comparable to the mean-life of 26Al [Taylor et al., 1993]. The Table 3 presents the dependence of the onset time of the initiation of core-mantle differentiation on the initial bulk composition of a 100 km sized planetesimal. We have considered both scenarios in order to study the effect of bulk initial compositions on the differentiation of the planetesimals. The results for the planetesimals with four distinct chemical compositions L, LL, CI, and CV and a sets of two binary combinations of compositions, one with (0.7LL + 0.3CI) and the other with (0.7L + 0.3CV) [Jurewicz et al., 1991, 1993, 1995; Taylor et al., 1993; Righter and Drake, 1997] have been compared with the H chondritic composition.
| TOnset (Ma) | H | L | LL | CI | CV | 0.7LL + 0.3CI | 0.7L + 0.3CV |
|---|---|---|---|---|---|---|---|
| 0 | 0.15 | 0.16 | 0.16 | 0.21 | 0.12 | 0.18 | 0.14 |
| 0.5 | 0.74 | 0.75 | 0.76 | 0.84 | 0.69 | 0.78 | 0.73 |
| 1.0 | 1.38 | 1.41 | 1.42 | 1.55 | 1.30 | 1.45 | 1.36 |
| 1.5 | 2.09 | 2.17 | 2.20 | 2.44 | 1.97 | 2.26 | 2.10 |
| 2.0 | 2.99 | 3.19 | 3.26 | 3.82 | 2.81 | 3.39 | 3.04 |
| 2.5 | 4.35 | 5.02 | 5.43 | none | 4.0 | 5.99 | 4.59 |
- a An accretion duration (tDuration) of 1 Ma and (60Fe/56Fe)initial = 2 × 10−6 were adopted.
| TOnset (Ma) | H | L | LL | CI | CV | 0.7LL + 0.3CI | 0.7L + 0.3CV | |||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| PM | RM | PM | RM | PM | RM | PM | RM | PM | RM | PM | RM | PM | RM | |
| 0 | 0.39 | 0.19 | 0.45 | 0.20 | 0.47 | 0.20 | 0.57 | 0.27 | 0.36 | 0.14 | 0.49 | 0.22 | 0.41 | 0.18 |
| 0.5 | 1.05 | 0.80 | 1.13 | 0.82 | 1.16 | 0.82 | 1.32 | 0.93 | 0.99 | 0.73 | 1.20 | 0.85 | 1.08 | 0.78 |
| 1.0 | 1.80 | 1.48 | 1.93 | 1.52 | 1.98 | 1.53 | 2.25 | 1.73 | 1.70 | 1.36 | 2.05 | 1.58 | 1.85 | 1.46 |
| 1.5 | 2.72 | 2.30 | 2.96 | 2.40 | 3.06 | 2.43 | 3.63 | 2.83 | 2.56 | 2.10 | 3.19 | 2.52 | 2.81 | 2.28 |
| 2.0 | 4.05 | 3.43 | 4.67 | 3.74 | 4.97 | 3.86 | 8.40 | 5.32 | 3.73 | 3.07 | 5.42 | 4.13 | 4.29 | 3.46 |
| 2.5 | 7.44 | 5.84 | none | none | none | none | none | none | 6.09 | 4.80 | none | none | none | 6.74 |
- a An accretion duration (tDuration) of 1 Ma and (60Fe/56Fe)initial = 2 × 10−6 were adopted.
4. Discussion
[31] Numerical simulations for the planetary differentiation of Vesta and other differentiated asteroids have been performed. With the present work, along with the earlier study [Sahijpal et al., 2007], a systematic theoretical framework has been developed to critically examine the various proposed planetary differentiation scenarios that differ from each other in terms of the relative timings between the core-mantle and mantle-crust differentiation besides the geochemical differences associated with various differentiation scenarios [Taylor et al., 1993; Righter and Drake, 1997].
[32] The dependence of the thermal evolution of the planetesimals on the various simulation parameters, specifically on the time of the initiation of the planetesimal accretion (tOnset) subsequent to the formation of the CAIs, the accretion duration (tDuration), the size of the planetesimals and the choice of (60Fe/56Fe)initial has been extensively studied in the earlier studies [Merk et al., 2002; Hevey and Sanders, 2006; Sahijpal et al., 2007]. The initiation and the growth of the Fe-FeS core, and the initiation of the melting of the bulk silicate critically depend upon these parameters, specifically on the onset time of the accretion of the planetesimals in the time interval ∼2–3 Ma. We avoid detailed discussion on these aspects. The Figures 5 and 6 along with the Tables 1–3 rather provide insight into these aspects for the presently adopted differentiation scenarios of the planetesimals with different compositions. In the present work, we focus on the outcome of the two presently simulated differentiation scenarios regarding the origin of the main group of basaltic achondrites along with the possible implications on the evolution of differentiated planetesimals. We have parametrically incorporated the influence of thermal convection in both the models. Due to thermal convection the heat transfers efficiently from the molten interiors to the overlying partially melted and un-melted exteriors of the planetesimal. This leads to a rapid heating and melting in the exterior regions (Figures 5 and 6). In the following text, we discuss the two differentiation scenarios and finally make assessments regarding the likely hood of these scenarios to be the cause of planetary differentiation of Vesta and other differentiated asteroids.
4.1. Partial Melt Origin of the Basalts: PM Model
[33] In the case of the PM simulations, the mantle-crust differentiation occurs prior to the core-mantle differentiation at 20% partial melting of the bulk silicate. Depending upon the onset time of the planetesimal accretion and the bulk composition of the planetesimal, the onset of the basaltic melt extrusion can occur any time from ∼0.15–6 Ma after the condensation of CAIs (Tables 1 and 2). Due to the uncertainty in the basaltic melt percolation velocity [e.g., Taylor et al., 1993] an additional time of up to a million years or more could be required for the crustal formation in case we consider the basaltic melt percolation velocity ≤0.3 m year−1 instead of the velocity of 300 m year−1 considered in majority of the simulations (Table 1). In general, the deduced onset time of the basaltic crustal formation is consistent with the recent compilation of the published ages of the basaltic meteorites within 4–10 Ma [e.g., Sanders and Scott, 2007, Figure 1; Nyquist et al., 2003, 2009; Kleine et al., 2009; Wadhwa et al. 2009], except for the scarcity of specimen within the initial couple of million years. The early dearth of the basaltic specimen could be due to the prolonged accretion of the planetesimals that could have obscured the early initial basaltic events by producing a chondritic shell around the differentiated body. It should be mentioned that in the present study we have performed simulations with planetesimal accretion duration up to one million years. A prolonged accretion of planetesimals [Weidenschilling, 2000] could have delayed the thermal processing. Further, as discussed earlier in case we consider the uncertainty in the basaltic melt percolation velocity the deduced temporal range of crustal formation would almost overlap with the published ages of the basaltic meteorites. However, it should be emphasized that the chronological records of the iron meteorites and achondrites seem to indicate that the core-mantle differentiation in general occurred prior to mantle-crust differentiation. This seems to contradict the partial melt scenario unless we speculate that the iron meteorite parent bodies formed and differentiated earlier compared to the prolonged accretion and differentiation of the distinct parent bodies of the achondrites. Nonetheless, the geochemical evidence of siderophile elements in eucrites [Righter and Drake, 1997] is difficult to reconcile with the partial melt origin of the basalts and hence remains one of the major dilemma for the scenario to explain differentiation of eucrite parent bodies.
[34] The extrusion of the basaltic melt to the surface in the case of PM simulations alters the subsequent thermal evolution of the planetesimal. The extrusion of the basaltic melt along with 26Al to the surface results in the retardation of the initiation of the formation and the growth of the core. This fact is clearly demonstrated in the case of smaller bodies, e.g., the 20 km sized planetesimals. In the case of simulations PM20-2-0.001–2(−6) (Figure 5a) and PM20-2-0.1–2(−6) (Table 1), the temperature at the center does not reach high enough even to initiate the core-mantle differentiation. The models with an identical set of parameters but with no 26Al rich basaltic melt extrusion would produce core-mantle differentiation. However, the influence of the basaltic melt percolation velocity on the rate of core-mantle differentiation is found to be insignificant for the melt percolation velocities attempted in the present work. As an example the simulations, PM100-2-1-2(−6) and PM100s-2-1-2(−6), with basaltic melt percolation velocities of 300 and 0.3 m year−1, respectively, infer almost identical rate of core-mantle differentiation.
[35] There are certain issues regarding sintering and differentiation of the planetesimals that need to be addressed in context to the simulation results. It should be mentioned that the onset time and the duration of the planetesimal accretion significantly influence the sintering process. As an example, in the case of the simulations PM20-2-0.001–2(−6) and PM50-2-0.1–2(−6) (Figures 5a and 5b, respectively), due to the rapid accretion, the process of sintering initiate after the accretion of the entire planetesimals. However, in the case of the simulations PM50-1-1-2(−6), PM100-1-1-2(−6) (Figures 5c and 5d), due to the prolonged accretion of the planetesimals over 1 Ma, the sintering in the deep interiors commence prior to the entire accretion. A sintering front propagates from the interior to the surface of the planetesimal. It should be noted that we have not presently worked-out the scenarios with the accretion timescales greater than 1 Ma as these scenarios could result in a complex situation where the planetary differentiation will commence in the interiors while the sintering and accretion in proceeding in the outer regions of the planetesimals.
[36] During the core formation, the radionuclide 60Fe moves inside the core. The further increase in the core temperature is partially due to its decay. The increase in the temperature of the core would primarily depend upon the initial 60Fe/56Fe. Due to the uncertainty in the initial 60Fe/56Fe, we have considered the (0.5–2) × 10−6 at the time of the CAIs formation. Apart from the 60Fe heat source, the heat is also transferred within the core due to the assumed convection in the molten iron core from the overlying mantle. We assumed the absence of any thermal gradient at the core-mantle boundary. This is reflected in almost uniform temperature isochrones in the simulations after the complete melting of the planetesimals. The assumption is supported by the fact that as long as the core and the mantle are both convective, the development of a strong thermal gradient at the core-mantle boundary will bring in thermal equilibrium among the two zones. It should be mentioned that it is difficult to handle the physical problems associated with discontinuities. Since the core-mantle discontinuity is essentially a step function in composition, the presence of any thermal gradient across the discontinuity will produce large thermal flow according to Fick's law.
[37] In order to understand the influence of the initial bulk composition on the differentiation of the planetesimals, we have run PM and RM simulations with seven different initial compositions, with the H chondritic composition as a benchmark (Tables 2 and 3). The differences in the bulk compositions significantly control the thermal evolution of the planetesimals. We have observed that the extrusion of the basaltic melt and the initiation of the segregation of the core occur fastest in the numerical simulations with CV chondrite composition, and slowest in the numerical simulations with CI chondrite composition (Tables 2 and 3). Further, the size of the core depends upon the mass abundance of the Fe that segregates to form the core.
4.2. Basalts From the Residual Melt in the Cooling Magma Ocean: RM Model
[38] The detailed studies of the siderophile element partitioning indicate that the basaltic achondrites were probably derived from the residual melts of a cooling magma ocean [Righter and Drake, 1997]. One of the main objectives of the present study was to understand the origin of the basaltic achondrites within the temporal constraints imposed by the observed ages of the basaltic meteorites. In the case of RM simulations, even in the presence of convection in the magma ocean, parametrically modeled by a three orders of magnitude high thermal diffusivity, we could not observe rapid and global cooling of the magma oceans in the planetesimals within 6–7 Ma for the planetesimals of sizes ≥50 km (Figure 6). The situation becomes even worse in the case of a 270 km sized asteroid, e.g., the Vesta that retains its internal heat for a long time due to its large size (Figure 6f). Hence, we do not anticipate major crystallization in the magma ocean of a 270 km sized planetesimal within the initial ∼6–7 Ma that would produce residual melts for the basalts on account of equilibrium crystallization as proposed by Drake [2001]. This is in contradiction with the published ages of the achondrites that are considered to have formed at least within the initial ∼7 Ma [Nyquist et al., 2003, 2009; Wadhwa et al. 2009; Sanders and Scott, 2007].
[39] The reason for the absence of rapid and global cooling of magma ocean is partially due to the presence of an un-melted chondritic ‘insulating’ crust of thickness few couple of kilometers that survives on the surface of the planetesimal, and partially due to the ongoing radioactive heating [Elkins-Tanton et al., 2008]. On the basis of a realistic treatment of convection, Elkins-Tanton et al. [2008] predicted the solidification of Vesta to occur on a timescale of the order of 100 years in the absence of outer crust and radiogenic heating. However, a comprehensive thermal model involving the cooling of magma ocean has not been developed so far.
[40] As mentioned earlier the chronological records of iron meteorites and achondrites broadly seem to indicate that the core-mantle differentiation occurred prior to mantle-crust differentiation. The scenario involving the residual melt origin of the basalt seems to have a more viable chronological explanation for various differentiated meteorites compared to partial melt scenario. The residual melt scenario is also supported by the geochemical evidence of siderophile elements partitioning in eucrites. The failure to rapidly cool the magma ocean seems to be the major short-comings of the residual melt scenario. In order to make the residual origin of the basaltic melt scenario viable by rapidly and globally cooling the magma ocean [Righter and Drake, 1997; Drake, 2001] it would be essential to frequently remove the insulating crust. The insulating crust is essentially the chondritic crust that survived the widespread melting of the planetesimal, and partially the cooling front of the magma ocean. As mentioned earlier, the variation in the viscosity by more than an order of magnitude can reduce the extent of convection in the outer regions of the magma ocean where we encounter the maximum thermal gradient [Marsh, 1989; Hort et al., 1999]. This would result in the formation of a rigid crust, mush and suspension that in turn would provide an insulating blanket to the cooling magma ocean. This thermal insulation along with the insulation provided by chondritic crust layer would eventually control the global cooling and the crystallization in the magma ocean. Rapid cooling could be achieved only if this insulating layer is efficiently removed. This could be achieved by continuous bombardment of the cooling crust by small planetesimals that would break the crust and maintain a convective magma ocean up to the surface. In the early bombardment era of the solar system it would be possible to speculate continuous bombardment of small planetesimals on Vesta and other planetesimals undergoing differentiation [Weidenschilling, 2000]. The cooling of the magma ocean will in turn produce fresh quenched crust that has to be again removed by planetesimal bombardment for efficient cooling. The wreckage crust would sink in the magma ocean and will be further replaced by the fresh quenched crust [Taylor et al., 1993]. These complex processes have to be tuned in a manner to produce basaltic achondrites within the initial 10 Ma of the solar system.
[41] Further, it should be mentioned that we could not observe rapid and global cooling of the magma ocean as anticipated by Drake [2001]. However, we cannot rule out the possibility of localized cooling of the magma ocean just beneath the outer insulating crust to produce localized regions where crystallization can produce basaltic melt residues. This scenario is much complex than the original proposed model of equilibrium crystallization. Hence, this scenario would require further scrutiny in terms of the geochemical evidence of the achondrites, and the treatment of convection in the magma ocean where the crystallization commence in the upper cooling regions.
5. Summary
[42] Numerical simulations have been performed for the differentiation of Vesta and other differentiated asteroids following two distinct planetary differentiation scenarios. These scenarios deal with the origin of the basaltic achondrites either by the partial silicate melting or from the residual melt left subsequent to the crystallization in a cooling magma ocean. The accretion of the asteroids within the initial ∼2 Ma resulted in the production of magma oceans. Convection was parametrically incorporated in the molten core and magma ocean. The partial melt origin model is consistent with the chronological records of the basalts except for the temporal sequence in the core-mantle and mantle-crust differentiation that seems to contradict the observed chronological records. Further, it is difficult to reconcile the partial melt scenario of basalt origin with the geochemical records of siderophile elements in eucrites. In the case of residual melt model that is compatible with the geochemical records of siderophile elements in eucrites, we could not observe rapid and global cooling of the planetesimals of sizes ≥50 Km within the initial 7 Ma. In order to make this model viable it would be essential to frequently and efficiently remove the insulating crust by continuous bombardment of small planetesimals. Alternatively, there is a possibility that the localized cooling of the magma ocean just beneath the insulating crust of the planetesimal produce regions where early crystallization produce residual basaltic melt plumes.
Acknowledgments
[43] We are grateful for the numerous comments made by two anonymous reviewers that led to significant improvement of the manuscript. We gratefully acknowledge the PLANEX (ISRO) research support and the associateship program of IUCAA, Pune.
