Recent seismic imaging of the mantle beneath western North America reveals complexities interpreted as structures ranging from plumes to lithospheric drips and slab fragments. A prominent high-velocity “curtain” beneath Idaho has been interpreted as a remnant of the subducted Farallon plate left dangling within the upper mantle since >40 Ma. Consequently, using numerical models, we explore the rheological, chemical, geometrical, and dynamic conditions under which a slab fragment might persist in the mantle for tens of millions of years. With thermal buoyancy alone, stalled slabs extending to 500 km depth tend to detach and sink vertically within ∼17 m.y. for the slab age and rheologic conditions explored here, and shorter slabs <300 km deep have the greatest impact on delaying detachment up to 28 m.y. Otherwise, we find that an unrealistic chemical density contrast of 90 kg/m3 with respect to the mantle is required for the stalled slab to remain attached to the lithosphere >40 m.y. An increase in upper- to lower-mantle viscosity contrast (1.4× to 100×) can slow sinking velocities and extend slab dangling time by up to 5 m.y. Dynamic effects such as those arising from active nearby subduction only slightly delay or do not affect stalled slab detachment timing but do affect the geometry of the slabs as they respond to suction pressures within the wedge. Overall, a combination of buoyant, viscous, geometric, and dynamic factors may allow cases of extended slab stalling, and conditions we explore here within realistic ranges can so far account for a delay of up to 28 m.y.
Increasingly detailed seismic imaging of the upper mantle beneath the western United States has been facilitated by the unprecedented data coverage and resolution of EarthScope’s USArray, revealing complexities of structures that motivate questions of their origins and dynamics. The complex heterogeneity of the region has been a topic of studies addressing the region as a whole and the potential overlapping effects of or coupled interactions among a number of processes (Xue and Allen, 2010), as well as separate investigations of features or processes such as Yellowstone plume dynamics (e.g., Manea et al., 2009; Smith et al., 2009), small-scale convection or mantle “drips” (e.g., Zandt et al., 2004), and slab remnants or gaps from the subduction history (Roth et al., 2008; Liu and Stegman, 2011).
A particular interest is the number of high-velocity features, such as those identified by Schmandt and Humphreys (2010a) as “lithospheric instabilities,” prominent in recent tomographic images. Such features may be either of a subducted slab nature or a delamination origin, as in the form of a “drip” or Rayleigh-Taylor instability. The ability to distinguish the two is important for understanding their origins, mechanisms, and dynamics. For instance, beneath the southern Sierra Nevada–Great Valley region of the United States, a roughly elliptical high-seismic-velocity feature >100 km wide and ∼200 km deep has been described as a large lithospheric “drip” interpreted as foundering of the dense batholithic root of the southern Sierra Nevadas (Biasi and Humphreys, 1992; Zandt et al., 2004). The feature and potential root-foundering origin have been addressed, for instance, in visco-elasto-plastic thermo-mechanical numerical models that have allowed the present-day images to be linked to observed uplift, magmatism, and extensional tectonics (Le Pourhiet et al., 2006). Another interpretation is that this same feature, also referred to as the Isabella anomaly, could be a remnant of the Monterrey slab now attached to the Pacific plate (Forsyth et al., 2011). It has also been proposed that the anomaly could be a composite of both batholithic root foundering and translation of a slab remnant.
Another prominent feature in tomographic models of the upper mantle beneath Idaho and eastern Washington is a large, vertically extending, high-velocity “curtain,” as addressed in a conceptual model by Schmandt and Humphreys (2011). The large volume of this anomaly, extending to depths between 230 and 600 km and with peak seismic anomalies of dVp = 1.5%–2.5% and dVs = 2.5%–4.0%, and its elongated shape with depth and slab-like (though curved) geometry in three dimensions have been interpreted by Schmandt and Humphreys (2011) as a remnant of subducted lithosphere rather than a lithospheric “drip.” Based on the geologic and magmatic history of the region, particularly a westward jump in subduction-related volcanism from the Idaho batholith to the Cascadia volcanic arc, Schmandt and Humphreys (2011) proposed that the high-velocity feature is a fragment of the Farallon plate that stalled between 50 and 40 Ma, after which subduction reinitiated farther west at what is currently the Juan de Fuca subduction zone (Fig. 1D). The timing constraints require the relict slab to have remained “dangling” in the upper mantle for the past 40–50 m.y. leading up to the present configuration seen in the tomography (Fig. 1D, last panel). The prolonged stalling might necessitate that such a slab became nearly neutrally buoyant with respect to the surrounding mantle, potentially explainable by loss of eclogotized crust, as suggested by Schmandt and Humphreys (2011), or even wholesale density segregation of the lithospheric mantle from the crust (Lee and Chen, 2007). The uncertainty regarding potential mechanisms for such extended stalling motivates an investigation into the buoyant, rheologic, and dynamic conditions under which slab fragments may persist in the upper mantle rather than detach and sink.
Conditions that might delay or prevent slab sinking include high lithospheric yield strength, high upper-mantle viscosity, and a positive chemical density contrast between subducting plate and mantle. Another possibility not previously addressed is the dynamic effect of flow from external processes, such as from nearby subduction. In this study, we use fully dynamic numerical models to explore the effect of slab geometry, rheology, and density to distinguish cases in which stalled lithosphere might detach or linger in the upper mantle for tens of millions of years. In particular, we present models that help test conceptual models such as that proposed for the Schmandt and Humphreys (2011) relict slab scenario in order to converge on a better understanding of the origins of high-seismic-velocity features imaged beneath the western United States.
We present time-varying two-dimensional (2-D) numerical models that include a visco-plastic rheology using the finite-element code CitcomS (Zhong et al., 2000; Tan et al., 2006). The software solves for thermal convection problems and allows for the use of complex boundary conditions, large variations in rheologic parameters, and mesh refinement. For a detailed description of the model setup, including implementation to simulate the weak region interface between plates, the viscosity law, and rheological parameters, see Burkett and Billen (2009).
Our model domain extends from the surface to 1350 km depth and 0°–29° in longitude, with resolution ranging from 1.4 km × 1.5 km (highest resolution in the mantle wedge) to 14 km × 10 km (bottom side boundaries); the domain is a thin (essentially 2-D) cutout of a sphere. Boundary conditions are free slip on all four boundaries. Temperatures range from 0 °C to 1400 °C, and the initial temperature structure includes surface plates from a half-space cooling model, with two slabs that have subducted to some depth according to a previous history and defined by the position of a spreading ridge positioned outboard of the trench.
We proceed from the 40 Ma frame shown in the motivating conceptual cartoon (Fig. 1B) by starting from an initial thermal structure that defines (1) a subducted plate that stalled when the slab was 500 km deep (referred to as the stalled or “dangling” slab), and (2) a reinitiated subduction zone shifted 450 km west of the terminated subduction zone (see Fig. 2A). The overriding plate initially has an age of 50 Ma (such that the 1000 °C temperature contour has a depth of 66 km).
We test cases in which there is no subduction of the westward plate, as well as cases in which only the western slab (not the stalled slab) is allowed to actively subduct. As in previous models (Burkett and Billen, 2009), subduction is fully dynamic in that we do not impose plate velocities; subduction occurs due only to slab pull along a weakened shear zone between the subducted and overriding plate, which is prescribed by a region of lower viscosity. In order to isolate the stalled slab from effects of ongoing subduction for some cases presented here, we do not modify the thermal structure but merely prevent subduction of the westward plate by deactivating the shear zone weakening effect.
To test cases of varied stalled slab chemical buoyancy in some cases, we define a density structure by tracers that fill a region defined below 80 km depth and within the 800 °C contour of “dangling” slab. The tracers are assigned a density anomaly, in cases here, a lower density by 0–90 kg/m3 than the reference density of 3300 kg/m3.
We include a composite upper-mantle rheology, a Newtonian-only lower-mantle rheology, and a plastic rheology that overrides the composite rheology when stress exceeds the yield strength (depth dependent to a maximum of either 500 or 1000 MPa). The rheology is implemented as in Burkett and Billen (2010) with the following exceptions: (1) We use a stress exponent of 3.0 instead of 3.5 for dislocation creep of olivine, which has been the preferred value in global models most consistent with observed plate motions and plateness (Stadler et al., 2010), and (2) we now explicitly account for the use of the second invariant of the stress and strain rate in the numerical models compared to the use of the axial strain rate and differential stress in the experimental flow law following Gerya (2010).
We present a summary of model results in Figure 3, showing the dependence of stalled slab residence time (the “dangling” time of the stalled slab before detachment) on slab maximum yield strength, buoyancy (chemical or dependent upon slab length), and dynamics (e.g., active nearby subduction). Table 1 also lists the initial parameters and detachment results for all models presented here.
General Model Dynamics
For models in which there is no active subduction, the dynamic flow is driven only by sinking and detachment of the stalled slab when buoyancy forces (thermal and chemical) overcome slab strength and viscous resistance from the surrounding mantle. Detachment of the dangling slab occurs for all cases explored here, as characterized by stretching and “necking” of the slab, occurring at ∼150 km depth in most cases, until separation at the 1050 °C temperature contour (see Fig. 4 for examples of detachment).
No Nearby Subduction
For our reference case, we start with a 500-km-deep dangling slab, no chemical density contrast with respect to the mantle, and no active nearby subduction. We find that detachment occurs within 15.6 m.y. (Fig. 4D), which is much earlier than the >40 m.y. residence time predicted by the Schmandt and Humphreys (2011) interpretation shown in Figure 1. The negative thermal buoyancy of the dangling slab causes it to detach by yielding and necking of the slab at ∼150 km depth.
No Nearby Subduction, Increased Chemical Buoyancy
A dangling slab chemical density contrast of −90 kg/m3 with respect to the mantle can partially offset the thermal negative buoyancy and delay slab residence time to 40–50 m.y. (Fig. 2B). However, such a slab versus mantle density contrast (∼3%) is unrealistic, and detachment still eventually occurs, suggesting chemical buoyancy alone cannot easily explain extended dangling times. The overall trend of delayed dangling time with increasing chemical buoyancy can be seen in Figure 3A, despite extending the value up to the unrealistic range of 90 kg/m3.
Effects of Maximum Yield Strength
For cases in which additional chemistry is not included, a greater slab strength of 1000 MPa versus 500 MPa delays timing of detachment by ∼6 m.y. (model 11 vs. model 12; Table 1; Fig. 3A). However, when the slab has additional chemical buoyancy from a density contrast greater than 15 kg/m3, there is no difference in detachment timing between yield strengths in the 500–1000 MPa range (Fig. 3A). This effect of the additional nonthermal buoyancy overrides the influence of plastic yielding by sufficiently supporting the slab long enough for thermal weakening to control detachment rather than plastic yielding (for further discussion of thermal- vs. yielding-controlled detachment, see Andrews and Billen, 2009).
Effects of Mantle Viscosity
Greater upper- or lower-mantle viscosities provide more viscous support at depth to a sinking slab and thus cause sinking velocities to decrease. We test models with a range of lower-mantle viscosities, from 1.4 to 100 times that of the upper-mantle viscosity, as might be a possible range based on geoid constraints for subducted slab density models (Hager, 1984). We find that the time it takes for the slab to detach increases with increasing lower-mantle viscosity. For the full two-order-of-magnitude viscosity jump, the delay of detachment time extends almost 5 m.y. (see models 15–18; Table 1; Fig. 3A). Additionally, an increase in the overall minimum background viscosity from 1 × 1019 Pa-s to 1 × 1020 Pa-s delays the detachment by only 1 m.y. (models 6 and 7 vs. model 1).
Effects of Active Nearby Subduction
Subduction is allowed to occur by inclusion of a weak shear zone of weakened viscosity that allows the western slab to respond to pull from its negative buoyancy and subduct dynamically (Figure 4A with ongoing subduction vs. Figure 4D without subduction). We find that the presence of active nearby subduction results in only a minor delay of <0.5 m.y. The more notable effect is a slight westward deflection of the stalled slab toward the active subduction zone (Fig. 4A vs. Fig. 4D) in response to suction pressures in the wedge (Fig. 5D). When the rate of nearby subduction is faster, ∼5–13 cm/yr versus ∼3–5 cm/yr, as allowed by a weaker shear zone viscosity, there is no change in the time the slab takes to detach, but there is a notable shallowing of the subducting plate dip (Fig. 4A vs. Fig. 4B). For an offset of 350 (vs. 450) km between the subducting and stalled slabs, detachment is delayed ∼1 m.y., and shallow part of subducting plate dip shallows, but the deeper slab remains vertical (less room to respond horizontally due to the proximity of the stalled slab; Fig. 4C). Deflections and shallowing of the subducting slabs occur in response to suction pressures in the wedge (Figs. 5C and 5D) and are similar to effects also seen in modeling of subduction wedge dynamics by Manea and Gurnis (2007).
Short Slabs and Combined Buoyancy Effect
For a stalled slab initially at a depth of 300 km versus 500 km (less thermal negative buoyancy since it is shorter), the detachment time is delayed by ∼10 m.y. (model 8 vs. model 1 with nearby subduction, or model 9 vs. model 2 without subduction). The shorter and less buoyant slab not only resides longer before detachment, but also is more easily deflected toward the active subduction zone (Fig. 4E vs. Fig. 4A) in response to the wedge suction pressure (Fig. 5B). A case composed of a combination of multiple conditions (within reasonable ranges) that delay detachment (e.g., 100× jump in lower-mantle viscosity, 15 kg/m3 chemical density contrast, no nearby subduction, short 300 km initial stalled slab) results in a detachment time of 28.3 m.y., though the major factor in delaying time to detachment is the minimization of slab negative buoyancy by shorter slab remnants.
Slab detachment is an expected result of stalled subduction scenarios when the dynamics depend only on thermally determined density anomalies since subducted plates are consistently colder and more negatively buoyant than the surrounding mantle. Such detachment as controlled by the thermal negative buoyancy (slab pull) opposing slab strength and viscous resistance is demonstrated here and has been explored in greater detail in previous models by Andrews and Billen (2009). Subducting plates, however, exist in a mantle that is likely affected by flow fields not associated with the immediate surroundings, and subducted lithosphere may include chemical buoyancy variations. We build upon previous modeling studies and specifically test the Schmandt and Humphreys (2010) dangling slab scenario by addressing further questions concerning the effects of slab chemical variations and nearby subduction dynamics on the residence time of stalled slabs within the upper mantle. Although we do find that effects such as chemical buoyancy, slab and mantle rheology, and flow patterns effect the timing, the stalled slab length turns out to have the largest influence on reducing the negative buoyancy and extending slab residence times beyond 20 m.y.
In addressing the effect of variations in stalled slab chemical buoyancy, the necessary addition of an unrealistically extreme difference of 90 kg/m3 for lower slab versus mantle density (Fig. 2) to delay residence time to the 40–50 m.y. range suggests that it may not be reasonable to call upon chemical buoyancy alone to explain lingering slabs within the upper mantle. For instance, if the slab consists of lithosphere chemically depleted in a basaltic component, its density may be up to 1%–1.5% less dense than the surrounding mantle (Lee, 2003), which would at most correspond to a 30–50 kg/m3 lower density. Furthermore, cooling of surrounding mantle as the slab warms (contours widen beyond the chemically buoyant region of the slab) does allow the thermal negative buoyancy to eventually pull the slab from the surface (Fig. 2B), thereby illustrating the temporary, even if delayed, nature of the slab stalling.
Many observations for regions involving recent cessation of subduction indicate that their associated slabs have detached or may be in the process of detaching. For instance, detachment of young lithosphere, such as may have occurred along Baja California and beneath western and central Mexico upon ridge-trench collision, seems to be a process that occurs on the order of a few million years, based on dynamic numerical modeling (Burkett and Billen, 2009; Burkett and Billen, 2010) and observations of magmatism linked to the timing of the 2-D break-off process and 3-D tear propagation (Ferrari, 2004; Michaud et al., 2006; Pallares et al., 2007). Within the India-Asia collision zone, a currently detaching remnant of Indian lithosphere may have been imaged, as evident from the clustering of seismicity and fault-plane orientations beneath the Hindu-Kush, and would have initiated within the last few million years some time after continental collision and closure of the Tethyan Ocean (Lister et al., 2008). Throughout the Mediterranean-Carpathian region, other instances of detachment due to continental collision and ocean basin closure have been described as occurring with rapid tear propagation rates on the order of a few million years (Wortel and Spakman, 2000). Overall, the observations linked to slab detachment processes seem to indicate a rapid process occurring within a few million years, and cases of stalling for longer than tens of millions of years may require particular buoyancy or dynamic conditions.
Most reported instances of “stalled” slabs tend to refer to the brief transitory time of stalling during collision (e.g., cases mentioned previously resulting in detachment after a few million years) or cases in which the stalling occurs for a short slab remnant, particularly following a previous occurrence of slab detachment. For instance, along California and Baja California, seismic tomography putatively reveals evidence for fossil slabs attached to the remnant microplates imaged to at least 150 km depth that were left abandoned after ridge-trench collision (Forsyth et al., 2011). In this particular setting, which would involve translation along the Pacific–North American plate boundary, such short and young slab fragments might not only be sufficiently buoyant but also feasibly able to withstand viscous drag effects during translation (Pikser et al., 2012).
Previous numerical models similar to those presented here suggest that after stalling of subduction, slab detachment times may range between ∼5 and 25 m.y, with the longer residence times for strong (high maximum yield strength) and old or short slabs (Andrews and Billen, 2009). These cases include a nonlinear upper-mantle rheology, which is important for allowing lower viscosity where strain rates are high, although they assume thermal-only buoyancy to drive flow velocities. In other numerical modeling studies of stalled subduction, which also include composition and melting in addition to thermal density effects, the slab detachment times range somewhat longer to 36 m.y. (Gerya et al., 2004; Baumann et al., 2010). Overall, it is difficult to explain the lingering presence of slabs in the upper mantle for longer than 36 m.y. for the variations in rheology and thermal structures within realistic ranges so far explored in numerical modeling studies.
We realize the limitations of these models given the 2-D geometry, absence of a “mantle wind” or reference frame to better compare to western North America conditions, and inclusion of phase transitions such as mantle viscosity jumps without buoyancy alterations from changes in slab density. However, the models provide useful constraints on the first-order processes that may help explain dynamics in more complicated scenarios. For instance, despite model simplicities, even the dynamic result of stalled slab deflection toward the subduction zone (Fig. 4) is consistent with the apparent westward deflection of an observed high-velocity “curtain” (Fig. 1D).
We have chosen to focus on testing the stalled slab scenario, realizing the need to eventually distinguish such dynamics from the viable alternative of a lithospheric drip scenario. The current understanding of both processes is in need of improvement, as evidenced by the mixed interpretations of evidence for seismic anomalies interpreted in the context of either process. For instance, crustal constraints may provide evidence for active continental lithospheric delamination beneath Romania where previous studies have interpreted observations in the context of a slab detachment scenario (Knapp et al., 2005). On the other hand and as previously mentioned, the reported mantle “drip” beneath the Sierra Nevada Mountains (Biasi and Humphreys, 1992; Zandt et al., 2004), or Isabella anomaly, has more recently been interpreted as a possible Monterrey slab fragment (Forsyth et al., 2011).
Our models demonstrate that stalled slabs tend to detach and sink well before 50 m.y. when slabs are entirely thermal. For the case of dangling slabs in which buoyancy is solely thermally controlled, the detachment time is on the order of 15 m.y. For the conditions explored in these models, an additional positive buoyancy force from a slab 90 kg/m3 less dense than the surrounding mantle is necessary for the slab to stall and “dangle” within the upper mantle for up to 40–50 m.y. Given that such densities are likely unrealistic beyond 1.5% or a 50 kg/m3 contrast, which would approach depletion similar to cratonic mantle (Lee, 2003), model results suggest that chemical buoyancy alone is not sufficient to explain stagnant slabs lingering longer than ∼20 m.y. An increase in upper- to lower-mantle viscosity contrast, however, extends slab dangling time by up to 5 m.y. (across a range from 1.4× to 100× contrast) due to a slowing of sinking velocity through viscous support by a more viscous lower mantle. Effects of active subduction near a stalled slab tend not to greatly affect timing of detachment, although geometric effects from flow-induced pressures include trenchward deflection of the slabs or shallowing of subducting slab dip. Overall, we conclude that a combination of buoyant, viscous, geometric, and dynamic factors may be responsible, and conditions we explore here within realistic ranges can so far account for a delay up to 28 m.y., with the greatest delay effect being the stalling of short slab remnants. Any further delay may require more extreme buoyancy, dynamics, three-dimensional geometric effects, or consideration of other processes not explored here, such as lithospheric delamination.
This work was supported by the National Science Foundation (CMMI-1028978), and the Caltech Tectonics Observatory (by the Gordon and Betty Moore Foundation). The original CitcomS software was obtained from CIG, Computational Infrastructure for Geodynamics (http://geodynamics.org), and code modifications for the shear zone and composite viscosity were made by Magali I. Billen. This is contribution 205 of the Tectonics Observatory, Caltech. We thank T.V. Gerya and V.C. Manea for constructive comments and M. Billen for discussions.