In this work, we compile several seismic velocity models publicly available from the Incorporated Research Institute for Seismology (IRIS) Earth Model Collaboration (EMC) and compare subcrustal mantle velocities in the models to each other and to the timing of tectonism across the continent. This work allows us to assess the relationship between the time elapsed since the most recent thermotectonic event and uppermost mantle temperatures. We apply mineral- and physics-based models of velocity-temperature relationships to calculate upper-mantle temperatures in order to determine cooling rates for the lower-crust and uppermost mantle following thermotectonic activity. Results show that most of the cooling occurs in the ∼300–500 million years following orogeny. This work summarizes current estimates of upper-mantle shear velocities and provides insights on the thermal stabilization of continental lithosphere through time.

One of the most significant products of the EarthScope facility has been the development of new tomography models that take advantage of the consistent station design, regular 70 km station spacing, and wide aperture of the EarthScope Transportable Array (TA) network. These models have sharpened our knowledge of known mantle structure and led to the discovery and interpretation of additional compositional, thermal, and density anomalies throughout the continental United States, especially in regions previously thought to be tectonically quiescent (e.g., Biryol et al., 2016; Menke et al., 2016; Wolin, 2016; Mazza et al., 2017). Discoveries from the EarthScope experiment have led to a newfound appreciation for the complexity of the deep crust and mantle lithosphere and the role of these regions in influencing topography, surface deformation, and regional tectonic processes at all levels of the lithosphere.

The wide aperture of the EarthScope experiment, when combined with other publicly available data sets, allows us to better understand the long-term effects of orogeny on the continent (e.g., Porter et al., 2016). It is well established that orogenic events have profoundly shaped the surface topography and crustal geology of the continent. However, less is known about the long-term effects of these tectonic events on the state of the lithospheric mantle and the thermal state of the lithosphere. By compiling multiple data sources including seismic velocities and geochemical data, we are able to assess the response of the lithospheric mantle to tectonism at continental scales.

Tectonism across the Continent

The North American continent has a long history of tectonism beginning with the Proterozoic assembly of Laurentia, continuing with Paleozoic collisions in eastern North America during the assembly of Pangea, and extending to modern deformation and volcanism along the western margin of the continent. These tectonic events led to deformation throughout the lithosphere and significantly altered the thermal structure of the mantle, typically by bringing hot asthenosphere to the base of the lithospheric mantle, which warms the lithosphere and the overlying crust. Because thermal conditions have been shown to dominate seismic velocities within the mantle, tectonism is commonly associated with reduced seismic velocities due to heating (see discussion below). In order to place mantle seismic-velocity observations calculated from EarthScope’s USArray data in the context of continental tectonism, we briefly summarize and discuss the most recent orogenic events affecting the continental United States by region.

The oldest basement rocks in the continental United States are located between the Rocky and Appalachian mountain chains, where they represent the southern extent of the Archean–Proterozoic North American craton (Whitmeyer and Karlstrom, 2007). Archean cratons (e.g., Superior and Wyoming) stabilized and assembled into the cratonic core of Laurentia (Houston et al., 1993; Henry et al., 2000; Chamberlain et al., 2003) during Paleoproterozoic orogenic events that formed the supercontinent of Nuna (Hoffman, 1988), also known as Columbia (Rogers and Santosh, 2002). This north-central region of the United States was largely stabilized within tens of millions of years after the Trans-Hudson orogenic cycle, at ca. 1.82 Ga (Lewry and Collerson, 1990; Corrigan et al., 2009).

A series of juvenile terranes were sequentially emplaced onto the southern and eastern margins of the North American cratonic core in the Mesoproterozoic (Whitmeyer and Karlstrom, 2007). These terranes make up the present-day Yavapai, Mazatzal, and Granite Rhyolite (Shawnee) provinces. The Yavapai province, which extends NE-SW from the Great Lakes to Arizona-California, was emplaced by the end of the Yavapai orogeny at ca. 1.68 Ga (Bowring and Karlstrom, 1990). This block is juxtaposed against the Mojave terrane in southern Arizona and against the Wyoming craton at the Cheyenne belt in southern Wyoming (Karlstrom and Bowring, 1988; Crosswhite and Humphreys, 2003; Duebendorfer et al., 2006; Lund et al., 2015). The Mazatzal province parallels the southeastern boundary of the Yavapai province, extending from Labrador-Newfoundland in Canada to southern Arizona in the United States (Karlstrom and Bowring, 1988; Shaw and Karlstrom, 1999) and was emplaced by ca. 1.64 Ga during a complex series of continental collisions and suturing events (Holm et al., 1998; Romano et al., 2000; Amato et al., 2008). The Granite-Rhyolite (Shawnee) province was emplaced by 1.35 Ga and lies to the east of the Yavapai and Mazatzal provinces (Whitmeyer and Karlstrom, 2007). Significant portions of the region that extends from the Laramide front in the western United States to the Appalachian Mountains in the east were overprinted by deformation during the 1.3–1.0 Ga Grenville orogenic cycle (Donaldson and Irving, 1972; Eriksson et al., 2003; Rainbird et al., 2012; Craddock et al., 2017).

Subsequent rifting events include the ca. 1.1 Ga Midcontinent rift, which cuts across the southern Superior Province, the Penokean orogeny and Marshfield terrane, and the Yavapai and Mazatzal provinces (Green, 1983; Van Schmus, 1992; Miller and Nicholson, 2013). The Midcontinent rift produced extensive but narrow bands of dense volcanic rocks in the upper crust, as well as associated underplated material between the crust and the mantle (Tréhu et al., 1991; Zhang et al., 2016). Younger rifting events include the early Paleozoic Reelfoot rift (Ervin and McGinnis, 1975; Mooney et al., 1983; Thomas, 1991; Johnson et al., 1994) and the formation of the Oklahoma aulacogen (Gilbert, 1983; Keller and Baldridge, 2006).

Paleozoic orogenic events include the Pennsylvanian–Permian formation of the ancestral Rocky Mountains (Kluth and Coney, 1981; Kluth, 1986; Leary et al., 2017). These orogenic events impacted the south-central United States and were followed by a series of contractional events (in the western United States) that initiated in the Mesozoic with the Sevier orogeny (DeCelles, 2004; Dickinson, 2004) and continued through the Laramide orogeny until ca. 40 Ma (Grove et al., 2003; Ducea et al., 2015). Post-Laramide extensional tectonics in the western United States include formation of the Basin and Range province (Spencer and Reynolds, 1990; McQuarrie and Wernicke, 2005) and other structures, such as the Rio Grande rift (Ricketts et al., 2014). Additionally, modern tectonism associated with upper-mantle flow due to lithospheric density instabilities, edge convection, or active upwellings has been hypothesized at several locales within the region (e.g., Schutt and Dueker, 2008; West et al., 2009; van Wijk et al., 2010; James et al., 2011; Allison et al., 2013; Huang et al., 2015).

Eastern and southeastern North America experienced a protracted series of orogenic events during the middle to late Paleozoic, often categorized as the Taconic, Acadian, and Alleghenian orogenies (Rodgers, 1967; Hatcher, 2005). This sequence, collectively called the Appalachian orogenic cycle, completed the assembly of the Pangea supercontinent and formed the Appalachian Mountains. These events were driven by the collision of Laurentia with continental fragments, island arcs, and ultimately, the African and Eurasian continents (Williams and Hatcher, 1982; Dalziel et al., 1994). The breakup of Pangea at roughly 175 Ma resulted in extension, rifting, and ocean formation along the eastern U.S. margin. Recent geophysical and geochemical work has also shown that localized areas along the eastern margin of the continent may be tectonically active in the present day (van der Lee et al., 2008; Mazza et al., 2014; Schmandt and Lin, 2014; Menke et al., 2016).

It is important to note that different styles of tectonism may have different effects on the thermal state of continental lithosphere. For example, transform tectonism does not influence the thermal state of the lithosphere as significantly as rifting, subduction, or continental collision. Here we do not differentiate between styles of tectonism. Because we use geochemical ages primarily from volcanic and igneous rocks as a proxy for tectonism, our age estimates are confined to events that produce these types of rocks, which are typically associated with lithospheric heating.

Factors Controlling Seismic Velocities

Within the upper mantle, seismic velocities are controlled by a combination of temperature, pressure, grain size, the presence of partial melt and/or fluids, and composition (summarized in Table 1). Each of these factors has varying impacts on the seismic velocities measured within the Earth, and isolating any one is challenging (van der Lee and Wiens, 2006; Schutt and Lesher, 2006). Previous seismic and mineral-physics work has quantified many of the effects that change in these variables would have on observed seismic velocities, and these works provide a framework for our calculations of temperature from shear velocity measurements (Goes et al., 2000; Cammarano et al., 2003; Faul and Jackson, 2005; Priestley and Mckenzie, 2006, 2013; Jackson et al., 2008; Jackson and Faul, 2010; McCarthy et al., 2011).

Previous work has shown that, in the absence of significant quantities of partial melts or fluids, temperature is a dominant factor in controlling upper-mantle seismic velocities at a given pressure (e.g., Goes et al., 2000; Cammarano et al., 2003; Faul and Jackson, 2005; Jackson and Faul, 2010), including within North America (Goes and van der Lee, 2002). Elevated temperatures lower elastic moduli within materials, and thus temperature is inversely correlated with seismic velocity. The effects of temperature variation stem from both anharmonic and anelastic causes. In most approximations of seismic velocity, anharmonic drops in the elastic moduli are linearly related to increases in temperature and are commonly calculated relative to a reference temperature and seismic velocity. The anelastic component reflects the dissipation of elastic energy from wave propagation into attenuation, its inverse represented by a quality factor (Q), and is controlled by the activation energy, activation volume, grain size, temperature, and pressure of a material. Its relationship to temperature is nonlinear, and at temperatures below ∼900 °C, mantle rocks typically fall within the elastic regime where anelasticity is not a significant factor in controlling seismic velocities (Goes et al., 2000; Faul and Jackson, 2005). At higher temperatures, anelasticity has an increasing impact on mantle seismic velocities.

The estimated effects of temperature variations on mantle seismic velocities vary with depth (pressure) and temperature, but at 50 km depth, Goes et al. (2000) calculate an ∼100 m/s change in seismic velocity for every 100 °C change in temperature. Within the North American upper mantle, temperatures have previously been estimated to vary by many hundreds of degrees Celsius between the eastern and western United States. Much of the variation in seismic velocity across the continent can be explained by these variations in temperature (Black and Braile, 1982; Goes and van der Lee, 2002; Schutt et al., 2011; Kaban et al., 2014). Given its large impact on mantle seismic velocities relative to other factors, temperature variations can be approximated from variations in absolute P- and S-wave velocities. In this work, we use S-wave velocities to estimate mantle temperatures.

Grain size strongly influences anelasticty and therefore seismic velocities within the mantle. Larger grain sizes are typically associated with faster seismic velocities (e.g., Faul and Jackson, 2005). However, the impact of grain size is challenging to calculate, because the relative effects of grain size on seismic velocity are influenced both by the frequency of the propagating wave and the temperature of the mantle (Jackson and Faul, 2010), resulting in a nonlinear relationship between grain size and seismic velocity. At higher temperatures, when anelastic effects are greater, grain size is thought to have a larger influence on velocity than at cooler temperatures (Jackson and Faul, 2010). Complicating our understanding of the relationship between grain size and mantle velocities is uncertainty in the actual grain sizes within the mantle, which can vary by several orders of magnitude (e.g., Hirth and Kohlstedt, 2003). Mantle grain sizes are partially controlled by stress, and regions of higher stress typically exhibit smaller grain sizes than lower-stress regions (Hall and Parmentier, 2003; Bürgmann and Dresen, 2008). In our calculations, we utilize grain-size estimates of 1 mm for our temperature estimates for all models, while exploring the effects of other grain sizes with the Jackson and Faul model (2010). For a more detailed discussion of mantle grain size, refer to Hirth and Kohlstedt (2003).

The presence of partial melt within the mantle is associated with a decrease in seismic velocity. However, it is hard to quantify this effect, because it not only depends on the percentage of melt but also on the mineralogy and melt geometry. Therefore, seismic velocities depend nonlinearly on melt percentages (Hammond and Humphreys, 2000), and the form of this nonlinear relationship is debated. While effects of large quantities of melt are substantial, the effects of small amounts of melt on shear moduli are debated by different authors (Anderson and Sammis, 1970; van der Lee, 2002; Dunn and Forsyth, 2003; Kreutzmann et al., 2004; Schutt and Dueker, 2008; Priestley and Mckenzie, 2013). We do not account for melt in our calculations and note that, in regions of high temperature where melt may be present, our temperature estimates likely represent a maximum possible temperature.

The presence of water is also associated with a decrease in seismic velocity. Although water is likely abundant in the uppermost mantle beneath the western United States (Dixon et al., 2004), constraining its effects on seismic velocity is challenging, and the influence of water on seismic velocities of the lithosphere-asthenosphere system is debated (Karato, 1995, 2003; Karato and Jung, 1998; Aizawa et al., 2007; Cline et al., 2018). Given the uncertainty, we do not account for hydration effects in our analysis.

The impact of composition on seismic velocities is thought to be secondary relative to temperature and the other factors listed above. Melt extraction (depletion) has been shown to have minor effects on density and little to no effect on seismic velocities for reasonable mantle compositions (<1%) (Cammarano et al., 2003; Priestley and Mckenzie, 2006; Schutt and Lesher, 2006). A limited exception to this is in regions where serpentinite is thought to be stable (i.e., subduction zone forearcs). Serpentinite has a very low shear velocity, which can produce decreased amplitude Moho conversions or even “inverted Mohos” in scattered wave imaging. An “inverted Moho” signifies that the lower crust exhibits higher seismic velocities than the underlying mantle (e.g., Bostock et al., 2002; Hyndman and Peacock, 2003). However, these conditions are only likely in a few places within the continent (i.e., the Cascadia forearc). Given its typically small impact on bulk velocities relative to the uncertainties in tomographic models and in the relationships used to convert velocity to temperature, we also do not account for compositional changes in our models.

The EarthScope program has led to the generation of an ensemble of seismic-velocity models at a variety of spatial scales, ranging from small-scale regional or local models that address a focused tectonic problem to large-scale, continent-wide models that span the breadth of geology described above. Here we utilize continent-scale velocity models that provide estimates of absolute velocity. The wide aperture for these models is useful, because it leaves no gaps beneath non-targeted regions and minimizes artifacts from variations in data regularization and processing techniques that would arise from a mosaic of multiple regional-scale models. Absolute velocity measurements can readily be tied to mineral-physics–based temperature models. Based on these criteria, we utilize shear-velocity models inferred from the inversion of earthquake and ambient noise shear waves with crustal thickness constrained by receiver functions (Schmandt and Lin, 2015; Porter et al., 2016; Shen and Ritzwoller, 2016). For contrast, we also include a pre-EarthScope shear-velocity model inferred from full-waveform fitting of earthquake shear waves and additional constraints from receiver functions (van der Lee and Frederiksen, 2005). For the data processing and regularization utilized in each model, the reader is referred to original works.

Because the goal of this work is to examine the thermal evolution of the mantle lithosphere, we utilize average velocities calculated between 6 and 12 km beneath the Moho as defined in each of the shear-velocity models and at 100 km depth for the pre-EarthScope model. These model data were downloaded from the IRIS Earth Model Collaboration (EMC), and velocities from 6 to 12 km below each model-defined Moho were extracted from the EarthScope models (Figs. 1 and 2) (IRIS DMC, 2011; Trabant et al., 2012). Mean upper-mantle velocities for the tectonic provinces of the continental United States were plotted to show the relationship between mantle velocities and physiographic provinces in each model (Fig. 2). We then interpolated the three EarthScope shear-wave velocity models to a uniform spacing and calculated the average and standard deviation of these velocities across the footprint of EarthScope’s USArray (Fig. 3). Given that pressure, which impacts the relationship between seismic velocities and temperature, is dependent on depth, the mean Moho depth from the three models and the standard deviation were also calculated (Fig. 3). Though the regularization schemes from the three different models vary, this exercise allows us to focus on regions of similarity in the three models, which are inferred to have the most robust velocity and crustal thickness data, and therefore provide the best constraints on lithospheric evolution.

We utilize age data from the Interdisciplinary Earth Data Alliance (IEDA): Earthchem Database to produce a model of most recent thermotectonic activity for the continental United States (Fig. 4). While this database is not a complete record of tectonic activity and sampling biases exist within the database, it still provides a reasonable approximation of the ages of tectonic activity across the continent. With a few exceptions, primarily at the model edges, these ages agree well with ages from the orogenic events described above, as well as with global-scale estimates of age (Artemieva, 2009).

To calculate an age model from these data, we first downloaded all geolocated and dated sample data from the IEDA database from within the footprint of USArray and surrounding regions and culled the data set to remove data points deemed non-representative of thermotectonic activity, leaving behind data from locally derived metamorphic, volcanic, and igneous rocks. Samples that were removed included duplicate entries, ages younger than 0 Ma, ages from sedimentary units that may not represent recent tectonic activity, and ages from volcanic units that may represent tectonic activity distal from the sample location (i.e., tuffs) (Fig. 4). After the data were culled, we generated a grid across the continental United States with 0.25° grid-point spacing. At each grid point, we identified the minimum age of the samples located within 1° and assigned that value to the point. Areas with no samples within 1° were not assigned values. In order to extrapolate this model across the continent to regions that were geochemically undersampled, we fit a surface to these grid points via a locally weighted linear regression using the nearest 5% of data points from each undersampled 0.25° grid point. Grid points where the fit produced ages less than 0 Ma, which were primarily at the edges of the continent, were removed (Fig. 4). The accuracy of this continent-wide age model is limited by the spatial resolution of the EarthChem database. However, this model largely reflects the tectonic activity described above, and the interpolation to uniform spacing allows a straightforward comparison to seismic-velocity models, which are smoothed via model regularization.

Temperature Calculations

We use mineral-physics–based models to quantify the anharmonic and anelastic effects of temperature and pressure on shear velocities in the continental mantle (Goes et al., 2000; Jackson and Faul, 2010; McCarthy et al., 2011). Because USArray is a continental array, we did not utilize other successful approaches to convert seismic velocities to temperatures that were calibrated with thermal models of oceanic lithosphere (Priestley and Mckenzie, 2006, 2013). In order to convert shear velocities to temperature, estimates of activation energy, activation volume, grain size, and pressure are required and are listed in the Supplemental Material1 for the different models.

Given the many factors that influence mantle seismic velocities, isolating any one variable is inherently an underdetermined problem. As discussed above, temperature is the dominant control on mantle seismic velocities, and thus, we take a first-order approach by associating all variations in velocity with temperature variations. However, these temperature estimates may represent maximum temperatures for mantle regions that are more likely to be affected by water or partial melt (i.e., high-temperature regions, such as the margins of the Colorado Plateau, Cascadia, Yellowstone, Basin and Range, etc.).

The three mineral-physics–based models utilized for our calculations were; (1) Jackson and Faul (2010); (2) Goes et al. (2000) utilizing updated inputs from Cammarano et al. (2003); and (3) McCarthy et al. (2011) using the unrelaxed shear modulus for olivine from Isaak (1992). Input values for the three models are given in the Supplemental Material (footnote 1), and readers are referred to the original works for specific details on the methodologies. Pressure was estimated using the mean crustal thickness from the three models. A crustal density of 2750 kg/m3 and an upper-mantle density of 3300 kg/m3 were utilized for this calculation. To investigate the effects of grain size on the seismic velocity, we utilize three different grain sizes (10 μm, 100 μm, and 1 mm) in our temperatures calculated for the Jackson and Faul (2010) model.

Seismic-Velocity Comparisons

The four shear-velocity models utilized in this work all exhibit similar broad trends found in previous continent-wide work (e.g., Bedle and van der Lee, 2009). These models include reduced velocities in the west and high velocities in the central and eastern United States (Figs. 1 and 2). The lowest velocities within the models are observed in the vicinity of Yellowstone, the margins of the Great Basin, and beneath the margins of the Colorado Plateau. Additional low velocities are observed along the Cascadia margin in the Shen and Ritzwoller (2016) model (Fig. 2). All of these regions are tectonically and/or volcanically active. The highest velocities are observed beneath the North American craton, where the oldest lithosphere is present. Intermediate velocities are along the eastern margin of the continent (Figs. 1 and 2).

The three shear-velocity models that used USArray data broadly agree with each other, especially in the western interior and eastern margins of the continental United States, though variations exist (Figs. 2 and 3). The three models also broadly agree with the pre-Earthscope model (NA04) (van der Lee and Frederiksen, 2005) (Figs. 1 and 2) but allow for much higher resolution results. The pre-EarthScope model misses smaller features, such as low velocities beneath Yellowstone and the Snake River Plain and around the margins of the Colorado Plateau.

There are increased variations in the three USArray velocity models within the central United States; these increased variations may be partially due to uncertainty in estimated crustal thicknesses between the models, and within Cascadia, where serpentinization and/or partial melting of the upper mantle may affect results. In the Cascadia region, serpentinization likely affects crustal thickness estimates and estimates of upper-mantle velocity, an issue that has been noted in previous seismic work (Bostock et al., 2002; Brocher et al., 2003; Gilbert, 2012). All three of the models represent significant improvements in resolution when compared to previous continent-wide velocity models (e.g., van der Lee, 2002; Bedle and van der Lee, 2009) due to the large aperture and high resolution of the EarthScope transportable array; these improvements in resolution were not available when these previous models were calculated.

Velocity versus Age

Using the seismic-velocity models and the tectonic age model described above, we produce plots of sub-Moho mantle seismic velocities versus timing since the last thermotectonic event for the continental United States. We first plot data from the EarthChem database versus the three S-wave velocity models and the upper-mantle P-wave velocity model of Buehler and Shearer (2016) (Fig. 5). This does not include any interpolation of the age data, and velocities are taken from the nearest grid point to the sample site. The plots show the same broad trend as our interpolated model. We then plot tectonic ages from the interpolated age model versus EarthScope-derived seismic velocities (Fig. 6). Finally, we used the average shear-velocity model (Fig. 3) calculated from the three shear-velocity data sets, and we plotted these velocities versus ages (Fig. 7).

These comparisons show an expected positive correlation between thermotectonic ages and upper-mantle seismic velocities (Figs. 57). Seismic velocities exhibit a sharp increase in the first 500 million years following tectonism; this increase becomes more gradual with time, which is consistent with previous models of lithospheric cooling (Fischer, 2002). Scatter in the data likely reflects errors in the velocity or age data or other factors, such as varied cooling rates, grain-size variations, lithospheric and crustal thicknesses variations, the presence of partial melts and/or fluids, compositional variations, and differences in orogenic heating impacting our measurements.

Temperature Estimates

Given that our temperature estimates are controlled by seismic velocities, similar patterns are observed in both our temperature and shear-velocity maps, though in our temperature calculations, temperature is inversely correlated with age (Figs. 810). Regardless of the mineral-physics models used, all temperature calculations show similar trends with the highest temperatures observed in the western United States, the coldest in the central United States, and intermediate at the eastern edge of the continent (Figs. 8 and 9). In these models, upper-mantle temperatures in the eastern and central United States range from ∼600 to 1000 °C, while western U.S. temperatures range from ∼700 to 1450 °C. The trend in temperatures within the western United States is similar to temperatures calculated by Schutt et al. (2018) using Pn waves, though our absolute temperatures are slightly higher. Temperatures calculated using the three models are largely similar but exhibit slight variations, especially at higher temperatures. These variations are primarily due to differences in how the anelastic components of the temperature–shear-velocity relationship are calculated, as well as the choices in inputs such as the anharmonic velocity model, activation energy, activation volume, etc. These inputs are discussed in greater detail in the Supplemental Material (footnote 1). Overall, our temperature calculations are in broad agreement with previous estimates of attenuation from the continental United States, with higher attenuation in the western United States, lower in the central United States, and intermediate in the U.S. interior (Dalton et al., 2008, 2017; Lin et al., 2012; Bao et al., 2016; Bowden et al., 2017). Direct comparisons between previous attenuation results and our measurements are challenging given that surface-wave attenuation results are presented in terms of surface-wave periods, which integrate measurements over a broad depth range.

Grain Size

Using the Jackson and Faul (2010) model, we demonstrate the effects of grain size on temperature estimates for the continent (Fig. 8). In the models, differences in calculated temperatures are most pronounced in the western United States where higher temperatures and greater anelastic effects are observed. The model with the smallest grain size (10 μm) exhibits the highest temperatures, while the largest grain size (1 mm) exhibits the lowest. Grain size in the mantle has been estimated using seismic data, xenoliths, and mineral-physics data but is still not well constrained and may vary considerably both vertically and laterally (Hirth and Kohlstedt, 1996; Armienti and Tarquini, 2002; Hirth and Kohlstedt, 2003; Behn et al., 2009; Priestley and Mckenzie, 2013). We assume a grain size of 1 mm for our plots of age versus temperature. Better constraints on grain size will be important for future temperature calculations and for determining both the mode of deformation and rheology of the mantle.

Temperature versus Age

We plot our estimates of the upper-mantle temperatures versus our age model in order to better understand lithospheric cooling following tectonism (Fig. 10). These plots show rapid cooling in the first ∼300–500 million years following tectonism; cooling slows as the lithosphere ages. This result is consistent with simple lithospheric cooling models.

To better understand this postorogenic cooling, we fit a curve based on lithospheric cooling models to our temperature and age data. In order to only fit the most robust data in our plots of temperature versus age, we utilize our average shear-velocity model (Fig. 3) and only fit data points where standard deviation of the velocities from the three tomographic models was less than 0.1 km/s and the crustal thickness standard deviation was less than 5 km. We then calculate a least-squares fit to these temperature and age data using a simple equation for cooling of continental lithosphere; this equation does not account for preexisting lithosphere thickness or heat-producing elements in the lithosphere (Turcotte and Schubert, 2014):
where t is time; T0 is the surface temperature; T1 is the temperature throughout the lithosphere at t = 0 and the temperature at the base of the lithosphere when t > 0 ; y is depth; yl is the maximum lithospheric thickness; and κ is the thermal diffusivity of the lithosphere. This model uses a maximum lithospheric thickness and does not account for radioactive heating within the lithosphere. We fix y at 46 km, which is the mean depth of the Moho in our models, and we use a least-squares fit to solve for the other variables, which are listed in Table 2. Though we only fit the most robust data, all temperature and age data are plotted in Figure 10. T0 values are expected to be close to 15 °C, the average temperature at Earth’s surface. The calculated values for T0 are higher than expected, which is possibly due to the fact the radioactive heating has not been taken into account. Estimated κ values are near the range expected from laboratory estimates, which commonly are between 4 and 22 × 10−7 m2/s−1 (e.g., Gibert et al., 2003).

Model Uncertainty

Further uncertainty in our first-order calculations stem from uncertainty in the seismic models and in the mineral-physics models used to calculate temperature from shear velocity. Shen and Ritzwoller (2016) conducted a detailed error analysis and calculated an average non-systematic model error of ∼0.4% in S velocity that increases near the Moho and does not account for trade-offs or sub-resolution structure. In the mineral-physics models utilized, a ∼0.001 km/s change in velocity would correspond to a ∼1 °C change in estimated temperature. Using these values, we would expect an average error in temperature of up to 20 °C based solely on seismic model uncertainty. Additionally, uncertainties in crustal thickness impact pressure calculations, and differences in regularization techniques for the different seismic models affect temperature calculations. We attempt to minimize these differences by focusing our analyses on areas where the models are in good agreement.

The mineral-physics measurements used to relate temperature to seismic velocities are typically scaled to seismic frequencies and upper-mantle conditions. This is because these measurements often were made at higher frequencies, lower pressures, and on materials that exhibit smaller grain sizes than expected within Earth’s mantle. To account for these variations, scaling factors are used to predict the effects of temperature on observed velocities within the mantle. Significant improvements have been made recently on these scaling factors (e.g., Jackson and Faul, 2010; McCarthy et al., 2011); nevertheless, they still introduce uncertainty into the model. Furthermore, uncertainties in grain size can introduce error into these calculations. Finally, fluids and/or melts, if present in sufficient quantity, can significantly impact seismic velocities. This is likely a concern in regions where recent volcanism has occurred, such as along the margins of the Colorado Plateau, the Yellowstone Region, and Cascadia. Due to the numerous sources of potential uncertainty and the nonlinear relationship between each input and temperature, we do not formally quantify uncertainty within our temperature models. However, Figure 10D gives a rough estimate of how the use of different mineral-physics models affect temperature estimates, and R-squared values for lithosphere cooling curves are given in Table 2.

Effects of Orogeny

Within orogenic belts, lithospheric heating is often associated with advection and convection, both of which involve flow of hot asthenosphere upward toward the base of the overriding lithosphere. This has been identified as corner flow in subduction zones (e.g., Stern, 2002), as small-scale convection in continental interiors (e.g., Zandt and Carrigan, 1993; van Zandt et al., 2004; Hales et al., 2005; Schmandt and Humphreys, 2010; Wijk et al., 2010), and as vertical plumes in a variety of tectonic settings (Campbell and Griffiths, 1990; Geist and Richards, 1993; Miller and Nicholson, 2013). Within subduction zones, this flow is driven by the motion of converging plates (Kincaid and Sacks, 1997; Stern, 2002). In continental interiors, this flow is associated with sharp gradients in lithospheric thickness and/or lithospheric downwellings (Elkins-Tanton, 2007; King, 2007; van Wijk et al., 2010). Plumes are commonly associated with deep-seated heat sources (Montelli et al., 2004; French and Romanowicz, 2015) or slab-rollback–induced flow (James et al., 2011), which lead to buoyancy-driven upwelling of mantle material due to thermal expansion. Because asthenospheric flow appears to be an important mechanism for lithospheric heating, orogenic events that produce favorable conditions for vertical flow may be more likely to produce significant heat signatures.

Across the United States, the differences in observed upper-mantle temperatures correlate well with time since the most recent tectonism, but it is likely that the style of orogeny played a role as well. Orogenic events that result in increased vertical mantle flow and lithospheric shortening or extension will likely lead to increased warming and therefore will impact on the thermal state of the lithosphere. Deformation from lateral motion at transform boundaries typically results in little vertical movement within the lithosphere, emplaces few heat-producing elements, and only results in localized frictional heating (Zoback et al., 2011). Transform boundaries should produce small thermal signatures, as evidenced by the lack of volcanism associated with strike-slip faulting. Within subduction zones, elevated temperatures are observed within the arc and backarc due to mantle flow (Syracuse et al., 2010). Within continental-collision zones, the thermal conditions are affected by crustal thickening, which can increase the quantity of heat-producing elements in the crust, by frictional heating and by variations in thermal diffusivity that can affect heat flow (Molnar et al., 1983; McNamara et al., 1997; Whittington et al., 2009). Continental collision is commonly preceded by subduction, which can lead to heating due to asthenospheric flow around the slab remnants after continental collision occurs (Kind et al., 2002).

In global seismic models, the western United States exhibits the lowest velocities for continental lithosphere on Earth even when compared to other orogenic zones such as the Himalaya (Debayle et al., 2016; Moulik and Ekstrom, 2016). This suggests that the thermotectonic events within the United States may have resulted in especially high or widespread elevated temperature within the lithosphere. These high temperatures are likely due to unique conditions within the region including warming of a hydrated continental upper mantle (Dixon et al., 2004) associated with sinking of the Farallon slab and upwelling due to extension, mantle plumes, return flow from downwellings, and possible small-scale asthenospheric convection (Farmer et al., 2008; Schutt and Dueker, 2008; West et al., 2009; van Wijk et al., 2010; Hansen et al., 2013; Porter et al., 2014; Porter et al., 2017). Even though the short-term thermal effects of these processes may be substantial, the long-term (>300–500 Ma) thermal conditions should appear similar to other orogenies, as shown within our data, because cooling is characterized by progressively smaller (exponential) decay in temperature (Fig. 10).

Heating and Cooling Rates

Because lithospheric warming is driven by the movement of warm material, through convection or advection and by concentrating heat-producing elements within the crust, time scales for lithospheric warming are likely much shorter than for lithospheric cooling, which primarily occurs through conduction. The Colorado Plateau provides an example of rapid heating. Prior to the Laramide orogeny, the Colorado Plateau was located at sea level and was relatively unaffected by nearby orogeny. If we use modern temperatures calculated for the relatively undeformed interior of the Colorado Plateau as a maximum temperature at that time, the plateau uppermost mantle was likely 900 °C or less during the latest Cretaceous. Since that time, the region has been uplifted 2 km, and temperatures within the mantle along the margins are possibly in excess of 1200 °C. This minimum heating rate of ∼5 °C per million years is higher than expected for conduction alone using the values from Equation 1 and, as such, requires vertical flow within the mantle or other mechanisms to drive heating.

Our estimates of lithospheric cooling show that most cooling happens within the first ∼300–500 years following orogeny, which is consistent with temperature predictions from our cooling model. Thermal models of cooling at Moho depths are also consistent with these observations (Fischer, 2002). In Figure 10, the mean standard deviation of the estimated temperatures at given age is about ∼100 °C in the Phanerozoic and ∼50 °C in the Precambrian. These standard deviations reflect uncertainties in parameters used as well as effects that we did not take into account. Our temperature curves show a decrease of approximately four times the Phanerozoic standard deviation from 541 to 0 Ma and a decrease of one standard deviation of the Precambrian region temperatures from 3500 to 541 Ma. Therefore our cooling rates for Precambrian regions are less robust than the Phanerozoic regions and may reflect data uncertainty. Future work will focus on incorporating heat-flow data (e.g., Goes and van der Lee, 2002) and crustal heat production (e.g., Mareschal, 2004) into these estimates to better understand these processes.

Temperature calculations from this work show a range of thermal signatures across the continental United States; these signatures correlate well with the time expired since the most recent thermotectonic events. Because these temperatures were calculated at depth near the Moho, these values can be used as an approximation of Moho temperature and allow us to examine the long-term effects of cooling on the lower crust. These thermal estimates are consistent with work showing densification of the crust over time due to metamorphism associated with cooling (Fischer, 2002). This cooling process likely has a long-term effect on surface topography and the support for mountain ranges, as well as, on Moho topography, as illustrated for South America by Lloyd et al. (2010). Future analyses of the thermal state of the lithosphere will help differentiate between the effects of cooling and fundamental differences in orogeny (e.g., Schmandt and Lin, 2015) in controlling the long-term evolution of topography.

This work summarizes some uppermost-mantle seismic imaging results from the EarthScope program. Broad trends in seismic-velocity variations imaged in these studies had been found in previous surface-wave work (e.g., Black and Braile, 1982; Grand, 1994; van der Lee and Nolet, 1997; van der Lee, 2002; Bedle and van der Lee, 2009), but the resolution is dramatically improved owing to the station density of the EarthScope Transportable Array. This improved resolution allows for the identification of fine features in the upper mantle, which has led to the refinement of models of continental lithosphere evolution. By using data products from the EarthScope program, we are able to examine the modern thermal state of the North America lithosphere and show that temperatures in the upper mantle vary dramatically across the continent with estimates ranging from >1300 °C in parts of the tectonically active western United States to <500 °C beneath the North American craton. By combining these seismic results with additional publicly available geochemical data, we are able to approximate cooling rates for continental lithosphere following thermotectonism and show that empirically measured cooling rates are consistent with models of lithospheric cooling.

The inspiration for this work came out of the EarthScope synthesis workshop “Synthesizing EarthScope Results to Develop a New Community Model for the 4-D Evolution of North America,” held at James Madison University in November 2016.

1Supplemental Material. Temperature-velocity models. Please visit or access the full-text article on to view the Supplemental Material.
Science Editors: Raymond M. Russo, David E. Fastovsky
Guest Associate Editors: John Hole, Michael L. Williams
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.