Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38
https://doi.org/10.1186/s40645-022-00497-0

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Hadean/Eoarchean tectonics and mantle
mixing induced by impacts: a three‑dimensional
study
Xavier Borgeat1*   and Paul J. Tackley2   

Abstract
The timing of the onset of plate tectonics on Earth remains a topic of strong debate, as does the tectonic mode thatpreceded modern plate tectonics. Understanding possible tectonic modes and transitions between them is also
important for other terrestrial planets such as Venus and rocky exoplanets. Recent two-dimensional modelling studies
have demonstrated that impacts can initiate subduction during the early stages of terrestrial planet evolution—the
Hadean and Eoarchean in Earth’s case. Here, we perform three-dimensional simulations of the influence of ongoing
multiple impacts on early Earth tectonics and its effect on the distribution of compositional heterogeneity in the
mantle, including the distribution of impactor material (both silicate and metallic). We compare two-dimensional
and three-dimensional simulations to determine when geometry is important. Results show that impacts can induce
subduction in both 2-D and 3-D and thus have a great influence on the global tectonic regime. The effect is particularly strong in cases that otherwise display stagnant-lid tectonics: impacts can shift them to having a plate-like regime.
In such cases, however, plate-like behaviour is temporary: as the impactor flux decreases the system returns to what
it was without impacts. Impacts result in both greater production of oceanic crust and greater recycling of it, increasing the build-up of subducted crust above the core-mantle boundary and in the transition zone. Impactor material
is mainly located in the upper mantle, at least at the end of the modelled 500-million-year period. In 2-D simulations,in contrast to 3-D simulations, impacts are less frequent but each has a larger effect on surface mobility, making thesimulations more stochastic. These stronger 2-D subduction events can mix both recycled basalt and impactor material into the lower mantle. These results thus demonstrate that impacts can make a first-order difference to the earlytectonics and mantle mixing of Earth and other large terrestrial planets, and that three-dimensional simulations areimportant to obtain less stochastic results, and also to not over- or under-predict the amount of impactor material
mixed into the mantle and the time during which a specific tectonic regime acts.
Keywords: Initiation of plate tectonics, Hadean and Eoarchean tectonics, Late heavy bombardment, Mixing of
impactor material, Three-dimensional
1 Introduction
The time of onset of modern-style plate tectonics is still
heavily debated. Large uncertainties exist over its time of
appearance, and what tectonic mode(s) occurred before(for reviews see Korenaga 2013; Palin and Santosh 2020;*Correspondence: xborgeat@student.ethz.ch
1
Department of Earth Sciences, ETH Zurich, Zurich, Switzerland
Full list of author information is available at the end of the article

Condie 2018; Cawood et al. 2018). Other planets may
also undergo changing tectonic mode—of particularinterest here is Venus, which is proposed to have experienced episodic resurfacing, perhaps caused by global
lithospheric overturn/subduction events (e.g. Turcotte
1993; Armann and Tackley 2012).
The tectonic mode of the early Earth is often describedas stagnant lid (e.g. O’Neill and Debaille 2014), with theidea that the lithosphere was basically rigid for a long

© The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

time until some mechanism managed to break it andinitiate plate tectonics (e.g. Bercovici and Ricard 2014;
Tang et al. 2020). However, recent works have questioned
this view of the Hadean. Hot mantle temperatures would
have caused extensive melting and igneous intrusion,
as indicated by widespread reworking of the Hadeanand Eoarchean crust inferred from zircons (Kirkland
et al. 2021), and this intrusion would have weakened thelithosphere and allowed it to deform, accommodating
the lateral motion of sections of it. Such behaviour has
been demonstrated in both regional models (Fischer andGerya 2016; Sizova et al. 2010; Piccolo et al. 2019), where
it was named “plume lid” and global models, where it
was named “plutonic squishy lid” (Lourenço et al. 2020,2018), occurring for intrusion fractions of ≥ 80%. Thus,
while there is general agreement that any subduction thatstarted under early Earth conditions would not result in
long-lived subduction zones due to the weakness of thelithosphere (van Hunen and van den Berg 2008; Sizova
et al. 2010; Moyen and Van Hunen 2012), considerable
crustal deformation might well have taken place. This isalso consistent with evidence for underthrusting from >4
Gyr zircons (Hopkins et al. 2008). Large plumes mightalso have had the capability to initiate transient subduction events (Gerya et al. 2015; Piccolo et al. 2020).
During that time, the Earth was accreting around 0.5%of its current mass, the so-called late veneer, throughimpacts. Currently, there is much debate over whether
there was a late peak of impacts, the so-called Late Heavy
Bombardment (LHB) (e.g. Bottke and Norman 2017;
Marchi et al. 2014), which is proposed to have beencaused by a destabilization of the E-belt and asteroid belt
caused by a change in Jupiter’s orbit (Morbidelli et al.
2012), or whether the impact rate declined monotonically(e.g. Michael et al. 2018; Morbidelli et al. 2018; Brasser
et al. 2020). In either case, the lunar cratering record
shows that there was a flux of large impactors extending
throughout the Hadean and Eoarchean eons, and thus,
the Earth should also have been hit by such impactors
during this time period.
A number of studies have addressed the influences ofimpacts on relatively small stagnant-lid bodies such as
the Moon (Rolf et al. 2017; Ghods and Arkani-Hamed
2007) and Mars (e.g. Reese et al. 2002, 2010; Roberts
et al. 2009; Monteux and Arkani-Hamed 2014), but relatively few have studied impacts on larger planetary bodies such as the Earth and Venus, despite their potential
importance (e.g. Maruyama et al. 2018). Gillmann et al.
(2016) showed that a large impact can have a substantial
influence on the tectonics of a stagnant-lid Venus-likeplanet, causing an episode of subduction that rolls backfrom the impact location. Multiple impacts continuing
for some period of time could thus potentially mobilise

Page 2 of 19

the lithosphere for a corresponding time period. Thiswas first demonstrated by O’Neill et al. (2017), who presented two-dimensional calculations with an impact flux
based on Marchi et al. (2014) to show the importance of
impacts in early Earth tectonics. The importance of such
multiple impacts on inducing outgassing on Venus wasdemonstrated also using two-dimensional models by
Gillmann et al. (2020), although their influence on tectonics was not analysed.
However, impacts may have an exaggerated effect on
two-dimensional models because they are effectively infinite in the out-of-plane direction. Thus, we here analyse
three-dimensional simulations and compare and contrast
them to identical two-dimensional simulations in orderto determine whether the same conclusions apply in 3-D.
Additionally, we analyse the effect of multiple impacts on
mantle compositional structure, focussing on the radial
distribution of recycled basaltic crust and on the distribution of impactor material.

2 Methods summary
Full details of the physical and numerical model are givenin the later Methods section; here is a brief overview. We
couple 3-D or 2-D simulations of mantle convection witha parameterized impact model, with the rate of impacts
vs. time based on Marchi et al. (2014), which is itself
based on the sawtooth bombardment time history proposed by Morbidelli et al. (2012), with a clear spike in theimpact rate starting 4.1 Gyr ago. This impact history has
a lower peak impactor flux than the classical LHB as well
as a commonly used monotonically decreasing model(see Fig. 1 in Hopkins and Mojzsis (2015)). We assume atotal mass addition of 1.0 ∗ 1020 kg, a power-law distribution of impact sizes with an exponent -3.5 and a maximum size of 1000 km diameter. Impacts are assumed
to add heat according to a commonly used parameterisation. Additionally, the metal core of the impactor ismodelled, but other direct mechanical effects are notincluded.
Convection and lithospheric dynamics is simulatedusing StagYY (Tackley 2008) with assumed physicalparameters similar to those in recent papers (Tackley
et al. 2013; Lourenço et al. 2020). The model includes
strongly temperature-, pressure- and stress-dependent
viscosity that combines diffusion creep at low stress, dislocation creep at intermediate stress and plastic failure
at high stress. Physical properties vary with depth under
the compressible anelastic approximation. Compositionalvariations between the endmembers basalt and harzburgite are included, with partial melting that can produce
basaltic crust.
We simulate a time period of 500 Myr, from 4.2 to 3.7
Gyr before present, starting with an initial temperatureBorgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

structure approximating what is expected for the Hadean(for details see the later Methods section).

3 Results
In this section, we first present the effects of a single
impact, and then analyse the long-term evolution of thesimulations, contrasting cases with impacts and without
impacts, and in 3-D versus 2-D. Both visual aspects andquantitative measures (such as surface mobility and heatflux) are presented. Finally, we focus on compositional
aspects, in particular the radial distribution of basaltic
material and the radial distribution of impactor material.
3.1 Effect of an impact

Figure 1 shows the effect of a 400 km radius impact ina two-dimensional simulation. The heat addition generates some melt, some of which produces more crust anddepleted harzburgitic material beneath. The buoyancyof the hot, depleted, partially molten material produceshigh stresses and low viscosity, which allows subduction
to start at the edge of the impact site. The buoyant material upwells and spreads laterally, while the subduction
zones roll back. This process occurs very quickly; the first
third columns span only 31,000 years. This sequence ofPage 3 of 19

events is in line with previous simulations (O’Neill et al.
2017; Gillmann et al. 2016). Additionally, we track themetallic impactor core as well as all impactor material.
The impactor’s metal core drops very rapidly to the coremantle boundary, while other impactor material spreads
below the lithosphere.
Figure 2 shows the effect of a 400 km radius impactor on surface fields of a three-dimensional simulation.
A cross section of this major impact is shown in Fig. 3.
Somewhat different to 2-D cases, subduction does notoccur at the edge of the impactor material, but instead
within the impactor material. Some of the subductedmaterial flows along the 660km discontinuity, while somereaches the lower mantle. In contrast to the two-dimensional simulations, no continuous subduction events lasting for several million years were observed in 3-D. Thisis the first main difference observed between the twogeometries.
3.2 Long‑term evolution

Simulations with impacts can evolve quite differently
from simulations without impacts, as we first explore
visually for cases with a yield stress of 85 MPa. Figures 4,
5 and 6 show, respectively, the long-term evolution ofFig. 1 Effect of an impact in a two-dimensional simulation with yield stress 100 MPa. Plotted are (top to bottom) basalt fraction, temperature, melt
fraction, impactor fraction and metal fraction. The time after the impact is displayed at the top

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Page 4 of 19

Fig. 2 Effect of an impact in a three-dimensional simulation with yield stress 100 MPa. Plotted are (top to bottom) basalt fraction, temperature andviscosity just below the surface. The time after the impact is displayed at the top

two-dimensional simulations, three-dimensional simulations and cross sections of three-dimensional simulations. Cases with (top panels) and without (lower panels)
impacts are compared and contrasted. Furthermore, twomovies are attached showing the time-evolution of 3-D
cases.
In two dimensions, the no-impacts case displays a stagnant lid over the illustrated time. The basalt field showsa gradually thickening crust and some basalt accumulation in the transition zone due to the “basalt barrier”
mechanism (Davies 2008; Yan et al. 2020). In contrast,
the impacts case shows much evidence of short-lived
subduction events, with many slab segments in the lower
mantle. As a result, basaltic material accumulates above
the core-mantle boundary (CMB) as well as in the transition zone. The crust is thinner than in the no-impacts
case but there is a thicker layer of depleted material
below the crust.
In three dimensions the effect of impacts on tectonics
is dramatic (Fig. 5). The no-impacts case again displays
a stagnant lid the entire time, while the impacts case displays considerable tectonic activity. In the basalt field,

subduction due to impacts results in crust-free regions.
The crustal thickness is on average much lower than in
the no-impacts case. The viscosity field displays a complex mixture of low-viscosity impact features and atsome times, linear low-viscosity features corresponding
to zones of extension or subduction, i.e. tectonic plateboundaries.
In the three-dimensional cross sections (Fig. 6), it isclear that the no-impacts case is stagnant lid, while in theimpacts cases there is substantial subduction, although
the lower mantle does not contain as many slab fragments as in the 2-D case. The impacts case again shows
more basalt accumulation above the CMB and above the660 km discontinuity.
3.3 Effect of the yield stress

The key parameter controlling the tectonic mode invisco-plastic mantle convection simulations is the yieldstress, as has been much studied in the past (e.g. Tackley
2000; Moresi and Solomatov 1998; Lourenço et al. 2020).
Here we study the influence of the yield stress on surface
mobility, which is the ratio of the rms. surface velocity

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Page 5 of 19

Fig. 3 Cross sections of the frames in Fig. 2 at the same times, showing (top to bottom) basalt fraction, temperature, melt fraction, impactor
fraction and metal fraction

to the volume-averaged rms. velocity (Tackley 2000). A
ratio larger than about 1 means that the surface is mobile
and that subduction is occurring, while a value close tozero means that the surface is immobile, indicating thestagnant-lid regime (O’Neill et al. 2017; Tackley 2000).
Figure 7 shows the relationship between the percentage of time that the surface is mobile (Mobility >1) andthe yield stress for cases with both geometries and withand without impacts. Simulations without impacts display a clear and quite sharp transition from mobile lid(plate) tectonics to a stagnant-lid regime at a yield stressof around 50-80 MPa. (There is some complexity in the2-D case due to random effects.) In strong contrast, caseswith impacts maintain a mobile surface over a muchwider range of yield stress, with some periods of mobility even at yield stress exceeding 150 MPa in 2-D. Thus,
it is clear that impacts promote surface mobility for yieldstress values that would normally cause a stagnant lid
regime.
Comparing 2-D and 3-D, the trends are very similarboth qualitatively and quantitatively, but there is morestochastic variation (randomness) in 2-D. This is something that is also apparent from other diagnostics discussed later. Another difference is that the mobile-lid tostagnant-lid transition occurs at higher yield stresses in
2-D cases.
Figure 8 shows the surface mobility vs. time for thethree different yield stresses in the both dimensions, withand without impacts. The sawtooth bombardment ramps
up at 0.1 Gyr. For the lowest yield stress (45 MPa), allcases are generally mobile, although the no-impacts cases
revert to stagnant lid for part of the time. The higher
yield stress cases (85 and 125 MPa) display a clearer distinction. In these, the no-impacts cases are stagnant lid
the whole time. The with-impacts cases, in contrast, aremostly mobile after 0.1 Gyr when the impact-rate sawtooth kicks in. Comparing 2-D and 3-D cases, the 2-D
cases display more variability, with mobility sometimes
being 0 and sometimes above 2. The 3-D cases are alsoquite time-dependent but the peaks are lower than thoseof the 2-D cases. Particularly in the 85 MPa 3-D case, themobility gradually decreases with time as the impact rate
decreases, approaching 0 by the end of the simulation.
This indicates that once the impacts stop, the tectonicmode returns to what it was before the impacts - there
is no permanent kick-starting of plate tectonics. The 125
MPa 3-D case has only two bursts of mobility, emphasising that mobility is more likely in 2-D.

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Page 6 of 19

Fig. 4 Long-term evolution of two-dimensional simulations with and without impacts with a yield stress of 85 MPa. Basalt fraction, temperatureand viscosity fields are shownBorgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Page 7 of 19

Fig. 5 Long-term evolution of three-dimensional simulations with and without impacts with a yield stress of 85 MPa. Basalt fraction, temperatureand viscosity just below the surface are plotted. Animations of the long-term evolution are provided in Additional files 1, 2

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Fig. 6 Cross sections of the long-term evolution for the same cases and times as Fig. 5

Page 8 of 19

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Page 9 of 19

causing mobility that lasts for several tens of millions ofyears. When there hasn’t been an impact for some time,
the system returns to a stagnant lid state. These twopoints: fewer impacts each with more influence, explain
the difference in mobility time series between 2-D and
3-D models.
3.4 Radial distribution of basalt

Fig. 7 Per cent of time when the surface mobility is higher than 1
versus yield stress for two- and three-dimensional cases with andwithout impacts. The spread of values for 5 different impact histories
in 2-D is shown

The relationship between periods of surface mobility
and individual large impacts is shown in Fig. 9. Spikes inmobility are often associated with impacts, as expected.
Not all the ∼ 12,000 impacts occurring during the simulations are shown; instead only the largest ones as indicatedby the right-hand vertical axes. In 2-D, only the impactors that influence the equatorial annulus, i.e. occur with
their own generated heat anomaly of the equator, areshown. Clearly, the number of impacts influencing the2-D cases is far lower than the total number of impacts
in 3-D. On the other hand, each impact has a stronger
influence on the 2-D slice than it does in a 3-D domain,

As basaltic crust is produced at the surface, the extent towhich basalt is present at various depths in the mantle
gives a clear indication of crustal recycling and mantle
mixing. Figure 10 shows radial profiles of the azimuthally
averaged basalt fraction of the mantle at the end of eachsimulation. A composition of 1 is basalt and 0 is harzburgite. All simulations start homogeneous with basalt fraction 0.2.
The compositional profiles of the mantle are similar
between all the cases except above the CMB, around the660 km discontinuity and near the top. At the CMB and660 km discontinuity, the amount of basalt is dependent
on the amount of subducted basalt and thus on the yieldstress and the presence or absence of impacts: impacts
increase the amount of subducted basalt, particularly
above the CMB. The lower the yield stress, the morebasalt accumulates above the 660 km discontinuity. The
amount of basalt above the CMB is less dependent onyield stress except for the highest yield stress in 3-D, forwhich it is lower. 2-D and 3-D cases are somewhat similar but 2-D models display much more basalt accumulation above the CMB. The no-impacts case is also shownfor 85 MPa yield stress. In both geometries, compared tothe equivalent impacts case this displays much less basalt
accumulation above the CMB and somewhat less above
the 660 km discontinuity. It also displays a thicker crust

Fig. 8 Surface mobility versus time for 2-D and 3-D cases with and without impacts for three different yield stress values. In 2-D, the plots show thecases with a percentage of surface mobility > 1 equal to 93.28%, 57.14% and 48.29%, from the left to the right panel

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Page 10 of 19

Fig. 9 Surface mobility and impactor size versus time for cases with yield stress 125 MPa in 2-D and 3-D. Only the larger impacts are shown aslabelled on the right-hand vertical axis. In 2-D only impacts generating a thermal anomaly that intersects the simulation plane are plotted. In 2-D,
the case with impacts is the same as that in the right panel of Fig. 8

Fig. 10 Azimuthally averaged basalt fraction versus radius at the end of the simulated time for various yield stresses in 2-D and 3-D. A basalt
fraction of 1 is pure basalt, while a basalt fraction of 0 is harzburgite

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

and less depletion (harzburgite) below the crust. Impacts
generate melt in this region, thereby depleting it.
The top panels of Fig. 11 show the time-evolution of
the mean basaltic composition around the 660 km discontinuity (between 640 and 690 km depth) over thesimulated 0.5 Gyr. The basalt fraction clearly increases
with time, even in the no-impacts cases. In no-impacts
cases, dripping/erosion is responsible for removing thebase of the crust. Sudden increases and decreases are due
to a slab flowing for a while along the discontinuity andthen sinking into the lower mantle, which can cause thebasaltic composition to drop suddenly. As the yield stressincreases and surface mobility decreases, such suddenchanges become rarer. The 45 MPa 3-D cases indicatethat impacts can increase crustal recycling and radial
mixing even in cases that are anyway mobile lid.
The lower panels of Fig. 11 show the time evolutionof the mean basalt fraction of the bottom 225 km of the
mantle. The simulations with impacts generally containmore basalt above the CMB than the simulations without impacts, particularly when the no-impacts cases arestagnant lid. This is particularly pronounced in the 2-D
impact cases with 85 and 125 MPa yield stress. This

Page 11 of 19

indicates that impact-induced subduction is more likely
to happen in 2-D and that more material is subductedduring such an event. Simulations without impacts
generally have less basalt in the deep mantle. Sudden
increases of the deep mantle basalt fraction are linked toa sudden decrease of the basaltic composition at the 660
km discontinuity, indicating a causal linkage.
To sum up, the amount of recycled basalt at various
depths in the mantle is dependent on the surface mobility
and on the geometry of the model.
3.5 Distribution of impactor material

Whether or not late veneer material was mixed throughout the mantle has implications for geochemical observations. Thus, we track impactor material, starting with thesimple assumption that impactor material initially exists
at the impactor site. In reality, impactor material might
be distributed in a more complex manner that is dependent on the impact angle (e.g. Golabek et al. 2018) but thatis beyond the scope of the present study. We also track
the impactor’s metallic core, which has a much higher
density.

Fig. 11 Time-evolution of the azimuthally averaged basalt fraction for various yield stress values and with or without impacts for both 2-D and 3-D
cases. Top panels show the basalt fraction near the 660 km discontinuity (between 640 and 690 km depth), while the lower panels show the basalt
fraction of the bottom 245 km of the mantle. On the vertical-axis, 1.0 is basalt, while 0.0 is harzburgite. The spread of values for 5 different impact
histories in 2-D is shownBorgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

The radial distribution of impactor material at the endof the simulated period (Fig. 12) shows that most impactor material remains in the upper mantle; the concentration in the lower mantle is almost an order of magnitude
lower. There is a peak at the CMB. Impactor core material is mostly above the CMB, although some, particularly from small impactors, is stuck in the lithosphere.
Inside the mantle there are only trace amounts of metallic impactor core material. There are some differences
between 2-D and 3-D models. In 2-D simulations more
impactor material reaches the deep lower mantle thanin 3-D simulations, particularly for higher yield stresses.
This is consistent with earlier observations that in 2-D,
impact-induced subduction events are stronger and easily
reach the lower mantle, whereas in 3-D they are weaker
and tend to get trapped by the 660 km discontinuity.
The spatial distribution of impactor material is shownin Fig. 13, which helps to explain some observations
from Fig. 12. In neither geometry could the material be
described as well mixed—the distributions are still very

Page 12 of 19

heterogeneous. In most cases, there is an accumulationof impactor material above the 660 km discontinuity
just as there is with subducted crust: the two often travel
together.
3.6 Time series of surface heat flow

Impacts can also greatly influence the thermal evolution.
Here we analyse their influence on the surface heat flow,
including both conductive heat flow and magmatic heat
flow. Magmatic heat flow, which is the sum of heat loss
by cooling of erupted magma (from its initial temperature to the surface temperature) and latent heat release
due to its solidification (Armann and Tackley 2012; Nakagawa and Tackley 2012) is a potentially important heatloss mechanism during the Hadean even when impacts
are not present (Moore and Webb 2013; Nakagawa andTackley 2012; Lourenço et al. 2018). Magmatic heat flux
is directly related to crustal production rate, so we do notseparately plot that here.

Fig. 12 Azimuthally averaged impactor fraction vs. radius at the end of simulations with impacts for 2-D and 3-D cases with three different yieldstress values (45, 85 and 125 MPa). The top row shows silicate impactor fraction, while the bottom row shows impactor metal core fraction. The
spread of values for 5 different impact histories in 2-D is shownBorgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Page 13 of 19

Fig. 13 Slices of the distribution of impactor material at the end of simulations with impacts in 2-D and 3-D and for three different yield stressvalues. In 3-D, two perpendicular pole-to-pole slices are plotted. In 2-D, the same cases as in Fig. 8 are shown

The conductive heat flux is plotted in the top panelsof Fig. 14. For no-impacts 85 MPa yield stress cases, theheat flux is similarly low (10-20 TW) in both 2-D and
3-D. Impacts greatly increase the conductive heat flux
for the cases that would normally be stagnant lid (85and 125 MPa yield stress). As with other diagnostics
considered earlier, the effect is higher and more timedependent in 2-D cases. The 45 MPa yield stress caseswith impacts display a conductive heat flux that is low
before for the first 100 Myr then increases to around 60
TW once frequent impacts begin, in both geometries.
Magmatic heat flux (lower panel of Fig. 14) shows
similar trends. For the no-impacts 85 MPa cases it isinitially of similar magnitude to the conductive heat
flux, decreasing over time. Impacts bring all cases into
the same range as for mobile-lid cases, with much morevariability for the 2-D cases, where there are some very
high peaks.
Some additional heat is added by impacts, so higher
surface heat flux does not necessarily mean faster mantle cooling. We have, however, also tracked the heat
added by impacts and it is only a small fraction of thetotal heat budget, so not significant.

4 Discussion
A robust finding from the above analyses is that theeffects of impacts is different in 2-D models than in 3-D
models. In 2-D models the effect of each impact is larger
because the impact is effectively infinite in the out-ofplane direction. Induced subduction is stronger, more
likely, and more able to penetrate the 660 km discontinuity into the lower mantle. Acting against this, the numberof major impacts that influences the 2-D equatorial plane
is much less than the total number hitting the planet. The
net result is that 2-D models display more time-dependence with larger variability in mobility, heat flux anderuption (crustal production) rate, and display more deep
recycling of basaltic crust and of impactor material.
In 3-D, impact-induced subduction tends to be weakerand less material penetrates the 660 km discontinuity toreach the lower mantle. This has consequences for theradial distribution of impactor material and recycledbasaltic crust. The lower mantle is, in 3-D simulations,somewhat lacking in these two types of material compared to the 2-D simulations. The main consequence isthat the impactor material is not as well mixed between
upper and lower mantles in 3-D, with a lot of impactorBorgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Page 14 of 19

Fig. 14 Time-evolution of conductive and magmatic surface heat fluxes for impact and no-impacts cases in 2-D and 3-D with various yield stresses.
The upper panels show surface conductive heat flux, while the lower panels show the magmatic heat flux. The spread of values for 5 different
impact histories in 2-D is shownmaterial staying in the upper mantle. Thus, the choice ofmodel geometry has an influence on the results.
Our 2-D models are broadly consistent with those presented by O’Neill et al. (2017) (the paper mentions a 3-D
model, but no 3-D results are presented). It is encouraging that similar results are obtained using two completely
different numerical codes with somewhat different model
setups. For example, O’Neill et al. (2017) included Peierl’s
creep, which we do not; while we include the metallic
cores of impactors and have a more detailed parameterisation of phase transitions. There are several additional
aspects to our current study: we analyse three-dimensional models, include metallic impactor cores and melting and crustal production, and study compositional
mixing, focussing on the distribution of recycled basaltic
crust and of impactor material.
In recent years there has been much debate over thetime history of impactor flux, with the previously oftenaccepted idea of a Late Heavy Bombardment (e.g. Morbidelli et al. 2012; Bottke and Norman 2017) losing
ground to that of a monotonically decreasing impactor flux (e.g. Michael et al. 2018; Morbidelli et al. 2018;
Brasser et al. 2020; Zhu et al. 2019). Therefore, it isimportant to consider what a different impactor history
might imply for these calculations. Firstly, we note thatthe peak impact rate of the “Lunar Sawtooth Bombardment” (Morbidelli et al. 2012) assumed here is alreadylower than that of the classical LHB concept (see Fig. 1
in Hopkins and Mojzsis 2015) and that the impact ratecoincides with that of the often-used exponential decay
model of Neukum and Ivanov (1994) after 4.1 Gyr before
present (Fig. 3 in Morbidelli et al. 2012). Thus, if the simple exponential decay model was assumed, there would
be more impacts during the first 100 Myr of the simulation and the same number otherwise. A recent probabilistic estimate of impact flux vs. time (Brasser et al. 2020)
suggests lower impact fluxes by a factor of several at earlier times. This would result in less frequent induced subduction in the simulations.

5 Conclusions
Here we have examined the influence of Hadean to Eoarchean impacts on subduction and plate tectonics andmantle mixing using three- and two-dimensional numerical models. The simulations show that impacts can havea huge effect on global tectonics, inducing subduction
and resulting in a higher surface mobility, greater crustal
production and greater recycling of crust. In both geometries, impacts can induce subduction in cases that would
otherwise have a stagnant lid, facilitating lid mobility

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

over a large range of yield stress values. In such models,
subduction is a direct consequence of impacts and after
the impact flux dies down, subduction stops and the system returns to being stagnant lid. Thus, impacts do not
influence the long-term tectonic mode: they cannot permanently start plate tectonics.
It is important to highlight the differences between thetwo- and three-dimensional numerical models. Firstly,
while the average surface mobility behaves similarly in
2-D and 3-D models, the implications for the deep mantle are quite different. In 2-D models, subduction ofbasaltic crust and impactor material into the lower mantle occurs over a much wider range of yield stress valuesin 2-D compared to 3-D. In 2-D, subduction is easier for
an impact to trigger, lasts longer and tends to reach thelower mantle, whereas in 3-D, impact-induced subduction is weaker, short-lived and tends to reach only theupper mantle.
Secondly, there is less recycled basaltic crust in thelower mantle in three-dimensional cases than in twodimensional cases with identical parameters. This islinked to the absence of continuous subduction in 3-D.
Thirdly, the radial distribution of impactor material
differs with geometry. In 2-D, the impactor material
spreads throughout the mantle, allowing the mantle to beenriched in highly siderophile elements (HSEs) and otherelements coming from the impactor material, whereas in
3-D a larger fraction of it stays in the upper mantle. Thisis again a consequence of the absence of continuous subduction in 3-D. Thus, the consequences of an impact are
more important in 2-D.
To sum up, the results show that Hadean and Eoarchean impacts delivering the late veneer are very
important for Hadean and Eoarchean tectonics, crustal production and mantle mixing on Earth. In somecases, the addition of impacts can temporarily shift thetectonic regime from a stagnant lid towards a plate-likeregime. Thus, in order to understand the first billion
years of Earth’s history, it is essential to consider impacts.
Such impacts were also likely to have been important on
Venus, as was already argued for Venus’ outgassing history (Gillmann et al. 2020). However, further studies areneeded for Venus to understand the role of impacts withhigher surface temperatures.

6 Methods
First, we detail the impact model, then the convection
model and finally the cases run.
6.1 Impact history model

As discussed earlier, several chronologies for impactor
flux versus time have been proposed. Here, for comparison with previous results, we assume the sawtooth-likePage 15 of 19

impact history proposed by Morbidelli et al. (2012) for
the Moon. This model predicts a sudden increase inimpact rate 4.1 Gyr before present. Before the start of the
LHB, Morbidelli et al. (2012) use two different curves tobracket their estimate of N20, which is the number of craters with a diameter larger than 20 km. The equation forthe lower bound is:
4.5−t 0.5
dN20
= 0.02 ∗ e−( 0.01 )dt(1)

and the equation for the upper bound is:
4.5−t 0.34
dN20
= 0.025 ∗ e−( 0.003 )dt(2)After the increase in impactor flux at 4.1 Gyr, the equation becomes:

dN20= 2.7 ∗ 10−16 ∗ e6.93∗t + 5.9 ∗ 10−7dt(3)where t is given in Gyr.
This lunar impact history was re-scaled for the Earthby Marchi et al. (2014). The key parameters are the total
mass of LHB impactors and the size distribution. Spatially, impacts are assumed to be equally distributed over
Earth’s surface. Additionally, a model for the effect of an
impact on Earth’s interior is necessary, which depends onthe velocity and density of the impactor and its impact
angle.
6.1.1 Total mass of impactors

The total mass accreted by the Earth after initial formation, or “Late Veneer”, has been estimated using geochemical constraints. Marchi et al. (2014) estimatedthis mass using HSEs and found that it lies between∼ (0.7 − 3.0) ∗ 1022 kg , representing up to 0.5% of Earth’s
mass. Here we assume that a value at the lower end of
this range: a total mass of 1 ∗ 1022 kg (only 0.17% ofEarth’s mass), was added by the later stage of bombardment. For simplicity we assume that the density of the
impactors is equal to the density of the planet.
6.1.2 Size distribution of impacts

We assume a simple power-law for the number of impactors N as a function of impactor radius r:dNr= a( )−α
dr
r0

(4)

with α = −3.5, a commonly cited value (de Pater andLissauer 2015), and the constant a calculated to give theneeded total mass. For this exponent, most of the mass isdelivered by the largest impactors, so the largest impactor size is an important parameter. We assume a largest

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

impactor diameter of 1000 km, roughly the size of Ceres,the largest asteroid, following O’Neill et al. (2017). The
smallest impactor diameter is 20 km, arbitrarily chosen
in order to avoid having to process a large number ofvery small impacts that have negligible effect. Given the
total mass, size distribution and time history, impacts are
generated randomly at each time step. Each simulation
has a different random impact history. The normalized
size distribution of 30 simulations is shown in Fig. 15.
Impacts larger than 100 km radius are rare; therefore, thedistribution above this size is rather stochastic.
6.1.3 Spatial distribution of impacts

It is assumed that each point on Earth’s surface has anequal probability of being hit by an impact. Thus, theprobability distribution is flat with respect to longitude
and sinusoidal with respect to colatitude. Impact locations are generated randomly over the full sphere evenfor 2-D simulations; the 2-D domain is assumed to be anannulus around the equator and thus only impacts thatare close to this plane will influence it.
6.2 Effects of an impact

As in previous studies on the influence of impacts on the
mantle, we consider only the thermal effect of an impact,
not its mechanical effect such as the formation of craters,
ejecta and the redistribution of mass. Our model is basedon that of Senshu et al. (2002) and later used by Monteux
et al. (2007) and others (e.g. Golabek et al. 2011; Gillmann et al. 2016). An impact creates a shock wave that
propagates spherically. Around the impactor position,
the pressure rises and is homogeneous, forming an isobaric core.

Fig. 15 Normalized histogram of the size of impactors. Each pointof the histogram corresponds to a bin spanning 10 km in size. Thishistogram contains impactors of 30 different simulations. The orange
dashed curved has a slope of 10−3.5

Page 16 of 19

The ratio, γic , between the radius of the isobaric core
( ric ) and the impactor size ( rimp ) is given by the following formula from Senshu et al. (2002):

γic =

1
ric
= 3 3 ≈ 1.44
rimp

(5)

The pressure increase caused by the impact decreases
away from the isobaric core as the square of the distanced from the centre of the isobaric core (Monteux et al.
2007). The temperature increase in the isobaric core (T0)is then given by Monteux et al. (2007):

T0 =4π γ ρGR2
9 f (m) Cp(6)

where G is the universal gravitational constant, ρ is thedensity, Cp is the specific heat capacity, R is the radius of
the impacted body, γ is the fraction of kinetic energy converted into thermal energy (around 30%) and the function f(m) describes the volume that is heated and is equal
to 2.7 in our case. This equation is valid if the impact
velocity vimp is equal to the escape velocity vescape. Whenvimp
the ratio vescapeis larger than one the above equation mustbe multiplied by the square of this ratio (Gillmann et al.
2016). Here we assume an impact velocity of 25 km/s
(Marchi et al. 2014), which is 2.235 times Earth’s escapevelocity. So, according to this model, the temperature
increase does not depend on the size of the impactor. The
temperature increase decreases as a function of the distance, d, from the centre of the isobaric core as given by
Senshu et al. (2002); Golabek et al. (2011):

T (r) = �T (ric 4.4)d(7)

For an Earth-sized target, typical post-impact temperatures exceed the liquidus, causing a magma pond
(Monteux et al. 2009). Rather than attempting to model
magma pond evolution, we instead assume that vigorous
convection and surface heat loss cools it cools rapidly tothe rheological transition, thus effectively placing a limit
the post-impact temperature.
The impactor’s metallic core is also tracked. Even if
the impactor is initially undifferentiated, differentiation would likely occur in a magma pond caused by the
impact (Monteux et al. 2009). Thus, the impactor core,
which is assumed to have a radius of half the impactor’sradius, is placed in the lower part of the isobaric core.
Metal has a zero-pressure density of 7000 kg/m3 , a surface bulk modulus of 210 GPa and a bulk modulus pressure derivative of 3.9.

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

6.3 Convection model

Simulations are run for a period of 500 Myr, covering 4.2
Gyr to 3.7 Gyr before present. The exact start time does notmake a difference; 4.2 Gyr was chosen to allow convection
to start before the peak in impact rate occurs. The mantle
temperature is initialised to an adiabat with 1900 K potential temperature plus 50 km thick thermal boundary layers
at top and bottom and random temperature perturbationsof amplitude 20 K. This initial potential temperature isbased on petrological evidence that the mantle was 300 K
hotter then (Herzberg et al. 2010). The initial CMB temperature is 4200 K and the surface temperature is 300 K.
Additionally, plastic yielding is included in order to failure of the lithosphere, which can result in plate-like behaviour, episodic overturn or stagnant lid tectonics dependingon the assumed value of the yield stress (Moresi and Solomatov 1998; Tackley 2000). The yield stress is thus one ofmain parameters that was changed in the different simulations. Reference viscosity η0 was changed in some.

is the harmonic sum of diffusion creep and dislocation
creep viscosities, as usual.
6.3.2 Yield stress

Plastic yielding is assumed, which can lead to plate tectonics-like behaviour (e.g. Moresi and Solomatov (1998);
Tackley (2000)). Based on rock mechanics experiments(Kirby 1980; Kohlstedt et al. 1995) and as used in earlier
studies (e.g. Schierjott et al. (2020)) the yield stress σy has
a brittle component σB and a ductile component σD:σy = min(σB , σD )The deformation mechanisms diffusion creep and dislocation creep are included, with plastic failure occurring at
higher stresses. An Arrhenius law is used to described thebasic temperature-, pressure- and stress-dependence ofviscosity:

E + pV (p)E
σ0 n−1
) exp
−
σ
RT
RT0η

(8)
where E is the activation energy, p is the pressure, T is
the absolute temperature, n is the power-law index (1 fordiffusion creep and 3.5 for dislocation creep), σ0 is the
crossover stress between diffusion creep and dislocation creep, and V is the activation volume, which can bepressure-dependent as in Tackley et al. (2013). η0 is thereference viscosity, which is the viscosity at the reference
temperature T0η of 1600 K and zero pressure. E and V may
be different for each phase. For diffusion creep we use therheological parameters in Tackley et al. (2013), in which
all upper mantle phases have the same activation parameters based on Karato and Wu (1993), all lower mantle
phases except post-perovskite have the same activation
parameters based on Ammann et al. (2009), there is no
intrinsic viscosity jump between upper and lower mantles, but a 0.1 viscosity jump going from bridgmanite topost-perovskite, which has activation parameters basedon Ammann et al. (2010). Dislocation creep parameters
are for the lower mantle are the same as diffusion creep
parameters with a crossover stress of 3 MPa, while for theupper mantle E = 532 kJ/mol, V = 14 × 10−6 m3 and thecrossover stress is 0.3 MPa. The overall creep viscosity(9)

σB is proportional to pressure P as,
σB = cf P(10)

where cf is the friction coefficient in Byerlee’s law taken
to be 0.5. σD is given by:

σD = σD0 + σD′ P

6.3.1 Viscosity law

η(T , p, σ ) = η0 (Page 17 of 19(11)

where σD0 is the surface yield stress (at p=0) and σD′ isthe rate of increase of σD with pressure. σD0 is the main
parameter varied, while σD′ is assumed to be 0.001 to
avoid yielding in the deep mantle.
The effective viscosity is the minimum of the creep viscosity and the yielding viscosity:σyηeff = min η,(12)2ε̇
where ε̇ is the second invariant of the strain rate tensor.
The final viscosity value is truncated between 1018 Pa s
and 1026 Pa s.
6.4 Geometries

3-D spherical models use the yin-yang grid with 64 radial
cells and 192 × 64 × 2 azimuthal cells. As a resolution
test, one case was repeated with 288 × 96 × 2 azimuthal
cells and diagnostics such as the time evolution of mobility and radial distribution of basalt were found to bealmost the same. The two-dimensional simulations use afull spherical annulus with the same radial resolution of
64 cells and a horizontal resolution of 512 cells. Radialgrid refinement gives higher radial resolution near surface, the 660 km discontinuity the CMB (Schierjott et al.
2020). This annulus is a 2-D slice from a spherical shelldomain taken at the equator and neglecting all terms that
depend on latitude (Hernlund and Tackley 2008). Some
impacts can happen outside the plane of simulation but
still have an effect on the simulation itself. Higher resolution simulations with 1024 × 128 cells display similar
results, indicating that the lower resolution is sufficient.

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Table 1 Main parameters of the simulations

Page 18 of 19

of StagYY have benefitted from funding from the ETH/CSCS Platform for
Advanced Scientific Computing (PASC) program.

Start time of the simulations

4.2 Gyr

Start time of the impacts

4.1 Gyr

Duration of the simulations

0.5 Gyr

Reference viscosity, η0

1.0e20 Pa s

Yield Stress

25–195 MPa

Dimensions

2-D & 3-D

Initial Temperature1900 K

Surface Temperature

300 K

CMB Temperature

4200 K6.5 Cases

We ran more than 350 two-dimensional simulations and around 60 three-dimensional simulations,changing mainly the yield stress. Table 1 gives a shortoverview of the range of each parameter used for thedifferent simulations.
The simulations ran on the ETH Euler cluster on 32
cores for two-dimensional simulations and on 64 coresfor the three-dimensional simulations.
Abbreviations
CMB: Core-mantle boundary; LHB: Late heavy bombardment; YS: Yield stress.
Supplementary Information
The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s40645-​022-​00497-0.
Additional file 1. Animation 1: Animation showing the long-term evolution of three-dimensional simulations withand without impacts from
Fig. 5. The basalt fraction field is plotted. Time is relative to the onset ofheavybombardment.
Additional file 2. Animation 2: Animation showing 400 Myr evolution ofa three-dimensional simulation withimpacts and a yield stress of 85 MPa,showing temperature, strain rate and basalt at the surface and a crosssectionshowing the distribution of impactor material in the mantle. Timeis relative to the onset of heavy bombardment.
Acknowledgements
We thank Antoine Rozel for being a co-examiner on XB’s Bachelor’s thesis andthe ETH SEC and CSCS/PASC for their financial support.
Author contributions
PJT proposed the topic, conceived and designed the study. XB carried out theexperimental study, analysed the data and wrote the first manuscript draft. PJT
helped in the interpretation and in writing the final manuscript. Both authors
have read and approved the final manuscript.
Authors’ information
Xavier Borgeat is studying for a Masters Degree in Earth Sciences at ETH
Zurich and performed this study as his Bachelors thesis.
Paul J. Tackley has been a professor in the Institute of Geophysics at ETH
Zurich since 2005 and has studied mantle convection using numerical simulations since the early 1990s.
Funding
Open access funding provided by Swiss Federal Institute of Technology
Zurich. Simulations were run on shares of the Euler cluster bought by funding from the ETH Scientific Equipment Commission. Recent developments

Availability of data and materials
The datasets used and/or analysed during the current study are available fromthe corresponding author on reasonable request.

Declarations
Competing interests
The authors declare that they have no competing interest.
Author details
1
Department of Earth Sciences, ETH Zurich, Zurich, Switzerland. 2 Instituteof Geophysics, Department of Earth Sciences, ETH Zurich, Zurich, Switzerland.
Received: 31 October 2021 Accepted: 30 June 2022

References
Ammann MW, Brodholt JP, Dobson DP (2009) DFT study of migration enthalpies in MgSiO3 perovskite. Phys. Chem. Min 36:151–158. https://​doi.​org/​10.​
1007/​s00269-​008-​0265-z
Ammann MW, Brodholt JP, Wookey J, Dobson DP (2010) First-principles constraints on diffusion in lower-mantle minerals and a weak D’’ layer. Nature
465:462–465
Armann M, Tackley PJ (2012) Simulating the thermo-chemical magmatic andtectonic evolution of Venus’ mantle and lithosphere: two-dimensionalmodels. J. Geophys. Res. 117:12003–1010292012004231. https://​doi.​org/​
10.​1029/​2012J​E0042​31
Bercovici D, Ricard Y (2014) Plate tectonics, damage and inheritance. Nature508(7497):513–6. https://​doi.​org/​10.​1038/​natur​e13072
Bottke WF, Norman MD (2017) The late heavy bombardment. Ann Rev
Earth Planet Sci 45(1):619–647. https://​doi.​org/​10.​1146/​annur​
ev-​earth-​063016-​020131
Brasser R, Werner SC, Mojzsis SJ (2020) Impact bombardment chronology of
the terrestrial planets from 4.5 Ga to 3.5 Ga. Icarus. https://​doi.​org/​10.​1016/j.​icarus.​2019.​113514
Cawood PA, Hawkesworth CJ, Pisarevsky SA, Dhuime B, Capitanio FA, Nebel O(2018) Geological archive of the onset of plate tectonics. Philos Trans A
Math Phys Eng Sci. https://​doi.​org/​10.​1098/​rsta.​2017.​0405
Condie KC (2018) A planet in transition: the onset of plate tectonics on Earthbetween 3 and 2 Ga? Geosci Front 9(1):51–60. https://​doi.​org/​10.​1016/j.​
gsf.​2016.​09.​001
Davies GF (2008) Episodic layering of the early mantle by the “basalt barrier’’mechanism. Earth Planet Sci Lett 275:382–392
de Pater I, Lissauer JJ (2015) Planetary sciences, 2nd edn. Cambridge University
Press, Cambridge. https://​doi.​org/​10.​1017/​CBO97​81316​165270
Fischer R, Gerya T (2016) Early Earth plume-lid tectonics: a high-resolution 3D
numerical modelling approach. J Geodyn 100:198–214. https://​doi.​org/​
10.​1016/j.​jog.​2016.​03.​004
Gerya TV, Stern RJ, Baes M, Sobolev SV, Whattam SA (2015) Plate tectonics
on the Earth triggered by plume-induced subduction initiation. Nature
527(7577):221–225
Ghods A, Arkani-Hamed J (2007) Impact-induced convection as the main
mechanism for formation of lunar mare basalts. J. Geophys. Res.
112(E3):03005. https://​doi.​org/​10.​1029/​2006j​e0027​09
Gillmann C, Golabek GJ, Raymond SN, Schönbächler M, Tackley PJ, Dehant V,
Debaille V (2020) Dry late accretion inferred from Venus’s coupled atmosphere and internal evolution. Nat Geosci 13(4):265–269. https://​doi.​org/​
10.​1038/​s41561-​020-​0561-x
Gillmann C, Golabek GJ, Tackley PJ (2016) Effect of a single large impact on the
coupled atmosphere-interior evolution of Venus. Icarus 268:295–312
Golabek GJ, Emsenhuber A, Jutzi M, Asphaug EI, Gerya TV (2018) Coupling
SPH and thermochemical models of planets: methodology and example of a Mars-sized body. Icarus 301:235–246. https://​doi.​org/​10.​1016/j.​
icarus.​2017.​10.​003

Borgeat and Tackley ﻿Progress in Earth and Planetary Science(2022) 9:38

Golabek GJ, Keller T, Gerya TV, Zhu G, Tackley PJ, Connolly JA (2011) Origin of
the Martian dichotomy and Tharsis from a giant impact causing massive magmatism. Icarus 215(1):346–357
Hernlund JW, Tackley PJ (2008) Modeling mantle convection in the sphericalannulus. Phys Earth Planet Inter 171(1–4):48–54
Herzberg C, Condie K, Korenaga J (2010) Thermal history of the Earth and its
petrological expression. Earth Planet Sci Lett 292(1–2):79–88. https://​doi.​org/​10.​1016/j.​epsl.​2010.​01.​022
Hopkins M, Harrison TM, Manning CE (2008) Low heat flow inferred from>4 Gyr zircons suggests Hadean plate boundary interactions. Nature
456(7221):493–6. https://​doi.​org/​10.​1038/​natur​e07465
Hopkins M, Mojzsis SJ (2015) A protracted timeline for lunar bombardmentfrom mineral chemistry, Ti thermometry and U-Pb geochronology ofApollo 14 melt breccia zircons. Contrib Miner Petrol 169(3):30
Karato S, Wu P (1993) Rheology of the upper mantle—a synthesis. Science
260(5109):771–778
Kirby SH (1980) Tectonic stresses in the lithosphere: constraints provided by the experimental deformation of rocks. J. Geophys. Res.
85(B11):6353–63
Kirkland CL, Hartnady MIH, Barham M, Olierook HKH, Steenfelt A, Hollis JA(2021) Widespread reworking of Hadean-to-Eoarchean continents during Earth’s thermal peak. Nat Commun 12(1):331. https://​doi.​org/​10.​1038/​s41467-​020-​20514-4
Kohlstedt DL, Evans B, Mackwell SJ (1995) Strength of the lithosphere—
constraints imposed by laboratory experiments. J. Geophys. Res.
100(B9):17587–17602
Korenaga J (2013) Initiation and evolution of plate tectonics on Earth: theories
and observations. Ann Rev Earth Planet Sci 41:117–151
Lourenço DL, Rozel AB, Ballmer MD, Tackley PJ (2020) Plutonic-squishy lid: a
new global tectonic regime generated by intrusive magmatism on Earthlike planets. Geochem Geophys Geosyst 21(4):2019–008756
Lourenço DL, Rozel AB, Gerya T, Tackley PJ (2018) Efficient cooling of rocky
planets by intrusive magmatism. Nat Geosci 11(5):322–327. https://​doi.​
org/​10.​1038/​s41561-​018-​0094-8
Marchi S, Bottke W, Elkins-Tanton L, Bierhaus M, Wuennemann K, Morbidelli A,
Kring D (2014) Widespread mixing and burial of Earth’s Hadean crust byasteroid impacts. Nature 511(7511):578–582
Maruyama S, Santosh M, Azuma S (2018) Initiation of plate tectonics in the
Hadean: Eclogitization triggered by the ABEL Bombardment. Geosci Front
9(4):1033–1048
Michael G, Basilevsky A, Neukum G (2018) On the history of the early meteoritic bombardment of the Moon: was there a terminal lunar cataclysm?
Icarus 302:80–103. https://​doi.​org/​10.​1016/j.​icarus.​2017.​10.​046
Monteux J, Ricard Y, Coltice N, Dubuffet F, Ulvrova M (2009) A model ofmetal-silicate separation on growing planets. Earth Planet Sci Lett
287(3–4):353–362
Monteux J, Arkani-Hamed J (2014) Consequences of giant impacts in earlyMars: core merging and Martian dynamo evolution. J Geophys Res Planets 119(3):480–505. https://​doi.​org/​10.​1002/​2013j​e0045​87
Monteux J, Coltice N, Dubuffet F, Ricard Y (2007) Thermo-mechanical adjustment after impacts during planetary growth. Geophys Res Lett 34(24).
https://​doi.​org/​10.​1029/​2007G​L0316​35
Moore WB, Webb AAG (2013) Heat-pipe Earth. Nature 501(7468):501–505.
https://​doi.​org/​10.​1038/​natur​e12473
Morbidelli A, Nesvorny D, Laurenz V, Marchi S, Rubie DC, Elkins-Tanton L,
Wieczorek M, Jacobson S (2018) The timeline of the lunar bombardment:
revisited. Icarus 305:262–276. https://​doi.​org/​10.​1016/j.​icarus.​2017.​12.​046
Morbidelli A, Marchi S, Bottke WF, Kring DA (2012) A sawtooth-like timeline
for the first billion years of lunar bombardment. Earth Planet Sci Lett
355:144–151
Moresi L, Solomatov V (1998) Mantle convection with a brittle lithosphere—
thoughts on the global tectonic styles of the Earth and Venus. Geophys.
J. Int. 133(3):669–682
Moyen J-F, Van Hunen J (2012) Short-term episodicity of Archaean plate
tectonics. Geology 40(5):451–454
Nakagawa T, Tackley PJ (2012) Influence of magmatism on mantle cooling,surface heat flow and Urey ratio. Earth Planet. Sci. Lett. 329–330:1–10.
https://​doi.​org/​10.​1016/j.​epsl.​2012.​02.​011
Neukum G, Ivanov BA (1994) Crater size distributions and impact probabilities on earth from lunar, terrestrial-planet, and asteroid cratering data.

Page 19 of 19

In: Gehrels T (ed) Hazards Due to Comets and Asteroids. University of
Arizona Press, Tucson, AZ, USA, pp 359–416
O’Neill C, Debaille V (2014) The evolution of Hadean–Eoarchaean geodynamics. Earth Planet Sci Lett 406:49–58
O’Neill C, Marchi S, Zhang S, Bottke W (2017) Impact-driven subduction on the
Hadean earth. Nat Geosci 10(10):793–797
Palin RM, Santosh M (2020) Plate tectonics: What, where, why, and when?
Gondwana Res. https://​doi.​org/​10.​1016/j.​gr.​2020.​11.​001
Piccolo A, Palin RM, Kaus BJP, White RW (2019) Generation of Earth’s earlycontinents from a relatively cool Archean mantle. Geochem Geophys
Geosyst 20(4):1679–1697. https://​doi.​org/​10.​1029/​2018g​c0080​79
Piccolo A, Kaus BJP, White RW, Palin RM, Reuber GS (2020) Plume—Lid interactions during the Archean and implications for the generation of earlycontinental terranes. Gondwana Res 88:150–168. https://​doi.​org/​10.​1016/j.​gr.​2020.​06.​024
Reese CC, Orth CP, Solomatov VS (2010) Impact origin for the martian crustal
dichotomy: half-emptied or half filled. J Geophys Res. https://​doi.​org/​10.​1029/​2009j​e0035​06
Reese CC, Solomatov VS, Baumgardner JR (2002) Survival of impact-induced
thermal anomalies in the Martian mantle. J Geophys Res Planets
107(E10):5082
Roberts JH, Lillis RJ, Manga M (2009) Giant impacts on early Mars and thecessation of the Martian dynamo. J Geophys Res. https://​doi.​org/​10.​1029/​
2008j​e0032​87
Rolf T, Zhu MH, Wünnemann K, Werner SC (2017) The role of impact bombardment history in lunar evolution. Icarus 286:138–152. https://​doi.​org/​10.​1016/j.​icarus.​2016.​10.​007
Schierjott J, Rozel A, Tackley P (2020) On the self-regulating effect of grain size
evolution in mantle convection models: application to thermochemical
piles. Solid Earth 11(3):959–982
Senshu H, Kuramoto K, Matsui T (2002) Thermal evolution of a growing Mars. Jf
Geophys Res Planets 107(E12):1–1
Sizova E, Gerya T, Brown M, Perchuk LL (2010) Subduction styles in the Precambrian: insight from numerical experiments. Lithos 116(3–4):209–229
Tackley PJ (2000) Self-consistent generation of tectonic plates in timedependent, three-dimensional mantle convection simulations. Geochem
Geophys Geosyst 1(8). https://​doi.​org/​10.​1029/​2000G​C0000​36
Tackley PJ (2008) Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang
grid. Phys Earth Planet Inter 171(1–4):7–18
Tackley PJ, Ammann M, Brodholt JP, Dobson DP, Valencia D (2013) Mantle
dynamics in super-Earths: post-perovskite rheology and self-regulation ofviscosity. Icarus 225(1):50–61
Tang CA, Webb AAG, Moore WB, Wang YY, Ma TH, Chen TT (2020) Breaking
Earth’s shell into a global plate network. Nat Commun 11(1):3621. https://​doi.​org/​10.​1038/​s41467-​020-​17480-2
Turcotte DL (1993) An episodic hypothesis for Venusian tectonics. J. Geophys.
Res. 98(E9):17061–17068
van Hunen J, van den Berg AP (2008) Plate tectonics on the early Earth: limitations imposed by strength and buoyancy of subducted lithosphere.
Lithos 103(1–2):217–235. https://​doi.​org/​10.​1016/j.​lithos.​2007.​09.​016
Yan J, Ballmer MD, Tackley PJ (2020) The evolution and distribution of recycled
oceanic crust in the Earth’s mantle: insight from geodynamic models.
Earth Planet Sci Lett. https://​doi.​org/​10.​1016/j.​epsl.​2020.​116171
Zhu MH, Artemieva N, Morbidelli A, Yin QZ, Becker H, Wunnemann K(2019) Reconstructing the late-accretion history of the Moon. Nature
571(7764):226–229. https://​doi.​org/​10.​1038/​s41586-​019-​1359-0

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.