Hadean/Eoarchean tectonics and mantle mixing induced by impacts: a three-dimensional study

The timing of the onset of plate tectonics on Earth remains a topic of strong debate, as does the tectonic mode that preceded 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 the simulations 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 early tectonics and mantle mixing of Earth and other large terrestrial planets, and that three-dimensional simulations are important 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.

time until some mechanism managed to break it and initiate 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 Hadean and Eoarchean crust inferred from zircons (Kirkland et al. 2021), and this intrusion would have weakened the lithosphere and allowed it to deform, accommodating the lateral motion of sections of it. Such behaviour has been demonstrated in both regional models (Fischer and Gerya 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(Lourenço et al. , 2018, occurring for intrusion fractions of ≥ 80% . Thus, while there is general agreement that any subduction that started under early Earth conditions would not result in long-lived subduction zones due to the weakness of the lithosphere (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 is also consistent with evidence for underthrusting from >4 Gyr zircons (Hopkins et al. 2008). Large plumes might also 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, through impacts. 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 been caused 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 of impacts 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. 2002Reese et al. , 2010Roberts 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-like planet, causing an episode of subduction that rolls back from the impact location. Multiple impacts continuing for some period of time could thus potentially mobilise the lithosphere for a corresponding time period. This was 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 was demonstrated 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 order to 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.

Methods summary
Full details of the physical and numerical model are given in the later Methods section; here is a brief overview. We couple 3-D or 2-D simulations of mantle convection with a 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 the impact 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 a total mass addition of 1.0 * 10 20 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 is modelled, but other direct mechanical effects are not included.
Convection and lithospheric dynamics is simulated using StagYY (Tackley 2008) with assumed physical parameters 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. Compositional variations 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

Results
In this section, we first present the effects of a single impact, and then analyse the long-term evolution of the simulations, contrasting cases with impacts and without impacts, and in 3-D versus 2-D. Both visual aspects and quantitative measures (such as surface mobility and heat flux) are presented. Finally, we focus on compositional aspects, in particular the radial distribution of basaltic material and the radial distribution of impactor material. Figure 1 shows the effect of a 400 km radius impact in a two-dimensional simulation. The heat addition generates some melt, some of which produces more crust and depleted harzburgitic material beneath. The buoyancy of the hot, depleted, partially molten material produces high 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 of events is in line with previous simulations (O'Neill et al. 2017;Gillmann et al. 2016). Additionally, we track the metallic 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 not occur at the edge of the impactor material, but instead within the impactor material. Some of the subducted material flows along the 660km discontinuity, while some reaches the lower mantle. In contrast to the two-dimensional simulations, no continuous subduction events lasting for several million years were observed in 3-D. This is the first main difference observed between the two geometries.

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 of  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, two movies 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 shows a 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 at some times, linear low-viscosity features corresponding to zones of extension or subduction, i.e. tectonic plate boundaries.
In the three-dimensional cross sections (Fig. 6), it is clear that the no-impacts case is stagnant lid, while in the impacts 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 the 660 km discontinuity.

Effect of the yield stress
The key parameter controlling the tectonic mode in visco-plastic mantle convection simulations is the yield stress, 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 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 to zero means that the surface is immobile, indicating the stagnant-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) and the yield stress for cases with both geometries and with and 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 stress of around 50-80 MPa. (There is some complexity in the 2-D case due to random effects.) In strong contrast, cases with impacts maintain a mobile surface over a much wider 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 yield stress values that would normally cause a stagnant lid regime.
Comparing 2-D and 3-D, the trends are very similar both qualitatively and quantitatively, but there is more stochastic variation (randomness) in 2-D. This is something that is also apparent from other diagnostics discussed later. Another difference is that the mobile-lid to stagnant-lid transition occurs at higher yield stresses in 2-D cases. Figure 8 shows the surface mobility vs. time for the three different yield stresses in the both dimensions, with and without impacts. The sawtooth bombardment ramps up at 0.1 Gyr. For the lowest yield stress (45 MPa), all cases 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, are mostly 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 also quite time-dependent but the peaks are lower than those of the 2-D cases. Particularly in the 85 MPa 3-D case, the mobility 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 tectonic mode 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. Page 9 of 19 Borgeat and Tackley Progress in Earth and Planetary Science (2022) 9:38 The relationship between periods of surface mobility and individual large impacts is shown in Fig. 9. Spikes in mobility 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 indicated by 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, are shown. Clearly, the number of impacts influencing the 2-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, causing mobility that lasts for several tens of millions of years. When there hasn't been an impact for some time, the system returns to a stagnant lid state. These two points: fewer impacts each with more influence, explain the difference in mobility time series between 2-D and 3-D models.

Radial distribution of basalt
As basaltic crust is produced at the surface, the extent to which 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 each simulation. 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 the 660 km discontinuity and near the top. At the CMB and 660 km discontinuity, the amount of basalt is dependent on the amount of subducted basalt and thus on the yield stress and the presence or absence of impacts: impacts increase the amount of subducted basalt, particularly above the CMB. The lower the yield stress, the more basalt accumulates above the 660 km discontinuity. The amount of basalt above the CMB is less dependent on yield stress except for the highest yield stress in 3-D, for which 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 shown for 85 MPa yield stress. In both geometries, compared to the 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   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 the simulated 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 the base of the crust. Sudden increases and decreases are due to a slab flowing for a while along the discontinuity and then sinking into the lower mantle, which can cause the basaltic composition to drop suddenly. As the yield stress increases and surface mobility decreases, such sudden changes become rarer. The 45 MPa 3-D cases indicate that 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 evolution of the mean basalt fraction of the bottom 225 km of the mantle. The simulations with impacts generally contain more basalt above the CMB than the simulations without impacts, particularly when the no-impacts cases are stagnant lid. This is particularly pronounced in the 2-D impact cases with 85 and 125 MPa yield stress. This indicates that impact-induced subduction is more likely to happen in 2-D and that more material is subducted during such an event. Simulations without impacts generally have less basalt in the deep mantle. Sudden increases of the deep mantle basalt fraction are linked to a 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.

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 the simple 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 that is beyond the scope of the present study. We also track the impactor's metallic core, which has a much higher density. The radial distribution of impactor material at the end of 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 than in 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 shown in 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 heterogeneous. In most cases, there is an accumulation of impactor material above the 660 km discontinuity just as there is with subducted crust: the two often travel together.

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 heat loss mechanism during the Hadean even when impacts are not present (Moore and Webb 2013; Nakagawa and Tackley 2012;Lourenço et al. 2018). Magmatic heat flux is directly related to crustal production rate, so we do not separately plot that here. The conductive heat flux is plotted in the top panels of Fig. 14. For no-impacts 85 MPa yield stress cases, the heat 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 (85 and 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 cases with 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 is initially 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 more variability 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 the total heat budget, so not significant.

Discussion
A robust finding from the above analyses is that the effects 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 number of 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 and eruption (crustal production) rate, and display more deep recycling of basaltic crust and of impactor material.
In 3-D, impact-induced subduction tends to be weaker and less material penetrates the 660 km discontinuity to reach the lower mantle. This has consequences for the radial distribution of impactor material and recycled basaltic 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 is that the impactor material is not as well mixed between upper and lower mantles in 3-D, with a lot of impactor  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 the time 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 is important to consider what a different impactor history might imply for these calculations. Firstly, we note that the peak impact rate of the "Lunar Sawtooth Bombardment" (Morbidelli et al. 2012) assumed here is already lower than that of the classical LHB concept (see Fig. 1 in Hopkins and Mojzsis 2015) and that the impact rate coincides 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.

Conclusions
Here we have examined the influence of Hadean to Eoarchean impacts on subduction and plate tectonics and mantle mixing using three-and two-dimensional numerical models. The simulations show that impacts can have a 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  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 the two-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 of basaltic crust and impactor material into the lower mantle occurs over a much wider range of yield stress values in 2-D compared to 3-D. In 2-D, subduction is easier for an impact to trigger, lasts longer and tends to reach the lower mantle, whereas in 3-D, impact-induced subduction is weaker, short-lived and tends to reach only the upper mantle.
Secondly, there is less recycled basaltic crust in the lower mantle in three-dimensional cases than in twodimensional cases with identical parameters. This is linked 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 be enriched in highly siderophile elements (HSEs) and other elements coming from the impactor material, whereas in 3-D a larger fraction of it stays in the upper mantle. This is 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 some cases, the addition of impacts can temporarily shift the tectonic regime from a stagnant lid towards a plate-like regime. 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 are needed for Venus to understand the role of impacts with higher surface temperatures.

Methods
First, we detail the impact model, then the convection model and finally the cases run.

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-like impact history proposed by Morbidelli et al. (2012) for the Moon. This model predicts a sudden increase in impact rate 4.1 Gyr before present. Before the start of the LHB, Morbidelli et al. (2012) use two different curves to bracket their estimate of N 20 , which is the number of craters with a diameter larger than 20 km. The equation for the lower bound is: and the equation for the upper bound is: After the increase in impactor flux at 4.1 Gyr, the equation becomes: where t is given in Gyr.
This lunar impact history was re-scaled for the Earth by 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 on the velocity and density of the impactor and its impact angle.

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) estimated this mass using HSEs and found that it lies between ∼ (0.7 − 3.0) * 10 22 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 * 10 22 kg (only 0.17% of Earth'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.

Size distribution of impacts
We assume a simple power-law for the number of impactors N as a function of impactor radius r: with α = −3.5 , a commonly cited value (de Pater and Lissauer 2015), and the constant a calculated to give the needed total mass. For this exponent, most of the mass is delivered by the largest impactors, so the largest impactor size is an important parameter. We assume a largest (1) dN 20 dt = 0.02 * e −( 4.5−t 0.01 ) 0.5 (2) dN 20 dt = 0.025 * e −( 4.5−t 0.003 ) 0.34 (3) dN 20 dt = 2.7 * 10 −16 * e 6.93 * t + 5.9 * 10 −7 (4) dN dr = a( r r 0 ) −α 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 of very 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, the distribution above this size is rather stochastic.

Spatial distribution of impacts
It is assumed that each point on Earth's surface has an equal probability of being hit by an impact. Thus, the probability distribution is flat with respect to longitude and sinusoidal with respect to colatitude. Impact locations are generated randomly over the full sphere even for 2-D simulations; the 2-D domain is assumed to be an annulus around the equator and thus only impacts that are close to this plane will influence it.

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 based on 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.
The ratio, γ ic , between the radius of the isobaric core ( r ic ) and the impactor size ( r imp ) is given by the following formula from Senshu et al. (2002): The pressure increase caused by the impact decreases away from the isobaric core as the square of the distance d from the centre of the isobaric core (Monteux et al. 2007). The temperature increase in the isobaric core ( T 0 ) is then given by Monteux et al. (2007): where G is the universal gravitational constant, ρ is the density, C p 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 v imp is equal to the escape velocity v escape . When the ratio v imp v escape is larger than one the above equation must be 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 escape velocity. 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): 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 to the 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's radius, is placed in the lower part of the isobaric core. Metal has a zero-pressure density of 7000 kg/m 3 , a surface bulk modulus of 210 GPa and a bulk modulus pressure derivative of 3.9.

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 not make 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 perturbations of amplitude 20 K. This initial potential temperature is based 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 depending on the assumed value of the yield stress (Moresi and Solomatov 1998;Tackley 2000). The yield stress is thus one of main parameters that was changed in the different simulations. Reference viscosity η 0 was changed in some.

Viscosity law
The deformation mechanisms diffusion creep and dislocation creep are included, with plastic failure occurring at higher stresses. An Arrhenius law is used to described the basic temperature-, pressure-and stress-dependence of viscosity: where E is the activation energy, p is the pressure, T is the absolute temperature, n is the power-law index (1 for diffusion 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 be pressure-dependent as in Tackley et al. (2013). η 0 is the reference viscosity, which is the viscosity at the reference temperature T 0 η of 1600 K and zero pressure. E and V may be different for each phase. For diffusion creep we use the rheological 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 to post-perovskite, which has activation parameters based on 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 the upper mantle E = 532 kJ/mol, V = 14 × 10 −6 m 3 and the crossover stress is 0.3 MPa. The overall creep viscosity (8) η(T , p, σ ) = η 0 ( σ 0 σ ) n−1 exp E + pV (p) RT − E RT 0 η is the harmonic sum of diffusion creep and dislocation creep viscosities, as usual.

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 : σ B is proportional to pressure P as, where c f is the friction coefficient in Byerlee's law taken to be 0.5. σ D is given by: where σ D0 is the surface yield stress (at p=0) and σ ′ D is the 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: where ε is the second invariant of the strain rate tensor.
The final viscosity value is truncated between 10 18 Pa s and 10 26 Pa s.

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 be almost the same. The two-dimensional simulations use a full spherical annulus with the same radial resolution of 64 cells and a horizontal resolution of 512 cells. Radial grid 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 shell domain 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.

Cases
We ran more than 350 two-dimensional simulations and around 60 three-dimensional simulations, changing mainly the yield stress. Table 1 gives a short overview of the range of each parameter used for the different simulations. The simulations ran on the ETH Euler cluster on 32 cores for two-dimensional simulations and on 64 cores for the three-dimensional simulations.
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 of heavybombardment.
Additional file 2. Animation 2: Animation showing 400 Myr evolution of a 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. Time is relative to the onset of heavy bombardment.