Impact-induced melting during accretion of the Earth

Because of the high energies involved, giant impacts that occur during planetary accretion cause large degrees of melting. The depth of melting in the target body after each collision determines the pressure and temperature conditions of metal-silicate equilibration and thus geochemical fractionation that results from core-mantle differentiation. The accretional collisions involved in forming the terrestrial planets of the inner Solar System have been calculated by previous studies using N-body accretion simulations. Here we use the output from such simulations to determine the volumes of melt produced and thus the pressure and temperature conditions of metal-silicate equilibration, after each impact, as Earth-like planets accrete. For these calculations a parametrised melting model is used that takes impact velocity, impact angle and the respective masses of the impacting bodies into account. The evolution of metal-silicate equilibration pressures (as defined by evolving magma ocean depths) during Earth's accretion depends strongly on the lifetime of impact-generated magma oceans compared to the time interval between large impacts. In addition, such results depend on starting parameters in the N-body simulations, such as the number and initial mass of embryos. Thus, there is the potential for combining the results, such as those presented here, with multistage core formation models to better constrain the accretional history of the Earth.

Geochemical models of core formation can be used to constrain the early melting history of the Earth.
The segregation of liquid metal from silicate during formation of the Earth's core depleted siderophile (metal-loving) elements in the Earth's mantle by transporting them into the core. Thus elements such as Fe, Ni, Co, Mo, W, V and Cr are depleted in the mantle by up to 2 orders of magnitude relative to their abundances in CI chondrites (e.g. Rubie et al., 2015a). The degree of the depletion of moderately siderophile elements is variable and depends on the metal-silicate partition coefficient which, for element M, is defined as the ratio of the molar concentration of M in the metal to its concentration in the silicate: The partition coefficient M met sil D − depends on pressure (P), temperature (T), oxygen fugacity (f O2 ) and, for some elements, the compositions of the metal and silicate phases. Because M met sil D − depends on P and T, the conditions of core formation have been estimated by matching experimentally-determined metal-silicate partition coefficients with core-mantle partition coefficients determined from McDonough (2003). For example, Li and Agee (1996) showed that to match M met sil D − for Ni and Co with core-mantle partition coefficients requires an equilibration pressure of ~28 GPa. Subsequent studies have derived conditions of core formation that range from 27 to 60 GPa and 2200 to 4200 K (see Table   3 in Rubie et al., 2015a). However, this approach is based on the concept of "single-stage core formation" whereby all the metal of the core is considered to equilibrate with all silicate of the mantle at some mid-mantle depth (e.g. Corgne et al., 2009;Righter 2011;Walter and Cottrell, 2013). In reality, the Earth accreted through a series of high-energy impacts with smaller bodies consisting of multi-km-size planetesimals (possibly hundreds of km in diameter) and Moon-to Mars-size embryos (e.g. Chambers and Wetherill, 1998). Most impacts, while delivering the energy that caused melting, added Fe-rich metal that segregated to the Earth's proto-core. Thus core formation was multistage and occurred as an integral part of the accretion process (Wade and Wood, 2005;Rubie et al., 2011).
The preliminary multistage core formation model of Rubie et al. (2011) has been integrated recently with N-body accretion simulations (Rubie et al., 2015b). Hundreds of impacts and associated coreformation events were simulated for all embryos as well as for the final terrestrial planets. A mass balance/element partitioning approach was used, which required the bulk compositions of all accreting bodies to be defined. Elements were assumed to be present mostly in Solar System (CI) relative abundances and the oxygen content was the main compositional variable. Metal-silicate equilibration pressures and bulk compositional variables were refined by least squares minimization in order to fit the calculated composition of Earth's mantle to that published for the primitive mantle (Palme and O'Neill, 2013); equilibration temperature was set midway between the peridotite solidus and liquidus at the equilibration pressure. A simplifying assumption is that metal-silicate equilibration pressures were a constant fraction of the core-mantle boundary pressure at the time of each impact, irrespective of the mass of the impactors. As constrained by the least squares fits to Earth's mantle composition, this fraction is 0.6-0.7. Thus, pressures increased as the Earth's mass increased during accretion and reached a maximum of ~80 GPa.
It has been assumed, in almost all studies of siderophile element partitioning, that the pressure of metalsilicate equilibration during core formation corresponds to the pressure at the base of the magma ocean, with temperature being defined by the equivalent peridotite liquidus or solidus. In the case of dispersed metal droplets that continuously re-equilibrated as they sink in a magma ocean, more realistic fractionation models have been proposed according to which the effective equilibration pressure differs significantly from the pressure at the base of a magma ocean (Rubie et al., 2003). However, accretion of Earth mostly involved collisions with bodies that previously had undergone core-mantle differentiation (e.g. Kleine et al., 2009). In this case, it is likely that the core of the impactor would sink through the magma ocean entraining silicate liquid to form a descending, expanding high-density plume of mixed metal and silicate (Deguen et al., 2011;Rubie et al., 2015b). The hydrodynamic model formulated by Deguen et al. (2011) to describe this plume has two important consequences for modelling core formation: (1) Only a small fraction of the silicate mantle of the target body equilibrates with the accreted metal during each impact (Rubie et al., 2015b). (2) The pressure of final metalsilicate equilibration corresponds to the pressure near the base of the magma ocean where the descending metal-silicate plume comes to rest. Based on this latter point, it is clear that understanding the evolution of magma ocean depth during accretion is of fundamental importance for modelling coremantle differentiation.
Complete chemical equilibrium between metal and silicate during core formation requires that the iron emulsifies into small (cm size) droplets (Rubie et al., 2003). In the case of small impacting bodies (planetesimals), emulsification is likely to be complete but for embryos the extent of emulsification is much less certain (Dahl and Stevenson, 2010;Samuel, 2012;Deguen et al., 2014;Kendall and Melosh, 2016). The extent of emulsification and the fraction of metal that equilibrates is of major importance for understanding the origin of planetary tungsten isotopic anomalies (Nimmo et al., 2010;Rudge et al., 2010) and also affects the final distribution of siderophile elements (Rubie et al., 2015b).

Aims of the study
The aim of this contribution is to present the results of preliminary calculations of melting depths caused by each individual impact in N-body accretion simulations and to develop a complete history of melting and evolving magma ocean depth throughout Earth's accretion history. Because magma ocean depths are related to the conditions of metal-silicate equilibration, this approach could eventually be used to input equilibration pressures and temperatures into the multistage accretion/core formation model of Rubie et al. (2015). This would enable the assumption that equilibration pressures are a constant fraction of core-mantle boundary pressures to be superseded. The approach may also enable the most realistic accretion models to be identified.

Model description
The volume of melt produced during an accretional collision depends on the kinetic energy of the impactor, as determined by its mass and velocity, and the impact angle which determines how much of the kinetic energy is transferred into the target body. These three parameters (impactor mass, impact velocity and impact angle) and the mass of the target body are determined for each impact in the Nbody accretion simulations. However, since the simulations involve 1000 to 2000 impacts, it would be much too time consuming to perform full three-dimensional hydrocode calculations of the melt volume produced by each impact (Marinova et al. 2011). Two-dimensional models also cannot be used for nonvertical impacts due to their assumed symmetry in the third dimension. We have therefore based our calculations on the analytical model of Bjorkman and Holsapple (1987), who constructed a simple relationship between melt volume, impactor mass and impact velocity for vertical impacts using dimensional analysis. Based on scaling theory and hydrocode simulations (Pierazzo and Melosh, 2000), a dependence on the impact angle has been included to cover the whole range of impact parameters that occur during planetary accretion (Abramov et al., 2012). Although our approach is undoubtedly simplified, it reveals general trends (such as the importance of magma ocean freezing time) that are unlikely to change significantly if more sophisticated approaches were to be implemented. The simple approach adopted here also makes the dependence of the results on input parameters -which are often poorly-known -more transparent.
We have used this model to calculate the melt volume produced by 1000-2000 impacts that occur during each of three "Grand Tack" N-body accretion simulations (Walsh et al., 2011; and have determined the corresponding metal-silicate equilibration pressures at the base of the resulting melt pools and global magma oceans.  (2014) and Rubie et al. (2015) and are chosen because they reproduce well the mass and orbital characteristics of Earth (see also below).

Melt volume
According to the analytical model of Bjorkman and Holsapple (1987), the mass of melt M melt produced by a vertical impact is given by: (1) Here m p is the mass of the projectile, v is the impact velocity, k and µ are scaling parameters and E m is the energy required to melt the target upon decompression (see also Table 1). The values of constants in this relationship have been constrained through laboratory experiments (e.g. Schmidt and Housen, 1987) and through numerical simulations (e.g. Elbeshausen et al., 2009;Barr and Citron, 2011).
By including a dependence on impact angle (Pierazzo and Melosh, 2000), the analytical model has been modified by Abramov et al. (2012) to give: Here V melt is the melt volume, ρ p and ρ t are the densities of the projectile and target respectively, D p is the projectile diameter, γ is a scaling parameter and θ is the impact angle (which for a vertical impact is defined as 90°). For simplicity, we assume that the mean densities of the bodies scale linearly with their mass, which is consistent with an increased self-compression for a larger body (Rubie et al., 2011).
If the target is already at elevated temperature (e.g. due to impacts), the energy required for melting is reduced. To take this into account, E m is calculated when the target is already at high temperature from: (adapted from Abramov et al., 2012), where T s is the surface temperature prior to an impact and T l0 is the liquidus temperature at the surface (1950 K); the temperature, T, and liquidus temperature, T l , are calculated at the maximum depth of melting, d m . For the liquidus temperature this is done through a temperature increase as a function of pressure, P. C p denotes specific heat and L m is the latent heat of melting. In the calculations described below, the surface temperature is set to 1750 K (Table 1) for a solid body before impact and the temperature increase with depth (0.1 K/km, based on a lower-mantle temperature gradient - Monnereau and Yuen, 2002), is assumed to be linear. Values of the material properties used here are listed in Table 2.
The depth of an impact-induced melt pool is calculated by solving Eqs. 2 and 3 numerically based on assuming that the melt pool has a specific geometry. Here we assume a spherical geometry (see Fig.   1a) which, in contrast to a hemispherical geometry, takes account of the free surface and closely approximates actual melt pool geometries . The pressure and temperature at the base of the melt pool can then potentially be used as the conditions of metal-silicate equilibration during an episode of core formation.
Each body in the N-body simulations is assumed to be differentiated into an iron core and a dunite mantle, where the core size is calculated from the mean planetary density, assuming a metal mass fraction of 0.34 and a core density that is a factor 2.5 larger than the mantle density. The core temperature of the target body depends on the deep mantle temperature with a temperature jump ΔT at the core-mantle boundary (CMB) of 100 K. Equations are first solved numerically using mantle parameters. If the depth of melting extends into the core, the equations are solved again in the core domain, with core values for all parameters, to determine the extent of core melting. Due to the spherical geometry of the melt pool, the amount of silicate mantle melting increases when the depth of melting extends deeper into the core.

Magma ocean evolution
If sufficient melt is produced as a spherical melt pool by an accretional impact, isostatic readjustment can cause the melt to spread over the surface of the target body to form a global magma ocean as depicted in Fig. 1 (Tonks and Melosh, 1992). This requires extensive solid-state deformation of the crystalline mantle material that is located below and adjacent to the melt pool (Fig. 1a). Reese and Solomatov (2006) have shown that the formation timescale of a global magma ocean is controlled by radial relaxation of solid mantle, because this timescale is much longer than the timescale of lateral melt flow. The radial relaxation timescale for a Mars-size body has been estimated by Reese and Solomatov (2006) and varies by 5-6 orders of magnitude depending on the volume of the impactinduced melt pool and the rheology of the underlying crystalline mantle. For small amounts of melting, the timescales are on the order of 10 4 -10 5 years (Reese and Solomatov, 2006). It is therefore unlikely that planetesimal impacts result in a global magma ocean because the melt would crystallize before it could spread over the surface, as discussed below. The timescale of radial relaxation decreases significantly as the volume of the melt pool become large (as in the case of a giant impact), especially in the case of a non-Newtonian rheology for the solid part of the proto-planet, which makes the rapid formation of a global magma ocean by isostatic re-adjustment following a giant impact highly likely.

Cooling times of magma oceans
Some of the impacts as calculated from the N-body accretion models may occur on a molten surface if a global magma ocean persists for a longer period of time. To determine whether this is the case, a simple model is used to estimate the cooling time of a global magma ocean. This model is described by the following equation: which balances the heat flux, F, through the planet's surface on the left hand side, with the latent heat, L m , that is released during crystallisation (first term on the right) and the heat released through secular cooling (second term on the right) (Elkins-Tanton, 2008). For symbol explanation, see Table 1. Here the heat flux is approximated by a constant value, which is certainly a simplification. However, what really matters is whether the solidification timescale is large or small compared to the time interval between giant impacts. The latter parameter is provided by the accretion simulations.
If cooling occurs via radiation from a deep magma ocean, without an insulating atmosphere or protocrust, the Rayleigh number is extremely high (10 27 -10 32 ), partly because of the extremely low viscosity of peridotite melt (Liebske et al., 2005). Thus, convection is in the hard turbulent regime with convective velocities of several meters per second. The heat flux is extremely high and the solidification timescale, at least for the deeper (lower-mantle) part of the magma ocean, is on the order of 10 3 years (Solomatov, 2000(Solomatov, , 2015Rubie et al., 2003) -much shorter than the time intervals between impacts. On the other hand, if a solid crust forms, the cooling rate is limited by conduction through this lid and mantle solidification could then take on the order of 10 8 years -much longer than the intervals between impacts. An intermediate case arises when a thick steam atmosphere is present, in which case the heat flux is limited to 300-500 W/m 2 and is set by the thermodynamic properties of the atmosphere (Sleep et al. 2014). In this case the solidification timescale can be comparable to the impact interval and here the assumption of a constant heat flux is quite good, because it depends on atmospheric and not surface or interior properties. Below we present the results of calculations assuming either a thick steam atmosphere (F=475 W/m 2 ) or a radiating magma ocean (F=2×10 5 W/m 2 ).
The former value was determined by fitting the cooling times of Lebrun et al. (2013) with our cooling model and is of the same order of magnitude as values determined for the maximum heat flux through a dense atmosphere (Sleep et al., 2014;Hamano et al., 2013). The latter value was determined using the same fitting procedure for a body without an atmosphere; the result is similar to the radiation of a black body with a surface temperature of ~1000 °C. Note that the constant heat fluxes used in our models represent time-averaged heat flows. In reality, heat fluxes will initially be relatively high and will decrease with time.

Melt production
For each impact in an N-body accretion simulation, the impact angle, impact velocity and the masses of the two bodies involved are used to calculate melt production using the equations listed above. The impact histories from three different N-body simulations are used here. These are 4:1-0.25-7, 4:1-0.5-8, and i4:1-0.8-4 Rubie et al., 2015b) and were chosen because they result in realistic Earth-mass planets at ∼1 AU (see Fig. 1 in Rubie et al., 2015). Here "4:1" designates the ratio of the total mass of all the embryos to the total mass of all the planetesimals at the beginning of the simulation. The middle term (e.g. "0.25") is the initial mass of embryos as a fraction of one Mars mass. For each set of starting parameters, as defined by these two terms, 10 simulations were run with very slight variations in the initial orbital characteristics of the starting bodies. The last term is the run number within the set of 10 simulations. The "i" in the third simulation number indicates that the starting mass of embryos increases with increasing heliocentric distant (from 0.7 to 3.0 AU) -see Table 3. The main difference between the three simulations used here is the initial mass of the embryos which ranges from 0.2 to 0.8 times the mass of Mars. Table 3 provides an overview of the initial parameters of the N-body models. Each simulation produces several final terrestrial planets and we define "Earth-like" planets as those that form close to 1 AU and have a mass close to one Earth mass (Rubie et al., 2015b).
Once the depth of melting for each impact has been calculated, the pressure and temperature at that depth is determined. Our ultimate aim is to use such conditions, and those at the base of subsequent global magma oceans, to model each episode of core formation and the resulting siderophile element partitioning throughout Earth's accretion. In the future, such models can potentially distinguish which N-body models most closely resemble Earth's accretion history if the calculated melt pressures for different N-body models are significantly different. Planetesimal impacts nominally create shallower melt pools than embryo impacts with maximum pressures that range from 0.4 to 1.2 times the core-mantle boundary pressure. A few high-velocity impacts create even deeper melt pools. However, the results shown in Fig. 2 for planetesimals are not realistic. This is because each "planetesimal" in the N-body simulations is actually a tracer that is used to represent a swarm of much smaller bodies in order to reduce calculation times (O'Brien et al., 2006).
Actual planetesimals likely have diameters of the order of 100 km. Impacts of bodies of this size will create only small amounts of melting and this is taken into account in the following discussion and the results that are presented below.
There are two end-member scenarios for planetesimal impacts. (1) They impact a solid surface due to the rapid crystallisation of a preceding magma ocean. In this case there will be little melting and the cores of such bodies will only segregate to the Earth's core when the next giant impact occurs and produces extensive melting.
(2) Small bodies impact a pre-existing global magma ocean, created by an earlier giant impact, and metal-silicate equilibration takes place near the base of this magma ocean. In this case, the unrealistically large size of planetesimals will not significantly change the results, because the effect of a large number of small bodies impacting in a magma ocean around the same time will be similar to the effect of one large body with the same total mass impacting the magma ocean.

Magma ocean crystallisation
To determine if some planetesimal impacts occur onto a pre-existing global magma ocean, we use the simple magma ocean crystallisation model described above (Eq. 4). If the time until the next impact on the proto-Earth is shorter than the crystallisation time of the magma ocean, the next impact will occur on a molten surface. In general, at the start of accretion the time intervals between impacts are relatively short. Most of these impacts are small and will cause only small amounts of melting.
However, some of these collisions are between two similar-sized bodies resulting in a deep magma ocean that may not crystallise before the next impact.
The cooling model for magma oceans described above (Eq. 4), based on a constant surface heat flux, is used to estimate the crystallization time of a global magma ocean. To calculate the magma ocean depth from the volume of the spherical impact-induced melt pool, the silicate melt volume is spread out geometrically over the entire surface of the new body after collision and accretion. Since the cooling of a magma ocean depends strongly on the presence or absence of an insulating atmosphere, the results of the two end-member cooling models described above are presented here.
Magma ocean crystallization times are plotted as a function of the time until the next impact in Figure  3. Only embryo-embryo collisions that form the final Earth-like planet are shown. In the absence of an insulating atmosphere, almost all magma oceans crystallise before the next impact (almost all points plot above the dotted 1:1 line). When a dense insulating atmosphere is present, the majority of magma oceans would still be present when the next impact occurs.

Metal-silicate equilibration depths
The magma ocean crystallisation model described above is used to determine the depth of the crystallizing magma ocean at the time of the next planetesimal impact. Based on the hydrodynamic model of Deguen et al. (2011), the metal from the impacting planetesimal can then be assumed to equilibrate with silicate liquid at this depth (Rubie et al. 2015b). In the case of an embryo impact onto a pre-existing magma ocean, the volume of the calculated melt pool that extends below the base of the magma ocean into previously solid mantle is added to the volume of the magma ocean to calculate a new magma ocean depth. Metal-silicate equilibration of the embryo material will likely occur at the bottom of the melt pool prior to isostatic readjustment.
Assuming that metal-silicate equilibration and segregation of metal to the proto-core occur only when a large volume of melt is created by an embryo-embryo collision (giant impact), we have the following possibilities for equilibration: • Embryo-embryo collisions: Equilibration takes place at the bottom of the melt pool (or at the coremantle boundary if the melt pool extends into the core) before a global magma ocean is created by isostatic readjustment and lateral spreading. In this case the equilibration pressure will generally be close to the core-mantle boundary pressure (Figure 4, large symbols).
• Planetesimal impact in a magma ocean that survives from a previous embryo-embryo collision: Equilibration will take place at the base of the magma ocean. The depth of the magma ocean is dependent on the time between the embryo-embryo collision and the planetesimal impact as well as the cooling rate of the magma ocean (Figure 4, points contained in the black ellipse).
• Planetesimal impact on a solid surface before the last giant impact: In this case, there is insufficient melt at the time of impact for significant equilibration and metal-silicate equilibration.
However, the next giant impact will mix this material with the target material. Therefore, equilibration takes place at the bottom of the deep melt pool created by the next giant impact ( Figure 4, points contained in the magenta ellipse).
• Planetesimal impact on a solid surface after the last giant impact (or in the case of no giant impacts): in this case there is no subsequent significant melting event available to equilibrate the material and little mixing of the impactor and target material will occur, the impactor material forms a late veneer (Figure 4, points contained in the red ellipse).
The evolution of metal-silicate equilibration pressures for the Earth-like planets in the three N-body accretion simulations and, in each case, for the two end-member cooling/atmosphere models are shown in Figure 5. Embryo-embryo collisions (large symbols) are visible as maxima, whereas the planetesimal impacts in a residual magma ocean result in a decreasing equilibration pressure with increasing time. Planetesimal impacts onto a solid surface result in a plateau at a pressure corresponding to the depth of the next magma ocean or a zero pressure when the material is late veneer and there is no subsequent giant impact.
As expected, the fast-cooling models without an insulating atmosphere mainly show plateaus because magma oceans crystallize before the next planetesimals impact the surface of the planet. The equilibration pressures for almost all collisions are then dominated by the initial depths of the melt pools that are created by embryo-embryo collisions (i.e. giant impacts). These pressures are mostly higher than the equilibration pressures in the slow-cooling models, where many planetesimals equilibrate in partially crystallised magma oceans.
The differences in equilibration pressures between the slow-cooling and fast-cooling models are greatest when the time between embryo-embryo collisions is just long enough for the magma ocean to crystallise in the slow-cooling models. This results in many small impacts in a partially crystallised magma ocean as can be seen in Figure 5b. In models with smaller embryo masses (model 4:1-0.25-7, Figure 5a) there are more embryos and therefore more embryo-embryo collisions. This results in embryo impacts into partially-crystallised magma oceans, resulting in magma oceans that reach almost to the core-mantle-boundary. Equilibration pressures become similar to equilibration pressures in the fast cooling model where most of the equilibration takes place at or near the core-mantle boundary.
The dotted lines in Figure

Discussion
The planetesimals in the N-body accretion models used here are unrealistically large and are actually tracers that represent swarms of much smaller bodies. However, since small planetesimal impacts do not create enough melt to cause a significant equilibration event, these bodies equilibrate either in a pre-existing magma ocean or at the bottom of the melt pool created by a subsequent giant impact. The exact size and timing of the planetesimal impactors therefore have little or no influence on the equilibration pressure.
Because of a number of simplifying assumptions, the results presented here must be considered to be preliminary. For example, the simple analytical model of melt production may not be accurate (Marchi et al., 2014) and does not take into account impacts onto existing magma oceans. In addition, the cooling models, based on constant time-averaged heat fluxes, are extremely simple. At this stage, our main intention is to emphasize the potential of this novel approach for modelling multistage core formation.

Conclusions
Our results show that metal-silicate equilibration pressures during core formation are likely to depend strongly on the rate of magma ocean crystallisation compared to the time interval between impacts. On a planet with a dense atmosphere, or one which develops a solid crust, magma oceans will cool slowly and numerous small impactors will equilibrate at the bottom of this magma ocean. If the liquid magma radiates heat directly to space, most magma oceans will have crystallised before the next impact and then planetesimal material generally only equilibrates when the next giant impact creates a large and deep melt volume.
If Earth lacked an insulating atmosphere and/or a solid crust during accretion, metal-silicate equilibration would have mostly taken place at core-mantle boundary pressures. This is unlikely because, based on the partitioning of siderophile elements between the core and mantle, average equilibration pressures must have been considerably lower (e.g. Wood et al., 2006;Righter, 2011;Rubie et al., 2015aRubie et al., , 2015b. Thus, based on the present results, the Earth most likely had magma oceans with lifetimes comparable to or longer than the interval between impacts. A single, long-lived magma ocean is a possibility, especially if a solid crust developed. This possibility, however, is not favoured based on noble gas isotopic arguments which have been used to argue for several generations of magma oceans (Tucker and Mukhopadhyay, 2014). Conversely these isotopic results are consistent with an Earth possessing a thick, insulating atmosphere for much of its accretion history.
The accretion history also strongly influences the equilibration pressures for slow-cooling models. The          Table 3. Starting parameters used for the three N-body accretion simulations (number and masses of embryos and planetesimals respectively). Collision histories of these three simulations are used for calculations of impact-induced melting in this study.