Using Stochastic Point Pattern Analysis to Track Regional Orientations of Magmatism During the Transition to Cenozoic Extension and Rio Grande Rifting, Southern Rocky Mountains

The southern Rocky Mountains in Colorado and northern New Mexico hosted intracontinental magmatism that developed during a tectonic transition from shortening (Laramide orogeny, ca. 75 to 40 Ma) through extension and rifting. We present a novel approach that uses stochastic weighted bootstrap simulations of a large set of new and historical geochronology data to better understand how regional anisotropies responsible for focusing magma emplacement evolved through time. This technique can detect subtle trends in directional distributions, including multi‐modal orientations, and can be filtered from regional to local scales. Our results indicate that magmatism followed first the northeast trend of the Colorado mineral belt between 75 and 40 Ma and deviated afterward. These deviations vary depending on the scale of the analysis. At the smallest scale we evaluated (<75 km), the orientation of magmatism from 45 to 30 Ma rotated counter‐clockwise before aligning with the north‐south trend of the modern Rio Grande rift. Larger, regional‐scale analyses indicate magma centers between 40 to 35 Ma and 25 to 20 Ma were dominantly oriented southwest‐northeast, whereas magmatism between 35 and 25 Ma had north‐south orientation. The large areal footprint of magmatism and shifting regional patterns suggest that ancient zones of weakness in the North American lithosphere accommodated magma flow at different moments in time, rather than controlled by a retreating interface of the Farallon and North American plates.

batches (Coleman et al., 2004;Eddy et al., 2022;Gaynor, Coleman, et al., 2019;Samperton et al., 2015).Consequently, some studies have visualized space-time patterns of points of magmatism that represent rock samples with known coordinates and ages (Armstrong & Ward, 1991;Chapin et al., 2004;Glazner, 2022).However, identification of any aligned magma systems in these space-time point patterns is usually left to visual examination.Detecting directional patterns with the naked eye can be challenging, which introduces uncertainty in understanding how regional-scale spatial-temporal magma orientations shift through time.
The Laramide orogeny (ca.75 to 40 Ma) and the subsequent transition to continental extension are of significant interest for understanding continental, magmatic, and mineral deposit evolution in the western United States.A significant body of research focuses on syn-and post-Laramide processes in the southern Rocky Mountains of Colorado and New Mexico because they are associated with deformation and magmatism that evolved deep within a continental interior (Chapin, 2012;Jones et al., 2011;Lawton & McMillan, 1999;Mutschler et al., 1987;Tweto, 1975), the development of the Rio Grande rift (Abbey & Niemi, 2020;Burke, 1980;Chapin & Cather, 1994;Tweto, 1979), and historically important mineral deposits (Bookstrom, 1990;Lovering & Goddard, 1950;Tweto & Sims, 1963).Although previous studies have focused on deciphering the timing and regional stress states of Laramide compression (Caine et al., 2006;Erslev, 2005) or extensional forces during formation of the Rio Grande rift (ca.30 to 20 Ma) (Abbey & Niemi, 2020;Chapin & Cather, 1994;Ricketts et al., 2016), the regional stress state during the post-orogenic transition to extension (ca.45 to 30 Ma) remains unclear.Some studies suggest that magmas that formed during this transitional period have chemical compositions consistent with early lithospheric extension (Gibson et al., 1993), but other studies across the western Cordillera have argued that voluminous magma systems that formed during the post-orogenic transition are associated with only minor, localized extension that predates the formation of the modern Basin and Range and Rio Grande rift provinces (Best & Christiansen, 1991;Dickinson, 2002;Elston & Bornhorst, 1979;Gans et al., 1989).
We test if any regional-scale anisotropies that may have focused magmatism can be tracked during the transition from compression to extension in the southern Rocky Mountains.We compiled existing geochronological data with new U-Pb zircon dates, and analyzed regional orientations in magma crystallization ages through time with Fry analysis (or "Fry method"; e.g., Fry, 1979).The Fry method was originally developed for analyzing strain in rock samples, but it has received considerable attention in evaluating structural features that control mineral systems, specifically because it can detect subtle, potentially multi-modal patterns in directional distributions in point patterns (Carranza, 2009;Lisitsin, 2015;J. Vearncombe & Vearncombe, 1999;Zheng et al., 2020).Fry analyses are often applied once for a given point pattern, but applying this approach to a compilation of geochronological samples yields biased results with no uncertainties to assess model precision.Biases in geochronological data sets can be spatial (e.g., non-representative sampling bias) or temporal (e.g., inaccurate dates due to incomplete isotopic system preservation).For example, these data sets typically include analyses from a variety of geochronological methods, each with their own considerations in terms of accuracy, precision, and the geological context of their age (e.g., igneous crystallization vs. uplift).This aspect creates challenges for handling older, often imprecise, geochronological data: they provide some information, but should they be entirely removed from the analysis, or can we still use them in some way?In this study, we use weighted bootstrap resampling to help mitigate temporal and spatial biases (C.B. Keller & Schoene, 2012;C. B. Keller et al., 2015;Rosera et al., 2021), and to help quantify uncertainties in our analysis of magmatic orientations through time.

Paleoproterozoic Through Paleozoic History
The pre-Laramide (>75 Ma) history of the southern Rocky Mountains includes multiple episodes of continental crust formation, accretion, and deformation that resulted in distinct regional structural trends, and these inherited structural fabrics influenced subsequent magmatism and deformation (Frothingham et al., 2022).The oldest basement rocks exposed in the southern Rocky Mountains of Colorado and northern New Mexico are composed of juvenile continental provinces, the Yavapai and the Mazatzal, that accreted onto the southern margin (present coordinates) of the Archean Wyoming craton during the Paleoproterozoic (Figure 1a) (Shaw et al., 2001;Whitmeyer & Karlstrom, 2007).The accretionary provinces have generally east-northeast striking boundaries that are characterized by deformation fabrics ranging from relatively narrow suture zones to more widespread regions of deformation (Karlstrom & Humphreys, 1998;Shaw & Karlstrom, 1999;Whitmeyer & Karlstrom, 2007).Paleoproterozoic northwest-trending structural fabrics have also been documented in Colorado 10.1029/2023TC007902 3 of 26 (Tweto & Lovering, 1977), and interpreted to represent complex subduction zone geometries during growth of juvenile continental material (Figure 1a) (Jessup et al., 2005).Voluminous plutonism during the Mesoproterozoic (ca.1.4 Ga) was associated with reactivation of Paleoproterozoic northeast-trending structures (Shaw et al., 2001).Deep-seated north-to northwest-trending structures overprinted Paleoproterozoic fabrics during an approximately 1.1 Ga extensional episode (Karlstrom & Humphreys, 1998).Pennsylvanian uplift of the Ancestral Rocky Mountains may have occurred along north-striking zones of weakness formed during this Proterozoic extensional episode (Marshak et al., 2000).

Late Cretaceous to Quaternary History
Late Cretaceous to Quaternary igneous rocks occur throughout the southern Rocky Mountains as upper crustal plutons, subvolcanic porphyries and plugs, and volcanic rocks that range from local flows and domes to regional ignimbrites (Figure 1; e.g., Lipman & Bachmann, 2015;Mutschler et al., 1987;Steven, 1975).These magma systems developed during a well-documented, broad change in plate dynamics throughout western North America, including the ca.75 to 40 Ma Laramide orogeny and transition to post-Laramide extension (Chapin, 2012;Rosera et al., 2021;  Jessup et al. (2005), G. R. Keller et al. (2005), and Whitmeyer and Karlstrom (2007).Rio Grande rift outline modified after Hudson and Grauch (2013), Colorado mineral belt outline modified after Chapin (2012) and Tweto and Sims (1963), and igneous rock units modified after Horton et al. (2017).(b) Time-slice panels showing the spatial distribution of samples from rocks formed during the Laramide orogeny (80 to 40 Ma), post-Laramide transition (40 to 30 Ma), and extension and rifting (30 to 5 Ma).Area of points is proportional to the weight each sample received in our bootstrap resampling method.Time-slices for each 5-Myr age bin are presented in Figure S5 of Supporting Information S1.Tweto, 1975Tweto, , 1979)).Early work suggested that the subduction of the Farallon plate along the western margin of North America could be tracked by an inward progression of magmatism into the continent interior, followed by westward retreat at the end of the Laramide when the subducting slab "rolled back" or foundered (e.g., Coney & Reynolds, 1977).However, Glazner (2022) recently argued that interpretations relating magmatism to the advance and "rollback" of the Farallon slab are incompatible with large data sets of igneous rocks in the western United States.
Magma systems broadly associated with the Laramide orogeny range in age from approximately 75 to 40 Ma (Bookstrom, 1990;Mutschler et al., 1987;Rosera et al., 2021) and are mostly shallow porphyries and plutons that are almost entirely confined to a northeast-trending (ca.043 azimuth) swath referred to as the Colorado mineral belt (CMB) (Figure 1a).Near the northeast limit of the CMB, Laramide magmas are dominated by alkaline intrusions and minor calc-alkaline plutonic rocks (Bookstrom, 1990;Mutschler et al., 1987).Laramide calc-alkaline quartz monzonites and granodiorites are the dominant igneous rocks in the central CMB, including multiple generations of sills in Paleozoic sedimentary rocks.The ca.48 to 45 Ma Swan Mountain complex (Bryant et al., 1981;J. W. Keller et al., 2017) is composed of monzonites and a younger quartz monzonite porphyry, and contains some of the largest sills that intruded in the region during this timeframe.The largest Laramide plutons in the CMB include the Twin Lakes batholith, which was episodically assembled between 60 to 57 Ma and 43 to 40 Ma (Feldman, 2010) and the 38.8 Ma Montezuma pluton (Rosera et al., 2021).The Montezuma pluton is cut by numerous northeast-trending trachyandesite dikes that have isotopic compositions comparable to the host pluton (Rosera et al., 2023), but are of unknown age (Figure S1b in Supporting Information S1).
Magmatism between approximately 40 to 30 Ma was transitional between earlier Laramide magmatism and younger systems associated with Rio Grande rifting.The eruption of large volume regional ignimbrites represents a major shift in the style and power of magmatism relative to that of the Laramide orogeny.Ignimbrite activity and caldera-forming magma systems initiated around 37.3 Ma (Wall Mountain Tuff; Zimmerer & McIntosh, 2012) in the north-south "Sawatch trend" of calderas (Lipman & Bachmann, 2015).A variety of intermediate to silicic plutons intruded during this transitional episode, including the Mount Princeton batholith (Mills & Coleman, 2013), sills and laccoliths, such as the Cimarron intrusive suite in northern New Mexico (Figure S1c in Supporting Information S1) (Armstrong, 1969), and silicic dikes and small plutons related to magmatic-hydrothermal alteration and low-grade porphyry Mo mineralization in the Sawatch Range (Rosera et al., 2021).Mafic to intermediate local volcanic centers developed during this period as well, including the central Colorado volcanic field (McIntosh & Chapin, 2004), the Raton-Clayton volcanic field in northeast New Mexico (Scott et al., 1990), and at the Two Buttes magma center in southeastern Colorado (Davis et al., 1996).
The northern Rio Grande rift began to form between approximately 30 to 25 Ma (Abbey & Niemi, 2020;Landman & Flowers, 2013;Ricketts et al., 2016).Large normal faults coeval with Rio Grande rifting developed from central Colorado to near the Wyoming-Colorado border, and some of the rift-related uplift was accommodated along ancient structures, suggesting that structural inheritance influenced the shape and dynamics of the Rio Grande rift (Tweto, 1979).Estimates of the magnitude of extension in the northern Rio Grande rift range from approximately 1 to 13 km, or 8%-17% extension (Chapin & Cather, 1994;Cordell, 1982).Magmatism that developed during this time period is compositionally diverse, including a suite of broadly alkaline magmas that intruded along the Front Range and Sangre de Cristo mountains in areas such as Cripple Creek, the Spanish Peaks, and Cerrillos Hills (Bachman & Mehnert, 1978;Kelley et al., 2020;Penn & Lindsey, 2009).Large volumes of calc-alkaline, intermediate to silicic magmas intruded the upper crust resulting in the formation of numerous caldera-forming magma centers in the San Juan mountains (Curry et al., 2021;Lipman & Bachmann, 2015;Steven & Lipman, 1976;Wotzlaw et al., 2013).Uplift along extension-related faults and differential erosion have exposed batholiths in the Never Summer igneous complex of northern Colorado, and the Questa-Latir complex in northern New Mexico (Gaynor, Coleman, et al., 2019;Jacob et al., 2015;Johnson et al., 1989).High-silica, F-rich magmas erupted as topaz rhyolites, including the Nathrop and Sugar Loaf domes in the northern Rio Grande rift, and also crystallized in the subsurface where they formed leucogranites and porphyry Mo mineral systems (Bove et al., 1990;Christiansen et al., 1983;Gaynor, Rosera, & Coleman, 2019;Mills & Coleman, 2013;Rosera et al., 2021).

Zircon U-Pb Geochronology
Samples were collected from magma systems near the Rio Grande rift that lacked existing U-Pb zircon data, including the Swan Mountain intrusive complex, Sugarloaf dome from the Nathrop topaz rhyolite dome complex, the Cimarron intrusive suite, and the Cerrillos Hills (Figure 1; Figure S1 in Supporting Information S1).We also sampled a trachyandesite dike that cuts the Montezuma pluton in the CMB (Figure S1b in Supporting Information S1), because there is an overall lack of precise geochronology for intermediate intrusions in the area.All samples were crushed in a tungsten mill, sieved to less than 300 μm, and then concentrated to heavy minerals using a Wilfley table, Frantz magnetic separation and heavy liquids.All zircon samples analyzed in this study were subsequently annealed in a muffle furnace at 900°C for 48 hr (Mundil et al., 2004).
These samples were analyzed for LA-ICPMS U-Pb zircon geochronology following the methods described in Ulianov et al. (2012), with additional information provided in Text S1 of Supporting Information S1.For these analyses, natural reference zircon GJ-1 crystal #67, which has a reported 230 Th disequilibrium-uncorrected 206 Pb / 238 U age of 600.28 ± 0.16 Ma (Schaltegger et al., 2021), was used as a primary standard for the determination of the relative sensitivity factors.Natural reference material zircon Plešovice was used as a secondary standard for accuracy control (336.373± 0.074 Ma, 2s uncertainty; Widmann et al., 2019), and analyses of this zircon yielded a value of 336.93 ± 0.43 Ma during the period of analyses (n = 48; MSWD = 1.12).
A subset of samples was analyzed using CA-ID-TIMS U-Pb geochronology using the methods detailed in Text S1 of Supporting Information S1.EARTHTIME ET100 Ma standard solution analyzed during the period of these analyses yielded a 206 Pb/ 238 U date of 100.176 ± 0.006 Ma (MSWD = 2; n = 22), within uncertainty of the recently reported inter-laboratory calibrated value of 100.173 ± 0.003 Ma for this solution (Schaltegger et al., 2021).For all analyses, uncertainty of weighted mean 206 Pb/ 238 U ages is given in the X/Y/Z or X/Z notation, where X = analytical error only/Y = with uncertainty of spike calibration added/Z = with uncertainty of decay constants added (Schoene et al., 2006).

Geochronology Database Compilation
A database of previously published ages for Late Cretaceous to Miocene intrusive rocks, volcanic plugs, and silicic calderas in the southern Rocky Mountains was compiled for this study (Table S1 in Supporting Information S1).We selected samples dated with U-Pb zircon TIMS, LA-ICPMS, SIMS (or SHRIMP), K-Ar, Ar-Ar, Rb-Sr, and fission track zircon.We excluded fission track apatite dates from the data set because the low closure temperature of the system generally does not reflect igneous crystallization ages.An attempt was made to preserve as much data as possible, but some samples were excluded because no uncertainties were reported.In order to ensure that this database was internally comparable, we recalculated the data so that all samples of individual isotope systems used the same decay constants, parental isotope compositions, the same monitor mineral values (when applicable), and finally included propagation of uncertainty associated with decay constant uncertainties.Further details are provided in Text S2 of Supporting Information S1.

Overall Approach and Assumptions
We analyzed the geographic distribution of Late Cretaceous to Miocene magma centers with the Fry method (Fry, 1979;Hanna & Fry, 1979).A single Fry analysis consists of taking a set of n original points, copying that entire set of points, and then transposing that copied set such that the centroid (or origin) is centered over one of the original points (Figure 2, see steps 1 and 2).A new copy and translation are generated for all n original points in the data set, resulting in a point cloud of n 2 − n points (Figure 2, step 3).The resulting Fry plot enhances subtle directional patterns and/or clustering in the original data set.The distance and orientation between any two points (point-pair) in the Fry plot can be measured, and measuring the orientations and frequency of orientations for all combinations of point-pairs is commonly used to construct a rose diagram.Each sector of the rose diagram represents the frequency of azimuths measured between point-pairs in that sector (Figure 2b).
Uncertainties introduced by imprecise or thermally perturbed geochronological data, as well as those from sampling and preservation bias, can be partially mitigated with a weighted bootstrap resampling approach (C.B. Keller & Schoene, 2012;C. B. Keller et al., 2015;Rosera et al., 2021).This approach uses random resampling with replacement (bootstrapping) to generate inferences about the original samples.The original sample data can be assigned weights that represent the probability that a given sample will be drawn in a given bootstrap sample.In our study, each sample in the data set was assigned a weight, w, following the same procedure as Rosera et al. (2021).That is, w is inversely proportional to (a) the spatial density of samples in the data set and (b) the analytical uncertainty of each sample.This approach effectively reduces the influence of data from areas that are oversampled or associated with imprecise geochronological data, respectively.
We use a weighted bootstrap approach to resample points representing potential igneous crystallization ages (i.e., geochronological samples) in different time slices.In each time slice, we draw a bootstrap sample and conduct a Fry analysis to generate rose diagrams representing the directional distribution of all point-pairs in the Fry plot.The bootstrap technique allows us to repeat this process to get a sense of uncertainty in the directional distribution of points in each age bin, and therefore can be used to better understand signals of shifting spatial magma patterns.This stochastic approach to Fry analysis is rarely applied in geological studies, but a similar approach was introduced by McNaught (2002) to estimate strain parameters.Our methodology adds sample weighting and focuses on detecting directional distributional patterns by comparing orientations of all point-pairs in the simulated Fry plots.
Sample coordinates from our geochronological compilation and new U-Pb zircon ages were reprojected to Universal Transverse Mercator zone 13N before processing.Stochastic directional distribution models of magmatism through time were constructed as follows: 1.A synthetic age was drawn for each sample in the full data set by drawing from a Gaussian distribution with the mean analysis age and one standard deviation (total uncertainty).2. The samples were placed into 5-Ma age bins between 80 and 5 Ma according to their synthetic ages.3. Within each bin, a subset of samples was selected where the probability of a sample being included was equal to sample weight w.Subsets were selected to match the size of each bin following grouping in the second step, and sampling was done with replacement.4. The subset for each 5-Ma age bin was passed into a Fry analysis, the angle between every point-pair in each bin was calculated and frequency of point-pair angles was recorded into a histogram of 15° sectors (i.e., a rose diagram).5. Steps 1-4 were repeated 1,000 times.6. Summary statistics are calculated for each 15° sector by using the frequency data for point-pair orientations.
The spatial density was calculated in ArcGIS 10.3 and attached to the sample data to calculate weights.All other analyses were completed in R v.4.1.1 (R Core Team, 2021).Fry points were determined using the SpatStat package (Baddeley et al., 2015) by incorporating the tool into a custom weighted-bootstrap script developed in R (Rosera & Gaynor, 2023).
A critical assumption of our approach is that the geochronological point sampling technique is representative of magmatism through time.The weighted bootstrap approach attempts to mitigate these issues by limiting spatial biases that could be related to local over-sampling or variations in erosion and preservation.If an intrusion is on a map but it has no corresponding geochronological data in our data set, then it cannot be included in the point sampling model.Extrapolating the points to areas on a geological map is a possible solution, but it is not a straightforward process.Digital local-scale maps (e.g., 1:24,000 to 1:48,000) are scarce, so we would have to rely on regional-or state-scale maps that overly lump intrusive units that are actually discrete magma bodies.Furthermore, texturally and compositionally similar intrusive rocks that might be combined for mapping purposes are not required to be the same age.This is exemplified by the "Lincoln porphyry" in the CMB, which is a label for megacrystic quartz monzonite to granites that range in age from ca. 67 to 43 Ma (Rosera et al., 2021;Wallace, 1995), and the term has even been applied to 35.6 Ma porphyries near Turquoise Lake (Craig, 1980;Rosera et al., 2021;Ruleman et al., 2020).Consequently, it should be kept in mind that the point pattern approach we present represents a snapshot of new and historical geochronological sample points, and we emphasize that there are many opportunities to further build out this data set.

Varying Models by Filtering Analytical Methods and Interpoint Distance
We constructed multiple models that deviate from the total data set in two ways: (a) some models excluded zircon and titanite fission track and Rb-Sr analyses, and (b) other models explore different maximum distance (2) Fry analysis is conducted by iteratively copying all points, and translating the copied points such that they are centered on a point in the original data set (green data set).(3) This process is repeated until every point in the original data set (shown in black) has been the center of a copied and translated set of points.The resulting point cloud is referred to as a Fry plot, and it can accentuate spatial anisotropies in the original point data.An example of one point-pair and the orientation between those two points is highlighted.(b) Rose diagrams showing the angle made for all possible point-pairs within the Fry plot shown in panel (a) Four rose diagrams show the effect of filtering the Fry plot at different scales, reflecting differences between regional (i.e., no distance filter) and more local point orientations (<75 km).Extent of spatial filters shown as dashed circles on Fry plot in panel (a).filters on the Fry plots before constructing the rose diagrams.Our first model used the entire geochronological data compilation.We nominally refer to this approach as "All Methods" because it includes samples from all analytical methods reported in Table S1 of Supporting Information S1.Fission track dates tend to have relatively high uncertainty and their dates represent relatively low closure temperatures.Similarly, Rb-Sr isochron ages, although included in our full data set for historical purposes, are not very reliable ages for crystallization (e.g., Bell & Blenkinsop, 1978;Tichomirowa et al., 2019), and therefore we generated a model that excluded these geochronological methods (hereafter, the U-Pb, K-Ar, Ar-Ar model).However, the fission track data represent a relatively large portion of the final data set (30%, 191/629), and in many places are the only geochronological information available (Figure 3b).We present both model types in this study to discuss the tradeoffs.
The directional distribution of magmatism and geological structures can vary depending on the spatial scale of the analysis.One benefit of Fry analysis in comparison to other directional distribution methods, such as standard deviational ellipses (e.g., Rosera et al., 2021), is that the point-pairs represented in a rose diagram constructed from a Fry plot can be filtered at various distances to differentiate between local and regional directional distribution patterns (S.Vearncombe & Vearncombe, 2002).This approach can lend further insight into how orientations of magmatism could vary with spatial scale.We explored filters ranging from 500 km down to 75 km, and ultimately selected 250-, 150-, and 75-km filters for our models.Filters greater than 250 km did not vary significantly from the unfiltered models.The smallest scale we consider, 75 km, was selected to approximately coincide with the scale of major uplifts in the study area, and a 150-km filter was selected as an intermediate value.We used the dmax function in SpatStat's Fry points tool to apply maximum distance filters (Baddeley et al., 2015; Figure 2).

Zircon Geochronology
All three samples from the Cimarron intrusive suite (CP15-01, CP15-02, CP15-03) contain distinct and abundant cores in CL-imagery, and have complicated internal textures (Figure S2 in Supporting Information S1).Zircon from the Cerrillos Hills (CH16-12) commonly have discrete cores, oscillatory zoning and occasional sector zoning.Zircon from sample MZ18-06, the trachyandesite dike that cuts the Montezuma pluton, commonly host apatite mineral inclusions and exhibit complicated sector zoning and rare zircon micro-inclusions.The Swan Mountain sill (SM17-11) contains some zircon with inherited cores, but more commonly there are no obvious cores in CL and the zircon have oscillatory and sector zoning.Zircon from the Sugarloaf dome of the Nathrop domes (SL19-04) were highly fractured, relatively cloudy, and did not yield significant CL response.Inclusions within zircon from the Sugarloaf dome included Th-phosphate, thorite and apatite (Figure S3 in Supporting Information S1), and zircon crystals analyzed by EDS commonly yielded measurable, high U and Th concentrations.

LA-ICPMS
LA-ICPMS spot data are compiled in Table S2 of Supporting Information S1, summarized in Table S3 of Supporting Information S1, and representative analytical spots are shown in Figure S2 of Supporting Information S1.All uncertainties are reported and presented as 2s (Figure 3a).Individual 206 Pb/ 238 U dates from the three samples from the Cimarron intrusive suite ranged from 30.3 ± 2.4 to 1704.8 ± 50.6 Ma (CP15-01), 29 ± 0.8 to 1656.2 ± 25.8 Ma (CP15-02) and 29.1 ± 1.0 to 943 ± 177 Ma (CP15-03).All three of these samples yielded a cluster of lower Oligocene 206 Pb/ 238 U dates, and all analyses with 206 Pb/ 238 U dates older than Oligocene-Eocene yielded Proterozoic 207 Pb/ 206 Pb dates (Figure S2 and Table S2 in Supporting Information S1).Sample CH16-12 from the Cerrillos Hills yielded 206 Pb/ 238 U dates ranging from 30.8 ± 6.6 to 1,883 ± 51 Ma.All of the Oligocene dates overlap within uncertainty, and analyses of cores all yielded Proterozoic 207 Pb/ 206 Pb dates ranging from 1,402 ± 80 to 1,808 ± 56 Ma.Analyses from sample MZ18-06 yielded 206 Pb/ 238 U dates ranging from 35.6 ± 2.6 to 38.9 ± 2.6 Ma, and did not include any spot dates older than Oligocene.Finally, 206 Pb/ 238 U dates from sample SM17-11 of the Swan Mountain quartz monzonite sill range from 41.0 ± 1.8 to 1,731 ± 29 Ma.

Database Compilation
The final data set we compiled contains 629 samples ranging in age from 81.2 ± 4.06 Ma to 5.53 ± 0.71 Ma (Figure 3b).A significant portion of the geochronological data has 2σ absolute uncertainties >3 Ma (approximately 26% of the samples).The most abundant type of analytical method is fission track (n = 191, or 30%; Figure 3b).Nearly half of the samples are K-Ar or 40 Ar/ 39 Ar dates, and U-Pb zircon methods represent only approximately 18% (115/629) of the data set (Figure 3b).High-precision TIMS analyses represent slightly less than half of the U-Pb zircon data (53/115), and about 8.5% of the total data set.Samples with dates between 80 and 40 Ma are almost entirely located in Colorado (Figure 1; Figure S5 in Supporting Information S1).Most of the northern New Mexico samples are <40 Ma, and the data set contains some dated samples in northeast Arizona and southeast Utah between 30 and 20 Ma (Figure S5 in Supporting Information S1).

Weighted Bootstrap Fry Analysis
We generated five models and the results for each were organized into 15 bins that span 5 Myr (Figures 4-7).For each model simulation, the azimuth for all point-pairs was binned into 15° angle sectors and a frequency of orientations was calculated.Hence, each angle sector contains the frequency of orientations determined from 1,000 simulations.Rose diagrams were constructed by calculating the median, 25th, and 75th percentile frequencies.In the following we refer to rose diagram frequencies as their northern-hemisphere azimuths (e.g., 030 to 045 refers to a northeast-southwest trending sector).

Effect of Removing Fission-Track and Rb-Sr Ages
To facilitate the comparison of model results with and without fission track and Rb-Sr geochronology, we isolated the frequencies of point-pairs that track different angle sectors and present them as boxplots (Figure 4).Although the median frequencies for angle sectors slightly vary, the relative frequencies are similar in both models: magmatism from 80 to 40 Ma has the highest mode oriented 030-045 (northeast-southwest), whereas Boxplots comparing frequency of point-pair orientations for three different angle sectors: 285 to 300 (west-northwest; green); 345 to 360 (north-northwest; blue); and 030 to 045 (northeast; red).The red sector includes the approximate regional trend of the Colorado mineral belt (azimuth 043).Note that both models show similar relative relationships for the frequencies through time, and both show that the 030 to 045 trend has the highest frequency for ages >40 Ma.The similar results of the two models indicate that fission track data do not significantly change our model, because their influence has been lowered in the weighted bootstrap simulations.magmatism younger than 40 alternates between modes oriented 285 to 300 (southeast-northwest) and 345 to 360 (approximately north-south).The most noticeable difference between the models is the lower number of points in the model that excludes samples dated with fission track and Rb-Sr methods.The 80 to 75-, 55 to 50-, and 45 to 40-Ma age bins have few samples (n < 10) in the model that excludes fission track and Rb-Sr data (Figure 4a).

Effect of Filtering Data at Different Spatial Scales
Fry analyses can be used to visualize spatial anisotropies at varying spatial scales by limiting the distance allowed between points (e.g., Figure 2).Hence, rose diagrams showing the frequency of point-pair orientations can be constructed to vary in time and by spatial scale (Figure 5).Note that the total number of points in the spatially filtered Fry plot is lower than what it would be without filtering (e.g., Figure 4a).Hence, in Figure 5 we report the total number of points in the filtered Fry plot, defined here as n f , rather than the input n.If no maximum distance is used, then: The highest numbers of points (n f ) are observed between 40 and 15 Ma (Figure 5a), corresponding to the Late Eocene to Miocene pulse of magmatism in the southern Rocky Mountains.A secondary peak in n f occurs between 70 and 60 Ma, during the early stages of the Laramide orogeny.
Decreasing the maximum distance of the Fry analysis for ages between 80 and 45 Ma results in a change in the directional distribution.Without any distance filter, samples between 80 and 45 Ma are mostly oriented 030 to 045, but at smaller spatial scales the distribution shifts, and in some cases becomes multi-modal (Figures 5b and 6; Figures S6-S9 in Supporting Information S1).Between 80 and 65 Ma, a more dominant east-west trend emerges at the smallest spatial scales analyzed here (<75 km; Figure 5; Figure S9 in Supporting Information S1).A north-south component emerges at smaller spatial scales for ages between 65 and 50 Ma.The 80 to 75-Ma bin contains only a few points (Figure 1b), and as a result, angle sectors are highly variable, where the frequency of orientation for some angle sectors ranges from 0 at the 25th percentile to 0.5 (i.e., the maximum) at 75th percentile (Figure S9 in Supporting Information S1).
The dominant orientation of magmatism in the 45 to 40-Ma age bin does not change as the search distance decreases.However, a minor component oriented 330-345 (northwest-southeast) does disappear at the smallest spatial scale (Figure 7).The directional distribution observed in the 40 to 35-Ma age bin rotates from 285 to 300 at the largest scale to 015 to 030 for scales limited to 250 km.When the 75-km limit is applied the dominant orientation in this age group is 000-015 (Figure 7).Rose diagrams in the 35 to 30-Ma age bin show remarkably little directional variation at the intermediate spatial scales we modeled (<250 km; Figure 7).
The frequency of orientations between 30 and 25 Ma does not exceed 10% for any angle sector, regardless of spatial scale considered (Figure 7).However, there are still detectable differences in median frequency of orientations, such as the north-northwest orientation observed in the <75 km model for this age bin (Figure 7).Rose diagrams between 25 and 20 Ma are oriented toward 105-120 if no distance filter is applied, but restricting the analysis to smaller scales highlights a bimodal distribution of orientations during this time frame, and the bimodal orientations differ in the <250-and <75-km spatially filtered models (Figure 7, top row).
Magmatism younger than 20 Ma also shows a tendency toward multi-modal directional distributions at smaller spatial scales.At search distances <250 km, the 20 to 15 Ma directional distribution has three distinct modes at 000 to 015, 045 to 060, and 300 to 315 (Figure 5, Figures S6-S9 in Supporting Information S1), which differs from a more unimodal directional distribution in the unrestricted model (Figure 5; Figure S7 in Supporting Information S1).

U-Pb Zircon Ages of Pre-Rio Grande Rift Intrusions in Colorado and New Mexico
Most of the geochronological samples analyzed in this study yielded protracted U-Pb age spectra far beyond what is possible by igneous crystallization, and therefore determining a crystallization age for the studied samples requires interpretation of the spectra.Protracted zircon age spectra can be caused by prolonged magmatic zircon growth (autocrysts), inclusions of zircon from previous pulses of a similar magmatic locus (antecrysts) or inclusions of ancient zircon material unrelated to the sampled pulse of magmatism (xenocrysts) (J. S. Miller et al., 2007).The ability to assess the cause of this dispersion is commonly driven by the level of analytical precision used for the U-Pb analysis (Gaynor et al., 2022).
Most of the samples analyzed by LA-ICPMS (5 out of 6) have overlapping clusters of Oligocene-Eocene zircon dates and anomalously older, and largely discordant, ages (Table S3 in Supporting Information S1).Many of these 207 Pb/ 206 Pb ages match those of exposed Proterozoic basement rocks in the region (i.e., mostly 1.7 to 1.4 Ga; Whitmeyer & Karlstrom, 2007), and are clearly xenocrystic.We interpret the crystallization ages of the studied intrusions using the youngest plateaus of LA-ICPMS ages.However, these young plateaus yield weighted mean ages with MSWD values >1 (1.3-7.8) and are over-dispersed.The relatively low temporal precision of LA-ICPMS spot analyses is not sensitive enough to discern between models to explain this dispersion, such as subtle Pb-loss or antecrystic inheritance (Gaynor et al., 2022).As a result, we calculated weighted mean ages for the younger plateaus with conservative uncertainties expanded by  √ MSWD , which is a simple model for accounting for overdispersion of ages (e.g., Vermeesch, 2018;Table S3 in Supporting Information S1).
The Swan Mountain quartz monzonite sill (SM17-11) yielded a relatively simple spectrum of ages determined via CA-ID-TIMS, with five of the seven grains overlapping within uncertainty, and two older grains which likely reflect xenocrystic inheritance.For this sample, we interpret that a weighted mean of the overlapping plateau of ages reflects crystallization of the intrusion at 43.446 ± 0.023/0.026/0.053Ma (n = 5; MSWD = 0.62).The age spectrum determined via CA-ID-TIMS for the trachyandesite dike from Montezuma is significantly more protracted, and the lack of a single obvious plateau precludes the ability to use a simple calculation for weighted mean age and uncertainty, even with a  √ MSWD expansion factor.Instead, we calculated a weighted mean and use a conservative uncertainty from the spread of dates (Δt/2, e.g., Rosera et al., 2021) to get a 206 Pb/ 238 U age of 38.123 ± 0.246/0.249/0.269Ma for sample MZ18-06.For all following discussion, we will use the CA-ID-TIMS ages for these two samples rather than the LA-ICPMS ages, because of the superior accuracy and precision afforded by CA-ID-TIMS for U-Pb measurements.
Zircon from the Sugarloaf dome topaz rhyolite lack concordant ages and yielded Th/U values (9.82-11.35;Table S4 in Supporting Information S1) far above typical igneous zircon (Schaltegger & Davies, 2017).These grains do not represent crystallization age and are probably affected by inheritance of Th-rich accessory minerals that were not completely removed by chemical abrasion, which would result in excess 230 Th relative to U, and because 230 Th is an intermediate daughter product of the 238 U-206 Pb decay chain, this leads to an excess of 206 Pb and reverse discordance.We do not use these data further and instead rely on the previously published Ar-Ar geochronology from the Nathrop domes for geostatistical analyses (McIntosh & Chapin, 2004;Nelson, 2021).

Comparison to Existing Geochronological Data
Quartz monzonite sills around Swan Mountain were previously dated with Rb-Sr geochronology (biotite-whole rock isochron; 45 Ma, no uncertainty provided; Simmons & Hedge, 1978), and more recently a sample collected east of Breckenridge yielded an 40 Ar/ 39 Ar biotite age of 45.27 ± 0.12/1.02Ma (J.W. Keller et al., 2017).The sample we collected is from the northernmost exposure of the sill complex and is distinctly younger than previous 40 Ar/ 39 Ar ages (43.446 ± 0.023/0.053Ma; Figure 3a), which suggests it crystallized after sills to the south cooled below Ar closure.The new age for the northern part of the Swan Mountain sill overlaps recent high-precision zircon ages for a granodiorite sill collected near Leadville (43.55 ± 0.12/0.14Ma; Rosera et al., 2021); however, radiogenic isotopic data indicate that these two magma systems are not directly related (Rosera et al., 2023).
New CA-ID-TIMS U-Pb zircon ages for trachyandesite dikes that cut the 38.8 Ma Montezuma pluton indicate that they intruded 400 to 900 ka after pluton assembly (38.123 ± 0.246/0.249/0.269Ma; Figure 3a).These data confirm that the dikes are more closely related to assembly of the Montezuma pluton during the late stages of the Laramide magmatic activity, rather than representing mafic to intermediate magmatism related to younger extension.The dikes are propylitically altered, which suggests that hydrothermal alteration that is observed in parts of the Montezuma pluton is younger than 38.1 Ma.
Previous K-Ar geochronology suggested the Cimarron complex crystallized at 34.7 Ma, while other reports suggested the intrusions are as young as 26 Ma, although the origins of that age are not clear (Armstrong, 1969;Kish et al., 1990).New LA-ICPMS U-Pb zircon dates indicate the complex was assembled between 30.48 ± 0.22 Ma to 31.48 ± 0.26 Ma (Figure 3a), indicating that magmatism within the Cimarron complex predates the formation of the neighboring Questa-Latir magmatic center (ca.28 to 19 Ma; Gaynor, Coleman, et al., 2019).
Andesite porphyry sills in the Cerrillos Hills were previously interpreted to have intruded between approximately 39.1 to 34.3 Ma (Ar-Ar hornblende dates, but no uncertainties provided; Maynard, 2014).New U-Pb zircon geochronology from an andesite porphyry sill in the Cerrillos Hills yielded an age of 33.05 ± 0.39 Ma, suggesting that previous ages are either too old, or that andesite sills formed over a longer time period than previously described.

Comparison of Models With and Without Fission-Track and Rb-Sr Data
Stochastic point pattern analysis of a data set of standardized igneous crystallization ages provides a new way to visualize how the regional orientation of magmatism in the southern Rocky Mountains shifted from the Laramide through younger continental extension.Removing fission track and Rb-Sr ages from the input data has little effect on the directional distribution patterns observed in the "All Methods" model (Figure 4b).This is except for time periods where the number of samples falls below approximately 10, and then uncertainty of our analysis increases.This impact is most notable for the 55 to 50-Ma age range (Figure 4), perhaps because the numerous zircon fission track ages in this age range likely represent cooling related to epiorogenic uplift, which occurred between approximately 54-46 Ma (Abbott et al., 2022), or by partial thermal resetting of Late Cretaceous to Paleocene intrusions by younger magmatism (e.g., Rosera et al., 2021).Regardless, we conclude that removing the fission track data has little effect on the directional distribution results for time periods with well-documented magmatic activity in our study area (i.e., approximately 70 to 55 Ma and 45 to 15 Ma; Figures 5a and 6a).Consequently, we focus our discussion on the "All Methods" model and its subsets (i.e., filtered at 75, 150, and 250 km) because it allows for analysis of a larger data set, capable of better detecting subtle shifts in magma orientations through time.

Orientation of Magmatism During the Laramide Orogeny
The distinct northeastern trend of magmatism before 40 Ma is most apparent in our model when no distance filter is applied; it is a regional magma pattern that, in part, defines the CMB (Figure 1).Rose diagrams constructed from Fry plots filtered below 250 km show a shift in magma orientations (Figures 5 and 6).This is most notable in the early Laramide peak in magmatic activity between 70 and 60 Ma.In the model restricted to distances below 75 km, the directional distribution of magmatism shows modes oriented between 015 to 030 and 075 to 090 for magmatism between 70 and 65 Ma, and it shows no strong anisotropy from 65 to 60 Ma (Figure 6).The Laramide magmatism in the CMB is locally complicated at smaller scales, but aligned along a northeast trend at regional scales.This is in agreement with previous research, that observed local structural control of magmatism in the CMB by fault zones (Chapin, 2012;Lovering & Goddard, 1950;Tweto, 1975;Tweto & Sims, 1963).Identification of these smaller-scale anisotropies in Laramide magmatism differs from previous directional distribution modeling that only identified the large-scale, approximately 043 trend of the CMB, with standard deviational ellipses (Rosera et al., 2021), and demonstrates the utility of using the Fry method for spatial analysis of magmatic patterns.
Recently, Abbott et al. (2022) proposed that shortening related to Laramide compression in the southern Rocky Mountains ceased by ca.67 Ma, and that an epiorogenic pulse of uplift between 54 and 46 Ma occurred in response to a mantle drip event.Our data set indicates that, if such an event occurred, it happened during a time period of relatively low magmatic activity (Figures 1 and 4), and that it does not correspond to any regional shifts away from the dominant northeasterly trend of magmatism that existed before and immediately after (i.e., 45 to 40 Ma) the proposed drip event (Figure 4).
Various hypotheses exist for explaining the northeast orientation of magma systems that formed within the CMB during the Laramide orogeny.These include magmatism following discontinuous shear zones (Tweto & Sims, 1963), the formation of secondary mantle convection cells between the Farallon and North American plates (Jones et al., 2011), or fluid ± melt leakage along a northeast tear or hinge in the Farallon plate as it subducted (Chapin, 2012).More recent seismic tomography has identified north-dipping structures in the North American lithosphere that are sub-parallel to the CMB (MacCarthy et al., 2014).There is also a growing body of literature that supports widespread aqueous metasomatism and conductive cooling of the overriding North American lithosphere during flat-slab Farallon subduction (Dumitru et al., 1991;Farmer et al., 2020;Humphreys et al., 2003).These recent observations, in conjunction with new geochronological data led Rosera et al. (2021) to suggest that the northeastern orientation of magmas between ca.75 and 40 Ma in the CMB was the result of crustal anatexis driven by channelization of slab-derived fluids along ancient, translithospheric structures in the North American plate (Figure 8a).The stochastic Fry analyses presented in this study suggest that any major lithospheric structures that focused magmatism during Laramide compression could have an en échelon pattern, because the directional distribution of magmas at smaller scales (<75 km) is offset from the northeast trend of the CMB (Figure 6).This could indicate that magmatism at this time followed complex translithospheric structures, but more high-resolution geophysical data are required to better understand these structures.It is evident that major northwest-southeast trending lithospheric structures in the southern Rocky Mountains were not major pathways for magmatism during Laramide compression, perhaps because northeastern movement of the Farallon plate relative to the fixed North American plate sealed ancient preexisting northwest-southeast structures, thereby inhibiting upward ascent of slab-derived fluids (Figure 8a; Rosera et al., 2021).

Transition to Rifting, 45 to 20 Ma
It is well accepted that post-Laramide magmatism in the region is characterized by the formation of voluminous calderas, sub-volcanic plutons and batholiths, and eventually syn-rifting magmatism (e.g., Chapin, 2012;Lipman & Bachmann, 2015;Mutschler et al., 1987).However, examinations of how ancient regional structures may have focused magmatism during the post-Laramide transition and subsequent Rio Grande rifting are scarce.Furthermore, some tectono-magmatic models suggest that rollback of the Farallon plate exerted a control on spatio-temporal patterns of magmatism, but examination of the directional patterns predicted by these models is usually left to the naked eye.Importantly, this model predicts that there should be a detectable directional distribution of magmatism sub-parallel to the retreating North American-Farallon lithosphere interface, such as in backarc settings.In the following section, we examine the implications of our analyses on predictions made by tectono-magmatic models that attempt to explain post-Laramide magmatism, extension, and Rio Grande rifting.

Regional-Scale Processes
Between 45 and 40 Ma, magmatism in the southern Rocky Mountain region was largely confined to central Colorado and oriented along the northeast trend of the CMB (ca.043; Figure 7).However, the orientation of magmatism shifted dramatically between 40 and 35 Ma, and again between 35 and 30 Ma, regardless of what scales are considered (Figure 7).At the largest scale we analyzed, magmatism between 40 and 35 Ma was broadly aligned toward 285 to 300, or southeast-northwest, along a transect that follows the modern northern limits of the San Juan volcanic field and extending into southeast Colorado and northeast New Mexico (Figure 1).This orientation of magmatism becomes dominant again between 25 and 20 Ma, after Rio Grande rifting initiated (Abbey & Niemi, 2020).Importantly, regional magmatism between 35 and 25 Ma lacks this southeast-northwest alignment, which suggests that processes controlling the regional orientation of magmatism fluctuated through time and do not follow a simple long-term pattern.
The shifting patterns of magmatism from 40 to 20 Ma could have important implications for models that rely on Farallon plate "rollback" to explain broader magma patterns in the region.For example, Davis et al. (1996) suggested magmatism at ca. 37 Ma, including the Two Buttes center, was in part related to deformation or foundering of the subducting Farallon plate beneath thickened lithospheric keel in eastern Colorado.This model is similar to many others that relate Late Eocene-Oligocene magmatism to the "leading edge" of the subducting Farallon plate as it advanced deep into the North American continent and eventually "rolled back" (Best et al., 2016;Bookstrom, 1990;Chapin, 2012;Coney & Reynolds, 1977); however, the retreating edge model becomes problematic if it is applied across the entire area affected by the Laramide orogeny in the western United States (Glazner, 2022).Rollback of a southeast-northwest oriented front of magmatism is also unlikely given our directional analysis because this regional southeast-northwest alignment of magmatism is only evident between 40 and 35 Ma, and again between 25 and 20 Ma, after the onset of Rio Grande rifting (Figure 7, left column).Regional magmatism between 35 and 25 Ma was dominantly oriented north-south (Figure 7, left column), and there is no obvious southeast-trending directionality that reflects a retreating "edge" of the Farallon plate as it decoupled from the North American plate.Gibson et al. (1993) suggested that a K-rich magma system distal to the modern Rio Grande rift axis (up to 300 km) formed in response to an early phase of deep lithospheric, but not necessarily crustal, extension.Examples of these systems include the ca.37 Ma Two Buttes, as well as in the Spanish Peaks and the Raton-Clayton volcanic field (Davis et al., 1996).However, these magma centers are part of a regional southeast-northwestern trend of magmatism that is oblique to north-south structures of the modern Rio Grande rift (Figures 1 and 7).Thus, if this time period does represent early extension, it was not likely in a similar stress field as dominantly east-west oriented extension that formed the younger Rio Grande rift.
Hence, slab rollback and early east-west oriented extension models cannot fully explain fluctuations in the regional orientation of magmatism from ca. 45 to 20 Ma.A simpler model suggests that the observed variations are due to widespread magmatism that episodically ascended along pathways related to ancient lithospheric structures.A growing set of observations suggest that flat-slab subduction of the Farallon plate during the Laramide orogeny may Figure 8. Maps comparing regional stochastic Fry analysis (rose diagrams with no distance filter applied), samples, and major structural and lithological features from Figure 1a.Paleoproterozoic structures are shown in dark gray-green and Mesoproterozoic major normal faults are shown in red.To the right of each map is a simplified cartoon block model showing the geodynamic framework for the time period.Center of rose diagrams is plotted near mean center of samples in each time slice, but adjusted slightly for ease of visualization.Colored points on block models are hypothetical and shown to highlight where upper crustal magma systems form in each time slice.(a) Laramide compression.Arrows on map show possible range of Farallon plate motion relative to fixed North America after Doubrovine and Tarduno (2008).Block model shows shallow subduction of Farallon plate that simultaneously conductively cools and metasomatizes the base of the North American lithosphere.Melting primarily occurs where fluids ascend deep lithospheric structures (blue lines) and into shallower levels.Fluid movement is inhibited along structures oriented nearly perpendicular to the direction of movement of the Farallon plate.(b) Post-Laramide transition and early ignimbrite flare-up.Replacement of the Farallon plate with hot asthenosphere drives widespread melting of metasomatized lithospheric mantle.Partial melts in the lithospheric mantle are transported laterally and then up major structures.Removal of compressive forces related to interactions between the Farallon and North American plates allows magmas to ascend along northwest-southeast oriented structures (red lines).Local high flux of magmatism drives extension oriented approximately northeast-southwest. (c) Continued partial melting of mantle source volume and lithospheric heating lead to Rio Grande rift formation.Magmas continue to ascend along major northwest-southeast structures, as well as more north-south structures related to rifting.
have simultaneously conductively cooled and metasomatized the base of the overriding North American lithosphere.Foundering of the Farallon plate may have exposed large regions of the North American metasomatized lithospheric mantle to hot asthenosphere, thereby driving spatially extensive partial melting (Farmer et al., 2008;Lake & Farmer, 2015).Volumetric estimates suggest that approximately 306,000 km 3 of lithospheric mantle-derived basaltic magma was required to drive the 37 to 23 Ma ignimbrite flare-up in the San Juan volcanic field, which corresponds to a capture area with a radius of 220-350 km, depending on the fraction of partial mantle melting involved (Lake & Farmer, 2015).This model requires lateral transport of magmatism to explain the formation of voluminous magma centers, such as the San Juan volcanic field's caldera clusters (Farmer et al., 2008;Lipman & Bachmann, 2015).This model varies from the "leading edge" model in that partial melting of the mantle may have occurred over large regions, and variability in local magma activity could be a function of the degree of metasomatism of mantle source rocks, local dynamics between upwelling asthenosphere and sinking Farallon plate material, and the availability of permeable pathways to transport magmas both laterally and vertically through the North American lithosphere.
We suggest that regional magmatism between approximately 40 to 20 Ma in the southern Rocky Mountains represents the combined effects of melting driven by prolonged exposure of a large area of metasomatized lithospheric mantle to hot asthenosphere, followed by magma movement along pre-existing zones of weakness in the North American lithosphere (Figure 8).The regional southeast-northwest alignment of magma centers between the early ignimbrite flare-up (40 to 35 Ma; Figure 8b) and again during Rio Grande rift formation (25 to 20 Ma; Figure 8c) is subparallel to Paleo-and Mesoproterozoic structures that were reactivated during the Laramide orogeny (Tweto, 1975).Recent structural analyses by Frothingham et al. (2022) observed that surficial fault traces and ductile fabrics in Colorado have a dominantly southeast-northwest orientation, and that these fabrics continue at depth where they can act as focal planes for magmatism.It is not clear why this southeast-northwest magmatic fabric was not present between 35 and 25 Ma (Figure 7, left panels).However, it could be related to some combination of: (a) dynamics between far-field stresses in the North American lithosphere; (b) increased permeability along ancient north-south zones of weakness, perhaps as northeast-directed compressive forces related to Farallon underplating waned (Figure 8); and/or (c) the shape of the asthenosphere/North American lithosphere interface where partial melting occurred, including potential vertical and lateral variability in the degree of mantle metasomatism that pre-conditioned melt generation (Lake & Farmer, 2015).

Local Scale Processes, <75-km Distances
At the smallest scale we evaluated (<75 km), magmatism between 45 and 30 Ma appears to rotate counterclockwise, beginning with the 030-045 trend of the CMB, and ending in a north-northwest orientation closer to that of the modern Rio Grande rift (Figure 7, highlighted in purple).We suggest that these data represent a gradual transition away from syn-, or late-orogenic magmatism in the CMB through early stages of local extension leading up to Rio Grande rifting.Regional studies of western North America suggest that many areas with post-orogenic voluminous magma production, such as the Cordilleran ignimbrite flare-up (ca.40 to 20 Ma), experienced local extension (Armstrong & Ward, 1991;Gans et al., 1989;Meyer & Foland, 1991).In the southern Rocky Mountains, the start of the ignimbrite flare-up corresponded to heating and anatexis of the deep crust, the development of calderas and porphyry Mo-F deposits, and rapid uplift within the north-south trend of the Sawatch Range (Gaynor, Coleman, et al., 2019;Lake & Farmer, 2015;Mills & Coleman, 2013;Rosera et al., 2021Rosera et al., , 2023)).We hypothesize that the observed counter-clockwise rotation away from the northeast trend of the CMB between from 45 to 30 Ma in <75-km filtered models reflects local extension and uplift driven by deep crustal heating by mantle-derived magmas (Figure 8b).Similar to uplift elsewhere along the flanks of the Rio Grande rift, uplift may have been partially accommodated by reactivation of Proterozoic structures that may have also experienced movement during subsequent tectonic events, including the Laramide orogeny (Marshak et al., 2000;Tweto, 1975Tweto, , 1979)).Hence, the directional distribution analysis in this study detects magmatism related to early, local, extension that occurred millions of years before the development of modern Rio Grande rift.

Data Fidelity and Future Applications
Tracking the progression of intracontinental extension leading up to rift basin development is challenging because: (a) it is difficult to date early extensional structures, (b) there could be limited sediment accommodation during early extensional phases, and (c) uplift along rift flanks erodes pre-rift rocks.The transition from extension to rifting therefore suffers from data fidelity challenges; preservation could be limited to a few isolated areas over a much larger system, and geochronological data for building absolute timing of events can vary in precision.
The stochastic point pattern analysis of magma systems presented in this study attempts to limit local oversampling and adds weight to high-temporal precision data.We suggest that such approaches are necessary to better quantify uncertainties and assess the "state" of our knowledge regarding key tectono-magmatic conditions leading up to rifting.A benefit of this approach is that it allows us to include historical geochronological data by decreasing its weight in the analysis, rather than only focusing on a small subset of high-precision analyses (e.g., U-Pb zircon TIMS ages account for only 8.5% of the data set we used, Figure 3b).In many places, lower precision historical geochronological data may be the best or only information we have, and therefore our method provides a way to retain those data and spatial coverage, albeit with lower weights on the final statistics.Furthermore, analytically demanding, high-precision data are usually acquired to test questions relevant to local magma systems.For example, many of the high-precision geochronological data available in the southern Rocky Mountains have been used to aid in parsing out subtle differences in pluton assembly (e.g., Gaynor, Coleman, et al., 2019;Jacob et al., 2015;Tappa et al., 2011), testing volcano-plutonic connections for local magma systems (e.g., Curry et al., 2021;Mills & Coleman, 2013), or to bracket mineralizing events in mining districts (Rosera et al., 2021).The approach presented in this study demonstrates a new way to leverage these high-precision data sets without eliminating historical information.In essence, it attempts to retain as much useable data as possible to address data "fidelity" issues, and it can be used to visualize or test inferences regarding spatial-temporal patterns in magma systems during transient tectonic episodes.Future studies can continue to build larger, modern geochronological data sets to further refine our understanding of regional tectono-magmatic processes.
Future directional distribution analysis of magmatism can further refine the methodology presented in this study.For example, extension of the northernmost Rio Grande rift was relatively minor (1-13 km; Cordell, 1982), and consequently our study did not factor in palinspastic restoration of the point data.Application of this tool to other terranes with regional-scale deformation may require integrating kinematic models to accurately capture regional orientations of magmatism at a given time period.Finally, the tool we presented relies on point locations of samples and does not extrapolate that information in space, such as within the extent of a mapped pluton, or along the strike of dikes.Integrating information from geological maps into these stochastic models could potentially enhance directional distribution analyses, especially at more local scales.

Conclusions
This study presents a new way to integrate modern and historical geochronological data to quantify and characterize regional anisotropies in magma patterns.An advantage of this technique is that it can resolve subtle and/ or multi-directional orientations in regional magma patterns, and therefore is an improvement over simple visual examination of two-dimensional point patterns for magmatism.Understanding the alignment of magma systems can enhance models that characterize Earth processes such as shifting regional stresses during tectono-magmatic transitions, the control of lithospheric structures on magmatism, and structural fabrics that could host unexposed magmatic-hydrothermal mineral systems.Stochastic Fry analyses reveal distinct shifts in the orientation of magma centers at various spatial scales from 80 to 5 Ma in the southern Rocky Mountain volcanic field.Magmatism was aligned along the regional northeasterly trend of the CMB before 40 Ma, and smaller-scale analyses indicate that the individual magma centers that define the mineral belt have locally more complicated orientations.The orientation of regional magmatism changed abruptly between 40 and 35 Ma: at regional scales, it was largely oriented southeast-northwest, but at smaller scales a north-south component emerged.The reemergence of southeast-northwest oriented magmatism during active Rio Grande rift development (25 to 20 Ma) challenges models that require a retreating southeast-northwest oriented interface between the Farallon and North American plates.Instead, we suggest that major regional zones of weakness acted as pathways for magma systems that were generated during spatially extensive melting of the lithospheric mantle following exposure to hot asthenosphere.At local scales (<75 km), magmatism between 45 and 30 Ma rotated counterclockwise away from an initial northeast trend along the CMB, and into the north-south alignment of the Rio Grande rift.Our data suggest that local extension and uplift along north-south structures occurred millions of years before development of the modern Rio Grande rift.Mélissa Ruiz assisted with mineral separation, Jean-Marie Boccard helped with mount preparation of the zircons and SEM imagery was supported by Agathe Martignier.Thanks to Max Webb, Tyson Smith, and an anonymous reviewer for providing constructive feedback.Conversations with attendees of "Records of continental rifting and their fidelity," Session T024 of the AGU 2021 meeting helped in the presentation of this work.This research was supported by the Swiss National Science Foundation project No. 200020_182007 to Urs Schaltegger.Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Figure 1 .
Figure 1.Maps of the southern Rocky Mountain region showing locations of Late Cretaceous to Miocene geochronological samples used in this study (literature and new U-Pb zircon ages) and simplified geological features.Literature samples include samples dated with U-Pb zircon (SIMS/SHRIMP, LA-ICPMS, TIMS), fission track zircon, Rb-Sr, K-Ar, and Ar-Ar.See text for further details.(a) Overview of study area, simplified geological features, and geochronological data used in this study.Proterozoic structures modified after Jessup et al. (2005), G. R.Keller et al. (2005), andWhitmeyer and Karlstrom (2007).Rio Grande rift outline modified afterHudson and Grauch (2013), Colorado mineral belt outline modified afterChapin (2012) andTweto and Sims (1963), and igneous rock units modified afterHorton et al. (2017).(b) Time-slice panels showing the spatial distribution of samples from rocks formed during the Laramide orogeny (80 to 40 Ma), post-Laramide transition (40 to 30 Ma), and extension and rifting (30 to 5 Ma).Area of points is proportional to the weight each sample received in our bootstrap resampling method.Time-slices for each 5-Myr age bin are presented in FigureS5of Supporting Information S1.

Figure 2 .
Figure 2. Examples showing how Fry plots and all point-pair rose diagrams are constructed.(a) Example workflow for generating one Fry plot (i.e., one realization) for geochronological samples with mean ages between 35 and 40 Ma.(1) Point locations of samples are extracted from and plotted in Universal Transverse Mercator (zone 13N) coordinates.(2)Fry analysis is conducted by iteratively copying all points, and translating the copied points such that they are centered on a point in the original data set (green data set).(3) This process is repeated until every point in the original data set (shown in black) has been the center of a copied and translated set of points.The resulting point cloud is referred to as a Fry plot, and it can accentuate spatial anisotropies in the original point data.An example of one point-pair and the orientation between those two points is highlighted.(b) Rose diagrams showing the angle made for all possible point-pairs within the Fry plot shown in panel (a) Four rose diagrams show the effect of filtering the Fry plot at different scales, reflecting differences between regional (i.e., no distance filter) and more local point orientations (<75 km).Extent of spatial filters shown as dashed circles on Fry plot in panel (a).

Figure 3 .
Figure 3. Rank-order plots showing results from new U-Pb zircon geochronology and database compilation.(a) U-Pb zircon results.Samples MZ18-06 and SM17-11 were analyzed with both LA-ICP-MS and CA-ID-TIMS methods, and results from both methods are shown for comparison.(b) Summary data for the entire geochronological data set used in this study.Vertical extent of bars represents 2s total uncertainty; bar color represents the analytical method.Pie chart in lower left shows the relative proportion of samples from the various analytical methods and isotopic systems.

Figure 4 .
Figure 4. Comparison of results for all point-pair orientations determined by stochastic Fry method described in text.Two model types are shown: those that include ("All Methods model") or exclude fission track and Rb-Sr ages ("U-Pb, K-Ar, Ar-Ar model").No distance filter was applied to the results.a. Boxplots comparing the number of simulated points used for each 5-Myr age bin.Dashed vertical line nominally marks n = 10, below which the point-pair orientations become noisy.b and c.Boxplots comparing frequency of point-pair orientations for three different angle sectors: 285 to 300 (west-northwest; green); 345 to 360 (north-northwest; blue); and 030 to 045 (northeast; red).The red sector includes the approximate regional trend of the Colorado mineral belt (azimuth 043).Note that both models show similar relative relationships for the frequencies through time, and both show that the 030 to 045 trend has the highest frequency for ages >40 Ma.The similar results of the two models indicate that fission track data do not significantly change our model, because their influence has been lowered in the weighted bootstrap simulations.

Figure 5 .
Figure 5. [Static placeholder image of Movie S1] Animation showing the effect of limiting the maximum point-pair distance for the "All Methods" model described in the text.(a) Boxplots showing the number of points included in the Fry plot for each age bin, across all four models.(b) Rose diagrams showing frequency of orientations for all point-pairs within the Fry plots.Animation progresses from models with no maximum distance filter (i.e., regional-scale) down to <75 km (local-scale).The changes in orientations of magmatic point-pairs highlight scale-dependent variations in the orientations of magma systems.Results are plotted at the 25th, 50th (median), and 75th percentiles.Each ring represents a 10% increment in frequency.Static images of the animation are available in Figures S6-S10 of Supporting Information S1.The animated version of this figure is available at (Movie S1).

Figure 6 .
Figure 6. Rose diagrams showing model results for age bins between 70 and 60 Ma.This timeframe represents the Late Cretaceous to Paleocene pulse of magmatism during the Laramide orogeny.Each ring represents a 10% increment in frequency.

Figure 7 .
Figure 7. Rose diagrams showing model results for age bins between 45 and 20 Ma.This timeframe corresponds to the Late Eocene stage of magmatism that was coeval with the Laramide orogeny, as well as post-orogenic transition and development of the Rio Grande rift.Each ring represents a 10% increment in frequency.Four panels highlighted in purple from the <75-km model show a counter-clockwise rotation from the northeast trend of the Colorado mineral belt to a north-northwest trend similar to the Rio Grande rift.