A Synoptic View of Mantle Plume Shapes Enabled by Virtual Reality

The shapes of mantle plumes are sensitive to mantle viscosity, density structure, and flow patterns. Increasingly, global tomographic models reveal broad plume conduits in the lower mantle and highly tilting conduits in the mid and upper mantle. Previous studies mostly relied on 2D slices to analyze plume shapes, but fully investigating the complexity of 3D plume structures requires more effective visualization methods. Here, we use immersive headset‐based virtual reality (VR) to visualize the full‐waveform global tomographic models SEMUCB‐WM1 and GLAD‐M25. We develop criteria for the identification of plume conduits based on the relationship between the plume excess temperature and the VS anomaly (δVS). We trace 20 major plume conduits, measure the offsets of the conduits in azimuth and distance with respect to the hotspots, calculate the tilt angle, and evaluate the δVS along all traced conduits. We compare our traced conduits with the conduits predicted by global mantle convection models and vertical conduits. The wavespeed variations along conduits traced from each tomographic model are slower than modeled or vertical conduits, regardless of which tomographic model they are evaluated in. The shapes of traced conduits tend to differ greatly from modeled conduits. Plume ponding and the emergence of secondary plumes, which could result from a combination of compositional variations, phase transitions, small‐scale convection, and variations in viscosity, can contribute to the complex observed plume shapes. The variation of δVS along the traced conduits and complex plume shapes suggest a thermochemical origin of many plumes.


Introduction
Deep mantle plumes originating from the Core Mantle Boundary (CMB) are thought to have a broad head, which generates Large Igneous Provinces (LIPs), and a narrower tail, which forms long-lived hotspots (e.g.Richards et al., 1989).The geochemical diversity of hotspot lavas, which are also known as Ocean Island Basalts (OIBs), reflects the entrainment and transport of different mantle materials by ascending plumes.Hence, understanding the shapes of mantle plumes is important for linking the rock record with deep mantle structures, including the Large Low Shear Velocity Provinces (LLSVPs) and Ultra Low Velocity Zones (ULVZs).Plume shape is influenced by the global pattern of mantle circulation as well as the intrinsic buoyancy and viscosity variations within ascending plumes and the ambient mantle.Seismic tomography is the only geophysical method that currently resolves plume-scale features at all mantle depths.Tomographic models shape our understanding of mantle plumes and naturally become a constraint on numerical models that aim to understand their structure and evolution.These comparisons bridge our theoretical models to tomographic images of mantle plumes and help advance our understanding of the physical and chemical properties of mantle plumes.Here we analyze the shapes of mantle plumes using immersive 3D visualization based on two recent global tomographic models and consider the implications of plume shape for the pattern of global mantle circulation and the variation of mantle viscosity.
Mantle plumes that rise to the surface have previously been described conceptually as primary and secondary plumes (Courtillot et al., 2003) on the basis of their buoyancy fluxes, upper mantle seismic signature, and the isotopic variations in OIBs.Primary plumes rise directly from the CMB, whereas secondary plumes rise from the superswells or broad primary plumes that pond below the upper mantle.State-of-art global tomographic models show patterns of slow shear velocity (Vs) resembling both types of plumes, although the plume shapes revealed by tomographic models have more complexities than what is proposed by the schematic plume model of Courtillot et al. (2003).
There has been considerable debate about whether hotspots are preferentially located at the edges of the Pacific and African LLSVPs (Torsvik et al., 2006;Steinberger & Torsvik, 2012) or whether they are associated with the LLSVP edges and interiors (Austermann et al., 2014;Davies et al., 2015;Doubrovine et al., 2016).These two hypotheses have different geodynamics implications: whether plumes rise from the edge of the pile-like LLSVPs (Tan et al., 2011;Hassan et al., 2015), or the LLSVPs are cluster of plumes (Davaille & Romanowicz, 2020).
Two complementary approaches have been taken to understand the evolution of mantle plumes.First, some numerical models of thermal and thermochemical plume ascent focus on idealized plumes and incorporate a high degree of physical realism at the expense of describing the geologic context of specific plumes within Earth's mantle (Dannberg & Sobolev, 2015;H. Liu & Leng, 2020).A second class of numerical models focuses on the influence of the global mantle flow associated with Earth's tectonic history on plume conduits at the expense of a complete treatment of mantle rheology, phase transitions, and plume buoyancy (e.g., Steinberger & O'Connell, 1998).For the first class of studies, the comparison between the shape of the observed and the modeled plume conduits is only qualitative because idealized models do not attempt to reproduce the detailed dynamics of specific plumes.The second class of models does make testable predictions of plume shape that can be qualitatively and quantitatively compared with plumes resolved in tomographic models but only in terms of the wavespeed variations (Boschi et al., 2007).
Plume shapes depend on both the inherent properties of a plume and the surrounding mantle conditions, so they provide information about the composition and dynamics of plume and mantle.For example, the amount and direction of shear of plume conduits reflect the large-scale mantle flow pattern.Changes in the conduit radius could indicate the viscosity variations across the mantle.The stagnation of plumes helps to reveal the influence of the pressure-induced phase transitions on mantle convection.It is crucial to measure the shapes of plume conduits quantitatively to make more appropriate and meaningful connections between numerical models and tomographic observations.Measuring plume shapes from tomographic models requires effective visualization of what are three-dimensional (3D) datasets, but most approaches to their visualization have involved two-dimensional (2D) slicing or the rendering of isosurfaces (surfaces defined by a constant value) on a 2D medium such as a computer screen or a paper (French & Romanowicz, 2015;Tsekhmistrenko et al., 2021;Celli et al., 2021).The understanding and insight gained from 2D visualizations of 3D data may be different than that gained through immersive 3D visualization.For example, the 2D cross-section of a plume cluster associated with the Pacific LLSVP seems to imply that the conduits of plume Samoa, and Tahiti are not resolved above 660 km depth (Figure 1a).However, the conduits of these plumes extend out of the vertical cross-section plane, as shown in Figure 1b.Selecting an isosurface with a specific negative δV S to represent the boundary of a plume reveals plume shapes better than 2D cross sections and allows us to "see through" the non-negative δV S that obscures our view.However, these approaches may fail if the shape of a plume is best represented by different isosurface values at different depths or when many plumes are clustered.In the first case, visualizing plumes requires observing the structures of many different δV S isosurfaces simultaneously.In the second case, the isosurfaces representing boundaries of conduits usually obscure each other, making it tricky to identify an individual conduit if the observer is outside the cluster.This is the scenario for the plumes feeding Pitcairn, Macdonald, Marquesas, Tahiti, Samoa, and Easter, which are located close together within the Pacific LLSVP (Figure 2).
Visualizing seismic tomographic models in a virtual reality (VR) environment can help to overcome these barriers.Immersive visualization allows an observer to explore mantle structures from within and view them quickly from arbitrary vantage points.Immersive 3D visualization is not new in geoscience research but has not seen widespread adoption due to the lack of commodity VR hardware and related software.Previously, the usage of VR environments centered on large, immobile, and expensive "cave" environments (e.g., Billen et al., 2008).As VR headset devices have become more preva-lent, immersive 3D visualization is becoming more accessible due to its lower cost and greater portability, presenting the potential to enable new discoveries.
The remainder of the paper is structured as follows.We establish a quantitative procedure to define mantle plume conduits and discuss the advantages and limitations of our conduit-choosing criteria.We present our traced conduits for well-resolved plumes in SEMUCB-WM1 (French & Romanowicz, 2014) and GLAD-M25 (V S ) (Lei et al., 2020) and the quantitative measurement of these conduits.We demonstrate that our traced conduits are more consistent with the distributions of slow seismic velocities than geodynamic model predictions.We discuss the implications and potential applications of this study.

Methods
The two tomographic models analyzed in this study, SEMUCB-WM1 and GLAD-M25, are state-of-art global tomographic models based on full waveform inversion (FWI).SEMUCB-WM1 inverts for 3-D variations in Voigt-average isotropic V S and radial anisotropy parameter ξ and parameterizes them radially using (continuous) cubic b-splines and laterally using spherical splines.Its starting model is SEMum2 (French et al., 2013) above 800 km and SAW24B16 (Mégnin & Romanowicz, 2000) below.The crust is approximated by a smooth anisotropic layer to account for the crustal effects on wave propagation and dispersion.GLAD-M25 inverts for the bulk sound speed and vertically and horizontally polarized V S in the mantle above 660 km.Its starting model is S362ANI (Kustowski et al., 2008) for the mantle and Crust2.0(Bassin et al., 2000) for the crust.As in the starting model S362ANI, GLAD-M25 uses a parameterization that includes first-and secondorder discontinuities in the radial direction, permitting abrupt changes in the pattern of heterogeneity across the mantle transition zone (MTZ).Both of the global tomographic models resolve broad plumes rising from the CMB to the upper mantle beneath many hotspots (French & Romanowicz, 2015;Lei et al., 2020).These enforced vertical discontinuities in GLAD-M25 could introduce artifacts to the resolved plume shapes around the MTZ, but plume structures resolved in the lower mantle should remain robust, discussed later.
We define plume conduits based on three considerations.First, we require plume conduits to be continuous pathways from the lithosphere to the CMB.Second, we require that plume conduits be slower than average across all mantle depths (i.e., having a negative δV S ).Third, we seek plume conduits for which the temperature anomaly implied by wave speed variations is consistent with petrological constraints on plume excess temperature.The third criterion may not always be satisfiable due to limitations in tomographic modeling, discussed later.
Following our criteria, we manually traced the conduits of 20 plumes (listed in Table S1), of which the buoyancy flux is larger than 1000 kg/s (Jackson et al., 2021) and are well-resolved in both SEMUCB-WM1 and GLAD-M25.We exclude the Yellowstone plume as it is only well-resolved in GLAD-M25.We include the Canary and St. Helena plumes, of which the buoyancy flux is only 800 and 500 kg/s, respectively, because similar plume shapes are clearly resolved in both tomographic models.Moreover, the OIBs associated with both hotspots display isotopic signatures supporting a deep mantle origin.
The plume conduits are traced in a headset-based immersive 3D visualization environment.We use the Valve Index VR headset and controllers and the Paraview 5.10.0 (Ahrens et al., 2005) visualization software.The identification of plume conduits was carried out using the following steps: 1.The traced conduit (TC) of each plume can be divided into an upper-mantle, a mid-mantle, and a lower-mantle part.We first identify candidate conduits (CCs) -conduit-like vertical negative δV S structures -that extend vertically across the mid mantle near each surface hotspot.There may be multiple candidate conduits for each hotspot, and we seek conduits that are closer to the hotspot's surface expression.
2. We use pipelines (control points connected by line segments) to represent the pathway of the traced conduit, where the control points are assigned every 200 km from 250 to 2450 km depth.We seek an upper-mantle TC, which connects the surface hotspot with the upper-end of the mid-mantle TC, and a lower-mantle TC, which starts from the lower-end of the mid-mantle TC.Where there is ambiguity, we prefer more vertical plume conduits.
3. After tracing the plume conduits, we validate our TCs according to two criteria.
First, the δV S along a TC should not be positive.Second, we use the plume and ambient mantle potential temperature calculated from olivine-liquid equilibria (Putirka, 2008) to estimate the excess temperature of plumes.We then calculate the profile of d(ln V S )/dT (Figure S1) assuming that the plume has a pyrolitic composition and use the profile of d(ln V S )/dT to calculate δV S corresponding to the petrologicallyestimated excess temperature at all depths for each plume that has an estimation.
δV S along the TC should be comparable to δV S converted from the petrologicallyestimated excess temperature at some depths above 1250 km.The second criterion is not hardwired because the variable resolution, parameterization, and regularization of global tomographic models can all contribute to modeled V S variations.

Results
We describe the general properties of the traced plume conduits (Figure 2), starting from describing the slowness of the traced conduits.We then describe overall trends in the amount of offset from the surface location, the tilt (measured in degrees away from the vertical) of plume conduits, and the depths at which large offsets or tilts occur.We describe the shapes of individual plume conduits in greater detail later.

Slowness along plume conduits
The δV S along conduits traced from SEMUCB-WM1 and GLAD-M25 is generally between 0% and -2%, comparable with each other (Figure 3-4).We find that plumes originating from the African LLSVP are faster than plumes stemming from the Pacific LLSVP above ∼ 1250 km depth in SEMUCB-WM1 and at all depths in GLAD-M25 (Figure 5b,c,g,h).We also evaluate the average δV S of conduits traced from SEMUCB-WM1 in GLAD-M25 as well as conduits traced from GLAD-M25 in SEMUCB-WM1 (Figure 5d, e, i, j).When plumes traced in one tomographic model are evaluated in the other tomographic model, the average δV S along TCs around the Pacific LLSVP remains negative at all depths, while it is negative only in the lower mantle for TCs around the African LLSVP.

Observed morphology
Tilt angles along the traced conduits generally remain smaller in the lower mantle (usually < 60 • ) than in the upper mantle with a few exceptions (Figure 6).For example, the Louisville and Azores plumes have a tilt angle (60 -70 • ) below 2000 km in SEMUCB-WM1.A comparison of the tilt angles of plumes (Figure 6) and the offsets of plume conduits (azimuth and distance, shown in Figure 7) shows that large tilt angles are associated with abrupt changes in offset distances and/or azimuths of TCs.Changes in offset azimuths and distances are small where the tilt is closer to vertical.The azimuth of a conduit is measured by assuming its hotspot as the origin, 0 degree at the north, and counting clock-wise.Due to the manual process of conduit tracing, the uncertainty in tilt of TCs is at least 5 • .Hence TCs with tilt less than this should be interpreted as nearly vertical.We do not report the average tilt angle of each conduit because these values do not accurately describe the shape of conduits.For example, in SEMUCB-WM1, the TC of Samoa has a similar average tilt angle (16.9 • ) to the TC of Pitcairn (16.1 • ).
However, the TC of Samoa appears to be ponded and deflected at 660 and 410 km depth, while the TC of Pitcairn tilts gently across the whole mantle.
Plume conduits traced in SEMUCB-WM1 and GLAD-M25 usually root at locations offset from their surface hotspots by 5 -10 • and most of the offset occurs in the upper mantle.A few plume conduits show larger offsets.The TCs of Galapagos, San Felix, and Tahiti root at locations offset from their surface hotspots by more than 10 • in both tomographic models (Figure 2 and 7).The offsets of conduits traced from SEMUCB-WM1 in the upper mantle can easily exceed 5 degrees (Figure 7), which converts to > 500 km offsets, while those of conduits traced from GLAD-M25 appear to be much smaller.

Paired plumes
In SEMUCB-WM1, the MacDonald and Pitcairn plumes seem to branch from the same conduit in the lower mantle and the Macdonald plume is significantly deflected at ∼ 1250 km depth (Figure 8a).The Canary and Cape Verde plumes also appear to share the same conduit from the CMB to at least ∼ 1250 km depth and branch into two conduits separated by ∼ 15 • in the upper mantle (Figure 8b).
In GLAD-M25, we identify CC with a similar shape as what is observed in SEMUCB-WM1 below the Canary and Cape Verde hotspots.We interpret Canary and Cape Verde as two adjacent plumes rising parallel to each other though this CC could be interpreted as either two separate conduits or one broad plume branching into two secondary plumes as it crosses the 660 km discontinuity.CCs of the Pitcairn and Macdonald plumes look less like those in SEMUCB-WM1.These two plumes seem to emerge from different locations at the CMB and merge into a broad plume conduit between 660 and 2000 km depth and branch again above 660 km depth.
The San Felix and Juan Fernandez plumes are another potential paired plumes.
These two plumes generally share the same CC in the mid-mantle in both tomographic models (Figure S2).We interpret it as two adjacent plumes rising parallel to each other and trace their conduits based on this interpretation.The conduit of San Felix is not resolved between 1250 and 660 km in SEMUCB-WM1 and above 660 km in GLAD-M25.
The conduit of Juan Fernandez is generally well resolved at all depths in both tomographic models.

Iceland
The Iceland plume is generally vertical in both tomographic models, but the detailed shape of the plume is different.Starting from the surface hotspot, the traced conduit from SEMUCB-WM1 is offset towards the northeast above ∼ 350 km and then offset back towards the hotspot at ∼ 660 km.The conduit remains generally vertical below 660 km and slightly tilts towards the east below ∼ 2000 km (Figure 6,7, and 8c).
Its TC from GLAD-M25 is vertical above 660 km, tilts first towards the east between 660 and 1000 km depth then towards the west between ∼ 1250 and 1500 km depth, and remains vertical below 1500 km.

Hawaii
The Hawaii plume appears to be mostly vertical in SEMUCB-WM1, while it appears to largely tilt towards the southeast in GLAD-M25.Its conduit is well resolved in SEMUCB-WM1 but not well resolved between 410 and 660 km depth in GLAD-M25 (Figure 8d).Although the TCs from SEMUCB-WM1 and GLAD-M25 are not consistent, both tomographic models resolve a similar CC between 660 and 1250 km depth below the surface hotspot location and a similar CC location at the CMB (Figure 8d).

Samoa, St Helena, Reunion, and Caroline
Similar CCs are identified in both tomographic models for the Samoa, St Helena, Reunion, and Caroline plumes.These plumes remain nearly vertical or slightly tilt in the lower mantle and tilt more heavily in the upper mantle (Figure 9a-c).We noticed that amplitudes of negative δV S along these TCs from SEMUCB-WM1 vary smoothly and reach a maximum near 660 km.Amplitudes of the negative δV S along these TCs from GLAD-M25, however, decrease abruptly above the 660 km discontinuity.These negative δV S amplitudes are larger (slower) than those of conduits traced from SEMUCB-WM1 by 0.5-1.0 % δV S below ∼ 2000 km (Figure 3 and 4).

Azores, Easter, Galapagos, Kerguelen, Marquesas, and Tahiti
We notice that for the Azores, Easter, Galapagos, Kerguelen, Marquesas, and Tahiti plumes, similar CCs are resolved in the two tomographic models but different TCs are identified (Figure 9d

Discussion
We first demonstrate the reliability of our traced conduits to justify that our TCs represent seismically slow paths through the mantle.We then compare our TCs with modeled conduits and discuss the reasons for their differences.Next, we discuss the implications for mantle and plume dynamics from our observed plume shapes and slowness along conduits.We conclude our discussion by proposing some applications of our TCs in future studies of plume dynamics.

Reliability of traced conduits
Seismic tomography is a mixed-determined inverse problem, and there exist many possible Earth structures that are equally compatible with seismic observables.The shapes of plumes could vary between different regional and global tomographic models due to different parameterization/regularization choices and different earthquake events used to constrain the tomographic models (French & Romanowicz, 2015;Wamba et al., 2021Wamba et al., , 2023)).Hence, one might question the veracity of mantle plume shapes determined on the basis of seismic tomography.Several lines of evidence suggest that the imaged and traced plume conduits are likely representative of real mantle structures.First, the slow V S structures near many hotspots are similar between the two models, suggesting that the imaged features are robust.Second, the average slowness along TCs is much greater than the average slowness along modeled or vertical conduits (Figure 5a-c, f-h).To further assess the robustness of the traced plume conduits, we evaluate the slowness along Pacific TCs obtained from SEMUCB-WM1 and GLAD-M25 in other P-and S-velocity tomographic models.We find that our Pacific TCs traced from GLAD-M25 are slower than the MCs and vertical conduits (VCs) in the lower mantle (below ∼ 660−1000 km depth) when they are evaluated in most of the other models (Figure 5i and S2 g-k).Our Pacific TCs traced from SEMUCB-WM1 are slower than the MCs and vertical conduits (VCs) but in a more restricted depth range between ∼ 1250 and 2100 km depth .(See Text S2 in the Supporting Information for more details.)This suggests that both sets of traced conduits, especially TCs from GLAD-M25, are more compatible with many other tomographic models than the modeled and vertical conduits in the mid to lower mantle.

Comparison between traced and modeled conduits
Simplified numerical models of mantle plume shapes have been used widely in geodynamics to understand the mobility of deep mantle hotspots and to establish the moving hotspot reference frames necessary for absolute plate reconstructions (e.g., Matthews et al., 2016).We compare modeled conduits (MCs) from (Steinberger & Antretter, 2006) with our traced conduits.These numerical models of plume dynamics start with a mantle buoyancy structure based on a tomographic model filtered to long wavelength.The buoyancy structure is reconstructed backwards in time through the reversal of buoyancy forces and the application of time-reversed plate reconstructions at the surface while ignoring the effects of thermal diffusion, which cannot be time-reversed due to non-uniqueness.This yields a model of long-wavelength (much longer wavelength than the widths of plumes) mantle flow in space and time.Then, initially vertical plume conduits are advected by the flow field forward in time.Previous studies demonstrated that the shapes of MCs are not very sensitive to the tomographic model used to compute the mantle flow field, the details of the plate reconstructions used, or the detailed mantle viscosity structure (Steinberger & O'Connell, 1998;Steinberger, 2000;Steinberger & Antretter, 2006;Williams et al., 2019).
The tilt angles and offsets of MCs show that most of MCs slightly tilt (tilt angle < 30 • ) below 660 km.This is likely because the deformation rate is slow due to the high viscosity of the lower mantle.Larger tilt angles (up to > 90 • ) of MCs observed above 660 km (Figure 6) are mainly due to the oscillations of the tightly spacing conduit elements in the lower-viscosity upper mantle.The offsets of modeled conduits (shown in Figure 7) show that MCs in fact tilt gently at these depths.Our TCs suggest that plumes generally slightly tilt in the lower mantle, but large tilt angles in the mid-mantle below 660 km are observed for many TCs from both tomographic models (e.g., Macdonald, Samoa, St Helena, and Tristan) (Figure 6).TCs generally have more complex shapes than MCs, especially in the mid-mantle.
Although the paths of TCs and MCs are generally not in very good agreement (Figure 2, Table S1), there are a couple of exceptions.TCs of plumes located at the edge of LLSVPs (Canary, Juan Fernandez, San Felix, St Helena, and Reunion)(Figure 7) seem to agree with their MCs better than TCs of plumes located near the center of LLSVPs.
TCs of these plumes share similar offset directions with their MCs, while the MCs have 5 -10 • more total offset distances than the TCs.These plumes have relatively simple plume shapes, that is, the offset direction of a TC does not change with depth.TCs of plumes located around the center of LLSVPs are usually vertical in the lower mantle but meander in the middle and upper mantle.Because of the physics included in the models, all MCs only have simple plume shapes (without stagnation or meandering).They are always smooth curves extending from the LLSVPs to the surface hotspots.We discuss this difference more in the next section.
The average seismic velocities of the TCs, MCs, and VCs are slower than the ambient mantle at all depths.However, TCs from SEMUCB-WM1 are up to 6 times slower than MCs and 3.7 times slower than VCs in the upper mantle, while they are 1.2-3 times slower than MCs and VCs in the lower mantle.TCs from GLAD-M25 are 1.1-3 times slower than MCs and VCs across the mantle.The average velocities of MCs are slower than the those along VCs only in the lower mantle (Figure 5a, f), which is consistent with the analysis of MCs and VCs done using older tomographic models (Boschi et al., 2007).
The δV S along MCs is often close to 0% or even positive in the upper mantle (Figure 3 and 4), while the δV S along TCs is negative in most cases.There are a few exceptions in SEMUCB-WM1 (Cape Verde and San Felix) and GLAD-M25 (Azores, Canary, Hawaii, San Felix, Tahiti, and Tristan).In these cases, no CC can be identified at some depths in the upper mantle.This may indicate that the global tomographic model does not resolve the plume conduit at these depths.It is expected that plume radius can significantly decrease as a plume rising from the more viscous lower mantle to the less viscous upper mantle (Leng & Gurnis, 2012).

Implications of the slowness along plume conduits
The excess temperature of a purely thermal plume conduit is not expected to change significantly with depth since plumes rise rapidly relative to the thermal diffusion timescale and mantle heat production is negligible on the timescale of material ascent through a plume conduit.For example, the exothermic phase transition (olivine to wadsleyite) at 410 km depth, and shear heating may be able to increase the temperature of a plume, but they are secondary effects compared with the plume's inherent excess temperature.This implies that if a mantle plume is purely thermal, the amplitude of its δV S should generally vary following the thermodynamically determined d(ln V S )/dT profile with depth.
Our observations from both tomographic models, however, show that the variation of δV S along plume conduits almost never strictly follow the d(ln V S )/dT profile, which suggests that non-thermal variations are present in plume conduits.
Non-thermal variations in mantle plumes include differences in intrinsic composition, water content, grain size, and melt fraction.At the 410 km discontinuity, the phase transition from wadsleyite to olivine may result in water release when plume materials rise and cross this boundary because wadsleyite has a higher water-bearing ability than olivine (W.Wang et al., 2019).Increasing water content can reduce V S (C.Liu et al., 2023) and may cause partial melting in this region, further reducing V S (Chantel et al., 2016).Isotopic measurements of OIBs and numerical models suggest that LLSVPs may be composed of a variety of different materials, ranging from primordial materials that get preserved at the CMB since the differentiation in early Earth's evolution (Labrosse et al., 2007;Deschamps et al., 2012) to piles of recycled oceanic crusts (Olson & Kincaid, 1991;Brandenburg & van Keken, 2007).For many of the traced conduits, we find that δV S in the lowermost mantle is slower than expected on the basis of d(ln V S )/dT .The incorporation of compositionally-distinct material within the lowermost mantle is one possible explanation for the slower than expected velocities (Figure 3 and 4).
The systematically faster plumes (in the upper-and mid-mantle) originating from the African LLSVP than those originating from the Pacific LLSVP (Figure 5b, c, g, h) are consistent with previous estimates of plume excess temperature based on upper mantle wavespeed variations (Bao et al., 2022).Y. Wang and Wen (2007) and He and Wen (2009) also show that the two LLSVPs have different shape and topology.They may indicate that the two LLSVPs have different origins, but we cannot rule out the possibility that the faster plumes from the African LLSVP are caused by different seismic data coverage between the Pacific and the Atlantic regions.

Implications of diverse plume shapes
The shape of a plume conduit depends on both the plume's properties and its interaction with its surrounding mantle.Buoyancy, which is determined by ∆ρ, the difference between the effective density of a plume and the density of its surrounding mantle (∆ρ = ρ plume − ρ mantle ), controls the behaviours of a plume as it rises.The buoyant ascent of plume material and its interaction with the large-scale mantle flow will re-sult in different plume conduit shapes.The composition of the plume, the pressure induced phase transitions, and the excess temperature (temperature difference between the potential temperature of a plume and the ambient mantle) together determine ∆ρ.
When a plume has a positive buoyancy (∆ρ < 0), it will rise, and it will start sinking when it has a negative buoyancy.When ∆ρ is close to or slightly smaller than 0, a plume could be ponded or develop a variety of complex shapes (Kumagai et al., 2008;Xiang et al., 2021).
The mantle viscosity structure and flow patterns of the ambient mantle also affect plume shapes.The mobility of a plume, that is how easily it gets deformed, is expected to be smaller in a more viscous than in a less viscous region (H.Liu & Leng, 2020).Largescale mantle flows driven by thermal convection, surface plate motion, and subduction could shear plume conduits or largely deflect the secondary plume stemming from a ponding primary plume (Steinberger, 2000;Farnetani & Samuel, 2005).
The more complex shapes of our TCs than the MCs suggest that the mantle convection models used to determine MCs may not consider all major factors affecting plume shapes, especially in the mid-mantle across and below the MTZ, where plume ponding and large tilt angles are only observed in TCs.
First, the mid-mantle below the MTZ could have significant viscosity variations (Marquardt & Miyagi, 2015;Rudolph et al., 2015;Shim et al., 2017), which indicates a more complex rheology than the numerical models' assumption that only a few deformations occur and diffusion creep is predominant at these depths (Ferreira et al., 2019).Furthermore, the transition from ringwoodite to bridgmanite at 660 km, which can lead to plume ponding at this depth, is not considered neither.In return, the numerical models lack the ability to produce plumes that are ponded and deflected at different depths due to their simplified physics, which does not consider the composition variations, phase transitions, nor a temperature-dependent or strain-rate-dependent viscosity.
Second, the mantle flow field converted from the global tomographic model (Steinberger & O'Connell, 1998) may not be accurate at a smaller scale due to our current incomplete understanding of mantle dynamics.MCs are determined based on the assumption that a plume rose to the surface vertically within a short time and left a vertical 100kilometer-radius conduit that gets passively advected by the large-scale mantle flows later.
However, this assumption is only valid if mantle plumes are purely thermal.Recent seismic tomographic models have imaged plume conduits with a radius of ∼ 500 km (French & Romanowicz, 2015) and much more complex morphology (Tsekhmistrenko et al., 2021;Celli et al., 2021;Wamba et al., 2023).Such broad plumes may not only be passively advected, but also influence the mantle flow field.Plumes with such large radius would have buoyancy fluxes that are much higher than previous estimations (Sleep, 1990;King & Adam, 2014).Together with the complex plume shapes, they suggest that many, if not all, mantle plumes are thermochemical rather than purely thermal.For example, a plume that incorporates an eclogitic component has a lower buoyancy flux and a larger radius than a purely thermal plume, which is more consistent with observations (Dannberg & Sobolev, 2015).
At ∼ 410 km depth, previous numerical models suggest that plumes with some eclogitic component will have a buoyancy barrier due to the different phase transitions that occur in pyrolitic and eclogitic materials.This buoyancy barrier can result in plume ponding and the emergence of a secondary plume (Farnetani & Samuel, 2005;Dannberg & Sobolev, 2015).It can potentially explain the ponding of Samoa, a large tilt angle, and a large change in offset distance observed in SEMUCB-WM1 at this depth (Figure 6,9a).
Large tilt angles at 660 km depth mostly reflect plume ponding, which could be caused by the combined effect of the ∼ 30-fold viscosity increase from above to below 660 km suggested by many geophysical studies (Hager, 1984;Mitrovica & Forte, 1997) as well as the endothermic phase transition from ringwoodite to bridgmanite (Faccenda & Dal Zilio, 2017).The phase transition can cause plume ponding as the hotter plume materials undergo this phase transition at a shallower depth, hindering ascent.Several scenarios may happen after a primary plume is ponded at this depth.First, the primary plume could penetrate the 660-discontinuity broadly while some plume materials are ponded.
These ponding materials become so hot that there is a significant viscosity reduction, allowing the conduit to be laterally deflected by hundreds of kilometers (Tosi & Yuen, 2011).This scenario is observed for St. Helena and Tristan in both tomographic models (Figure 6, 7, 9b and S6).
When the primary plume cannot penetrate the 660-discontinuity in the first place, significant amount of plume materials will accumulate at this depth.The ponding materials will spread like a pancake and secondary plumes can develop from anywhere above the ponding zone.As a result, the offset distance between an upper-mantle secondary plume and a lower-mantle primary plume is not large, while the offset azimuth can be irrelevant to the flow patterns (Caroline in GLAD-M25, Azores, Iceland, Reunion in both tomographic models) (Figure 6, 7, 8c, 9c, S3 c and d).They may resemble the "plumetree" model proposed in Liu and Leng (2020), which requires a thin low-viscosity layer beneath the 660 km ponding depth and a low-viscosity upper mantle to allow secondary plume(s) develop from any part of the ponding materials.
At a greater depth ∼ 1250 km, large tilt angles observed of Tahiti in both tomographic models, Hawaii in GLAD-M25, and Kerguelen in SEMUCB-WM1 (Figure 6) could arise if the viscosity is higher around this depth than in the mantle above and below it.
Owing to the higher viscosity, conduits tilt less around this depth, so the conduit above this depth could be preferentially deflected by mantle flow.Some inversions of geophysical data suggest that there exists a viscosity hump, a one-to-two-order of magnitude viscosity increase, between 800 and 1200 km depth (King & Masters, 1992;Mitrovica & Forte, 1997;Rudolph et al., 2015).Studies on mineral physics also suggest that the increasing strength of ferropericlase (Marquardt & Miyagi, 2015;Deng & Lee, 2017) and decreasing the iron-enrichment in bridgmanite (Shim et al., 2017) at the mid-mantle depth can both result in this mid-mantle viscosity hump.
Another mechanism that may produce large tilt angles at ∼ 1000−1250 km (Canary and MacDonald in SEMUCB-WM1) is plume ponding and secondary plumes emerging.This mechanism is proposed by Wamba et al. (2023) to explain alternating vertical conduits and horizontal ponding zone observed for the Reunion and Comores plumes from ∼ 1000 km depth to the top of the asthenosphere in the latest tomographic models.There is no known endothermic phase transition, which could cause plume ponding, at these depths.However, a denser mantle below ∼ 1000 km depth due to its higher basalt content (Ballmer et al., 2015) could cause plume ponding at this depth if the thermal expansion effect is not strong enough to reduce the plume effective density to be smaller than the mantle density above ∼ 1000 km (Xiang et al., 2021).Seismic observations imply a not-global discontinuity presenting at 1000 km depth (Zhang et al., 2023), which may indicate a compositional layered mantle.
Other than these various behaviours of a single plume conduit, plume merging may further complicate the observed plume shapes.For example, we identify two CCs for Galapagos in the mid-mantle that merge into one CC with < 1% δV S above 660 km in SEMUCB-WM1.It may represent that two adjacent conduits are ponded at 660 km and the ponding zones of them merge into one conduit or these two resolved CCs are caused by a lack of resolution in SEMUCB-WM1 as they are only observed in SEMUCB-WM1.The TCs of Macdonald and Pitcairn from GLAD-M25 suggest these two plumes merge in the midmantle and branch above 660 km.Merging of two adjacent plumes has been demonstrated by both lab experiment (Moses et al., 1991) and numerical models (e.g., Lewis-Merrill et al., 2022;Brunet & Yuen, 2000), and the branching of the merged conduit could be explained by secondary plumes emerging from a ponding plume.
Given all these uncertainties in our interpretations of plume dynamics from observed plume shapes, our TCs are useful for future numerical modeling.For example, idealized plume models can explore under which geodynamics setting, the observed plume shapes can be reproduced.Our TCs can also provide a better schematic model for future studies to interpret the geochemical heterogeneity of OIBs from different hotspots.For example, previous studies have tried to interpret the heterogeneous isotopic signals of OIBs from neighbouring hotspots by correlating them with the vertical projection of the hotspots onto the CMB (Huang et al., 2011;Harpp & Weis, 2020) or interpreting these isotopic signals under simplified schematic plume models (Williams et al., 2019;Cordier et al., 2021).Our TCs can provide information about potential inter-plume interactions and the ascent history of plumes, which can be critical to the interpretation of geochemical observations.

Conclusion
Broad plumes clustering around LLSVPs have been recognized from the latest global tomographic models.Our study presents a systematic analysis of the pathways of these plume conduits.We carried out an analysis of the shapes of plume conduits in an immersive headset-based virtual reality (VR) environment.The wavespeed variations along the traced conduits from SEMUCB-WM1 and GLAD-M25 generally appear to be slower than the conduits predicted by geodynamic models and vertical conduits in the mid to lower mantle depth regardless of which tomographic models they are evaluated in.The traced conduits are 1.1 -3 times slower than either modeled or vertical conduits.This suggests that our manually-traced conduits are more consistent with the locus of slow seismic velocities within the mantle than either the vertical conduits that some authors have assumed when relating surface observables to deep mantle structures or the shapes of plume conduits predicted using physically simplified geodynamic models.Moreover, our traced conduits are more consistent with the petrologically-determined excess temperature than either of the other types of conduits.
In our manually traced conduits, the total amount of offset from the surface to the deep mantle is comparable between many traced and modeled conduits (usually smaller than 10 • ), while the offset direction of traced and modeled conduits usually differ.Some traced conduits of plumes stemming from the edge of the LLSVPs (Canary, Juan Fernando, Reunion, San Felix, and St Helena) tend to be 5 -10 • less offset than their modeled conduits, but the traced and modeled conduits share similar offset directions.Our traced conduits reveal a tendency for plumes to stagnate or to be offset at mid-mantle depths (660 -1250 km), a behavior that is not captured in modeled conduits.Previous geophysical studies, mineral physics studies, and geodynamics modeling provide multiple mechanisms that could contribute to plume ponding or deflection, including the buoyancy barrier induced by phase transitions and the viscous decoupling of conduits.The large variations of V S anomaly along plume conduits and the complex observed plume shapes together suggest that many plumes are thermochemical.Our analysis of plume conduit shapes provides a dataset that can be of value across multiple disciplines including geodynamic modeling, geochemistry, and mineral physics.

SEMUCB-WM1
Modeled GLAD-M25        March 30, 2024, 11:35pm and S3-5).One of the main causes is the poor inter-model agreement above 660 km and below 2000 km.The other main cause is that the δV S of CCs with similar shapes can amplify at different depths in different tomographic models.It can result in very different interpretations of the most-reasonable conduit path.

Figure 1 .
Figure 1.(a) Cross section of Pacific plumes in SEMUCB-WM1 and the location of the crosssection on the map, and (b) the 3D image of -2%, -1.2%, and -0.75% δVS isosurfaces taken from the same region.

Figure 2 .Figure 5 .
Figure 2. Traced and modeled (Steinberger & Antretter, 2006) plume conduits in SEMUCB-WM1 (top) and in GLAD-M25 (bottom).The colorful dots represent modeled conduits, while black-white dots represent traced conduits.The green circles represent the location of hotspots.The background shows δVS at 2850 km depth.Plate motions in the spreading-aligned mantle reference frame of Becker et al. (2015) are shown with gray arrows.

Figure 6 .
Figure 6.The depth profile of tilt angle along 20 plume conduits.Blue represents the conduits modeled in Steinberger and Antretter (2006).Red represents the traced conduits in SEMUCB-WM1.Yellow represents the traced conduits in GLAD-M25.The gray line marks the 60 • angle.

Figure 7 .
Figure 7. Azimuth and offset distance of model-predicted conduits and conduits traced in SEMUCB-WM1 and GLAD-M25 with respect to hotspots.Blue represents the azimuth of a conduit at different depths.Green represents the angular offset between a conduit and its hotspot.

Figure 8 .
Figure 8. Cross section and map view of the traced conduits of a) Macdonald and Pitcairn, b) Cape Verde and Canary, c) Iceland, d) Hawaii in SEMUCB-WM1 and GLAD-M25.From top to bottom, the dash lines represent 410, 660, and 1250 km depth.

Figure 9 .Figure S1 .Figure S2 .
Figure 9. Cross section and map view of the traced conduits of a) Samoa, b) St Helena, c) Reunion, d) Easter similar to Figure 8.

Figure S3 .
Figure S3.Cross section and map view of the traced conduits of San Felix and Juan Fernandez

Figure S4 .:
Figure S4.Cross section and map view of the traced conduits of a) Caroline, b) Louisville,