Onset of Plate Motion in the Presence of Chemical Heterogeneities in the Mantle and the Effect of Mantle Temperature

Chemical heterogeneities of various origins have been observed in the Earth's mantle. Their assumed higher density acts on the style of mantle convection and therefore, as tectonic plates form the highly viscous upper boundary layer of mantle convection, the onset of surface motion is also affected by chemical heterogeneities. We perform 2D basally heated thermochemical mantle convection models which initially include layers of dense material at different mantle depths. In comparison to purely thermal convection models, the onset of first plate motion is delayed. This delay is even more pronounced, the deeper the location of dense material, the larger its volume and the higher its density. Additionally, we varied the starting temperature of our models. In an initially hot mantle, a strong temperature‐dependent viscosity contrast between the mantle and cold surface leads to a cold, highly viscous plate (stagnant lid). For an initially cold mantle, rising plumes first need to transport basal heat upwards to create a sufficient viscosity contrast. Consequently, the formation and subsequent breaking of the stagnant lid is delayed. Considering a hotter mantle after Earth's magma ocean phase, we find that plate motion can occur within approximately the first 0.5 Gyr of solid‐state convection if no chemical structures are present or for dense material situated at the surface. Deep chemical heterogeneities delay the onset by at least one order of magnitude. Furthermore, the early surface motion will have been more erratic than todays stable plate motion, probably with a few single subduction events.


Introduction
The initiation of plate motion has widely been of interest in geodynamical and geochemical studies, but is still under debate (for reviews see Korenaga, 2013;Cawood et al., 2018;Condie, 2018;Wyman, 2018;Hawkesworth et al., 2020;Palin & Santosh, 2020).The onset time, for example, is placed either late in Earth's evolution, that is, in the Neoproterozoic (e.g., Hamilton, 2011;Stern, 2020), or the Archean (e.g., Brown & Johnson, 2018;Condie, 2018) or even as early as in the Hadean (e.g., Borgeat & Tackley, 2022;Maruyama et al., 2018).Discrepancies as to when plate motion started are partly due to definition.Typically, studies favoring late onset of plate tectonics exclude the modern-style plate motion with deep subduction and slab pull as primary driving force in the hotter early Earth, because the lithosphere in a hotter mantle would have been thinner and therefore not (or less) negatively buoyant (cf.Stern, 2007).However, plate tectonics probably operated differently at first.Evidence of crustal growth and many geodynamical studies suggest that an initial episodic behavior is likely (e.g., Bedard, 2018;Debaille et al., 2013;Kreielkamp et al., 2022).Short-lived subduction episodes (e.g., Moore & Webb, 2013;O'Neill & Debaille, 2014;van Hunen & Moyen, 2012) or delamination in form of lithospheric drips as the dominant process, like for the regional plume-induced subduction (plume-lid tectonics) or global intrusive magmatism (plutonic squishy lid) are discussed (e.g., Gerya et al., 2015;Fischer & Gerya, 2016;Lourenço et al., 2020;Baes et al., 2021, and references therein).Further mechanisms that have been proposed are, for example, drip-like events due to grain damage, leading to often but locally confined subduction (Foley, 2018), or impact-induced subduction that only occurs temporary after Hadean or Eoarchean impacts (e.g., Borgeat & Tackley, 2022;O'Neill et al., 2017).A reason for transient subduction is seen in the weakness of the Archean Earth's lithosphere (Moyen & van Hunen, 2012), where slab break-offs and re-initiation happens every few million years fitting geochemical observation of alternating arc-and non-arc-style volcanism.
The strength of the lithosphere is strongly controlled by the mantle temperature (Stein et al., 2013), which is assumed to have been maybe as much as ∼1500-1600°C about 3 Gyr ago (compared to ∼1350°C today), and decreased since the Archean due to secular cooling and a reduction in radiogenic heat production (Brown, 2014;Herzberg et al., 2010;Labrosse & Jaupart, 2007;and references therein).The higher temperature is consistent with petrological estimations of the Archean Barberton Greenstone belt in the Kaapvaal craton (Abbott et al., 1994).The thinner lithosphere forming in a hotter Earth would take longer to become negatively buoyant, which led Stern (2007) to the conclusion that plate motion only occurred sporadically with a more magmatic style in the Archean, while the modern-style only emerged later in the Neoproterozoic after the Earth has cooled.Likewise, O'Neill et al. (2016) describe the change from stagnant lid convection to episodic and finally steady plate motion when cooling the mantle with time.This is similar to the observation of Stein et al. (2013), where a decrease of internal heating resulted in a weakening of the stagnant lid and yielded episodic or even mobile lid behavior.The weak plate is characterized by a low viscosity contrast (due to the reduced temperature difference) between the interior and the surface.In this way, only small stresses are necessary to break the plate, either by a plume impinging or by subduction due to gravitational instability.
The higher mantle temperature also leads to basalt formation by decompression melting resulting in a thick positively buoyant mafic protocrust (Sizova et al., 2015;Zegers & van Keken, 2001) which has been considered as an argument against early subduction (Vlaar et al., 1994).However, eclogitization of a part of the lower crust could have caused delamination (Zegers & van Keken, 2001), possibly triggering a large resurfacing event (van Thienen and van den Berg, 2004).As seen in zircons, the increased amount of melting will have enhanced the recycling of the Hadean crust (Kirkland et al., 2021).For example, zircons from the Jack Hills region (Australia) have been used to argue for subduction in the Hadean (Hopkins et al., 2008).Mixing times of short-lived 142 Nd-isotopes, however, suggest a stagnant lid mode for the most part with only limited tectonic episodes (Debaille et al., 2013).
All these various events obscure Earth's earliest geological record, explaining why the initiation of plate motion remains poorly understood.The accessible rock record of the unaltered early Earth is sparse.For example, the earliest formed continental crust, largely made up by the tonalite-trondhjemite-granodiorite (TTG) suite (Sizova et al., 2015), is preserved in only two pristine terranes, the Pilbara craton (Australia) and Kaapvaal craton (South Africa), both essentially undeformed until ∼3.2 Gyr ago (Zegers & van Keken, 2001).The earliest rock record is limited to gneiss complexes, such as the 3.55-3.85Gyr-old Itsaq gneiss complex in Greenland (Moore & Webb, 2013), while for even earlier times (>4 Gyr) only zircons can be used to constrain the dynamics of the Earth.The oldest zircons from the Jack Hills in Australia are nearly 4.4 Gyr old (Moore & Webb, 2013).
An alternative method of investigation is therefore given by numerical simulations that consider conditions relevant for the early Earth after the magma ocean period (e.g., Borgeat & Tackley, 2022;Debaille et al., 2013;Fischer & Gerya, 2016;Foley, 2018;van Hunen & Moyen, 2012).In this context, dense accumulations, that could be remnants of the magma ocean period (Ballmer, Lourenço, et al., 2017;Caracas et al., 2019;Labrosse et al., 2007), will also play an important role, as the material interacts with the convective flow.Additional sources for chemical heterogeneities could be dense material penetrating the mantle either by core-mantle interaction (Stein & Hansen, 2023) or by impacts (e.g., Borgeat & Tackley, 2022).On the one hand, the dense material moves with the flow, forming dense piles beneath mantle plumes and being largely absent beneath subducted slabs (e.g., Lassak et al., 2007;McNamara & Zhong, 2005;Stein & Hansen, 2023;Steinberger & Torsvik, 2012), on the other hand, dense material has a restoring effect on rising plumes because the higher density counteracts with the thermal uplift and reduces the buoyancy (e.g., Dannberg & Sobolev, 2015).Consequently this affects plate motion, the surface expression of mantle convection.Large density contrasts or extensive volumes of dense material can strongly delay stable surface motion leading to initial episodic behavior (Kreielkamp et al., 2022;Stein et al., 2020) or even prevent mobilization completely leading to stagnant lid convection (Trim et al., 2014).
Here, we focus on the initiation of the surface motion in a convecting mantle under the influence of dense material and for different initial mantle temperatures.The model's starting condition is the post-magma ocean period.The two end-member scenarios of magma ocean crystallization are fractional and equilibrium (batch) crystallization (e.g., Elkins-Tanton, 2012;Solomatov, 2015;Zhang et al., 2022).In fractional crystallization of the magma ocean, crystals settle and the residual melt becomes more enriched in incompatible elements (e.g., iron).Therefore, in the subsequent solid-state convection a prescribed dense basal layer is considered as a remnant of a basal magma ocean (Ballmer, Lourenço, et al., 2017;Labrosse et al., 2007).However, magma ocean models revealed that dense material can be kept in suspension under the influence of strong rotation (e.g., Maas & Hansen, 2019) or in the equilibrium crystallization (e.g., Solomatov, 2015;Zhang et al., 2022).An example for equilibrium crystallization is given by Caracas et al. (2019), who find that a density crossover between solids and liquids can lead to crystal accumulations at the crossover depth, yielding a surficial and a basal magma ocean.The density crossover depth is the depth, where the densities of the melt and bridgmanite crystals are equal and the system reaches neutral buoyancy.Due to the vigorous convection in the magma ocean the crystals remain entrained.Another possibility for less deep chemical heterogeneities are therefore strong bridgmanite-enriched ancient mantle structures (BEAMS) which (due to their high viscosity contrast of 30-50) prevent efficient mantle mixing and persist at depths of 1,000-2,200 km beyond 4.5 Gyr (Gülcher et al., 2020(Gülcher et al., , 2021)).Additionally, material introduced to the Earth by later impacts (Borgeat & Tackley, 2022;O'Neill et al., 2017) might be an explanation for less deep chemical heterogeneities.
We therefore analyze a number of model scenarios featuring variations in mantle temperature, thickness, density and depth of a layer of dense material.In particular, we determine the first collapse of the stagnant lid as onset time of surface motion.This value serves as a quantitative measure for the comparison between different simulations and gives estimates on the time, surface motion initiated on Earth.

Model
In this study, we apply an incompressible thermochemical mantle convection model that solves the energy equations (for temperature T and composition C), the continuity equation and the momentum equation: ∇p + ∇ ⋅ (ηϵ) + Ra(T BC)e z = 0.
(4) t is the thermal diffusion time, v the velocity, ρ the density, H the internal heating rate, p the dynamic pressure, η the effective viscosity (cf.Equation 9), ϵ the strain-rate tensor and e z the unity vector in vertical direction.
The thermal Rayleigh number, buoyancy ratio and the Lewis number comprise the following system and fluid parameters: the thermal expansivity α, the gravitational acceleration g, the density ρ 0 , and the viscosity η 0 (all values defined at the surface), and the superadiabatic temperature drop over the system ΔT, the system height d, the thermal diffusivity κ T and the chemical diffusivity κ c .Δρ is the density difference between dense accumulations and the ambient mantle.It can be computed by applying α = 1.2 × 10 5 K 1 , the reference density of ρ 0 = 3,300 kg/m 3 at the surface and a temperature drop of ΔT = 2500 K (Trim et al., 2014) from the buoyancy ratio as: Δρ is the maximum density difference that appears at the beginning of a simulation before the piles erode.Alternatively, we can compute the density difference Δρ = ρ ρ 0 at each point in time from the equation of state: With respect to the reference density ρ 0 , the normalized density contrast is δρ = Δρ respectively).Therefore a buoyancy ratio of B = 1 yields a density difference of Δρ = 99 kg/m 3 between dense accumulations and the ambient mantle, which means the density of the accumulations is δρ = 3% of the reference density.
As the main topic of this work is the initiation of surface motion, we use a rheological model to allow for surface motion according to the convective flow.The effective viscosity depends on temperature, depth and stress: with a temperature/depth-and stress-dependent viscosity: Δη T and Δη p are the viscosity contrasts with respect to temperature and depth (with z the vertical coordinate).σ is the yield stress and Ė is the second invariant of the strain-rate tensor (Stein et al., 2004;Tosi et al., 2015).

Model Setup
We employ the 2D Cartesian version of the finite volume code "plaatjes", where the conservation equations are solved with the multigrid method and the fully implicit time-stepping Crank-Nicolson method (Tosi et al., 2015;Trompert & Hansen, 1996).
The initial and boundary conditions are summarized in Figure 1 and key parameters with their dimensional values are given in Table 1.
The boundary conditions for the temperature are isothermal at the top (T = 0) and bottom (T = 1) to represent the cold surface and the hot core-mantle boundary, respectively.Following Trim et al. (2014), in dimensional values we have T 0 = 0°C at the surface and T CMB = 2500°C at the core-mantle boundary (Table 1).At the sides we choose reflective boundary conditions so that there is no heating or loss of heat from the sides.In order to ensure that these side conditions do not affect the results, we tested different box widths ranging from a square box (aspect ratio 1) to an aspect ratio 4 model geometry.The small box allows for a large number of simulations, while the larger box allows for more plates and subduction zones to form.The overall system behavior is the same in all tested box widths.Similarly, we varied the mesh resolution from 64 × 64 up to 512 × 256 and all parameter sets apply a Chebychev grid refinement toward the upper and lower box boundary to resolve the plate and pile thicknesses (cf.Stein et al., 2004).
As the initial mantle temperature is an important factor for surface motion (O'Neill et al., 2016), we vary the starting temperature T init in our systems.We either employ a linear temperature profile as shown in Figure 1a or a temperature which is a fraction of the core-mantle boundary temperature: T init = f • T CMB , with f varying between 0 and 1.The latter results in a constant initial temperature throughout the box that mimics an initially cold or hot mantle.A cold mantle will not be likely for the Earth, but is shown here to cover the full physical range of the system.The constant hot mantle is an approximation to the temperature profile of mantle convection featuring stagnant lid convection or radioactive heating, where the interior temperature is constant apart from temperature gradients in the boundary layers (e.g., Grasset & Parmentier, 1998;O'Farrell & Lowman, 2010;Tosi et al., 2015).In our simulations, these gradients develop with the start of the convection, thus prior to the onset of plate motion.Additionally, we analyze models with an initial linear temperature profile throughout the entire system (cf.Table 2).
The chemical boundaries are insulating at all sides so that here the dense material represents primordial material from early fractionation or differentiation processes.In this study, we neglect dense material being introduced into the solidified mantle from the core by core-mantle interaction (Stein & Hansen, 2023) or at the surface by impacts (Borgeat & Tackley, 2022).
We prescribe this primordial layer with an overall initial composition of C max = 1 (cf.Figure 1b), thus a given density difference Δρ according to the chosen buoyancy ratio (cf.Equation 6).The thickness d L (or volume V) of this layer is given as percentage of the mantle thickness d (mantle volume).z L gives the depth at which we initially place the base of a dense layer (cf. Figure 1b).For example, for z L = 1,740 km a layer with V = 10% (d L = 0.1 • d = 290 km) extends from 1,450 km down to 1,740 km.
The employed free-slip boundary conditions for the velocity allow for the free movement of surface plates.

Styles of Surface Motion and Plate Diagnostics
Due to a strong temperature dependence of the mantle viscosity, the cold surface layer of the system is highly viscous (rigid) and can completely cover the convecting mantle without taking part in the convection itself (Solomatov, 1995).This so-called stagnant lid can break by deformation processes, which are simulated by a reduction in lithospheric strength (controlled by the viscosity, cf.Stein et al., 2013) as high convective stresses are inserted to the lid.Here, we use the yield stress (cf.Equation 10) as the limit at which the lid breaks, if the stresses become too high (e.g., Moresi & Solomatov, 1998;Stein et al., 2004;Tackley, 2000).
The lid can either break by a strong plume impinging the lid from below (bottom-up mode) or by a thickened stagnant lid causing gravitational instabilities (top-down mode).After the breakup, small parts or the whole surface layer start to move episodically or continuously.In the case of episodic surface motion, the system shows long periods of stagnant lid convection that are interrupted by short bursts of subduction events (e.g., Moresi & Solomatov, 1998).In the case of continuous surface motion, the surface layer either moves fluid-like (surface viscosity similar to the interior viscosity) or as a plate with a higher viscosity than the interior and with little deformation (e.g., Stein et al., 2004).
Besides these three regimes (mobile, episodic and stagnant lid) two further subregimes of the stagnant lid have been classified: the heat-pipe model for a high extrusion efficiency and the plutonic-squishy lid for a low extrusion efficiency (Lourenço & Rozel, 2023).In the latter, melt is intruded in the crust as plutons because the melt is below the depth of negative buoyancy, whereas it is higher in the heat-pipe model and erupts at the surface.
For simplicity we neglect melting in our models, but will provide a further discussion after presenting our results.
We employ two diagnostic values to analyze the surface behavior in our models with respect to mobility and the plate-likeness (cf.Rolf et al., 2022;Stein et al., 2004;Tackley, 2000).
The surface motion is evaluated by means of the surface mobility M, which is the ratio of average surface velocity v s to the average system value v rms (Tackley, 2000).The mobility for the stagnant lid is indicated by M < 0.1 and for surface motion by M > 0.9, meaning that the surface layer moves about as fast as the interior (cf.Kreielkamp et al., 2022).
Additionally, the plateness P gives the amount of surface defined as a plate (Stein et al., 2004;Tackley, 2000) and is computed as the surface region with less than 20% of the maximum surface strain . In this way, some deformation is allowed within a plate and also small numerical inaccuracies are permitted.
Regions with more than 20% of the maximum surface strain are plate boundaries and regions with ϵ S ≈ 0 are regarded as stagnant lid.Additionally, a stagnant lid has a viscosity of at least two orders of magnitude higher than the interior and a mobility below 0.1.A fluid-like surface is given, if the surface viscosity is similar to the interior viscosity but mobility is higher than 0.9, whereas plate-like motion also shows a higher viscosity contrast (high plateness).For mobilities between 0.2 and 0.8 and a high plateness, we classify the system as being in the sluggish lid regime (cf.Lenardic, 2018).
As onset time of surface motion we define the time, when the mobility M of a simulation rises above 0.9 for the first time.In this way we take into account that at first a rigid plate (stagnant lid with a low mobility) forms that breaks either from the top by gravitational instabilities or from below by a plume.By choosing 0.9 as critical value, we guarantee that a mobile surface fully developed.The lowest limit for the onset time would be defined by M ≥ 0.1 as critical value.This is the point where the stagnant lid starts to break, which correlates with the moment when the plateness starts to change from 0% (stagnant lid) to a higher plate-like value.Typically, at the lower bound the plateness has not reached its final value, so that we have chosen M ≥ 0.9 as critical value.Throughout the results we will only give this upper bound, but will also present this lower bound of the onset time in the final discussion.

Parameter Sets
For the exclusive analysis of the effect of chemical heterogeneities and the mantle temperature, we neglect further complexities (e.g., no internal heat sources or melting effects in our model).
We will begin with one exemplary model calculation featuring a high Rayleigh number (parameter set 1 in Table 2), that is followed by simulations in which we vary one specific parameter at a time in order to understand its influence.Therefore, we run a number of simulations with low Rayleigh numbers to reduce the computational effort.As a consequence of the low Rayleigh number, the onset time of surface motion is increased beyond realistic values, but still shows the tendency for a specific parameter.
Table 2 summarizes the parameters which we used in the simulations.The Rayleigh number of Ra = 5000 in parameter set 1 is defined at the surface and converts to a bottom Rayleigh number (≈Ra • Δη T /Δη p ) of about 1 • 10 8 , which is in the range assumed for the Earth's mantle (Bercovici & Mulyukova, 2021).As the Rayleigh numbers are about 1-2 orders of magnitude lower in the other parameter sets, we also decrease the yield stress in these sets (set 2-set 5) by 1-2 orders of magnitude to be in the same flow regime as in set 1. A higher yield stress would result in stagnant lid convection (cf.Stein et al., 2004).
The Lewis number in the field approach is limited due to numerical issues and thus much lower than assumed for the Earth (∼10 13 , Kellogg & Turcotte, 1987), but has given similar flow results as the tracer approach which uses an infinite Lewis number (Stein & Hansen, 2023).The consequence of the lower Lewis number is that our thermochemical structures will erode more quickly than in the Earth's mantle.
For today's mantle, it is assumed that thermochemical structures have a volume of V ≈ 8% (Cottaar & Lekic, 2016), but due to erosion the structures will have been larger in the past.A density of about 0.5% higher than the surrounding density is proposed (Lau et al., 2017), with some structures (ultra-low velocity zones, ULVZs) even having a density increase of up to 10% (Rost et al., 2005).In numerical studies a high density contrast is needed to explain the observed stability of thermochemical structures (e.g., McNamara et al., 2010), but Koelemeijer (2020) debates if the structures are indeed denser than their surroundings.Due to these uncertainties and in order to fully understand their physical effect, we vary the density contrast and volume over a wide range to cover the full spectrum of flow results.Similarily, we vary the interior temperature over the full range (0°C like the surface temperature to 2500°C like the core-mantle boundary temperature).
We will provide a more detailed discussion on the implications of this choice of parameters, neglected complexities and model limitations in Section 4, after presenting the results of our models.

Results
We first consider the scenario where a uniform basal layer of dense material has formed after the magma ocean phase (Ballmer, Lourenço, et al., 2017;Labrosse et al., 2007).As shown in the initial model setup in Figure 1b we therefore prescribe a dense layer at z L = 2,900 km.This is a commonly assumed setup in solid-state mantle convection (e.g., Lassak et al., 2007;Li et al., 2014;Citron et al., 2020;Li et al., 2022;Langemeyer et al., 2022).
In agreement with the earlier studies, that focussed on the formation and stability of dense structures, we also find that material from the dense layer is entrained by rising mantle plumes (cf. Figure 2a at t = 1.15 Gyr) and pushed aside by sinking mantle currents to form dense accumulations, so called thermochemical piles (cf. Figure 2a at t = 1.48 Gyr).Due to the entrainment, piles erode over time and the density constrast between piles and ambient mantle decreases (cf. Figure 2b).
Here, we study thermochemical structures with respect to their influence on the surface motion.Trim et al. (2014), for example, ad hoc added a basal dense layer to a convecting system with moving plates and found that this could stop plate motion.In the self-consistent approaches by Stein et al. (2020) and Kreielkamp et al. (2022), it was also shown that dense structures delay the start of (continuous) plate motion as we observe on Earth today.While  2020) on the dynamic topography beneath dense structures, we now focus on the first initiation of surface motion (cf.Section 2.2).In addition to the previous works, we here also investigate the effect of shallow dense structures and the mantle temperature.3).The marked onset time gives the time, when the stagnant lid breaks and the surface starts to move (M ≥ 0.9 after the initial stagnant lid phase).For model parameters see simulation set 1 in Table 2.

Evolution of Thermochemical Structures and Surface Motion
At first an exemplary thermochemical model calculation (set 1 in Table 2) showing the surface mobility, along with temperature, composition and density fields for illustration, is given in Figure 2. The temperature and density fields for the corresponding purely thermal model are displayed in Figure S1 in Supporting Information S1, and the plateness for both model calculations in Figure S2 in Supporting Information S1.
The thermochemical model starts with vigorous convection beneath a stagnant lid.This is reflected by the topmost cold and dense layer being only deformed at the base by small downwellings, but no large cold slab sinks deep into the modeled mantle.At the same time first plumes start to rise and entrain some of the dense material from the primordial layer (see marked points 1-5 at t = 1.15 Gyr).These plumes are much smaller than the one observed in the purely thermal model (Figure S1 in Supporting Information S1 at t = 0.15 Gyr) which reaches up to impinge and break the stagnant lid (Figure S1 in Supporting Information S1 at t = 0.3 Gyr).
At t = 1.48 Gyr the thermochemical model shows a cold, dense slab sinking down to the base of the mantle (Figure 2a), where it pushes the denser primordial material aside and forms a large accumulation at the left part of the model.New small plumes form in the area that heats up after slab break-off (see marked points 1-5 at t = 1.95 Gyr).Additionally, a large plume has formed that is anchored by a large pile (point 0 at t = 1.95 Gyr).
From t = 1.95 to 3.33 Gyr the lid grows thicker by conductive cooling at the surface (cf.Stein et al., 2004).In the snapshots at t = 3.33 Gyr the areas of the slab remnants have heated up and the surface is again completely covered by a cold lid.At the same time the dense, primordial material distributed to an almost uniform layer covering the whole base of the system.The formation of a flat dense layer in stagnant lid convection agrees with the findings of Trim et al. (2014).
Further periods of subduction events and stagnant lid convection follow (t = 10.16 and 11.52 Gyr).
The initial episodic behavior is also reflected in Figures 2b and 2c in the corresponding surface mobility M (and in Figure S2 in Supporting Information S1 in the plateness P).The red lines show the results for the thermochemical model simulation displayed in the snapshots in the top part of Figure 2.For comparison, the black lines give the results for the corresponding purely thermal convection model in Figure S1 in Supporting Information S1 (B = 0, Le = 1 but otherwise the same parameters).
The model calculation of thermochemical convection reflects a large period (up to 15 Gyr) of initial episodic behavior (red line in Figure 2b).The mobility M varies between high values (during subduction events, e.g. at t = 10.15Gyr) and low values that represent stagnant lid convection (e.g., at t = 1.95 Gyr).After the episodic phase, the simulation shows a period of erratic, sluggish behavior (M > 0.5): the mobility still strongly varies but does not reach down to the value of the stagnant lid mode of convection.Finally, at t ≈ 40 Gyr the statistically steady state, with a mean mobility value (and plateness) equal to that of the thermal model, is reached (Figure 2b; Figure S2 in Supporting Information S1).As described in detail in Kreielkamp et al. (2022), this state is reached much later in the presence of dense material in the mantle.Here, we additionally find that the onset time (first breaking of stagnant lid) for the thermochemical model (t ≈ 1.3 Gyr) is delayed compared to the thermal model (Figure 2c).
In the model of purely thermal convection, the stagnant lid quickly breaks and surface mobility increases.At first the movement is erratic, but becomes more stable after approximately t = 15 Gyr (black line in Figure 2b).The onset time of surface motion (M ≥ 0.9) is t ≈ 0.3 Gyr (black dot in Figure 2c).
The results suggest that initially plumes are weakened by the presence of chemical heterogeneities and insert less stress to the base of the lid.In the thermal model, we find a large plume impinging the stagnant lid and breaking it in a bottom-up manner.This is supported by the strain-rate field in Figure S3a in Supporting Information S1.In the thermochemical model, however, the lid becomes gravitationally unstable and breaks in a top-down manner, before a plume reaches the lid (cf.strain-rate field in Figure S3b in Supporting Information S1).The dense material at the base of the system hinderes plumes to rise to the surface.
The blue lines in Figures 2b and 2c give the density contrast between piles and the ambient mantle, which displays a decrease as dense material is entrained.This entrainment is also visible in the composition fields in Figure 2a as light blue streaks at the position of plumes (e.g., point 0 marked at time t = 1.95 Gyr).Consequently, the composition C of the piles is reduced from C = 1 at t = 1.48 Gyr to about C = 0.5 at t = 10.16Gyr.More clearly visible is the entrainment in the density fields.The density of the pile at x = 0-5,800 km has decreased by more than 70 kg/m 3 from t = 1.48 Gyr to t = 10.16Gyr.
With the end of the pile period (i.e., no piles left due to complete entrainment), the system enters the period of sluggish surface motion where plumes can largely rise unhindered (note the larger plumes at t = 11.52 Gyr compared to t = 1.15 Gyr in Figure 2a).The system, that initially cools with time due to thermal diffusion at the surface and by the cold slab mixing into the interior (cf.Figures S4a and S4b in Supporting Information S1), can then heat up with the rising plumes (cf. Figure S4a in Supporting Information S1).It takes a while for the system to reach the steady state at approximately 40 Gyr.At this point, the thermochemical model calculation shows a modeled mantle that is homogeneously mixed (cf. Figure S4c in Supporting Information S1).The observation that thermochemical plumes carry less heat into the mantle, because the thermochemical layer acts as a thermal insulator, agrees with the studies of Farnetani (1997) and Li and McNamara (2018).
These exemplary simulations show that the first onset of surface motion and the onset of steady surface motion can strongly be delayed in the presence of dense mantle material.Additionally, we observe that the presence of dense piles leads to initial episodic or erratic surface motion, but as piles erode over time, a change to continuous surface motion appears at a later time.
In purely thermal convection, a similar change from episodic to continuous surface motion is achieved by increasing the buoyancy in form of an increased Rayleigh number (cf.Stein et al., 2004).In this study, the buoyancy of the rising plumes is at first reduced, compared to within purely thermal convection, due to the dense material entrained by plumes (cf.Figures S4d and S5 in Supporting Information S1).This is similar to the observation of Dannberg and Sobolev (2015).With progressing entrainment (reduction of density contrast), the buoyancy increases and we observe the change from erratic to continuous motion as in the purely thermal model.
In the following, we study the effect of various starting scenarios in more detail by separately analyzing different densities and volumes.Additionally, different depths of a dense layer and different initial system temperatures are systematically investigated.

Effect of Thermochemical Piles (Density Contrast and Pile Volume)
We start by elucidating the initial density contrast and the initial volume of a basal primordial layer (Figure 3a) considering model calculations from parameter set 2 (cf.Table 2).We clearly find that the first collapse of the stagnant lid (i.e., the onset time of surface motion) is more strongly delayed if the density contrast is increased.The same we observe when increasing the volume of the initial dense layer.For example, for the density contrast of δρ = 9%, the onset time changes from 10 to 1,000s of Gyr when increasing the volume from 3% to 20%.For comparison, the purely thermal reference case for this parameter set gives an onset time of approximately 4.3 Gyr.Note again, that here onset times are much higher than realistic due to the low Rayleigh number, which we have chosen from hereon in order to perform a large number of model calculations to understand the underlying physics.
As discussed in Section 3.1, plumes are weaker when they entrain dense material.The increased restoring force (i.e., higher density) hinders the ascent of plumes.However, we observe that not only rising plumes are suppressed in thermochemical convection, but also downwelling instabilities.For the thermochemical model, the strength of convection is reduced in the whole system (cf.Figure S6 in Supporting Information S1).This also hinders the start of surface motion by the top-down manner.With increasing density or volume of chemical heterogeneities this effect is more pronounced and the onset of surface motion is delayed.
Interestingly, this seems only valid up to a certain density contrast.For parameter set 4 in Figure 3b, a late onset for high density contrasts (here >10%) is also visible, but the delay is smaller than for the low density contrasts and seems to be almost constant, independent of the density contrast or starting temperature.

Variation of Initial Temperature
Like O'Neill et al. ( 2016), we observe that the mantle temperature plays an important role.Chemical heterogeneities indirectly affect the mantle temperature by hindering the rising of plumes, therefore inhibiting mantle heating (cf.Section 3.1).Varying the initial temperature between 1% and 100% of the core-mantle boundary temperature has an even bigger effect than the heating by plumes.We find that for low density contrasts (<10%), a hot mantle (T init = 2500°C) leads to an earlier onset of surface motion than the cold mantle (Figure 3b).The same we observe for the purely thermal model (δρ = 0%): for T init = 2500°C the onset time is about 10 Gyr, which is increased to almost 100 Gyr for T init = 250°C.
In a hot mantle, the temperature-dependent viscosity contrast over the surface is sufficient to quickly form a stagnant lid, whereas in a cold mantle the system first has to heat up to establish a strong viscosity contrast (cf. Figure S7 in Supporting Information S1).For high density contrasts, the strength of convection is suppressed so strongly, that diffusion is the main heat and composition transport and late onset of surface motion is obtained independent of the initial temperature (cf. Figure S7a in Supporting Information S1 right).
Figure 4a reveals the importance of the initial temperature in a regime diagram for a large number of models (parameter set 4 in Table 2).Here, for δρ ≤ 6% the onset time decreases strongly, the hotter the interior.The same is true for the purely thermal model with no density contrast.Furthermore, all simulations independent of the volume of dense material show this behavior (orange and blue line in Figure 4b, set 2 in Table 2).For example, for the density contrast of δρ = 3.6% the volumes of 3% and 5% are shown and for both we find that a low starting temperature delays the onset time.This is even more pronounced with increasing layer volume (light blue line).We observe that the initially hot systems typically evolve from the top, that is, by conductive cooling the stagnant lid grows thicker and becomes negatively buoyant.The initially cold systems instead first need to heat up by plumes transporting and mixing hot material from the base of the mantle into the interior.Consequently, it takes longer for the stagnant lid to form.The larger the volume or the higher the density contrast, the more strongly plumes are hindered and therefore the more pronounced is the delay in onset time.In Figure 4b we present two temperature field snapshots for simulations with the same density contrast and layer volume (δρ = 3.6% and V = 3%) but different initial temperatures to illustrate this behavior.The upper snapshot for T init = 2300°C shows a hot interior with one plume reaching below the thin cold lithosphere which starts to sink at the left side of the model domain.The lower snapshot is taken from a model calculation that uses 25°C as starting temperature.In this case, the initially cold system first needs to be heated up (cf. Figure S7 in Supporting Information S1).Here, two hot rising structures are visible below a much thicker cold top area.As the system only slowly heats up from below (even slower in the presence of a thicker or denser primordial layer), the formation of a stagnant lid takes longer and as a result the onset of surface motion is later.The same behavior is observed in the purely thermal models in Figures 4a and 4b (cases with δρ = 0%).Here, the onset time also increases more strongly below an initial temperature of about 50% of the CMB-temperature, that is, ∼1250°C (Figure 4b).
For very strong density contrasts (δρ > 6% in Figure 4a) the differences in the onset times of cold and hot systems is less pronounced.These unrealistically high contrasts lead to an initial period of conduction because the restoring force of the primordial layer is too strong.All systems fall into a linear stratified temperature and composition field (cf.also Figure S7 in Supporting Information S1).When the density contrast is diminished due to chemical diffusion, convection sets in and allows for the onset of plate motion.

Variation in the Depth of an Initial Dense Layer
Typical mantle convection models assume a dense primordial layer at the base of the mantle (e.g., Kreielkamp et al., 2022;Lassak et al., 2007;Li et al., 2014).This represents a remnant of the early magma ocean that occurred after the giant impact (Labrosse et al., 2007;Solomatov, 2015).However, magma ocean models featuring strong rotation reveal that dense material only sinks to the bottom of the magma ocean at the poles, whereas particles at the equator can be kept in suspension (Maas & Hansen, 2019).These results suggest an inhomogenous solidification with respect to depth and latitude.Thus an (at least regionally) elevated layer could be possible.The idea of chemical heterogeneities at different mantle depths has also been matter of investigation in other studies.For example, Ballmer, Houser, et al. (2017) and Gülcher et al. (2020Gülcher et al. ( , 2021) ) propose the existence of Bridgmaniteenriched ancient mantle structures (BEAMS) at mid-mantle depths, which are a result of batch crystallization (Caracas et al., 2019).Furthermore, there are assumptions that the chemical heterogeneities within today's mantle are recycled crustal material (Brandenburg & van Keken, 2007;Frost & Rost, 2014;Jones et al., 2020;Niu, 2018;Trønnes et al., 2019), meaning that initially the dense material was located at shallower depths.Tusch et al. (2022), for example, argue for Hadean crustal restite or O'Neill et al. ( 2017) and Borgeat and Tackley (2022) have followed the idea of dense material being introduced into the Earth by smaller impacts after the giant moonforming impact.We therefore also analyze how chemical heterogeneities at shallower depths affect the initiation of plate motion (Figure 5).
O'Neill et al. ( 2017) have already demonstrated that impacts induce short-lived subduction events.Dense material due to impacts, however, would only act very locally.For comparison to the primordial basal layer scenario, we consider bands of dense material.In Figure 5a we display the onset time for models (set 5) with an initial dense layer (V = 10%, d L = 290 km) which is placed at different mantle depths.As explained in Figure 1, here the given depth z L defines the layer's lowest depth.Thus, the 290 km thick layer at z L = 1,740 km actually expands from 1,740 km up to 1,450 km.
Generally we observe that the onset is earlier, the higher the initial dense layer.For example, for δρ = 3% the onset time is roughly 537 Gyr for a deep dense layer (z L = 2,610 km), 374 Gyr for a layer at 1,740 km, and 237 Gyr for a shallow dense layer (z L = 580 km).
Interestingly, a slightly earlier onset is observed if the dense layer is placed right at the base of the system (e.g., for δρ = 3% the onset time is 458 Gyr for z L = 2,900 km, compared to 537 Gyr for z L = 2,610 km).Here, the thermal boundary layer is thicker than the basal composition layer (cf. Figure S8a in Supporting Information S1).
Therefore the dense material only partly influences the plumes which rise above this dense layer.In other words, the volume of dense material acting on the plumes is lower, which reduces the onset time as we have already discussed (Figures 3a and 4b).
The earliest onset is obtained for a dense layer located at the very top of the system.For δρ = 3% and z L = 290 km, the onset time is 6.2 Gyr which here is even earlier than for the purely thermal model (about 8.2 Gyr).This early onset can be explained by the dense material lying within the cold thermal boundary layer (cf. Figure S8b in Supporting Information S1).The strong negative buoyancy allows for earlier initiation of sinking slabs (top-down scenario).
We analyzed the bottom-up and top-down scenarios further by means of temperature field snapshots, here with additional corresponding viscosity-depth profiles.Figures 5b-5i give examples of the top-down (left figures) and bottom-up tectonics (right figures), respectively.In the top row, the model calculations of parameter set 3 in Table 2 feature a dense layer at the depth z L of 327 km with a hot starting field (left snapshot) and a cold starting field (right snapshot).
In case of an initially hot system, a thin cold layer forms that subducts (Figure 5b), in case of a cooler system no such thin cold layer is observed (Figure 5c).Rather a vast area of cold material is obtained that is punctuated with a single plume in the center of the system.The corresponding viscosity-depth profiles reveal that for the initially hot system, due to conductive cooling from the top, a strong stagnant lid (highly viscous lid) has formed that extends to a depth of roughly 725 km (Figure 5d), while for the cooler system the viscosity of the lid is less high (Figure 5e).The smaller viscosity jump over the top layer represents a weaker lid.
Likewise we observe for a dense layer located at the depth of 2,350 km.The hot system yields a stagnant lid that starts to subside at the left (Figure 5f + 5h), whereas the cold system does not show a stagnant lid before plumes reach the surface (Figure 5g + 5i).In fact, in this scenario no pronounced viscosity jump in the top is obtained.Due to a viscosity contrast of less than one order of magnitude compared to the interior, we find that the surface also behaves more as a sluggish lid.
The snapshots were taken at the time when the first surface motion initiates.For the two cold systems the depth of the dense layer seems to have little effect on the onset time.For the hot systems, however, a larger difference is visible.The system with the shallow dense layer shows an early onset (Figure 5b), but the system with a deep dense layer instead displays a later onset (Figure 5f).
The earlier onset of plate motion in models with a shallower dense layer can be explained by two observations.On the one hand, shallow dense material quickly drips down due to negative buoyancy (Figure S9 in Supporting Information S1).This can trigger an early subduction event.On the other hand, rising plumes are less hindered by dense material located at mid-depth (Figure S10 in Supporting Information S1).

Effect of Mantle Temperature
In order to understand some observations when analyzing the presence of chemical structures within the convecting mantle, which is the focus of this work, we also briefly compare initially cold and hot systems.Changing the initial mantle temperature yields a similar, but more pronounced effect, as the mantle heating by plumes.As reference also the model of purely thermal convection with a density contrast of 0% is shown.Panels (b) and (c) show snapshots of the color-coded temperature field for a layer extending from 240 km to a depth of 327 km for the initial temperatures of 2250°C and 750°C, respectively, with (d) and (e) corresponding viscosity-depth profiles.Panels (f) and (g) show snapshots of the color-coded temperature field for a layer ranging from the depth of 2,263-2,350 km for the two initial temperatures with (h) and (i) corresponding viscosity-depth profiles (parameter set 3 with δρ = 3.6%).
From our results we conclude that the initial interior temperature plays an essential role which agrees with earlier studies (e.g., Moyen & van Hunen, 2012;O'Neill & Debaille, 2014;Zegers & van Keken, 2001), where a hotter mantle is related to melting or delamination.Here, we simply compare initially hot and cold systems and neglect the additional aspects of melting and delamination.Nevertheless, as in the former studies, we obtain an earlier onset of surface motion if increasing the mantle temperature (cf.models 'temp' in Table 3).Like in Stein et al. (2013), who found that an increase in radiogenic, internal heating promotes the formation of a stagnant lid, the strong temperature difference between the cold surface and initially hot mantle quickly leads to the formation of a highly viscous lid in our thermoviscous convection models.The interior, however, only slowly heats up by plumes transporting the basal heat to shallower depths in our initially cold systems.At the onset time of surface motion, the viscosity contrast between the surface layer and the interior is still less than two orders of magnitude, so that we only have a very weak lid in initially cold systems.In some cases, we even obtain a contrast of less than one order of magnitude, meaning the surface layer behaves as fluid-like as the interior (Figure 5i).Additionally, the surface motion is very slow (cf. Figure S7c in Supporting Information S1) as the overall convective strength is strongly reduced (Figure S6a in Supporting Information S1).Consequently, the time for the formation of a stagnant lid and the subsequent surface motion is delayed compared to a hot mantle (Figure 4).
Considering that the Archean Earth's mantle was hotter than today (e.g., Jaupart et al., 2015), an early onset of surface motion seems therefore likely.Our model with a realistic Rayleigh number for today's mantle suggests that the onset of plate motion is possible within the first 0.5 Gyr after magma ocean solidification (Figure 2 and model "TC" in Table 3).Especially, if we take into account that the neglected factors (melting and delamination  Note.Here, we give the parameter sets with the varied parameter highlighted separately.Additionally, we give the figures in which the results are shown.As onset time we present the (upper) limit that we use in all the figures.This is when the mobility M exceeds 0.9 for the first time.Additionally, the lower limit (M ≥ 0.1), where the stagnant lid starts to break, is given.Furthermore, we provide the plateness P at the time of the onset with M ≥ 0.1.The plateness typically further increases until the steady state is reached (cf. Figure S2 in Supporting Information S1).
due to eclogitization) further promote the breaking of the stagnant lid (e.g., Moyen & van Hunen, 2012;O'Neill & Debaille, 2014).Also, the Rayleigh number was presumably higher than today and declined with the cooling of the core (cf.Davies & Greenwood, 2023;Labrosse et al., 1997).A higher Rayleigh number would also make early mobilization more like (cf.Stein et al., 2004).This is visible in our models "density" and "depth" with δρ = 0% in Table 3, where the onset time increases from 4.09 to 8.32 Gyr, when reducing the Rayleigh number by a factor of 4.
The convective strength in thermochemical convection is strongly affected by the presence of dense material.Observed chemical heterogeneities in Earth's mantle should therefore be taken into consideration when analyzing the onset time of plate motion.

Temporal Evolution
In this study we investigate the effect of dense basal structures on the onset of plate motion.We observe that at first the style of plate motion is more episodic or sluggish compared to the modern-style plate motion.This is either characterized by short episodes of plate motion during long phases of stagnant lid behavior or by a surface layer that shows slow movement (M ≈ 0.5 ± 0.3).Initial episodic behavior has been favored by many modellers (Condie, 2018;Debaille et al., 2013;O'Neill et al., 2016;van Hunen & Moyen, 2012) and was also suggested by studies of zircon data (e.g., Bedard, 2018;Kirkland et al., 2021).
During these episodes, dense material is entrained upwards, reducing the density contrast or volume of dense accumulations, which finally allows for a more stable plate motion (cf. Figure 2b) with deep subduction rather than cold drippings detaching from the stagnant lid (cf. Figure 2a).Accordingly, the time until stable modernstyle plate motion is observed, is prolonged the larger the density contrast and volume of a dense layer.For a more detailed discussion on the evolution toward stable plate motion, we refer to Kreielkamp et al. (2022).

Onset Time of Plate Motion
Here, we evaluate the moment, when the stagnant lid breaks the first time and the surface starts to move (cf.Section 2.2).Our models "TC", "density", and "volume" in Table 3 all suggest that deep dense material delays the onset of surface motion compared to the purely thermal model simulations with no density contrast.
Generally, deeply rooted dense material reduces the buoyancy of ascending plumes and thus hinders their rise (e.g., Figure S5 in Supporting Information S1), and with that the effective heating of the mantle (e.g., Figure S4 in Supporting Information S1).This thermal insulating effect of dense structures was also described by Farnetani (1997) and Niu (2018).Likewise, Baes et al. (2021) argue that plume-induced tectonics is difficult for strong piles.We observe that the presence of piles leads to weaker convection in the whole mantle.Too large values of the density can even reduce the strength of convection as much as to prevent convection completely (cf.Figure S7 in Supporting Information S1).A conductive phase results, in which the dense material slowly diffuses before plate motion is possible.
Assuming realistic values for the density contrast of LLSVPs of about 0.5% (Lau et al., 2017) and a volume of 8% (Cottaar & Lekic, 2016) or the higher density contrast of 10% but smaller volume for ULVZs (Rost et al., 2006) would increase the onset time by roughly one order of magnitude compared to the purely thermal model (cf.models "density" and "volume" in Table 3).Taking into account that the Rayleigh number was higher than in parameter set 2 and further effects such as melting were neglected, would nevertheless place the onset of plate motion late in the Archean or even as late as the Proterozoic with steady plate motion only occurring in the very recent time.
An earlier onset could then be in favor of less dense structures at the base of the mantle as argued by Koelemeijer (2020).Alternatively, the presence of less deep chemical heterogeneities (e.g., Ballmer, Houser, et al., 2017;Tusch et al., 2022) will affect the onset time differently.Due to dense material acting on the buoyancy of convective instabilities, the location of dense material strongly matters.While deep chemical heterogeneities increase the density of rising plumes and therefore reduce their buoyancy, chemical heterogeneities located at the surface increase the negative buoyancy and therefore enhance downwellings.We thus turned our attention in this study on the investigation of dense material at different depths to give estimates on the onset of plate motion.
As visible from the models "depth" in Table 3, the onset time strongly decreases for chemical heterogeneities placed higher in the mantle.Deeply situated dense structures hinder plumes in rising and heating the mantle, while dense structures at mid-depth allow plumes to rise higher (Figure S10 in Supporting Information S1).
Most notably is the strong reduction in onset time, if the dense material lies near the surface (Figure 5a and Figure S9 in Supporting Information S1).Similar to the impact-induced models (Borgeat & Tackley, 2022;O'Neill et al., 2017), we find that dense material located within the lithosphere leads to an early subduction event.This is due to an increase in negative buoyancy when the density of the surface layer is increased.For parameter set 4, we even find that the onset time can be earlier than in the purely thermal model (cf.models "depth" in Table 3).

Model Limitations
In our simple approach, which excludes phase transitions, shallow dense material sinks to the base of the mantle, where it then acts on the ascending plumes in a similar way as for a primordial basal layer (cf.Figures S9 and S10 in Supporting Information S1).Consequently, further plate motion is hindered.Like in Trim et al. (2014) who added a dense basal layer to convecting flow featuring active plate motion, we observe that after the very early subduction event a phase of sluggish and erratic behavior can follow.In Comeau et al. (2021) we find, that a phase transition hinders dense eclogite to sink and rather causes a slowly delaminating surface slab that buckles above the phase transition.Therefore, we assume that additional phase transitions might reduce the sluggish period, giving a more continuous surface motion.This will have to be studied in future works.
Further limitations of our models are the neglect of additional radiogenic heating or the usage of the 2D Cartesian box.Apart from the curvature effect, these limitations might not have a big impact as the geotherm in a Cartesian box is larger than in the spherical shell, which is compensated by a low or no internal heating rate in Cartesian models (O'Farrell & Lowman, 2010).
As mentioned above, melting will have been a relevant process in the hotter Archean mantle (e.g., Moore & Webb, 2013;Moyen & van Hunen, 2012).In our models, showing stagnant (or sluggish) lid behavior in the early evolution, this effect was neglected.However, the shielding from effective cooling in systems with a highly viscous lid increases the interior temperature, which enhances melting.Consequences could then be volcanic outgassing and the growth of a secondary crust as melts in the upper mantle rise, cool and solidify.Extrusive volcanism and magmatic intrusions were relevant forms of magmatic transport in the Hadean and Archean mantle (e.g., Lourenço et al., 2020;Moore & Webb, 2013).In the so-called plutonic squishy lid, magmatic intrusions weaken the lower part of the crust and fragment the lid allowing for slow surface motion (Lourenço et al., 2020).
Often eclogitic dripping or lithospheric delamination appear close to the weak areas.As such, the plume-lid tectonics (e.g., Fischer & Gerya, 2016;Gerya et al., 2015) closely resembles the plutonic-squishy lid scenario with the difference of being more regional and showing a weaker crust (lower plateness) than the plutonic-squishy lid regime (Lourenço & Rozel, 2023).Less movement occurs in the so-called heat-pipe model, where melt diapirs reach the surface (Moore & Webb, 2013;Stern et al., 2018).Both, extrusion efficiency and lithospheric strength (viscosity), determine the resulting tectonic regime (Lourenço et al., 2020).For example, a thick and strong crust could prevent extrusive magmatism and could therefore be in favor of plutonism rather than volcanism (Rolf et al., 2022).An earlier onset of plate motion than observed in our models might be a result of melt weakening the early stagnant/sluggish lid obtained in our models.

Model Parameters
The exact onset times depend on a number of controlling parameters, which are not well constrained for the early Earth.For example, a higher Rayleigh number promotes the movement of the surface plates (Stein et al., 2004), which we could confirm by our different parameter sets.Our thermal model with a high Rayleigh number (parameter set 1, model "TC" in Table 3) reveals an onset time of less than 0.5 Gyr, while the simulations of purely thermal convection with a lower Rayleigh number (which we have chosen here for reasons of computational time) lead to a later onset time (∼3-4 Gyr for parameter set 2, model "density" in Table 3).
Due to the low Rayleigh number, the yield stress in our models is lower than assumed Earth values (cf.Katayama & Karato, 2008;Kohlstedt et al., 1995), but was chosen to yield surface motion in all cases.The yield stress in our models scales to dimensional values of σ • κ T η 0 /d = 0.3-30 MPa (with parameters from Table 1 and κ T = 2.5 • 10 6 m 2 /s).Generally, an increase in yield stress reduces the surface mobility (cf.Stein et al., 2004), which would shift the onset of surface motion to a later time.For example, from model set 2 to model set 4 we reduced the yield stress by one order of magnitude (and the Rayleigh number by a factor of 4), which roughly doubles the onset time from ∼4-8 Gyr (cf.Table 3).Besides the strength of the lithosphere (viscosity contrast and yield stress), also the Lewis number will play a role.A higher (more realistic) Lewis number should lead to an earlier onset than in our simulations as there will be less entrainment of the dense material (Stein & Hansen, 2023).

Conditions for Archean or Proterozoic Plate Motion
In our numerical study on the onset time of plate motion we find that the mantle temperature plays an important role.In an initially cold mantle, the formation of the stagnant lid is very late and consequently the break-up into single rigid plates is delayed.For an initially hot mantle, as more likely for the Archean Earth (e.g., Herzberg et al., 2010), the onset time of plate motion is one order of magnitude lower.In our simulation with a Rayleigh number in the realistic range for today's Earth's mantle, we obtain an onset time of 0.15-0.3Gyr after the magma ocean solidification.
A basal layer with a density of up to ∼10% higher than the surrounding, as assumed for some dense structures at the core-mantle boundary (cf.Rost et al., 2005), increases the onset time by at least one order of magnitude.For a volume of about 3%, the onset time could then be in the range of 2-3 Gyr after magma ocean solidification.Assuming a volume of about 8% (cf.Cottaar & Lekic, 2016) or even higher in the past (due to erosion), would strongly increase the onset time further.This might hint at less dense structures as suggested by Koelemeijer (2020).
Another possible scenario could be chemical heterogeneities initially located near the surface (at the base of the stagnant lid).Our models show that shallow dense material sinks down to also form thermochemical piles.Basal structures with a volume and density as observed today could therefore additionally be explained by some of the dense material dripping down to the core-mantle boundary after solid-state convection has started.Drip-like events have been described, for example, by Gerya et al. (2015), Foley (2018) and Baes et al. (2021) and the idea that the chemical heterogeneities within Earth's mantle stem from recycled crust is favored by Mulyukova et al. (2015), Niu (2018) or Trønnes et al. (2019).In the case of initially shallow dense material, we observe that the onset of plate motion is again in the same range or can even be earlier than in the purely thermal model.

Summary
We conducted numerical models of solid-state mantle convection featuring a dense primordial layer as a remnant of the magma ocean period.In an extensive parameter study we investigated the influence of the depth, the density and the volume of this dense layer as well as mantle temperature on the initiation of plate motion.
We observe that a low mantle temperature and the presence of chemical heterogeneities lead to a late onset of plate motion, which is further delayed if increasing the volume, density and mantle depth of dense structures.
The onset of plate motion in the hot Archean mantle is therefore only feasible for basal structures with a low density contrast/volume or shallow chemical heterogeneities.In line with several other studies (e.g., Bedard, 2018;Debaille et al., 2013;van Hunen & Moyen, 2012), the early plate motion is then characterized by single subduction episodes interrupting periods of sluggish or stagnant lid behavior.With progressing entrainment of dense material by plumes, the erratic surface motion fades out and results in steady plate motion as we observe today.
has also been benchmarked against several other open-source) mantle convection codes.All data supporting the findings of this study are available in the TRR 170 repository Stein and Hansen (2024).

Figure 1 .
Figure 1.(a) Boundary conditions for all simulations, with T the temperature, C the composition, (u, w) the velocity components in horizontal (x) and vertical (z) direction.Here, the initial temperature profile is linear with zero at the top and 1 ( =2500°C, cf.Table 1 for dimensional values) at the bottom of our system (T init = (d z)/d • T CMB , with d the mantle thickness, z the vertical coordinate (height) and T CMB the core-mantle boundary temperature).(b) The initial composition field is zero (purely thermal model) plus a layer of dense material (thermochemical models).This layer has the composition of 1 ( = density difference of Δρ), a width of d L with its base located at depth z L .

Figure 2 .
Figure 2. (a) Temporal evolution of the color-coded temperature, composition and density field.(b, c) Density contrast δρ (blue line) and surface mobility M for the thermochemical model (red line)and the purely thermal model as reference (black line) for a modeled time of 50 Gyr and the first 4.5 Gyr, respectively.The density contrast is computed as the average over the piles (defined as areas with C > 0.3).The marked onset time gives the time, when the stagnant lid breaks and the surface starts to move (M ≥ 0.9 after the initial stagnant lid phase).For model parameters see simulation set 1 in Table2.

Figure 3 .
Figure 3. (a) Onset time as function of the initial density contrast and volume of the dense material (model parameter set 2 with T init = 2300°C).(b) Onset time as function of the density contrast for different initial temperatures (model parameter set 4).The simulations with a density contrast (and volume) of zero mark the corresponding purely thermal models.

Figure 4 .
Figure 4. Onset time as function of the initial temperature.(a) Regime diagram spanned by the initial system temperature and the density contrast (model parameter set 4).(b) Model calculations with different density contrasts and volumes of the initial dense layer along with a purely thermal reference case.The snapshots show the temperature field for two simulations with δρ = 3.6% and V = 3% and the initial temperatures of 2300°C (top) and 25°C (bottom).Here, model parameter set 2 has been used.In both figures purely thermal models (δρ = 0%) are given for comparison.

Figure 5 .
Figure 5. (a) Onset time as function of the depth z L of the dense layer and the density contrast (parameter set 5).As reference also the model of purely thermal convection with a density contrast of 0% is shown.Panels (b) and (c) show snapshots of the color-coded temperature field for a layer extending from 240 km to a depth of 327 km for the initial temperatures of 2250°C and 750°C, respectively, with (d) and (e) corresponding viscosity-depth profiles.Panels (f) and (g) show snapshots of the color-coded temperature field for a layer ranging from the depth of 2,263-2,350 km for the two initial temperatures with (h) and (i) corresponding viscosity-depth profiles (parameter set 3 with δρ = 3.6%).

Table 1
Key Model Parameters (in Bold Those Varied Within This Study) and Their Dimensional Value (cf.Trim et al., 2014)

Journal of Geophysical Research: Solid Earth
Kreielkamp et al. (2022)concentrated on the evolutionary path (number of episodes) toward this continuous plate motion andStein et al. (

Table 3
Key Results for Some Exemplary Runs