Coupled dynamics and evolution of primordial and recycled heterogeneity in Earth's lower mantle

The nature of compositional heterogeneity in Earth's lower mantle remains a long-standing puzzle that can inform about the long-term thermochemical evolution and dynamics of our planet. Here, we use global-scale 2D models of thermochemical mantle convection to investigate the coupled evolution and mixing of (intrinsically dense) recycled and (intrinsically strong) primordial heterogeneity in the mantle. We explore the effects of ancient compositional layering of the mantle, as motivated by magma ocean solidification studies, and of the physical parameters of primordial material. Depending on these physical parameters, our models predict various regimes of mantle evolution and heterogeneity preservation over 4.5 Gyr. Over a wide parameter range, primordial and recycled heterogeneity are predicted to co-exist with each other in the lower mantle of Earth-like planets. Primordial material usually survives as medium- to large-scale blobs (or streaks) in the mid-mantle, around 1000–2000 km depth, and this preservation is largely independent of the initial primordial-material volume. In turn, recycled oceanic crust (ROC) persists as large piles at the base of the mantle and as small streaks everywhere else. In models with an additional dense FeO-rich layer initially present at the base of the mantle, the ancient dense material partially survives at the top of ROC piles, causing the piles to be compositionally stratified. Moreover, the addition of such an ancient FeO-rich basal layer significantly aids the preservation of the viscous domains in the mid-mantle. Finally, we find that primordial blobs are commonly directly underlain by thick ROC piles and aid their longevity and stability. Based on our results, we propose an integrated style of mantle heterogeneity for the Earth involving the preservation of primordial domains along with recycled piles. This style has important implications for early Earth evolution and has the potential to reconcile geophysical and geochemical discrepancies on present-day lower-mantle heterogeneity.

1 Introduction

The lower mantle is the largest geochemical reservoir in the Earth's interior and controls the style of mantle convection and planetary evolution. Despite efficient stirring by vigorous convection over billions of years, the Earth's lower mantle appears to be chemically heterogeneous on various length scales (e.g., Allegre and Turcotte, 1986; van Keken and Ballentine, 1998). Knowing the distribution of this heterogeneity is key for assessing Earth's bulk composition and thermochemical evolution, but this remains a scientific challenge that requires cross-disciplinary efforts.

On relatively small scales (mm to km), the concept of a “marble cake” mantle has gained wide acceptance, emphasizing that much of the mantle is made out of recycled oceanic lithosphere, deformed into narrow streaks of depleted and enriched compositions (Fig. 1a). Partial melting in the upper mantle creates heterogeneity between basaltic (magma) and harzburgitic (residue) end-members, forming a physically and chemically layered oceanic lithosphere. Subsequent injection into the mantle during subduction causes a non-equilibrated mantle that is a mechanical mixture of basalt and harzburgite (Allegre and Turcotte, 1986; Christensen and Hofmann, 1994; Morgan and Morgan, 1999; Xu et al., 2008). Additional small-scale heterogeneity is introduced into the mantle by recycling of continental material through subduction and/or delamination. Such influxes of felsic material to the deeper Earth have been used to explain the heterogeneous isotopic compositions of mantle-derived rocks, particularly in terms of trace elements (Hofmann, 1997; Kawai et al., 2009; Stracke, 2012).

Figure 1 Several conceptual models of mantle compositional structure with heterogeneity on various scales. For explanations regarding these heterogeneity styles, the reader is referred to the text: (a) “marble cake” mantle, modified following Xu et al. (2008); Woodhead (2015); (b) thermochemical piles, modified after Deschamps et al. (2015); and (c) viscous “blobs”, modified after Becker et al. (1999); Ballmer et al. (2017a).

On larger scales (10–100 s of km), compositional heterogeneity may be preserved by delayed mixing of this marble cake with either intrinsically dense or intrinsically strong materials. Intrinsically dense materials may accumulate as piles atop the core–mantle boundary (CMB), as illustrated in Fig. 1b. In particular, the two large low-shear velocity provinces (LLSVPs) in the deep Earth are commonly thought to have resisted mantle mixing due to their thermochemical origin, with an excess density of a few percent compared to pyrolitic material (Hernlund and Houser, 2008). Hypotheses regarding the composition of these dense piles involve recycled oceanic crust (ROC) (e.g., Christensen and Hofmann, 1994; Hirose et al., 2005; M. Li et al., 2014), ancient dense material (e.g., Labrosse et al., 2007; Tackley, 2012; Bower et al., 2013; Y. Li et al., 2014), or a combination of the two (e.g., Tackley, 2012; Ballmer et al., 2016).

On the other hand, intrinsically viscous domains may survive mantle convection and be preserved as “blobs” in the mid-mantle over long timescales (Manga, 1996; Becker et al., 1999; Ballmer et al., 2017a; Gülcher et al., 2020), such as plums in the mantle “plum pudding” (Fig. 1c). The physical properties (e.g., high viscosity) required for long-term preservation of these blobs are thought to be caused by an enrichment in the strong lower-mantle mineral bridgmanite (Mg, Fe)SiO₃ (i.e., stabilized by an enrichment in silica). During crystallization of a deep magma ocean at high pressures, MgSiO₃ bridgmanite is the relevant liquidus phase over a wide range of conditions and may hence be fractionated, potentially giving rise to such ancient heterogeneity in the mantle (Elkins-Tanton, 2008; Boukaré et al., 2015; Xie et al., 2020). These viscous domains in the mid-mantle may potentially reconcile recent seismic observations of mid-mantle heterogeneity (e.g., Fukao and Obayashi, 2013; Jenkins et al., 2017; Waszek et al., 2018), as well as cosmochemical and geochemical observations that indicate that the lower mantle hosts an ancient primordial reservoir that is potentially enriched in SiO₂ with respect to the upper mantle (e.g., Jackson et al., 2010; Mukhopadhyay, 2012). This regime has been successfully reproduced in 2D spherical annulus convection models with composition-dependent rheology (Gülcher et al., 2020). Various styles of primordial heterogeneity preservation were found as a function of its physical parameters, and the survival of sharp-to-diffuse primordial domains in the mid-mantle was proposed for planet Earth (Gülcher et al., 2020). However, the employed rheology was simplified, and the models did not produce Earth-like tectonic behavior. Moreover, none of the models explicitly predicted the formation of thermochemical piles in the lowermost mantle, which are perhaps the most evident heterogeneities in Earth's lower mantle (e.g., Solomatov and Stevenson, 1993; Christensen and Hofmann, 1994; Li and Romanowicz, 1996).

Many studies have explored the formation and preservation of either ROC (intrinsically dense) or primordial (intrinsically dense and/or viscous) heterogeneity, but only a few (if any) have quantified mantle dynamics in the presence of different types of heterogeneity with distinct physical properties (e.g., Nakagawa and Tackley, 2014; Y. Li et al., 2014). Understanding the interplay between recycled and primordial heterogeneity is critical to validate, falsify, and/or integrate the various views on chemical heterogeneity in Earth's mantle (Fig. 1). The goal of the present contribution is to investigate mantle dynamics in the presence of different types of heterogeneity with distinct physical properties (i.e., intrinsically viscous and dense) using Earth-like numerical models of mantle convection. We establish multiple regimes of chemical heterogeneity that are dependent on the rheological properties of primordial material. The coexistence of viscous, primordial blobs in the mid-mantle with dense recycled (and possible ancient) piles in the lowermost mantle is robustly predicted in many of our experiments. This coexistence is surprisingly largely independent of the initial compositional layering set-up. Based on these results, we propose that Earth's mantle is in a hybrid state between marble cake and plum pudding heterogeneity styles. Finally, the results are put into context with recent discoveries from geochemistry and seismology and the long-term evolution of Earth.

2 Methods

2.1 Numerical technique and boundary conditions

In this study, we model mantle convection in two-dimensional spherical annulus geometry using the finite-volume code StagYY (Tackley, 2008; Hernlund and Tackley, 2008). Free-slip and isothermal boundary conditions are employed at the top and bottom boundaries, with the surface temperature fixed to 300 K and that of the CMB fixed to 4000 K. The numerical experiments are purely bottom heated (i.e., no internal heating). StagYY solves for the conservation equations for mass, momentum, energy, and composition on a staggered grid for a compressible fluid with an infinite Prandtl number. The computational domain is decomposed in 512 × 96 cells, in which around 1.2 million tracers, tracking composition and temperature, are advected (25 tracers per cell). We further impose vertical grid refinement near the boundary layers and near 660 km depth. The size of grid cells therefore varies between 40 and 80 km in the horizontal direction and between 10 and 35 km in the vertical direction. Resolution tests were carried out with up to double the number of grid cells in each direction and up to 50 tracers per cell (Supplement S2). In these tests, we did not observe any significant changes in the dynamics and heterogeneity styles of our models. In fact, lower-mantle heterogeneity is slightly to moderately increased with increasing resolution (as in Gülcher et al., 2020). Our estimates of heterogeneity preservation thereby remain conservative.

2.2 Rheology

In our models, the assumed viscous deformation mechanism is diffusion creep. Viscosity is governed by a simplified temperature-, pressure-, and composition-dependent Arrhenius-type viscosity law (Newtonian rheology):

$$\begin{array}{}\text{(1)}& \mathit{\eta }(T,P,c)={\mathit{\eta }}_{\text{0}}\cdot {\mathit{\lambda }}_{\text{c}}\cdot \exp \left({\displaystyle \frac{E_{\text{a}}+P\cdot V_{\text{a}}}{R\cdot T}}-{\displaystyle \frac{E_{\text{a}}}{R\cdot T_{\text{0}}}}\right),\end{array}$$

where η₀ is the reference viscosity at reference temperature T₀ (=1600 K) and zero pressure; E_a and V_a are the activation energy and volume, respectively; T is the absolute temperature; P is the pressure; and R is the gas constant (8.314 J mol⁻¹ K⁻¹). The composition dependency of viscosity is considered through prefactor λ_c, which can be phase dependent (≠1 for primordial and post-perovskite compositions; see Table 1 and Sect. 2.5). All physical and rheological parameters used in this study are listed in Table 1. In cases where a grid cell contains a mixture of compositions, the overall cell viscosity is calculated as a mass-weighted geometrical average of the viscosities of different compositions (as in, e.g., Tackley and King, 2003).

Table 1 A summary of the physical properties used in the numerical simulations of this study. UM stands for upper mantle, LM stands for lower mantle, and PPV stands for post-perovskite. As we solve for compressible convection, the adiabatic temperature, thermal conductivity, thermal expansivity, density, and heat capacity are pressure dependent, following a third-order Birch–Murnaghan equation of state (Tackley et al., 2013). The asterisk (*) notes the parameters that are varied between the models.

The effective rheological properties of lower-mantle materials may be well represented by our relatively low effective activation energy and volume (e.g., Yang and Gurnis, 2016). For upper-mantle materials, the applied rheological parameters are significantly smaller than lab measurements (Karato and Wu, 1993; Hirth and Kohlstedt, 2003). In our Newtonian rheology formulation, low activation energies can account for the effects of grain size or stress-dependent rheology (e.g., Christensen and Hofmann, 1994; van Hunen et al., 2005; Yang and Gurnis, 2016). Low values of activation volume may result in reduced mechanical decoupling of the lithosphere and asthenosphere in our models, although our choice of yield stress (Table 1) allows for the development of stiff, mobile plates in our models (see below).

In order to obtain plate-like behavior at the surface, we assume that the material deforms plastically once a critical pressure-dependent yield stress is reached (as in Tackley, 2000; Crameri and Tackley, 2014). Plate-like behavior is determined using the same diagnostics as in Tackley (2000), measuring plateness P (the degree to which surface deformation is localized) and mobility M (the extent to which the lithosphere is able to move). For details regarding these diagnostics, the reader is referred to Supplement S3. Plate-like behavior occurs for P close to 1 and M close to or larger than 1 (Tackley, 2000).

2.3 Mantle composition, phase changes, and melting

The composition in our modeled mantle has three lithological end-member components: harzburgite, basalt, and primordial material. Each tracer carries either a primordial-material composition or a mechanical mixture of harzburgite and basalt (hz and bs). Harzburgitic and basaltic mantle materials are treated as a mixture of olivine and pyroxene–garnet systems (for details, see Nakagawa et al., 2010), following a parametrization based on mineral physics data (Irifune and Ringwood, 1993; Ono et al., 2001). Primordial material is not defined in terms of a specific mineral composition but solely through its material properties. It corresponds to a bridgmanite-enriched material with a (Mg + Fe) $/$ Si ratio of ≈ 1.0 and is thought to represent cumulates resulting from the fractional crystallization of a magma ocean in the early Earth (Ito and Takahashi, 1989; Elkins-Tanton, 2008; Xie et al., 2020). Phase transition parameters and physical properties for each mineral system and for primordial material are given in Table 2. The density profiles of the mantle materials are shown and discussed in Supplement S1.2. Our reference primordial material roughly corresponds to Mg_0.88Fe_0.12SiO₂ (Mg number of 0.88) or any other material with the same density profile (Supplement S1.2). The modeled mantles are initialized with a bottom primordial layer below a ≈ pyrolitic layer (see Sect. 2.4 below). The initial chemical density difference between such a primordial and pyrolitic layer stabilizes chemical stratification, while the thermal density difference destabilizes chemical stratification. The competition between the two is expressed as the non-dimensional buoyancy ratio B (Hansen and Yuen, 1988; Davaille, 1999):

$$\begin{array}{}\text{(2)}& B={\displaystyle \frac{\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{C}}}{\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{T}}}}={\displaystyle \frac{\mathrm{\Delta }{\mathit{\rho }}_{\mathrm{C}}}{\mathit{\rho }\mathit{\alpha }\mathrm{\Delta }T}},\end{array}$$

where Δρ_T and Δρ_C are the appropriate thermal and compositional density contrasts, α is the thermal expansivity, ρ is the lower-layer density, and ΔT is the super-adiabatic temperature contrast between surface and CMB. B is calculated for relevant lower-mantle depths, taking depth-dependent parameters ρ (density of primordial material), Δρ_C, and α at 1500 km depth. Accordingly, our reference model with a Mg number of 0.88 corresponds to a buoyancy number of 0.28.

Table 2 Parameters for the phase transitions for the olivine, pyroxene–garnet, and primordial systems. The primordial system is parameterized to fit the density profile of a mechanical mixture of 40 % basalt and 60 % harzburgite from Xu et al. (2008) in the upper mantle (see Supplement S1.2). These parameters include the depth and temperature at which the phase transition occurs, the density jump across the phase transition (Δρ_pc), and the Clapeyron slope (γ). The olivine and pyroxene–garnet system parameters are similar to previous studies (e.g., Tackley et al., 2013). Finally, the reference bulk modulus for each individual layer of a system is shown as K₀, which is increased to 230 GPa for primordial material in the lower mantle, as discussed in Supplement S1.2.

Changes in composition, carried on tracers, arise from melt-induced differentiation. Tracers in the basalt–harzburgite space undergo partial melting as a function of composition, pressure, and temperature, which generates basaltic oceanic crust and a complementary depleted residue (for details, see Nakagawa et al., 2010). Tracers of primordial composition undergo “melting” once they cross the relevant solidus of pyroxenite melting (Pertermann and Hirschmann, 2003). This solidus is used as pyroxenes are the low-pressure polymorphs of bridgmanite. During such a “melting” episode, the primordial tracer is converted into a tracer with 40 % basalt and 60 % harzburgite, and its primordial flag is removed. Any melting and related degassing would likely reduce, or even destroy, the ancient isotopic fingerprint of the previously primordial material (Gonnermann and Mukhopadhyay, 2007). The 40:60 basalt $/$ harzburgite ratio is used as it corresponds to a (Mg + Fe) $/$ Si ratio of ≈ 1.0, such as that in bridgmanite. Due to this tracer conversion, the final composition of the non-primordial mantle (from here on termed “ambient mantle”) $f_{\text{bs,amb}}^{\text{final}}$ evolves over time as primordial material is processed (melts) in the upper mantle:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{f}_{\text{bs,amb}}^{\text{final}}& =\left(f_{\text{bs,amb}}^{\text{ini}}\cdot V_{\text{amb}}^{\text{ini}}+c_{\text{bs,prim}}\cdot V_{\text{prim,molten}}^{\text{final}}\right)/\\ & \left(V_{\text{amb}}^{\text{ini}}+V_{\text{prim,molten}}^{\text{final}}\right),\end{array}\end{array}$$

where $f_{\text{bs,amb}}^{\text{ini}}$ is the initial bulk composition of the ambient mantle with volume $V_{\text{amb}}^{\text{ini}}$, c_bs,prim is the conversion factor of primordial material into bs–hz space (=0.4) and $V_{\text{prim,molten}}^{\text{final}}$ is the total volume of primordial material that has melted during model evolution ($=\left(\mathrm{1}-{\mathit{\chi }}_{\text{prim}}^{\text{pres}}\right)\cdot V_{\text{prim}}^{\text{ini}}$ with a primordial preservation factor of ${\mathit{\chi }}_{\text{prim}}^{\text{pres}}=\frac{V_{\text{prim}}^{\text{final}}}{V_{\text{prim}}^{\text{ini}}}$).

2.4 Initial set-up of the models

The initial mantle temperatures are calculated from an adiabat with a potential temperature of 1900 K together with the top and bottom boundary layers and small random perturbations. The resulting temperature profile is ≈ 300 K warmer than that Earth's present-day geotherm. Such high initial mantle temperatures are applied to roughly mimic the thermal evolution of the mantle, with high vigor of mantle convection and near-surface melting that likely occurred in the early Earth (e.g., Davies, 2007; Herzberg et al., 2010).

In terms of composition, our models are initialized as simplified two-layered (and for selected cases, three-layered) profiles, motivated by a fractional-crystallization sequence of the magma ocean (Ito and Takahashi, 1989; Elkins-Tanton, 2008; Boukaré et al., 2015; Xie et al., 2020). Similar to previous work (Gülcher et al., 2020), we impose a bridgmanitic “primordial” material layer in the lower mantle with layer thickness D_prim (1568–2230 km). By including randomly distributed 5 % pyrolitic “noise” in this primordial layer, the layer does not represent a pristine fractional-crystallization end-member cumulate. The upper mantle is initially a mechanical mixture of 15 % basalt and 85 % harzburgite, i.e., close to pyrolitic composition (Supplement S1.1). For some models, a third layer is added just above the CMB (150, 200, or 250 km thick) for which we impose the same physical properties as for basalt. This intrinsically dense layer is thought to represent an iron-enriched material, e.g., originating from the last cumulates from a crystallizing basal magma ocean (e.g., Wang et al., 2021) or from ancient crust that settled at the CMB (Elkins-Tanton et al., 2003).

2.5 Parameter study

In this study, we systematically investigate the styles of chemical heterogeneity in the mantle as a function of the physical properties of primordial material. The physical properties explored are the intrinsic density (FeO enrichment in (Mg, Fe)SiO₃) and viscosity contrast (silica enrichment) relative to pyrolite (as in Ballmer et al., 2017a; Gülcher et al., 2020). The density profile of primordial material is shifted throughout the mantle, which intrinsically changes the initial buoyancy ratio B between the models (varying between 0.07 and ≈ 0.57, as per Eq. 2). The approximate Mg number of the primordial material thereby ranges from 0.9 to 0.86 (Supplement S1.2). We further impose a viscosity contrast λ_prim (λ_c in Eq. 1) between primordial and pyrolitic mantle material in the lower mantle that is motivated by the high intrinsic strength of bridgmanite relative to ferropericlase (Yamazaki and Karato, 2001; Girard et al., 2016).

Additionally, we varied the initial conditions of chemical layering in the models (Supplement S1.1). The initial layer thickness of primordial material D_prim was changed such that the final composition of the ambient (i.e., convecting) mantle is similar for all cases, with a target value of a 25 % basalt + 75 % harzburgite mechanical mixture ($f_{\text{bs,amb}}^{\text{final}}=\mathrm{0.25}$). This value is appropriate for the Earth, particularly if considering that the ambient lower mantle may be moderately enriched in ROC compared to the upper mantle in the present day (e.g., Murakami et al., 2012; Mashino et al., 2020; Yan et al., 2020). As $f_{\text{bs,amb}}^{\text{final}}$ for any given model depends on the extent of primordial-material processing over time, which creates basaltic materials through near-surface melting (Eq. 3), D_prim is set to a different value for a given combination of parameters B and λ_prim. The applied D_prim ranges from 1569 to 1844 km (i.e., 60 vol %–75 vol % of the lower mantle) between the models. These choices for the initial conditions of our models lead to final ambient-mantle bulk compositions of $f_{\text{bs,amb}}^{\text{final}}=\mathrm{0.24}$–0.26, i.e., very close to our target value for all cases. In order to calculate these appropriate values for D_prim in the first place, we ran an additional suite of test models with a constant initial volume of primordial material. From the predictions of these test models and Eq. (3), the appropriate values for D_prim for each case of the main model suite were computed. For more details on this approach and the results of the test suite, see Supplement S4.

Finally, we compare model outcomes of the selected cases which include a thin, dense FeO-rich layer just above the CMB (three-layered models) with those of corresponding models in the main model suite (two-layered model).

3 Results

We have conducted a total of 119 numerical experiments for this study. The relevant model parameters and selected output variables of each case are summarized in Supplement Tables S1 and S2. We first describe the various heterogeneity styles identified in the models (Sect. 3.1) and their dependence on composition-dependent rheological parameters (Sect. 3.2). We then discuss how initial compositional layering affects mantle dynamics and chemical heterogeneity preservation (Sect. 3.3). Finally, we focus on a specific heterogeneity style and analyze the characteristics of the distinct chemical mantle domains and how they are linked to mantle flow (Sect. 3.4).

3.1 Mantle heterogeneity styles

Our results reveal different styles of long-term convection and mixing in the mantle. We label these regimes in line with previous work on variable styles of primordial heterogeneity preservation (Gülcher et al., 2020) as regimes I–III. For all main models, primordial and ROC heterogeneity can coexist in a fully convecting mantle (partial heterogeneity preservation, regime III), with several distinct preservation styles (Fig. 2). A subset of the test experiments predicted various styles of stable chemical layering (regime II), which are summarized in Sect. 3.3 and discussed in more detail in Supplement S4. In contrast to previous work (Gülcher et al., 2020), we do not observe any models with pervasive convective mixing of ROC and primordial mantle materials (regime I) in this study.

Figure 2 The left columns show mantle sections for all sub-styles in regime III (partial chemical heterogeneity) at 4.5 Gyr model time; λ_prim, B, and D_prim are as labeled. The right columns show the corresponding profiles of primordial fraction, basaltic fraction, temperature, and viscosity. These profiles are radially averaged and time-averaged (4.25–4.5 Gyr). The shaded region indicates the range for all models in any given sub-style. Dashed lines refer to the case shown in panel (a).

3.1.1 General geodynamic trend

All experiments in the first model suite are characterized by whole-mantle convection and moderate compositional heterogeneity preservation in the lower mantle after 4.5 Gyr of model evolution. Plate-like tectonic behavior is displayed for most of their evolution (Supplement S3). Their evolution commonly starts as a two-layered convective system that is sustained for several 100 Myr. Almost immediately, weak downwellings develop from the cold top thermal boundary layer and are deflected at the compositional interface in the mid-mantle. The onset of convection in the lower layer is delayed with respect to the upper layer. However, progressive heating and cooling of the lower and upper layers promotes a whole-mantle overturn at ≈ 0.2–1 Gyr. During this overturn, much of the primordial material in the lower mantle may reach the upper mantle and can subsequently be processed by extensive near-surface melting. Consequently, basaltic crust is formed that can re-enter the mantle through subduction. Moreover, the strong primordial material is separated into several blobs and/or streaks that are slowly eroded and entrained by the convecting pyrolitic mantle over time. The final extent of heterogeneity differs between the models, and the following four sub-regimes are established.

3.1.2 Viscous blobs (regime III.B)

In many numerical experiments, primordial material is ultimately preserved as blobs in the mid-mantle. Here, the primordial material is neutrally buoyant (Supplement S1.2) and relatively stable. Deformation due to mantle convection is mainly focused along narrow upwelling and downwelling conduits (Fig. 2a, b), and blobs are variably preserved in poorly mixed regions in between. These blobs periodically rotate and/or are displaced laterally as they are passed by sinking slabs. Occasionally, two blobs coagulate or are separated again as the convection patterns reorganize through time. When slabs reach the CMB they push away basaltic piles to the side, disconnecting the otherwise continuous basaltic layer at the base of the mantle. The final lower-mantle primordial (${\mathit{\chi }}_{\text{prim}}^{\text{LM}}$) and ROC (${\mathit{\chi }}_{\text{ROC}}^{\text{LM}}$) heterogeneity range from 8 vol % to 40 vol % and 1.5 vol % to 6 vol %, respectively. This style of material preservation is similar to that described in Ballmer et al. (2017a); Gülcher et al. (2020) and is here established with the presence of ROC (as piles atop the CMB and marble cake streaks throughout the mantle).

3.1.3 Marble cake mantle (regime III.M)

In this style, the final mantle composition is characterized by the preservation of a relatively large amount of ROC (4.5 vol %–6.5 vol %) and low primordial heterogeneity (<5 vol %). In addition to large ROC piles in the lowermost mantle, chemical heterogeneity takes the form of thin streaks of recycled and primordial material in an otherwise well-mixed mantle. Ultimately, the mantle is mostly made up of a marble cake with both ROC and primordial streaks (Fig. 2c).

3.1.4 “Diffuse” primordial domains (regime III.D)

Here, primordial heterogeneity is not preserved as coherent blobs or streaks (e.g., with sharp boundaries as in regime III.B or III.M) but is rather diffusely distributed throughout much of the lower mantle (Fig. 2d). Even though no discrete primordial domains are sustained in the convecting mantle, significant amounts (2 vol %–7 vol %) of primordial heterogeneity can survive in the lower mantle for 4.5 Gyr model time. As in regime III.B, ROC is preserved as piles in the lowermost mantle, as well as streaks throughout the mantle. Aside from these additional structures, regime III.D is similar to the diffuse-domain regime detected in Gülcher et al. (2020).

3.1.5 Metastable primordial piles (regime III.P)

In some numerical experiments, primordial material is preserved as large primordial blobs and streaks that are mostly confined to the lowermost mantle (Fig. 2e). Initial model evolution differs from the main geodynamic trend as hot upwellings form in a delayed manner, and even when these upwellings reach the compositional interface, a double-layered convective system is maintained for several 100 Myr. Gradually, the two-layered system breaks down as subducting slabs break through the lower layer, and in return hot upwellings, including primordial material, reach upper-mantle depths. Large viscous blobs are confined to, or settle near, the CMB due to their relatively large negative buoyancy. Some of the blobs are occasionally pushed up from the CMB by convective currents and intermittently float through the mid-mantle. In contrast to other partial heterogeneity styles, short intervals with low mobility occur in this sub-regime (Supplement S3). With the addition of the presence of dense ROC heterogeneity (1.5 vol %–4 vol %) throughout the lower mantle as thin streaks, this regime is similar to previously described as metastable piles (regime III.P) in Gülcher et al. (2020).

3.2 Influence of physical parameters of primordial material

The regime diagram shown in Fig. 3 pinpoints the parameter space where the various heterogeneity styles discussed above are identified. The viscous blobs heterogeneity style (III.B) occurs for the majority of the parameter space (≈ λ_prim>20 and variable B), while the total volume of the preserved blobs strongly depends on the intrinsic viscosity contrast λ_prim. High λ_prim results in a delayed overturn, more localized deformation, and slow entrainment of primordial material by the convecting mantle (Figs. 3, 2a). In contrast, low λ_prim causes an earlier overturn, increases convective vigor and mixing, and prevents large primordial domains from being sustained in the convecting mantle: a marble cake mantle (III.M) occurs for low B, while primordial material is diffusely distributed throughout the mantle (III.D) for moderate to large B. In a narrow parameter space (high B and λ_prim of approximately >0.45 and >300, respectively), metastable, primordial “piles” settle due to their intrinsic negative buoyancy (regime III.P).

Figure 3 Summary of key results as a function of physical parameters of primordial material for model suite 1 (see text). The vertical axis represents the initial buoyancy ratio B, defined at 1500 km depth (Sect. 2 and Supplement S3), or the corresponding Mg number of primordial material. The horizontal axis gives the viscosity contrast λ_prim between primordial and pyrolitic material in the lower mantle. The fraction of primordial heterogeneity (${\mathit{\chi }}_{\text{prim}}^{\text{LM}}$, red) and ROC (${\mathit{\chi }}_{\text{ROC}}^{\text{LM}}$, blue) in the lower mantle is averaged between 4.25 and 4.5 Myr model time (for definitions, see Supplement S1.3). The circle denotes the reference model in this study. Regime boundaries are established based on the style of mantle evolution and chemical heterogeneity: (III) partial heterogeneity preservation as marble-cake mantle (M), marginally stable piles (P), viscous blobs (B), or diffuse domains (D).

3.3 Influence of initial layering

3.3.1 Primordial layer thickness

The model suite discussed above was also run with a fixed primordial layer thickness (D_prim=2230 km, i.e., extending from the CMB up to the top of the lower mantle, as in Gülcher et al., 2020). The details regarding this model suite are described in Supplement S4. As a summary, the final ambient mantle composition greatly varies between the models, as very different amounts of primordial material are processed (as expected from Eq. 3). There is also more variety in geodynamic behavior, as two additional sub-styles of stable chemical layering (regime II) are identified (see Fig. S7): two-layered convection with an undulating primordial layer preserved for large B (II.T) and post-overturn ROC layering for low B and λ_prim (II.R). Although these regimes are not relevant for present-day Earth, they may be relevant for mantle dynamics within other (exo)planets with slightly different bulk composition, e.g., enhanced in silica. Comparisons of models with equal B and λ_prim in both model suites allows us to identify the effect of a different primordial layer thickness on model behavior. Selected model outputs are presented in Fig. 4. Final ROC heterogeneity is increased for most models with more primordial material initialized (Fig. 4a–b), and sub-regime boundaries are slightly shifted as primordial domains are more coherent in comparison with lower D_prim models (Supplement S4). Most notably, the fraction of initialized primordial material that is preserved after 4.5 Gyr of mantle stirring (${\mathit{\chi }}_{\text{prim}}^{\text{pres}}$) is similar in both cases or even lower in models with more primordial material initialized (Fig. 4c–d). Moreover, the final internal temperature of all experiments in the test model suite are systematically lower than those the main model suite.

Figure 4 Output quantities time-averaged between 4.25 and 4.5 Gyr for the main model suite (a, c, e) and test model suite (b, d, f): (a–b) primordial preservation factor ${\mathit{\chi }}_{\text{prim}}^{\text{pres}}$ in vol % of initialized primordial material, (c–d) average mantle temperature, and (e–f) convective vigor. The horizontal axis represents the primordial viscosity contrast λ_prim, while colors denote different B used in the models.

Figure 5 Temporal evolution of the radially averaged compositional profiles of selected models. Models with variable D_prim (main model suite) are on the left, and their corresponding model with fixed D_prim (=2230 km, test model suite) is on the right for given B and λ_prim: (a–b) B=0.28, λ_prim=100; (c–d) B=0.35, λ_prim=50; (d–e) B=0.21, λ_prim=30. The dotted black line indicates the onset of the compositional overturn. In the test model suite 0 (but not in the main suite), significant amounts of primordial material immediately reach the upper mantle after the overturn.

Figure 6 Snapshots during the compositional overturn phase for three models with equal B (0.28) and λ_prim (100) but a different initial compositional layering (as shown on the bottom left). The panels provide the potential temperature field and a zoom-in on the compositional field.

Figure 5 shows the temporal evolution of the radially averaged compositional profiles for selected experiments in both model suites. The compositional overturn (vertical dashed black line in Fig. 5) occurs later for models with a thicker primordial layer, as convective instabilities form less rapidly in the upper thinner layer to drive the overturn. Moreover, with a primordial-dominated lower mantle, primordial material is immediately swept to upper-mantle depths during the overturn while being displaced by the slabs that sink into the lower mantle (left panels of Fig. 5). Subsequent melting of primordial material and crustal recycling leads to a quick buildup of basaltic piles in the lowermost mantle. In contrast, a thinner primordial layer (left panels of Fig. 5) detaches from the CMB during the overturn and is swept to mid-mantle depths without immediately reaching the upper mantle. At these depths, the primordial material remains relatively stable. Melting of primordial material, crustal recycling, and ultimately the formation of ROC piles underneath is thereby delayed. This discrepancy between mantle dynamics and the extent of melting affects the subsequent thermal evolution of the models: melting is an efficient mechanism to remove heat from the mantle by transporting heat from the interior to the surface and enhancing the plates' mobility (Nakagawa and Tackley, 2012; Lourenço et al., 2016). Thereby, the lower extent of melting in models with a thinner bridgmanitic layer can explain their higher internal temperatures (Fig. 4c). The dynamic behavior described above can also be seen in Fig. 6, which shows detailed dynamics of selected models during the compositional overturn.

3.3.2 Additional basal FeO-rich layer

We further explore several additional cases that include an additional FeO-rich layer just above the CMB as an initial condition (Supplement S1.1). The results of these cases are summarized in Supplement Table S1. The imposed material properties of the FeO-rich layer are the same as for basalt (i.e., intrinsically dense; see Supplement S1.2). Notably, the ancient dense basal layer aids the preservation of bridgmanitic primordial material over time (higher ${\mathit{\chi }}_{\text{prim}}^{\text{pres}}$). During the compositional overturn, the FeO-rich layer remains at the base of the mantle, while the bridgmanitic material is swept to mid-mantle depths (Fig. 6c). Significant amounts of heat can be retained in this negatively buoyant layer that acts an insulating boundary layer above the core. Plumes readily form on top of this layer and are less hot than those otherwise directly forming above the CMB, resulting in a less vigorous compositional overturn and thus an increased coherency of primordial domains. Downgoing slabs separate the initially laterally connected dense layer into distinct piles, and they also transport ROC to the lowermost mantle, which is then added to the piles of ancient FeO-rich material (see last the two panels of Fig. 6c). The incoming ROC material is cooler and thus denser than the hot ancient FeO-rich material. Accordingly, it sinks through the ancient layer to settle directly above the CMB. Over time, much of the ancient FeO-rich material is progressively entrained by mantle upwellings and ultimately processed by near-surface melting. Only a small fraction is preserved after 4.5 Gyr, preferably located in the uppermost parts of some thermochemical piles (Fig. 7). The general trend in these piles is aging upwards: the youngest, most recently subducted (cool) ROC material is at the bottom, with progressively older ROC material progressing upwards, and the ancient (hot) FeO-rich material at the very top. Thereby, the piles are thermally stratified, and internal small-scale convection occurs in the lowermost, dense layer. Ultimately, the coexistence of large, coherent bridgmanitic primordial domains with basaltic piles (ROC and ancient) is robustly predicted, even for low λ_prim (Supplement S5).

Figure 7 Final snapshots (at 4.5 Gyr) of composition for selected models, λ_prim and B, as stated in the panels. The models are initialized with a 200 km thick FeO-rich basal layer (pink) overlain by a (a) 1546 and (b) 1514 km thick primordial layer. The panels also show a zoom-in on the formation age of ROC within the piles (defined as f_bs>0.6 and T_pot>2000 K). The formation age quantifies the time since the last melting episode.

Figure 8 Each panel shows a 2D histogram that relates several selected output quantities of a representative model displaying both primordial blobs and dense piles underneath (regime III.B, model M_300dD, Fig. 2a). For six time steps (4.25–4.5 Gyr with steps of 0.05 Gyr), selected quantities were calculated for each radial column in the model: pile height D_pile, formation age of pile material age_pile, mid-mantle primordial-material fraction f_prim, mid-mantle radial velocity v_z, and mid-mantle temperature anomaly dT. Thermochemical pile thresholds are f_bs>0.6 and T_pot>2000 K. Mid-mantle is defined as the depth range 1000–2000 km. The formation age of pile material is volume-averaged over the full pile thickness and weighted by the basaltic fraction. Each row in the histograms is normalized by the total counts of the row. For rows in greyscale, the number of counts in these rows represents less than 1 % of the total counts.

3.4 Relating dense piles, viscous blobs, and mantle dynamics

Finally, we analyzed the potential relationship between dense piles, viscous blobs, and mantle dynamics for a model that displays both piles and blobs (model M_300dD, regime III.B, illustrated in Fig. 2a). Figure 8 shows 2D histograms relating selected model output quantities with one another, including pile characteristics (pile height, formation age of pile material) and mid-mantle characteristics (primordial-material fraction, radial velocity, and temperature anomaly over a 1000–2000 km depth range). Several relations between these parameters can be inferred: first, the presence of primordial material in the mid-mantle is related with little radial dynamics (near-zero v_z, Fig. 8a). Similarly, most piles are stable at near-zero radial velocity, although thick (>100 km) piles are also stabilized in regions of upwellings. Second, the absence of primordial domains (f_prim=0, bottom rows in Fig. 8e, f) and piles (D_pile=0, bottom rows in Fig. 8e, f) is associated with columns of cool downwellings (v_z<0 and dT<0). Third, the presence of a large amount of primordial material in the mid-mantle is associated with positive temperature anomalies (Fig. 8b), and thus large primordial blobs (high f_prim) are typically slightly warmer than average (also see Fig. 2a, b). There is also a slight positive correlation between pile height and mid-mantle temperature anomaly: large piles are often related to overlying warm primordial blobs (see below). Although the histograms remain scattered, there is a slight positive trend in formation age of pile material versus mid-mantle primordial fraction and formation age versus pile height (Fig. 8c, g), indicating that piles containing older material are more often overlain by primordial domains than younger piles and that old ROC tends to be preserved near to top of the piles (see above). Finally, Fig. 8d and h indicate a spatial relationship (in radial direction) between the presence of primordial material in the mid-mantle and the height of the underlying ROC material. Even though this relationship is not strict, it indicates that the highest ROC piles tend to be centered just below the largest primordial blobs. This configuration away from upwellings and downwellings (Fig. 8b, f) appears to promote the coupled preservation of both these thermochemical domains.

4 Discussion

4.1 Styles of mantle mixing

Many of the identified dynamic styles of primordial-material preservation in our prior study (Gülcher et al., 2020) are here established in the presence of plate-like behavior and recycling of oceanic crust into the deep mantle. This indicates that primordial-material preservation robustly occurs for significant rheological contrasts, regardless of tectonic behavior or the presence of dense, recycled materials in the lower mantle. This study even shows the robustness of primordial-material preservation with respect to initial layering set-up. Despite no detection of models in which both primordial and recycled materials are efficiently mixed throughout the mantle (regime I in Gülcher et al., 2020), we expect such a behavior to occur for small rheological contrasts of primordial material in combination with lower-density anomalies of ROC than modeled here (Nakagawa and Tackley, 2014). Indeed, different ROC properties may affect the segregation, accumulation and survival of ROC in the deep mantle (e.g., Christensen and Hofmann, 1994; Hirose et al., 2005; Nakagawa et al., 2010; Mulyukova et al., 2015).

All of our numerical experiments are run in 2D geometry. Admittedly, the planform of mantle flow is different for 3D than for 2D geometry (Ferrachat and Ricard, 1998; Davies, 1990). An additional toroidal component of mantle flow in 3D can enhance mixing relative to mantle flow in 2D with its rather simple geometry of convection. In terms of the preservation of intrinsically strong heterogeneity, however, we expect competing effects to occur: on one hand, coherent blobs may be difficult to preserve in the presence of toroidal and poloidal flow components (Ferrachat and Ricard, 1998); on the other hand, mantle flow may be more efficiently guided around viscous blobs in 3D than in 2D geometry (Merveilleux du Vignaux and Fleitout, 2001). In addition, sheet-like upwellings in 2D can enhance mixing. Indeed, several studies showed that for high Rayleigh number convection models, the different geometry of boundary layer instabilities between 2D and 3D does not significantly modify the timescale or the regime of mixing (e.g., Coltice and Schmalzl, 2006; O'Neill and Zhang, 2018), although it is uncertain whether this would be the case in the presence of composition-dependent rheology as modeled here.

Figure 9 Sketch of our proposed style of present-day mantle heterogeneity in the Earth. Recycled crustal material (ROC, dark blue) is present throughout the mantle as marble cake streaks, with enhancement in the mantle transition zone (between 410 and 660 km depth). Moreover, the dense ROC material settles at the CMB to form thermochemical piles. Any ancient basal layer may be preserved within these piles, here proposed to be on the roof of the piles as predicted by our models (see Fig. 8). Finally, primordial (Mg, Fe)SiO₃-enriched material is preserved as blobs and streaks that reside in the mid-mantle (dark green) due to its intrinsic strength (BEAMS: bridgmanite-enriched ancient mantle structures, Ballmer et al., 2017a).

The relevance of each regime for terrestrial planet evolution depends on the actual initial compositional profile of the mantle at the onset of long-term solid-state convection. The initial conditions of our models are mainly motivated by an enrichment of a thick, basal layer in (Mg, Fe)SiO₃ bridgmanite (and hence in silica with respect to the overlying mantle). This may occur due to fractionation during magma ocean or basal magma ocean (BMO) crystallization (e.g., Labrosse et al., 2007). Indeed, bridgmanite is the liquidus phase in a magma ocean over a wide compositional and pressure range (e.g., Ito and Takahashi, 1987; Caracas et al., 2019). Alternatively, incomplete core–mantle equilibration during the last stages of planetary accretion (Deng et al., 2019; Kaminski and Javoy, 2013) or disproportionation of FeO during core formation (Frost and McCammon, 2008; Armstrong et al., 2019) may have promoted the formation of (bridgmanite-enriched) heterogeneity in the lowermost mantle. Moreover, chemical transfer between the outer core and the (basal) magma ocean in a hot early Earth may have caused delivery of SiO₂ to the BMO and extraction of iron from the BMO into the core, stabilizing bridgmanite as the primary liquidus phase during prolonged periods of BMO crystallization (e.g., Trønnes et al., 2019). Indeed, various studies have proposed SiO₂ exsolution from the core during secular cooling of the Earth or Earth-like planets (e.g., Badro et al., 2016; Hirose et al., 2017; Helffrich et al., 2018a; Rizo et al., 2019).

4.2 A recipe for Earth's lower mantle: marble cake plus plum pudding

For Earth, geodynamic styles that predict the coupled preservation of primordial material and ROC, such as regimes III.D (diffuse primordial domains), III.M (marble cake mantle), and III.B (viscous blobs), can potentially explain a wide range of geophysical, geochemical, and geological observations. First, these geodynamic styles are consistent with whole-mantle convection operating on Earth, which is motivated by seismic tomographic images of recently (<200 Myr) subducted lithosphere in the lowermost mantle (e.g., van der Hilst et al., 1997) and deep-rooted plumes that rise through the entire mantle (French and Romanowicz, 2015). Second, they all include the existence of thermochemical piles in the lowermost mantle, which are robustly featured in seismic tomographic studies as two antipodal, anomalously slow domains (LLSVPs) just above the CMB. (e.g., Li and Romanowicz, 1996; Ritsema et al., 1999; Trampert et al., 2004). Third, coupled preservation of primordial and recycled materials in the models agrees with the geochemical record of igneous rocks carrying the footprint of both (old) recycled crust (e.g., Chauvel et al., 1992; Hofmann, 1997; Cabral et al., 2013) and primordial mantle sources, such as ¹⁸²W $/$ ¹⁸⁴W and ¹⁴²Nd $/$ ¹⁴⁴Nd isotopic anomalies and/or high ³He $/$ ⁴He ratios (e.g., Jackson et al., 2010; Touboul et al., 2012; Rizo et al., 2016; Peters et al., 2018). The specific isotopic fingerprints of the different mantle reservoirs predicted by our models (e.g., primordial and recycled) and the subsequent sampling thereof by mantle plumes remain key topics for future research. Finally, the plate-like behavior displayed by the numerical models agree with various geological indicators of plate tectonics operating on Earth for at least several billions of years before present day (Korenaga, 2013, and references therein).

Out of the likely geodynamic styles (diffuse domains, marble cake mantle, and viscous blobs), the primordial blobs style (III.B) is most commonly predicted (Fig. 3). Model predictions for this sub-regime are also consistent with the stagnation of some slabs in a depth range that is similar to primordial-domain roofs, while other slabs readily sink into the deep mantle (Fukao and Obayashi, 2013), as seen in Figs. 2, 7, and videos in the Supplement (Gülcher, 2021). Furthermore, the presence of viscous domains in the mid-mantle may explain the observation of sharp seismic velocity contrasts in the uppermost lower mantle (850–1100 km depth), usually occurring away from major upwellings and downwellings (Jenkins et al., 2017; Waszek et al., 2018). Moreover, such bridgmanite-rich regions would suppress the average effects of the pressure-induced spin transition of iron in ferropericlase on seismic properties (Crowhurst et al., 2008; Murakami et al., 2012). Indeed, a recent study noted the absence of the expected signal of this spin transition for pyrolitic compositions at mid-mantle depths in global seismic profiles (Shephard et al., 2020), corroborating the presence of bridgmanite-rich regions in a depth range as predicted here. Finally, in terms of seismic tomography, the lack of clear evidence for primordial domains in the mid-mantle may be related to a trade-off between the seismic effects of their compositional (intrinsically fast bridgmanite, Wentzcovitch et al., 2004) and thermal (slightly warmer than the ambient mantle) anomalies.

Although various alternative hypotheses have been put forward to explain each of the observations above, a hybrid style of present-day Earth mantle heterogeneity as predicted here can provide a unified explanation for these various geochemical and geophysical observations. This hybrid style is illustrated in Fig. 9 and includes the coexistence of viscous primordial material and dense ROC material in a convecting mantle, with several primordial streaks and/or blobs forming a larger composite bridgmanitic structure (100 s–1000 km) around which convective flows are organized. Except for very large λ_prim, this prediction is somewhat different from that in Ballmer et al. (2017a), who predict large coherent blobs. ROC material is present as marble-cake streaks throughout the mantle with local enhancement in the mantle transition zone (MTZ) and at the CMB (i.e., forming thermochemical piles) (as in Nakagawa and Buffett, 2005; Yan et al., 2020). Piles in this regime are mostly of recycled origin but may include ancient FeO-enriched material that is preserved near their top if a FeO-rich basal layer ever existed in the early Earth. In this scenario for mantle composition and structure, the upper mantle is on average pyrolitic, while the lower mantle is significantly more enriched in silica, in line with previous studies (e.g., Murakami et al., 2012; Ballmer et al., 2015; Mashino et al., 2020; Yan et al., 2020).

4.3 Thermochemical pile layering

Thermochemical piles formed in our models comprise of ROC with varying formation ages. For models including an ancient FeO-rich basal layer, piles consist of both recycled and ancient materials (e.g., a “basal melange”, Tackley, 2012). Towards the roof of the piles, the formation age of ROC increases, while the overall ROC content decreases (Fig. 7). A radial distinction into regions of a high- and low-ROC content within piles was also found by Mulyukova et al. (2015). For a high buoyancy ratio of basaltic material, the thermochemical piles consist of a high-density basal layer covering nearly the entire CMB, overlain by high-topography piles with a much lower fraction of ROC (Mulyukova et al., 2015), as is the case for the thermochemical piles in our models. Moreover, in our experiments, ancient FeO-rich material survives for several billion years at the roof of some piles (Fig. 7). Compositional layering of ancient versus recycled domains within LLSVPs has been previously proposed (Ballmer et al., 2016; Trønnes et al., 2019), albeit with an inverted maturity gradient as found here.

Chemical layering strongly depends on the density profile and rheological properties of the material that makes up the thermochemical piles, which remain poorly constrained. We modeled the ancient FeO-rich basal layer to have the same physical properties as ROC. Final preservation of an ancient basal layer would be advanced, and the maturity gradient potentially inverted, if the layer were enriched in FeO or SiO₂ relative to ROC. Such an enrichment may or may not occur depending on the formation scenario of the ancient basal layer, e.g., the style of magma ocean crystallization (e.g., Labrosse et al., 2007; Ballmer et al., 2017b) or core–mantle interaction (e.g., Hirose et al., 2017; Trønnes et al., 2019). Moreover, ancient crustal rocks may have different chemistry and hence density at lower mantle depths than young basalts (e.g., Herzberg and Rudnick, 2012), which may also affect the layering sequence as seen in our models.

4.4 Linking recycled piles, primordial blobs, and mantle dynamics

The slight positive correlation between pile height and overlying primordial material (Fig. 8) differs from the previously conceptional model proposed in Ballmer et al. (2017b), in which the authors suggested that dense mantle material ideally piles up in zones of convergence atop the CMB, free from overlying viscous blobs (BEAMS). As we find piles to comprise of older material when overlain by primordial blobs (Fig. 8d, h), we infer that primordial domains can shield piles from any incoming subducted material and aid their longevity. The presence of primordial domains in the mid-mantle may therefore constrain the shapes and spatial fixity of LLSVP piles (Torsvik et al., 2010; Dziewonski et al., 2010).

Our results further contribute to the ongoing debate about whether thermochemical piles are intrinsically stable features in the deep mantle, which spatially determine mantle convective patterns (Dziewonski et al., 2010; Torsvik et al., 2014) or are instead pushed around by subduction zones (McNamara and Zhong, 2005; Zhang et al., 2010). In our models, piles and blobs both mostly reside well away from lower-mantle downwellings, as they are pushed away from them (videos in the Supplement, Gülcher, 2021). Indeed, as noted in previous studies (e.g., McNamara and Zhong, 2005; Zhang et al., 2010; Schierjott et al., 2020), we observe that downgoing slabs are mainly controlling the spatial distribution of piles (as well as viscous blobs) and not the other way around.

4.5 Early Earth shift in mantle dynamics, tectonics, and chemical sampling

An interesting prediction of our models involves the breakdown of ancient compositional layering (“overturn”), marking the beginning of whole-mantle convection (see Fig. 6, Supplement S3). The timing of the overturn varies between 0.5 and 2.5 Gyr of model evolution, mainly depending on the intrinsic strength of the primordial layer. This overturn triggers a shift in surface-tectonic style, as plate-tectonic behavior (P and M) robustly occurs in the models following the event. For Earth, a major change in geodynamic style is commonly proposed to mark the onset of “modern-style” plate tectonics (e.g., Bédard, 2018; Condie, 2018, and references therein), shifting from either an (episodic) stagnant lid (e.g., O'Neill et al., 2016; Bédard, 2018; Lenardic, 2018) or a sluggish lid with mostly vertical tectonics (e.g., Rozel et al., 2017; Foley, 2018; Lourenço et al., 2020) to a mobile-lid style with horizontal plate-like tectonics. Estimates for the timing of this event vary between studies (e.g., 2–3.2 Ga, Korenaga, 2013, and references therein). According to our model predictions, the breakdown of ancient compositional layering may be related to such a tectonic shift in early Earth. The related stirring of the mantle during the overturn may have drastically changed upper-mantle composition (e.g., Stein and Hofmann, 1994), which is likely followed by the oxygenation of the atmosphere, with implications for the evolution of life (Andrault et al., 2017).

Another implication of this geodynamic shift is that deep-rooted mantle plumes can reach the upper mantle for the first time. Accordingly, the geochemical signatures of mafic rocks of this age should carry stronger lower-mantle signatures compared to older mafic rocks. Many geochemical studies indicate a rapid, widespread change in isotopic signatures in basalts around ≈ 3 Ga (e.g., Gamal El Dien et al., 2020, and references therein). For example, μ¹⁸²W anomalies in basalts show a steady trend of mostly positive values for rocks of before ≈ 3 Ga but display negative values in the modern mantle. This change has previously been attributed to inner-core segregation (Rizo et al., 2019) or the onset of deep slab subduction and recycling of crustal material and sediments (e.g., Liu et al., 2016; Rizo et al., 2019). Here, we propose that it may (additionally) be related to the onset of primordial-material entrainment by mantle plumes via whole-mantle convection. The primordial material could have obtained a negative μ¹⁸²W signal through SiO₂ exsolution from the core, which has a strongly negative μ¹⁸²W signature (e.g., Kleine et al., 2009), to the (basal) magma ocean during rapid initial cooling of the planet (e.g., Helffrich et al., 2018b; Trønnes et al., 2019; Rizo et al., 2019).

It must be noted that the modeled thermal evolution of the mantle and related early convective vigor somewhat depends on our parameter choices. Although our models do not take into account heating from radioactive decay or core cooling (Nakagawa et al., 2010), the modeled mantles cool over time with internal potential temperatures from 1900 to ≈1600–1650 K, consistent with petrological studies on Earth's mantle temperatures over time (Herzberg and Rudnick, 2012). A prior study found that early model dynamics differ as a different initial mantle temperature or core temperature is applied, although the ultimate thermal and chemical structure after 4.5 Gyr does not depend much on these parameters (Nakagawa and Tackley, 2012). In particular, the scale of early mantle flow and the geometry of boundary layer instabilities are different when evolving core temperatures and internal heating sources are incorporated (Nakagawa and Tackley, 2012). With our constant core temperature of 4000 K, early Earth's CMB heat flow is likely underestimated and present-day heat flux slightly overestimated. The presence of internal heating from the decay of radioactive nuclides (or heat-producing elements, HPEs) increases convective vigor and decreases the length scales of mantle flow. On the other hand, it suppresses active hot instabilities and mantle plumes, which instead become diffuse and passive return flows (e.g., Mckenzie et al., 1974). Keep in mind that erosion of viscous blobs and the entrainment of primordial material is primarily accommodated by plumes. Nevertheless, considering uniform internal heating is likely to impede preservation, as blobs should heat up and rise into the upper mantle, where they would be efficiently processed. However, HPEs are unlikely to be uniformly distributed. In fact, they are unlikely to be incorporated into bridgmanitic material for a wide range of formation scenarios. For example, bridgmanitic magma ocean cumulates will not incorporate any significant levels of highly incompatible elements, including HPEs (Corgne et al., 2005; Brown et al., 2014). Instead, if most HPEs are partitioned into the dense, final MO cumulates (i.e., such as the ancient FeO-rich basal layer discussed in Sect. 3.3.2), they would be “trapped” in this insulating layer.

4.6 Model limitations and outlook

As in any geodynamic study, the numerical models remain a simplified approximation of nature. For example, our model resolution only allows us to simulate relatively large heterogeneities (on scales of several kilometers or longer). In addition, we do not include any felsic material in our models, which may contribute to lower mantle heterogeneity by recycling of continental material (Hofmann, 1997; Stracke, 2012). Furthermore, future investigations are needed to clarify how the heterogeneity styles operate in 3D geometry or in the presence of internal heating sources and core evolution.

The work presented here is primarily a numerical study that establishes various feasible regimes of mantle dynamics and heterogeneity mixing and preservation as a function of initial conditions and primordial-material properties. In the future, a thorough integration with inter-disciplinary observations on mantle heterogeneity is needed to further establish the applicability of some of our detected regimes for the Earth's (lower) mantle. For example, model results should be quantitatively compared with seismic constraints: piles are predicted to be anomalously hot and dense, which would leave a clear seismic signature (e.g., Wang et al., 2020). On the other hand, bridgmanitic blobs in our models are usually only slightly warmer than typical ambient mantle temperatures (Fig. 8). Therefore, the related seismic anomalies may be much less evident in tomography, particularly when considering that bridgmanite is slightly faster than ferropericlase (Tsuchiya et al., 2020). Indeed, seismic tomographic images do not show obvious large-scale seismic anomalies in the mid-mantle, and future studies should aim at establishing the presence or absence of these domains in the present-day Earth, e.g., focusing on seismic anisotropy or out-of-plane reflections.

5 Conclusions

We performed a numerical study to investigate the coexistence of different chemical heterogeneity with distinct physical properties in Earth's lower mantle. Our results demonstrate the following conclusions.

• Primordial and recycled heterogeneity may coexist in various styles in the mantles of Earth-like planets, depending on their rheological properties (intrinsic density and viscosity).

• The final volume of primordial heterogeneity preserved in the mantle only very weakly depends on the initial thickness of the primordial layer in the range of values explored here (35 vol %–65 vol % of the total mantle).

• The preservation of coherent viscous blobs in the mid-mantle is enhanced by the existence of an ancient FeO-rich basal layer underneath.

• The coupled existence of viscous, primordial blobs in the mid-mantle with dense, recycled (and partially ancient) piles in the lowermost mantle is a robustly predicted over a wide range of parameters.

• The presence of large viscous blobs in the mid-mantle stabilizes the underlying dense thermochemical piles, contributing to their preservation and longevity.

Our results provide a quantitative and dynamic framework for the coupled evolution of recycled and primordial materials in a convecting mantle. For planet Earth, we suggest that the lower mantle may be in a hybrid state between marble cake and plum pudding style of chemical heterogeneity. Such a hybrid mantle structure can reconcile geochemical evidence for ancient rock preservation in the convecting mantle through the present day with seismic evidence for whole-mantle convection. Finally, the breakdown of a stratified system with double-layered convection towards a whole-mantle convective system with plate-tectonic behavior as predicted by our models may be relevant for a major geodynamic shift during the Archean, as is indicated by the geological and geochemical record.