Plume‐Lithosphere Interaction and Delamination at Yellowstone and Its Implications for the Boundary of Craton Stability

Delamination of the lower lithosphere has occurred under different tectonic settings, thereby inducing strong lateral lithospheric strength heterogeneities. Here, we examine a recent case of plume‐induced lithosphere delamination associated with the Yellowstone plume. We propose that interaction between the Yellowstone plume and western North American lithosphere is what led to delamination of the lower lithosphere beneath the Columbia River Basalt (CRB) roughly ∼17 Myr ago. The excess melting that occurred when this “hole” was being filled by Yellowstone plume material was the primary trigger for the CRB excess melting event. The delaminated lower lithosphere currently lies to the northeast of the CRB; it can be clearly identified in recent seismic tomographic models. We suggest that both weak zones, for example, lithosphere‐crossing faults or suture zones, or vertical weaknesses associated with hotspot tracks, and strong zones (regions not underlain by weak mid‐lithosphere discontinuity material) can become the lateral boundaries of lithosphere delamination.


10.1029/2021GL096864
2 of 11 beneath the eastern part of the North China Craton (Zhu et al., 2012(Zhu et al., , 2021;;and references therein) and the western part of Wyoming Craton (Foster, et al., 2006;Hansen et al., 2015), with both of these events apparently related to subduction.Such nonuniform modification of cratons will shape lateral variations in cratonic lithosphere thickness (Hoggard et al., 2020).
We first discuss a possible location for where delaminated lower CRB lithosphere should lie in the mantle, based on the geometry of the Yellowstone hotspot track and the present-day GPS-based deformation field for western North America (Kreemer et al., 2014; https://www.unavco.org/software/visualization/GPS-Velocity-Viewer/GPS-Velocity-Viewer.html).Basaltic eruptions that mark the Snake River Plain/Yellowstone (SRP/Y) track are clearly visible starting at ∼12-13 Ma.These are located just southeast of the Steens Mountain Basalts that erupted at ∼16 Ma (Camp & Ross, 2004;Christiansen et al., 2002;Pierce & Morgan, 2009;Pierce et al., 1992), contemporaneous with CRB volcanism.
The guiding assumptions in our search are that the delaminated lithosphere sank quasi-vertically in the mantle, and that local deformation and surface motions have not changed over the past 16 Ma.We will see that these simple assumptions fit quite well with observations.To reconstruct how northwest North America has moved and internally deformed in the ∼16 Ma since the CRB delamination event, we first back-track the location of the Steens Mountain basalt based on GPS measurements of ongoing deformation, that is, using the present estimate (stations P739 and P391) ∼2.8 mm/year relative motion with respect to stable North American Plate to "undeform" the SRP/Y track (see the green stars and the green dashed lines in Figure 1).Then, we correct the CRB location (white stars in Figure 1; stations P020, P447, P451, and WALA, ∼1.7 mm/year) to account for the effects of its internal plate deformation in the last 16 Myr, again assuming that this occurred at modern measured rates (Note this deformation pattern also typifies the general subduction-linked compression to the north of this 10.1029/2021GL096864 3 of 11 region [Kreemer et al., 2014]).Finally, we add the corrected SRP/Y track (green dashed lines in Figure 1) to the corrected CRB location (filled white circle in Figure 1).This yields our best guesstimate for the surface location above potential delaminated CRB lithosphere (Note that not including these estimated corrections for the internal deformation of the North American Plate will not strongly change the estimated backtrack location for the CRB keel.)This location for the delaminated CRB lithosphere is marked by the empty white circle in Figure 1b.(Ernst, 1988).Also shown are the Steens Mountain, Blue Mountain Province, Monument and Chief Josesph.(b)-(d) S-wave velocity perturbations at 100 km, 200 km and 600 km depths (Each % perturbation scale is shown beneath its corresponding cross-section).The white star and white filled circle are the CRB location before and after correction.The white empty circle is the location of the delaminated lithosphere.The region of thinned lithosphere beneath the CRB and the location of the possible delaminated CRB lithospheric keel are marked with dashed "oval outlines".
Interestingly, at ∼400-600 km depths, a strong high velocity anomaly exists below this region with a shape similar to the low velocity anomaly at depths of 60-200 km beneath the current surface location of the CRB.This reconstruction supports the hypothesis that the delaminated lithosphere linked to the formation of the CRB currently lies in the mantle transition zone to the east-northeast of the surface location of the CRB; and that this lithospheric keel sank essentially vertically at an average rate of ∼2 cm/year over the last 17 Myr (e.g., ∼350 km over 17 Myr).
The proposed region of the possible delaminated CRB lithospheric keel lies completely to the north of the contemporaneous Yellowstone hotspot trail.This may seem curious compared with the idea that the CRB formation is directly linked to the Yellowstone hotspot trail.Clearly, a delaminated "hole" in the lithosphere, when being filled by plume materials that are hotter than the ambient underlying mantle, has the potential for the plume material to more extensively pressure-release melt, thereby creating more extensive melting consistent with measured CRB-like basalt compositions (McKenzie & O'Nions, 1991).However, flood basalt formation only rarely occurs-CFBs are geologically rare events (Courtillot et al., 1999).At first glance, delamination does not provide an obvious argument for why the hotspot track would lie along the southern margin of the flood basalt.
Here, we are inspired by the seismic evidence that only localized thinning of lithosphere is present above the younger portions of the Yellowstone track (see Figure 1).Our hypothesis is that lithospheric faults, shear zones, and hot spot tracks can create long-lived linear belts of weakness in continental lithosphere.Such linear weak zones can act as the horizontal boundaries to where lithospheric delamination would take place.When plume material flows beneath a region, there is the enhanced potential for the delamination of one (or more) lithospheric fragments bounded by vertical weak zones in the lithosphere, so long as a semi-continuous MLD is present above the delaminating region (Shi et al., 2021).

Initial Model Design
To explore this mechanical hypothesis, we have performed a suite of large-scale two-dimensional (2D) numerical experiments (4,000 × 670 km).The goals of these experiments were both to further study the dynamics of recent plume-lithosphere interaction in the Yellowstone region (albeit in an idealized 2D geometry), and also to make an initial exploration of the factors that can potentially shape the lateral extent of plumed-enhanced lithosphere delamination (Figure 2a).The spatial resolution of the model is 2 × 2 km in the central area, gradually increasing to 30 × 30 km at margins far away from the region of interest (e.g., Gerya & Yuen, 2003a, 2003b;Shi et al., 2021).
The initial model configuration includes an assumed 200 km thick lithosphere with a hot mantle plume beneath it.Multiple vertical weak zones with a spacing of 400 km, similar to the characteristic distance between Paleozoic orogens and basins in the U.S.A. (Ham & Wilson, 1967), are prescribed in the initial continental lithosphere.
Here, a weak zone can represent the effects of either a fault zone, a fossil suture zone, or the persistent mechanical weakness created by a prior plume track.For this continental lithosphere, the initial material field is assumed to be a 20 km thick felsic upper crust, a 15 km thick mafic lower crust, and a mantle layer of anhydrous peridotite (Detailed flow laws and properties of the compositionally variable lithologies are shown in Tables S1 and S2 in the Supporting Information S1; Gerya & Yuen, 2003a, 2003b;Li, 2020;Li et al., 2016;Ranalli, 1995).For the 2D hot mantle plume material in the reference model, the initial plume "radius" is taken to be 25 km and the plume material has a temperature 100-300 K higher than its surrounding mantle (Camp et al., 2015;Colón et al., 2018;Leeman et al., 2009;Smith et al., 2009;Steinberger et al., 2019).The relatively weak MLD layer is set at a depth of 100 km with a thickness of 20 km in the reference model (Griffin et al., 1999(Griffin et al., , 2004;;Selway et al., 2015;Shi et al., 2020).It is assumed to have a constant viscosity of 3 × 10 19 Pa s.We recently used the same parameters when exploring both subduction-and plume-induced delamination (Shi et al., 2020(Shi et al., , 2021)).Once a mechanical flow connection would occur through a pre-existing vertical weak zone between uprising weak material from the mantle and the weak MLD layer, then a favorable condition for delamination would exist that would lead to lower continental lithosphere delamination as initially proposed by Bird (1979).Detailed numerical methodologies (i.e., governing equations, rheological flow laws, boundary conditions, and the implementation of partial melting) are given in the Supporting Information S1.

Models With Multiple Weak Zones
The reference model contains three vertical weak zones, each 12 km wide (Figure 2a).When hot and buoyant plume material rises beneath the continental lithosphere, part of the plume material moves upwards through the pre-existing weak zone, thereby widening the weak zone.The rest spreads horizontally beneath the lithosphere.At 10.5 Myr, the plume is directly connected to the MLD layer, which results in a favorable condition for lithosphere-delamination to occur (Bird, 1979;Lei et al., 2019;Shi et al., 2020).Because of the relatively small "radius"/volume of plume material (i.e., 25 km), this material only moves along the left side of the MLD layer (Figure 2b).The lithospheric root below the MLD layer between the left weak zone and the middle weak zone rapidly delaminates within 5 Myr, to sink into the mantle transition zone by 20.0 Myr.
High strain rates (>10 −13 s −1 ; Figure 2d) and low stresses will surround the delaminated lithosphere (<10 6 Pa; Figure 2e).However, the stress stays much higher (>10 8 Pa) within the lithosphere to the right side of the middle weak zone and at the left side of the left weak zone during the entire experiment.The MLD layer thins in this region and is not filled by plume material.This leads to "locking" of the MLD layer; locking that stops the delamination process.If we introduce two more weak layers 400 km away from the existing weak layers (i.e., five weak zones in total), this results in an evolution identical to that of the reference model (Figure S2 in the Supporting Information S1).This suggests that a vertical weak zone can indeed act as an internal boundary that controls the lateral extent of the lithosphere delamination process.
While seismic tomographic models reveal the deep plume origin of the Yellowstone hotspot, the radius and the volume flux of this (and other!) plumes remain uncertain (French & Romanowicz, 2015;Nelson & Grand, 2018).The actual volume of the hot upwelling mantle and partial melting materials could significantly affect the spatial extent of lithosphere thinning (Gorczyk et al., 2018;Hu et al., 2018;Yang & Leng, 2014).Therefore, we also explore models with a larger plume radius/volume of 50 km instead of 25 km with a total initial area ("volume") of plume material of 28,000 km 2 .In this case, both sides of the lower lithosphere between weak zones will delaminate (Figure S3 in the Supporting Information S1) instead of the single side delamination seen in the reference model.The stresses surrounding the delaminated portion of the lithosphere also decrease lower than 1 MPa during the delamination events (Figure S3 in the Supporting Information S1).
The location of the upwelling plume will often deviate from the location of pre-existing weak zones (Note that the plume will also be creating a new vertical weak zone as the lithosphere migrates over it).Therefore, we next explored a series of models with different distances between the plume (50 km radius) and pre-existing weak zones, exploring three different scenarios for plume-assisted lithosphere delamination (Figure S4 in the Supporting Information S1).When the plume-weak zone offset is relatively small (e.g., ≤60 km), then both sides of the lithosphere can delaminate.However, once the offset reaches 100 km, only the lithospheric root that is closer to the plume will be removed.An even larger lateral offset (e.g., ≥120 km) only results in a widening of the nearest weak zone, but no delamination.

Models With a Discontinuous MLD Layer
The distribution of the MLD layer itself is likely be discontinuous and complicated (Abt et al., 2010;Chen et al., 2014;Hansen et al., 2015;Hopper et al., 2014;Liu et al., 2018aLiu et al., , 2018b;;Selway et al., 2015).We also tested whether the boundary of a regional MLD layer could control the lateral extent of lithosphere delamination.We first explored a model with an 800 km wide MLD layer and a single weak zone located above the plume material (Figure 3a).Similar to the previous models, the lithosphere will delaminate along the weak layer.While the lithosphere keel below the MLD layer sinks continuously, once plume material reaches the edge of laterally discontinuous MLD material, then the lithosphere keel breaks off, and delamination stops (Figure 3b).High stresses Figure 2. Initial configuration of the numerical model and snapshots of composition, viscosity, the second-invariant of strain rate, and stress fields for the reference model.Model size is 4,000 × 670 km.White lines show isotherms in °C.In panels (a) and (b), colors indicate the possible rock types that can arise during a numerical experiment, specified as: 1: air/sticky air (Gerya & Yuen, 2003a, 2003b); 2: water; 3: upper continental crust; 4: lower continental crust; 5: lithosphere mantle; 6: sub-lithospheric mantle; 7: hydrated mantle (i.e., the initial weak zones); 8 and 9: sediment; 10: partially molten sediment; 11: partially molten upper continental crust; 12: partial molten lower continental crust; 13: mid-lithosphere discontinuity layer; 14: mantle plume.Note that molten materials (10, 11, and 12) do not form in the experiments.
at the lateral edge of the MLD layer arrest further delamination (Figure 3e).With a narrower MLD layer that is 400 km wide, while a 12 km thick weak zone will no longer lead to lithosphere delamination, a 20 km thick weak zone leads to similar evolution as in the model with an 800 km wide MLD layer (Figures 3a'-3e').These results suggest that a laterally discontinuous MLD layer could also act as a barrier to limit lithosphere delamination.

Discussion and Conclusion
While most continental lithosphere appears to be stable due to its high viscosity and relatively low density (Carlson et al., 2005;Wang et al., 2014;Jordan, 1988;Sleep, 2003;Lee et al., 2011), some regions have experienced lithospheric thinning and removal under different tectonic settings (Bird, 1979;Griffin et al., 1998;Lee et al., 2011;Menzies et al., 2007;Sleep, 2003;Wu et al., 2014).Many mechanisms have been proposed to explain the origin and dynamics of cratonic thinning and destruction (Hu et al., 2018;Liu et al., 2018aLiu et al., , 2018b;;Liu & Li, 2018;Shi et al., 2020Shi et al., , 2021;;Windley et al., 2010;Wu et al., 2014;Yang & Leng, 2014;Zhu et al., 2012;Wang & Kusky, 2019).Recent numerical models indicate that the interaction between subduction zones and/or plumes and a weak MLD layer can result in the delamination of the lower lithospheric keels (Hu et al., 2018;Liu et al., 2018aLiu et al., , 2018b;;Shi et al., 2020Shi et al., , 2021;;Wang & Kusky, 2019; and references therein).Another typical characteristic of the cratons that have experienced thinning and destruction is that only part of those cratons has been modified while the rest remains stable.For example, lithosphere thinning and removal only appear at the eastern part of the North China Craton with the Trans-North China Orogenic Belt and the North-South Gravity Anomaly Belt as the boundary (Zhu et al., 2012(Zhu et al., , 2021;;and references therein), the western part of Wyoming Craton with the Bighorn Mountains, Casper Arch, and Laramie Range as the boundary (Foster, et al., 2006;Hansen et al., 2015), and the eastern part of the Dharwar Craton with the Chitradurga schist belt and Closepet granite belt in as the boundary (Maurya et al., 2016).
Here, we propose that the lithosphere beneath the Columbia River Flood Basalts experienced delamination, while the lithosphere to its east did not delaminate.The Yellowstone plume, which led to the eruption of the Columbia River Basalts starting at ∼17 Ma (e.g., Brueseke et al., 2007;Camp & Hanan, 2008;Camp & Ross, 2004;Reidel et al., 2013), is also linked to lithospheric thinning beneath the CRB.We propose that the interaction between rising and melting plume material and the weak MLD layer in the lithosphere, limited by vertical weak zones, is what led to localized lithosphere delamination beneath the CRB.To the east of the CRB, the boundary of which is marked by the 87 Sr/ 86 Sr = 0.706 line (Ernst, 1988), the oldest Archean lithosphere has retained its thickness (Figure 1b).Because of the overall west-southwest absolute motion of the North American Plate, the delaminated lithospheric root that originated from beneath the CRB currently lies beneath the northern Rocky Mountain Range, where it can be clearly identified in recent seismic tomographic models (Figure 1; James et al., 2011).
Our numerical models imply that cratonic lithosphere heterogeneities-either in the form of pre-existing vertical weak zones and/or a laterally discontinuous weak MLD layer are what controls the lateral extent of lithospheric delamination.A vertical weak zone, in the form of an active fault zone, a fossil suture zone (e.g., the Trans-North China Orogenic Belt in North China Craton Zhu et al., 2012Zhu et al., , 2021; Chitradurga schist belt and Closepet granite belt in the Dharwar Craton Raval & Veeraswamy, 2003), or a mechanical weakness created by a prior plume track (the SRP/Y hot spot track in this study or other hotspot tracks Bird, 1979;Murphy et al., 1998;Sleep, 1997), are ubiquitous features in continents.The strength of the lithosphere is significantly affected by a horizontal weak MLD layer; when yield stresses are reduced below ∼1-10 MPa, then delamination can easily happen (Figure 2) if weak material can reach and fill a delaminating MLD.The uneven distribution of the MLD layer is also widely seen beneath cratons (Selway et al., 2015), such as cratons in the US (Abt et al., 2010;Hansen et al., 2015;Hopper et al., 2014), the North China Craton (Chen et al., 2014), and the Dharwar Craton (Kumar et al., 2013).Lithospheric delamination can also be shaped by strong lateral "edges" to the weak MLD layer, where local yield strengths can be much higher (≥100 MPa; Figure 3), and so halt keel delamination.

Figure 1 .
Figure 1.Topography and S-wave tomographic models of the Yellowstone hotspot province and its surrounding region (Modified from James et al., 2011).(a) Topography and geological setting of the Yellowstone hotspot province.The red shaded area indicates the Columbia River (CRB) and Steens Mountain flood basalts.Green circles indicate the Snake River Plain-Yellowstone hotspot track with the main eruption times shown.The dark pink dashed line is the 87 Sr/ 86 Sr = 0.706 line(Ernst, 1988).Also shown are the Steens Mountain, Blue Mountain Province, Monument and Chief Josesph.(b)-(d) S-wave velocity perturbations at 100 km, 200 km and 600 km depths (Each % perturbation scale is shown beneath its corresponding cross-section).The white star and white filled circle are the CRB location before and after correction.The white empty circle is the location of the delaminated lithosphere.The region of thinned lithosphere beneath the CRB and the location of the possible delaminated CRB lithospheric keel are marked with dashed "oval outlines".

Figure 3 .
Figure 3. Initial configuration of numerical experiments with different initial mid-lithosphere discontinuity (MLD) layer widths.Shown are snapshots of the compositional, viscosity, second-invariant of strain rate, and stress fields in the experiments.Left: 800 km wide MLD layer; right: 400 km wide MLD layer.