Insight into the Evolution of the Eastern Margin of the Wyoming Craton from Complex, Laterally Variable Shear Wave Splitting

Dense seismic arrays such as EarthScope’s Transportable Array (TA) have enabled high-resolution seismic observations that show the structure of cratonic lithosphere is more heterogeneous and complex than previously assumed. In this study, we pair TA data with data from the Bighorn Arch Seismic Experiment and the Crust and lithosphere Investigation of the Easternmost expression of the Laramide Orogeny (CIELO) to provide unprecedented detail on the seismic anisotropic structure of the eastern margin of the Wyoming Craton, where several orogens emerged from nominally strong cratonic lithosphere during the Laramide Orogeny. In this study, we use the splitting of teleseismic shear waves to characterize fabrics associated with deformation in the Earth’s crust and mantle. We constrain distinct anisotropic domains in the study area, and forward modeling shows that each of these domains can be explained by a single layer of anisotropy. Most significantly, we find a fast direction in the southern part of the Powder River Basin, which we refer to as the Thunder Basin Block (TBB), that deviates from absolute plate motion (APM). This change in splitting behavior coincides with changes in other modeled geophysical observations, such as active source P-wave velocity models, potential field modeling, and seismic attenuation analysis, which all show a significant change moving from the Bighorn Mountains to the TBB. We argue that these results correspond to structure predating the Laramide Orogeny, and most likely indicate a Neoarchean boundary preserved within the lithosphere.The Wyoming Craton is an Archean to Proterozoic block of lithosphere situated in the center of the North American continent (Figure 1) often divided into three main subregions: in decreasing order of age, the Montana metasedimentary province in the northwest, the Beartooth–Bighorn magmatic zone across the middle, and the southern accreted terranes in the southeast [1]. By ~2.5 Ga, all three subregions were cratonized and assembled as a distinct block of lithosphere [2]. Following cratonization, the Wyoming Craton’s interaction with the other cratons of Laurentia is debated. There is general agreement that terminal collision in Northern Laurentia between the Superior Craton, Hearne-Rae Craton (for location, Figure 1), and Slave Craton of northwestern Canada began earlier (~1.815 to 1.780 Ga) than in southern Laurentia between the Wyoming and Superior cratons (~1.750 to 1.700 Ga). It is unclear whether this represents one orogenic event (i.e. the Trans-Hudson Orogeny [THO]) or several discrete events, with some authors referring to the southern portion of the THO as the Black Hills or Dakotan Orogeny [3-5]. Following the formation of Laurentia, the Wyoming Craton was tectonically quiescent until ~80 Ma, when flattening of the Farallon slab initiated the Laramide Orogeny [6-9]. Cratons are assumed to be stable, thick lithosphere resistant to deformation or destruction under most circumstances [10], but across the Wyoming Craton, the crust deformed in basement-cored uplifts during the Laramide (i.e. the Wind River Range, Granite Mountains, Owl Creek Mountains, Bighorn Arch, Laramie Mountains, and Black Hills; Figure 1).The location of the eastern edge of the Wyoming Craton is still debated: some studies suggest that the craton extends through the Black Hills, based on the Archean age of rocks in the Black Hills and in the uplifts throughout the Wyoming Craton [2, 11], as well as the age of detrital zircons from the Black Hills [12]. However, there is an age difference of ~300 Myr between the Little Elk Granite in the Black Hills [11] and Archean volcanic rocks in the Bighorn Mountains [13]. Others posit a more westerly boundary along the eastern edge of the Bighorn Arch. This argument relies on the presence of a crustal-scale west-dipping seismic reflector that may reflect a Precambrian suture zone coupled with a magnetic contact east of the Bighorn Mountains [14]. Magnetotelluric studies image a highly conductive anomaly that extends from northern Canada through the Cheyenne Belt [15]. Bedrosian and Finn [15] term this the North American Central Plains Anomaly and posit a connection to the THO and assembly of Laurentia. This interpretation places the eastern margin of the Wyoming Craton east of the Black Hills. Given the disagreement between interpretations in these studies, key outstanding questions regarding the Wyoming Craton are as follows: Where is the eastern edge of the craton? Did the Laramide Orogeny affect the physical state (e.g. strong and intact or weakened and destabilized) of the crust and mantle lithosphere in the Wyoming Craton? How does the current state of the lithosphere relate to the craton’s formation and evolution? In this study, we utilize shear wave splitting (SWS) to investigate these questions.SWS can provide crucial insights into dynamic and static mantle. For instance, in the tectonically active western United States, a pattern of swirling fast directions indicates either a lithospheric drip under the region [16] or toroidal flow around the Gorda-Juan de Fuca slab edge [17]. In cratonic Australia, SWS fast directions spatially correlate with ancient deformational event strain orientations [18]. Fast directions and delay times can also vary at or across tectonic terrane boundaries [19], allowing for determination of a more precise boundary location relative to other methods. Previous splitting results in and around the Wyoming Craton indicate complex anisotropy that generally does not match asthenospheric deformation predicted by North American plate motion [20, 21]. In addition, small delay times within the craton increase eastward into the Proterozoic THO [20, 22]. SWS analysis near the Yellowstone hotspot shows plate motion in parallel fast directions to the northwest of the craton, with more variable fast directions in the Archean lithosphere of the craton [21]. Studies using different methods suggest that anisotropy varies with depth across the region [23, 24]. Preliminary SWS work focusing on the Bighorn Mountains found that fast directions do not generally align with plate motion but instead align with boundaries between crustal aeromagnetic anomalies, implying coherence between crustal and mantle structures and an Archean origin to observed fabrics [25]. That study additionally observed a change in splitting parameters east of the Bighorn Mountains, corresponding with the more westerly boundary of the Wyoming craton posited in Worthington et al. [14].In this study, we present SWS results from three networks (Figure 2). First, the Crust and lithosphere Investigation of the Easternmost expression of the Laramide Orogeny (CIELO), deployed from 2017 to 2019 [26]. Second, the Bighorn Arch Seismic Experiment (BASE), which was deployed for a year and a half starting in 2009 and formed the basis of the preliminary EarthScope SWS work in this region [25]. Finally, we use data from the EarthScope’s Transportable Array (TA), which was deployed across the United States in 2-year intervals. Combining data from these studies allows us to generate the most laterally detailed view of anisotropy beneath the eastern half of the Wyoming Craton to date. Our study suggests the existence of heterogeneous cratonic lithosphere that records complex deformation patterns.After cratonization, the Wyoming Craton was assembled with other cratons to form the composite craton Laurentia. However, there is a debate regarding terminal collision between the Wyoming and Superior cratons: using geologic exposure and structural data, some have argued that the suturing of the two occurred ~1.750 Ga [4, 27], whereas others have relied on dating of dike swarms and paleomagnetic data to suggest it occurred much later at ~1.715 Ga [5, 28]. Regardless, by ~1.7 Ga, Wyoming was part of the Laurentian core. The eastern boundary of the Wyoming Craton is the suture with the Superior Craton via the THO, though the exact location of this boundary is debated (Figure 1). To the north of the Wyoming Craton is the Great Falls Tectonic Zone, the boundary between the craton and the Great Falls Tectonic Zone is indicated by a northeast trending high conductivity anomaly [29] and a north-dipping seismic reflector [30]—which may be the remnants of the Proterozoic suturing of the Wyoming Craton with cratons to its northwest. To the southeast of the Wyoming Craton is the Cheyenne Belt (Figure 1), where cratonic Archean and Paleoproterozoic rocks abut younger Paleoproterozoic rocks of the Yavapai province [31, 32].During the Laramide Orogeny, the crust deformed in basement-cored uplifts across the Wyoming Craton. There are several hypotheses for how this deformation occurred. A common argument is that Precambrian faults or other preexisting weaknesses within the crust were reactivated, facilitating orogenesis deep within the continent [33-35]. Other models for Laramide orogenesis include Moho-penetrating faults defining discrete blocks [36, 37], buckling of the upper crust accommodated by thickening of the lower crust [38], upper crustal buckling forced by crustal detachment [39], and buckling of the entire lithosphere [40]. Recent seismological examination of the Moho under the Bighorn Arch via receiver function and deep crustal refraction analysis disproves hypotheses involving Moho-penetrating faults or thickening under the range by pure shear deformation [14, 41]. Mechanisms of crustal detachment and lithospheric buckling still stand. How these uplifts occurred has implications for the current state of the Wyoming Craton. For instance, if the lithosphere was heavily deformed during orogenesis through buckling, faulting, or metasomatism, the craton may still record evidence of that deformation, or it may have destabilized or delaminated—similar to the North China Craton, which saw subduction and potential destabilization in the Mesozoic [42].Regional and continental tomographic models [43-49] reveal several broad features in the Wyoming Craton: first, remarkably low velocities occur in parts of the western Wyoming Craton—likely due to interaction with the Yellowstone Plume. Second, the central and eastern portions of the craton feature a high-velocity anomaly extending below 200 km depth. Third, a swath of low velocities to the southeast of the craton occurs along the Cheyenne Belt. Several models also estimate a low-velocity anomaly beneath the Black Hills [43-45, 48, 49] extending to at least 100 km depth. Bedle et al. [43], Burdick et al. [44], Schmandt and Humphreys [47], Schmandt and Lin [48], and Shen and Ritzwoller [49] observe a NE–SW trending high-velocity region in the center of the craton, bounded by velocity lows to the west and southeast, with a clear slowing of velocity to the east below 100 km depth. Several studies suggest that this high-velocity block extending to ~300 km depth may represent a remnant of the Farallon slab [43, 47, 48]. Humphreys et al. [50] develop this hypothesis further, arguing that because xenoliths from northern Wyoming indicate a warmer geothermal gradient than predicted for Archean lithosphere, this high-velocity block relates to the emplacement of depleted oceanic mantle in the form of the Shatsky Conjugate beneath the Wyoming Craton. Dave and Li [45] argue that tomographic models show interaction of the Wyoming craton with the Yellowstone hotspot and suggest that flat-slab subduction initiated small-scale convection at the base of the lithosphere as opposed to the emplacement of depleted mantle. Both hypothesize, though, that the lithospheric root of the craton has been severely deformed and modified as a result of the Laramide Orogeny.A recent attenuation study has confirmed a thick block of lithosphere (upwards of 200 km) beneath the central Wyoming Craton [35], in agreement with the high-velocity block imaged by tomography. These results also show higher attenuation beneath mountain ranges than beneath basins in this region. The magnitude of the signal requires thicker and less attenuating lithosphere beneath the Bighorn and Powder River basins than the Bighorn Mountains and Black Hills. These variations in thickness may provide a mechanism for localizing stress during the Laramide, as stronger blocks could resist deformation and transfer stress to weaker adjacent blocks. Zhu et al. [35] refer to the high-velocity anomaly roughly coincident with the southern Powder River Basin as the Thunder Basin Block (TBB), see Figure 2—we utilize the same terminology here.This study utilizes data from the CIELO, the BASE, and the TA (Figure 2). CIELO was a 2-year (2017 and 2019), linear array deployed across the eastern portion of the Wyoming Craton, from the southern terminus of the Bighorn Mountains, through the TBB, and ending to the east of the Black Hills [26]. CIELO aimed to image the structure of the easternmost Wyoming Craton and better understand the relationship between mantle lithosphere and Laramide orogenesis. It utilized twenty-four Nanometrics Meridian Compact Posthole seismometers. Refer Ford et al. [26] for a detailed review of CIELO installation procedures and data quality.BASE was a year-and-a-half (2009 and 2010) experiment centered on the Bighorn Mountains, with stations to the west in the Bighorn Basin and the east in the Powder River Basin. The goal of this experiment was imaging the Moho of the Bighorn Arch and understanding the structural mechanism for uplift of foreland basement-cored arches. The broadband component of BASE had thirty-eight Guralp CMG-3T seismometers [41]. Previous studies that examined BASE data include the following: receiver function analysis by Yeck et al. [41], reflection and P-wave analysis by Worthington et al. [14], and local seismicity analysis by O’Rourke et al. [51]. Here, we provide the previously calculated, unpublished splitting results from BASE [25] in conjunction with new results from CIELO. The BASE and CIELO arrays overlap at the southern tip of the Bighorn Mountains and the westernmost portion of the TBB. Combined, the two experiments provide a higher spatial resolution view than previous studies of lateral splitting variations in the eastern half of the Wyoming Craton, which relied on TA stations or even fewer for pre-TA studies.Seismic anisotropy (i.e. directional wave speed dependence) results from macroscopic alignment of intrinsically anisotropic minerals. Shear waves passing through an anisotropic medium are split into a fast and slow phase that separate in time, with the time lag proportional to the strength of anisotropy and/or thickness of the anisotropic layer; this is known as SWS and is an unambiguous indication of anisotropy [52].We calculate splitting parameters using similar methodologies for all three arrays in an updated Splitlab [53, 54] using the rotation correlation method [55]. This is a grid search method which rotates the two horizontal seismogram components in 1° increments and time shifts them in 0.1-second increments. At each orientation and lag time, a cross-correlation between components gives a measure of the similarity of the waveforms and the maximum indicates the most probable fast direction and delay time of that split. This method does produce a systematic 45° misorientation at null and near-null backazimuths, but this misorientation is well understood and can be used to forward model single layers of anisotropy [56, 57]. We require signal-to-noise ratios of 5.0 or greater and a correction from initially elliptical particle motion to corrected rectilinear particle motion. Because agreement between different splitting methods is often used as a quality control metric, we calculate splits with two additional methods: the minimum energy and eigenvalue methods [58]. The degree of agreement between methods aids in determining the quality of splitting. We verify nulls—results with no evidence of splitting—with the splitting intensity method [59], which has values near zero for nulls.We use SKS, SKKS, and PKS phases, which each convert from P waves to S waves at the core–mantle boundary and thus isolate anisotropic signatures on the receiver side. Events recorded on the BASE and some TA stations were limited to Mw ≥ 5.75 and epicentral distances of 85–145°; for CIELO and other TA stations, events were limited to Mw ≥ 5.5 and epicentral distances of 90–130°. Results for BASE stations include 519 SKS, 109 SKKS, and 81 PKS phases. At CIELO stations, results include 282 SKS and 14 SKKS phases. TA stations had 716 SKS, 108 SKKS, and 58 PKS phases. BASE-affiliated researchers calculated splits using a filter of 0.02 to 1.0 Hz. CIELO-affiliated researchers used a range of passbands between 0.02 and 1.0 Hz to optimally isolate the signal. We note that splitting results were calculated by two different groups at different points in time. However, agreement between the results of the two groups (discussed in Section 3.1) suggested that a full reprocessing of results was not necessary.There are 709 splits and 180 nulls at the 35 BASE stations. At the 24 CIELO stations, there are 296 splits and 241 nulls. Finally, at the 39 TA stations, there are 882 splits and 331 nulls; because both groups use 9 TA stations in common, 47 splits and 10 nulls are calculated twice, but only one version is used in our analysis. Most measurements are made on data from events with westerly backazimuths (180–360°) corresponding to Pacific subduction zones (Figure 3). All splitting results are provided in online supplementary Data Sheet S1.In general, stations outside of the TBB have more splits. In online Supplementary Material Figure S1, we show the percent of nonnull splits for individual stations. The most obvious trend is that the southern TBB into the Black Hills has a lower percentage of nonnull splits than elsewhere. One possible explanation for this is an increase in seismic noise: Ford et al. [26] noted that the United States Geological Survey (USGS) catalog had 428 mine blasts during the CIELO deployment (primarily in the eastern TBB and likely an undercount given the USGS minimum threshold of ML 2.5), indicating substantial mining activity. O’Rourke et al. [51] used BASE and TA data to find 1563 mine blasts in the region. The combination of basin structure [41] and unspecified mining noise may increase the levels of ambient noise in the TBB.In online Supplementary Material Figure S2, we show fast directions and delay times according to latitude and longitude. Fast directions (online Supplementary Material Figures S2(a)–S2(b)) do not show much geographic dependence. There is a decrease in the number of splits east of 106°W, consistent with our observation that splits within the TBB have fewer splits. Most fast directions (52%) are within 20° of 90° or −90°, though there is significant scatter outside of that. Delay times (online Supplementary Material Figures S2(c)–S2(d)) also do not show any strong geographic dependence. The vast majority of splits (90%) have delay times of less than 1.0 second, which holds true across the entire region.We find good agreement in observations from different groups. The splits come from independent analysis by BASE and CIELO-affiliated researchers, with slight differences in events, magnitude ranges, epicentral distances, and passbands. To ensure these differences do not impact our interpretation of splitting, we compare splits from the different groups in two ways. First, we directly compare splits at the nine TA stations used by both groups (Figure 4). This direct comparison of splits at TA stations indicates good agreement. For instance, both groups observe a change in splits with backazimuth at I21A (Figure 4): fast directions for events from the southwest trend ~NW–SE, whereas those from the northwest trend more ~E–W. One of the key observations of our study is the marked variation in fast direction and delay time with backazimuth (online Supplementary Material Figure S3). This is consistent across networks. For example, I21A (TA), TB010 (CIELO), and BH3E (BASE) all have a transition in fast direction from NW–SE at southwesterly backazimuths to more E–W at northwesterly backazimuths (online Supplementary Material Figure S4); fast directions also rotate to a more NE–SW orientation close to 0° backazimuth. Overall, there is good agreement between results from the two research groups at TA stations and minor disagreement across networks. We conclude that the separate sets of splitting observations are consistent and may be interpreted together.Splits are often reported as station averages, which many researchers correlate with plate motion or geographic features. However, we prefer not to focus on station averages since averaging splits obscures backazimuthal variation (for reference, averages are shown in online Supplementary Material Figure S5). As shown in online Supplementary Material Figure S3 and Figure 4, our splits have significant variation with backazimuth. For clarity, and given the large number of splits, we also plot results from all stations in arbitrary 90° bins from 0° to 360° (Figure 5) for the initial analysis and depiction of major trends in the data. In the 0–90° bin, most splits trend NE–SW, with some variation near the Bighorn Mountains (Figure 5(a)). From 90° to 180°, CIELO stations shift to a nearly W–E trend, whereas BASE stations maintain a NE–SW trend (Figure 5(b)). This bin has the fewest splits from all networks, and hence, this difference between BASE and CIELO likely results from waves coming from different backazimuths and epicentral distances (as discussed in Section 3.3), rather than a fundamental disagreement between networks. Splits from 180° to 270° trend NW–SE to the west of the Bighorn Mountains and east of the Black Hills and a subset shifts to NNE in the eastern TBB. Splits in this bin at stations north of the Bighorns have a smaller delay time and moderate heterogeneity. Finally, the bin with the most splits (270° and 360°) displays the most variability in fast direction and delay time. In particular, we note that there are different fast directions between backazimuths of ~300° and 360°—we address this in Section 3.3. While fast directions show considerable variation with backazimuth, delay times do not. Delay times are clustered between 0.3 and 0.6 seconds, with a mean of 0.4 seconds (online Supplementary Material Figure S6).For a more robust analysis of variations within the larger dataset, we attempt to find trends within smaller backazimuthal bins. We begin by identifying bins that are substantially different but have enough data to compare (marked in Figure 3(b)): 130–170°, 220–270°, 290–330°, and 340–20°. For each bin, we identify the fast direction that occurs most frequently (the modal fast direction). All splits are plotted in Figure 6, color coded according to their deviation from the modal fast direction for that bin.From 130° to 170° (Figure 6(a)), the modal fast direction is ~20°, which is 47° counterclockwise from APM, with significant variation existing between the CIELO and BASE stations. However, we note that for both networks, there is only one split within the corresponding backazimuth range for most stations and these events vary in backazimuth and epicentral distance, with events for BASE coming from ~150° backazimuth and ~85° epicentral distance, whereas events for CIELO come from ~140° backazimuth and ~120° epicentral distance. For BASE stations, the modal fast direction is 28°, whereas for CIELO stations, it is 97°. Different ray paths most likely cause the variation, not a systematic difference between networks or analyses. In no other backazimuthal range, we do observe a systematic difference between CIELO and BASE results, as shown in Section 3.1 and online Supplementary Material Figure S4.The modal fast direction from 220° to 270° is 120° (Figure 6(b)), rotated approximately 60° from APM. Of the 541 splits in this bin, 145 splits have a fast direction that deviates more than 15° from the mode. Fast directions that deviate most from the mode are located within the TBB in the eastern part of our study region with a few to the west of the Bighorns. A similar deviation is observed in the 290–330° backazimuth bin (Figure 6(c)), where splits have an average fast direction of 90°, 23° clockwise from APM. This bin has the most splits, with 705. Of those, 207 had a deviation from the modal fast direction of 15° or more. As with the 220°–270° bin, the splits with the largest deviations are within the TBB, east of the Bighorn Mountains in a NW–SE trending elongated region. In the 340°–20° bin (Figure 6(d)), the modal fast direction is 50°—17° counterclockwise from APM. Fast directions for 61 of 153 splits deviate 15° or more from the average, with a general consistency of fast direction across the region. Fast directions in this bin are the closest to APM.As with our initial binning of splits, we examine delay times for the same bins as for fast directions (online Supplementary Material Figure S7). Again, we find no significant variation with backazimuth for our delay times, with a similar clustering between 0.3 and 0.6 seconds. Because of this, we will generally restrict further analysis to fast directions.We begin by highlighting the most important features of our results: first, there is a significant variation in fast direction with backazimuth, although there is general consistency at any given backazimuth (Figures 4-6). While backazimuthal variations can result from layered anisotropy, we will demonstrate below that the systematic backazimuthal variations are consistent with a single layer of anisotropy. Second, fast directions are generally not parallel to APM. Third, splits found in the TBB, deviate from other splits in the region (primarily between 130° and 330° backazimuth) but are locally consistent with each other. As we argue below, this indicates inhomogeneous lithosphere within the Wyoming Craton. Importantly, the anomaly in the TBB is only observed with the inclusion of datasets with finer spatial resolution than TA (which has an average station spacing of 70 km), afforded to us with the addition of BASE (average station spacing of 32 km) and CIELO data (average station spacing of 19 km).The fast directions we observe display clear variation with backazimuth, exhibiting a sawtooth pattern (Figure 7(a)). The rotation correlation method yields a systematic 45° misorientation in fast direction with 90° periodicity for single layers of anisotropy [56, 57]. Eakin et al. [56] developed an algorithm based on equations laid out in Wüstefeld and Bokelmann [57] for generating models that utilize this systematic misorientation. We use this algorithm to test for single layers of anisotropy that might explain the backazimuthal variations seen in our splitting results. The misfit between the sawtooth pattern and each observation is calculated, then summed. This summed misfit is normalized by the maximum summed misfit for the whole suite of models. We note that between backazimuths of 230° and 360° (where 89% of our results are), an obvious feature is a change in fast direction from ~125° to ~50° around a backazimuth of 275° (Figure 7(a))—any successful model should explain this observation.An obvious starting point is to test whether fast directions can be explained by plate motion. In this region, plate motion is oriented at roughly 247° [60]. A fast direction parallel to plate motion (gray line, Figure 7(a)) fits 74% of our fast directions within 20° (summed normalized misfit of 0.18); however, the plate motion model fails to capture the change in fast directions at ~275° backazimuth. We therefore rule it out for further consideration.Testing the full range of fast directions from 0° to 180° in 1° increments yields a best-fitting model with a fast direction of 88°, a delay time of 0.6 seconds, and a summed normalized misfit of 0.12 (blue line, Figures 7(a) and 7(b)). The best-fitting model is a global minimum, but between 65° and 110°, the misfit curve has a relatively low slope, which steepens around a misfit of 0.2 (Figure 7(c)). While the fit for much of the data is good (82%), a subset of the data is fit poorly (18%). We therefore divide the results into those that are within 20° of the 88°-Fast Direction (FD) model (1345 splits) and those outside that range (297 splits) and model both groups separately. A virtually identical model to the initial test fits the first group (87° and 0.6 seconds with a summed normalized misfit of 0.17), whereas a model with a fast direction of 52° and 0.8 seconds fits the second group (red line, Figures 7(a) and 7(b); summed normalized misfit of 0.13). Few results are not fit well by either model, for instance near 300° backazimuth and 160° fast direction. These results have slightly different epicentral distances (online Supplementary Material Figure S8), so this may be explained by anisotropy elsewhere along the raypath. However, a large majority of our results are fit well by one of the two best-fitting single-layer models.Multiple single-layer models can explain the observed backazimuthal variation in our results. Variations in epicentral distance are not sufficient to explain differences between the data best fit by the 88°-FD and 52°-FD models (online Supplementary Material Figure S8). Certain features seem epicentrally dependent: for instance, results near 150° backazimuth and many of the results with fast directions near 180°. These features are less well defined than the overall sawtooth pattern and do not explain the bulk of our observations. Results between 230° and 360° backazimuth (which clearly conform to a sawtooth pattern) come from several different epicentral distances.A better explanation for differences between the two populations of data (those fit by an 88°-FD or a 52°-FD model) is lateral variation in anisotropy across the region. We divide our region into 0.5° by 0.5° cells and determine whether the 88°-FD or 52°-FD model fits the data in each cell better (Figure 8). Smaller cells would result in many cells with no data, whereas larger ones would smooth over short length-scale variations. We model each cell regardless of the number of splits; the minimum number of splits per cell is 5, whereas the maximum is 134. We consider neither model a good fit for the data if the summed normalized misfit for each of them is larger than 0.2 (i.e. outside the region of low slope in Figure 7(c) and the point at which misfits for all cells display a sharp drop off in online Supplementary Material Figure S9). Otherwise, the model with the smaller misfit is used for that cell. Thirty-four cells (1097 splits) are fit well by the 87°-FD model—these are scattered throughout the study area and include the orogens. Those fit by the 52°-FD model (16 cells and 306 splits) are almost entirely in the basins. Basin depth does not appear to be a primary control on this distribution, as these fast directions occur at various sediment thicknesses (online Supplementary Material Figure S10). Only one cell with 50 splits was poorly fit by either model, so we consider these results overall to be robust.While these findings are robust, we note that SWS alone may not be sufficient to describe complex anisotropic regions. SWS is a path-integrated effect and therefore has no depth sensitivity. Previous work examining synthetic SWS has found that even in cases where there are multiple layers of anisotropy, SWS tends to mirror the uppermost layer of anisotropy [61]. Further, the expected relationship between multilayer anisotropy and the systematic misorientation of the rotation correlation method has not been well explored, as pointed out by Eakin et al. [18]. In addition, it has been demonstrated that SWS and surface wave studies do not always agree in regions with complex and vertically varying anisotropy [18, 62]. However, as discussed further in Section 4.3, our SWS results are congruent with other geophysical methods examining lateral changes in the lithosphere of the Wyoming Craton.While our finding of laterally variable seismic anisotropy agrees with previous studies of SWS with coarser station spacing, our key finding of single-layer anisotropy does not. Station-averaged splits show fast directions deviating from APM to the east of the craton, which has been used to argue multiple layers of anisotropy further east in the THO and Superior Craton [22]; our modeled 88° fast direction is close to trends seen to the east of the TBB, where splits are more variable. Another SWS study [63] shows station-averaged fast directions transitioning from E–W in the TBB to NW–SE near the Black Hills, then back to E–W in the THO. The E–W orientation in Yang et al. [63] matches well with our overall preferred single-layer fast direction of 88°, though our other single-layer model with a fast direction of 52° disagrees with their NW–SE splitting orientations seen along the margin of the Wyoming Craton. Though these previous SWS studies [22, 63], respectively, used the multichannel analysis [59] and minimum energy method [58], their reported fast directions within the Wyoming Craton are similar to our modeled fast directions, particularly the 88° model.Other methods have indicated multiple layers of anisotropy in our study area. Ps receiver function analysis shows evidence of multilayered anisotropy within the mantle lithosphere [64]. Surface wave analysis also indicates multilayered anisotropy throughout the North American continent [23]. Our SWS results do not require multiple layers of anisotropy. There is a persistent disagreement between SWS and other methods in regard to multilayer anisotropy. For instance, in cratonic Australia, surface waves [65] and Ps receiver functions [66] image multiple layers of anisotropy, whereas SWS can be fit with a single layer of anisotropy [18, 66]. Similarly, Yang et al. [62] attempted to model SWS results in the THO using the three-layer model of Yuan and Romanowicz [23] but found the three layers unnecessary to explain their results.SWS is a path-integrated observation: previous studies suggest that most anisotropy originates in the upper mantle, but the crust [67] and lowermost mantle [68] can contribute as well. Schulte-Pelkum and Mahan [69] used Ps receiver functions to estimate crustal anisotropy by removing contributions to the signal from horizontally layered isotropic structure: they found stronger crustal anisotropy in the western US transitioning to weaker anisotropy east of the Rocky Mountains. In our study area, this transition to weaker anisotropy occurs along the eastern edge of the Black Hills [69], indicating that crustal anisotropy may contribute a nonnegligible amount of splitting to our results. The overall effect of crustal anisotropy is likely small but requires further analysis (i.e. using Ps receiver functions from BASE and CIELO). However, we do not think that our modeled delay times (0.6 and 0.8 seconds) could be generated in the crust alone. Regarding the lower mantle, we did not have a sufficient number of SKS–SKKS pairs (which have divergent paths in the lowermost mantle) to determine whether there is a significant contribution to SWS from lowermost mantle depths. However, as demonstrated in Section 4.1, backazimuthal dependence of fast directions is well explained by a single-layer anisotropy. We thus consider it unlikely that a contribution from the lowermost mantle is significant to our results.Because our modeled fast directions do not match APM, we explore the possibility that there is instead anisotropy within the cratonic lithosphere. The western portion of the Wyoming Craton has been modified through interaction with the Yellowstone Plume: a previous SWS study [21] found roughly plate motion parallel splitting near the plume, with more complicated patterns of fast directions further away. It is possible that the pattern in the western end of our study region (with fast directions predominantly fit with a 52° model, 15° away from APM) may result from interaction with the plume and plate motion. However, this would not explain the variability in modeled fast directions we see further to the east. The penultimate tectonic event in the eastern Wyoming Craton was the Laramide Orogeny, which occurred ~60 Ma. Analysis of structural trends in exposed rocks suggests that Laramide shortening throughout the craton was northeast directed [70, 71]. The link between shortening and SWS is largely dependent on the deformation style. In the case of uniaxial compression resulting in pure shear, the olivine b-axis (the slow axis) aligns with the orientation of shortening [72, 73]. In the case of large strain, fast directions should be parallel or within 30° of ductile shear [74-76]. Therefore, fast directions are expected to be parallel to the strike of tectonic features, such as mountain ranges [72, 73, 77]. The Bighorn Mountains and Black Hills both roughly trend NW–SE, which is oblique to our modeled fast direction of 88° and nearly perpendicular to the 52° modeled fast direction; therefore, pure shear under uniaxial compression does not fit our observations. This is in line with prior studies that suggest pure shear thickening does not support Wyoming mountain ranges [41]. If a simple shear model is assumed, fast directions should be aligned perpendicular to the strike of tectonic features. Our modeled fast direction of 52° is broadly consistent with northeast-directed shortening. However, these modeled fast directions are scattered throughout the study region and do not generally appear under orogens, where such shearing would be most expected.Another candidate process that may have produced our modeled fast directions is shear between the Shatsky Conjugate portion of the subducting Farallon and the overriding North American Plate. There is some debate regarding the trajectory of the Shatsky Conjugate around the time of the Laramide Orogeny, with some arguing for motion roughly N15°E in a relatively narrow corridor [8, 50, 78, 79]. Others argue for a broader corridor with motion oriented at N65°E [80-82]. In either case, the Shatsky Conjugate would have been beneath eastern Wyoming around 50 Ma. Importantly, if the subducting plate were the main cause of anisotropy in the region, then there should be very little lateral variation in fast directions (in contrast to our modeling). We therefore argue that subduction of the Farallon slab and shear at the lithospheric base of the Wyoming Craton is not the primary cause of the anisotropy we observe.The results of our modeling can be divided into three broad regions (Figures 8(b) and 9): first, to the west of the Bighorn Mountains, most modeled fast directions are 52° (67%); second, between the Bighorn Mountains and the Black Hills (the TBB), both modeled fast directions are scattered throughout; finally, to the east of the Black Hills, modeled fast directions are almost entirely 88° (86%). Previous studies using BASE data [14, 41] observed a transition in the structure of the crust along the eastern edge of the Bighorn Arch, roughly coincident with our observed change in modeled fast directions (Figures 8(b) and 9). More specifically, Worthington et al. [14] observed a west-dipping reflection boundary and modeled a sharp magnetic contact immediately east of the Bighorn Arch (along the black, model line in Figure 9), which they hypothesized was a suture between the Wyoming Craton and Proterozoic orogens. Building on this interpretation, Kilian et al. [5] presented a model based on paleomagnetic data in which the Black Hills are an exotic block within the Black Hills/Dakotan Orogen, with the eastern boundary of the Wyoming Craton immediately east of the Bighorn Arch.However, Hrncir et al. [27] used previously published geologic observations and new age data from detrital zircons to argue the Black Hills occupied a position in close proximity to the rest of the Wyoming Craton throughout the Proterozoic and that the suture imaged by Worthington et al. [14] more likely resulted from a Neoarchean collisional event. This explanation of the suture observed by Worthington et al. [14] would be consistent with the North American Central Plains Anomaly (NACP in Figure 1), indicating the remnants of oblique subduction zones formed during the THO [15, 29]. As noted earlier, our change in modeled splitting behavior from NE–SW under the Bighorn Basin to more mixed fast directions under the Powder River Basin occurs in roughly the same area as the crustal transition argued for by Worthington et al. [14] and Kilian et al. [5]—refer Figure 9.A recent attenuation study [35] found variations in upper-mantle, path-integrated attenuation of teleseismic P phases recorded within the eastern Wyoming Craton as part of the CIELO study: generally higher attenuation (low Q) occurs under regional orogens (i.e. the Bighorn Mountains and Black Hills), and generally lower attenuation (high Q) occurs under the basins directly west of these orogens (i.e. the Bighorn and Powder River basins). We see a transition from NE–SW-oriented fast directions in the Bighorn Basin to mixed fast directions in the TBB, roughly coincident with the changes in attenuation (see especially the “Low” transition line trending WNW–ESE across the Powder River basin in Figure 9). Zhu et al. [35] modeled path-integrated attenuation beneath the TBB to argue for thick, strong lithosphere—in agreement with previously published seismic tomography studies [45, 47-49] which show high velocities extending to 300 km.Several hypotheses for these high velocities from seismic tomography exist. First, the Wyoming Craton may be decratonizing (i.e. losing its long-term stability); in this model, high velocities would correspond to cratonic lithosphere that is downwelling and undergoing erosion [45]. Second, the eastern half of the craton may have been destabilized and then restabilized following the Laramide Orogeny—here, the high velocities would correspond to the restabilized block [50]. And third, the high velocities may represent a locally thickened region original to the craton. Supporting the model of decratonization (now or in the past) is the path of the Farallon slab, predicted to have been beneath the eastern portion of the Wyoming Craton by ~65 Ma [50], providing a means to hydrate and destabilize the lithosphere of the Wyoming craton. Xenolith data from the northern Wyoming Craton indicate Eocene kimberlites sampled metasomatized mantle lithosphere, possibly from interaction with the Farallon slab during the Laramide Orogeny [83]. In the decratonization model, evolution of the Wyoming Craton since the Mesozoic reflects that of the North China Craton, which was destabilized during the Mesozoic by interaction with a subducting slab, as evidenced by thinned lithosphere, high heat flow, and Mesozoic–Cenozoic volcanism [42]. Decratonization of the Wyoming Craton would imply that the base of the lithosphere has either been deformed (e.g. a drip) or metasomatized. However, the pattern of splitting we observe in this study does not match models of dripping lithosphere, which predict a reduction in delay time at the center [16] and a radial pattern of fast directions mirroring mantle flow [84]: as noted, calculated delay times are consistent throughout the craton and while there are backazimuthal variations in fast direction, they do not exhibit a radial pattern. Zhu et al. [35] argue that the region of low attenuation they model beneath the TBB is coincident with high seismic wavespeeds imaged by tomography, indicating strong lithosphere, and by extension implying that the lithosphere of the TBB is not currently delaminating or dripping (i.e. weak lithosphere). Taken together, this evidence suggests that active decratonization is less likely than other explanations for our results.In the recratonization hypothesis, a layer of oceanic lithosphere known as the Shatsky Conjugate was attached to the base of the Wyoming Craton beneath the TBB around ~70 Ma after the original cratonic keel was thinned [50]. As noted earlier, the combination of low attenuation and high velocities in the TBB implies strong lithosphere that is not actively decratonizing, leaving open the possibility of recratonization (i.e. strong lithosphere). If the allochthonous and original cratonic lithospheres each had their own anisotropic fabrics, with fast axes offset from each other, the observed path-integrated splitting results would vary as a function of backazimuth due to the two layers of anisotropy but with a sinusoidal pattern different than the sawtooth one we observe. The predicted areal extent of the Shatsky Conjugate would be quite large, covering and extending beyond our study area. If the recratonization model were correct, we would expect to see two layers of anisotropy without many local variations and low attenuation across the region. We therefore argue that the recratonization model is not the best explanation for our results.Alternatively, the thickened, strong lithosphere indicated by low attenuation beneath the basins could be original to the craton, and the observed lateral changes in single-layer splitting the result of structures inherited during the craton’s formation. This model matches well with Bader [34], who presents evidence from measured Precambrian fabrics and the structural trends of Laramide orogens that preexisting Precambrian heterogeneity within the lithosphere-facilitated Laramide orogenesis. Hopper and Fischer [85] reported multiple seismic discontinuities from receiver functions within the Wyoming Craton and at its boundaries, which they speculate are related to the formation of the craton itself—bolstering the case for ancient structure preserved within the lithosphere. Given the dearth of tectonism in the region between the Trans-Hudson and Laramide orogenies [86], the rheological contrast between lithosphere beneath basins and lithosphere beneath orogens suggested by Zhu et al. [35] may date to the Proterozoic, unless subduction of the Farallon slab prior to the Laramide Orogeny weakened the locations of future orogens. In either case, variations in lithospheric strength between thicker and thinner lithosphere may have then led to the formation of the Black Hills and Bighorn Mountains [35]. The transition from high attenuation in the Bighorn Mountains to low attenuation in the TBB is roughly coincident with our modeled change from more NE–SW fast directions west of the Bighorn Mountains to mixed fast directions in the TBB (Figure 9). The transition in our modeled fast directions coincides with the location of a suture zone suggested by Worthington et al. [14]—refer Figure 9. Using crustal-scale seismic tomography and gravity modeling of the Bighorn Arch, they image a high-density, high-velocity layer (7.xx) under the Bighorn and Powder River basins that is missing under the arch. Worthington et al. [14] suggest that the 7.xx layer may have strengthened the crust below the basins, focusing strain under the arches. Additional evidence for an intact cratonic block comes from Bedrosian and Frost [29], who image clear high-conductivity belts around the margins of the Wyoming Craton, which they argue correlate to the major tectonic events of Laurentian formation. There are no corresponding high-conductivity belts within the craton, which suggests a relatively stable block. In summary, our splits exhibit patterns characteristic of single layers of anisotropy not aligned with APM and do not fit with predicted patterns of delamination. We argue that (given other evidence presented here) the lateral changes in anisotropy we observe are likely original to the craton and date to the Precambrian.In this study, we generate SWS observations from three separate temporary networks (BASE, CIELO, and TA), two of which are newly reported, to image variations in anisotropy along the eastern margin of the Wyoming Craton at unprecedented scale and resolution. We would not be able to image this heterogeneity without the smaller station spacing (generally less than 50 km) of these temporary deployments. Our results confirm the presence of complex, backazimuthally dependent splitting. We can adequately explain our observations with a model consisting of multiple single-layer blocks of anisotropy. Our modeled fast directions deviate from APM, ruling out asthenospheric flow as the primary cause of anisotropy in the region. In addition, observed delay times are larger than expected for crustal anisotropy alone, which suggests preservation of fabrics in the mantle lithosphere of the Wyoming Craton. Critically, our dense station spacing reveals lateral variability in the TBB. A region of mixed fast directions is coincident with the re-emergence of a high-density lower crustal layer moving from the Bighorn Mountains to the TBB and a transition from high to low attenuation in the same location. Taken together, these observations indicate inhomogeneous lithosphere within the Wyoming Craton. We argue that the variations in splitting are the result of preserved structure within the Wyoming Craton that predates the Laramide Orogeny, as other hypotheses predicting syn- or post-orogenic lithospheric deformation predict anisotropy different from what we observe. This implies a substantial portion of the Archean-aged cratonic lithosphere is preserved within the Wyoming Craton beneath the TBB. While the zone of mixed fast directions could correlate to the eastern edge of the Wyoming Craton, other studies suggest that it is more likely a suture internal to the craton inherited during cratonic formation. The Wyoming Craton (and other cratons) cannot be treated as monolithic, homogeneous blocks. Rather, it is important to consider how discontinuities and fabrics within the lithosphere may influence or be reactivated during later deformation.All seismic data were accessed through the Seismological Facility for the Advancement of Geoscience (SAGE) Data Management Center (DMC) and are freely, publicly available. All shear wave splitting measurements are included in the Supplemental Material.The authors declare no conflicts of interest.A. Birkey received funding as a postdoctoral researcher through the University of Delaware and as a graduate student through the University of California, Riverside. H.A. Ford was funded as faculty through the University of California, Riverside. M. Anderson was funded through the Washington Geological Survey. J.S. Byrnes received funding through the Northern Arizona University. M.J. Bezada was funded as faculty through the University of Minnesota, Twin Cities. M. Shapovalov worked on this project as an undergraduate at the University of California, Riverside, and did not receive funding.The authors thank all those who were involved in the deployment of the CIELO and BASE seismic experiments. The authors also acknowledge the helpful feedback of one anonymous reviewer and the assistance of the editors.All shear wave splitting results are provided in a spreadsheet. DOIs for all networks are in Supplemental Text 1. Additional figures for shear wave splitting results are shown in Figures S1–10.