Sluggish thermochemical basal mantle structures support their long-lived stability

Domain: earth_science
Published: 2024-11-19
Editorial outcome: minor_revision

Abstract
Large low shear-wave velocity provinces (LLSVPs) in the lowermost mantle are the largest geological structures on Earth, but their origin and age remain highly enigmatic. Geological constraints suggest the stability of the LLSVPs since at least 200 million years ago. Here, we conduct numerical modeling of mantle convection with plate-like behavior that yields a Pacific-like girdle of mantle downwelling which successfully forms two antipodal basal mantle structures similar to the LLSVPs. Our parameterized results optimized to reflect LLSVP features exhibit velocities for the basal mantle structures that are ~4 times slower than the ambient mantle if they are thermochemical, while the velocity is similar to the ambient mantle if purely thermal. The sluggish motion of the thermochemical basal mantle structures in our models permits the notion that geological data from hundreds of millions of years ago are related to modern LLSVPs as they are essentially stationary over such time scales.

Introduction
The makeup and origins of Earth’s deep and inaccessible interior have long been mysterious. Ever since the detection by seismic tomography of the two antipodal large low shear-wave velocity provinces (LLSVPs) — regions in the lowermost mantle in which the shear-wave velocity is lower than that of the ambient mantle1 — their origin has been extensively debated. Although a purely thermal origin of LLSVPs cannot be excluded2, multiple sources of geophysical data indicate that the LLSVPs are thermally hot and should also be chemically distinct3–6. Whether the LLSVPs derive most or part of their density from a primordial layer of dense materials (such as from the crystallization of magma ocean7 or moon-forming giant impact8) and/or a layer of subducted oceanic crust growing over time, the antiquity of their current morphology into two piles most likely developed at some later point in Earth’s history1. Yet despite great advances in recent decades, many features of the LLSVPs remain enigmatic. There is a great debate about the stability and antiquity of the LLSVPs. In general, there are two end-member interpretations where the LLSVPs are fixed, ancient structures9,10 or dynamic, young structures11–13.

Various types of geological data appear to suggest the African LLSVP has been essentially fixed since 200–300 million years ago (Ma) (i.e., supercontinent Pangaea age), including large igneous provinces (LIPs)10, palaeogeography14, and net characteristics of plate motions15. Some studies have even argued for LLSVP stability back to >500 Ma or even since early in Earth’s history9,16. Alternatively, recent work has challenged LLSVP stability as far back as ca. 200–300 Ma Pangaea age. Numerical modeling of basal mantle structures (BMSs; referring to LLSVPs created in models) parameterized by plate kinematics and thus the mobility of subduction zones at the surface and slabs sinking into the mantle suggest that the lower mantle is more dynamic than previously considered12. Although BMSs are suggested to be very actively dynamic in recent studies, this conclusion is based on imposed plate reconstructions whose accuracy decreases back in time12. Therefore, a fully dynamic approach provides a viable alternative approach.

In this study, we conduct numerical modeling of BMSs to assess their mobility under self-consistent mantle convection with plate-like lithospheric behavior as the top boundary condition. We explore the parameter space that does or does not permit the creation of two antipodal BMSs akin to the areal extent of the LLSVPs. Those results that successfully satisfy such criteria are able to constrain the relative velocity of the BMSs and the ambient mantle, thereby testing the validity of LLSVP stability on 100 million years (Myr) and 1 billion years (Gyr) timescales as both have been proposed by competing models.

Results
Role of plate tectonics in the topology of basal mantle structures
Three-dimensional thermochemical mantle convection models in spherical geometry on the Yin-Yang grid17,18 were conducted with plastic yielding to generate plate boundaries at the top19,20 and an initially even basal layer of dense materials in an attempt to generate BMSs over an evolutionary period of 4.56 Gyrs (comparable to the history of the Earth; “Methods”). As our goal is to investigate the stability time scale of LLSVPs, we first investigate the reaction of the topology of BMSs to plate tectonics by varying the strength of the oceanic lithosphere. We then analyze the mobility of BMSs in an Earth-like model where a Pacific-like subduction girdle and antipodal degree-2 BMSs form.

In our models, slabs penetrate into the lower mantle21 and push aside the BMSs. When plastic yielding is taken into account, two antipodal BMSs are split by a great subduction girdle (Fig. 1d), which is consistent with the present-day Earth where the two LLSVPs are separated by the Pacific “ring of fire” of subduction at the surface as well as with previous 3D models12,22. In 2D convection models with plastic yielding, two antipodal slabs could cut the layered BMSs and then push them into two piles22. In our work, we achieve a 3-D Earth-like model with quasi-stationary BMSs that maintain degree-2 structure for about 630 Myrs when yield stress is set to 100 MPa (Fig. 1d).

To evaluate the influence of plate tectonics on BMSs, we vary the yield stress, hence the strength of the oceanic lithosphere from 10–200 MPa, which controls the plate tectonics configuration20. Figure 1 shows snapshots of the surface viscosity, BMSs, plumes, and subduction. In our series of models, the observed planforms can be divided into four groups. (i) Ultra-low yield stress (30 MPa) yields four subduction systems that stir BMSs, and no piles are observed on the CMB. (ii) Low yield stress (50 MPa) yields a triple-junction-like subduction zone that is concentrated in a semi-sphere and cannot migrate because of the hindering of intense upwelling from the margin or interior of the one-pile BMS with high topography. (iii) Intermediate yield stress (70–100 MPa) yields more than one subduction girdle that cuts the BMSs and pushes them into two piles, where the distance between them increases with increasing yield stress. (iv) High yield stress (200 MPa) yields one circular pile that is pushed against by the downwelling. We further find that successful Earth-like models with two antipodal piles covering the CMB with areas of ~25% similar to LLSVPs occur in our study within a narrow range of moderate yield stress (90–100 MPa) (Fig. 2).

More sluggish thermochemical LLSVPs than thermal LLSVPs
To quantify the mobility of LLSVPs, we calculate the horizontal velocities of BMSs and the ambient mantle in our simulations (Fig. 3). The horizontal velocities of the BMSs are ~4 times slower than that of the ambient mantle at the same depth in the lowermost mantle. These sluggish BMSs could be maintained for several hundreds of million of years (Supplementary Fig. 1). For comparison, we formulate an isochemical convection case to investigate the mobility of purely thermal LLSVPs2. A great subduction girdle akin to the Pacific ring of fire subduction girdle and plume bundles can be obtained with the same yield stress as the thermochemical convection models (Fig. 4a). We define plumes (purely thermal LLSVPs) as the areas of temperature higher than Tp(z) (see “methods”), and find the horizontal velocity of plumes is similar to that of the bulk mantle in the lower mantle (Fig. 4b). Thus, comparing Figs. 4, 3, the mobility of purely thermal LLSVPs is higher than thermochemical LLSVPs.

To further clarify such mobility contrast, we calculate the vote maps (Methods) of BMSs and plumes for the preferred thermal and thermochemical models over the last 630 Myrs (Fig. 5). The main parts of the BMSs at the bottom of the mantle remain stable during this period of time, while the margins of BMSs are deformed and migrate depending on the push from the sinking slabs (column 1 of Fig. 5a, b), which is consistent with the temporal evolution of the BMSs (Supplementary Fig. 2). In contrast, plumes in the deep mantle of the thermal model are so mobile that almost none of them can be maintained for such a period of time compared to the BMSs in the thermochemical model (Column 3 of Fig. 5a, b). Our vote maps, therefore, indicate that the thermochemical LLSVPs in the lower mantle are more sluggish in their lateral motion than thermal LLSVPs (Fig. 5a, b). These differences between the thermal and thermochemical models in the deep mantle may also affect the temporal and spatial distributions of plumes in the shallow mantle. Since the roots of these plumes are generally associated with the hot basal mantle structures, the plumes in the shallow mantle are temporally less intermittent and spatially more stationary in the thermochemical model than those in the thermal model (columns 2 and 4 of Fig. 5a, b).

Influence of physical properties of basal mantle structures on model space
The topology of BMSs in our models (i.e., 1, 2, or no piles) is affected by the distribution of the subduction girdle (Fig. 1). However, in turn, the pattern of BMSs may also affect the plate tectonic configuration. We thus adjust the buoyancy ratio, which previous work has shown to be the first-order factor dominating the configuration of BMSs12,18,23 regardless of the vigor of mantle convection to change the pattern of BMSs. We vary the buoyancy ratio (B) from 0.18–0.30, corresponding to the chemical density contrast from 74–124 kg/m3 (Supplementary Fig. 3), which is consistent with the estimate of density anomalies in the lower mantle4, and keep yield stress and chemical viscosity contrast constant (100 MPa and 10, respectively). We find the BMSs are fully entrained by ambient mantle materials, and subduction forms a great-circle girdle in small buoyancy ratio cases (B = 0.18 and 0.20). In large buoyancy ratio cases (B = 0.22, 0.26, and 0.30), subduction also forms a great-circle girdle, but BMSs form in 2 piles with different topographies and shapes. Thus, the adjustment of the buoyancy ratio in large buoyancy ratio cases has effects on the topography of BMSs but limited effects on the position of subduction zones, hence the plate tectonic configuration, and the stability of BMSs (Supplementary Figs. 3, 4).

Tomography study indicates that LLSVPs contain more bridgmanite than the ambient mantle4. MgSiO3 bridgmanite is more viscous than MgO periclase by up to 1000 times24,25. A chemical viscosity contrast of 10 for BMSs in our preferred model can decrease the viscosity difference between BMSs and the ambient mantle because the BMSs are hot and thermally less viscous26. We vary the chemical viscosity contrasts between the BMSs and the ambient mantle from 1–100 to test the influence of chemical viscosity on the pattern of BMSs (Supplementary Fig. 5). With increasing chemical viscosity contrast, the topography of BMSs rises, and the area of BMSs decreases, due to the decreasing efficiency of mantle convection and hence decreasing heat loss efficiency26. When the chemical viscosity contrast is >30, BMSs form three piles, and the position of subduction changes considerably. BMSs form two piles when the chemical viscosity contrast is <30. Note that the horizontal velocity of BMSs in these cases remains notably slower than that of the bulk mantle, indicating these BMSs are relatively sluggish (Supplementary Fig. 6).

Methods
Thermochemical mantle convection models
We assume an anelastic and compressible fluid with an infinite Prandtl number, and we solve the non-dimensional equations of mass, momentum, heat, and chemical conservation in a 3D spherical shell on the Yin-Yang grid40 using the finite difference/volume code StagYY17. The main model setup used here is similar to that published previously18, and we further take plastic yielding into account to generate plate boundaries19,20. The reference Rayleigh number Ras is expressed as:

Ras = \rho s g\alpha \eta s \Delta T s D3 / \kappa 0 s, (1)

where \rho s is the surface density, g the gravity, \alpha the thermal expansion, \eta s the surface viscosity and \kappa 0 s the surface thermal diffusivity. In this study, Ras is fixed to 1 \times 108.

Temperature-, depth-, and composition-dependent viscosity coupled with plastic yielding is used. Meanwhile, we imposed a viscosity jump of 30 at 660 km depth. The viscosity is described fully by

\eta eff = 1 / (1/\eta(T,d,C) + 1/\eta Y), (2)

with

\eta(T,d,C) = \eta 0 (1+29H(d-660)) exp(VaD) exp(EaT+T off / (T+T s)) + KaC, (3)

and

\eta Y = (\sigma 0 + \sigma_ z) / 2\epsilon_ \Pi, (4)

where \eta Y is the yield viscosity, \sigma 0 the surface yield stress, \sigma_ z the yield stress gradient, which is set to 0 to introduce transform faults as suggested by previous study41, and \epsilon_ \Pi the second invariant of the strain rate tensor. H is the Heaviside step function, Va the non-dimensional activation volume, Ea the non-dimensional activation energy, T off the offset temperature, which is fixed to 0.88\Delta T s to reduce the temperature-induced viscosity variation at the surface, C the composition field, and Ka the parameter controlling the chemical viscosity contrast \Delta \eta C = exp(Ka). No prescribed plate boundaries are imposed and the yield stress help to generate plate-like behavior at the top of the domain. To avoid numerical problems, the viscosity is truncated between 10-3 and 105 of the reference viscosity.

Initial condition and model setup
To introduce basal mantle structures (BMSs), we initially setup an even dense layer above CMB whose thickness is 0.06 (173 km in dimensional unit)42,43. The grid resolution is 64 \times 192 \times 64 \times 2 with grid refinement at both the surface, the bottom, and at 660 km depth (the ringwoodite phase change to bridgmanite). The boundary conditions of the bottom and the top are free-slip, and the temperature of the top is set to 0.12 (300 K), and the bottom is 1.5 (3750 K) (Supplementary Fig. 8a). The model includes both internal and basal heating. The internal heating rate is 7.08 \times 10-12 W/kg, as suggested in ref. 44. Given the LLSVPs possibly enriched in heat-producing elements, we assume the heating rate of the BMS is ~7 times that of the ambient mantle materials26,45. The superadiabatic thermal boundary layers are located at the surface and the bottom, the thickness of which is ~87 km. Random temperature perturbation is implemented at both boundaries. The curvature of the shell f defined as the ratio between the radius of the core and the total radius is set to 0.55 corresponding to the Earth’s value. We use a tracer ratio method to track the compositional field46. Two types of tracers are used to represent the compositional field, where 0 corresponds to the ambient mantle, and 1 corresponds to BMS material. The number of total tracers of all models is 24 million, corresponding to about 15 tracers per cell, which is sufficient for simulating material entrainment46. The composition field is calculated from the ratio between BMS tracers and total tracers. The initial thickness of the BMS material is 0.06 (corresponding to 173 km), and the initial volume of the BMS material is ~3% of the mantle42,43. The density contrast between BMSs and the ambient mantle is controlled by the buoyancy ratio

B = \alpha s \rho s \Delta T s / \Delta \rho c, (5)

where \Delta \rho c is the density contrast between BMSs and the ambient mantle.

Detection of plumes and slabs
To determine plumes and subduction zones in our models, we calculate threshold temperatures of plumes and slabs following those of ref. 47, which assumes that a point is located within a plume (subduction zone) if the temperature at this point is higher (lower) than:

T p(d) = T avg(d) + fp(T max(d) - T avg(d)), (6)

T s(d) = T avg(d) + fs(T min(d) - T avg(d)), (7)

where T avg(d), T max(d) and T min(d) are the horizontally averaged, maximum, and minimum temperature at depth d, and fp and fs are fixed to 0.5 and 0.1, respectively, in this study.

Horizontal velocities of BMSs and plumes
The horizontal velocities of BMSs and plumes are defined as the root-mean-square averaged velocity within the BMSs and plumes at each depth by

V HBMSs(d) = \sqrt{ (\int C>0.5 V H2 dA) / (\int C>0.5 dA) }, (8)

V HPlumes(d) = \sqrt{ (\int T>Tp(d) V H2 dA) / (\int T>Tp(d) dA) }, (9)

where V H = (V\theta, V\phi) is the horizontal velocity at the point of (d, \theta, \phi), \theta colatitude, \phi longitude, dA infinitesimal area. The temporal variation of the horizontal velocity of BMS is calculated at each timestep by

V HBMSs = \sqrt{ (\int C>0.5 V H2 dV) / (\int C>0.5 dV) }, (10)

where dV is infinitesimal volume.

Vote maps of BMSs and plumes
To quantify the stability of BMSs, we calculate vote maps of BMSs following a similar strategy of ref. 48. For a given depth of the composition field, a binary threshold of 0 (no BMSs) or 1 (BMSs) is applied to every latitude-longitude pixel in each time slice within the period of time we selected, followed by a pixel-wise summation of all N time slices within the period selected. The pixels of the vote map hence take values between 0 and N. We then divided the value of each pixel on the vote map by N to get a normalized vote map on which all pixels take values between 0 and 1. A value of 0 on a vote map indicates no BMSs have existed in this pixel during the whole selected period of time, while a value of 1 indicates that BMSs have existed in this pixel at some point during the whole selected period of time, and a value between 0 and 1 indicates the normalized time length that BMSs have existed in this pixel compared to the selected period of time. Vote maps of plumes are calculated by using the same method in the temperature field.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of Sluggish basal mantle structures support their long-lived stability, submitted to Nature\nCommunications by Shi et al. 2023\nThe two antipodal tomographically-observed large low-shear velocity provinces (LLSVPs) are the most\ndistinct structures in Earth's interior. Although they could still be purely thermal, most of the community\nexplain them as intrinsically dense materials that has higher density than the background mantle. Given\nthis assumption and also based on paleomagnetic results from large igneous provinces, hotspots and\nkimberlites, it has been hypothesized that the LLSVPs may be temporally stationary over hundreds of\nmillions of years. But a stable LLSVPs is hard to reconcile with many observations and models, which\nreveals the mobile subduction zone and passive and deforming patterns of LLSVPs. Thus whether LLSVPs\nremain stable or dynamic is one unsolved question in Earth's interior community.\nIn this manuscript, the authors performed global mantle convection simulations to assess whether the\ntwo LLSVPs were stationary over relatively short geological time scales. They found that the average\nhorizontal velocity of the thermochemical piles (LLSVPs in models) in their preferred Earth-like models is\n~4 times slower than the background mantle, indicating the LLSVPs can be stationary over 200-300 Myr\nto allow the connection of surface volcanisms to the deep mantle structure as indicated by geological\nevidence. I commend the authors for a clear written and distinct approach to understand the mobility of\nLLSVPs–which is one of the hottest debates on LLSVPs, and for undertaking global 3D models which is\ncomputationally expensive. However, I find that the calculations are not appropriate for supporting their\nconclusion. Some nontrivial issues are ignored, which would greatly change the results.\nFirst, one of the most critical issues is their models did not include plate reconstruction history, which is\nessential in the validation of numerical models and a key and necessary component in related studies for\ndecades (e.g., Hans-Peter Bunge et al., 2002). Including realistic surface plate motions is crucial as\nindicated by two very recent papers, one in nature (Flament et al., 2022) and another one in nature\ngeoscience (Bodur and Flament, 2023). These two studies clearly show that the mobile LLSVPs can\nstatistically explain associated large igneous provinces and kimberlite eruptions at the surface for the\npast 320 Myr, without a need for stationary LLSVPs. Although the authors’ models with reasonable\nlithosphere yield stress can form two antipodal LLSVPs dividing by the Pacific-like girdle, similar results\nhave been already shown more than 15 years ago (McNamara & Zhong, 2005, Nature). The previous\nstudy further made a detailed comparison of their geodynamic calculation and tomography models,\nwhich is important to validate numerical models, and this comparison also was not considered in the\nstudy.\nAnother major comment is a technical problem: whether it is useful to compare the average velocity of\nLLSVPs with that of the ambient mantle. It is well known that the surface and upper mantle have higher\nvelocity compared to the lower mantle, where the two LLSVP reside. Thus, without any experiments, I\nexpect to see a lower velocity of LLSVPs compared to the average velocity of the bulk mantle. This\nactually can be clearly seen in their figure 3. It might be helpful to compare the LLSVP velocity with the\nambient mantle only in the lower mantle, but from their figure the average velocity between them is\nsimilar (0.3 vs 1 cm/year). I appreciate the authors presented another isothermal calculation that shows\na purely thermal LLSVP have a similar average velocity to the bulk mantle, but how the results will\nchange if you use a different threshold of defining the plumes? A short summary for this comment is that\nI think it is much more meaningful to explore the average lateral motion of plumes from LLSVPs, rather\nthan the horizontal velocity of LLSVPs themselves, to help constrain the mobility vs fixity of LLSVPs.\nWhen exploring the velocity of LLSVPs, the viscosity is one of the most critical parameters. This study\nonly considered cases where LLSVPs are of higher chemical viscosity contrast (1–100 times) given one\nstudy suggesting it may include more bridgmanite than the ambient mantle. However, we have to admit\nthat the viscosity of LLSVPs are still highly uncertain. Because it is widely accepted they are hotter than\nthe ambient mantle, you might want to explore cases with a lower chemical viscosity contrast, such as\n0.01. For a lower viscosity of LLSVPs, I expect the average velocity would be similar or even faster than\nthat of the average lower mantle.\nAnother specific comment is I am a bit concern of the time scale in your reference 3D model. It feels like\nthe whole simulation is modelled for only a few hundred million year, rather than 4.5 Gyr. I am not a user\nfor StagYY, but I suggest doubling check the time scale during dimensionalization. The time scale in other\n3D StayYY models looks more realistic and dynamic. If there is not an issue, it may be related to the\nrheology structures or initial conditions. In terms of visualizations, the plume structures in your uploaded\nvideo seems too massive, even bigger than the downwelling structures. This is in contrast with many\nnumerical models (they also used plasticity), and using a different threshold for both slabs and plumes\nmay help.\nGiven these, I cannot recommend publication of the manuscript in its current form. Despite my non-\nsupportive review, I’d like to make the positive note that I believe the motivation of this study is very well\nselected. I encourage them to find a way to bring more appropriate dynamic models, particularly with\nplate reconstruction history (now the past 1Ga history is available at Müller et al., 2022), and improved\nalgorithm to address this problem.\nI have also several detailed comments as below that may be helpful.\n1. Line26: inaccessible，not largely inaccessible\n2. Line35, a debate or debates\n3. In the introduction you should first state that the LLSVPs could be purely thermal\n4. Line42, these work challenge the assumed “mobility” or “stability” of LLSVP, not their longevity\n5. Line46-48: the main question pointed out here is indeed nicely explained by the work by Flament et\nal. last year and also the one in NG 2023\n6. Line77-78: this might be wrong, Don Bower’s phd work show quasi-stationary pile in both 2D and 3D\n7. Line 78, most Earth-like?\n8. Line88-89, it is hard to believe that the upwelling can hinder the subduction zones in 3D, why is that?\nThis was not found in previous models to my knowledge.\n9. Line94, 90-100 MPa or 70-100 MPa\n10. Line104, how does it influence the result if you change the threshold defining plumes, like\nsmaller or larger threshold\n11. Line115, why only 0.18-0.3, what happen for larger buoyancy, such as 0.4 or 0.5?\n12. Line120-122: the limited effect is because of the buoyancy is not high, and it might be very different\nusing 0.4 or 0.5\n13. 124-130, as I mentioned, you need to explore LLSVP of lower viscosity, e,g., 0.01 contrast\n14. Line154-156: I am very surprised to see the subduction did not fully divide the pile until 3.7 billion\nyears. This might suggest the viscosity structure may be problematic, and previous 3D model takes 200\nMyr from a layer to degree two structures.\n15. 166-167, long-term stable pattern of mantle convection does not directly indicate the stability of\nLLSVPs\n16. Line171-172, show clear evidence that it is easier to move slab than BMSs\n17. Line177-178: this might be your slabs are not that strong\n18. Line180-181: too strong, and not consistent with many previous models\n19. Fig2. How is the case with 200 MPa? Why it gets smaller area?\n20. Fig3. As you can see the velocity in lower mantle is much lower compared to the upper, and note the\nLLSVP velocity is similar to the average velocity of lower mantle.", "reviewer_2": "LLSVPs probably represent a combination of primordial Fe-rich material and dense recycled oceanic\ncrust, and exploring their potential stability over hundreds of million years is of vital importance of\nunderstanding the workings of our planet. Testing their stability by numerical modelling is highly\ndependent on parameterization and I welcome this study that explore their stability under different\nparameter space. I am not an expert on StagYY modelling but I find their parameterization and\nconclusions to be very sound. The paper is well written, the rationale is very clear, it is well referenced\nand I truly enjoyed reading this paper.", "reviewer_3": "This article presents numerical models of mantle convection to examine the stability and longevity of\nbasal mantle structures (BMSs) which serve as a proxy for LLSVPs in the real Earth. This could be an\nexciting result if additional evidence can be provided to support claims made in the paper.\nThe authors run a series of models with various parameters to determine the longevity of the BMSs.\nHowever, more evidence is needed to justify their claim of BMS stability over hundreds of millions of\nyears. As presented, the work is not so convincing. There are some modeling details which should be\nspecified to further understand the model setup. There are also some grammatical and style issues\nwhich could be fixed to enhance readability.\nThe authors show that 2 piles can be formed at the CMB with various parameters. However, there should\nbe more convincing evidence that these results show longevity or stability of these BMSs.\nThe authors claim that the velocity in the BMSs is lower than the ambient mantle and surface plate\nvelocities. However, the authors did not show that this lower velocity is consistent over any length of\ntime. The velocity measured in the BMSs could be a measure the internal convection rate (see >>>>) and\nnot a measure of how the BMS are moving across the CMB. The paper does not provide convincing\nenough evidence that the velocities measured reflect the stability of the BMSs. The authors should\nexplain in more detail how the velocity is being calculated. It would be helpful to include a plot of how\nthe velocities are changing over time since that would be a better indication of pile stability.\nThe authors should show snapshots of a model with degree-2 structure at various times to show the\nstability of the BMSs. It is difficult to tell exactly from video how long the 2 piles are lasting in the same\nplace at CMB.\nThere are many model details that are missing or poorly described to be able to be confident in the\nresults as presented. The models are described as similar to previous research but what exactly does this\nmean? Are the model parameters the same? the resolution? Are the plate boundaries prescribed in the\nmodel and then plastic yielding allows for subduction ? How are the boundaries determined? What is\nthe viscosity profile in the model that would allow for slabs to easily penetrate into the lower mantle and\npush aside the BMSs. Are all models run for the same dimensional time?\nThe authors quantify the mobility of the BMSs by calculating the horizontal velocities of the BMSs (is that\nvelocity inside the BMS or their edge motion as whole pile?). They compare these with a purely thermal\nmodel. They compare the motion of thermal plumes to the ambient mantle but previous work has\nshown plumes to be influenced by mantle motion. A better comparison would be thermo plume to\nthermochemical plumes rather than the plumes to piles. Or more description is needed to justify how\nthermo plume to thermochemical pile comparisons support stability.\nFigure images are clear but captions and divisions could be made more clear. In particular there are\nmultiple references to Figure 1d but there is no \"d\" labelled on figure 1. In fact Figures 1, S1, S2 are all in\ncolumns and should referenced as such in the paper (or labeled with a, b, c, etc).\nLine by line comments:\nLine 62: add \" \"even\" basal layer\" as it could imply 2 piles are setup to start\nLine 70. more like this like: \"like a model where a Pacific-like ... \"\nLine 72: explain how the slabs can penetrate and push aside BMS is there no viscosity jump at 660 ?\nLine 80: no Figure 1d - no part d\nLine 100: typo: \"formulate an\" instead of formulate another\nLine 103: change \"compared to\" to \"as.. the thermochemical convection models).\nLine 104; change \"As we\" to just \"We\"\nLine 104: define what is the value of Tp(z) = ??\nLine 106: refer to the figure you want me to look at\nLine 112: Figure 1 refer to the columns 1 and 4 you want me to look at\nLine 113: change \"which from previous work is ..\" to \"which previous work has shown\"\nLine 118 to 120: For the buoyancy ratios, list the values, a range wasn't tested just specific values (B=0.18\nand 0.20 or B=0.22, 0.26, 0.30)\nLine 119 to 120 move \"form in 2 piles\" so it reads \"but BMSs form in 2 piles with different topographies\nand shapes\".\nLine 124-127:: Change order sentences to enhance clarity\n\"We vary the chemical viscosity contrasts of the BMSs from 1-100 to test the influence of chemical\nviscosity on the pattern of BMSs (sup Fig 2). Tomography indicated that LLSVPs contain more\nbridgmanite than the ambient mantle. MgSiO3 bridmangite is more viscous than the MgO periclase by\n10-1000 times. With increasing chemical viscosity ..... \"\nLine 143 \"are in this range\" list the range here\nLine 170 change \"In consideration of \" to \"Considering \"\nLine 175: change \"in compensation of \" to \"after\" Pangaea\nLine 178: change \"Considering also the less\" to \"Considering the lower\""}, "author_response": "We would like to thank the reviewers for their constructive and helpful comments. This is our\nresponse to the reviews of manuscript NCOMMS-23-28774, entitled “Sluggish basal mantle\nstructures support their long-lived stability”. Comments from the reviewers are in black and our\nresponses are in blue. Note that the line numbers correspond to the clean version without track\nchanges.\nResponse to Reviewer #1:\nReview of Sluggish basal mantle structures support their long-lived stability, submitted to Nature\nCommunications by Shi et al. 2023\nThe two antipodal tomographically-observed large low-shear velocity provinces (LLSVPs) are the\nmost distinct structures in Earth's interior. Although they could still be purely thermal, most of the\ncommunity explain them as intrinsically dense materials that has higher density than the\nbackground mantle. Given this assumption and also based on paleomagnetic results from large\nigneous provinces, hotspots and kimberlites, it has been hypothesized that the LLSVPs may be\ntemporally stationary over hundreds of millions of years. But a stable LLSVPs is hard to reconcile\nwith many observations and models, which reveals the mobile subduction zone and passive and\ndeforming patterns of LLSVPs. Thus whether LLSVPs remain stable or dynamic is one unsolved\nquestion in Earth's interior community.\nIn this manuscript, the authors performed global mantle convection simulations to assess whether\nthe two LLSVPs were stationary over relatively short geological time scales. They found that the\naverage horizontal velocity of the thermochemical piles (LLSVPs in models) in their preferred\nEarth-like models is ~4 times slower than the background mantle, indicating the LLSVPs can be\nstationary over 200-300 Myr to allow the connection of surface volcanisms to the deep mantle\nstructure as indicated by geological evidence. I commend the authors for a clear written and distinct\napproach to understand the mobility of LLSVPs–which is one of the hottest debates on LLSVPs,\nand for undertaking global 3D models which is computationally expensive. However, I find that\nthe calculations are not appropriate for supporting their conclusion. Some nontrivial issues are\nignored, which would greatly change the results.\nResponse: Thank you for acknowledging the potential significance of our study and for the\nconstructive suggestions and comments that helped us to improve the manuscript. Below are our\npoint-by-point responses to your suggestions and comments.\n1) First, one of the most critical issues is their models did not include plate reconstruction history,\nwhich is essential in the validation of numerical models and a key and necessary component in\nrelated studies for decades (e.g., Hans-Peter Bunge et al., 2002). Including realistic surface plate\nmotions is crucial as indicated by two very recent papers, one in nature (Flament et al., 2022) and\nanother one in nature geoscience (Bodur and Flament, 2023). These two studies clearly show that\nthe mobile LLSVPs can statistically explain associated large igneous provinces and kimberlite\neruptions at the surface for the past 320 Myr, without a need for stationary LLSVPs. Although the\nauthors’ models with reasonable lithosphere yield stress can form two antipodal LLSVPs divided\nby the Pacific-like girdle, similar results have been already shown more than 15 years ago\n(McNamara & Zhong, 2005, Nature). The previous study further made a detailed comparison of\ntheir geodynamic calculation and tomography models, which is important to validate numerical\nmodels, and this comparison also was not considered in the study.\nResponse: We thank the reviewer for expanding our perspective and for reminding us of another\nadvanced simulation strategy for surface conditions. This question deserves careful clarification.\nThe reviewer is correct in that plate reconstructions are important for validating the style of mantle\nconvection, including the work presented in this study. Here, we want to clarify that the conclusion\nof our study—the long-term stability of LLSVPs with a suitable lithospheric yield stress—is not\ninconsistent with our knowledge of ancient plate reconstructions. For example, the 200–300 Myr\nstability of our modeled LLSVPs is similar to results directly based on plate reconstructions\n(McNamara & Zhong, 2005). We hope to emphasize that the LLSVPs in our models are not\nentirely fixed as one may envisage from some end-member models (e.g., Torsvik et al., 2010).\nInstead, we showed that the bulk LLSVPs are largely stationary, but their edges are deforming and\nmigrating depending on the push of the sinking slabs. This is also what Flament et al. (2022)\ndemonstrated in their model. Another independent verification from plate reconstructions comes\nfrom the study of Conrad et al. (2013) where they found that global plate motions obeyed a\nstationary degree-2 pattern during the past ~250 Myr.\n[Figure Redacted]\nOn the other hand, we also want to make it clear that models with imposed plate kinematics at the\nsurface cannot objectively conclude on the stability of deep mantle structures like LLSVPs, due to\nthe forced nature of the surface boundary conditions. In contrast, a freely moving surface like the\none incorporated in our models is more appropriate, as demonstrated in earlier studies (e.g., Zhong\net al., 2007). Because of this philosophy, the cited works by the reviewer (e.g., McNamara &\nZhong, 2005; Flament et al., 2022) cannot be used to claim the true stability of the LLSVPs, as\nthey all depend on the style of the imposed plate motions. The fact that our result resembles theirs\nreveals that the reconstructed plate motion since as early as 250 Ma is likely reasonable, which is,\nagain, consistent with the finding of Conrad et al. (2013).\nThis brings up another question regarding the pre-250 Ma plate reconstruction: should one trust\nthe recent plate reconstructions managed to extend beyond the supercontinent Pangea or fully\ndynamic simulations like the ones we propose here? Many people realize that these plate\nreconstructions were done with limited constraints, especially from the ocean basins which are the\nkey constraints for the reconstructions for the past 250 Myrs reconstructions. Indeed, the\nreconstruction experts themselves clearly warned about this caveat (e.g., Müller et al., 2022) by\nstating that these early reconstructions enforced a constant mean global plate motion over time\nsuch that oceanic subduction never stopped, even after the assembly of a supercontinent;\nconsequently, published plate reconstructions all feature a dominant degree-2 pattern over time,\nwith a circum- supercontinent subduction girdle reminiscent of the degree-2 convection pattern\nfirst proposed in Zhong et al. (2007). No matter if this assumption is valid or not, such a long-term\ndegree-2 pattern is again similar to our finding of LLSVPs stability for up to 600 Myr. We therefore\nconclude that our result is consistent with and supported by approaches using plate reconstructions,\nand that our approach is more appropriate for examining the long-term stability of LLSVPs than\nthose with imposed plate motions—or at least provides a viable alternative approach.\n2) Another major comment is a technical problem: whether it is useful to compare the average\nvelocity of LLSVPs with that of the ambient mantle. It is well known that the surface and upper\nmantle have higher velocity compared to the lower mantle, where the two LLSVP reside. Thus,\nwithout any experiments, I expect to see a lower velocity of LLSVPs compared to the average\nvelocity of the bulk mantle. This actually can be clearly seen in their figure 3. It might be helpful\nto compare the LLSVP velocity with the ambient mantle only in the lower mantle, but from their\nfigure the average velocity between them is similar (0.3 vs 1 cm/year). I appreciate the authors\npresented another isothermal calculation that shows a purely thermal LLSVP have a similar\naverage velocity to the bulk mantle, but how the results will change if you use a different threshold\nof defining the plumes? A short summary for this comment is that I think it is much more\nmeaningful to explore the average lateral motion of plumes from LLSVPs, rather than the\nhorizontal velocity of LLSVPs themselves, to help constrain the mobility vs fixity of LLSVPs.\nResponse: We thank the reviewer for this comment. This comment deserves a detailed\nclarification.\nFirst, we agree with the reviewer that the surface and upper mantle have higher velocities compared\nto the lower mantle. There might be some misunderstanding in comparing the velocities between\nLLSVPs and the ambient mantle. In our study, we compared the averaged horizontal velocity\ndifferences between LLSVPs and the bulk mantle at the same depth. We found that the\nthermochemical LLSVPs have a lower velocity of about 1/4 that of the bulk mantle as shown in\nFigure 3 (0.3 vs ~1.2 cm/y), meanwhile, pure thermal LLSVPs have a similar velocity compared\nto the bulk mantle as shown in Figure 4 (~1.4 vs ~1.5 cm/y). This difference between thermal and\nthermochemical LLSVPs suggests that thermochemical LLSVPs are sluggish structures in the\nlower mantle. We have clarified the comparison at the same depth in the revised manuscript (line\n95) and removed the comparison with surface velocity to avoid misunderstanding (line 95).\nSecond, we thank the reviewer for the suggestion to test “a different threshold of defining the\nplume”. We tested different thresholds (0.4, 0.5, and 0.6) to define plumes (Fig. R2), and we found\nthat although the structure of plumes becomes less broad with increasing threshold, the averaged\nhorizontal velocities of plumes are comparable to each other, and no significant differences were\ncaused by changing the threshold. And this observation also applies to the thermochemical plumes.\nFollowing the reviewer’s suggestion, we use a threshold of 0.5 in the revision (instead of 0.4 in\nthe original manuscript) to improve the visualization of our results.\nFig. R2. Thermal convection model with yield stress of 100 MPa. (a) Viscosity at the surface (top row) and\nplume and slabs (bottom second to fourth row). Right column shows the opposite hemispheres of the left column.\nThe plumes are detected by setting fp=0.4 (second row), 0.5 (third row), 0.6 (fourth row), respectively. (b)\nHorizontal velocities of plumes and the bulk mantle.\nThird, we appreciate the reviewer’s suggestion to explore the average lateral motion of plumes.\nBecause an individual plume has a much shorter life cycle compared to the length of time in our\nstudy and is therefore not easy to track, we have made new calculations to plot the vote maps (see\nMethods) of plumes as well as that of the thermochemical LLSVPs for the reference case and the\nthermal convection case to show how the locations of these features change with time (Fig. 5) in\nthe last 630 Myrs of evolution. We can observe that the main parts of the LLSVPs at the bottom\nof the mantle remain stationary during this period of time, while the margins are deforming and\nmigrating with time depending on the push of the sinking slabs. We also plot the vote maps of\nplumes at 300 km, and consistent with our main conclusion, the plumes in the thermochemical\ncase are more stationary compared to those in the thermal convection model (Fig. 5, columns 2\nand 3). The vote maps also indicate that the thermochemical piles in the lower mantle are sluggish\nin their lateral motion compared to the plumes above them, which again can serve as one line of\nevidence to support our main conclusion that these BMSs are sluggish structures in the lower\nmantle. We have added the description of the features mentioned above in the revised manuscript\n(lines 105-113).\n3) When exploring the velocity of LLSVPs, the viscosity is one of the most critical parameters.\nThis study only considered cases where LLSVPs are of higher chemical viscosity contrast (1–100\ntimes) given one study suggesting it may include more bridgmanite than the ambient mantle.\nHowever, we have to admit that the viscosity of LLSVPs are still highly uncertain. Because it is\nwidely accepted they are hotter than the ambient mantle, you might want to explore cases with a\nlower chemical viscosity contrast, such as 0.01. For a lower viscosity of LLSVPs, I expect the\naverage velocity would be similar or even faster than that of the average lower mantle.\nResponse: The viscosity of LLSVPs is affected by both temperature and composition, as we\ndetailed in the Methods. Since LLSVPs are hotter than the ambient mantle, the thermal viscosity\nis smaller than for the ambient mantle, and for the case where chemical viscosity contrast equals\n1, LLSVPs are already less viscous than the ambient mantle, as shown below (Fig. R3).\nFig. R3. Snapshots for viscosity (a) and composition (b) fields with chemical viscosity contrast equals 1.\nLess viscous LLSVPs lead to more mixing with the ambient mantle. If we further decrease the\nchemical viscosity contrast to 0.1(i.e., by one order of magnitude), the stable part of LLSVP\ndecreases, nonetheless overall, the velocity of the stable part remains smaller than that of the bulk\nmantle (Fig. R4). The primordial piles cannot be preserved during the entire run if the chemical\nviscosity contrast decreases to 0.01(i.e., an additional order of magnitude). To summarize, the\nchemically distinct LLSVPs cannot preserve in the lower mantle during the entire run if their\nchemical viscosity contrast is relatively small (e.g., 0.01); and the velocity of stable LLSVPs can\nbe still smaller than that of the bulk mantle even if the overall viscosity is smaller than the bulk\nmantle.\nFig. R4. Viscosity, basal mantle structures, plumes, slabs, and horizontal velocities with chemical viscosity\ncontrast of 0.1. (a) Viscosity at the surface (columns 1 and 4), basal mantle structures (isosurface of the\ncomposition with contour level 𝐶 =0.5, columns 2 and 5), plume (isosurface of non-dimentional temperature\nwith contour level 𝑓 =0.5, columns 3 and 6, orange), and slabs (isosurface of non-dimensional temperature\np\nwith contour level 𝑓 =0.1, columns 3 and 6, blue). Columns 4–6 show the opposite hemisphere of columns 1-\ns\n3. (b) 1-D averaged horizontal velocities of basal mantle structures (BMSs; solid line), and the ambient mantle\n(dotted line).\n4) Another specific comment is I am a bit concern of the time scale in your reference 3D model.\nIt feels like the whole simulation is modelled for only a few hundred million year, rather than 4.5\nGyr. I am not a user for StagYY, but I suggest doubling check the time scale during\ndimensionalization. The time scale in other 3D StayYY models looks more realistic and dynamic.\nIf there is not an issue, it may be related to the rheology structures or initial conditions. In terms\nof visualizations, the plume structures in your uploaded video seems too massive, even bigger than\nthe downwelling structures. This is in contrast with many numerical models (they also used\nplasticity), and using a different threshold for both slabs and plumes may help.\nResponse: We checked the dimensionalization of time, and the model time is 4.56 Gyr (D2/κ, D\nthe mantle thickness, κ the reference thermal diffusivity). In the revised manuscript, we have added\nthe profiles of the temperature and viscosity both at the initial and final steps (Supplementary\nFig.S8), which are comparable to previous studies (e.g., Langemeyer et al., 2022; Nakagawa et\nal., 2010; Li et al., 2014). For the visualizations, plumes and slabs in our plots are defined by\nexcess temperatures (Labrosse, 2002), and we think these hot and cold anomalies are generally\nconsistent with those in previous studies (e.g., Nakagawa and Tackley, 2008; Nakagawa and\nTackley, 2009). The same as our response to the comment above, by following the reviewer’s\nsuggestion, we have tested different thresholds for plumes, and we think the visualizations could\nbe improved by using a bigger threshold (0.5 in the revised manuscript compared to 0.4 in the\noriginal manuscript) while the horizontal velocities of plumes don’t show significant differences.\nWe show the plots for the thermal convection as an example of major comment 2 (Fig. R2), and\nwe have updated all figures in the revised manuscript.\nGiven these, I cannot recommend publication of the manuscript in its current form. Despite my\nnon-supportive review, I’d like to make the positive note that I believe the motivation of this study\nis very well selected. I encourage them to find a way to bring more appropriate dynamic models,\nparticularly with plate reconstruction history (now the past 1Ga history is available at Müller et al.,\n2022), and improved algorithm to address this problem.\nResponse: We appreciate the reviewer for providing the constructive and helpful comments. We\nhope our responses can satisfy the reviewer.\nI have also several detailed comments as below that may be helpful.\n1. Line26: inaccessible, not largely inaccessible\nWe have modified this in the revised manuscript (line 26).\n2. Line35, a debate or debates\nWe have modified this in the revised manuscript (line 36).\n3. In the introduction you should first state that the LLSVPs could be purely thermal\nThank you for this suggestion. In the revised manuscript, we have modified the Introduction to\nstate that the LLSVPs could be purely thermal (lines 29-30).\n4. Line42, these work challenge the assumed “mobility” or “stability” of LLSVP, not their\nlongevity\nWe have replaced ‘longevity’ with ‘stability’ in the revised manuscript (line 44).\n5. Line46-48: the main question pointed out here is indeed nicely explained by the work by\nFlament et al. last year and also the one in NG 2023\nWe acknowledged that the work by Flament et al.(2022) offers an alternative explanation to the\ngeological observations, but we also think it could not exclude the stable LLSVPs such as those\nproposed by previous works. And we have revised this to state that a fully dynamic approach is\nmore appropriate (lines 48-50).\n6. Line77-78: this might be wrong, Don Bower’s phd work show quasi-stationary pile in both 2D\nand 3D\nBower et al. (2013) uses a plate-reconstruction history back to 250 Ma and an initial condition of\ntwo existing blobs which we think may not be able to represent quasi-stationary piles as long as\nthose that we reveal in this study.\n7. Line 78, most Earth-like?\nWe have revised this wording as ‘an Earth-like model’ to refer to the reference model in our work\nwhich contains a subduction girdle similar to the ‘ring of fire’ of subduction and degree-2 antipodal\nBMSs that take up an area of >25% of CMB akin to LLSVPs (line 75).\n8. Line88-89, it is hard to believe that the upwelling can hinder the subduction zones in 3D, why\nis that? This was not found in previous models to my knowledge.\nIn our study, we find that because the LLSVPs at the CMB are stationary, most of the plumes\ngenerate atop these LLSVPs, and these regions usually relate to warm or the divergent plate\nmargins where large-scale subduction zones may not form.\n9. Line94, 90-100 MPa or 70-100 MPa\nWe have clarified this stress range in the revised manuscript. Although intermediate yield stress\n(70–100 MPa) could yield a subduction girdle to form two piles, our best-fit cases are within the\nrange of 90–100 MPa, which also leads to the coverage of the CMB area similar to the observations\nof the LLSVPs (lines 88-90).\n10. Line104, how does it influence the result if you change the threshold defining plumes, like\nsmaller or larger threshold\nWe have tested the effects of the threshold of plumes by varying its value (0.4, 0.5, and 0.6), and\nwe find the horizontal velocities of plumes defined by different thresholds are similar to each other\nwithout significant differences (Fig. R2, see our response to major comment 2).\n11. Line115, why only 0.18-0.3, what happen for larger buoyancy, such as 0.4 or 0.5?\nA buoyancy ratio between 0.18 and 0.3 corresponds to a density contrast between 74 and 128\nkg/m3, which is consistent with geophysical estimates (e.g., Trampert et al., 2004; Mosca et al.,\n2013). Larger values such as 0.4 or 0.5 lead to larger chemical density contrasts between BMSs\nand the ambient mantle, which may not fit the geophysical constraints well.\nWe have run simulations with B = 0.4 and 0.5 as the reviewer mentioned in the comment, and the\nresults are shown below (Fig. R5). Consistent with previous studies (e.g., Deschamps and Tackley,\n2009; Li et al., 2014), a large buoyancy ratio of primordial material leads to more stable and flat\npiles above the CMB without large topography which is not comparable to the LLSVPs. The\nprimordial material is almost layering flatly above the CMB without large topography which is\nquite different from seismological observations, and the CMB area covered by these piles is much\nlarger (>50%).\nFig. R5. Viscosity, basal mantle structures, plumes, and slabs with different buoyancy ratios. Viscosity at\nthe surface (columns 1 and 4), basal mantle structures (isosurface of the composition with contour level 𝐶 =0.5,\ncolumns 2 and 5), plume (isosurface of non-dimentional temperature with contour level 𝑓 =0.5, columns 3\np\nand 6, orange), and slabs (isosurface of non-dimensional temperature with contour level 𝑓 =0.1, columns 3 and\ns\n6, blue). Columns 4–6 show the opposite hemisphere of columns 1-3. From (a) to (b), the values of buoyancy\nratios correspond to 0.40 and 0.50 (other parameters the same as Supplementary Fig. 3).\n12. Line120-122: the limited effect is because of the buoyancy is not high, and it might be very\ndifferent using 0.4 or 0.5\nThe reviewer is correct: a larger buoyancy ratio would strongly alter the lower mantle conditions,\nand therefore affect the mantle convection pattern. But as mentioned in our response to minor\ncomment 11, a large buoyancy ratio of 0.4 or 0.5 would avoid the formation of BMSs comparable\nto the observed LLSVPs, although it may have a more pronounced effect, it is not suitable for our\nstudy (Fig. R5).\n13. 124-130, as I mentioned, you need to explore LLSVP of lower viscosity, e,g., 0.01 contrast\nPlease see our response to your major comment 3.\nLLSVPs are high-temperature regions that have low viscosities due to thermal viscosity contrast.\nAn additional lower chemical viscosity contrast (e.g., 0.01) made the viscosity even much lower,\nwhich makes the primordial material more easily entrained by plumes and mantle flow and\ntherefore more easily mixed with the ambient mantle. In these cases, large reservoirs formed by\nsuch primordial material cannot survive the entire run in our simulations.\n14. Line154-156: I am very surprised to see the subduction did not fully divide the pile until 3.7\nbillion years. This might suggest the viscosity structure may be problematic, and previous 3D\nmodel takes 200 Myr from a layer to degree two structures.\nWe thank the reviewer for pointing this out. In our preferred model with a yield stress of 100MPa,\nthe most recent subduction girdle occurred about 1Gyrs ago, and it takes another 200~300 Myrs\nbefore the formation of two stationary BMSs, and then these BMSs are stable for about 630 Myrs.\nWe have modified this part in the revised manuscript (lines 163-164). In the revised manuscript,\nwe have also added a new plot (Supplementary Fig. 8) to show the 1-D initial viscosity profile and\nthe viscosity profile at the end of the simulation for our preferred case (yield stress of 100MP).\nThese viscosity profiles are generally comparable to previous studies (Deschamps & Tackley,\n2009; Li et al., 2014; Nakagawa et al., 2010; van Heck & Tackley, 2008).\n15. 166-167, long-term stable pattern of mantle convection does not directly indicate the stability\nof LLSVPs\nWe agree with the reviewer that the long-term stable pattern of mantle convection does not directly\nindicate the stability of LLSVPs. However, if they are long-term stable, it is favored by core\ndynamics studies as we indicated in our manuscript (Gubbins et al., 2011; Labrosse et al., 2001;\nLasbleis et al., 2020; Tkalčić, 2015; Wang & Vidale, 2022; Waszek et al., 2011).\n16. Line171-172, show clear evidence that it is easier to move slab than BMSs\nIn the revised manuscript, we have added vote maps for BMSs to show that their main parts are\nstationary in our Earth-like model during the last 630 Myrs (column 1 of Figs.5a and 5c). We also\nnote that although the main parts of BMSs remain stationary, the sinking slabs significantly deform\nand shift the margins of BMSs (column 1 of Figs. 5a and 5c). We think that when the slabs move\nthe margins of BMSs toward the main parts of BMSs, the resistance of BMSs increases and get\nlarger than the ambient mantle. Therefore, the slabs have to move to the slab graveyard away from\nBMSs. In seismic tomography models, LLSVPs are surrounded by such slab graveyards which are\nhigh-velocity anomalies.\n17. Line177-178: this might be your slabs are not that strong\nThe slabs can easily push small piles merging into larger ones (e.g., in the 50 MPa case); but the\nresistance of piles to the slab push increases with the size of piles, such as those in our reference\ncase, where the resistance can prevent the subduction girdle from pushing the two large piles\ntogether.\n18. Line180-181: too strong, and not consistent with many previous models\nThank you for pointing this out. We have removed this sentence and revised this statement to avoid\noverstatement. A subduction girdle may have formed and simply shifted during the last billion\nyears as observed in our self-consistent model, and consistent with some geological observations\n(e.g., Doucet et al., 2019).\n19. Fig2. How is the case with 200 MPa? Why it gets smaller area?\nFor cases with large yield stress, the area covered by BMSs varies slightly varies with different\ncases and with time (Fig. R6), but mostly between 40% and 50% because of the margins deforming\nand shifting with time. We have updated Fig 2 in the revised manuscript accordingly.\nFig. R6. Spatial and temporal distributions of BMSs in cases with different yield stress. (a) 1-D distribution\nof BMSs materials in the mantle. (b) Temporal variation of the area covered by BMSs at the core-mantle\nboundary.\n20. Fig3. As you can see the velocity in lower mantle is much lower compared to the upper, and\nnote the LLSVP velocity is similar to the average velocity of lower mantle.\nThe velocity in the lower mantle is much lower compared to the upper mainly due to the viscosity\njump at 660km. In the lower mantle, the velocity of thermo-chemical LLSVPs is much smaller\nthan the average velocity of the lower mantle at the same depth (e.g., 0.3 vs 1.2 cm/y). In\ncomparison, the velocity of such anomalies in a pure thermal case is similar to that of the average\nvelocity of the lower mantle (Fig. 4).\nResponse to Reviewer #2:\nLLSVPs probably represent a combination of primordial Fe-rich material and dense recycled\noceanic crust, and exploring their potential stability over hundreds of million years is of vital\nimportance of understanding the workings of our planet. Testing their stability by numerical\nmodelling is highly dependent on parameterization and I welcome this study that explore their\nstability under different parameter space. I am not an expert on StagYY modelling but I find their\nparameterization and conclusions to be very sound. The paper is well written, the rationale is very\nclear, it is well referenced and I truly enjoyed reading this paper.\nWe greatly appreciate the reviewer for making this critical recommendation.\nResponse to Reviewer #3:\nThis article presents numerical models of mantle convection to examine the stability and longevity\nof basal mantle structures (BMSs) which serve as a proxy for LLSVPs in the real Earth. This could\nbe an exciting result if additional evidence can be provided to support claims made in the paper.\nThe authors run a series of models with various parameters to determine the longevity of the BMSs.\nHowever, more evidence is needed to justify their claim of BMS stability over hundreds of millions\nof years. As presented, the work is not so convincing. There are some modeling details which\nshould be specified to further understand the model setup. There are also some grammatical and\nstyle issues which could be fixed to enhance readability.\nResponse: We greatly appreciate the reviewer’s constructive and supportive comments. We have\nrevised the manuscript with additional analysis. We also have added more modeling details to\nshow the model setup and carefully revised the manuscript to solve the grammatical and style\nissues as suggested to enhance readability, as suggested. And we hope our response and the\nrevision can satisfy the reviewer.\n1) The authors show that 2 piles can be formed at the CMB with various parameters. However,\nthere should be more convincing evidence that these results show longevity or stability of these\nBMSs.\nResponse: We have revised our manuscript and added more evidence to hopefully make our study\nmore convincing to show the stability of these BMSs.\nTo show the longevity and stability of the degree-2 BMSs in our preferred case, we further analyze\nthe harmonic spectrum in the lower mantle in Supplementary Figure 7. We find the dominant\nharmonic structure of degree-2 can be maintained even as long as almost 1Gyr after its formation\n(associated with the degree-2 subduction girdle at the surface). We also plot the vote maps for both\nthermal and thermochemical cases for plumes and piles (Fig. 5), which show that the main part of\nthe chemical piles remains stable during the last 630 Myr as in the reference case in which degree-\nwe can observe that the plumes in the thermal model are more mobile than those in the\nthermochemical model. These observations together indicate that the 2 piles can be maintained for\na long period of time.\n2) The authors claim that the velocity in the BMSs is lower than the ambient mantle and surface\nplate velocities. However, the authors did not show that this lower velocity is consistent over any\nlength of time. The velocity measured in the BMSs could be a measure the internal convection rate\n(see >>>>) and not a measure of how the BMS are moving across the CMB. The paper does not\nprovide convincing enough evidence that the velocities measured reflect the stability of the BMSs.\nThe authors should explain in more detail how the velocity is being calculated. It would be helpful\nto include a plot of how the velocities are changing over time since that would be a better indication\nof pile stability.\nResponse: We thank the reviewer for this suggestion. In the revision, we try to clarify the sluggish\nmotion of thermochemical BMSs and the stability of the BMSs over time.\nFirst, following the reviewer’s suggestion, we have added a plot of the temporal variations of the\nhorizontal velocity of BMSs and the bulk mantle which shows how these velocities are changing\nover time (Supplementary Fig. 1). The same chemical viscosity contrast is used in all the cases\nshown in this figure, and the horizontal velocities of BMSs are different. We therefore think the\nhorizontal velocity we calculated can reflect the mobility of BMSs; and although the internal\nconvection rate can not be excluded, it may not play a major role here. This plot also shows that\nthe thermochemical BMSs are more sluggish than the bulk mantle most of the time after their\nformation in the lower mantle. (lines 95-96)\nSecond, following the reviewer’s suggestion, we have added the description of the details that the\nvelocity is being calculated in the Methods section (lines 255-264).\nThird, to further reflect the stability of the BMSs, we have added a plot (Fig. 5) to show the vote\nmap of plumes and BMSs during the last 630 Myrs in our Earth-like model and the thermal\nconvection model as in our response to the above major comment 1. We think that all these can\nserve as pieces of evidence together to support the stability of these piles and that they are sluggish\nstructures in the lower mantle.\n3) The authors should show snapshots of a model with degree-2 structure at various times to show\nthe stability of the BMSs. It is difficult to tell exactly from video how long the 2 piles are lasting\nin the same place at CMB.\nResponse: Thank you for this suggestion. In the revised manuscript, we have added snapshots of\nthe temporal evolution of BMSs in the final 1.5 Gyrs (Supplementary Fig. 2) and the temporal\nevolutionary histories of the dominant degree of the harmonic spectrum of the lower mantle\n(Supplementary Fig. 7). These plots, together with the vote maps for the thermochemical piles at\nthe CMB (Fig. 5), can serve to support the claim of our study that the two piles can last for a long\nperiod of time, and further indicates that the main bodies of the BMSs remain stable while the\nmargins vary with time (Fig. 5). We think all these figures can together give a clearer view than\nthe video of the original manuscript.\n4) There are many model details that are missing or poorly described to be able to be confident in\nthe results as presented. The models are described as similar to previous research but what exactly\ndoes this mean? Are the model parameters the same? the resolution? Are the plate boundaries\nprescribed in the model and then plastic yielding allows for subduction ? How are the boundaries\ndetermined? What is the viscosity profile in the model that would allow for slabs to easily penetrate\ninto the lower mantle and push aside the BMSs. Are all models run for the same dimensional time?\nResponse: In the revision, we have added more detailed descriptions of the model in the Methods,\nincluding physical parameters, initial condition, model setup, and details in calculating horizontal\nvelocities and vote maps of BMSs and plumes.\nThe main parameters and 3-D model setup used in this study are similar to a previous study (Li et\nal., 2014) and we further introduce yield stress in this study to generate plate boundaries (Tackley,\n2000; van Heck & Tackley, 2008, 2011). The 3D spherical shell is discretized by Yin-Yang\ngrid(Kageyama & Sato, 2004) and the resolution is 64 × 192 × 64 × 2. No plate boundaries are\nprescribed and the plate boundaries are only generated by rock yielding (Tackley, 2000; van Heck\n& Tackley, 2008). The viscosity profile of 100 MPa case is shown in Supplementary Figure 9b.\nAll the models are run for the same dimensional time of 4.56 Gyrs.\nWe have modified the revised manuscript to clarify the details mentioned above.\n5) The authors quantify the mobility of the BMSs by calculating the horizontal velocities of the\nBMSs (is that velocity inside the BMS or their edge motion as whole pile?). They compare these\nwith a purely thermal model. They compare the motion of thermal plumes to the ambient mantle\nbut previous work has shown plumes to be influenced by mantle motion. A better comparison\nwould be thermo plume to thermochemical plumes rather than the plumes to piles. Or more\ndescription is needed to justify how thermo plume to thermochemical pile comparisons support\nstability.\nResponse: We calculated the horizontal velocity inside the BMSs in the lowermost mantle, and\nwe have added the details of how we calculate the horizontal velocity of BMSs and plumes in the\nMethods of the revised manuscript.\nWe thank the reviewer for the suggestion to compare the motions between thermal and\nthermochemical plumes. One of our main purposes in the original manuscript was to show that\nthermochemical LLSVPs are sluggish basal mantle structures that are more stable than the ambient\nmantle, and even more stable than the plumes, so that they can serve as a geological framework as\nsuggested by previous studies (e.g., Torsvik et al., 2014). In the revised manuscript, instead of\ntracking individual plumes, we plot the vote maps for both plumes and the BMSs (Fig. 5) to make\ncomparisons for how likely the BMSs and plumes can be maintained in the same regions. Again,\nthese vote maps support our conclusion that the thermochemical BMSs are more sluggish than the\nambient mantle, and also more sluggish than the plumes above. We also observe that the plumes\nin the thermochemical case are likely to be more stationary than those in the purely thermal case\n(Fig. 5).\n6) Figure images are clear but captions and divisions could be made more clear. In particular there\nare multiple references to Figure 1d but there is no \"d\" labelled on figure 1. In fact Figures 1, S1,\nS2 are all in columns and should referenced as such in the paper (or labeled with a, b, c, etc).\nResponse: We have updated the captions and divisions in the revised manuscript for the figures\nboth in the main text and the supporting information to make them more clear.\nLine by line comments:\n1) Line 62: add \" \"even\" basal layer\" as it could imply 2 piles are setup to start\nThank you for this suggestion. We have modified this in the revised manuscript (line 64).\n2) Line 70. more like this like: \"like a model where a Pacific-like ... \"\nWe have modified this in the revised manuscript (line 67-68).\n3) Line 72: explain how the slabs can penetrate and push aside BMS is there no viscosity jump at\n660 ?\nWe have a viscosity contrast of 30 between the upper and lower mantles. The phase change at\n660km leads to the density increase; and when the phase change occurs within the cold slabs, they\npenetrate into the lower mantle (as in Tackley et al., 1993) and push aside BMSs.\n4) Line 80: no Figure 1d - no part d\nWe have corrected this in the revised manuscript.\n5) Line 100: typo: \"formulate an\" instead of formulate another\nWe have corrected this in the revised manuscript (line 97).\n6) Line 103: change \"compared to\" to \"as.. the thermochemical convection models).\nWe have modified this in the revised manuscript (line 99).\n7) Line 104; change \"As we\" to just \"We\"\nWe have modified this in the revised manuscript (line 100).\n8) Line 104: define what is the value of Tp(z) = ??\nThank you for pointing this out. Tp(z) is defined in the Methods (as in the original version)\n𝑇 (𝑑) = 𝑇 (𝑑)+𝑓 (𝑇 (𝑑)−𝑇 (𝑑)),\np avg p max avg\nwhere 𝑇 (𝑑), 𝑇 (𝑑) and 𝑇 (𝑑) are the horizontally averaged, maximum and minimum\navg max min\ntemperature at depth d. And we have clarified this in the revised manuscript (lines 100-101).\n9) Line 106: refer to the figure you want me to look at\nWe have modified this in the revised manuscript (line 102).\n10) Line 112: Figure 1 refer to the columns 1 and 4 you want me to look at\nWe have modified this in the revised manuscript (Figure 1).\n11) Line 113: change \"which from previous work is ..\" to \"which previous work has shown\"\nWe have modified this in the revised manuscript (line 118).\n12) Line 118 to 120: For the buoyancy ratios, list the values, a range wasn't tested just specific\nvalues (B=0.18 and 0.20 or B=0.22, 0.26, 0.30)\nWe have modified this in the revised manuscript (lines 123-124).\n13) Line 119 to 120 move \"form in 2 piles\" so it reads \"but BMSs form in 2 piles with different\ntopographies and shapes\".\nWe have modified this in the revised manuscript (line 125).\n14) Line 124-127: Change order sentences to enhance clarity. \"We vary the chemical viscosity\ncontrasts of the BMSs from 1-100 to test the influence of chemical viscosity on the pattern of\nBMSs (sup Fig 2). Tomography indicated that LLSVPs contain more bridgmanite than the ambient\nmantle. MgSiO3 bridmangite is more viscous than the MgO periclase by 10-1000 times. With\nincreasing chemical viscosity ..... \"\nThank you for this suggestion, and we have rearranged the order of these sentences, as suggested,\nfor clarity (lines 130-133).\n15) Line 143 \"are in this range\" list the range here\nWe have listed the range (between 10 MPa and 200 MPa) in the revised manuscript (line 150-151).\n16) Line 170 change \"In consideration of \" to \"Considering \"\nChanged (line 177).\n17) Line 175: change \"in compensation of \" to \"after\" Pangaea\nChanged (line 183).\n18) Line 178: change \"Considering also the less\" to \"Considering the lower\"\nChanged (line 184)."}, {"round": 2, "reviewer_comments": {"reviewer_3": "The authors have provided evidence and new analysis/figures to address the concerns of reviewers 1\nand 3, but I would like to see a few more issues addressed or clarified in the paper before recommending\nfor publication.\nIn particular, I believe the authors failed to address major comment #1 by reviewer 1, about the\ninfluence of plate reconstructions on the results.\nIn Lines 48-50, the authors states that plate reconstructions create very mobile BMSs but they do not\ninclude the plate reconstructions in their models. If a plate reconstruction model is not possible given\nthe long time scale they wish to evolve the models, then some clarity about this omission is needed in\nthe text. It should be addressed that this could change the results at least over the last 200 Myr when\nplate reconstructions are more accurate. Specifically, lines 57 states that the authors test the stability of\nthe BMSs over a shorter period of 100 Million years, but surely plate reconstructions would influence\nthese results.\nThere needs to be some more details about when stability of the BMSs was established in the models\n(was it after 1billion years of model run time, was it closer to 3.5 billion years?). Is the BMSs stability\nalways for the last 630 Myr and was there any stability in the piles before that? How long does it take\nafter the subduction girdle is formed to push the BMSs into two piles?\nMajor comment #3 from reviewer 1 asks the authors to address the viscosity in the BMSs by looking at a\nsmaller chemical viscosity contrast. I commend the authors on providing us with the results of such\nmodel runs. However, I believe the results presented in Reponse figures 3 and 4 emphasize the\nimportance this unknown viscosity of the LLSVPs plays in the stability question. Given these new runs,\nthe results of slower BMSs (compared to ambient mantle) seems contingent on them having larger\nviscosity contrast. In particular, there is some tradeoff between the chemical viscosity of LLSVPs allowing\nfor entrainment and the yield stress altering the number of piles and entrainment of material. To\ninvestigate this fully would require more model runs and analysis which is beyond the focus of this paper.\nI suggest instead, the authors include qualifications to the viscosity assumption that must be present to\nuphold these results (rather than run an entirely new suite of models).\nReviewer #3 comment 1 asks for clarity or more evidence about the stability of these two BMSs. The\nauthors included Supplemental Figure 7 but this could be improved. There is much overlap between the\ndifferent yield stress runs and it is difficult to discern the values of the 100MPa (their preferred model)\nover the entire run time. I suggest using different line types to help provide some clarity or ensuring the\nappear to show that degree 2 convection is stable over long periods for any of the models, but this could\nbe overlap in lines in the figure. Related to this figure is Supplemental Figure 1 which shows the velocity\nin the BMSs for the entire 4.5 Gyr model run. Is this velocity for the BMSs being calculated in the large\nstable piles only and after they formed? Is the RMS velocity being calculated for all of that material at\neach time step, or in the center of the pile? The stability of the velocity in the BMSs for the entire model\nrun seems to contrast the large variation in convection degree from Supp Fig 7 and the amount of BMS\nmaterial left in piles as show in figure 2.\nLine 102 Add \"Figure 4\" so it reads \"Thus, comparing Figure 4 to Figure 3\" to help with clarity\nLine 107 refers to \"The main parts of BMS\". How much area is the \"main part\" versus the margins that\ndeform? Is it just the center of the BMS that remain stable or most? Is there a way to quantify how much\nis remaining stable?\nLine 144 compares BMS velocity to \"surface plate motions\" but surface plate motions are not mentioned\nelsewhere in the paper. Either the authors need to include some surface plate motion details or they\nmean to say \"ambient mantle motion at the same depth\".\nFigure 2 needs some clarification either in the text or caption. It would help to highlight (with box/circle\nor color) the main model runs that are deemed the preferred models.\nDoes the amount of BMSs material and thus area being covered decrease or change over time in the\ndifferent models? Are the piles higher for some models or does smaller area imply more entrainment?\nIs this the final value in the model run or averaged area being covered during the \"stable\" time?\nFigure 4: Do the plumes not rise above 700-800km depth?\nSupplemental Figure 1: Include in caption that the solid lines are BMSs and dashed lines are bulk mantles\nto enhance clarity. It appears the velocities are fairly constant for the entire 4 billion year. Why is that? In\nparticular, shouldn't the BMSs have higher velocities earlier on when more entrainment is happening\nand then stable velocities when the piles appear stable? It also seems like the bulk mantle velocities are\nconverging for all model runs, whereas the BMS velocities are increasing slightly towards the end of the\nmodel run.\nThe authors present interesting work and a bit more clarity and revision would improve the paper to be\nready for publication.", "reviewer_4": "See comments in attached Word document.\nComments to Author\nThe composition of the two continent-sized large low velocity provinces (LLVPs), that are consistently\nimaged by seismic tomographic inversions atop the core-mantle boundary, has been a topic of debate\nsince their discovery four decades ago. This study aims to address this controversy using a suite of\nbillion-year-duration global mantle convection models that include plastic yielding to reproduce plate-\nlike behaviour. The authors find that forming two antipodal LLVPs (or Basal Mantle Structures [BMSs])\nat the present-day requires intermediate values of mantle yield strength (~100 MPa), as well as denser\nand (~+100 kg m-3) and slightly more viscous material within the BMSs. The result is that the BMSs\nmove a factor of ~4 slower than their surroundings, producing relative stability over 100-Myr\ntimescales.\nI am evaluating this paper after it has already undergone a thorough review, with specific instructions\nto comment on the degree to which the main concerns raised by Reviewers 1 and 3 have been\nresolved (i.e., the suggested need to use plate reconstructions to set surface boundary conditions and\nthe level of support given for the main conclusions regarding LLVP mobility). The authors have\nsubstantially revised and augmented their manuscript and overall, I think they have carefully\nconsidered and more than adequately addressed these key criticisms (see below), as well as the\nreviewers’ other concerns. In my view, this is an interesting piece of work that provides further\nindependent evidence to support the view that LLVPs are chemically distinct from the ambient mantle\nand move more sluggishly. I have some relatively minor additional concerns and suggestions I would\nlike the authors to address, but I would be happy to see this manuscript published in once these are\ndealt with.\nComments on Response to Key Concerns Raised in Previous Reviews\n1) Suggested need to use plate reconstructions as a surface boundary condition (Reviewer 1)\nI disagree with Reviewer 1’s contention that adopting a plate-reconstruction-based approach is\n“crucial” to addressing the question of LLVP stability in general. I would agree with them if the authors\nwere aiming to accurately reconstruct LLVP motion over the past 200-300 Myrs and link it directly to\nthe geological record. However, this work—as I read it—aims only to place bounds on the length of\ntime over which a degree-two LLVP structure, like that presently observed, can remain relatively\nstationary. That is, they are testing whether relative stationarity over a ~250 Myr period (as advocated\nby some) is possible given the geodynamic conditions required to form such structures in the first\nplace, not whether such limited lateral movement actually occurred in the last 250 Ma of Earth history.\nI therefore agree with the authors that their work provides a useful counterpoint to previous plate\nreconstruction-based approaches for investigating LLVP mobility. As highlighted by the authors in their\nrebuttal, the application of plate reconstructions also has its downsides, since reconstructions are\nlikely to be highly inaccurate prior to ~150 Ma due to a lack of oceanic constraints and, even at younger\nages, imposed surface forcing may not be dynamically consistent with the underlying flow field,\nintroducing artefacts. Both methods (fully dynamic and plate-reconstruction-based) have their\nadvantages and shortcomings, but when simulating LLVP formation and mobility on very long\ntimescales (>250 Ma)—as is the focus here—it is reasonable to argue that the fully dynamic approach\nmay give comparably reliable results. In any case, their results are broadly compatible with certain of\nthose obtained from plate-reconstruction-based techniques, with substantial, albeit slow, reshaping\nof LLVPs over time (see changes in BMS location in Figure S5 and variation in the dominant degree of\nthe lower mantle temperature power spectrum in Figure S7).\n2) Disputed level of support for the conclusions regarding LLVP mobility (Reviewers 1 & 3)\nReviewer 1’s other qualms about the level of support for the conclusions are, I think, mainly the result\nof certain misunderstandings. These have now been clarified by the authors, and I feel that these\nprevious misgivings have now been demonstrated to be unfounded. First, the authors have always\nbeen comparing LLVP and ambient mantle velocities at the same depth (i.e., not comparing bulk values\nfor the whole mantle). Secondly, cases where LLVP viscosity is lower than that of the ambient mantle\nhave already been tested by the authors, since the high temperatures obtained in these regions,\ncoupled with strongly temperature-dependent viscosity, outweighs the effect of the imposed\nchemical viscosity contrast when it is set to low values. Finally, there appear to be no problems with\nthe dimensionalisation of the calculations, and the choice of threshold for identifying plumes is shown\nto have negligible impact on the conclusions.\nReviewer 3’s comments centre on the degree to which LLVPs (or BMSs) have been demonstrated to\nremain relatively stationary over 100-Myr periods. In response, the authors have added much needed\nmethodological detail regarding how horizontal velocities inside the BMSs are calculated, as well as\nseveral new figures to support their interpretation (e.g., Figure 5). Overall, these clearly show that,\nwhile their margins deform, the bulk of the modelled BMSs move relatively little over a 630 Myr\nperiod, addressing Reviewer 3’s point. That said, I do have some queries regarding the vote maps used\nto conduct the thermal versus thermochemical plume mobility comparison (see below).\nFurther Comments\n1) Figure 1: What time in the model run is plotted here? The final timestep in each case? Are they\nrandom snapshots (and if so, how are they chosen)? This information should be made clear in the\ncaption.\n2) Figure 2: What time in the model run is plotted here? Or is the area of BMS and 1 pile/2 pile\ncategorisation being calculated based on average behaviour? This information needs to be made\nclear in the caption.\n3) Figure 5: In the thermal and thermochemical plume vote map comparison in the last two columns\nof this figure (maps and histograms), it appears as if there is a larger overall proportion of plumes\nthrough time in the thermochemical case; i.e., if you integrated the voting values over the surface\nof the 300 km depth slice in each case it does not look like they would give the same total, a\nconclusion which seems to be confirmed by the area under the histograms being different. It is\nalso not clear to me, even after accounting for the apparently smaller number of plumes present\nover the 630 Myr time period, that those in the thermal models are more laterally mobile than\ntheir thermochemical counterparts. Instead, it looks like regions with non-zero plume voting\nvalues are distributed over a smaller area in the thermal case, and that the lower maximum voting\nvalues (generally <0.2) are a result of individual plumes being shorter-lived rather than more\nlaterally mobile.\nCould it be that the apparently larger time-integrated proportion and greater spatial coverage of\nplumes in the thermochemical models is simply an illusion caused by the adopted colour scheme\nand histogram y-axis? What I mean is, could it be that thermal plumes are so mobile that many\npixels have voting values that are non-zero yet so small that the adopted colourbar makes them\nappear as if they are equal to zero? If that is the case, the vote map colourbar should be modified\nto differentiate between a voting value of zero (i.e., plume not present at a given location at any\ntime in the 630 Myr model time period) and a very small value (plume locally present for at least\none time step). In this case, the y-axis range on the relevant histogram would also need to be\nextended to highlight the larger number of near-zero values in the thermal case and to\ndemonstrate that the area under the histogram remains the same as that in the\nthermochemical case, despite the distribution being more positively skewed.\nIf there really is a greater proportion of plume material through time in the thermochemical\nmodels, then this result requires some explanation, especially given that the definition applied\nhere uses a relative rather than absolute plume temperature threshold. As a result, despite the\nextra radiogenic heating from the BMSs (~7 times ambient) in the thermochemical models, it’s\nnot clear that there should be a greater number of features classified as plumes since higher\naverage temperatures due to elevated radiogenic heating values won’t inherently lead to more\nplume material. Similarly, if the white regions in the vote maps really do have voting values equal\nto zero then, as stated above, it looks to me like there are fewer non-zero values in the\nthermochemical case, indicating that these plumes are potentially more laterally mobile, not less.\nThat said, the higher maximum voting values do suggest that individual plumes are less\nintermittent in the thermochemical case, which if confirmed, is an interesting result.\nIf the apparently smaller time-integrated proportion of plumes in the thermal model is a real\nfeature, then a physical explanation should be given, and it should be clarified that rather than\nnecessarily being less laterally mobile, the thermochemical model plumes are longer-lived. To\ndetermine relative lateral mobility, average horizontal velocity of the plume regions over time\nneeds to be assessed instead. Alternatively, if the apparently smaller time-integrated proportion\nof plumes it is a plotting artefact, then this issue should be rectified through modification of the\nmap colourbar and histogram y-axis, as suggested above.\n4) Line 31: I would suggest adding references to Lu et al., 2020, G-Cubed and/or Richards et al.,\n2023, EPSL here, as these studies recently showed that LLVPs are not purely thermal (i.e., their\nbasal portions are chemically anomalous) whilst reconciling a wider range of geophysical\nconstraints than the currently cited papers, making them more convincing in my view. That said,\nI must declare a potential bias/conflict here as I am an author of one of the suggested papers.\n5) Lines 75-76: Do you mean “can maintain degree-2 structure for about 630 Myr” or “can maintain\ndegree-2 structure for up to 630 Myr”? It’s not clear to me whether the focus here is on the\nmaximum length of time for which piles are stable or the average length of time. Both would be\nhelpful to state somewhere.\n6) Lines 88–90: Check grammar and sentence construction, currently very hard to parse this\nsentence.\n7) Lines 111–113: I don’t see how the referenced vote maps (Figure 5) confirm that thermochemical\nLLVPs move more sluggishly than thermal LLVPs since thermal LLVPs are not plotted, only thermal\nplumes at 300 km depth. This sentence should either be removed or better justified (e.g., by\nshowing the plume vote map at 2800 km depth to illustrate mobility of the thermal LLVP over\nlast 630 Myrs of the model run).\n8) Paragraph starting on Line 115: You should mention somewhere in this paragraph that yield\nstress and chemical viscosity contrast are being held constant (100 MPa and 10) in these variable\nbuoyancy number runs.\n9) Lines 133–135: I don’t fully understand why the increases in BMS topography and decreases in\nBMS area that occur for large chemical viscosity contrasts are the result of decreasing heat loss\nefficiency. Can this argument be fleshed out a little more?\n10) Line 154–155: Not clear what kind of “end-member” these models are argued to represent from\na yield stress point of view, i.e., are they a low-strength or high-strength end-member? They have\nlow yield stress compared to some estimates for anhydrous oceanic lithosphere, but high yield\nstress compared to others (e.g., estimates for continental lithosphere and oceanic lithosphere at\nplate boundaries), so it is not clear which extreme is being argued for here. Can this point be\nclarified?\n11) Line 177–179: Not clear how the vote maps show that it’s easier for slabs to move the slab\ngraveyard surrounding the BMSs rather than the BMSs themselves. This point was picked up on\nby Reviewer 1 too. Comparisons of the horizontal velocity of the slabs in the lower mantle\ncompared to those of the BMSs would help to substantiate this claim.\n12) Lines 184–187: I don’t really follow the argument that “the present antipodal LLSVPs should have\nbeen cut by a subduction girdle one time from the age of their formation more easily than there\nafter”. Why does the ease with which LLVPs can be cut by a subduction girdle change over time?\n13) Lines 187–189: Struggling to understand what is being proposed in this sentence and how it links\nwith the previous one. Is the idea that the girdle cuts a single BMS into two and these two BMSs\nare then shifted to their present antipodal configuration by continued slab push over the course\nof a supercontinent cycle? Or is the argument that two distinct BMSs exist prior to girdle\nformation but are not pushed into their present configuration until they interact with it? This\nsentence could do with rewording to clarify these points.\n14) Line 189–191: Ending the paper with these unresolved issues seems a slightly negative way to\nfinish. I would suggest concluding by reiterating what has been firmly established by this study\nand its implications after the need for further study to resolve these open questions has been\noutlined.\n15) Line 224: Why is 173 km set for the initial BMS layer thickness? This value stands out as being\nsurprisingly specific so a short justification would help here. I assume this choice is related to the\n3% volume fraction being assumed (Line 239), but even so, some reasoning should be given for\nwhy these numbers are adopted.\n16) Worth asking a mother-tongue English speaker to do a final language check. I found nothing that\naffected my comprehension of the text (beyond the sentences specifically highlighted above),\nbut there are several instances of incorrect/unusual syntax and word choice errors. For example,\nthe word order should be flipped for “volume infinitesimal” and “area infinitesimal” in lines 260\nand 264; “No prescribed plate boundaries are determined…” should be “No prescribed plate\nboundaries are imposed…” or “No plate boundaries are prescribed…” in lines 217–218).", "reviewer_author_response": "We thank again the reviewers and the editor for their valuable contribution to the article.\nThis is our response to the reviews of manuscript NCOMMS-23-28774A, entitled\n“Sluggish basal mantle structures support their long-lived stability”. Reviewer\ncomments are shown in black and our responses are in blue. Note that the line numbers\ncorrespond to the clean version without track changes.\nResponse to Reviewer #3:\nThe authors have provided evidence and new analysis/figures to address the concerns\nof reviewers 1 and 3, but I would like to see a few more issues addressed or clarified in\nthe paper before recommending for publication.\nReviewer #4.\n2) There needs to be some more details about when stability of the BMSs was\nestablished in the models (was it after 1 billion years of model run time, was it closer\nto 3.5 billion years?). Is the BMSs stability always for the last 630 Myr and was there\nany stability in the piles before that? How long does it take after the subduction girdle\nis formed to push the BMSs into two piles?\nReviewer 1’s other qualms about the level of support for the conclusions are, I think,\nmainly the result of certain misunderstandings. These have now been clarified by the\nauthors, and I feel that these previous misgivings have now been demonstrated to be\nunfounded. First, the authors have always been comparing LLVP and ambient mantle\nvelocities at the same depth (i.e., not comparing bulk values for the whole mantle).\nSecondly, cases where LLVP viscosity is lower than that of the ambient mantle have\nalready been tested by the authors, since the high temperatures obtained in these regions,\ncoupled with strongly temperature-dependent viscosity, outweighs the effect of the\nimposed chemical viscosity contrast when it is set to low values. Finally, there appear\nto be no problems with the dimensionalisation of the calculations, and the choice of\nthreshold for identifying plumes is shown to have negligible impact on the conclusions.\nReviewer 3’s comments centre on the degree to which LLVPs (or BMSs) have been\ndemonstrated to remain relatively stationary over 100-Myr periods. In response, the\nauthors have added much needed methodological detail regarding how horizontal\nvelocities inside the BMSs are calculated, as well as several new figures to support their\ninterpretation (e.g., Figure 5). Overall, these clearly show that, while their margins\ndeform, the bulk of the modelled BMSs move relatively little over a 630 Myr period,\naddressing Reviewer 3’s point. That said, I do have some queries regarding the vote\nmaps used to conduct the thermal versus thermochemical plume mobility comparison\n(see below)."}, "author_response": ""}]

Decision letter