A Giant Impact Origin for the First Subduction on Earth

Hadean zircons provide a potential record of Earth's earliest subduction 4.3 billion years ago. It remains enigmatic how subduction could be initiated so soon after the presumably Moon‐forming giant impact (MGI). Earlier studies found an increase in Earth's core‐mantle boundary (CMB) temperature due to the accumulation of the impactor's core, and our recent work shows Earth's lower mantle remains largely solid, with some of the impactor's mantle potentially surviving as the large low‐shear velocity provinces (LLSVPs). Here, we show that a hot post‐impact CMB drives the initiation of strong mantle plumes that can induce subduction initiation ∼200 Myr after the MGI. 2D and 3D thermomechanical computations show that a high CMB temperature is the primary factor triggering early subduction, with enrichment of heat‐producing elements in LLSVPs as another potential factor. The models link the earliest subduction to the MGI with implications for understanding the diverse tectonic regimes of rocky planets.


Introduction
As planetary exploration continues, Earth remains unique as the only planet that is known to operate in the plate tectonic regime.In plate tectonics, subduction is the major process affecting the geodynamic and geochemical evolution of our planet, but when and how the first subduction occurred remains contentious.Although their interpretation is debated (Kemp et al., 2010), analyses of Hadean detrital zircons show geochemical signals consistent with subduction as early as 4.3 Ga (Valley et al., 2002;Watson & Harrison, 2005).Thus, the earliest subduction may have followed less than ∼200 million years after the Moon-forming giant impact (MGI) (Canup & Asphaug, 2001).The collision of a protoplanet called Theia with the proto-Earth at ∼4.51 Ga is expected to largely reset the initial conditions of Earth's evolution, creating a need for systematic and quantitative exploration of the planetary state established by MGI and its lasting influence on Earth's tectonic evolution.Here we propose that the thermal and chemical structure established by a canonical MGI creates favorable conditions for the initiation of subduction in Earth's earliest period, exemplifying the pivotal influence of initial conditions set by giant impact processes for the tectonic evolution of terrestrial planets.
A number of mechanisms have been proposed to initiate early subduction (Stern & Gerya, 2018), but all except three of these models require a pre-existing weakness in the lithosphere.The exceptions are grainsize evolution (Foley et al., 2014), impact-driven (Borgeat & Tackley, 2022;O'Neill et al., 2017) and plume-induced (Gerya et al., 2015;Ueda et al., 2008) subduction initiation.However, large (>700 km) and energetic impacting bolides capable of initiating subduction are rare and only occurred very early in Earth's history (Marchi et al., 2014; • The mantle thermochemical structure left by the Moon-forming impact triggers strong plumes that may have initiated the first subduction • The strong mantle plumes may rise from a heated core or from large lowshear velocity provinces enriched in heat-producing elements • Early giant impact events may have a profound influence on the diverse tectonic regimes of rocky planets

Supporting Information:
Supporting Information may be found in the online version of this article.O'Neill et al., 2017), which may coincide with the period of magma ocean solidification.Recently, plumeinduced subduction has been increasingly proposed as a viable mechanism for subduction initiation in the early Earth (Baes et al., 2021).Nevertheless, it is unclear how to generate strong mantle plumes in a vigorous convecting Hadean mantle (Korenaga, 2008).
The MGI is typically thought violent enough to melt a substantial portion or even the entire mantle, especially in scenarios involving high energy and high angular momentum (Lock & Stewart, 2017;Nakajima & Stevenson, 2015).Using two different giant impact computational methods with improved equations of state (Stewart et al., 2020) and at high (Deng et al., 2019) and ultra-high (Kegerreis et al., 2022) resolution, our recent work (Yuan et al., 2023) demonstrates that the lower half of Earth's mantle would have remained mostly solid after a canonical MGI.It also suggests that large intact domains of Theia's likely Fe-rich mantle are candidates of the two seismically-observed large low-shear velocity provinces (LLSVPs).Furthermore, the concomitant core formation process during the MGI may substantially heat the core to raise the temperatures at the core-mantle boundary (CMB) (Canup, 2004;Deng et al., 2019).This unique thermochemical structure, featuring a solid lower mantle and a hotter CMB resulting from the MGI, provides ideal conditions for the formation of strong mantle plumes, a scenario less likely when the mantle and core are in equilibrium.Accretion scenarios built from only smaller-scale impacts dissipate more of their gravitational potential energy in the mantle and yield a more stable initial thermal structure.Here, we test the hypothesis that such strong mantle plumes can be developed in the early Hadean, weakening the lithosphere and eventually causing subduction initiation (Figure 1).We model this scenario using 2-D and 3-D whole-mantle thermo-chemical-mechanical models with a visco-elasto-plastic rock rheology.

Governing Equations and Material Description
The computations solve the conservation equations of mass, momentum, and energy:  (Deng et al., 2019;Yuan et al., 2023).After solidification of the upper magma ocean within several million years (Hamano et al., 2013;Nikolaou et al., 2019), we assume a proto-lithosphere gradually formed over 50-100 million years (c), which was then destabilized and segmented by a mantle plume from an large low-shear velocity province made of the Theia mantle remnant (Yuan et al., 2023) (d).
where u is velocity, σ′ is deviatoric stress, P is pressure, T is temperature, ρ is density, g is gravitational acceleration, ẑ is the vertical unit vector, c p is heat capacity, k is thermal conductivity and H is heat production.
An incompressible Maxwell-rheological model is adopted to reflect viscoelastic behavior: where εij is the strain rate tensor, σij is the stress rate tensor, K e is the shear modulus for elasticity, and η is the viscosity.Non-Newtonian and temperature-dependent rheology is specified by: where E is the activation energy, R is the gas constant, and n is the exponent for shear thinning in dislocation creep.Material yielding is included by imposing an upper limit on the stress, thus the effective viscosity becomes where τ y is the yield stress and εII is the second invariant of the strain rate tensor.The effective viscosity is restricted to the range between 10 18 and 10 24 Pa s.The yield stress follows the Drucker-Prager yielding criteria with an upper cut-off: where τ y is the yield stress, C and μ y are cohesion and friction coefficient, p is pressure, and τ y0 is the maximum yield stress.
Similar to previous studies (Li & Gurnis, 2022), weakening processes are approximated by reducing the yield stress with plastic strain, following a two-stage process: prior to a strain saturation, C and μ y linearly decrease as follows: where the variables are the accumulated plastic strain ε P , the reference plastic strain ε P0 that determines the rate of weakening (thus, the inverse of ε P0 is the weakening rate), minimum yield stress τ yf , initial cohesion C 0 , and initial friction coefficient μ y0 .
Three materials are present within the model domain: a 30-km-thick mafic crust on top that can metamorphose to eclogite at depth (Hacker et al., 2003), LLSVP material, and background mantle.

Model Setup in 2D
Most 2D models have a domain that is 5,780 × 2,890 km with 256 × 128 linear, quadrilateral elements, with 30 Lagrangian particles in each element.We have verified the validity of 2D results with a higher-resolution grid of 512 × 256 elements, and in longer aspect ratio of 14,450 × 2,890 km with 640 × 128 elements.A free-slip boundary condition is applied to all model domain boundaries.The temperature at the upper boundary is maintained at 273 K, and the bottom boundary is fixed at a varying temperature (5273 K down to 3273 K).The lithosphere's maximum yield strength in 2D is set between 80 and 180 MPa.As seismic observations require that LLSVPs have a distinct bulk modulus (Tan & Gurnis, 2005), we use a depth-dependent density profile within the LLSVP.The values of physical parameters are listed in Table S1 of Supporting Information S1.

Model Setup in 3D
Three-dimensional plume-induced subduction computations are used to verify the 2D results and consist of a Cartesian domain 2,890 km deep, 2,000 × 2,000 km horizontally.We use a resolution of 128 × 80 × 80 elements with a uniform grid spacing in each coordinate direction.Each element initially has 30 material points randomly located in each element.Free slip conditions are applied to all the boundaries.Maximum yield stresses of lithosphere between 60 and 100 MPa were tested in 3D.

Results
Following the MGI, a surface magma ocean extending to a certain depth is expected.Tracking magma ocean solidification is beyond the scope of this study, but it is generally accepted that the solidification process would be mostly complete within several million years (Hamano et al., 2013;Nikolaou et al., 2019).For simplicity, our models begin 50 Myr after the magma ocean solidification; the initial state is derived from a 2D thermomechanical model with a proto-lithosphere defined by a half-space cooling model with a prescribed age of 50 Ma.
Most model runs include an intrinsically dense layer of Theia mantle materials above the CMB.Current estimates of the present-day CMB temperature that account for a lower boundary layer range from ∼3473 to 4573 K (Deschamps & Cobden, 2022), but because of Earth's secular cooling, the Hadean CMB temperature should have been higher.Moreover, a hot CMB is anticipated after the MGI due to the deposition of Theia's iron core, as shown in previous MGI simulations (Canup, 2004).Consequently, for most cases we used a higher CMB temperature (5273 K) and a boundary-layer contrast of 1273 K between the CMB and the initial LLSVP temperature.Such high CMB temperature may lead to partial or complete melting in the boundary layer, but temperature should drop quickly below the solidus of likely lower mantle compositions above the CMB.A thin partial melting zone within the mantle does not substantially alter the dynamics, since the fully molten core sets the lower free slip boundary condition in all cases.Furthermore, we note that similar or higher CMB temperatures were used in previous convection simulations (Nakagawa & Tackley, 2010), and we confirm our results with lower CMB temperatures (4273 and 3773 K).We mostly used a 80 MPa yield strength for the proto-lithosphere (Rudi et al., 2022).
Model runs with this initial configuration consistently show that a mantle plume nucleates atop the LLSVP-like thermochemical pile, ascends through the mantle, ruptures the lithosphere, and initiates subduction after ∼120 Myr (Figure 2).In detail, typical model evolution can be divided into several stages.First, the initial layer of intrinsically dense material (Figure 2a) at the CMB concentrates into discrete piles due to a combination of strong CMB heating and its depth-dependent density (Tan & Gurnis, 2005).Meanwhile, the lithosphere grows thicker due to inefficient cooling while in a stagnant-lid convection, although the mode of magmatism may affect a planet's tectonic regime (Lourenço et al., 2018).At 99 Myr, a strong mantle plume develops from the top of the central pile (Figure 2b).This plume quickly rises to the upper mantle, becoming narrower due to the viscosity decrease in the transition zone.As the plume approaches the surface, the lithosphere is locally weakened through yielding above the plume and, at 109 Myr, the persistent localized thinning leads to plume penetration through the lithosphere, forming a high-temperature, low-viscosity wedge structure (Figure 2c).Meanwhile, the proximal lithosphere thickens around the plume head, resulting in two lithospheric-scale shear bands on each side (Figure 2c).The plume then starts to spread near the surface, overcoming the yield strength of the surrounding lithosphere segments and forming shear zones along the subducting plate boundary that promote self-sustaining and retreating subduction at 123 Myr (Figure 2d).
Our reference case above shows that strong plumes can arise from the LLSVPs and initiate subduction when the lithosphere's strength is less than 100 MPa, a threshold yield stress consistent with previous geodynamic studies (Moresi & Solomatov, 1998;Van Heck & Tackley, 2008).Here, we aim to explore the maximum lithosphere strength that the MGI-caused mantle plumes can breach to initiate subduction and assess the influence of LLSVPs on this capability (Figure 3).Initially, we increase the lithosphere yield strength from 80 to 100 MPa and find that it becomes more challenging to induce subduction; subduction only occurs on one side of the plume head (Figure 3a).When the lithosphere yield stress is further increased to 150 MPa, plumes from LLSVPs do not penetrate through the lithosphere to induce subduction but instead spread at the base of the lithosphere (Figure 3b).Next, for the same calculation after removing LLSVPs materials, we find that strong plumes rising from the CMB can initiate subduction, as shown in Figure 3c for a case with a yield stress of 150 MPa, and in another case with a yield stress of 180 MPa (Figure 3d).
As studied in previous research on plume-induced subduction, temperature, size, and buoyancy of plumes play a major role in subduction initiation.Therefore, we systematically explore the influence of CMB temperature, which significantly affects all these factors in models where plumes are self-consistently generated.We lower the CMB temperature to 4773 K (Figure 3e) and 3773 K (Figure 3f), and find an increase in the time delay for plumeinduced subduction initiation (Figures 3i and 3j) due to the lower buoyancy of the plumes.Nonetheless, subduction still occurs as long as CMB temperature is ≥3773 K.When CMB temperature is 3273 K, no subduction is initiated (Figures 3g and 3k).However, if we assume that LLSVP materials are enriched in heat-producing elements for the same calculation (Citron et al., 2020), the results indicate that plumes rising from hotter LLSVPs can still generate subduction (Figures 3h and 3l).Computations with different rheologies, compositions, lithosphere thicknesses, and model geometries showed that these factors have a second-order influence on outcomes compared to the maximum yield strength and CMB temperature.Cases with different rheologies show that a smaller viscosity of the LLSVPs (Figures S1a in Supporting Information S1), a larger reference viscosity (Figures S1b in Supporting Information S1), and a more strongly temperature-dependent viscosity (Figures S1c in Supporting Information S1) do not significantly alter the plume-induced subduction process, except for the onset time (77-320 Myr)-in general, higher mantle viscosity leads to later plume-induced subduction initiation.We considered cases with lower (Figures S1d in Supporting Information S1) and higher (Figures S1e in Supporting Information S1) LLSVP density than the reference case, as in our earlier work (Yuan andLi, 2022a, 2022b), as well as a more compressible (more depthdependent density) LLSVP case (Figures S1f in Supporting Information S1).We find the overall dynamics of these computations closely resemble those of the reference case.Different ages of the lithosphere from 25 to 100 Myr (Figures S2a-S2c in Supporting Information S1) and different weakening rates (Figures S2d and S2e in Supporting Information S1) were also considered, and their results remain similar to the reference case.Models with higher resolution and periodic boundary conditions (Figure S2f in Supporting Information S1), and a different aspect ratio (5:1) (Figure S3 in Supporting Information S1) do not significantly differ in the results.
We further test our model using 3D whole-mantle thermomechanical models, as plume-induced subduction necessitates a three-dimensional perspective.As in the 2D geometry, we observe that strong mantle plumes can be generated from the top of an LLSVP-like pile (Figures 4a and 4b).As the plume rises and pierces the lithosphere (Figure 4c), it spreads over broken lithosphere segments in a circular area, pushing them downward into the mantle by overcoming their yield strength.As the plume continues to spread, the surface shear zones extend deeper along the subducting plate boundary, which facilitates subduction and retreat of the circular subduction zone (Figure 4d).The evolution resembles previous 3D models of plume-initiated subduction in which the plume For panels (e-g), the CMB temperature is 4773 K (e), 3773 K (f), 3273 K (g).Panel (h) has the same parameters as panel (g) except that the LLSVP materials (Radi_LLSVP) are enriched in heat-producing elements.
was imposed (Gerya et al., 2015).In another 3D case, we find that a slightly higher yield strength for the lithosphere (80 MPa) does not significantly influence the results (Figure S4 in Supporting Information S1), but subduction is not initiated if the lithosphere strength is 100 MPa as shown in another case (Figure S5 in Supporting Information S1).However, if we simply remove LLSVP materials, plumes rising from the CMB can still cause subduction initiation (Figure S6 in Supporting Information S1), similar to what we identified in our 2D models.In the next 3D case, we lower the CMB temperature to 3773 K, the plume does not penetrate the lithosphere or initiate subduction (Figure S7 in Supporting Information S1).Thus, the systematic trends we identified in 2D computations-that a higher CMB temperature is crucial in plume-induced subduction when the lithosphere strength is not very strong-remain unchanged in 3D, though the threshold values may differ due to effect of geometry, plume size and morphology.

Discussion and Conclusion
Hadean detrital zircons and anomalies in isotope ratios derived from short-lived radioactive decay systems support the notion that subduction operated in the Hadean (Hyung & Jacobsen, 2020;Valley et al., 2002;Watson & Harrison, 2005).Lithospheric dripping, suggested as an alternative mechanism in early Earth tectonics, may arise from plume-lid interaction (Fischer & Gerya, 2016) or high magmatic intrusion/extrusion ratios (Lourenço et al., 2020).However, lithospheric dripping represents non-plate-tectonic subduction, and a recent experimental study (Hastie et al., 2023) lends support to subduction over the dripping mechanism for the formation of Earth's earliest continental crust.Their results indicate that the absence of garnet in the lower crust would complicate the initiation of the gravitational instability typically necessary for crustal drips and delamination, though not necessarily in grain-damage scenarios.Other studies suggest an early surface ocean could thermally weaken the lithosphere, aiding subduction initiation (Korenaga, 2007).Yet, advanced models on crack formation and serpentinization effects leading to subduction are still lacking.Plume-induced subduction presents an appealing mechanism for the earliest subduction initiation, as it does not necessitate pre-existing weak zones or other external forces (Baes et al., 2021).It also holds the potential to trigger plate tectonics (Gerya et al., 2015), but developing such strong mantle plumes in an expected hot Hadean Earth is difficult.Our computations show that is a path to developing such plumes, and it is closely linked to the presumed MGI, which happened <150 Ma before the oldest dated terrestrial zircons.Following the earliest subduction, the proto-lithosphere can usually be recycled globally, depending on lithosphere strength and plume-slab interactions, with subsequent periods of sluggish lid convection where the subducting slab periodically tears off, similar to the tectonic mode inferred for early Earth (Foley, 2020).
The present model builds on results of high-resolution MGI simulations that yield an end-state with a largely solid lower mantle (Yuan et al., 2023) and an expected increase in CMB temperature due to the gravitational deposition of Theia's core during the agglomeration process (Asphaug, 2014).Although the exact CMB temperature increase varies in different MGI models, it should be higher than that resulting from the core-formation process during previous smaller impacts (Nimmo, 2022).Pebble accretion may also increase CMB temperatures, while recent studies suggest this mechanism is so efficient that it could melt the entire mantle (e.g., Johansen et al., 2023), leaving no solid mantle to produce strong mantle plumes.Our results show that the unique context caused by the MGI can give rise to strong mantle plumes capable of inducing subduction even when the maximum lithosphere yield strength is >150 MPa, which is noticeably higher than inversions based on present-day plate tectonics (Rudi et al., 2022).In addition, our recent work suggests parts of Theia's iron-rich mantle may persist over time and form the seismically-observed LLSVPs (Yuan et al., 2023).Here, we assess the influence of dense LLSVPs on subduction initiation.Our results suggest that plumes rising from the CMB, rather than from LLSVPs, are more capable of generating subduction and allow the mechanism to operate under a stronger lithosphere.In addition, the timing of subduction due to the presence of LLSVPs is shown to be slightly later compared to the case where plumes rise from the CMB, consistent with an earlier study (Kreielkamp et al., 2022).Nevertheless, our results suggest that plumes rising from LLSVPs are still able to generate subduction to match the Hadean zircon records as long as the CMB temperature is ≥3773 K.As the CMB temperature further decreases over time, plumes may not be strong enough to induce new subduction (Figure 3e).However, if LLSVP materials are enriched in heatproducing elements, LLSVP-sourced plumes may continue to trigger transient episodes of subduction (Figures 3h  and 3l).The enrichment of heat-producing elements in LLSVPs can be created by Earth's internal differentiation process (Citron et al., 2020).If LLSVPs are made of Theia's mantle (Yuan et al., 2023), they are likely also enriched in heat-producing elements, given that the canonical MGI suggests Theia resembles Mars, which has been estimated to have higher heat-producing elements than Earth (Table S2 in Supporting Information S1).Thus, exploring the potential link between ancient plume-induced subduction sites and the location of LLSVPs may provide a powerful means to reconstruct global plate motions beyond 200 Ma.
Earth's unique convective regime, mobile-lid plate tectonics, is distinctly absent on its twin planet, Venus.This is true even though there are potential sites of subduction, at least on modern Venus, thought to be plume-induced, which do not apparently mature into full mobile-lid tectonics.This may be due to the high surface temperature of Venus, which allows rapid healing of damage along shear zones (Bercovici & Ricard, 2014) and arrests the formation of continuous plate boundaries (Davaille et al., 2017).Alternatively, it may be harder to convert Venus crust to eclogite, hindering development of planetary-scale subduction (Chen et al., 2022).In addition to the various ongoing differences between the two planets, another significant distinction may be the initial conditions.The absence of a moon and a core-generated dynamo suggest that Venus never experienced late Moon-forming giant impacts, based on planetary accretion models (Jacobson et al., 2017).Our proposed link between the MGI and subduction initiation on Earth thus implies that the conditions for nucleation of early subduction on Earth may have been absent on Venus.Given that a planet's tectonic development is history-dependent (O'Neill et al., 2016), the early giant impact events may play a crucial role in shaping their divergent evolutionary paths.

Figure 1 .
Figure 1.A schematic illustration of the path from the Moon-forming giant impact (MGI) to plume-induced subduction in the early Hadean.A canonical MGI (a) leads to a two-layered mantle structure with a mostly solid lower mantle and an elevated core-mantle boundary temperature (b)(Deng et al., 2019;Yuan et al., 2023).After solidification of the upper magma ocean within several million years(Hamano et al., 2013;Nikolaou et al., 2019), we assume a proto-lithosphere gradually formed over 50-100 million years (c), which was then destabilized and segmented by a mantle plume from an large low-shear velocity province made of the Theia mantle remnant(Yuan et al., 2023) (d).

Figure 2 .
Figure 2. Temporal evolution for the temperature (left column) and viscosity (right column) fields of the reference case showing a LLSVP-sourced plume inducing subduction initiation.(a) Initial set up of the model; (b) Plume rising from the LLSVP; (c) Plume penetrates the lithosphere; (d) Slab retreat and melting of crust.The yellow contours in temperature field marks the boundary of different materials including crust, lithosphere, molten materials, and large low-shear velocity province materials.

Figure 3 .
Figure 3. Snapshots of viscosity and temperature field of eight representative cases with varying lithosphere yield stress (τ y0 ) and core-mantle boundary (CMB) temperatures (T_cmb).Each panel from (a-h) represents a different calculation, with panels (i-l) corresponding to the viscosity field for cases in panels (e-h).In panels (a-b), the lithosphere yield strength is 100, 150, 150, and 180 MPa, respectively.Note there are no large low-shear velocity province (LLSVP) materials in panels (c-d).For panels (e-g), the CMB temperature is 4773 K (e), 3773 K (f), 3273 K (g).Panel (h) has the same parameters as panel (g) except that the LLSVP materials (Radi_LLSVP) are enriched in heat-producing elements.