Yellowstone Plume Conduit Tilt Caused by Large‐Scale Mantle Flow

Mantle plumes are hot upwellings of rock thought to originate at the core‐mantle boundary. As they rise through the mantle, their conduits may become tilted due to lateral large‐scale mantle flow. Recent tomographic images have revealed a strongly tilted plume conduit starting at the core‐mantle boundary beneath northern Baja California rising toward the Yellowstone hot spot from the southwest. Here we perform numerical computations of plumes deflected in large‐scale mantle flow with the aim of finding if realistic model parameter ranges exist that yield a good fit with the tomographically observed conduit. We restrict ourselves to models that yield reasonable results for plume conduit tilt and hot spot motion globally. These models require high viscosity ≈1023 Pa s at some lower mantle depths. For a plume head reaching the surface 17 Ma, corresponding to the start of the Columbia River Basalts, our models require rise times ≈80 Myr or longer to match the tilt of the conduit observed by tomography. We used several tomography models to determine mantle density with almost all models predicting southwestward flow in the lowermost mantle beneath the western United States. Exact details of the shape of the predicted conduit's southwesterly tilt vary, depending on the density and viscosity structures we used. In many cases we find comparatively strong tilts in two depth ranges, in the upper and lower portions of the mantle, which is also a characteristic of the tomographically observed conduit. We expect that future models may help to constrain large‐scale flow by matching these corresponding depth ranges.


Introduction: Yellowstone-A Whole-Mantle Plume?
Yellowstone National Park, with its geysers and other geothermal activity, and three large volcanic eruptions during the past ≈2 Myr, is often regarded as an intraplate hot spot caused by a deep mantle plume (Pierce & Morgan, 1992). Older silicic caldera-forming volcanism with increasing age toward the southwest is found along the Snake River Plain (Armstrong et al., 1975), as would be expected if the North American plate moved over a stationary or slowly moving hot spot (Morgan, 1971). Eruption of the Columbia River Basalts (CRB)-a Large Igneous Province (LIP)-is thought to have been initiated when the corresponding plume head (Richards et al., 1989) reached the base of the lithosphere around 17 Myr ago. Alternatively, the Yellowstone plume may be older and the CRB caused by interaction of the plume conduit with the subducting slab (Wells et al., 2014).
Despite these features, shown in Figure 1, Yellowstone is a somewhat atypical hot spot candidate: The CRB are rather small compared to other LIPs, and the largest eruption centers of the CRB are located 400 km away from the start of the Yellowstone hot spot track at the McDermitt Caldera. More importantly, its location is unusual: Whereas most hot spots (Thorne et al., 2004) classified as deep-rooted plumes, based on tomographic images (French & Romanowicz, 2015)-and corresponding LIP eruption locations (Torsvik et al., 2006) are associated with the margins of Large Low Shear Velocity Provinces (LLSVPs), the Yellowstone hot spot is located far from any LLSVP and instead is located in a region of recent subduction where an upwelling plume would have to circumvent sinking slabs. The subduction has been either at the western margin of the westward moving north American continent, resulting in the Farallon slab (Grand, 1994;Grand et al., 1997) sinking beneath North America, or at intraoceanic subduction zones further west (Sigloch & Mihalynuk, 2013;van der Meer et al., 2012). An example for a rather simple plate reconstruction, where subduction basically occurs at the western margin of the North American continent, as it moves westward, is shown in Figure 2.  Nelson and Grand (2018), using TX2016 (Lu & Grand, 2016) as the starting model. Tickmarks are every 500 km in depth. Colored circles indicate the hot spot track along the Snake River Plain, from Yellowstone (YS) to McDermitt Caldera (McD). The area covered by the Columbia River Basalts (CRB) is shown in gray. The direction of recent North America absolute plate motion is indicated by the blue arrow. Plate boundaries in the region include the Cascadia subduction zone (indicated by line with solid triangles), the Pacific-North America plate boundary with strike-slip motion (indicated by half-arrows) along the San Andreas fault, and the spreading ridge and transform faults between the Pacific and Juan de Fuca plate. (right) Cross section through the tomography model of Nelson and Grand (2018) along the plume conduit as shown in the left panel, with angular distance along the track in degrees.
With the deployment of the USArray, aperture and resolution became sufficient to image a slow conduit-shaped anomaly to the southwest of Yellowstone down to ≈900 km using standard S wave tomography techniques (Obrebski et al., 2010;Schmandt et al., 2012). Recently, Nelson and Grand (2018) were able to image this slow anomaly extending down to the core-mantle boundary using core waves, which they interpreted as a whole-mantle plume. According to their image (Figure 1), the plume conduit has a diameter of about 350 km in the deep mantle and is tilted toward the south-southwest. Its source location in the lowermost mantle is beneath northern Baja California, about 1,500 km away from Yellowstone. The imaged conduit is not uniformly tilted, but the strongest tilts appear in two depth ranges, about 2,500-2,000 and about 800-1,300 km. Continuity of the conduit from the core-mantle boundary to the surface is not completely clear from this image: The method of detecting plumes using SKS and SKKS waves applied by Nelson and Grand (2018) starts to lose sensitivity around 1,000-km depth, due to the waves' vertical raypaths in the uppermost mantle, leaving a resolution gap between ≈900and 1,300-km depths, where continuity of the conduit is only weakly seen.
In the lowermost mantle, a negative seismic velocity anomaly in the presumed plume source region is also imaged by several global tomography models, both P wave (Amaru, 2007) and S wave whole mantle models (Auer et al., 2014;Houser et al., 2008) and D ′′ models (Castle et al., 2000;Kuo et al., 2000). This could be either a separate low-velocity anomaly, or a tendril extending from the Pacific LLSVP. The latter can be inferred from the map of the SAVANI (Auer et al., 2014) tomography model in the lowermost mantle, but the connection is not seen in most other tomography models. Because a plume conduit in the lowermost mantle has only recently been imaged and other characteristics of the hot spot system are not predicted by the classical plume hypothesis (see Foulger et al., 2015), alternative scenarios for the origin of Yellowstone volcanism and the CRB related to subduction processes have been suggested. Based on tomographic images, James et al. (2011) suggest that volcanism along the hot spot track results from slab fragmentation, trench retreat, and mantle upwelling at the tip and around the truncated edges of the descending plate. Liu and Stegman (2012) propose an origin of the CRB controlled by propagating rupture of the Farallon slab, and Zhou et al. (2018) devise models where western U.S. volcanism is due to intruding oceanic mantle driven by ancient Farallon slabs. Zhou (2018) finds an anomalous mantle transition zone beneath the Yellowstone hot spot track and suggests that reversed-polarity subduction generated passive upwellings from the lower mantle, which ascended through a water-rich mantle transition zone to produce melting and age-progressive volcanism. These hypotheses would be unlikely or at least need to be modified if the plume imaged by Nelson and Grand (2018) is correct.
Our motivation for this paper is to see if the large degree and direction of tilting of the proposed Yellowstone plume imaged by Nelson and Grand (2018) is supported by what we know about mantle flow and to place some constraints on mantle dynamics, especially the rise time of plume heads. We find that our results provide additional support for the whole-mantle plume scenario, because the geometry of the conduit as imaged by tomography can be generally matched by numerical model predictions using a variety of reasonable parameters and starting models. In section 2, we introduce our numerical method for computing the shape of a plume conduit in large-scale mantle flow. Our method is similar to what is used by Steinberger and O'Connell (1998) but with an important addition of considering conduit distortion during the time when the plume head rises from the core-mantle boundary to the surface. In section 3, we explain the input used by our model-especially the mantle density anomalies and viscosity structure used to compute mantle flow. These are not direct measurements but instead are based on models-such as from seismic tomography. Given the many assumptions and approximations, it is difficult to consider formal uncertainties. Therefore, our approach is to use various input models to assess which features of model predictions arise more consistently and therefore appear to be more robust. These results will be shown in section 4, discussed in the context of our existing knowledge in sections 5 and 6, and summarized in section 7.

Methods: Computing the Shape of the Plume Conduit Distorted in Large-Scale Mantle Flow
Since our approach is largely based on a method first developed by Steinberger and O'Connell (1998), we only briefly summarize it here. Instead, we mainly focus on describing the new feature introduced here of how to consider changes in plume conduit shape while the plume head rises through the mantle. In earlier work (Smith et al., 2009;Steinberger, 2000;Steinberger & Calderwood, 2006;Steinberger & O'Connell, 1998;Steinberger et al., 2004), conduits were taken as initially vertical, corresponding to the assumption that plume heads rise very fast such that the advection of the trailing plume conduit during their rise time can be neglected. The new model development has also been applied and described in Torsvik et al. (2019). First, we compute large-scale flow following the spherical harmonic-based method of O'Connell (1979, 1981), extended to consider compressibility (Panasyuk et al., 1996), with gravity g = 10 m/s 2 . How we obtain input models of mantle density anomalies, surface plate velocities, and radial mantle viscosity structure will be described in the next section. The boundary condition at the CMB is free slip except for spherical harmonic Degree 1 toroidal flow (i.e., the net rotation component) where we use no slip such that a net rotation of the lithosphere does not lead to net rotation in the deep mantle. Mantle flow is time dependent, because we considered both time-variable plate motions and backward advection of mantle density anomalies from the computed flow field. We limit backward advection to the most recent 68 million years because for older times, the effect of disregarding thermal diffusion becomes increasingly severe. To check the effect of backward advection, we compare models of plume conduit shape made using both time-dependent flow and constant present-day flow and show that the results are quite similar; that is, the backward advection is not an essential component to understand or produce the general nature of our results.
In the second step, we compute the shape of the plume conduit through time by assigning to each point in the conduit a velocity, which is the vector sum of ambient flow and vertical buoyant rising velocity. We treat the plume conduit as a chain of points, reaching from the source depth z b near the base of the mantle (see next section for specific values depending on mantle viscosity structure) to the base of the lithosphere at z t = 100 km depth. How we compute the buoyant rising velocity of the conduit is also explained in the next section.
To model rising of the plume head, we assume for simplicity that the vertical velocity component is inversely proportional to mantle viscosity at that depth; that is, v(z) = C∕ (z) whereby (z) is the radial viscosity structure (see next section). C takes into account the density differences between the plume and mantle and also the size and geometry of the plume head. But since we do not want to make explicit assumptions on these parameters, we rather prescribe a total plume rise time T, which will be varied. C can then be Geochemistry, Geophysics, Geosystems In contrast to the plume conduit, we only consider buoyant rising, not ambient mantle flow to compute the vertical velocity component of the plume head through equation (1). Ambient flow is only used for its horizontal components of motion. This approximation is more appropriate for the plume head than for the conduit, because the rise velocity is much faster due to the large amount of hot buoyant material in the head. We use this approximation such that the total rise time T, and hence, the initiation time of the plume head can be prescribed through equation (1), rather than also finding it iteratively. The modeling procedure is sketched in Figure 3; the large-scale flow field is shown in Figure 4. The cross sections show that, through all times, the flow field can be approximately represented by a simple convection cell, as in the cartoon. However, especially close to present day, flow in the lower part of the mantle is not nearly horizontal but shows a wavelike pattern with ups and downs, which can also be seen in map view. The conduit becomes sheared between different flow directions at different depths, and tilted conduit elements may rise buoyantly, aided by large-scale upward flow. As the plume head rises faster than the conduit, intermediate conduit elements have to be introduced, such that the conduit remains connected to the head, as is the case in experiments and probably also in reality. As a limiting case for T → 0 we also consider an initially vertical conduit. The plume starting position is varied in order to closely match the present-day hot spot position at the surface. Figure 5 shows the rising speed of both the plume head (left) and the conduit (right). The width (corresponding to ≈300 km) of the line for plume head rising is chosen arbitrarily; the actual size of the plume head is probably larger. Colors for conduit rising are shown on the right to 0 Ma, since the conduit continues to connect lithosphere and D ′′ after the head has reached the lithosphere until present day. The plume head rises much faster-note the nonlinear color scale-but nevertheless remains connected to the conduit, which is newly created in the wake of the rising head. In our modeling procedure, the rising speed of the plume head is independent of the tomography model and only depends on the viscosity structure we use and a prescribed rise time, while the conduit rise speed depends on the density and mantle flow structure inferred from the tomography model, in addition to those factors. For the conduit, ambient flow velocity may dominate buoyant rising (see next section), especially in the lower mantle, and conduit rise speed varies with both depth and time. Conduit rise speed tends to increase with time as the conduit is advected southward toward a region with a stronger upward mantle flow component. The strong dependence of rise speed with depth for both conduit and plume head is largely due to the strong viscosity variation with depth, which is also shown. Depth of the plume head versus time follows the line between the colors indicating conduit rising and head rising.

Model Input: Density Anomalies, Viscosity Structure, Plate Velocities, and Plume Parameters
Density models are based on Tomography Models TX2016 (Lu & Grand, 2016), S40RTS (Ritsema et al., 2011), SAVANI (Auer et al., 2014), and SEMUCB-WM1 (French & Romanowicz, 2014). In some cases we combine the global models with a regional modification of the mantle beneath the western United States from Nelson and Grand (2018). Since this regional model also features the plume conduit, one might consider to somehow "remove" the conduit before computing large-scale flow in order to avoid circularity. However, our results will show that this regional refinement does not change results by much; we hence chose not to attempt such a modification. To further assess the robustness of our results, we include a number of other density models in the supporting information.
To convert from seismic velocity-to-density anomalies, we use the thermal conversion factor Model 2 from Steinberger and Calderwood (2006) but scaled down by a factor 0.5 in the uppermost 220 km as a simplified way to account for compositional anomalies in the lithosphere (Forte & Perry, 2000;Jordan, 1978). In some cases we have additionally scaled up the density anomalies by a factor 1.5 in order to compensate  for the underprediction of velocity anomalies by tomographic inversions due to damping. All models are provided and used in a spherical harmonic format (Becker & Boschi, 2002). Additional density anomalies due to deflection of the "410" and "660" phase boundaries are considered following Steinberger (2007) but disregarding the effect of latent heat release. The preliminary reference Earth model of Dziewoński and Anderson (1981) is used to convert relative to absolute density anomalies. Additionally, we tested scenarios with a more detailed treatment of the lithosphere, and thermochemical piles in the lowermost mantle, but find their effects on the plume to be minor (see supporting information).
If not specified otherwise, we use the preferred viscosity Model 2b from Steinberger and Calderwood (2006) (Figure 5). For comparison we also show ( Figure 8 and supporting information Figure S1) and discuss results for some other viscosity models. Surface plate velocities are adopted from Torsvik et al. (2010). In their plate reconstruction, absolute plate motions are determined separately in the Pacific and Indo-Atlantic hemispheres, based on separate sets of hot spot tracks (not including Yellowstone) and taking hot spot motion into account. Relative plate motions within these hemispheres are adopted from elsewhere, mainly based on seafloor magnetic anomalies.
We compute the distortion of plume conduits with time-dependent flow. In the supporting information, we also show results with constant present-day flow for comparison. A drawback of using time-dependent flow is that backward advection causes the model density structure in the mantle to become increasingly layered back in time, which artificially lowers flow speeds. However, the similarity of results with constant and time-dependent flow shows that this most likely does not affect our results in a major way. To determine the buoyant rising speed of conduit elements, we mostly (unless mentioned otherwise) follow the formulas given in Steinberger et al. (2004), with anomalous mass flux 1.5 · 10 3 kg/s (Sleep, 1990). Essentially, rising speed is computed with a Stokes formula that has been modified to account for the difference in shape (conduit vs. sphere). The conduit varies between radius 57 km at a depth where the ambient mantle viscosity is lowest and radius 110 km where it is highest. This variation occurs because we assume that higher ambient mantle viscosity also corresponds to higher viscosity inside the conduit at the same depth, and a larger conduit radius is required to achieve the same mass flux for higher viscosity (Steinberger et al., 2004). An alternative model of rising speed is considered in the supporting information.
The present-day surface position to be matched is 110.5 • W, 44.6 • N, located in the northeastern part of the third Yellowstone caldera. The conduit is modeled between source depth and the base of the lithosphere. Which source depth is chosen depends on the radial viscosity model: If there is a low-viscosity layer at the base of the mantle, we choose the top of this layer (depth 2,620 km in the reference case), corresponding to the concept that the plume is being fed by this low-viscosity layer. Without such a low-viscosity layer, we rather choose a source depth 2,800 km near the base of the mantle. Figure 6 shows results for different plume head rise times. For the left column, the computation is initiated with a vertical conduit, as in our previous work. However, this case yields a conduit tilt much smaller than inferred from tomography. We have therefore developed the above-described model extension: If the plume head takes some time to rise, the conduit will develop a tilt, while the plume head is rising as opposed to initially being completely vertical. As expected, predicted plume conduit tilt becomes larger the longer the plume head takes to rise through the mantle, and rise times required to approximately match tomography Figure 6. Hot spot tracks (row 1), and map views (row 2) and side view (row 3) of plume conduits computed for different rise times of the plume head. For column 1, the computation is initiated with a vertical conduit (corresponding to the limit of very fast plume head rising); rise times in columns 2, 3, and 4 are 40, 80, and 120 Myr; corresponding times when the plume head starts rising are indicated in each column. In all cases, the plume head reaches the surface at 17 Ma. The density model is based on the global model TX2016 (Lu & Grand, 2016) with velocity-to-density conversion scaled up by a factor 1.5 using the same viscosity model as in Figure 5. Brown lines in row 1 are for a fixed hot spot; black lines consider hot spot motion; tickmarks are every 5 Myr. Red lines in row 3 are projections of the plume conduit onto the same plane used by Nelson and Grand (2018)  are quite long, around 120 Ma (Figure 6, panel 2d), with this particular mantle flow model. In addition to the overall amount and direction of the conduit's tilt, there are further similarities to the conduit shape inferred from tomography, both here and (even more so) in some of the cases discussed below: In two depth intervals in the upper and lower part of the mantle the tilt is stronger than elsewhere. However, the stronger tilts do not occur at exactly the same depths as what is inferred from tomography. We also see that time dependence of flow does not have a large influence on the results, as the overall structure of large-scale mantle flow remains similar through the modeling time (supporting information Figure S1).

Results: Plume Conduit Rising From South-Southwest and Hot Spot Track From Southwest
In the top row of Figure 6, the predicted hot spot tracks 0-17 Ma for both the computed motion of the hot spot (i.e., where the plume conduit reaches the lithosphere) and for a fixed hot spot, and for North America absolute plate motion from Doubrovine et al. (2012) are shown. The observed hot spot track along the Snake River plain begins near the Oregon-Idaho-Nevada state corner around 15-16 Ma (Rytuba & McKee, 1984, McDermitt Caldera), and there is a mismatch with the predicted fixed hot spot track (brown line in top row), which begins further east. Because the plume conduit in the upper mantle is tilted over by eastward mantle flow (see also the supporting information Movie S1), an eastward hot spot motion is computed, yielding a better fit than assuming hot spot fixity. This result is independent of plume head rise time, as in all cases, the plume head rises very fast through the upper mantle, causing an initially almost vertical conduit in the upper mantle. In contrast, the trailing conduit rises much more slowly causing it to become tilted due to upper mantle flow, resulting in eastward movement of the hot spot.
Alternatively, we also consider an older plume age with initially vertical conduit instead of longer plume head rising time. Results in supporting information Figure S1 show that the amount of tilt that develops increases with age, as there is more time available to shear the conduit. For a plume age around 77 Ma (row 4) both direction and amount of tilt inferred from tomography are approximately matched. However, in this case, no eastward hot spot motion during the last 17 Myr is computed, as the conduit is already tilted over at 17 Ma. Hence, the observed hot spot track is fit less well. Also, in some cases the present-day location of Yellowstone cannot be exactly matched; that is, our iteration procedure sometimes cannot find a suitable initial plume location to match the present-day hot spot location, especially if the plume conduit gets strongly sheared. Figure 7 shows results for density based on various tomography models. Again, most of them share the same overall characteristic that the plume conduit rises from the southwest, with one notable exception for SAVANI (column c). The different tilt in this case corresponds to a different regional flow field in the lowermost mantle: This is the only model we tested where the flow is not predominantly toward the southwest. Globally, the flow field in the lowermost mantle computed using SAVANI looks similar to other tomography models, with flow toward large-scale upwellings above the two LLSVPs and away from downwellings driven by subducted slabs. However, for SAVANI the southwestward flow toward the Pacific LLSVP does not extend as far toward the northeast as for other models because in contrast to the other models there is not a sufficient quantity of downgoing slabs beneath the northern and eastern parts of North America that would push material southwestward in the lowermost mantle. Model TX2016 is upscaled by a factor 1.5 to compensate for lower amplitude due to damping, which would otherwise lower flow speeds. Differences in the detailed conduit shape correspond to variations in smaller-scale structure between tomography models.
Supporting information Figure S2 also shows results for several more density models. These results are similar to those shown in Figure 6: In particular, when a more realistic density model for the lithosphere is used, results remain almost indistinguishable compared to using the simple downscaling in the upper 220 km as for the other models. Hence, we regard this simplified procedure as appropriate, as details of lithospheric density hardly matter. Also, the effect of excess density in thermochemical piles on the Yellowstone plume conduit is minimal due to the distance from Yellowstone. Comparison of results in Figures 7 and supporting information Figure S2 based on identical tomography models shows that there is a trade-off between velocity-density scaling and plume head rise time, as an increase in predicted tilt can be achieved by increasing either of them. Results remain similar if the tomography model is upscaled and the plume head rise time is reduced at the same time.
The left column of Figure 8 differs from the left column of Figure 7 in that there is a regional refinement using the Nelson and Grand (2018) model. Results are very similar to each other suggesting that global structure is more impactful to the conduit's tilt than the local structure is. We also compared results with and without regional refinement for several other tomography models and found that differences are always rather small. The second column shows the result for a different, much simpler viscosity structure. Comparison shows that the predicted conduit shape is rather robust with respect to viscosity structure, but the amount of tilt depends on overall viscosity magnitude: Somewhat lower overall viscosity leads to somewhat stronger tilt in the right column. In the third column we use a four-layer model with an additional viscosity jump at 1,000-km depth, following the suggestions of Rudolph et al. (2015) with conduit source depth at 2,800 km near the base of the mantle. Again, results remain very similar. Some other viscosity structures that we tested gave unrealistic results such as plume conduits that are too strongly tilted even in the case with an initially vertical conduit at 17 Ma. Such unrealistic results-not only for Yellowstone but for many other hot spots as well-occur for models with comparatively low viscosities and hence too fast flow, or models with a thin layer of very low viscosity between upper and lower mantle, such that flow in the upper and lower mantle is largely decoupled. An example (Roy & Peltier, 2015, VM6) is shown in the fourth column of Figure 8. In contrast, the other viscosity models for which results are displayed here and in supporting information Figure S2 yield generally reasonable tilts for mantle plumes globally (e.g., Torsvik et al., 2019). Similar results were also found for the A-family viscosity model fromČížková et al. (2012), which has similar magnitude and shape as the model in Figure

Discussion: The Rather Long Rise Time of the Yellowstone Plume Head ≈100 Myr
The best agreement between model predictions and the tomographic image of the plume conduit is found for plume head rise times of ≈80 Myr or longer. For shorter rise times, with the plume head reaching the surface 17 Ma, corresponding to the age of the CRB, the plume conduit would remain more vertical than tomographically imaged. Therefore, we favor models with longer rise times to allow for sufficient time to acquire enough tilt to match the seismic observations. For some of our results we have scaled up tomography models by a factor 1.5, hence scaled up flow speeds. For some of the tomography models, this leads to flow speeds such that a plume head rise time of around 60 Myr is sufficient to achieve a conduit tilt comparable to observations. Faster flow speeds and hence lower implied rise times can also result for lower viscosities. These cases are not definitely excluded but tend to yield less realistic results for hot spot motion and plume conduit shapes for other plumes (Steinberger & O'Connell, 1998). Only with a fairly high viscosity in the lower mantle, and correspondingly slow flow speeds, are modeled hot spot motions sufficiently slow to be compatible with observations (Richards, 1991;Steinberger & O'Connell, 1998). Hence, we favor results to viscosity and density structures, which yield reasonable results for hot spots globally and which imply plume head rise times of ≈80 Myr or longer as opposed to ≈60 Myr or even less.
Alternatively, the plume head could have reached the surface earlier than 17 Ma. We find that rather similar conduit shapes can result for an initially vertical conduit (corresponding to a plume head rising very fast) at an earlier time, and a plume head reaching the surface at 17 Ma after a longer rise time. Presumably, intermediate cases combining somewhat earlier plume head arrival with a finite head rise time can again give similar results, because already results for these limiting cases are rather similar.
For an earlier plume head arrival, earlier parts of the hot spot track may have formed on the Farallon plate and subsequently been subducted or accreted. Wells et al. (2014) describe such a hypothesis, showing that-consistent with current plate tectonic reconstructions-the Yellowstone hot spot could have provided a 56-to 49-Ma source on the Farallon plate for the Siletzia LIP, which accreted to North America by 50 Ma and is now located in coastal Oregon, Washington, and southern Vancouver Island. A similar scenario is also discussed by McCrory and Wilson (2013). In this case, the CRB may be due to transient ponding of plume tail upwelling beneath a slab before it ruptured in the Miocene, which may also explain that CRB are smaller than most LIPs. If such a secondary pulse causes the plume conduit in the upper mantle to straighten up and become nearly vertical again, its subsequent behavior will be similar to after eruption of the plume head: the conduit gets tilted over by eastward upper mantle flow (see supporting information Movie S1).
Combining the corresponding eastward hot spot motion with absolute plate motion yields a better match with the observed track than in the case where the plume conduit is already tilted in the upper mantle at 17 Ma, such that no eastward hot spot motion after 17 Ma is modeled. Since our model does not include secondary pulses that may straighten the conduit in the upper mantle, we rather focus here on the scenario with a plume head taking longer time to rise through the mantle and reaching the surface about 17 Ma.
An important restriction of our model is that the rise of the plume head is not modeled fully dynamically, and, in particular, we do not consider the effect that a vertical ambient mantle flow component may have on it. Hence, it is an a priori assumption that the plume head can reach the surface in a region with lots of sinking slabs, and not part of the modeling.
As a further simple test to see if a long plume rise time is reasonable, we can use Stokes' law for the rising speed of a sphere: If we assume a plume head diameter r h = 1,000 km (Campbell, 2007) and a density contrast of Δ = 30 kg/m 3 , which is realistic for a thermal plume (Steinberger & Calderwood, 2006), we obtain a plume head rise time of somewhat more than 150 Myr if we use viscosity at a given depth to compute rise speed v at that same depth. This is approximately compatible with our result given the uncertainties in density contrast and plume head diameter.
In contrast, Torsvik et al. (2019) find that rise times for plumes in the vicinity of LLSVPs are probably only about 30 Myr or less. For longer rise times, plume heads would need to start lifting off the core-mantle boundary further away from the LLSVPs, in order to appear at the observed surface locations, whereas for short rise times they would in general start close to the margins of the LLSVPs, which are likely Plume Generation Zones (Burke et al., 2008;Torsvik et al., 2006). There could be a number of reasons for this discrepancy with our results, which are more in line with Stokes' law. The CRB are less voluminous than other LIPs and may therefore be caused by a comparatively small plume head. According to the Stokes formula, a smaller plume head rises more slowly than a faster one. Also, mantle above and in the vicinity of LLSVPs could be generally hotter and less viscous than elsewhere, allowing plumes to rise faster. Hotter mantle would also cause large-scale upward flow in these regions, further assisting the rise of plume heads.
Despite their apparently faster rising, the modeled amount of tilt of other plumes is often similar: This is, because they often appear to be much older, corresponding to longer associated hot spot tracks. Hence, their conduits also may have acquired a large tilt by being sheared in mantle flow during that time. In contrast, French and Romanowicz (2015) state that their tomographic images are in support of nearly vertical plume conduits. Hence, it remains unresolved whether the tilt of the Yellowstone plume conduit is unusual in a global context.
Our results also show that it could be possible to match the conduits tilt as a function of depth in addition to the overall tilt direction and magnitude. Conduit tilt is a result of shearing due to different flow velocities at different depths in the mantle. In most models, flow in the lower part of the mantle beneath western North America is toward the south to southwest and becomes stronger with depth, as the core-mantle boundary is free-slip and viscosity likely decreases in the lowermost mantle. This explains the overall tilt direction observed in the lower part of the mantle. This southwestward flow is likely a consequence of former subduction northeast of Yellowstone ( Figure 2): Although since ≈50 Ma, subduction has occurred southwest of Yellowstone, as North America moved westward, and its volume flux has been reduced, there are still vast amounts of slabs sinking beneath North America northeast of Yellowstone, as they probably take about 200 Myr to sink to the lowermost mantle. These slabs are imaged by tomography models, explaining the predominant directions in our flow models. More recent subduction southwest of Yellowstone could counteract, however it appears that the combined effect of flow driven by hot upwellings above the Pacific LLSVPs and subducted slabs beneath northeastern North America is dominant: For example, Figure 4 shows downward flow, possibly related to subduction, beneath the northwestern Pacific, but nevertheless, the horizontal flow component is southward to southwestward all across the downward flow.
In contrast, upper mantle flow in the vicinity of Yellowstone is likely approximately eastward, driven by subducted slabs sinking into the lower mantle further to the east, which can explain a different tilt direction in the upper part of the mantle. Tilted conduit segments can also rise due to their own buoyancy or aided by upward mantle flow. Tilt angle varies with depth, and-similar to the tomography model-we often find two separate depth intervals, in the upper and lower half of the mantle, respectively, where the tilt is strongest. By adjusting model parameters, it might be possible to tune the model such that the depths of strongest tilt match the tomography model observations.
A strong limitation of using backward advection in our models to compute past mantle flow is that diffusion is not included, as running diffusion backward in time would be inherently unstable, and the computed density structure becomes increasingly stratified further back in time. Accordingly, we limit backward advection to 68 Ma. For this time frame, the limitations of our approach are not too severe, as results computed with backward advection remain broadly similar to those with constant present-day flow. To improve the "retrodictions" of past mantle structure, the adjoint method (Bunge et al., 2003;Ismail-Zadeh et al., 2004; can be used instead of backward advection. This requires iterations and therefore more computation time.
Given the simplicity of our model, we cannot expect to model the plume conduit in every detail. We regard the similarity achieved between most geodynamic model predictions and the tomography model-similar amount and direction of tilt, with two depth intervals of stronger tilt, already as quite remarkable. It supports both the fact that Yellowstone is a lower mantle plume, which took on the order of 100 Myr from its initiation near the core-mantle boundary to the first volcanoes and corroborates direction and magnitude of modeled lower mantle flow. We leave it to future models to further improve the fit and possibly find implications about mantle properties that allow an improved fit.

Outlook and Speculation: The Cause of the Yellowstone Plume
The inferred source location of the Yellowstone plume in the lowermost mantle is beneath southern California or northern Baja California. Besides the model of Nelson and Grand (2018), this is also seen as a region of low seismic velocity in several global tomography models (Amaru, 2007;Auer et al., 2014;Houser et al., 2008) and D ′′ models (Castle et al., 2000;Kuo et al., 2000). This could be one of perhaps two Low Shear Velocity Provinces (LSVPs)-similar to, but smaller than, the two LLSVPs (Burke et al., 2008), but it could also be part of a northward extension from the eastern part of the Pacific LLSVPs, as suggested by the SAVANI (Auer et al., 2014) tomography model.
The cause of the Yellowstone plume at the core-mantle boundary is still unclear. Torsvik et al. (2006) and Burke et al. (2008) suggest that plumes tend to rise at the margins of LLSVPs. However, the source of the Yellowstone plume that we infer here is far from LLSVPs and rather close to a region where a lot of slabs have been subducted during the Mesozoic and Cenozoic.
These slabs are imaged more clearly than plumes in the lowermost mantle due to their planar compared to tubular geometries and the absence of wave front healing effects that are caused by slow anomalies (Sleep, 2006). In particular, tomography models show an apparent hole in the slabs toward the southwest of Yellowstone, in the region where the plume conduit rises from the lower mantle in our models. This hole could enable the plume to rise from the lowermost mantle, despite being located in a region with abundant subduction, or alternatively, it could have been caused by the plume penetrating downgoing slabs, but this has been suggested to be unlikely (Leonard & Liu, 2016).
In the model of Steinberger and Torsvik (2012) mantle flow is caused by slabs sinking into the mantle at the locations of subduction zones during the past 300 Myr. This causes a thermochemical layer at the base of the Figure 9. Map views (right) and cross sections (left; along the black lines in the right panels) for the model of Steinberger and Torsvik (2012) at times 120 Ma (top), 80 Ma (second row), and 40 Ma (third row) and present-day (bottom row) to illustrate a conceptual model of how subduction zones intruded onto the Pacific LLSVP, and subducted slabs split off a part of it, from which the Yellowstone plume rose. Slabs and plumes are only shown below 500-km depth. Orange stars indicate surface position of Yellowstone, and red circles plume position at CMB as seen by tomography (Figure 1). mantle to be mainly accumulated into two piles, one beneath the Pacific and one beneath Africa, resembling the LLSVPs. Flow driven by these slabs also causes plumes to be initiated mainly along the edges of these piles. Because the model of Steinberger and Torsvik (2012) does not consider lateral viscosity variations, it yields rather wide plumes, which may more actively drive mantle flow and lead to a different plume tilt than the treatment of plumes as essentially passive. If lower viscosity due to higher temperature in plumes is considered, their width is reduced but the overall pattern of plumes rising from the margins of LLSVPs remains the same (Hassan et al., 2015;Gassmöller, 2015). Figure 9 illustrates how slabs from a subduction zone, which migrated westward at the western margin of North America, may have split off part of the Pacific LLSVP and a plume may have risen from the part that was split off. In this process, slabs play an active role; hence, such a scenario should also, or even more so, be viable if plumes are less wide and hence presumably less actively influencing mantle flow.
Compared to the models that we focus on here (Figures 6-8), such an approach has the advantage that subduction, large-scale flow, and plumes are contained in the same self-consistent model, whereas here we compute large-scale flow and plume conduits in two separate steps. However, our two-step approach allows present-day hot spot positions to be exactly matched in most cases. In contrast, if the plumes self-consistently arise in the model, then they generally will not form at the exact locations of hot spots, although often in the general vicinity (Gassmöller, 2015;Hassan et al., 2015). Also, a forward approach based on subduction history generally does not match present-day mantle structure inferred from tomography, although some features are similar. Furthermore, our two-step approach is also computationally less demanding, because large-scale flow only needs to be computed once for various plume parameters. There is no need to resolve the dimensions of the plume conduit in fine detail because large-scale density anomalies are most relevant for driving large-scale flow.
The plume location in Figure 9 is several thousand kilometers southeast of Yellowstone, and the plume reaches the surface several tens of million years before the age of the CRB in that model, with the shape of the plume conduit at 40 Ma giving a better match to the observed one than the one modeled for present day. A location further west could possibly result if the plate reconstruction was modified to include Pacific intraoceanic subduction. A later rise time of the plume could possibly result with slower slab sinking, as the sinking speeds of slabs in the model is faster than what is inferred from matching predicted slab locations with tomography (Domeier et al., 2016;van der Meer et al., 2010). Further, more advanced modeling of mantle convection and plumes, with lateral viscosity variations similar to Dannberg and Gassmöller (2018) coupled to improved plate reconstructions, could help to support or refute these speculations. If such models could approximately fit the location and timing of Yellowstone, they could again be directly compared to tomography.

Conclusions: Numerical Models Provide Support for Whole-Mantle Plume Scenario
We have computed the shape of the Yellowstone plume conduit as the plume rises through the large-scale mantle flow field to see if the tilt magnitude and direction observed by Nelson and Grand (2018) is reasonable. In order to assess the variability and robustness of results, we varied a number of modeling parameters and assumptions: These include the age of the plume, the time it takes the plume head to rise through the mantle, buoyant rising speed of the plume conduit, and mantle density anomalies and mantle viscosity structure. We find that in most cases computations yield that the plume conduit rises from the south to southwest-as is also inferred from mantle tomography. The amount of tilt is more variable and can be adjusted to match tomography by varying plume age or plume head rise time, as well as the magnitude of mantle density anomalies and overall mantle viscosity. We find that, for models compatible with other constraints, a good fit can be obtained if the plume head took rather long-around 90 Myr or longer-to rise through the mantle. This is probably longer than the typical rise time of other plumes, which are mostly in the vicinity of LLSVPs. The tilt of the conduit varies with depth, as segments of the conduits get variably distorted through shear flow, and subsequently rise buoyantly, or aided by upward flow. The conduit shape inferred from tomography can therefore provide information about mantle flow. The unusual location of the Yellowstone plume in the vicinity of subduction zones could be due to parts of the Pacific LLSVP being separated off by intraoceanic subduction, and a plume rising from the separated part, triggered by slabs approaching the lowermost mantle. Whether, and under what circumstances, such a scenario is feasible could be the theme of future numerical models.