Primordial Earth mantle heterogeneity caused by the Moon-forming giant impact

TL;DR

This study proposes that the Earth's mantle's primordial heterogeneity, caused by the Moon-forming impact, explains the compositional similarities between Earth and Moon and aligns with geochemical and geophysical observations.

Contribution

It introduces a new hydrodynamical model showing impact-induced stratification of Earth's mantle, explaining the preservation of primordial heterogeneity and Earth-Moon compositional similarities.

Findings

Impact leads to a stratified mantle with distinct layers.

Primordial heterogeneity can be preserved over billions of years.

Earth-Moon system may have a chondritic bulk composition.

Abstract

The giant impact hypothesis for Moon formation successfully explains the dynamic properties of the Earth-Moon system but remains challenged by the similarity of isotopic fingerprints of the terrestrial and lunar mantles. Moreover, recent geochemical evidence suggests that the Earth's mantle preserves ancient (or "primordial") heterogeneity that predates the Moon-forming giant impact. Using a new hydrodynamical method, we here show that Moon-forming giant impacts lead to a stratified starting condition for the evolution of the terrestrial mantle. The upper layer of the Earth is compositionally similar to the disk, out of which the Moon evolves, whereas the lower layer preserves proto-Earth characteristics. As long as this predicted compositional stratification can at least partially be preserved over the subsequent billions of years of Earth mantle convection, the compositional similarity between the Moon and the accessible Earth's mantle is a natural outcome of realistic and high-probability Moon-forming impact scenarios. The preservation of primordial heterogeneity in the modern Earth not only reconciles geochemical constraints but is also consistent with recent geophysical observations. Furthermore, for significant preservation of a proto-Earth reservoir, the bulk composition of the Earth-Moon system may be systematically shifted towards chondritic values.

Tables

Table 1. Table 1: Parameters and results of impact simulations

Run	$M_{tar}$	$M_{imp}$	$b$	$\frac{v_{imp}}{v_{esc}}$	$\frac{L_F}{L_{EM}}$	$\frac{L_D}{L_{EM}}$	$\frac{M_D}{M_L}$	$\frac{M_{planet}}{M_E}$	$\frac{M_U}{M_{mantle}}$	$F_{U,tar}$	$F_{D,tar}$	$\frac{M_{Fe}}{M_D}$	$\frac{M_M}{M_L}$
1	0.85	0.12	0.71	1.00	0.91	0.16	0.87	0.94	0.46	0.84	0.38	0.16	0.73
2	0.85	0.12	0.73	1.00	1.18	0.35	1.48	0.95	0.50	0.86	0.40	0.04	1.48
3	0.85	0.16	0.71	1.00	1.30	0.32	1.70	0.98	0.49	0.80	0.30	0.04	1.5
4	0.85	0.16	0.73	1.00	1.30	0.23	1.22	0.98	0.43	0.79	0.39	0.12	1.08
5	0.90	0.16	0.71	1.00	1.32	0.32	1.98	1.0	0.52	0.81	0.47	0.01	1.19
6	0.90	0.16	0.73	1.00	1.53	0.38	1.64	1.0	0.52	0.80	0.34	0.05	1.64
\hdashline7	0.85	0.2	0.574	1.20	1.33	0.29	1.63	0.99	0.59	0.79	0.52	0.08	1.26
8	0.85	0.2	0.574	1.25	1.26	0.26	1.45	0.98	0.62	0.80	0.47	0.09	1.14
9	0.90	0.2	0.537	1.20	1.30	0.05	0.30	1.06	0.55	0.79	0.64	0.01	0.20
10	0.90	0.2	0.537	1.25	1.32	0.13	0.81	1.04	0.60	0.79	0.62	0.03	0.48
11	0.90	0.2	0.537	1.30	1.30	0.20	1.20	1.04	0.54	0.79	0.55	0.12	0.79
12	0.90	0.2	0.574	1.15	1.38	0.14	0.88	1.04	0.54	0.79	0.59	0.01	0.51
13	0.90	0.2	0.574	1.20	1.34	0.23	1.23	1.02	0.56	0.79	0.57	0.08	1.07
14	0.90	0.2	0.574	1.25	1.30	0.24	1.30	1.02	0.59	0.80	0.54	0.09	1.30
15	0.90	0.2	0.574	1.30	1.20	0.23	1.18	1.02	0.50	0.79	0.44	0.16	1.18

Full text

{affiliations}

Center for theoretical Astrophysics and Cosmology, Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland

Institute of Geophysics, ETH Zurich, Sonneggstrasse 5, 8092 Zurich, Switzerland

Institute of Geochemistry and Petrology, ETH Zurich, Clausiusstrasse 25, 8092 Zurich, Switzerland

As the only planet in our solar system, the Earth is orbited by a single and massive moon. The leading theory for the formation of the Earth-Moon system with its high angular momentum involves a giant impact followed by lunar aggregation from the impact debris disk[1]. The canonical giant impact model involves a graze-and-merge impact, in which a Mars-sized body (or “Theia”) collides with the proto-Earth at an oblique angle at roughly the escape velocity of the system[2]. In this model, however, Theia contributes a larger fraction of silicates ($\sim$70% by mass) to the proto-lunar disk than to the proto-Earth. Unless Theia and the proto-Earth had almost the same isotopic composition, this imbalance is at odds with the strong isotopic similarity of the Earth’s and lunar mantles, e.g., in terms of oxygen[15] and titanium[16].

One way to reconcile this compositional similarity involves the post-impact re-equilibration of the Earth and the Moon-forming disk[17]. This model, however, is unable to explain the isotopic similarity of the Earth and Moon in highly refractory elements, for example, titanium[16]. More recently, several alternative giant-impact models have been proposed. A near equal-mass “Sub-Earth” impact[18] or the disruption of a fast-spinning Earth (close to self-breakup) by a small impactor[19] can indeed explain the isotopic similarity. However, the proposed solutions are low-probability events: equal-mass collisions are highly unlikely in the late stage of planetary accretion[8], and planetary embryos hardly reach spins that are close to self-breakup in hydrodynamic simulations[20]. Both models also predict an angular momentum that is too high for the early Earth-Moon system, and the mechanisms proposed to remove the excess angular momentum exclusively work in a narrow tidal parameter range[21]. Alternatively, multiple-impact models have been suggested for lunar origin[20], certainly broadening the favourable parameter space compared to single-impact models, and giving rise to mixing through mergers of moonlets with different isotopic composition. However, the dynamics of moonlets are highly uncertain, and primitive moonlets might be lost during repeated impacts pre-dating the Moon-forming stage [22].

Here, we explore the mixing state of the Earth-Moon system in low-velocity, low-angular-momentum impact scenarios. We consider the canonical model[2] and the hit-and-run model[23](see Figure 1a, for these $v_{imp}<1.3v_{esc}$, where $v_{imp}$ is the impact velocity and $v_{esc}$ is the escape speed), which are both high-probablity impact configurations [24] and lead result naturally in the current angular momentum of the Earth-Moon system. We apply a new hydrodynamical Lagrangian method: Meshless Finite Mass (MFM) (see Methods subsection 2). This method is better at capturing fluid mixing and at resolving the core-mantle boundary than the widely-used smoothed-particle hydrodynamical approach[7].

In our simulations of such impacts, a strong shock propagates almost perpendicular to the line connecting the centers of the two impacting bodies after the first contact. In both the canonical and hit-and-run scenarios, the part of the impactor that can avoid direct collision is sheared into a spiral structure which, shortly afterwards, collapses into clumps. The clumps re-impact the highly distorted target to eject some additional target material into the circum-planetary debris disk (Supplementary Video; Figure 1b). Figure 1c shows the predicted fraction of target material in the post-impact body, $F_{tar}$ as a function of normalized enclosed mass (which corresponds to the radius, within which a given planetary mass fraction is enclosed). Even though mixing is more efficient with our new method than in previous studies, most of the impactor’s silicates still remain in the outer layer of the post-impact target mantle (low $F_{tar}$). This prediction is explained by inefficient transfer of angular momentum during the impact (Methods subsection 1). Focused shock heating in the outer layer (Supplementary Video) results in a steep entropy profile (Figure 1d) through most of the mantle. In particular, there is a distinct entropy jump at radius $R$, or at a normalized enclosed mass of $\sim$0.7 $M_{\oplus}$. In terms of mass, this $R$ corresponds to $\sim$1000 km depth in the present-day mantle.

For the metal core, our models also predict a steep entropy and compositional profile in the aftermath of the impact. The entropy of the deep core even drops below its initial condition of 1200 J/kg/K. This is well explained by the effects of phase transitions near the core-mantle boundary, which can be captured by the MFM method, and results in a redistribution of energy, and entropy, from the core to the mantle [7](Methods subsection 2, Extended Data Fig. 2, 6). The impactor’s metal mostly remains near the top of the core. Depending on the conditions of metal-silicate (impactor metallic core/target silicate mantle) equilibration, this prediction may provide an explanation for the $\sim$300-km-thick compositionally stratified layer that is seismically observed at the top of the present-day outer core[25].

To test our model predictions with geochemical observations, we estimate the unknown compositions of the impactor’s (Theia) and target’s (proto Earth) mantles from the known isotopic compositions of the accessible part of the Earth’s mantle and Moon[15] (See Supplementary material for calculation using alternative data). Figure 2a shows the allowed $\Delta{}^{17}$O-difference as a function of the mass of the mantle that remains poorly homogenized (or preserved) over the age of the Earth. If the present-day mantle is fully homogenized such that any primordial stratification is completely removed, a common assumption in previous studies, Theia and the proto-Earth must have been rather similar in composition. For example, only 30 ppm $\Delta^{17}$O difference between parent bodies are allowed in our best canonical model (run 5). On the other hand, if the assumption of full homogenization is relaxed, and the predicted mantle stratification can be (partially) preserved through the present day, larger compositional differences can be reconciled with the available data. For example, considering the preservation of a compositionally-distinct domain below $R$, differences of up to 54 ppm (run 13) in $\Delta{}^{17}$O are allowed (Figure 2b), particularly for the hit-and-run models, which display larger fractions of proto-Earth silicate material in the disk than the canonical models (table 1; Figure 1). These large values are consistent with the compositional difference between potential parent bodies of the Earth-Moon system in N-body simulations of planetary accretion[8]. For realistic parent-body compositional differences[26], inefficient mixing of the Earth’s mantle through time can indeed help to resolve the geochemical similarity of the accessible Earth mantle and Moon. Even for moderate mixing across $R$, our models can critically increase the allowed compositional difference between parent bodies (see Figure 2a).

Whether post-impact heterogeneity can persist through $\sim$1 Myr of magma-ocean and $\sim$4.5 Gyrs of mantle convection is controlled by the initial compositional and thermal profiles of the mantle. Our models predict the formation of a deep magma ocean due to the energy release of the giant impact[27] (see Extended Data Fig. 6). While major-element compositions of the post-impact mantle layers above and below $R$ depend on the unknown bulk compositions of Theia and the proto-Earth, respectively, an enrichment of the lower (proto-Earth) layer in FeO and SiO2 is generally consistent with the evolving of physical conditions of multi-stage core formation during progressive planetary accretion[28, 29]. Furthermore, FeO-enrichment of the deep proto-Earth’s mantle may have been promoted by compositional fractionation (and subsequent overturn)[30] during any magma-ocean episode(s) that predate(s) the Moon-forming impact. Even just a slight FeO-enrichment of the lower layer is sufficient to promote stable stratification through various magma-ocean stages. A long lasting stratification should be favoured by the entropy gradient across the post-impact mantle predicted by our models as the observed sub-adiabatic entropy gradient should prevent redistribution by convection.

After the final magma-ocean episode that follows the giant impact, mixing during long-term solid-state mantle convection is controlled by the density and viscosity contrasts between the two layers. Primordial heterogeneity can survive mantle stirring as blobs over a range of spatial scales (meters to terameters)[31, 32]. Intrinsically high densities and viscosities of the primordial deep-mantle layer as sustained by an enrichment in FeO and SiO2 of the proto-Earth’s mantle (see above) impede efficient across $R$. For example, ref. [28] predicts a molar Mg/Si of $\sim$0.98 for the proto-Earth’s (lower) mantle, corresponding to a predominant abundance of the high-viscosity mineral bridgmanite (Mg,Fe)SiO3, implying poor mixing. While some degree of whole-mantle mixing is indicated by the sinking of a subset of subducted slabs of oceanic lithosphere through the entire mantle, the stagnation of other slab segments[11] and the deflection of upwelling plumes[12] at $\sim$1000 km depth (i.e. about radius $R$) is indeed consistent with restricted mixing. Sharp seismic-velocity contrasts at similar depths support this interpretation, and provide direct evidence for large-scale compositional mantle heterogeneity[13, 14]. The preservation of primordial noble gases[10], e.g. the large missing budget of argon[9], provides complementary evidence for incomplete homogenization of primordial mantle reservoirs. 182W isotopic evidence[4, 5] requires that least a subset of the preserved heterogeneity predates the moon-forming impact, and thus reflects proto-Earth compositions.

According to our results, the preservation of a significant fraction of primordial mantle heterogeneity through the age of the Earth can explain the isotopic similarity between the Earth’s and lunar mantles for a relatively wide range of parent-body compositional differences (Figure 2). Preservation of SiO2-enriched heterogeneity further helps to balance the bulk-Earth’s silica budget relative to the chondritic compositional range. Future geochemical and geophysical studies of lower-mantle composition will contribute to constrain the chemistry and origin of Earth’s parent bodies, and thus, ultimately, of the inner solar nebula, providing the means to test effectively our scenario.

{methods}

0.1 Estimation of the penetration depth of the impactor’s silicates.

In the gravity dominated regime, the interaction between the impactor’s mantle and the target’s mantle behaves like a fluid collision. One fluid element can only deliver half its momentum/angular momentum to a roughly equal mass fluid element due to the completely inelastic nature of a fluid collision. The transport of angular momentum through shock waves, primarily the contact shock at their first contact, is also inefficient because the shock is almost symmetric to the line of centres. The impactor’s mantle keeps roughly half its initial angular momentum, $0.35(1-\gamma)L_{imp}$ (confirmed by our simulations), so that it cannot sink deep into the target’s mantle. To avoid rotational instability, the specific angular momentum of the post-impact target must not decrease as radius increases (Rayleigh’s criterion); in simulations, the outer part rotates faster. In the best possible case, we assume the impactor’s silicates concentrate in a thin shell and they rotate with materials from the target’s mantle, residing outside the shell, at constant specific angular momentum (Extended Data Fig. 1). The materials beyond the shell contain the rest of the angular momentum, $(0.65+0.35\gamma)L_{imp}$. Even for $\gamma=0.15$, the impactor’s mantle cannot penetrate half the target’s mantle. A hit-and-run collision is more complicated due to interaction with the escaping part and stronger oblique shocks, but the angular momentum transport is still inefficient with $R$ a little deeper (see table 1).

0.2 Simulations and analysis.

We simulate the giant impacts using the GIZMO code[33], which is a descendant of the GADGET code[34] and its SPH method widely used in previous impact simulations). GIZMO indeed contains the legacy SPH implementation of the GADGET code, newer SPH implementations such as PSPH, and, most importantly the novel Meshless Finite Mass (MFM). Gravity is coupled to all these different hydro methods using the same treecode scheme [33]. The MFM method is an improved hydrodynamics formulation that is fundamentally different from the Smoothed Particle Hydrodynamics (SPH). GIZMO MFM does effective volume partition according to the particle distribution and then solves the Riemann problem to update the fluid variables and can be regarded as a generalized moving mesh method. It employs no explicit artificial viscosity and thus shows better conservation property than SPH[35]. MFM captures shocks and subsonic turbulence in giant impact simulations accurately, so it can simulate the mixing properly[7]. We use about 500K particles in our simulations, which is comparable to present-day high resolution simulations. We run a 2M particles simulation as a convergence test (Extended Data Fig. 2). A comparison between the standard SPH simulation, which suppresses mixing, and MFM simulation is presented in Extended Data Fig. 3.

We apply the ANEOS/M-ANEOS equation of state[36, 37] with iron comprising the core and dunite comprising the mantle. We build the initial condition for planets (30 wt% iron, 70 wt% dunite) by solving the hydrodynamic equilibrium with an isentropic profile and place the particles (computational elements) in spherical shells to represent the equilibrium profile[38]. The temperature on the planet surface is about 2000K, corresponding to an entropy of 2700 J/kg/K in the mantle and 1200 J/kg/K in the core. We run a comparison study with higher initial entropy and our results are robust concerning the initial entropy value (Extended Data Fig. 4). Note that the same version of ANEOS/M-ANEOS is implemented in SPH and MFM, so that any difference that we observe between simulations carried out with these two methods will only stem from the underlying hydrodynamical solver.

We characterize the modeled impacts by determining the disk mass/angular momentum, the predicted Moon mass, planet mass/angular momentum and the internal structure of the post-impact target. The former is done following a standard approach, i.e., bounded particles with periapsis distance larger than the equatorial radius of the planet are classified as disk particles[39, 20]. This calculation is preformed at least 40 hours after the impact, i.e. when the system saturates to a quasi-steady state [39] (Extended Data Fig. 5). We calculate the entropy profile by arithmetically averaging the entropy across spherical shells (assuming spherical symmetry). The fraction of materials from the target is calculated analogously. The iron core is slightly oblique with an equatorial radius that is slightly ($<5\%$) larger than the polar radius; hence, the profiles is not messed up at the core-mantle boundary. A counter-intuitive result involves that that core entropy can drop below its initial value for MFM, as shown by the entropy profile (e.g., Figure 1d).

This prediction by our MFM models contrasts with that of our SPH models (Extended Data Fig. 3). The entropy drop in the core for MFM is caused by a phase transition and the associated redistribution of internal energy, and entropy, from the core to the mantle, hence it has a physical origin[7]. Indeed, as the outer core, near the core-mantle boundary, melts first, entropy locally increases, forcing the inner core to decrease its own entropy in order to maintain thermodynamical equilibrium. This entropy loss outweighs the entropy gain through shocks in the central core region. In contrast, in SPH simulations, the core and mantle are separated by an artificial tensional force [40, 41, 7] (see also Extended Data Fig. 6). This force largely isolates the core, causing the core to evolve nearly adiabatically[7]. The physical heating associated with shocks is somewhat overestimated by the usage of artificial viscosity[34, 7]; as a result, the entropy slightly increases in the core for SPH. That the entropy in the central core does not decrease for SPH, is thus mostly due to a numerical artifact. The marked difference between MFM and SPH, i.e. in terms of the temperature and entropy distribution near the core-mantle boundary, is visualized in Extended Data Fig. 6.

0.3 Code availability

The latest version of the GIZMO code is made available by its author, Philip Hopkins at http://www.tapir.caltech.edu/ phopkins/Site/GIZMO.html

0.4 Data availability

The data files that support our analysis will be made available upon reasonable request.

{addendum}

We thank Martin Jutzi, Prasenjit Saha, Romain Teyssier and Qing-Zhu Yin for stimulating discussions. We acknowledge support from the Swiss National Science Foundation via the NCCR PlanetS.

H.D. conceived the idea of linking the moon formation impact to the Earth mantle heterogeneity and planned the project. C.R. built the equation of state (EOS) library and F.B. prepared the EOS lookup table. H.D. and C.R. incorporated the EOS library into the hydrodynamical code and prepared the initial conditions. H.D. ran the simulations and did the visualisation and interpretation. C.R., L.M., J.S. also helped in the interpretation. M.B. and M.M. contributed to the geodynamic and cosmochemistry argument, respectively. H.D., M.B. and M.M. prepared the manuscript. L.M. reviewed the manuscript and all authors commented on it.

The authors declare that they have no competing financial interests.

Correspondence and requests for materials should be addressed to H.D.(email:hpdeng@physik.uzh.ch)