Mantle Thermochemical Variations beneath the Continental United States Through Petrologic Interpretation of Seismic Tomography

The continental lithospheric mantle plays an essential role in stabilizing continents over long geological time scales. Quantifying spatial variations in compositional and thermochemical properties of the mantle lithosphere is crucial to understanding its formation and its impact on continental stability; however, our understanding of these variations remains limited. Here we apply the Whole-rock Interpretive Seismic Toolbox For Ultramafic Lithologies (WISTFUL) to estimate thermal, compositional, and density variations in the continental mantle beneath the contiguous United States from MITPS 20, a joint body and surface wave tomographic inversion for Vp and Vs with high resolution in the shallow mantle (60-100 km). Our analysis shows lateral variations in temperature beneath the continental United States of up to 800–900°C at 60, 80, and 100 km depth. East of the Rocky Mountains, the mantle lithosphere is generally cold (350–850°C at 60 km), with higher temperatures (up to 1000°C at 60 km) along the Atlantic coastal margin. By contrast, the mantle lithosphere west of the Rocky Mountains is hot (typically >1000°C at 60 km, >1200°C at 80–100 km), with the highest temperatures beneath Holocene volcanoes. In agreement with previous work, we find that the predicted chemical depletion does not fully offset the density difference due to temperature. Extending our results using Rayleigh-Taylor instability analysis, implies the lithosphere below the United States could be undergoing oscillatory convection, in which cooling, densification, and sinking of a chemically buoyant layer alternates with reheating and rising of that layer. Supplemental Information for Mantle Thermochemical Variations beneath the Continental United States Through Petrologic Interpretation of Seismic Tomography William J. Shinevar, Eva M. Golos, Oliver Jagoutz, Mark D. Behn, Robert Van der Hilst 1 MIT/WHOI Joint Program in Oceanography/Applied Ocean Engineering 2 Brown University 3 Massachusetts Institute of Technology 4 Boston College now at University of Colorado Boulder This supplement contains two supplemental tables containing compiled xenolith and primary magma thermobarometry, two supplemental figures, and a data file containing the results plotted in Figures 3-6.

The first potential cause for the systematic difference in our results and magmatic temperature estimates is the scale over which each is measuring. Magmatic temperature estimates sample the temperature at which the melt was in equilibrium with the mantle and may only represent small regions (˜1 m) of the mantle, especially if melt is focused (e.g., Kelemen and Dick, 1995) and/or escapes the mantle on short (˜10 kyr) timescales (Feineman and DePaolo, 2003). Conversely, seismic tomographic inversions such as MITPS 20 are limited to resolving seismic anomalies greater than˜1.5°x1.5°with vertical resolution on the order of 10 km . Thus, non-pervasive, small-scale seismic anomalies due to thermal upwellings or melt may be smeared or not sensed. Furthermore, seismic inversions smooth their wave speeds in order to stabilize the inversion, though MITPS 20 corrects for some of this effect (see Methodology). Subduction zones are especially difficult to image due to any smearing of the cold, subducting lithosphere, which increases seismic wave speed, decreasing the temperature estimate. The fact that the 80-km temperatures are in better agreement with magmatic estimates may suggest that at 60 km, vertical smearing may increase the inverted seismic wave speeds as the tomography samples from starkly colder lithosphere along a steep geotherm. Regional, high-resolution seismic studies are necessary to understand these effects.
A second reason for the temperature discrepancy could be error in anelastic corrections of seismic wave speed. Anelasticity experiments on olivine and peridotite are difficult, with various experimental groups giving different results and sensitivities (Faul and Jackson, 2015;Karato and Park, 2018). The Behn et al. (2009) power-law formulation of anelasticity does not fit experimental data well at low quality factors (Q -1 >0.1, high temperatures or melt present. Jackson and Faul, 2010). Certain parameters in the both anelasticities are relatively unconstrained, like the activation volume that controls the pressure sensitivity (Faul and Jackson, 2015;Jackson and Faul, 2010). Other comparisons of high-quality seismic experiments and forward calculations of peridotite seismic wave speeds required altering the relaxation peak of the frequencies in order for the observations to be interpreted by the Jackson and Faul (2010) anelasticity (Ma et al., 2020). Furthermore, the effect of water content on anelasticity is currently debated (Aizawa et al., 2008; Cline et al., 2018; Karato and Park, 2018). Increasing the water content decreases V s at high temperatures, thus shifting all forward calculations above˜900°C to the left in Figure 2. As we assumed relatively dry water contents (C OH =50 H/10 6 Si), assuming an increased water content would decrease the interpreted temperatures. While grain size is an important parameter for anelasticities (Behn et al., 2009;Faul and Jackson, 2005), we have assumed a reasonable grain size near the upper bound observed in xenoliths. Any grain size reduction increases anelastic effects, thus reducing the modeled temperature. Therefore, variable grain size and its effect on anelasticity cannot reconcile the temperature discrepancy discussed here. Oxidation has been found to increase dissipation (Cline et al., 2018), not incorporated in our methods. This would also reduce our temperature results in arc settings as arc mantle is more oxidized (Kelley and Cottrell, 2012). Conversely, as long as melt and fluids are focused, oxidation would not decrease large-scale seismic wave speed as only a small portion of the mantle may be highly oxidized.
While increasing water content can drastically reduce the peridotite solidus, tholeiitic (dry) magmas are observed in the western United States (Till, 2017). Melting experiments on dry peridotite compositions require temperatures >1300°C at 60 km depth, greater than nearly all our temperature results (Hirschmann, 2000). Given the existence of Holocene age tholeiitic magmas, the reported uncertainty in magmatic thermobarometry (11-43°C, 0.1-0.4 GPa) cannot explain the temperature discrepancy present at 60 km.
Miscalculation in the forward calculation of seismic wave speeds in WISTFUL is also unlikely to explain the temperature estimate discrepancy. As we incorporate expected non-systematic uncertainty from our forward calculations into the error allowed for fitting and utilize current experimental moduli for most mineral endmembers, the only systematic error from this could be due to a difference in mixing assumptions, e.g. anisotropy. WISTFUL calculates the isotropic wave speeds, so comparing with the fast direction wave speeds would produce colder than realistic temperature estimates. As MITPS 20 inverts isotropic wave speed from combination of teleseismic body waves and surface wave arrival times, a systematic increase in recovered wave speed due to anisotropy beneath all regions with magmatic temperature estimates is unlikely.
Lastly, WISTFUL does not incorporate any effect of melt and hydrous phases, both of which would decrease the predicted temperature. Melt strongly reduces V s (e.g., Hammond and Humphreys, 2000), but the exact wave speed reduction is heavily dependent on the melt content and the melt connectivity (Zhu et al., 2011). Thus, incorporating the effect of melt would make supersolidus temperatures require even lower wave speeds than observed. Similarly, pargasitic amphibole, (NaCa 2 (Mg 4 Al)(Si 6 Al 2 )O 22 (OH) 2 ), the most common hydrous phase predicted for the shallow mantle (Dawson and Smith, 1982), has V s =3.85 km s -1 andV p /V s =1.83 at 60 km pressure and 1200°C assuming the same anelasticity described in Section 3 and moduli from Abers and Hacker (2016). At the same conditions, diopside (MgCaSi 2 O 6 ) hasV s =4.35 km s -1 andV p /V s =1.78. Replacing clinopyroxene at the same temperature with pargasite would decreaseV s and increaseV p /V s (shifting all forward calculations to the upper left in Figure 2). Therefore, the addition of melt and/or hydrous phases would decrease the predicted temperature for the same seismic wave speed. Despite amphibole dikes and peridotites being hypothesized as a significant source of volatiles for alkaline and ocean-island basalts (Harry and Leeman, 1995;Pilet et al., 2011), the volume required to have a geochemical impact is unlikely to have a noticeable impact on seismic wave speed of the shallow mantle on the scale we are can interpret with MITPS 20 (˜1.5°).
In summary, our results agree within error of recent (<10 Ma) xenolith compositions, predict temperature greater or equal temperature to spinel-bearing and garnet-bearing xenoliths, but underpredict magmatic temperature estimates, especially at 60 km. The systematic difference between our best-fit temperature estimates with magmatic temperature estimates are best explained by a difference in the scale of the estimates, smearing in the tomographic models at shallow depths along a steep geotherm, and/or errors in anelastic corrections. Further experimental work on anelasticity is required to better interpret high temperature mantle regions like the western United States. Vs from MITPS 20 at 60 km Vs80

Supplementary
Vs from MITPS 20 at 80 km Vs100 Vs from MITPS 20 at 100 km Figure S1: Temperature estimates from magmatic temperature estimates against average WISTFUL temperature within 0.5°arcdistance of the surface outcrop from Figure S3 for 60 (blue) and 80 km (green). Error bars depict a 1-sigma error. Figure S2: Compositional buoyancy as defined in section 6.2, equation 7. Boundaries as in Figure S1. States, is interpreted in terms of temperature, composition, and density. 34

References
• Our method predicts temperatures of 260-1430°C, compositions of Mg# 85-92, and 35 density between 3230-3370 kg m -3 between 60-100 km. 36 • These results imply that the mantle lithosphere has enough compositional buoyancy to 37 compensate for half the negative thermal buoyancy.

Introduction 39
The North American continent consists of amalgamated continental and arc fragments 40 originating over the course of Earth's history (Whitmeyer and Karlstrom, 2007  Godson, 1989). The coherence between Bouguer anomalies and topography predicts smaller 50 elastic thicknesses in the west than in the east (Steinberger and Becker, 2018), which is 51 consistent with the inference from seismology that the depth at which mantle seismic anisotropy 52 aligns with absolute plate motion is shallower in the west (~80 km) than in the east (~200 km) 53 (Yuan and Romanowicz, 2010). Mantle seismic wave speeds also differ between the two regions, 54 with slower wave speeds in the west . 55 In the upper mantle, temperature exerts the dominant control on rheology (Hirth and 56 Kohlstedt, 2003), density, and seismic wave speed (Shinevar et al., 2022). Thus, taken together 57 these observations suggest elevated shallow mantle temperatures in the western US compared to 58 those in the east (e.g., Goes and van der Lee, 2002). This variation in mantle temperature has 59 been related to the age of the lithosphere as the lithosphere cools and thickens with time 60 (Mareschal and Jaupart, 2013).
Re-Os dating and isotopic measurements suggest that the shallow cratonic mantle is as 62 old as the surface rocks (e.g., Pearson, 1999). In order for cratonic mantle to remain buoyant and 63 persist over billions of years, the density increase due to cooling has been hypothesized to be 64 balanced by a density decrease due to the depletion of the mantle through melting-the so-called 65 isopycnic hypothesis (Jordan, 1975). The density structure associated with the combined effects 66 of composition and temperature is therefore vital to the force balance within the North American and databases chosen to best-fit the mineral modes of well-studied mantle xenoliths, and 83 experimental calibrations of olivine anelasticity. In this study, we first briefly discuss the methodology behind MITPS_20 and WISTFUL. We then present maps of inferred temperature, 85 composition, and density for the continental United States and eastern Canada at 60, 80, and 100 86 km depth, where the MITPS_20 model is best constrained by data. To validate our methodology, 87 we compare our results at 60 and 80 km with estimates of temperature and composition from 88 recently erupted xenoliths and magmatism. Using the best-fit density, composition, and 89 temperatures, we then investigate the relative chemical and thermal buoyancy of continental 90 lithosphere. We find an imbalance in these buoyancy terms, suggesting that the continental 91 lithosphere is density unstable. We argue that these observations could be the result of oscillatory 92 convection, in which cooling, densification, and sinking of a chemically buoyant layer alternates 93 with reheating and rising-resulting in laterally harmonic perturbations to the interface between 94 the layers. 95  (Atwater, 1989). Since cessation of Laramide compression, the 100 Basin and Range has undergone large-scale extension (Parsons, 2006) and the Snake River 101

Geological Setting and Previous Work
Plain/Yellowstone Plateau was formed by the impingement of the Yellowstone Plume beginning 102 at 16 Ma (Leeman, 1982). In contrast, the plate interior last experienced internal deformation 103 during the Neoproterozoic due to the Mid-Continent Rift (1100 Ma) and Grenville Orogeny (ca.  (1) (2) After scaling, the relative variations are converted to absolute and ( / ) ( Figure  170 1a, b) using 171 Where and are from a modified version of the ak135 reference model (Kennett et al., 172 1995).
/ for ak135 is greater than interpretable values obtained from WISTFUL (e.g., 173 ak135 Vp/Vs is 1.793 at 60 km, whereas the range produced by WISTFUL is 1.73-1.77, Figure  174 2). Therefore  In order to interpret seismic wave speeds under mantle conditions, we correct for the 217 anelastic behavior of olivine at high temperatures using the power-law formulation of Behn et al. 218 (2009). We assume that olivine anelasticity applies to all minerals present. We apply anelastic 219 corrections using the periods that dominate surface wave sensitivity for MITPS_20 (38 s for 60 220 km, 45 s for 80 km, and 57s for 100 km). We assume a grain size of 1 cm, in line with grain sizes 221 for cratonic xenoliths (e.g., Baptiste    the probability that the relationship is due to randomness. This correlation between mantle 256 density and topography is in agreement with the hypothesis that the western United States is near 257 isostatic equilibrium and that mantle density variations due to temperature partially support 258 topography (e.g., Molnar et al., 2015). It is important to note that the temperature errors are small compared to the total predicted variation (60°C error for a variation of ~1000°C), but the errors 260 for Mg# (1.0 for a 3-7 variation) and density (20 kg m -3 for a 110-140 kg m -3 variation) can be a 261 significant fraction of the total variations. 262   To validate our approach, we compare our results with estimates of temperature and 286 composition from xenoliths, as well as melting P-T estimates calculated from primary magmas. 287 The benefit of comparison with well-studied young (<10 Ma) xenoliths is that they directly 288 constrain composition. In addition to composition, estimates of mantle temperature have been 289 calculated from xenolith using appropriate thermobarometry. Comparison with thermobarometry 290 using primary magma compositions is discussed in Supporting Information and Supplemental 291

Xenolith Compositions 296
We compiled 15 xenolith localities from recently erupted volcanos (<10 Ma, black 297 squares, Figures 5, 6), and which have at least four xenolith bulk rock compositions with 298 Mg#>85 and <2 wt % loss on ignition to avoid averaging xenoliths that have been refertilized or 299 altered before/during eruptive processes (references in Supplementary Table 1). We assume that 300 the average xenolith composition represents the mantle composition beneath that location and bearing xenolith localities for which the 60 km depth slice was used; the green square signifies a 314 garnet-bearing xenolith locality for which the 80 km slice was used. b) Temperature estimates 315 from xenoliths against average WISTFUL temperature within 0.5° arc distance of the surface 316 outcrop for 60 (blue) and 80 km (green). Error bars depict a 1-σ error. 317

Xenolith Thermobarometry 318
Mantle xenolith thermobarometry relies on using the relative mineral composition to 319 calculate the equilibrium P-T. We compiled 13 localities with at least one temperature (or P-T) 320 estimate (references in Supplementary Table 1). Because no reliable spinel barometer exists, 321 most spinel thermometry is calculated at 1.5 GPa as spinel-bearing xenoliths could originate 322 submitted to Earth and Planetary Science Letters anywhere in the spinel stability field (~0.7-2 GPa). Here, we compare all spinel thermometers 323 with our 60-km estimate (~1.7 GPa) and acknowledge that our temperature estimates should be 324 equal to or exceed the estimates from this thermometry as the spinel xenoliths could be sampling 325 a shallower mantle. Reliable barometers exist for garnet-bearing peridotites, and we compare all 326 temperatures from garnet-bearing peridotites within 0.3 GPa of the 80 km pressure (~2.4 GPa). 327 To compare these temperature estimates with ours, as above, we average the best-fit 328 temperatures from all results within 0.5° arc-distance from the locality. As with composition, we 329 take the WISTFUL temperature error to be the maximum of the average temperature uncertainty 330 and the standard deviation of averaged temperatures. We consider the temperature error in a 331 xenolith locality to be the maximum of the standard deviation of the calculated temperatures and 332 published thermometer uncertainty. 333 Our temperature estimates agree within error for 4 out of 11 spinel xenolith localities, 334 while overpredicting temperature by ~125°C for five xenolith localities (Figure 7b). Xenoliths 335 from Green Knobs in New Mexico and Vulcan's Throne in Arizona predict much lower 336 temperatures (750-775°C) at 60 km than our results (~1000°C). These values appear to be 337 anomalously low as there is recent (<1 Myr) magmatism nearby and recently exhumed granulite 338 facies lower crustal xenoliths (Cipar et al., 2020). It is therefore possible that these spinel 339 xenoliths sample significantly shallower mantle, potentially as shallow as the regional Moho 340 (~30 km, Schmandt et al., 2015). 341 Our results disagree within error for both garnet xenolith localities (green triangles, 342 Figure 7b). The more enigmatic garnet-bearing xenolith locality (Big Creek, Sierra Nevadas, 343 California) predicts low temperatures for 80 km (~750°C), in stark disagreement with our results. 344 It is unlikely that these xenoliths were not in equilibrium at high temperatures (~750°C). Instead, the xenoliths may record a no-longer-present thermochemical state, as they erupted in an 8 Myr 346 old diatreme (Chin et al., 2012) in a region hypothesized to be undergoing continental 347 delamination between 10-3 Ma (Zandt et al., 2004). 348

Discussion 349
Here we discuss regions with anomalous temperature and composition to compare with 350 regional geology and tectonic history. Subsequently, we present estimates of lithospheric 351 buoyancy, and consider implications for our understanding of cratonic lithospheric stability. 352 given the low regional isostatic gravity anomalies (Molnar et al., 2015). The excess topography 449 of the Rocky Mountains is instead more likely due to a steep decline in density (~30 kg m -3 , 450 ~1%), observed here beginning at the eastern front of the Rocky Mountains (~255° E, Figure 6), 451 as previously hypothesized and inferred (e.g., Levandowski et al., 2014). 452

Lithospheric Buoyancy and Stability 453
As WISTFUL predicts temperature, density, and composition, we have the unique 454 opportunity to investigate the relative importance of compositional and thermal buoyancy for the 455 stability of the mantle lithosphere beneath the entire continental United States. To investigate 456 this, we calculate a dimensionless buoyancy number, B, the ratio of the intrinsic (compositional) 457 and thermal buoyancies, 458 where Δ is the density variation attributed to compositional heterogeneity and Δ is the 459 change of density due to temperature differences (Shapiro et al., 1999b). Negative B values imply negative compositional buoyancy. A B value of 0 implies that no compositional buoyancy 461 exists. A B value of 1 implies that the compositional and thermal effects are equal (the isopycnic 462 hypothesis). A B value much greater than 1 implies that the compositional buoyancy is greater 463 than thermal effects. Thus, when B>1, the compositional effect on density is sufficient for the 464 lithosphere to remain positively buoyant, but when B<1 there will not be sufficient 465 compositional buoyancy to overcome the negative thermal buoyancy of the cold lithospheric 466 mantle. 467 We calculate B using the WISTFUL predictions of density for the compositions that are 468 within error seismically at the best-fit temperature (Table 1). We define 469 where is the average density of the best-fit compositions at the same pressure and a mantle 470 potential temperature (1350°C) weighted by the inverse of the seismic error at the best fit-471 temperature, and is the density calculated for depleted MORB mantle (DMM, Mg# 89.4, 472 Workman and Hart, 2005) at the same conditions as and utilizing the same thermodynamic 473 calculations as WISTFUL (3245 kg m -3 at 60 km, 3280 kg m -3 at 80 km, and 3305 kg m -3 at 100 474 km). We define 475 where the best-fit density at the best-fit temperature. At high temperatures approaching 476 the mantle potential temperature, Δ becomes very small, which results in extremely large 477 values of B. We therefore only calculate B at temperatures less than 1300°C (50°C below the 478 reference potential temperature). 479 Mountain Front compared to the mantle in the cratonic portions or east of the Grenville Front, as 517 subsolidus oscillatory convective mixing would fertilize the shallow depleted mantle. As the 518 LAB becomes significantly deeper as seen in the cratonic mantle (>200 km, Steinberger and 519 Becker, 2018; Yuan and Romanowicz, 2010), oscillatory convection may eventually become too 520 slow due to increasing mantle viscosity with pressure (Hirth and Kohlstedt, 2003). This process 521 would predict an overall decrease of cratonic mantle age with increasing depth, but also predicts 522 submitted to Earth and Planetary Science Letters the presence of some older outliers at depth due to convective mixing, as observed (Pearson, 523 1999