Narrow, Fast, and “Cool” Mantle Plumes Caused by Strain‐Weakening Rheology in Earth's Lower Mantle

The rheological properties of Earth's lower mantle are key for mantle dynamics and planetary evolution. The main rock‐forming minerals in the lower mantle are bridgmanite (Br) and smaller amounts of ferropericlase (Fp). Previous work has suggested that the large differences in viscosity between these minerals greatly affect the bulk rock rheology. The resulting effective rheology becomes highly strain‐dependent as weaker Fp minerals become elongated and eventually interconnected. This implies that strain localization may occur in Earth's lower mantle. So far, there have been no studies on global‐scale mantle convection in the presence of such strain‐weakening (SW) rheology. Here, we present 2D numerical models of thermo‐chemical convection in spherical annulus geometry including a new strain‐dependent rheology formulation for lower mantle materials, combining rheological weakening and healing terms. We find that SW rheology has several direct and indirect effects on mantle convection. The most notable direct effect is the changing dynamics of weakened plume channels as well as the formation of larger thermochemical piles at the base of the mantle. The weakened plume conduits act as lubrication channels in the mantle and exhibit a lower thermal anomaly. SW rheology also reduces the overall viscosity, notable in terms of increasing convective vigor and core‐mantle boundary heat flux. Finally, we put our results into context with existing hypotheses on the style of mantle convection and mixing. Most importantly, we suggest that the new kind of plume dynamics may explain the discrepancy between expected and observed thermal anomalies of deep‐seated mantle plumes on Earth.

To date, it is understood that mantle plumes can start as deep as the core-mantle boundary and rise all the way to the base of the lithosphere, where they sustain intraplate hotspot volcanism (Morgan, 1971). In the classical view (Howard, 1964;Richards et al., 1989), a rising mantle plume is characterized by a large head atop a narrow tail, although chemical complexities may result in deviations from such classical shapes (Davaille & Vatteville, 2005;Farnetani & Samuel, 2005;Lin & Van Keken, 2006). Plume shapes are difficult to clearly be imaged by seismic tomography due to wavefront healing (Ritsema et al., 2021). Although recent full-waveform tomography models hint at the presence of plume-like features associated with major hotspots (e.g., French & Romanowicz, 2015), ambiguity remains as to the vertical continuity of these features, as well as their shapes and stability (e.g., Davaille & Romanowicz, 2020;French & Romanowicz, 2015;Wolfe et al., 2009). Another controversy lies in the temperature excess of such mantle plumes. Based on petrologic thermo-barometry, excess temperatures of plumes have been estimated at 100-300 K (Boehler, 1996), with even lower estimates based on seismic observations (30-150 K, Bao et al., 2022). Both these estimated ranges are significantly lower than the expected CMB temperature difference of ∼1,000 K. While the dynamics and shapes of plumes are well studied in geodynamic models with Newtonian rheology, they strongly depend on the material properties of mantle rocks (e.g., Massmeyer et al., 2013). However, these properties, and in particular the rheology of the lower mantle, are ill-constrained.
The two main constituents of the lower mantle are bridgmanite (Br) and ferropericlase (Fp) (Hirose, Morard, et al., 2017). The viscosity of the strong mineral bridgmanite (η Br = 10 21 − 10 23 Pa⋅s) is several orders of magnitude larger then that of the weak ferropericlase (η Fp = 10 18 − 10 21 Pa⋅s) (Kaercher et al., 2016;Yamazaki & Karato, 2001). It has been suggested that for typical mantle rocks that comprise Br as well as Fp, the bulk viscosity of the rock decreases with ongoing deformation as the weaker ferropericlase crystals become elongated in the direction of strain and interconnect with each other (Girard et al., 2016) (see Figure 1a). This experimental result was further confirmed in numerical studies on the effective rheology of a lower-mantle two-phase medium during deformation (e.g., de Montserrat et al., 2021;Thielmann et al., 2020). The weakening of the bulk rock with accumulating strain ("strain weakening") implies that deformation may localize in the lower mantle, analogous to localized shear zones in crustal rocks. Such strain localization may potentially explain the isolation of large unmixed domains in the lower mantle, which may host primordial (or "hidden") geochemical reservoirs away from regions of localized deformation (Ballmer et al., 2017;Chen, 2016;Gülcher et al., 2020;Mundl et al., 2017) (see Figure 1b). However, the effects of strain weakening on lower-mantle convection patterns and mixing dynamics have not yet been studied using global-scale geodynamic models.
Here, we implement a macro-scale description of strain-weakening (SW) rheology in a global mantle convection model. We present 2D numerical models of thermochemical convection in a spherical annulus geometry that include a new implementation of tracking the strain ellipse at each tracer through time. We allow lower mantle materials to rheologically weaken to various degrees and investigate the effects of this rheological weakening on mantle convection dynamics. We particularly focus on the characteristics of mantle plumes in the models. We find that SW rheology has several effects on mantle dynamics, including on the (a) pattern of mantle flow, (b) thermal evolution of the mantle, (c) pile stability, and (d) mantle plume dynamics. We distinguish first-order effects (directly caused by spatial viscosity variations resulting from SW rheology) and second-order effects (indirectly caused by a changing Rayleigh number of the ambient mantle caused by SW), and link the results to the previously-proposed style of mantle convection (Figure 1b). The changing plume dynamics are of particular interest since weakened plumes could explain the discrepancy between expected and observed thermal anomalies of deep-seated mantle plumes on Earth.

Numerical Technique and Model Set-Up
In this study, we use the finite-volume code StagYY (Tackley, 2008) to model mantle convection in two-dimensional spherical annulus geometry (Hernlund & Tackley, 2008) for 5 Gyrs of model time. The conservation equations for mass, momentum, energy and composition are solved on a staggered grid for a compressible fluid with an infinite Prandtl number (find full description of these equations in e.g., Tackley, 2008). In the energy conservation, we include contributions of adiabatic, latent, and shear heating, but ignore radiogenic heating (see below). The computational domain is discretized by 1,024 × 128 cells, in which ∼2.5 million tracers, tracking composition, temperature, and strain, are advected (20 tracers per cell). Due to vertical grid refinement near the boundary layers and near 660 km depth, the size of grid cells varies between 4 and 25 km in the vertical direction. Free-slip and isothermal boundary conditions are employed at the top and bottom boundaries, with a fixed surface temperature of 300 K and CMB temperature of 4,000 K. The numerical experiments are purely bottom heated (i.e., no internal heating and constant core temperatures).
Initial mantle temperatures are calculated from an adiabat with a potential temperature of 1900 K, together with the top and bottom boundary layers, and superimposing small random temperature perturbations of ±10 K on the cell level. The resulting temperature profile is ∼300 K warmer than Earth's present-day geotherm. Such increased initial mantle temperatures are applied to crudely mimic the thermal evolution of the mantle (especially since our models are purely bottom-heated), with higher convective vigor and widespread near-surface melting in the early model evolution (e.g., Davies, 2007;Herzberg et al., 2010). Most modeled mantles in this study cool down to ∼1,600-1650 K, consistent with petrological studies on Earth's mantle temperatures over time (Herzberg & Rudnick, 2012).

Treatment of Mantle Composition, Phase Changes and Melting
The driving forces of mantle convection are related to rock density, which depends on temperature and composition (Nakagawa & Buffett, 2005;Tackley et al., 2013). Composition in our modeled mantle has two lithological end-member components: harzburgite and basalt, which are each treated as a mixture of olivine and pyroxenegarnet systems that undergo different solid-solid phase transitions (for details, see Nakagawa et al., 2010). Harzburgite is considered to be a mixture of 75% olivine and 25% pyroxene-garnet; basalt is considered as pure pyroxene-garnet (see their phase change parameters in Table 1). Each tracer carries a mechanical mixture of harzburgite and basalt, initially set to 80% harzburgite and 20% basalt (i.e., pyrolitic composition). The density profiles of harzburgite and mid-ocean ridge basalt are consistent with those from Xu et al. (2008), and they are plotted in Figure S3 of Supporting Information S1. Compositional anomalies carried on tracers evolve from the initial state due to melt-induced differentiation: they undergo partial melting as a function of pressure, temperature and composition to sustain the formation of basaltic crust, and leaving a harzburgitic residue (for details, see Nakagawa et al., 2010). fraction for the two end-member textures of "load-bearing framework" (LBF) and non-linear "interconnected weak layers" (IWL) (Handy, 1994). Image adapted from Ballmer et al. (2017). (b) Suggested mantle convection dynamics in which shear localization of weak Fp grains induces weak layers of IWL along slabs and plumes, and mixing is less efficient for the bridgmanitic LBF part of the lower mantle, potentially promoting the preservation of long-lived geochemical reservoirs. Reprinted from Chen (2016), with permission from The American Association for the Advancement of Science (AAAS).

Visco-Plastic Rheology
Viscosity is temperature-, pressure-, composition-(or phase-), and strain-dependent following an Arrhenius-type viscosity law (Newtonian rheology): where η 0 is the reference viscosity at zero pressure and reference temperature T 0 (=1,600 K), E a and V a are the activation energy and volume, respectively, T is the absolute temperature, P the pressure, and R is the gas constant (8.314 J⋅mol −1 K −1 ). Composition-dependency is considered through prefactor λ c : a viscosity increase of one order of magnitude is imposed along the 660-km depth boundary (consistent with e.g., Čížková et al., 2012), and a viscosity decrease of 10 −3 relative to the lower mantle is imposed at post-perovskite phase transition in the lowermost mantle (as suggested by Ammann et al., 2010). By considering a Newtonian rheology, the implicit dominant deformation mechanism is diffusion creep. Note, however, that we consider an activation energy that is smaller than experimental constraints for the upper mantle, an approach that can account for the effects of dislocation creep (e.g., Christensen & Hofmann, 1994;van Hunen et al., 2005). With our chosen activation energy and volume, the activation enthalpy in the lower mantle varies from 262 kJ/mol to 548 kJ/mol, in line with perovskite predictions (Yamazaki & Karato, 2001). An additional strain-dependency of viscosity is implemented through the weakening factor f w (see next section). All physical and rheological parameters used in this study are listed in Table 2.
In order to obtain plate-like behavior at the surface, we assume that the material deforms plastically when a critical pressure-dependent yield stress is reached (as in Crameri & Tackley, 2014;Tackley, 2000a) (see Table 2). Plate-like behavior is evaluated using diagnostic criteria from Note. The table shows the depth and temperature at which a phase transition occurs; △ pc and denote the density jump across the phase transition and the Clapeyron slope, respectively. The Clapeyron slopes for these phase changes are similar to those used in previous numerical studies (e.g., Tackley et al., 2013). In the olivine system, the 410 and 660 phase changes are made discontinuous, whereas all other phase changes in all systems are defined as hyperbolic tangent functions (tanh) that transition between the phases across a predefined phase loop width. Finally, 0 refers to the reference bulk modulus for the system for each individual layer (marked by the depth range), and ′ 0 refers to its pressure-derivative. The reference profiles of density, thermal expansivity, thermal conductivity, and specific heat of modeled mantle materials are given in Figure S3 of Supporting Information S1. Note. PV = perovskite; PPV = post-perovskite. Since we solve for compressible convection, the adiabatic temperature, density, thermal conductivity, and heat expansivity are pressure-dependent following a thirdorder Birch-Murnaghan equation of state . The reference profiles of density, thermal expansivity, and thermal conductivity of modeled mantle materials are given in Figure S3 of Supporting Information S1.  Tackley (2000b): plateness p (the degree to which surface deformation is localized) and mobility m (the extent to which the lithosphere is able to move). An ideal plate tectonic style gives p = 1 and m close to or larger than 1.

Finite Strain
Our finite deformation approach follows that of McKenzie (1979) and builds up on previous work (Xie & Tackley, 2004). The deformation tensor M is tracked on the advecting tracers. At each timestep, the velocity gradient tensor is calculated and interpolated to the position of each tracer, from which the additional deformation for that timestep can be retrieved. This additional deformation is then added to the existing deformation tensor. A second-order accurate in time approach is implemented following (equations 7-10 in McKenzie, 1979). This approach tracks stretching, rotation and advection in full tensor form ( Figure 2a). The eigenvectors of the matrix MM T give the principal directions of the strain ellipse, and the square root of the eigenvalues of this matrix give the amount of strain in each direction (i.e., the semi-major and semi-minor axes a and b). This implementation has been tested and analyzed in Text S1 of Supporting Information S1. The scalar (finite) strain ϵ is calculated from the semi-major and semi-minor axes of the strain ellipse (see Text S2 in Supporting Information S1 for more detail) and is defined as: The deformation matrix M is reset to a unit matrix (i.e., the strain ellipse is reset to a circle) when a tracer undergoes melting-related differentiation or passes the depth of 660 km. This represents the resetting of the micro-structure of the rock when it melts or passes through the major 660 solid-solid phase transition (e.g., Solomatov & Reese, 2008).
We have chosen to track deformation and apply weakening/healing (see below) in full tensor form, because this allows for analysis of full history-dependent deformation. This differs from other studies that only track the second invariant of strain rate (e.g., Fuchs & Becker, 2019;Tackley, 2000b). Tracking the deformation matrix (e.g., strain ellipse) also has many potentials for future research, including directional information of deformation related to mineral fabric, lattice-preferred orientation, and seismic anisotropy. (the maximum weakening factor f w ), the critical strain threshold ϵ th (the strain at which half of the maximum weakening has taken place after which weakening occurs), and the shape factor C. A comparison of these weakening curves with other studies is presented in Figure  S4b of Supporting Information S1 and discussed in Text S2 of Supporting Information S1.

Strain Weakening Parametrization
Lower-mantle material is rheologically weakened following a simplified strain-dependent weakening curve ( Figure 2b): where ϵ crit is the critical strain threshold parameter (the strain at which half of the weakening has taken place, see Figure 2b), max w is the maximum weakening factor (1 for no weakening), and C (=2) controls the shape of strain weakening curve. Equation 3 was chosen as such to mimic a smoothly evolving weakening curve that can be easily parameterized (i.e., tanh function). Weakening occurs almost instantaneously as material becomes deformed (in agreement with Girard et al., 2016;Thielmann et al., 2020;de Montserrat et al., 2021) and the general trend is in line with predictions from micro-scale numerical studies (with Thielmann et al., 2020;de Montserrat et al., 2021), see Text S2 in Supporting Information S1 for a more elaborated discussion. The maximum amount of weakening ( max w ) is a free parameter. For simplicity, we neglect the composition-dependency of strain weakening, as well as the anisotropy of viscosity according to the strain tensor, to establish the first-order effects of SW rheology on the lower mantle. Neglecting anisotropy assumes that the strain ellipse is well-aligned with the dominant shear direction at any finite strain.

Rheological Healing
Processes such as diffusion-dominated annealing and/or grain growth may lead to the relaxation of deformed grains (e.g., Solomatov & Reese, 2008). For part of the models, we approximate such rheological healing by relaxing the deformation matrix (strain ellipse) toward a unit matrix (circle) with a temperature-dependent (grain growth highly depends on temperature) and pressure-dependent (diffusion is limited at higher pressures) term. This term may represent diffusion creep-dominated healing of micro-structures. Healing lowers the local strain ϵ and subsequently reduces the effective weakening. The temperature-dependency of rheological healing allows for long-term memory of deformation at cold temperatures and faster healing at deep mantle temperatures, whereas the pressure-dependency may counteract fast rheological healing at high T-P regions. Rheological healing is implemented as: where H is the rheological healing rate (H = 1 h with reference timescale of rheological healing τ h ). By rewriting Equation 4, M is updated as follows: where M new is the updated deformation matrix, M old is the deformation matrix before healing, Δt is the timestep in seconds (=t new − t old ), and I is the unit matrix. This healing rate causes M to relax toward a unit matrix I. This relaxation produces a slight rotation of the strain ellipse, but this does not affect our model results since we assume isotropic properties. Text S1 in Supporting Information S1 shows a more detailed analysis of this healing implementation and an outlook on how to keep the orientation equal. The rheological healing rate H is temperature-and pressure-dependent following: where H 660 is the reference healing rate (s −1 ) at the top of the lower mantle (P = P 660 , T = T 660 , along the reference adiabat); this depth is chosen because weakening and healing only occur in the lower mantle. The activation energy and volume are assumed to be the same as those for diffusion creep because atomic (vacancy-) diffusion is the mechanism by which healing occurs. More information on this rheological healing rate, the chosen reference values, and how this parametrization ultimately affects the distribution of strain in the models is described in Text S3 of Supporting Information S1.

Automated Detection of Mantle Domains
We use the geodynamic diagnostics toolbox StagLab (Crameri, 2018) to automatically detect regional flows that are either self-driven (i.e., active) or induced (i.e., passive). Active regional flows represent mantle plumes (active upwellings) or active slabs (active downwellings). Passive slab remnants in the mantle that are not actively sinking, are not detected. Detecting mantle plumes using a purely thermal definition (e.g., Labrosse, 2002) or a purely dynamic definition (e.g., Hassan et al., 2015) does not work for our set of numerical models, as our modeled plumes significantly vary in terms of their anomalies in both temperature and radial velocity. We use a new combined thermo-dynamical approach of identifying mantle plumes and slabs. As a summary, active plumes and slabs are defined based on a combination of the horizontal residual fields of temperature and radial velocity: T res ⋅ |v r, res |. The residual fields T res and T res ⋅ |v r, res | at each timestep are statistically analyzed, and their percentiles are used for the definition of the plume/slab thresholds. More details on this approach is given in Text S4 of Supporting Information S1. Finally, the mantle diagnostics routine also detects thermochemical piles present atop the CMB based on composition and temperature. In terms of composition, pile material must consist of at least 60% basal (C bs ≥ 0.6), while the temperature constraint is defined using the average of a mid-mantle temperature of 3000 K and the CMB temperature: ≥ (3000 + CMB∕2) (K) (as in Schierjott et al., 2020). Using this routine, the physical properties within all different mantle domains can be separately explored.

Parameter Study
In this numerical study, we vary several parameters. Most of them are related to the newly-implemented SW parameterization: we systematically vary the maximum strain weakening factor max w from 1 (no weakening) to 0.01. The range of explored weakening factor (1-0.01) implies a viscosity contrast range of 1-100 between the "interconnected weak layers" (IWL) and "load-bearing framework" (LBF) viscosity mechanisms. For a ∼pyrolitic mantle, this would coincide with a viscosity contrast between bridgmanite and ferropericlase ranging from 1 to ∼200-300, which is consistent with mineral physics estimates (e.g., Kaercher et al., 2016;Yamazaki & Karato, 2001). The reference healing rate H 660 is varied in the range of 10 −14 -10 −16 s −1 , which corresponds to 1 healing time scales of 3-300 Myr for the uppermost lower mantle. The fastest healing time of 3 Myr is in line with studies on single-mineral chemical diffusivity (using e.g., diffusivities of 10 −18 m 2 s −1 and d = 1 cm in Ammann et al. (2010)). We expect such self-diffusion timescales to give lower bounds in terms of healing timescales for a complete mineral fabric. Finally, since these models including SW rheology show different final viscosity profiles, we select certain SW cases which we run again with an increased viscosity contrast between the upper and lower mantle (λ 660 > 10). The combined effect of such increased intrinsic lower-mantle viscosity and strain weakening behavior (lowering the lower mantle viscosity) gives more similar final viscosity profiles. As such, the effect of a changing Rayleigh number on mantle dynamics is eliminated and the direct effect of SW rheology can be determined.

Results
The relevant model parameters and output variables of all models run in this study are summarized in Tables S1 and S2 (Supporting Information S1). Videos related to the cases discussed in the text and/or figures can be found in the Supporting Information S1 related to this article. All models show plate-like behavior according to the diagnostics of Tackley (2000b), and final average viscosity profiles of most models approximately agree with estimates from literature (see Figure S4 in Supporting Information S1). We first describe the evolution of our reference model (Section 3.1). In this model, the strain field is tracked according to the newly-implemented finite strain approach without the application of strain-dependent weakening or rheological healing. In Section 3.2, we 8 of 21 describe the effect of the implemented SW rheology on model behavior. Two case studies are highlighted with various degrees of SW and rheological healing. Finally, in Section 3.3, we summarize the results of several case studies with various degrees of SW, but a similar final viscosity profile (Section 3.3).

Reference Model Evolution
The temporal evolution of our reference model (M 0 ), which does not include SW, is shown in Figure 3 and Video S1. Soon after the model onset, the thermal boundary layers grow in amplitude, and after ample growth of boundary layer instabilities, a mantle overturn initiates the onset of whole-mantle convection and plate-tectonic behavior (Table S1 in Supporting Information S1). The viscous flow associated with early model dynamics (at 1.0 Gyr) causes a lower-mantle strain field that is localized in regions of buoyant, hot upwellings, and areas which are deflected by incoming, strong lithospheric drips/slabs (Figure 3a). From the start of whole-mantle convection, the mantle gradually cools and the frequent occurrence of active mantle plumes and subducting slabs and active mantle plumes causes further complexity of the mantle strain pattern (Figure 3). As the deformation history is reset at the 660 km boundary layer, small-scale strain patterns in the upper mantle are not carried into the lower mantle. Instead, strain builds up in downwelling (e.g., around slabs) and upwelling regions (plumes) of the lower mantle ( Figure 3).
The radially-averaged temperature profile displays the typical signal of efficient whole-mantle convection with boundary layers superimposed on a mostly adiabatic geotherm (Figure 4a). The radial viscosity profile reflects the temperature-and depth-dependent rheology, as well as its compositional dependency expressed as a viscosity step toward higher values from the upper-to lower mantle (λ 660 ) (Figure 4b). Mantle velocity is highest in the bottom ∼150 km of the lower mantle and in the upper mantle ( Figure 4c). Finally, the compositional profile shows efficient basalt segregation in a thin region on top of the CMB and in the mantle-transition zone (Figure 4d).

Influence of Strain-Dependent Rheology
In this section, we separately describe the effects of SW rheology on several key aspects of mantle convection which were introduced in Section 1. First, the effect of SW rheology on convective flow patterns is described (Section 3.2.1), followed by its effect on the thermal evolution of the mantle (Section 3.2.2), on thermochemical piles formation (Section 3.2.3), and on the dynamics of mantle plumes (Section 3.2.4).

Global Mantle Convective Patterns
The radial profiles of viscosity and root mean-square (RMS) velocity of modeled mantles (Figures 4b and 4c), as well as their final averages (Figures 5b and 5c) show clear trends for SW rheology models. The average mantle viscosity is significantly lowered when SW rheology is applied, mostly in the lower mantle ( Figure 4b). Final convective vigor (∼v RMS ) is increased for most SW models, also mainly accommodated in the lower (most) mantle ( Figure 4c). Figure 6 shows the detected mantle domain field, that is, passive/active up -and downwellings, of three selected models with variable degrees of SW rheology at ∼4.5 Gyr. It further shows the temporal evolution of the lateral distribution of this mantle domain field at 1,800 km depth (i.e., in the middle of the lower mantle). In comparison to the reference model (Figure 6a), SW models with an increased convective vigor show a more chaotic planform of mantle flow (Figures 6b and 6c) with a larger number of small plumes present. The timescale of convection decreases with SW rheology as the convective vigor increases (overturn time ∼ 1 RMS , Miyagoshi et al. (2017)). The length-scale of convection also decreases with SW rheology as the lower mantle consists of convection cells with smaller aspect ratios, that is, more narrow regions of up-and downwellings (Figures 6b and 6c). The temporal evolution shows that the upwelling plumes in the lower mantle have shorter lifetimes than those in non-SW models. Figures 7a-7c show the distribution of selected quantities in the whole mantle domain for the same three models. The v RMS histograms for SW models are more skewed than that for the reference model, highlighting small domains in the mantle with very high velocities. This highlights the (albeit small) domains in which SW efficiently occurs. Interestingly, despite the changing pattern and vigor of mantle flow (described above), the statistical distribution of the age of mantle materials (defined by the time since a tracer last underwent a melting episode) in the whole mantle is similar for cases with and without SW (right panels of Figure 7).

Thermal Evolution
In models with SW rheology, the mantle cools down to higher final equilibrium mantle temperatures than in models with less or no SW rheology (Figures 4c and 5a). Even with a high healing rate of H 660 = 10 −14 s −1 , which causes strain in the lower mantle to heal on short geological timescales (Text S3 in Supporting Information S1), higher final internal temperatures are reached. This is also apparent by the relation between top and bottom Nusselt number (Nu top and Nu bot ). Nu top is not much affected by SW rheology, while Nu bot is on average significantly higher for SW models (Figure 5d). Likely, Nu top is not much affected since convective stresses at the base of the lithosphere are roughly similar for all cases due to similar upper-mantle viscosities (see Figure 4b). A higher Nuy bot implies that heat is more efficiently removed from the core. In our models (constant core temperature), this causes a higher final mantle temperature in models in which SW rheology is applied. How our assumption of constant core temperature may affect these results is discussed in Section 4.4.

Formation of Piles
Basalt segregation is most efficient in models with strong SW, because these models exhibit a low-viscosity, high-temperature mantle (see above and Figure 4d, and Videos). Due to efficient segregation in these models, basaltic material tends to settle more efficiently in the lowermost mantle, stabilizing the formation of large thermochemical piles that cover more of the CMB (Figure 6; Table S1 in Supporting Information S1). Also, the transition zone (410-660 km depth) becomes more enhanced in basaltic material (f bs ≈ 40%), while the regions around it are more depleted (Figure 4d) (as in e.g., Davies, 2008;Nakagawa & Buffett, 2005;Ogawa, 2003;Yan et al., 2020). Moreover, the thermochemical piles are internally convecting in SW models, affecting overall heat fluxes through the mantle: due to the high intrinsic pile density, heat can build up within the piles, forming another thermal boundary layer on top of them. Thus, these piles can act as plume generation zones in the models as active upwellings sporadically form on top of them, as well as along their edges.

Plume Dynamics
Strain mainly localizes in hot upwellings, where material is deformed soon after model onset and subsequently subjected to weakening. Even in the presence of rheological healing, plume channels are still weakened due to The horizontal axis represents the reference healing rate H 660 , and the icon shape stands for the implemented SW factor max w . The error bars and colored shaded areas indicate the standard-deviation of the parameter over that time period. For selected cases, outline-only symbols are also plotted, which represent the results for the additional cases with an increased viscosity jump at the 660 km discontinuity (see Section 3.3).

Figure 6.
strain localisation (see Text S3 and Figure S3 in Supporting Information S1). According to classical thermal plume formation theories (e.g., Howard, 1964), the temperature and width of a mantle plume is related to the time over which the boundary layer grows. A thermal instability occurs when the critical boundary layer Rayleigh number is reached (Ra local = Ra local,crit ) (Howard, 1964). If the viscosity of the material is large (low Ra local ), longer onset times occur, that is, the growth rate of the instability depends inversely on the local viscosity (e.g., Griffiths & Campbell, 1990;Howard, 1964;Olson, 1990). In our SW models, the lower viscosity of weakened plume materials would, according to the theory above, decrease the onset times of the instabilities. In these SW thermal instabilities, there is also a smaller temperature build-up compared to non-SW thermal instabilities. This precludes the growth of a wider, anomalously hot mantle plume with mushroom-shaped heads as seen in our reference model (Figure 6a). This absence of mushroom-shaped plume heads is evident for weakened plumes in both pure SW cases as well as SW + rheological healing models (Figures 6b and 6c; Videos S2-S3).
The positive feedback between weakening and strain localization causes a low-viscosity channel to form in and around plumes, allowing for rapid transport of mass and heat from depth toward the surface. The typical velocities in the plume conduit increase for more efficient SW (i.e., for lower max w or lower H 660 ), while the excess  temperatures are lower for more efficient SW (Figure 7, Table S2 in Supporting Information S1). These distinct plume dynamics caused by SW rheology are further apparent in the bottom Nusselt number (Nu bot ). Final Nu bot is, on average, higher for models with most efficient SW rheology, linked to the thinner thermal boundary layers and higher boundary layer Rayleigh numbers (Figure 5f). This prediction implies that heat is more effectively lost by convection (i.e., via mantle plumes) rather than conduction. The weakened conduits are easily deflected by background flow or by incoming slabs (see Videos S2 and S3). Typical timescales of plume lifetimes decrease from 500∼1,000 Myr for the reference case to few 100 Myrs for the extreme SW cases (lower panels in Figure 6).

Influence of the Mantle Viscosity Profile
Each model discussed above displays a distinct effective viscosity profile through time (Figure 4), which, in turn, controls convective vigor and thereby strongly affects model evolution. In order to distinguish the direct (first-order) effects of SW rheology on mantle dynamics from the indirect effects of SW rheology through the radial viscosity profile (second-order), we explore five additional SW cases with a higher intrinsic viscosity jump at 660 km depth (λ 660 ), such that the final viscosity profile is similar to that of the reference case M 0 , and to each other (see Figure 8b). Videos S4 and S5 show the time evolution of internal dynamics for two selected cases (those illustrated in Figures 6d, 6e and 9). These additional cases show a similar average thermal evolution as the reference case (Figure 8a). Despite the similarities in mantle viscosity, the convective vigor is still affected by SW due to localization of flow. The average values for v RMS are ∼50% higher in the additional SW models than in the reference model case (Table  S2 in Supporting Information S1; Figure 9b). In particular, the v RMS histogram is more skewed, with high values (>2 cm/yr) for SW models, and hardly any such very high velocity domains for the reference case (Figure 9b). Most of this faster flow is focused in narrow, weakened upwelling regions (Figures 9b and 9e and Videos S5). The higher convective vigor, in turn, affects the pattern of mantle flow (Figures 6d-6e), but much less so than the corresponding SW rheology models with a fixed λ 660 = 10 (Figures 6b and 6c). Despite more similarity of the length-and timescales of convection for the additional SW cases and the reference case, SW rheology still causes the formation of narrower convection cells and more mantle plumes with shorter lifetimes (Figures 9b and 9c). In terms of the distribution of the age of all mantle materials (Figure 9c), the mean age is indistinguishable between cases. Further, top Nusselt numbers are similar between cases. In turn, bottom Nusselt numbers are still slightly increased in the additional SW rheology models (Figures 5d and 5e), but much less than in the related SW rheology models with a fixed λ 660 = 10. These results strengthen our conclusion that CMB heat flux is preferentially accommodated via convection versus conduction, particularly in models with SW. However, plume dynamics and the size of thermochemical piles are still affected in the same way as in the previously described SW models. In fact, the localization of increased flow velocity in the narrow upwelling mantle plumes is even significantly more pronounced in these additional SW models. Moreover, thermochemical piles in the additional SW models are still substantially larger. Hence, we conclude that SW rheology is the critical ingredient for the weakening of plume channels, their narrow shapes and relatively low thermal anomalies, as well as the formation and stabilization of large thermochemical piles. Second-order effects, such as the higher final mantle temperature, and the significantly higher average mantle flow velocities, are caused by the modification of the viscosity profile through SW.
Based on the outcomes above, we conclude that SW rheology is the critical ingredient for (i) weak and (ii) narrow plumes with (iii) relatively low thermal anomalies, as well as (iv) large and long-lived thermochemical piles. Secondary effects of SW rheology, such as on CMB heat flux and mantle cooling, are caused by its effects on mantle dynamics by increasing convective vigor (locally). At a fixed average lower-mantle viscosity, CMB heat flux and mantle thermal evolution is still affected by SW, but much less so than in the models discussed in Section 3.2. Figure 9. Histograms of selected quantities for the reference model M 0 (black) and two additional cases with variable combinations of strain-weakening and healing plus an increased viscosity jump (λ 660 ) at the 660 km discontinuity (orange, blue). The color scheme corresponds to that in Figure 6, with the difference being the higher viscosity jump in the additional cases explored here. See Section 3.3 for details. (a-c) Distribution of the strain field, the root mean-square (RMS) velocity, and the material age within the whole mantle domain. (d-f) Distribution of the strain field, the RMS velocity, and the average horizontal temperature anomaly of the material within the detected active mantle upwellings (plumes). The vertical lines represent the mean values for each histogram. Note that the vertical axis of the strain histogram in panel (a) is logarithmic, whereas that in (d) is not.

Mantle Mixing and Geochemical Reservoirs
With an increasing convective vigor and decreasing length-scale of convection for SW models (Figure 6), one might expect mixing efficiency of mantle material to increase (e.g., Coltice & Schmalzl, 2006). However, in high convective-vigor SW models, basalt more easily segregates from harzburgite and thermochemical piles, which are in turn more stable over time (Section 3.2.3). Such a relation between lower mantle viscosity and more efficient basalt segregation is consistent with other studies (e.g., Yan et al., 2020). Hence, heterogeneity mixing turns out to be less efficient in models with SW rheology. Yet, the statistical distribution of mantle material age in the whole mantle is similar for all cases (Figures 7c and 9c). While a slightly higher proportion of very ancient material (>4 Ga) is preserved in SW models, a significant part of this material portion is accommodated in the larger thermochemical piles. The similarity of preservation in all our models, and particularly in the convecting mantle, is contrary to earlier suggestions that SW promotes the survival of primordial materials (e.g., Chen, 2016;Girard et al., 2016). In our SW models, convection patterns are not critically stabilized over time. This result may be attributed to the lack of effective strain weakening in the low-strain downwelling regions. Only if both upwellings and downwellings were significantly weaker than the regions in-between, efficient preservation may occur in these in-between regions (Ballmer et al., 2017;Gülcher et al., 2020). For example, it has been proposed that grain-size reduction in cold slabs that enter the lower mantle causes local weakening (Dannberg et al., 2017;Ito & Sato, 1991;Karato et al., 2001;Yamazaki et al., 2005). Such grain-reduction weakening in combination with SW plumes may cause a style of convection dynamics more akin to previously proposed (Figure 1b), with weakening occurring in both downwelling slabs and upwelling plumes. Future work should test if this is indeed the case.

Planetary Interior Evolution
Since SW rheology in the lower mantle affects CMB heat fluxes and their ratio to surface heat fluxes ( Figure 5, Table S1 in Supporting Information S1), it may have a substantial control on core dynamics as well as mantle cooling rates. The heat transfer from the core into the base of the mantle greatly affects the sustainability of a planetary dynamo through its control on the vigor of core convection, and the onset time of inner core crystallization (Lay et al., 2008;Stevenson, 2003). Moreover, the spatial pattern of (geo-)magnetic secular variations is commonly attributed to changes in CMB heat fluxes and mantle plumes (e.g., Biggin et al., 2012;Courtillot & Olson, 2007;Larson & Olson, 1991).
Modern estimates of CMB heat flux for the Earth range from several TW up to 15 TW (e.g., Frost et al., 2022;Lay et al., 2008;Nakagawa, 2020), that is, significantly lower than outcomes in our numerical models (Table  S1 in Supporting Information S1). However, these estimates may be underestimated since they do not consider additional CMB heat flux by advection due to cool plumes (subducted slabs) arriving at the base of the mantle (Labrosse, 2002). In turn, our models neglect internal heating (as they are not tuned to exactly match Earth thermal evolution), and thus overestimate plume activity and strength (e.g., Mckenzie et al., 1974). They also ignore core cooling (constant core temperature of 4,000 K), and hence early and present-day CMB heat fluxes are likely underestimated and overestimated (Labrosse, 2003), respectively. Nevertheless, in our models, mantles with SW rheology pull out heat more efficiently from the core, which would alter planetary thermal evolution. If core cooling were combined with SW rheology, the core would likely tend to cool faster which, in turn, would lower CMB heat flux and possibly final mantle temperatures. At a given present-day geotherm, this may imply hotter early-core and early-mantle temperatures, and/or more early mantle melting. Future studies should investigate the combined effect of SW rheology and core cooling, and assess whether the indirect effect of SW on mantle temperatures (i.e., only minor when comparing models with similar viscosity structures, see Section 3.3) still holds.
The relevance of SW rheology for (exo-)planets depends on their mineralogy and internal structure. Stars in the solar neighborhood show diverse Mg-Fe-Si compositions (Hinkel et al., 2014). Assuming stellar compositions as a proxy for rocky planet compositions (as in Spaargaren et al., 2020), planets in stellar systems with Mg/ Si < 1 likely have no ferropericlase in their mantle, hence no strain weakening is expected to occur. Rocky planets associated with Mg/Si ≫ 1.5 stars feature significant amounts of ferropericlase in their mantles, hence mantle rocks would form a IWL at very small (or even zero) strain. The majority of rocky exoplanets should have an Earth-like bulk composition with 1 < Mg/Si ≪ 1.5 (e.g., Putirka & Rarick, 2019;Spaargaren et al., 2020), where SW rheology potentially occurs in the mantle. Moreover, a recent study established the stability of a very weak B2-(Mg,Fe)O phase under extreme pressures (Coppari et al., 2021), which may dramatically affect the deep-mantle rheology of Super Earths.

Thermochemical Piles
The current degree-2 pattern of Earth's mantle flow, anchored by the two antipodal LLSVP piles (Dziewonski et al., 2010), has been suggested to be a stable energy configuration from the point of view of Earth's moment of inertia and to exist for at least 200 Myrs (Burke et al., 2008;Conrad et al., 2013;Torsvik et al., 2010Torsvik et al., , 2014). Yet, this configuration is only energetically stable if LLSVPs are (much) denser than slabs, which remains unclear (Koelemeijer et al., 2017a(Koelemeijer et al., , 2017bMcNamara, 2019). The intrinsic high density of recycled crustal materials (basalt) causes piles to form in our models (Figure 6), in agreement with various geodynamical studies (Li et al., 2014;Nakagawa & Buffett, 2005;Nakagawa et al., 2010;Tackley, 2012). Here, we show how SW rheology causes more efficient basalt segregation, and the formation of piles that cover a larger extend of the CMB (see Figure 6 and Table S1 in Supporting Information S1). Such thermochemical piles can act as a thermal insulator of part of the heat coming from the CMB (e.g., Lay et al., 2008;Nakagawa, 2020). Yet, the overall CMB heat flux is increased in SW rheology models. These increased values are accommodated by the much larger fluxes within the weakened plume channels (Table S2 in Supporting Information S1) as well as small-scale convective fluxes within the piles (e.g., Figure 6b). Plumes formed from this (secondary) thermal boundary layer have by default less heat available (Farnetani, 1997), since the temperature difference between the piles and ambient mantle is less than that between the adiabat and the CMB. This could partially explain lower temperature anomalies of these mantle plumes, although note that weakened plumes not only rise from the piles, but also from the CMB (see Figure 6; Figure S6 in Supporting Information S1, and Videos S2-S5).

Plume Formation
As described in the Section 3.2.4, the differences in plume dynamics between our numerical models agree well with scalings and relationships found in early classical thermal plume formation theories (e.g., Griffiths & Campbell, 1990;Howard, 1964;Olson, 1990). In our models, plumes weakened by SW rheology differ most from the classical head-and-tail plume structure (Davaille, 1999;Griffiths & Campbell, 1990;Richards et al., 1989;Sleep, 1990). Weakening of the narrow conduit provides a pathway (lubrication channel) through which hot material can readily rise. In such weakened plume conduits, transport of mass and heat occurs more efficiently. As relatively little thermal buoyancy needs to be built up to drive the plume, no head-and-tail geometry is formed. Indeed, a number of studies have shown that the conduit radius is proportional to η 1/4 , where η is the viscosity of the hot thermal boundary layer (Griffiths & Campbell, 1990;Olson et al., 1993). Such narrow weakened plume conduits have a shorter lifetime than non-weakened plumes ( Figure 6) and can be more easily diverted by large-scale mantle flow and rheological contrasts in the mantle, as is visualized in the Videos S1-S5. An alternative mechanism for the generation of non-classical (e.g., head-less) plumes is the entrainment of intrinsically dense materials. Thermochemical plumes have been found to display highly variable shapes and ascent styles, and to be generally wider than thermal plumes (Davaille & Vatteville, 2005;Farnetani & Samuel, 2005;Lin & Van Keken, 2006).
Most, if not all plumes in our models are also of thermochemical origin (i.e., due to the entrainment of basalt), with negligible compositional differences between SW and non-SW models (see Table S2 in Supporting Information S1). However, the entrained basalt fractions are rather small (f bs ∼ 0.25), and hence the effects of strain weakening on plume ascent geometry are dominant. At a reference depth of 1,800 km, the basalt fractions in our SW and non-SW models correspond to average plume buoyancy numbers B pl (as formulated in e.g., Davaille, 1999;Kumagai et al., 2008) of 0.1-0.21, see Table S2 in Supporting Information S1. Such low B values indicate that thermochemical effects alone cannot account for variable plume shapes, particularly as basalt fraction for SW and non-SW models are similar.

Mantle Plumes on Earth
On Earth, the mismatch between lower mantle and core adiabats implies a super-adiabatic temperature jump across the CMB of about 1,000-1500 K (Boehler, 1996;Jeanloz & Morris, 1986;Lay et al., 2008). Most mantle plumes on Earth, however, are inferred to have excess temperatures of only 100-300 K (e.g., Albers & Christensen, 1996). When extrapolating such excess temperatures to the lower mantle, temperature anomalies of about 500 K have been inferred at the CMB (Albers & Christensen, 1996), still lower than the expected CMB temperature difference (Boehler, 1996). It has been argued that this mismatch is an indication for the super-adiabatic ascent of plumes (Bunge, 2005), or that plumes rise from the top of a compositionally distinct layer at the base of the mantle (Farnetani, 1997). SW rheology of lower mantle materials could additionally help to explain the discrepancy between expected thermal anomalies and observed thermal anomalies of deep-seated mantle plumes, via the shorter onset times of thermal instabilities (see Section 3.2.4). On Earth, many deep-rooted plumes are thought to ascend within a few tens of million years to the base of the lithosphere (e.g., Torsvik et al., 2021). From our modeled mantle plume upwelling velocities, predicted average ascent times are reduced from ∼200 Myr (for non-SW models) to an ∼110 Myr (with SW). The fastest plumes in our SW models have ascent times of only 30 Myr (for 8 cm/yr rising speed), in contrast to 70 Myr for the fastest non-weakened plumes. Therefore, SW rheology could help to explain the plume rise speeds in the Earth's mantle. An alternative explanation for fast plume ascent involves stress-dependent non-Newtonian rheology in the lower mantle (van Keken, 1997), which is argued to produce significantly reduced plume ascent times -although in combination with larger temperature excesses (van Keken, 1997), which is not the case for weakened plumes in this study. Moreover, recently, several surprisingly "cool" hotspots have been identified on Earth (Bao et al., 2022), attributed to either a passive upwelling, shallow source, or entrainment and small-scale convection. We suggest that SW rheology may be a key ingredient for unconventional deep-sourced "cool" mantle plumes as identified by Bao et al. (2022).
Even though LLSVPs have been commonly linked to plume generation (e.g., Torsvik et al., 2010), this link remains controversial. French and Romanowicz (2015) used whole-mantle seismic imaging techniques to argue for the presence of broad, quasi-vertical plumes (i.e., ∼1,000 km in width) beneath many prominent hotspots. These broad plumes were inferred to be thermochemical in origin and rooted at the base of the mantle in patches of greatly reduced shear velocity (e.g., LLSVPs). Other studies propose that the same seismic structures are actually a collection of poorly-resolved narrow mantle plumes (Davaille & Romanowicz, 2020;Davaille & Vatteville, 2005;Schubert et al., 2004). In this scenario, LLSVPs are in fact a cluster of plumes rather than being made up of stable, wide thermochemical piles with broad plume structures atop (as argued by French & Romanowicz, 2015). In our models, SW rheology promotes the existence of multiple weak and thin plume channels, that could possibly be imaged as plume clusters. Yet, in contrast to the plume bundle hypothesis (Davaille & Romanowicz, 2020;Schubert et al., 2004), such narrow, weakened plumes occur in combination with stable, thermochemical piles in the lowermost mantle, and they rise mostly vertically throughout the mantle. It must be noted that even with state-of-the-art seismic methods, narrow mantle plumes (especially those weakened by SW) are difficult-even impossible-to be uniquely distinguished (e.g., Hwang et al., 2011;Kumagai et al., 2008); hence further methodological advances are needed to convincingly discriminate the effects of strain weakening. Moreover, time-dependency is a key factor when interpreting present-day tomographic images of mantle upwellings. The image of a conduit rising from the CMB all the way to the lithosphere is only valid during part of the plume's lifetime (Davaille & Vatteville, 2005). Therefore, plumes might not be easy to detect in tomographic images, particularly if they are weakened by SW rheology and/or deflected by mantle flow. Due to the SW effect on plume width and lifetime, SW rheology could help to explain the faint seismically slow anomalies -or even the absence of detectable anomalies -beneath several hotspots, such as Louisville, Galapagos, and Easter (Davaille & Romanowicz, 2020).

Future Studies
Several future scientific avenues may be carried out to advance this study on the effect of SW rheology on mantle dynamics. First of all, future studies could advance our implementation by making the SW rheology composition-dependent, causing strain-dependent weakening to mainly occur in ferropericlase-enhanced regions. Subducted oceanic crust at lower mantle conditions does not contain any ferropericlase, but instead, contains much cubic CaSiO 3 perovskite Tschauner et al., 2021;Wicks & Duffy, 2016), which may be intrinsically weak (Immoor et al., 2022). It is further interesting to test this composition-dependent weakening in combination with the existence of ancient bridgmanite-enhanced regions in the mid-mantle (as in Gülcher et al., 2020;Gülcher et al., 2021), which should not exhibit SW due to the absence of ferropericlase, hence promoting their preservation. Moreover, it remains to be explored how strain-dependent and grain size-dependent rheology interact with each other in the lower mantle (see Section 4.1). Finally, with the tracking of the deformation matrix in full tensor form, additional work can focus on the direction-dependency of the strain ellipse and weakening behavior and their effects on whole-mantle dynamics.

Conclusions
• We implemented a new strain-dependent rheology for lower mantle materials, combining rheological weakening and healing, in numerical models of global-scale mantle convection. • Strain particularly focuses in anomalously hot regions, such as piles atop the CMB and hot mantle plumes, also when rheological healing is applied. • SW rheology is the key ingredient for the weakening of plume channels as well as forming large thermochemical piles • Second-order effects of SW rheology, caused by the changing mantle dynamics due to a reduction of viscosity in the lower mantle, are higher equilibrium mantle temperatures and the significantly higher average mantle flow velocities. • Weakened mantle plumes form narrow lubrication channels in the mantle through which hot material readily rises, and they have shorter lifetimes. • This new kind of plume dynamics may explain moderate plume excess temperatures beneath hotspots (only up to 200-300 K), given the much larger temperature difference across the core-mantle boundary (∼1,000 K).

Acronyms
Myr million year Gyr billion year LBF load-bearing framework IWL interconnected weak layers CMB core-mantle boundary LLSVP large low shear-wave velocity province SW strain-weakening

Data Availability Statement
The open-source StagLab toolbox (Crameri, 2018) was used for detecting different mantle domains in the numerical models, creation of histogram data (Figures 6 and 7;Figures S6 and S7 in Supporting Information S1), and creating the Videos S1-S5. The new mantle domain detection scheme (as discussed in Supporting Information Text S4 in Supporting Information S1) is implemented in a STAGLAB 6.0 version, and it is available on https://doi.org/10.5281/zenodo.6801104 (STAGLAB-OS release for SW paper (G-cubed, 2022), v1.0.0 https://doi.org/10.5281/zenodo.6801104). Moreover, the open-source Python module StagPy (https://stagpy.readthedocs.io/en/stable/, last access: 17 July 2021) was also used for post-processing of the numerical data and production of radial and temporal profiles (Figures 4 and 5; Figures S2, S4, and S5 in Supporting Information S1). The numerical code is available by reasonable request to Paul J. Tackley. All the data corresponding to the Figures and Tables shown in this paper, as well as the input files for all model runs, are given on https://doi.org/10.5281/zenodo.6805415. As the raw data of the numerical experiments of this paper are too large to be placed online (several terabytes), they can be requested from the corresponding author (Anna J. P. Gülcher). Visitor's Program. We thank Lukas Fuchs for helpful discussions related to this work, Fabio Crameri for discussions related to visualizations and colormaps, and Albert de Montserrat for sharing data from de Montserrat et al. (2021) for Figure S4b in Supporting Information S1. We thank Juliane Dannberg for constructive comments that helped to improve the manuscript considerably. All numerical simulations were performed on ETH Zürich's Euler cluster. For the figures showing 2D snapshots of the models, we used the open-source software ParaView (http://paraview.org, last access: 2 September 2021), and the Videos S1-S5 were made using the geodynamics diagnostics toolbox StagLab (Crameri, 2018). Several perceptually uniform scientific color maps (Crameri, 2018, https://doi. org/10.5281/zenodo.1243862) were used to prevent visual distortion of the figures.