{'fraction_non_alphanumeric': 0.03005444991321302, 'fraction_numerical': 0.025417885250968924, 'mean_word_length': 4.920326576576577, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Continental Tectonics Inferred From High‐Resolution Imaging of the Mantle Beneath the United States, Through the Combination of USArray Data Types

A comprehensive North American upper mantle seismic‐tomographic model, NA13, was inferred from a combination of several decades of seismic data from early North American seismic networks and USArray data. This data‐driven modeling inverted a combination of regional waveform fits, teleseismic S delay times, and point constraints on Moho depths. The joint inversion model combines the contrasting and overlapping resolving power of the different data sets, and demonstrates enhanced resolution over regional models created with a single geophysical data set. Resolution in the upper mantle is achieved on the scale of ∼100 km, although travel time delay studies suggest that anomalies in parts of the western United States are smaller. In NA13, velocities beneath the Yellowstone plume are low enough to require the presence of partial melt. NA13 also models a variety of smaller scale velocity variations beneath the Central and Eastern United States that reveal remnant velocity anomalies related to lithospheric variations in the upper mantle between Proterozoic units.

focus of many local tomographic studies as USArray moved eastward (e.g., Bollmann et al., 2019;Levin et al., 2018;Long et al., 2019), and improved understanding of these anomalies may provide additional insight into current central and eastern earthquakes. Most large-scale tectonic episodes leave remnant signatures in the lithosphere and asthenosphere by altering the thermal and mineralogical structure. Improved tomographic modeling that employs multiple types of data sets, such as NA13, allow for the detection and resolution of these anomalies to aid in better constraining paleotectonics.
Further improvement in resolving North American mantle structure results from better accounting for crustal structure and the simultaneous inversion of multiple data types Obrebski et al., 2011;Tian et al., 2011). The model presented here, NA13, combines body wave travel times, constraints from regional S and Rayleigh waveform trains, and independent Moho depth constraints to jointly invert for a three-dimensional S-velocity (V s ) model for the North American mantle. The joint inversion leads to better velocity structure resolution than can be achieved by any of the individual data sets, as previously demonstrated by Schmid et al. (2008) and Chang et al. (2010). In addition, model NA13 demonstrates improved lateral resolution of upper mantle features on the order of 1° × 1°, compared to other joint-inversion models which demonstrate limits of 5° × 5° (e.g., Golos et al., 2018;).

Model Parameterization
The seismic velocity of the three-dimensional model is parameterized by 25 spherical shells of grid points at 25 depths (0, 2, 20, 35, 60, 90, 120, 155, 190, 230, 270, 320, 370, 420, 480, 540, 600, 670, 740, 820, 900, 1,000, 1,100, 1,300, and 1,500 km). An additional spherical shell of grid points represents the laterally varying Moho depth. The vertical spacing of the spherical shells within the mantle increases gradually with depth to accommodate decreasing resolution. A triangular tessellation method discretizes each spherical shell (Baumgardner & Frederickson, 1985;Van der Lee & Nolet, 1997b;Wang & Dahlen, 1995). The seismic velocity is defined at any point within the mantle or crust by linear interpolation of values at the six surrounding grid points, which are nodes of triangles in each of two spherical shells, right above and below the point. Each spherical shell consists of 176 triangles and 66,266 vertices, yielding 1,722,916 grid points. The triangle-side length ranges from 40 to 47 km with an average of 45 km on the surface, which is less than the TA station spacing of 75 km. The triangular grid is centered 45°N and 95°W and extends 55° in all directions from the center as shown in Figure 2a. Moho depth is defined using linear interpolation of values at the three nodes of its bounding triangle in the Moho shell. Figure 2 additionally displays the input parameters and data for NA13, with station coverage and events displayed in Figures 2b-2d. Coverage extends past the North American continent to accommodate the extent of the body and surface wave paths. In this study, 386,762 grid points are hit by teleseismic S waves, while 876,864 grid points are hit by the regional S and Rayleigh waveform fits. In total, the data samples 960,947 grid points when combining the body and surface waves.

Joint Inversion
Teleseismic S wave arrival times and regional S, multiple S, and Rayleigh waveform fits are jointly inverted for S velocity perturbations of a reference model. This joint inversion combines the resolution and sensitivities of body and surface waves to create a more accurate and detailed velocity anomaly model. Regional S and multi-mode Rayleigh waveform fits were measured by the method of partitioned waveform inversion (PWI; Nolet, 1990;Van der Lee & Nolet, 1997a, 1997b. Independent Moho depth measurements such as from receiver functions, active source surveys and gravity measurements are also incorporated, along with appropriate relative uncertainties, to improve resolution of the velocity structure. The joint inversion method is similar to previous work by Feng et al. (2007), Schmid et al. (2008), andChang et al. (2010) and described in detail for this inversion in Lou (2013). The joint inversion of differing data sets provides enhanced resolution of the velocity structure because each data set can improve the resolving power of other datasets (e.g., surface waves best constrain the uppermost mantle, and anomalies not accounted for in the upper mantle can be positioned in the lower mantle by the teleseismic waves). In addition, the joint inversion of these various datasets allows for the creation of a velocity model that has superior placement and detection of smaller scale velocity anomalies where data coverage allows.
A typical seismic tomography problem is to solve a large, sparse, and mixed determined system of linear equations in the form of Gm = d where m is the unknown model vector, d is the known data vector, and G is the modeled sensitivity kernel matrix. The model vector, m, includes the three-dimensional S velocity 4 of 21 perturbations to a one dimensional reference model, two dimensional Moho depth perturbations to a reference value and corrections for earthquake origin times and hypocenter locations. Residuals between the observed data and data predicted by the reference model are arranged in the data vector, d, with the contributions from the multiple data types. Partial derivatives of these data, with respect to the model unknowns, formulate the G matrix. The set of equations (Gm = d) is solved for the best fitting model (m) using the iterative conjugate-gradient method LSQR (Paige & Saunders, 1982a, 1982b with regularizations. Because wave paths are irregularly distributed in the model space, the inverse problem is partly overdetermined and partly underdetermined. Therefore, regularization is added to the inversion in terms of flattening. The flattening matrix limits the spatial gradient of the velocity model, which also limits the anomaly Figure 2. Map of (a) triangular grid (red points) of a spherical shell, which represents the model space of NA13. Note that the actual triangle-side length of the grid for NA13 is 45 km, and is too small to be properly visualized at this global scale. So, a grid with the same area coverage, but twice the horizontal spacing is plotted for illustration. (b) Seismic stations for the body wave travel time data set, (c) teleseismic earthquakes and stations for the body wave travel time data set, (d) map of seismic stations (red triangles) and regional earthquakes (green circles) with regional waveform constraints (from Bedle, 2008;Bedle & van der Lee, 2009). size to some extent. Because of this, we only use the flattening matrix in the tomographic joint inversion, rather than also adding damping equations. Details of investigation into the model regularization and resultant travel time misfit are presented in Lou (2013), but a preferred flattening multiplier of 0.5 was chosen as it balanced the trade-off between the waveform fit data's variance reduction and the model's smoothness. Investigations were also performed to determine the ideal weighting between the different data types in the inversion. Preferred weighting based on parameter testing resulted in a travel-time data weight of 1, a waveform fit weight of 2, and a weight of 8 for the point constraints on Moho depth (additional details see Lou [2013]). The RMS residuals for the combined data normalized over 500 iterations converged to about 0.25 after a few dozen of the iterations.

Input Data
As previously mentioned, different types and sets of seismic data are combined to jointly invert for seismic velocity structure beneath North America. The Crust 2.0 model (Bassin et al., 2000) was used to initially define the seven-layer one-dimensional crustal model that consisted of ice, water, soft sediments, hard sediments, upper crust, middle crust and lower crust on a 2° × 2° global grid. The models of crustal thickness and V p /V s from Lowry and Pérez-Gussinyé (2011) and Lowry et al. (2013), both whom jointly inverted seismic receiver functions, gravity and surface heatflow measurements, were used to correct Moho depth delay times at each station. For stations not in the grid of Lowry's model, we used the Moho depth from NA04 (Van der Lee & Frederiksen, 2005). The body wave travel time data are comprised of teleseismic P and S wave absolute delay times from previous IRIS PASSCAL arrays, EarthScope's USArray, and various permanent networks as displayed in Figures 2b and 2c. The delay times were measured using AIMBAT . Regional S and Rayleigh wave trains have been fitted using the method of Partitioned Waveform Inversion (PWI) as first introduced by Nolet (1990). PWI methodology and algorithms were further updated to include Moho depth as a variable by Das and Nolet (1995). Van der Lee and Nolet (1997b) performed the first waveform inversion for both continental scale variations in Moho depth and upper-mantle S-velocity structure beneath North America. Model NA13 draws on a total of 9,331 waveforms from 252 regional earthquakes that were used for three-dimensional model NA07 (Bedle & van der Lee, 2009), an update of NA04 (Van der Lee & Frederiksen, 2005). An additional 1,255 waveforms from five earthquakes recorded by USArray's Transportable Array in the western United States were used to refine knowledge of the western United States (Bedle, 2008) and are included in this inversion. Lou (2013) provides additional detail on stations, events, and delay time measurements. The stations, regional, and teleseismic earthquakes are plotted in Figure 2.

Results
A detailed study to determine the optimal regularization parameters for the full NA13 data set was performed by Lou (2013). In that study, weighting between the data types was determined through trial inversions to minimize variance in the data. NA13 extends deeper than previous PWI North American models NA04 (Van der Lee & Frederiksen, 2005) and NA07 (Bedle & van der Lee, 2009) due to the addition of teleseismic arrival-time data in the joint inversion.

S Velocity Model NA13
Model NA13 (Figures 3 and 4) is the joint S velocity model inverted simultaneously from teleseismic S wave delay times, regional S and Rayleigh waveform fits, and local Moho depth constraints. The velocities in this model are plotted as a percent deviation from the IAASPI91 velocity model (Kennett & Engdahl, 1991). As a result, relatively small scale (∼100 km) velocity variations in the upper mantle can be more clearly imaged than when using a single data set. In the western continent, NA13 displays narrow bands of relatively higher velocity anomalies running parallel to the coast along the Cascades subduction zone between depths of 100 and 300 km. This velocity anomaly is typically interpreted as corresponding to the subducting Juan de Fuca and Gorda slabs. The Snake River Plain and the Basin and Range province are characterized by low velocity anomalies between 100 and 300 km depths. The strongest low velocity anomaly at these depths is imaged beneath the Yellowstone hotspot in the northwestern corner of Wyoming. Both the Columbia River Plateau, in Washington and Oregon, and the Colorado Plateau exhibit relatively high uppermost mantle velocity anomalies at 100 km depth when compared to the rest of the western United States. Surrounded by lower velocities, the high-velocity Colorado Plateau's southern segment is of a higher velocity than the northern segment at 150 and 200 km ( Figure 3). A similar difference has been observed in regional models of the Colorado Plateau, including that of Liu et al. (2011), these variations lay just to the east of the proposed region of crustal delamination as proposed by Levander et al. (2011). In the Southern Rocky Velocities are plotted throughout as a percent deviation from the IASP91 1D velocity model (Kennett & Engdahl, 1991).
Mountains province, low velocities at 150 and 200 km are modeled in a north-south trend beneath the Rio Grande Rift on the eastern side of the Colorado Plateau ( Figure 3).
East of the Rocky Mountains, the tectonically stable Great Plains, Central Lowlands, and adjacent smaller provinces, display higher S velocity anomalies in the upper mantle, which reflect the relatively cool and rigid cratonic lithosphere of the North American Craton (Figures 3 and 4). These high-velocity anomalies extend from the tectonically active western United States to the Coastal Plains. Most of this lithosphere is Proterozoic in age. Velocities are higher in the northern Central Lowlands and Great Plains, including the Wyoming Craton, where the lithosphere is Archean in age. Velocities appear slightly less high beneath the Yavapai and Mazatzal provinces on the southwestern and southern edges of the North American Craton. Further to the south and east, the Paleozoic Ouachita and Appalachian orogenies display lower velocity anomalies, likely reflecting their more recent tectonic stabilization.
NA13's shallow, lithospheric portion of the United States agrees at very large scales with lithospheric models created via ambient noise tomography (Bensen et al., 2009), but additionally images smaller scale features, such as high velocities related to subduction in the Pacific Northwest, and higher velocities beneath northern Idaho. The Bensen et al. (2009) model better images intra-crustal structures than NA13. NA13 models lower velocities in the uppermost mantle beneath New England, which are not observed as prominently in the Vs ambient noise models of Bensen et al. (2009) potentially because they occur relatively deep within the cool lithosphere. Another ambient noise model by Zhao et al. (2017) also reveals similar velocity structures to the large scale velocity features of NA13 such as the slower velocities beneath the Basin and Range and Rio Grande Rift, although the velocity maps created therein are smoother and are limited to larger-scale features than NA13 due to the smoothing and damping requirements in their ambient noise method.
In the transition zone and lower mantle, the magnitudes of the velocity anomalies are smaller. Relatively high velocity structures are present beneath southeastern Washington, northeastern Oregon, eastern Nevada, and western Utah in the transition zone ( Figure 3). In the eastern United States, high velocity anomalies are also observed in a broad zone throughout the transition zone and down to 900 km depth. These structures are not as anomalous as those beneath the western continent. The high velocity feature dipping to the east can be observed on the left side of cross section A-A′ of Figure 4. At the large scale, this feature is similar to that modeled by Sigloch (2011) using teleseismic tomography to detect V p anomalies. Differences between the two models also exist in the uppermost mantle, such as lower velocities at the northern edge of the Coastal Plain in NA13, and significantly lower V s beneath the Columbia River Plateau and the eastern Basin and Range province. Higher velocities extend to the south beneath Texas in NA13. In the immediate region of the transition zone and 410 km discontinuity, studies have mapped the topography of this discontinuity (Cao & Levander, 2010;Tauzin et al., 2013). At 400 km (Figure 4), NA13 reveals a relatively high velocity beneath northern Nevada around 400 km, which is in congruence with a shallower 410 discontinuity as imaged by Cao and Levander (2010), suggesting that these high velocities could be due to the same cool thermal regime responsible for shallowing the discontinuity in this region. Tauzin et al. (2013) also present a western US 410 km topography model, noting a similar shallow region in northern Nevada, in line with Cao and Levander (2010) and the high velocities observed in NA13. A similar shallower topography in Tauzin et al.'s (2013) model also is imaged at the northernmost limits of their model, beneath northern Idaho, which are also in the same region as relatively high velocities in NA13. This alignment between 410 km discontinuity topography maps and the NA13 velocity model confirm that the dominant origin of seismic velocity variations is likely temperature.

Resolution Tests
To examine how well and to what scales the imaged velocity anomalies are resolved by the tomographic inversions, resolution tests are necessary. Because the applied regularization would render a biased model covariance matrix, the data's resolving power is tested for structures of general and particular interest. A typical approach to quantify resolution is the checkerboard test. A series of typical checkerboard resolution tests are conducted on hypothetical models using the same inversion parameters (including regularization and data type weighting) as for model NA13. As displayed in the left columns of Figures 5 and 6, two hypothetical checkerboard models are created by placing a ±200 m/s Vs perturbation into repeating rectangular boxes of sizes of 300 and 100 km. For every input model, the test outputs are plotted in three horizontal slices and four vertical cross sections in the right column of the corresponding figure. Two of the three horizontal slices are at depths to display the input velocity anomalies and for testing the data's ability to recover the presence of these velocity anomalies. The other horizontal slice is for a depth without any input velocity anomaly to test the data's ability to recover the absence of a velocity anomaly.
The first checkerboard test ( Figure 5) is for velocity anomaly boxes of size 300 × 300 × 300 km. These velocity heterogeneities are well recovered throughout the upper mantle, and through most of the lower mantle, particularly beneath the continental US, where most of the seismic stations are located. The horizontal resolution becomes less ideal to the east of the −90° meridian due to poor station and data coverage. Velocities are smeared across the US borders with Canada and Mexico, and across the east and west coasts. This is also due to relatively abrupt termination of USArray in these regions at the time. The resolution in the upper mantle is better than in the lower mantle because the teleseismic body waves have more crossing paths in the upper mantle than in the lower mantle. Some smearing of the velocity anomalies extends into the transition zone, where no input velocities exist in the resolution model. This smearing is minor and less than 0.5% of the total S velocity.
The second resolution test has boxes of the size of 100 × 100 × 100 km (Figure 6). At depths of 300 and 500 km, the models show strong smearing along the southwest-northeast direction. This is because the majority of the teleseismic earthquakes are from the back azimuth of northwest and southeast ( Figure 2). Despite the higher degree of smearing, this checkerboard test does reveal that velocity anomalies of the 100 km size should be detectable, and resolved, particularly in the upper mantle and beneath the western half of the United States. Additional resolution tests for NA13 are presented in Lou (2013), and demonstrate that the resolution for the lower mantle beneath the western and central US is ∼200 km. Due to data coverage for NA13 favoring the western and central USA, the eastern US has a resolution of 300-500 km in the upper mantle, and minimal resolution in the lower mantle of this region.

Predicted Delays Through NA13
Tomographic models can be assessed by predicting delay times of seismic data though the velocity model and comparing them to the observed delay times. Figure 7 shows the station-average observed delays for NA13. NA13 has the ability to recreate observed S delay time patterns within the western and central United States. In the eastern US, NA13 also performs well in the regions with good data coverage. Interestingly, despite high station and data coverage, NA13 has some difficulty in recovering the correct delay times in California, beneath northwestern Wyoming, and in the regions encircling the Colorado Plateau. While discussed in more detail in Lou and van der Lee (2014), this could be attributed to the idea that while NA13 has high resolution of ∼100 km in this region, the velocity anomalies are actually of an even smaller range, causing the misfit of delay times.

Discussion
S velocity anomalies are examined in more detail for the jointly inverted model, NA13, with a focus on the tectonic relevance. Figure 4 displays two horizontal slices at 200, 400 km depths, and two cross-sectional profiles across the North American continent. Figures 8 and 9 display additional depth slices and cross-sectional profiles through NA13 to highlight mantle velocity anomalies of particular interest.

Pacific Northwest Subduction
In the northwestern US, model NA13 images a north-to-south elongated strong high velocity band along the Cascades at 200 km depth. Profile AA′ through the southern tip of this anomaly shows to be an east-dipping velocity anomaly likely representing the subducting Juan de Fuca plate. This high-velocity anomaly extends to the transition zone where it appears to flatten beneath the Basin and Range. Profile BB′ cuts through these anomalies at a different azimuth, and suggests a possible continuity with high-velocity anomalies at about 800 km depth, in agreement with other regional tomographic models (e.g., Schmandt & Humphreys, 2010). This higher velocity features does not tend to be as clearly imaged in global scale S-velocity models, including Simmons et al.  are also imaged in the lower mantle beneath provinces east of the Great Plains (Figure 3). These high velocity anomalies are identified as subducted Farallon lithosphere as originally noted by Grand (1994). The high velocities are caused by thermal differences between the lithospheric slab and the surrounding lower mantle material. Our images of the subducted Farallon plate are largely consistent with previous tomographic studies (Burdick et al., 2008(Burdick et al., , 2012Grand, 2002;James et al., 2011;Obrebski et al., 2010Obrebski et al., , 2011Schmandt & Humphreys, 2010;Shearer & Buehler, 2019;Sigloch et al., 2008;Tian et al., 2009Tian et al., , 2011Van der Lee & Nolet, 1997a). Figure 7 additionally focuses in on slices and transects beneath the western portion of NA13 in the region of the subducting slabs. Profile CC′ shows the high-velocity slab signature beneath Washington, Oregon, and northern California. The southern edge of the high-velocity slab contrasts sharply with low upper-mantle velocities in the slab window to its south. The 4% velocity transition occurs over less than 100 km. Profiles DD′, EE′, and FF′ are perpendicular to CC′ and show more segmentation of the high-velocity slab anomalies, including complex fragments in the transition zone. Profile EE′ is located where subducting lithosphere is the youngest and where the Gorda Ridge nearly meets the trench. This profile barely shows a high-velocity slab-related anomaly in the mantle above the transition zone. Previous tomographic models have also imaged a gap in the high-velocity slab-related anomaly (Obrebski et al., 2011;. This gap has been proposed to be due to slab break-off (Obrebski et al., 2011;Roth et al., 2008), but more recent studies of receiver functions (Tauzin et al., 2016) do image reflections in the region, suggesting the subducted slab is present, but not easily detected in previous tomographic models, possibly due its young age upon subduction. Interestingly, the region of this weaker slab corresponds to a region of faster continental crust, and less seismic tremor based on a study by Delph et al. (2018). In this region, Delph et al. (2018) suggest that their data can be explained by a lower fluid content in the slab, or due to a lower fluid permeability that is not allowing the mantle wedge to be serpentinized. The lower seismic velocities observed in NA13 are slightly north of where the Gorda Ridge nearly meets the trench, which is consistent with the oblique convergence of the Juan De Fuca and Gorda plates with the North American Plate.
In addition, NA13 displays a high-velocity anomaly beneath the Wallowa Mountains of the Columbia Plateau, similar to that imaged by Schmandt and Humphreys (2011). It is possible that this anomaly is a western promontory part the Precambrian North American lithosphere, given its connection to cratonic lithosphere shown in the right half of profile DD′. Alternatively, this anomaly could represent a phase of Farallon or Kula subduction from before 40 Ma, given its connection to a northern transition zone high-velocity anomaly, as shown in the center of profile DD′, that extends upwards to a depth of ∼250 km. We interpret the southern transition zone high-velocity anomaly, seen in EE′ and FF″ as more recently subducted fragments of the Farallon Plate. More precisely, a ∼1% high-velocity anomaly resides between 300 and 500 km depth range beneath the Snake River Plain and the northern Basin and Range (Figure 3), and is imaged in the EE′, FF″, GG′, HH′, and II′ cross-sections. Other S velocity models, such as Van der Lee and Nolet's (1997a) and Sigloch et al.'s (2008) image a similar feature, interpreted as related to Farallon subduction.
First noted by Zandt et al. (2004), a small high velocity lithospheric anomaly is observed in southern California (to the east of C′ in Figure 7). This anomaly, referred to as the South Sierra Nevada anomaly is observed by Obrebski et al. (2011) in their joint inversion extending to 400 km. In NA13, this anomaly appears to extend to ∼350 km depth. The lateral resolution of the Golos et al. (2018) model is limited to features larger than 400 km, and indeed this model does not show this feature.

Yellowstone Region
The Yellowstone hotspot and its related thermal anomalies are another significant component of the recent modification western United States' crust and mantle (Humphreys et al., 2000). Profiles GG′, HH′, and II′ examine the upper mantle of the Yellowstone region. Here, strong low velocity anomalies are present in the upper mantle from the Basin and Range to the Snake River Plain. These low-velocity anomalies are strongest in the upper 200 km beneath Yellowstone as mapped in GG′. The maximum slow anomaly occurs right beneath the Yellowstone Caldera and exceeds a 10% anomaly at a depth of ∼60 km. This maximum perturbation value is similar to that of Schmandt and Humphreys (2010) (−11% for V s ) and significantly larger than the −7% of Tian et al. (2011)'s. At a depth of 30 km, model NA13 displays an even larger peak low velocity of −13.8%, some of which might be attributable to velocity modeling interference from crustal structure. According to Cammarano et al. (2003), a −14% S velocity anomaly requires a temperature anomaly of 700 K if it is purely caused by heat. Such high temperature is above the solidus of peridotite at 50-100 km depth range (Goes & van der Lee, 2002;Thompson, 1992). Therefore, it is likely that partial melt occurs beneath the Yellowstone Caldera and the easternmost Snake River Plain, which is linked to magmatism in the area. The presence of partial melt can be further noted by investigating P-velocity anomalies. As NA13 is a S-velocity model, Lou (2013)'s XL13P P-velocity model is queried, revealing a −4.5% V p anomaly in this region, just slightly less of an anomaly then the −5%-7% P-wave velocity noted by Huang et al.'s (2014) local Yellowstone regional study, which estimated a melt fraction of ∼9% in the lower crust, by following the methods of Chu et al. (2010). Calculating partial melt based solely on seismic velocities includes assumptions regarding density, temperature, and attenuation. Several studies have used seismic velocity data to attempt to quantify partial melt deeper in the earth (Goes & van der Lee, 2002;Schutt & Humphreys, 2004), as compared to Huang et al.'s (2014) crustal study. Both of these upper mantle partial melt calculations estimate up to 1% partially molten. The NA13 and XL13P velocity anomalies are in line with those of Goes and van der Lee (2002) and Schutt and Humphreys (2004), reaffirming a 1% upper mantle partial melt beneath Yellowstone and the Snake River Plain. Based on the GG′ profile in Figure 7, this strongly low-velocity anomaly in the uppermost mantle beneath Yellowstone continues at depth as a tall −1.5% S velocity anomaly that extends through the transition zone and into the lower mantle. Nelson and Grand (2018) suggest that the Yellowstone hotspot has a plume origin, based on shear wave tomography. In the Nelson and Grand (2018) model, their hypothesized plume extends almost vertically, with a slightly southward dipping structure though the upper mantle and down to almost 1,000 km, where it then strongly dips to the southwest. Model NA13 might support this inference, but more strongly suggest a northward dipping structure though the upper mantle. Either way, the low velocities do not extend perfectly vertically through the transition, which shows high velocities beneath Yellowstone in the upper mantle and transition zone (Figure 7 transect GG′ and HH′). However, the bottom half of the transition zone and the lower mantle down to 1,000 km indeed show similar low velocities as imaged by Nelson and Grand (2018). The non-vertical nature of the tall low-velocity anomaly connected to Yellowstone and the associated high velocities in the TZ beneath Yellowstone and the northern Basin and Range Province suggest that if Yellowstone has a plume origin, the plume may have been deflected by subducted lithosphere.

Southwestern United States
Figure 8 displays five profiles cross cutting the Rocky Mountain front that separates the tectonically active western United States from the relatively stable central and eastern United States. On close inspection, the velocity boundaries in the mantle, do not coincide directly along the Rocky Mountain front, and instead deviate to the west along the JJ′ and KK′ profiles, and to the east along LL′, MM′, and NN′. The intricate velocity pattern can be observed on a map at 150 km depth (Figure 8). At this depth, three high-velocity features extend from the stable Great Plains to the western US through the Columbia Plateau beneath Idaho, through the Wyoming craton and Colorado Plateau, and through northwestern and central Texas. The northern two high velocity features are separated along the Snake River Plain, in association with the reworking of the mantle in relation to the Yellowstone hotspot.
The southern two high V s features are separated by lower velocities beneath the eastern margin of the Colorado Plateau (Figure 8), where the San Juan Mountains are. Similar P velocity anomalies were observed by Karlstrom et al. (2012). This separation in the uppermost mantle of high velocities in the Colorado Plateau region appears to align along the SW-NE trending border between the Northern and Southern Yavapai provinces, that can further be defined to the east by the Colorado mineral belt, which trends SW-NE just north of the San Juan volcanic region as noted in Figure 8. A joint inversion for 1D Vs profiles by Bailey et al. (2012) also shows the faster anomaly beneath the Colorado Plateau. They suggest that this thick, cool lithospheric root is at least partially responsible for the buoyancy and topography of the Colorado Plateau. NA13 models a faster anomaly to the south, than the Bailey et al. (2012) model. In addition, the NA13 anomalies extend slightly deeper, to ∼200 km, as opposed to ∼150 km in the Bailey et al. (2012) model.
Beneath southern Colorado, S velocities are less high, in alignment with higher expected temperatures, as modeled by Reiter (2008). Higher temperatures in this region would result in lower S velocities. The lowest S-velocities in this region are beneath the Southern Rocky Mountains and the San Juan volcanic, as seen in Figure 8, profile MM′. Indentations into the high velocity anomalies beneath the Colorado Plateau on the eastern side align with the Aspen Anomaly in the north, south to the San Juan anomaly, then to the Jemez lineament to the southeast, as seen at 150 km depth in Figure 8. The velocities in NA13 agree with those modeled in a 2D transect by Gao et al. (2004), with higher velocities in the uppermost mantle beneath northwestern New Mexico, and lower velocities in central and southern New Mexico related to the Rio Grande Rift. Near the southern New Mexico and Texas border, high velocities are imaged again (transect N-N′), which Gao et al. (2004) associate with a downwelling in the mantle, associated with convection processes. As seen in profile NN′, the mantle velocity boundary between the tectonically active western US and the North American craton appears to align with the surficial boundary.

Central and Eastern United States
Anomalies in mantle structure are also observed beneath the eastern United States, which is less well resolved in NA13 as mantle structure beneath the western USA. Resolution in the eastern United States is not as high as in the central and western US, due to data density, but features ∼300 km diameter can still be resolved ( Figure 5).
Maps of model NA13 in Figure 9 reveal heterogeneity in the mantle beneath the tectonically stable central and eastern US. Along the OO′ profile through the Great Plains close to the Rocky Mountains, the S velocity undulates with a slow-fast-slow-fast pattern from north to south, indicating complex variation at the western edge of the North American craton. The PP′ profile exhibits an east dipping low-velocity anomaly from New Mexico and northern Texas into the Great Plains to western Oklahoma in the Central Lowland. The low velocity anomalies just above the transition zone could most likely be related to mineralogical remnants of the Oklahoma Aulacogen or the thickening of the sub-lithospheric asthenosphere in conjunction with the thickening lithosphere. South of the Oklahoma Aulacogen which possibly formed during rifting of the proto-Atlantic (Burke & Dewey, 1973), the majority of Texas displays fast anomalies to 200 km depth as seen in the map view in Figure 9, representative of a more stable lithosphere, but is bordered to the west by lower velocities coincident in region to the Delaware Aulacogen.
Low velocities extend beneath Arkansas (profile PP′) from the bottom of the 200-km thick lithosphere and the 410-km discontinuity, and beneath the New Madrid Seismic Zone at the northernmost point of the Coastal Plain up until the Interior Low Plateaus (Figure 1). The three-dimensional structure is further explored in QQ′, which cuts through the NMSZ to the NNE. Here, the eastern Texas and Louisiana portion of the Coastal Plains display low velocities throughout much of the upper mantle. A northerly elongated slow anomaly hoovers above the transition zone and underlies the Mississippi Embayment and the New Madrid Seismic Zone, up until the Central Lowlands.
The mantle beneath the eastern Coastal Plains is characterized by low velocity anomalies throughout the upper mantle and transition zone (profile TT′). Strong low velocity anomalies are additionally observed beneath New England (profiles RR′ and TT′) and could be attributed to the thermal Northern Appalachian Anomaly (Menke et al., 2016), although alternate studies link a similarly placed low velocity anomaly to remnants of the Great Meteor hotspot track (Taylor & Fitzgerald, 2010). Two profiles along the MOMA (Wysession et al., 1996) (RR′) and FLED (Wysession & Fischer, 2001) station arrays (SS′) demonstrate velocity contrasts across the Appalachian Mountains between Precambrian lithosphere on the left of the profiles and that beneath the Paleozoic margins in the right of the profiles. The contrast is more pronounced along MOMA in the north than FLED in the south, but both are weaker, with contrasts of 2%dV s /100 km and ∼1%dV s /100 km than the ones across the Rocky Mountains. Profile TT′ crosses the path of the Virginia upper mantle anomaly, which can be seen around 300 km beneath the West Virginia and Virginia border. This anomaly has been suggested to be related to the ∼48 Ma Highland County, Virginia Volcanics (Mazza et al., 2014;, possibly related to a lithospheric delamination. Low velocity anomalies are also imaged in the transition zone and lower mantle beneath the east coast. These features were proposed by Van der Lee et al. (2008) to represent a wet upwelling, as the velocity anomalies cannot be as easily explained through variations in temperature or composition other than water. These low velocity anomalies extend down to the top of the lower mantle and possibly originate from dehydration of dense hydrous magnesium silicates in the subducting Farallon Plate.

Conclusions
The joint inversion of teleseismic S wave delay times, regional waveform fits, and constraints on Moho depth, resulted in an updated North American tomographic model, NA13. This model combines several decades of seismic data used for early North American seismic models (NA95, NA04, and NA07) and includes a substantial portion of USArray data. The joint model combines the merits of body waves' substantive lateral resolution and regional wave trains' substantive vertical resolution, and is superior to models based on one of these data sets. Model NA13 better reproduces the observed station-average delay time patterns between the coasts than a variety of tomographic models (Lou & van der Lee, 2014). It also sharpened velocity contrasts between the western, central, and eastern US.
In the western US, the subducting Juan de Fuca and Gorda slabs are clearly imaged, consistent with previous body-wave tomographic models (e.g., Burdick et al., 2012Burdick et al., , 2014James et al., 2011;Obrebski et al., 2010;Schmandt & Humphreys, 2010;Sigloch, 2011). Beneath the Snake River Plain and Yellowstone, very strong low velocity anomalies are imaged. In NA13, the peak anomaly beneath the Yellowstone Caldera is increased to −14%, enhancing the strong low velocity anomaly imaged by Humphreys (2010, 2011) and Tian et al. (2011). This large anomaly is too large to be solely attributed to thermal anomalies (Cammarano et al., 2003;Goes & van der Lee, 2002) implying that partial melting is likely present beneath Yellowstone. These strong low velocity anomalies at shallow depths beneath Yellowstone are associated with weaker, deeper low-velocity anomalies, suggesting a continuous plume conduit down to the lower mantle. This conduit is not vertical and appears to be deflected by higher velocity structures related to slab segments or continental lithosphere. The contrast between low velocities in the uppermost mantle beneath Yellowstone and those of the bordering Wyoming craton are over 17%.
Between the western and central US, the Rocky Mountain Front is defined by sharp velocity contrasts. The neutral anomaly contour within the mantle does not coincide with surface geological boundaries, deviating westwards in the north and eastwards at central latitudes. In the southernmost portion of the USA, mantle and crustal boundaries appear more aligned.
The complexity of the interface between the tectonically active and stable portions of the continent suggest that different processes and interactions are occurring between the surface and mantle, and cannot always be directly linked geographically. The relationship between surficial and upper mantle deformation is complex and highly dependent on regional and local variations in lithosphere strength. Seismic velocity anomalies reveal that the interaction of the North American craton and the tectonically active western United States can be influenced by local features such as Yellowstone and the Aspen Anomaly, as well as deeper features such as the topography of the 410 km discontinuity.
In the central and eastern United States, mantle is also heterogeneous despite being tectonically stable compared to the western United States. For example, deep low velocity anomalies were identified from beneath eastern Texas and Louisiana portion of the Coastal Plain to the Mississippi Embayment and the New Madrid Seismic Zone. The New England anomaly (Tian et al., 2011;Van der Lee & Frederiksen, 2005) and the eastern US anomaly (Van der  are imaged with low velocities within the upper mantle and the transition zone. Overall, NA13 improves upon previous tomographic models by using both surface and body waves over a wide range of distances to constrain mantle structure. This improves the amplitude and location placement of velocity anomalies, and extends tomographic regional models deeper, with reliable results below the transition zone and into the lower mantle.