Catastrophic shear-removal of subcontinental lithospheric mantle beneath the Colorado Plateau by the subducted Farallon slab

The causes of Cenozoic uplift of the Colorado Plateau, southwestern USA, are strongly debated, though most hypotheses acknowledge the importance of northwest-directed subduction of the Farallon oceanic plate beneath North America since c. 100 Ma. Existing thermomechanical models suggest that the Farallon slab underthrust the proto-plateau region at ~200 km depth, removing the basal portions of its subcontinental lithospheric mantle (SCLM) root, although such small-volume subduction erosion cannot fully account for the degree of uplift observed today. Here we show via petrological modeling of lawsonite-bearing eclogite xenoliths exposed in diatremes in the center of the plateau that the Farallon slab surface penetrated through the proto-plateau SCLM at much shallower depths (~120 km) than these previous estimates, allowing shear-removal of ~80 km of SCLM – a volume up to three-times greater than previously suggested. This removal led to asthenospheric upwelling and isostatic rebound of the plateau region during the late Cretaceous to the Eocene. We posit that similar shear-removal of SCLM likely played a major role in inhibiting cratonic growth and stabilization in the Neoarchean and Paleoproterozoic – when low-angle subduction of oceanic lithosphere was more prevalent than today – accounting for the atypically thin roots existing below many ancient cratons worldwide.

A representative P-T phase equilibrium diagram for xenolith 17MSR09 is shown in Fig. 3a. The observed peak assemblage Grt-Omp-Lws-Rt-Coe-Ph is stable at P > 35 kbar and T < 640 °C. This field is limited at lower pressure by the stabilization of glaucophane -which is not observed within the stable assemblage -and at higher   www.nature.com/scientificreports www.nature.com/scientificreports/ temperatures by the consumption of lawsonite. Isolines of equal volume proportion of garnet, omphacite, and lawsonite constrain peak metamorphism at ~37 kbar and ~620 °C, with associated 1σ uncertainties of ~1 kbar and ~50 °C 38,39 . This pressure for the Farallon slab surface thus infers a minimum depth of equilibration of ~120 km, assuming lithostatic conditions, as outlined below. Petrological modeling results for xenoliths 17GR11 and 17MSR11 are shown in Supplementary Figs 6 and 7 and record equivalent P-T conditions for peak metamorphism within error. An independent estimate of peak metamorphic pressure for sample 17MSR09 was obtained via the garnet-omphacite-phengite barometer calibration for white mica-bearing eclogites 41 . Compositions for garnet rim domains and adjacent matrix omphacite and phengite in textural equilibrium (Supplementary  Table S1) produced a pressure of ~38 kbar at 620 °C, with an uncertainty of ±2 kbar (Fig. 3a). This equilibrium has a positive slope in P-T space with dT/dP = 40 °C/kbar and passes through the center of the interpreted peak assemblage field, as constrained by the observed phase assemblage and mineral proportions.
We limit the calculated peak conditions of metamorphism to lie outside of the 'forbidden zone' , defined by (Liou, Hacker & Zhang 2000) 42 as geothermal gradients colder than 5 °C/km. This geotherm represents the minimum rate of conductive heating that rocks may experience during descent into the Earth, and no known entirely mafic crust is documented to have experienced such P-T conditions. While some studies report metamorphic conditions that lie within this high-pressure/low-temperature domain (e.g. Zhang et al. 43 ), these occurrences are exclusively UHP terranes formed via deep subduction of continental crust, which provides positive buoyancy for rapid exhumation and preservation. Subducted oceanic crust or eroded mafic fragments from the base of an overlying arc subjected to thermal diffusion over periods longer than ~5-10 Myr are expected to thermally equilibrate, as the characteristic diffusion distance for this time scale at mantle conditions is ~10-15 km 41 . Thus, it is unlikely that the uppermost or middle portions of the Farallon slab crust could preserve such low thermal gradients, as geochronology of eclogite xenoliths indicates that the experienced prograde metamorphism for at least 45 Myr 33 .

Dip-angle of the Farallon plate and thickness of the sCLM. Existing thermomechanical models of
Farallon slab evolution over time suggest that flattening was promoted by the subduction of an oceanic plateau at c. 90 Ma 44 . The c. 80-30 Ma age range recorded by eclogite xenoliths [31][32][33] from the NVF suggests that metamorphic recrystallization was prevalent during the low-angle and flattened subduction stages of the Farallon plate. Thus, we interpret the peak P-T conditions determined from xenoliths 17MSR09, 17GR11, and 17MSR11 in this work ( Fig. 3 and Supplementary Figs 6 and 7) to represent the maximum depths reached by the top of the plate during these stages, assuming that it remained at approximately the same depth once it flattened and migrated inboard of the continent.   www.nature.com/scientificreports www.nature.com/scientificreports/ The results of thermobarometry shown here constrain the depth and dip-angle of the Farallon plate during this period, and so the thickness of the SCLM below the Colorado Plateau (see Supplementary Methods). In a two-layer model that considers representative densities for the upper and lower continental crust 45 , the continental Moho beneath the proto-plateau would have been at a pressure of ~11 kbar (Fig. 4a). The pressure difference of ~26 kbar between this value and that calculated here for the Farallon plate slab surface (i.e. ~37 kbar from eclogite 17MSR09; Figs 3a and 4a), combined with a representative mantle density of 3.34 g/cm 3 from mantle xenoliths recovered from the SCLM beneath the Colorado Plateau 46 , implies an intermediate SCLM root ~80 km in thickness (Fig. 4a). The Farallon slab-top was therefore located at ~120 km depth, in agreement with the present lithospheric thickness of the Colorado Plateau 7,25 (Fig. 4a). As the Colorado Plateau and adjacent Great Plains (Fig. 1a) share similar Mesozoic geological histories and overlie the same Proterozoic terranes 47 , geophysical measurements of the thickness of SCLM beneath the latter today (~200 km) 25 imply a similar thickness beneath the former prior to shear removal. Consequently, over ~80 km of SCLM must have been removed from the base of the proto-plateau's root due to northeastward progression of the subducted Farallon plate, which is greater than estimates provided by recent studies (~20-50 km) 20 .
If this depth is considered the maximum reached by the Farallon plate slab-top prior to and/or during flattening, the mean dip angle of the slab during the Late Cretaceous to the Eocene can also be estimated from our P-T data. Current geodynamic models suggest that the Farallon slab flattened ~1000 km inboard of the oceanic trench during the early stages of the Laramide orogeny [19][20][21] (Fig. 4a). If correct, during low-angle subduction and prior to flattening, the slab would have required a mean dip angle of ~7°. Sensitivity analysis considering geologically realistic uncertainties for these geometric data indicates that the mean dip angle may change by ±0.6° for every 10 km uncertainty in proto-plateau continental crust thickness and ±0.3° for every 50 km variation in absolute distance from the trench (Supplementary Fig. 8). Calculated SCLM thickness may vary by ±3 km for every 1 kbar of pressure variation and ±1 km for every 50 kg/m 3 change in density ( Supplementary Fig. 9), although the absolute magnitude of error associated with pressure determination via petrological phase equilibrium modeling is considered to be less than ~1 kbar for well-equilibrated parageneses 38 .  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Constraining the rates and styles of lithospheric-scale geodynamic processes has critical importance for validating our understanding of how the Earth has evolved throughout geological time. Furthermore, such P-T-t data obtained via petrological modeling are routinely used as benchmarks and/or constraints for numerical models of lithospheric evolution. If our knowledge of the rates and styles of metamorphism and tectonic deformation are poorly constrained, the results and implications drawn from thermo-mechanical simulations become less reliable.
While eclogite xenoliths from the NVF have been studied by many workers [31][32][33] , all previous P-T estimates for peak conditions have been determined via conventional thermobarometry, which is subject to large inter-calibration uncertainties and relies on measured mineral compositions accurately representing those attained at the metamorphic peak 38,48 . This is exemplified by the wide range of estimated peak pressures reported for eclogite xenoliths (~27-50 kbar) [31][32][33] , which necessarily must have been exhumed from the same paleo-slab top due to their spatially restricted occurrence in the northwestern NVF. However, over the past decade, petrological phase equilibrium modeling has become the preferred technique with which to perform reliable forward and inverse modeling of subduction zone processes 49,50 as it employs a multi-equilibrium approach combining internally consistent thermodynamic data and more geologically realistic a-x relations for minerals with solid solutions. The peak P-T conditions obtained here using this technique (~35-37 kbar and ~615-625 °C; Fig. 3 and Supplementary Figs 6 and 7) lie at the low-pressure end of the range defined by previous studies [31][32][33] , posing the question of what these alternative higher-pressure estimates mean. One possibility is that conventional thermobarometers in these studies were applied to minerals that were not in mutual chemical equilibrium at the time of metamorphism. For example, relatively minor retrograde diffusional equilibration of garnet rim compositions with adjacent matrix phases in mafic rocks has been demonstrated to cause true P-T conditions obtained via garnet-omphacite-phengite thermobarometry to vary by up to ±100 °C and ±4 kbar 51 , although this does not seem to be the case in 17MSR09, as the results of garnet-omphacite-phengite barometry match very well with the results obtained via inverse petrological modeling (Fig. 3a).
An alternative solution to these discrepancies is the influence of tectonic overpressure that complicates the conversion of pressure to depth within the Earth. Recent studies have shown that localized, non-lithostatic overpressure may occur to different degrees through a lithospheric column in convergent tectonic settings 52,53 . Further, relatively rigid rheological components preferentially act as foci for overpressure, such as cold and dense oceanic crust in direct contact with hotter and more malleable mantle lithosphere 54 . Modeled tectonic overpressure at the surface of subducted oceanic crust at ~120 km has been calculated to potentially reach magnitudes of 1-5 kbar 55 , yet this is too small to account for the absolute pressure differences reported between studies in this case. It is critical to note that any component of non-lithostatic overpressure within the total pressure calculated in our modeling procedure would have the effect of decreasing the depth of subduction of the Farallon slab, meaning that it would have penetrated through the proto-plateau's lithospheric root at even shallower depths and sheared away even larger mass of SCLM. Thus, the ~120 km reported here must represent a maximum depth, and that the calculated ~80 km of SCLM removed during its advance should thus be considered a minimum.
Prior study of mantle xenoliths from the NVF suggest that the hydrated SCLM beneath the Colorado Plateau has a density of ~3.34 g/cm 3 (Lee et al. 46 ) and thus is more buoyant than other SCLM beneath North America 56 , which has a density of ~3.38 g/cm 3 -similar to that forming the roots of many Archean cratons 57 (Fig. 4a). Most workers suggest that as the Farallon slab subducted northeast beneath the USA, fluid released from prograde devolatilization reactions -mainly driven by the consumption of chlorite, actinolite, talc and/or lawsonite 58promoted hydration of the overlying SCLM, promoting its chemical depletion via metasomatism 59,60 (Fig. 4a). Importantly, sheared-away SCLM would not have been exposed to this metasomatism, being positioned at depths below the advancing slab top (Fig. 4).
Thermomechanical modeling shows that a weaker (i.e. hydrated) SCLM that is more chemically depleted will be significantly more conducive to penetration by the Farallon slab during flat subduction 20 , which we suggest allowed this large-scale removal. The relatively dry underlying SCLM would then have a very slight negative buoyancy (Δρ = −0.01 g/cm 3 , assuming an asthenosphere density of ~3.37 g/cm 3 from Dziewonski and Anderson 61 ), and eventually sink into the underlying mantle (Fig. 4a). Our data therefore support the widely acknowledged hypothesis that lithospheric loss beneath the proto-Colorado Plateau during the Late Cretaceous and/or Eocene was driven by the inboard migration of the subducting Farallon slab 60 , leading to the atypically thin Colorado Plateau SCLM observed in geophysical profiles today 25 . These geophysical studies report a current lithospheric thickness between 120 km and 150 km, which directly correlates with our calculated estimate derived via thermobarometry, but conflicts with the results of some thermos-mechanical models, suggesting deep-seated subduction (Fig. 4b). In these cases, calculated metamorphic pressures of around 70 kbar would be expected from thermobarometric analysis, yet this is not the case for the samples investigated herein. Although the timing of removal of this mass of SCLM is uncertain, upper-mantle seismic tomography indicates that the Farallon slab began to break apart and sink in at least two stages at c. 86-60 Ma 62 and c.  . This shear-induced erosion and slab breakaway would then allow asthenospheric upwelling, which has been cited by some workers to be a critical factor in accounting for subsequent plateau uplift 16,23,24,62 . The proposed mechanism of lithospheric thinning by lateral shearing beneath the Colorado Plateau has been suggested by other studies 9,22 , and so our findings potentially provide the first direct constraints on this having occurred in North America.
What are the implications of this result for the efficacy of subduction erosion and removal of continental roots throughout Earth history? While the timing of the onset of subduction within the geological record is unresolved 63 , many studies agree that early plate tectonics involved shallow subduction within a hot Neoarchean-Paleoproterozoic mantle (c. 3.0-1.7 Ga) 64 , with occasional slab tearing and breakoff. This period of Earth history represents a convergence of many major geological transitions, including the widespread emergence of continental crust from beneath the oceans and the associated saturation of Earth's atmosphere with oxygen (the Great Oxygenation Event) 65 . Sub-horizontal subduction of the Farallon oceanic slab shows fundamental parallels with www.nature.com/scientificreports www.nature.com/scientificreports/ this early-Earth regime and may further account for the occurrence of some atypically thin Archean cratonic roots worldwide. Thermal estimates of continental lithospheric thickness show a consistent increase with age, from ~100 km in the Phanerozoic to ~250 km in the Early Proterozoic 66 . Archean cratonic lithosphere, however, has a bimodal distribution at ~350 km and ~220 km, despite thermal modeling suggesting an equilibrium thickness of >400 km 66 . Mechanisms such as thermo-mechanical erosion by secondary mantle convection 67 , erosion by mantle plumes 68 , delamination due to Rayleigh-Taylor-type gravitational instabilities in the lower lithosphere 69 , and erosion by basal drag 70 have all been suggested by previous workers, although our results suggest that the hotter ambient Neoarchean mantle may have alternatively promoted lithospheric shear-removal of these keels, akin to that documented here for the Colorado Plateau. As such, similarities between thinned SCLM in Neoarchean and Phanerozoic terranes provides support for these key geodynamic processes having operated similarly throughout much of geological time.

Methods
Petrological characterization of eclogite xenoliths. Electron microprobe analysis (EPMA). Majorelement compositional analyses of minerals in each xenolith were acquired using a JEOL JXA 8900 electron microprobe housed at the Denver Microbeam Laboratory at the United States Geological Survey (USGS), Colorado, USA. Both natural and synthetic materials were used as standards for calibration, and a ZAF correction routine was applied. Operating conditions comprised an acceleration voltage of 15 keV, a beam current of 20 nA, and a beam diameter of 5 μm for mica and 1 μm for all other minerals.
Automated mineralogy. Volume proportions of minerals in each xenolith were acquired via automated mineralogy, using a TESCAN VEGA-3 model LMU VP scanning electron microscope housed at the Department of Geology and Geological Engineering, Colorado School of Mines, USA. An electron beam was rastered across the surface of a thin section of each xenolith sample at a pixel resolution of 25 µm, using an acceleration voltage of 25 keV and beam intensity of 14.5 nA. Four energy-dispersive detectors simultaneously acquired a backscattered electron (BSE) image at each pixel and a chemical composition. Mineral characterization was performed using both BSE and compositional information, which were compared to spectra stored in an internal database. This procedure produced quantitative mineral abundance maps for each sample, with area proportions in thin section assumed to be representative of volume proportions throughout each xenolith. The calculated proportion of zoisite was considered to represent that for lawsonite at peak metamorphic conditions, as the fine-grained acicular zoisite pseudomorphs retain the original rhombohedral outlines of parent lawsonite, indicating direct replacement with no significant volume change during retrogression 31 . Automated mineralogy scans and obtained proportions for each eclogite xenolith are shown in Supplementary Figs 2-4. thermobarometry. Petrological modeling. Phase diagram construction was performed using the Gibbs free energy minimization software Theriak-Domino 71,72 and the internally consistent thermodynamic data set ds62 73 . Eclogite xenoliths 17GR11 and 17MSR11 were modeled in the nine-component Na 2 O-CaO-FeO-MgO-Al 2 O 3 -SiO 2 -H 2 O-TiO 2 -O 2 (NCFMASHTO) system, whereas petrological modeling of phengite-bearing xenolith 17MSR09 additionally considered K 2 O. The following a-x relations for solid-solution phases were used: clinopyroxene (diopside-omphacite-jadeite) and clinoamphibole (glaucophane-actinolite-hornblende) 74 ; garnet, biotite, muscovite-paragonite, and chlorite 75 ; epidote 73 ; plagioclase 76 ; and ilmenite 77 . Pure phases comprised talc, lawsonite, kyanite, zoisite, quartz, coesite, and rutile. Mineral abbreviations follow the guidelines of Whitney and Evans 78 . Effective bulk compositions for each xenolith were calculated using mineral proportions derived by automated mineralogy and representative EPMA-derived compositions 79 . Adjustments to measured mineral compositions were made using the "ideal analysis" approach of Powel and Holland 38 , where necessary, to reduce misfit between natural and modeled compositional systems. Individual bulk-rock compositions used to perform phase equilibrium modeling are shown in Supplementary Table 2.
Uncertainties related to the absolute positions of assemblage field boundaries on calculated phase diagrams have been shown to be less than ±1 kbar and ±50 °C at the 2σ level 38,39 , with this variation being largely a function of propagated uncertainty on end-member thermodynamic properties within the data set. However, as all phase diagrams were calculated using the same dataset and a-x relations, similar absolute errors associated with dataset end-members cancel, and calculated phase equilibria are relatively accurate to within ±0.2 kbar and ±10-15 °C 38,39 . The P-T conditions of peak metamorphism for each xenolith were determined by comparing mineral proportions calculated by automated mineralogy against those predicted in each petrological model ( Fig. 3 and Supplementary Figs 5 and 6).

Data Availability
All petrological data necessary to reproduce the results described herein are provided in Supplementary Information. Software enabling petrological calculations (Theriak-Domino) can be downloaded at no cost from the following web address: http://www.rocks.uni-kiel.de/theriakd/html/down_en.html.