Thermochronological transect across the Basin and Range/Rio Grande rift transition: Contrasting cooling histories in contiguous extensional provinces Open Access

The Rio Grande rift (RGR) and the Basin and Range (B&R) province of the southwestern United States (Fig. 1) are two of the most widely studied continental extensional provinces in the world (e.g., Van Wijk et al., 2018), and yet the timing of the onset of extension and cooling of exhumed blocks in the key area where the provinces meet in southern New Mexico has not been documented. Because of this, it is not clear whether they represent distinct regions of extension with a sharp boundary, or if there is a gradual transition zone between them, or if the entire region from the southeastern B&R to the RGR experienced coeval extension. This study was undertaken to evaluate the cooling histories of rocks along an E-W transect from southeastern Arizona to south-central New Mexico to address this question using low-temperature thermochronometric methods of (U-Th)/He dating of apatite and zircon and apatite fission track. These methods have been used to constrain cooling histories of exhumed crust in a wide variety of tectonic settings (e.g., House et al., 2003; Fitzgerald et al., 2006; Wildman et al., 2016).

There have been several different models for the causes of extension in the RGR, including extension synchronous with rapid growth of the San Andres transform, oblique extension related to rotation of the Colorado Plateau, collapse of over-thickened crust following the Laramide orogeny, and passage of a weak mantle plume beneath the RGR (e.g., Hamilton, 1981; Humphreys, 1995; Landman and Flowers, 2013; van Wijk et al., 2018; Abbey and Niemi, 2020). Existing models vary with regard to how the RGR and B&R provinces are related, and how, when, and where the RGR formed (e.g., Smith, 2004; Ricketts et al., 2016; Axen et al., 2018; van Wijk et al., 2018).

A critical constraint that is lacking is the timing of extension in these two provinces where they meet in southern New Mexico. Improved characterization of the relationship between the Basin and Range and Rio Grande rift in southern New Mexico is key to creating a robust, unified model for possible driving mechanisms and timing of widespread Neogene lithospheric extension. The three options for understanding the relationship between these two extensional provinces are (Fig. 2): (1) The Basin and Range and Rio Grande rift are treated as separate extensional provinces (e.g., Abbey and Niemi, 2018; Fig. 2A); (2) the Basin and Range and RGR are merely two expressions of the same extensional event separated by the Colorado Plateau and merging to the south into a coherent extensional province (e.g., Baldridge et al., 1984; Dickinson, 2006; Fig. 2B); (3) the RGR effectively ends at the southern margin of the Colorado Plateau, raising the possibility that the southern RGR is merely the eastern extension of the B&R, and the Basin and Range province effectively extends all the way to the Great Plains province in southern New Mexico (Fig. 2C).

The primary goals of this project are: (1) to constrain the timing of main phase extension in the southern RGR and southeastern Basin and Range; (2) to identify and interpret differences and similarities between the RGR and Basin and Range in southern New Mexico; and (3) to test conflicting models for the initiation and evolution of extension that rely on different driving forces. The data support a model in which the cooling histories of the southeastern Basin and Range are influenced by Paleogene magmatism and shallow crustal intrusions. In contrast, ranges in the southern RGR experienced fault-related exhumation and cooling that were coeval with fault-related cooling in the northern and central segments of the rift and lag behind the youngest ages from Basin and Range samples by 2–3 m.y.

To address these topics, spatial and temporal patterns in a robust data set of new single-grain low-temperature thermochronology dates are considered as a first-order constraint on delineating different times of cooling throughout southern New Mexico. Newly presented thermochronological data include apatite (U-Th)/He (AHe), apatite fission-track (AFT), and zircon (U-Th)/He (ZHe) dates and build on previous research by Ricketts et al. (2016) in the northern and central segments of the Rio Grande rift.

The Rio Grande rift and Basin and Range are two extensional provinces in the southwest United States (Fig. 1) that formed in the mid to late Cenozoic following the rollback and sinking of the Farallon slab at the end of the Laramide orogeny (e.g., Coney and Reynolds, 1977). Farallon plate subduction to the east under the western margin of North America initiated ca. 180 Ma (DeCelles, 2004). Circa 90 Ma, the subduction angle of the Farallon slab began to shallow, possibly related to subduction of the conjugate Shatsky Rise, a low-density ocean plateau (Dickinson and Snyder, 1978; Liu et al., 2010). Nevadan-Sevier arc magmatism waned at this time, and tectonic and magmatic activity migrated >1000 km eastward, heralding the onset of the Laramide orogeny (e.g., Coney and Reynolds, 1977; Saleeby, 2003; DeCelles, 2004). Despite the Laramide orogeny being known as amagmatic in the central United States, in the southwestern United States, magmatism occurred in three pulses, starting at 75 Ma (Amato et al., 2017), followed by early Paleocene and mid-Eocene phases, the latter of which began ca. 44 Ma (McMillan, 2004) and ended ca. 39 Ma (Creitz et al., 2018). During the middle-late Eocene (ca. 50–35 Ma), shallow subduction of the Farallon plate slowed as it began to tear and delaminate (e.g., Axen et al., 1993; Humphreys, 1995, 2009; Madsen et al., 2006; Ricketts et al., 2016). This exposed the asthenosphere to the hydrated and weakened lithosphere, leaving behind a keel of displaced continental mantle lithosphere that was removed through subduction erosion (Axen et al., 2018). Large-volume, caldera-style volcanism began in southern New Mexico ca. 36 Ma in the Mogollon-Datil volcanic field (McIntosh et al., 1992; Chapin et al., 2004) in what is known as the “ignimbrite flare-up” (Coney, 1978; Lipman, 1992).

Ignimbrite volcanism that migrated westward and was contemporaneous with migrating continental extension is attributed to rollback of the Farallon slab (Coney and Harms, 1984). Silicic volcanism lasted from ca. 60–55 Ma in the Great Basin (Fig. 1; Dickinson, 2009), ca. 35–18 Ma in San Juan volcanic field (Lipman et al., 1978), and from 40 to 24 Ma in the Mogollon-Datil and Boot Heel volcanic fields (McIntosh et al., 1992). Migrating volcanism coincided with the beginning of basin-range block faulting, formation of sedimentary basins of 1–3 km depth (e.g., Averill and Miller, 2013), and a sweeping trend in core complex development throughout the Basin and Range (Axen et al., 1993; Dickinson, 2002, 2009).

The southern RGR is adjacent to the Basin and Range in southern New Mexico. The RGR is a narrow zone of lithospheric extension that extends >1000 km from Colorado into northern Mexico (Baldridge et al., 1984). It trends roughly north-south, parallel to the axis of the Rocky Mountains (Keller et al., 1991; Pazzaglia and Hawley, 2004; Fig. 1.) The rift consists of a series of en echelon half-grabens separated by broad transfer zones. Extensional strain in the RGR increases to the south, from 8%–12% in the San Luis Basin (Kluth and Schaftenaar, 1994) to 50%–100% in the southern rift (Morgan et al., 1986). From the vicinity of Leadville, Colorado to central New Mexico, the RGR is a discrete, narrow feature that is only tens of kilometers wide (Morgan et al., 1986; Ricketts et al., 2016). Along this stretch, the rift is separated from the central Basin and Range province by the unextended Colorado Plateau. South of Socorro, New Mexico (Fig. 1), the RGR widens as it merges with the Basin and Range province to the west along the southern margin of the Colorado Plateau.

In the Basin and Range of southwestern New Mexico and southeastern Arizona (Fig. 3), the onset of extension is not well documented because many of the volcanic rocks associated with basins have yet to be dated, and exposures of Miocene and possibly older basin fill are limited (Mack, 2004). However, several inferences can be made regarding the initiation of extension. In the Boot Heel region of southwestern New Mexico, thick accumulations (250–1000 m) of rift fill suggest that much of it may be Miocene in age (Thompson et al., 1978; Heywood, 1992; Klein, 1995; Chang et al., 1999; Hawley et al., 2000; Shearer and Miller, 2000). Circa 26–28 Ma, ignimbrite outflow sheets covered large areas of the Boot Heel region, suggesting that at the time of their eruption, no significant extension-related topography existed (McIntosh and Bryan, 2000).

The timing of onset of extension in the southern RGR has been controversial. Seager and Clemons (1975) and Mack et al. (1994a) noted that a single half-graben in the Goodsight–Cedar Hills “depression” (Fig. 3) was filled with ca. 36 Ma ash-flow tuffs. These authors used this relationship to suggest that extensional faulting related to the RGR had initiated at this time, and that the tuffs were derived from the Mogollon-Datil volcanic field. The first voluminous mafic magmatism in the RGR was the eruption of the ca. 27 Ma Uvas basalt (Clemons, 1979), which is widespread in the southern rift and inferred to have resulted from lithospheric thinning and asthenospheric upwelling (McMillan et al., 2000). Both Ricketts et al. (2016) and Abbey and Niemi (2020) reported that the initiation of rifting in the northern and southern RGR occurred ca. 25 Ma and proceeded to ca. 15 Ma. Abbey and Niemi (2020) proposed that extension in the RGR was related to initial oblique strain, which was followed by rotation of the Colorado Plateau as originally proposed by Chapin and Cather (1994). Abbey and Niemi (2020) reviewed the occurrence of bimodal magmatism along the entire rift and noted that voluminous mafic magmatism did not occur until ∼10 m.y. after extension was detected using thermochronology, particularly when focusing on intra-rift magmatism.

Several competing geodynamic models account for the initiation of the Rio Grande rift. Some models suggest that RGR extension is driven by mantle forces, either through large-scale mantle upwelling based on mantle tomography (Moucha et al., 2008) or small-scale mantle convection at the edge of the stable Great Plains province (van Wijk et al., 2018). Other models rely on crustal forces to drive extension and include clockwise rotation of the Colorado Plateau (Hamilton, 1981; Chapin and Cather, 1994; Landman and Flowers, 2013), collapse of over-thickened continental crust (Cordell, 1978; Eaton, 1986), and initiation of a transform setting along the western margin of the North American plate (Dickinson and Snyder, 1978; Baldridge et al., 2006). Some authors attribute extension to detachment and foundering of a piece of the Farallon plate beneath the Rio Grande region in New Mexico, and this extension enhanced asthenospheric upwelling in the slab window (Ricketts et al., 2016). Abbey and Niemi (2018; 2020) suggest that initial extension in the rift involved oblique strain that was replaced by rotation of the Colorado Plateau to drive later stages of extension (e.g., Chapin and Cather, 1994).

Because the B&R and RGR are contiguous geographically, the boundary between them is not clearly defined based on surficial geology, and the nature and location of such a boundary between the two provinces remain poorly understood. The boundary between the RGR and Basin and Range south of the Colorado Plateau has been generally positioned east of the Cookes Range, Florida Mountains, and Deming basin, and west of the Mimbres Basin (Fig. 3), approximately aligned with the extension of the eastern edge of the Colorado Plateau to the north (Woodward et al., 1978; Mack, 2004; Villarreal-Fuentes et al., 2016; van Wijk et al., 2018). However, this position is somewhat arbitrary, and the regional relationship between the two provinces remains poorly understood. The clear rift boundaries that are present in the central and northern segments of the RGR become more diffuse in the southern RGR where it lies adjacent to another extensional province.

Apatite fission-track (AFT), AHe, and ZHe thermochronology data have been previously reported from uplifted fault blocks along the RGR in Colorado and New Mexico. AHe data combined with syn-rift stratigraphy suggest a main phase of unroofing beginning in the Late Oligocene–Miocene for the northernmost RGR in Colorado (Landman and Flowers, 2013). In the Albuquerque Basin >45 Ma ages on inactive rift flanks indicate Laramide activity (Ricketts et al., 2016), but a 22–16 Ma period of initial uplift has been recorded for the Sandia Mountains, an active rift-flank (House et al., 2003; Ricketts et al., 2016). Accelerated sedimentation was observed after 16 Ma in the Sandia Mountains, as well as a period of isothermal relaxation following volcanic activity in the Oligocene (House et al., 2003). Laramide AFT dates were also used to identify the eastern and western extents of the RGR in the Sierra Ladrones, Blanca Peak, the Sawatch Range, and the Joyita Hills (Kelley et al., 1992). Abbey and Niemi (2018) documented RGR rifting in the Sawatch Range in central Colorado initiating at 22 Ma.

Complex patterns in exhumation and uplift similar to those revealed by AFT data in the northern and central segments of the RGR have also been observed in the southern RGR, suggesting that the entire RGR may have experienced extension at roughly the same time (Kelley and Chapin, 1997; Ricketts et al., 2016). AHe data from new areas as well as from the AFT samples of Kelley and Chapin (1997) were found to constrain a period of rapid exhumation that occurred ca. 25–10 Ma along the entire >1000 km length of the Rio Grande rift (Ricketts et al., 2016).

In many low-temperature thermochronology studies, a single thermochronometer is used on multiple samples collected along an elevation traverse in the footwall of a normal fault (e.g., Kelley et al., 1992; Kelley and Chapin, 1997; Foster and John, 1999; Stockli et al., 2002). However, in the Rio Grande rift–Basin and Range transition of southern New Mexico, our study area (Fig. 3), topographic relief in the footwalls of normal faults is generally limited to less than 1.5 km such that acquiring age-elevation transects is not possible. In regions with limited footwall topography, an approach using multiple thermochronometric dating techniques from a single sample can be used to constrain its time-temperature path as it cooled during extensional activity (e.g., Foster et al., 1993; Harrison et al., 1995; Axen et al., 2000; Stockli, 2005). Our approach was to collect one or more samples in fault-block uplifts adjacent to major normal faults for AHe, AFT, and ZHe analysis. These three techniques are complementary because they cover a significant temperature range and have overlapping temperature-sensitivity windows, ranging from 30 to 90 °C for AHe (Flowers et al., 2009), 60–110 °C for AFT (Gleadow et al., 1986), and 50–210 °C for ZHe (Guenthner et al., 2013; Johnson et al., 2017).

The majority of the samples are from Proterozoic felsic intrusive rocks from basement rocks beneath the major nonconformity in the region overlain by Cambrian–Ordovician sandstone. These samples include granitoids, amphibolite, and felsic orthogneiss. Other samples include Cambrian granite from the Florida Mountains, Cambrian–Ordovician sandstone from the San Andres Mountains (Bliss Formation; Amato and Mack, 2012), Cretaceous sandstone from the Potrillo Mountains, and Paleogene granitic rocks from the Cookes Range and Peloncillo Mountains. For the AHe data, 13 samples were Proterozoic rocks, two samples were Cambrian granite, and two samples were Paleogene intrusive rocks. For the ZHe data, 16 samples were Proterozoic, one sample was Cambrian granite, one sample was Cambrian sandstone, two samples were Cretaceous sandstone, and two samples were Paleogene intrusive rocks. The AFT data include analyses of 16 Proterozoic samples, three Cambrian samples, and two Paleogene samples (see summary in Table S1¹).

Apatite and zircon were separated using standard mineral separation techniques, and individual grains were selected under a polarizing microscope. Suitable grains for AHe and ZHe analysis are >70 mm measured perpendicular to the c-axis, are euhedral, and do not have visible inclusions. Crystal lengths, widths, and depths were measured to apply a Ft age correction (Farley et al., 1996). Individual apatite and zircon grains were packed into niobium tubes for analysis at the Colorado Thermochronology Research and Instrumentation Laboratory at the University of Colorado, Boulder. Apatite fission-track analyses for this study were carried out at New Mexico Institute of Mining and Technology following the methods outlined in Kelley et al. (1992) and the equations of Hurford and Green (1983). Please see Table S2 (^{footnote 1}) for details of AFT data used in HeFTy modeling.

Individual AHe and ZHe dates obtained from a sample typically do not correspond with any single geologic event, and multiple grains from a single sample sometimes result in some spread of dates (e.g., Flowers and Kelley, 2011). Flowers et al. (2009) noted that AHe dates are governed by differential radiation damage in the crystal, resulting in differential He diffusion. The radiation damage accumulation and annealing model (RDAAM) predicts a correlation between AHe date and effective uranium of a crystal (eU = U + 0.235Th), where eU is used as a proxy for radiation damage. In cases where a sample cools slowly through the AHe partial retention zone (PRZ), AHe dates will be positively correlated with eU, and in cases where a sample cools quickly through this temperature window, AHe dates will be consistent across a wide range of eU values (Flowers et al., 2009). Similarly, Guenthner et al. (2013) described the zircon radiation damage accumulation and annealing model (ZRDAAM), which incorporates a relationship between ZHe dates and radiation damage in the zircon crystal. Similar to the AHe method, ZHe dates either show a positive or negative trend with eU values, where individual ZHe dates are the cumulative result of the overall thermal history and effects of radiation damage on helium diffusion. This model is used to explain intra-sample ZHe dates from a single sample that sometimes span hundreds of millions of years. For these reasons, multiple single-grain AHe and ZHe dates were obtained from each sample to detect possible intra-sample date variability.

The AFT method is complementary to AHe and ZHe because differential annealing of fission tracks in apatite can be used to constrain a continuous thermal history in the AFT partial annealing zone (PAZ) (60–110 °C; Gleadow et al., 1986). For example, long track lengths (15–16 mm) are indicative of rapid cooling through the AFT PAZ, whereas shortened tracks or a bimodal track length distribution suggest prolonged residence in this temperature window or reheating (e.g., Donelick, 2005; Ketcham, 2005; Stockli, 2005).

The software program HeFTy (v. 1.9.3; Ketcham, 2005) was used to generate inverse models that use AHe, ZHe, and AFT data sets to generate probable time-temperature paths that match the observed data and uncertainty. Complete modeling inputs, assumptions, and constraints used for each sample are shown in Table S3 in the Supplemental Material (see ^{footnote 1}), following methods outlined in Flowers et al. (2015).

The data selection process of Flowers and Kelley (2011) was implemented to select grains for inverse modeling inputs. If single-grain AHe dates are older than AFT dates for a sample, and/or eU relationships cannot explain dispersion in AHe or ZHe dates, those grains were not used for inverse modeling. For grains selected for modeling, U, Th, and Sm concentrations for each grain are included, as well as the grain radius (the radius of a sphere of equal volume to that of the crystal) (Ketcham, 2005). The radiation damage accumulation and annealing models of Flowers et al. (2009) and Guenthner et al. (2013) were used as helium diffusion models for AHe and ZHe grains, respectively.

The Chiricahua Mountains are a southwest-tilted block in southeast Arizona bounded to the east by a normal fault (Fig. 3). This buried fault separates the Chiricahua Mountains from the downthrown San Simon Valley block. In the late Cretaceous, the range was folded and faulted during Laramide shortening (du Bray et al., 1997). Major igneous activity occurred in the Chiricahua Mountains from 29 to 23 Ma, with volcanic material from the 29 Ma Turkey Creek caldera dominating the central part of the range. Normal faulting coeval with Cenozoic volcanism and plutonism is responsible for the modern topography and position of the Chiricahua Mountains block (Drewes et al., 1995).

A total of seven AHe dates (Table 1), six ZHe dates (Table 2), and two AFT dates (Table 3) were obtained from two samples collected in the Chiricahua Mountains (Biddle et al., 2018). AHe dates range from 14 to 21 Ma across an eU range of 40–183 ppm (Fig. 4A). ZHe dates are also consistent, ranging from 24 to 29 Ma with corresponding eU values ranging from 120 to 265 ppm. Both AHe and ZHe dates show nearly flat correlations with eU. AFT dates of these samples are 25.6 ± 1.7 Ma and 23.9 ± 1.7 Ma with mean track lengths of 13.2 ± 0.2 µm (n = 49) and 10.8 ± 0.4 µm (n = 15), respectively. These dates are older than accompanying AHe dates from the same samples and younger or similar to the range of ZHe dates.

The Peloncillo Mountains in southwestern New Mexico (Fig. 3) were uplifted and deformed in the mid to late Tertiary. This range contains a record of Late Cretaceous to early Tertiary andesitic volcanism but no known thrust faulting (Gillerman, 1958). Normal faults bound the mountains to the west and east, in the San Simon and Animas Valleys, respectively. The central part of the range consists primarily of southwest-tilted fault blocks. The central range includes Granite Gap, part of an east-west–trending, up-thrown block bounded by south-dipping normal faults. The relationship between the Granite Gap faults and the other faults in the range is unclear, but Granite Gap faults appear to offset the north-trending blocks.

A single sample was collected from the Peloncillo Mountains (Fig. 3). This sample yielded a total of two AHe dates, three ZHe dates, and one AFT date (Fig. 4B). AHe dates are 25 Ma and 46 Ma, with corresponding low eU values of 20 and 16 ppm, respectively (Table 1). ZHe dates range from 24 to 26 Ma and have higher eU values ranging from 1058 to 1236 ppm (Table 2). The AFT date of 33.8 ± 7.4 Ma (n = 7) (Table 3) is slightly older than ZHe dates but overlaps within the large error in the AFT date related to low track counts.

The geology of the Little Hatchet Mountains in southwestern New Mexico (Fig. 3) is dominated by Laramide thrust faulting and folding of Cretaceous, Jurassic, and minor Proterozoic rocks (Clinkscales and Lawton; 2015, 2017). Late Cretaceous volcanism associated with the Laramide orogeny is present throughout the Little Hatchet Mountains, and a second pulse of igneous activity in the Eocene resulted in a granitic intrusion in the southern part of the range (Channell et al., 2000). The modern topography of the Little Hatchet Mountains is the product of Cenozoic extension, where boundary faults strike north along the east and west sides of the range (Zeller, 1970; Clinkscales and Lawton; 2015).

AHe and ZHe data were obtained from three samples collected from the Little Hatchet Mountains (Biddle et al., 2018). However, these samples yielded poor-quality apatites, and only a single AHe date (23 Ma) was obtained from sample 17LH01. A total of nine ZHe dates were obtained from these three samples, ranging from 21 to 151 Ma across a considerable eU spread from 167 to 734 ppm (Fig. 4B).

The Burro Mountains in southwestern New Mexico (Fig. 3), considered part of the Basin and Range (Woodward, 1987), consist of the Big Burro Mountains to the north and the Little Burro Mountains to the south. Range-bounding faults strike NNW along the flanks of these uplifts (Gillerman, 1970). An unconformity exists between Proterozoic basement rock (ca. 1680–1220 Ma; Amato et al., 2008, 2011) and the mid-Cretaceous Beartooth quartzite (Lawton et al., 2020). This unit is overlain by Paleogene and Neogene volcanic and sedimentary rocks, as uplift and subsequent removal of Paleozoic rocks occurred during the Laramide orogeny (Hedlund, 1980; Ross and Ross, 1986; Kues and Giles, 2004).

A total of five samples were collected from the Burro Mountains. AHe dates from these samples range from 9 to 33 Ma across an eU spread from 4 to 204 ppm (Fig. 4C). ZHe data were only obtained from sample 10BM-196. These dates range from 23 to 41 Ma with eU values ranging from 312 to 1297 ppm. Two AFT dates of 26.7 ± 1.4 Ma (10BM204) and 17.3 ± 1.9 Ma (10BM219) with mean track lengths of 13.7 ± 0.2 mm (n = 100) and 13.2 ± 0.2 mm (n = 100), respectively, are very similar to corresponding AHe dates. One AFT date for 05BM183 is younger at 11.0 ± 3.8 Ma; the large error is due to low track density.

The Cookes Range (Fig. 3) is bounded by normal faults on three sides and exposes Paleozoic sedimentary rock units in the northern extent of the range. The majority of Proterozoic rock exposure is along Fluorite Ridge in the southern part of the range (Clemons, 1982). Cookes Peak, the high point of the range, is a ca. 38 Ma granodiorite pluton (McLemore et al., 2001) associated with latest Laramide subduction magmatism. Biddle et al. (2018) attempted (U-Th)/He thermochronometry on two samples from the Proterozoic rocks of the Cookes Range but were unable to generate a conclusive thermal model because of unexplained spread in AHe and ZHe ages.

AHe data were obtained from three samples from the Cookes Range, ZHe data from one sample, and AFT data from one sample (Fig. 4D). AHe dates from all three samples, including the 38 Ma Cookes Peak pluton and Proterozoic granite, range from 8 to 20 Ma across an eU range of 7–73 ppm (Biddle et al., 2018, and this study). ZHe dates collected from sample 13COOKS-2 range from 28 to 38 Ma across a limited eU span of 102–186 ppm, and the AFT date for this sample is 25.8 ± 3.0 Ma with a mean track length of 13.7 ± 0.2 (n = 76). Sample 12CR-04 yielded an AFT age of 10.1 ± 1.3 Ma and a mean track length of 12.7 ± 0.2 (n = 100). In contrast, AHe and AFT data from the Cookes Range are notably younger than previously reported ZHe dates, which range in age from 44 to 446 Ma and are largely older than Basin and Range and Rio Grande rift extension (Biddle et al., 2018; Reade et al., 2020).

The Florida Mountains (Fig. 3) form an east-tilted fault block with range-bounding faults that are mostly concealed by the expansive bajada surrounding the range (Clemons, 1998). The Florida Mountains fault is the major range-bounding fault and was likely active from ca. 29 Ma to Pleistocene time, producing ∼2000 m of vertical displacement and tilting the range 23° to the northeast (Clemons, 1998). Most of the crystalline rock in the range is Cambrian in age, but limited Proterozoic granite exists in the west-central region and is dated at ca. 1.63 Ga (Amato and Mack, 2012). Approximately 1300 m of Paleozoic section is exposed in the Florida Mountains and is overlain unconformably by Paleogene Lobo Formation (Clemons, 1998).

AHe data are presented from two samples, ZHe data from one sample, and AFT data from four samples collected in the Florida Mountains (Biddle et al., 2018, and this study). AHe dates range from 10 to 22 Ma with 5 to 36 ppm eU values (Fig. 4D). ZHe dates are similar (23–28 Ma), yet at much higher eU values of 537–737 ppm. AHe dates are 16.5 ± 2.4 Ma age with no track lengths, 28.6 ± 2.3 Ma with mean track length of 13.0 ± 0.2 mm (n = 39), 30 ± 2.1 Ma with mean track length of 13.4 ± 0.5 mm (n = 24), and 12.2 ± 4.4 Ma with no track lengths.

The Caballo and Fra Cristobal Mountains are two north-striking, east-dipping fault blocks separated by the Cutter Sag accommodation zone (Seager and Mack, 2003; Fig. 3). Both ranges expose Proterozoic granites, gneisses, and amphibolite overlain by 2–4 km of Paleozoic sedimentary strata. The Caballo–Hot Springs fault forms the major range-bounding structure in the Caballo Mountains and has accommodated an estimated ∼6 km of stratigraphic separation (Seager and Mack, 2003). The Fra Cristobal Mountains have similar stratigraphy, and their Proterozoic basement rock has been dated at 1.45 Ga (U-Pb; Amato et al., 2018). The Hot Springs fault borders the range on the west and has a minimum throw of 1100 m (Seager and Mack, 2003). Neither range has been the subject of any (U-Th)/He study, but AFT ages of 5–16 Ma have been reported from the Caballo Mountains (Kelley and Chapin, 1997). The range of AFT ages, which were collected from different points along a N–S transect, are interpreted to have been caused by progressive growth and/or propagation of the range-bounding fault. The Mud Springs Mountains are a northeast-tilted fault block north of Truth or Consequences, New Mexico, and Mud Mountain, where our samples were collected, is a spur off of the main ridge of the Mud Springs Mountains. Mud Mountain is cored by Proterozoic gneiss, meta-sandstone, quartzite, and granite gneiss (Maxwell and Oakman, 1990). A single AFT date of 8.4 ± 3.8 Ma indicates Miocene cooling (Kelley and Chapin, 1997).

Four samples were collected from the Caballo Mountains, three from the Fra Cristobal Mountains, and two from nearby Mud Mountain (Fig. 3), resulting in 34 single-grain AHe dates, 16 single-grain ZHe dates, and seven AFT dates (Fig. 4E). In the Caballo Mountains, AHe dates range from 10 to 36 Ma across eU values from 2 to 122 ppm. ZHe dates from three samples range from 7 to 271 Ma. Collectively, these data show a negative relationship with eU values, which range from 475 to 2708 ppm (Fig. 4E). A single AFT date of 10.6 ± 2.1 Ma with mean track length of 14.4 ± 0.4 mm is reported from the Caballo Mountains. This date is similar to previously reported AFT dates from the Caballo Mountains that range from 5 to 11 Ma (Kelley and Chapin, 1997). Other AFT dates determined during this study are older at 14.6 ± 3.9 Ma and 17.3 ± 1.8 Ma.

AHe dates from the Fra Cristobal Range vary from 8 to 13 Ma and eU values range from 11 to 139 ppm eU (Fig. 4E). ZHe dates are nearly identical (8–13 Ma) but at higher eU values of 657–1809 ppm. An AFT date from sample 15FC1 yielded a central date of 14.9 ± 1.8 Ma, with a mean track length of 13.0 ± 0.6 mm (n = 23). Two other AFT dates are 9.7 ± 1.2 Ma and 16.7 ± 1.2 Ma, with mean track lengths of 13.3 ± 0.3 µm (n = 66) and 13.1 ± 0.3 µm (n = 72), respectively.

Two samples from Mud Mountain yield AHe dates that collectively range from 9 to 18 Ma and ZHe dates that range from 26 to 288 Ma (Fig. 4F). AHe eU values range from 2 to 40 ppm, and ZHe eU values range from 388 to 851 ppm. One sample from Mud Mountain gave an AFT date of 10.3 ± 1.1 Ma and mean track length of 12.4 ± 0.5 µm (n = 33), within error of an 8.4 ± 3.8 Ma date measured previously by Kelley and Chapin (1997).

The San Andres Mountains are a 240-km-long west-tilted fault block bounded on the west by the Jornada del Muerto Basin and on the east by the Tularosa basin (Fig. 3; Seager, 1981). Proterozoic granites and schists of the San Andres Mountains (Howland, 2018) are overlain by ∼2600 m of Paleozoic strata (Seager, 1981). Existing AFT dates from the San Andres Mountains vary between 7 and 22 Ma and suggest a complex Cenozoic cooling history (Kelley and Chapin, 1997). Ricketts et al. (2016) obtained AHe dates for a Proterozoic sample that ranged from 17 to 9 Ma, younger than a published AFT age of 22 Ma of the same sample (Kelley and Chapin, 1997).

A total of 30 AHe dates, 19 ZHe dates, and six AFT dates were collected from ten samples in the San Andres Mountains from locations in the range not previously studied (Fig. 4G). Collectively, AHe dates range from 3 to 36 Ma, and ZHe dates show a much wider spread from 2 to 190 Ma. Six samples from the San Andres Mountains produced AFT dates between 10 and 14 Ma with mean track lengths between 12.8 and 14.1 mm (only one had n>31). These are consistent with previous AFT dates recorded in the southern San Andres Mountains by Kelley and Chapin (1997).

The East Potrillo Mountains in south-central New Mexico (Fig. 3) are a horst block bounded by the East Potrillo fault to the northeast and the Mount Riley fault zone to the southwest (Seager and Mack, 1994). The range contains a system of late Oligocene–early Miocene low-angle normal faults that crosscut late Cretaceous–early Tertiary compressional structures (Seager and Mack, 1994). These structures are Laramide thrust sheets and associated drag folds that increase in deformational severity toward the north of the range. Mid-Paleogene extension in the East Potrillo Mountains has been estimated to be 75%–100%, similar to estimates of extension across the central Rio Grande rift, indicating that the Potrillo Mountains were uplifted during this later period of extensional tectonism.

No suitable apatite was obtained from the East Potrillo Mountains for AHe dating, and a single sample was collected for ZHe analysis (Fig. 4G). Sample 16EP01 yielded three ZHe dates that range from 21 to 46 Ma from eU values ranging from 135 to 234 ppm. This sample also contained apatite shards that yielded an AFT date of 17.5 ± 2.3 Ma with a mean track length of 12.7 ± 0.7 mm.

The Franklin Mountains are a north-trending, west-tilted block in El Paso, Texas, east of the Mesilla Basin in the Rio Grande rift (Fig. 3). The mountains record both thin- and thick-skinned Laramide deformation and represent a transition between the two deformation styles in the Chihuahua trough (Carciumaru and Ortega, 2008). Rio Grande rift–associated uplift of the Franklin Mountains is thought to have occurred in the Miocene (Harbour, 1972). Quaternary uplift is evidenced by bowed and deformed terraces that flank the east and west sides of the Franklin Mountains (Armour et al., 2018), and based on fault scarp morphology, the latest movement along the East Franklin boundary fault occurred in the Pleistocene or early Holocene (Machette, 1987).

A total of 26 AHe dates were collected from five samples from the Franklin Mountains in western Texas (Fig. 4G). Collectively, AHe dates range from 8 to 40 Ma with very restricted eU values of 4–43 ppm. Previously reported ZHe data from the Franklin Mountains show a very large spread in ZHe dates, from 19 to 401 Ma and eU values from 63 to 828 ppm (Biddle et al., 2018; Reade et al., 2020). One sample yielded an AFT date of 17.7 ± 3.5 Ma; the large error is associated with low track density.

Valid interpretation of AHe and ZHe dates and the resulting thermal history of a sample sometimes requires knowledge of the spatial distribution of U and Th within the crystal. U and Th zonation (or lack of zonation) affects He diffusivity in the crystal by modifying the spatial distribution of radiogenic He (e.g., Farley, 2000; Meesters and Dunai; 2002; Hourigan et al., 2005). When combining multiple low-temperature thermochronometers, crystal zonation can sometimes explain seemingly incompatible data sets, such as AHe dates that are older than corresponding AFT dates (e.g., Flowers and Kelley, 2011). The lack of spread in AHe and ZHe dates in this study and the observation that AHe and ZHe dates are similar to corresponding AFT dates for many of the samples in this study suggest that crystal zonation most likely does not affect our interpretation of the thermal history.

However, etched fission tracks from some of the AFT samples show U zonation (Fig. 5), which may influence interpretation of the thermal models. To test the possible influence of crystal zonation, we include several example inverse models using AHe, ZHe, and AFT data from sample 16CH01 from the Chiricahua Mountains (Fig. 5). These end-member models include no zonation, U- and Th-rich cores, and U- and Th-rich rims. The final inverse models for each scenario are similar, and a summary of all the resulting best-fit and weighted-mean paths in the bottom plot are nearly identical. This simple analysis suggests that possible zoning in each apatite and zircon crystal most likely does not affect our final interpretations from the inverse models. Therefore, the inverse models presented below assume uniform distribution of U and Th. They are presented in two geographic areas outlined on Figure 3, and the average ages from each sample are plotted on maps showing the Basin and Range (Fig. 6A) and the Rio Grande rift (Fig. 6B).

For samples that yield AHe, AFT, and ZHe dates that are all younger than 100 Ma, models begin at 100 Ma and 280–300 °C, a temperature higher than the ZHe PRZ (Figs. 7–9). Some samples show ZHe dates older than 100 Ma, and ZHe dates for these samples show well-defined negative correlations with eU. For these samples, our models begin at the time of rock crystallization during the Proterozoic to explore the effect of long-term radiation damage on the resulting thermal history. In addition, we present two inverse models for samples from the Basin and Range (Fig. 7), including the transition zone areas of the Cookes Range and Florida Mountains (Fig. 8), to test for possible effects of reheating during Oligocene magmatism. One model includes monotonic cooling to surface temperatures, and the other allows for, but does not require, a pulse of reheating during 36–23 Ma volcanism associated with the Boot Heel volcanic field (Chapin et al., 2004). All inverse models are forced to surface temperatures at present-day conditions (15–25 °C). Ideally, models only incorporate data from a single sample, but in several instances, multiple samples were combined to provide enough data for a robust thermal model.

Samples 16CH01 and 16CH02 were modeled using AHe, AFT, and ZHe data (Fig. 7). Both models record cooling through the ZHe PRZ from ca. 30–20 Ma. Sample 16CH02 suggests the sample was cooled to near-surface temperatures by 20 Ma, whereas sample 16CH01 records a second pulse of rapid cooling beginning 5 Ma and continuing to the present. The initial cooling in both samples is contemporaneous with two regional geologic events. The Desert Southwest core complex extensional episode of Dickinson (2002) occurred from 30 to 12 Ma in southern Arizona and easternmost California. More locally, development of the Turkey Creek caldera within the Chiricahua Mountains occurred at 29 Ma (du Bray et al., 1997).

Apatite fission-track and ZHe dates for the Chiricahua Mountains both overlap in time with the eruption of the caldera, suggesting that this igneous system was most likely emplaced at depths shallow enough to reset both ZHe and AFT systems. If so, the main phase of extension to bring these rocks to near-surface levels occurred prior to eruption of the Turkey Creek caldera, and AHe, AFT, and ZHe data reflect relaxation of isotherms.

For sample 16PE01, we present one inverse model that includes AHe and ZHe data (Fig. 7). Thermochronologic data from this sample are all similar to the crystallization age of 33.2 ± 0.2 Ma (McLemore et al., 1995), and the resulting inverse model suggests rapid cooling to <50 °C by 30 Ma.

Two samples from the Little Hatchet Mountains (17LH01 and 17LH02) were combined to produce a single inverse model (Fig. 7). These two samples were collected relatively close together, and are not separated by any major faults, so should have very similar thermal histories. The final inverse model includes AHe and ZHe data and indicates rapid cooling beginning at 30–25 Ma and reaching near-surface temperatures <50 °C by 27–20 Ma.

Four samples from the Burro Mountains were modeled using a combination of AHe, AFT, and ZHe data (Fig. 7). Sample 10BM-196 was modeled using a combination of AHe and ZHe grains and recorded a period of rapid cooling (>20 °C/m.y.) beginning at 25 ± 5 Ma and reaching near-surface temperatures (<40 °C) by 15 ± 5 Ma. Sample 10BM-204 records a similar onset of cooling at 30 ± 5 Ma, but the model suggests slower cooling to near-surface temperatures between 25 Ma and 15 Ma. Sample SGBM-147 only produced grains for AHe, and the thermal history model produced records cooling below the AHe PRZ (>50 °C) by 20 ± 5 Ma. The thermal model alone does not provide enough information to suggest a cooling rate, but age-eU information suggests that sample cooled through the AHe PRZ relatively rapidly. Sample 05BM-183 also only used AHe grains and recorded a similar time of cooling.

Two calderas lie within an 80 km radius of the Burro Mountains: Twin Sisters and Schoolhouse Mountain (31.4 Ma and 34.9 Ma, respectively; Fig. 3; Chapin et al., 2004; Swenton, 2017). ZHe ages from the Burro Mountains are limited to two data points, but these can be used to infer the small effect of nearby magmatic activity on overall cooling history. The two grains are older (41 Ma) and younger (21 Ma) than both caldera systems, but the older age overlaps with other Mogollon-Datil volcanism (e.g., Chapin et al., 2004). Two AFT samples have mean track lengths of 13.2 µm and 13.7 µm, suggestive of a relatively short residency time for samples in the AFT PAZ. AFT ages from this area are also significantly younger than any magmatic event in the region, thus lending additional confidence to a limited reheating interpretation for this range. Despite inconclusive evidence for the driving force of rapid cooling observed in the Burro Mountains, it can be inferred based on geologic evidence that the basement had already been exhumed to shallow depths in the crust before the onset of Mogollon-Datil volcanism.

The Cookes Range and Florida Mountains lie in an area near the postulated boundary between the RGR and B&R (Fig. 3), and thus we have separated them into a “transition zone” region. The Proterozoic sample from the Cookes Range (12CR-04) (Fig. 8) was modeled with AHe and AFT data and records a period of rapid cooling below the AFT PAZ (∼100 °C) beginning at 15 ± 5 Ma. Sample 13COOKS-2 is from the ca. 38 Ma Cookes Peak pluton (Amato et al., 2018) and records a similar period of young, rapid cooling from ca. 15–10. AFT age and track length data from 12CR-04 (Proterozoic granite) suggest prolonged residency time in the AFT PAZ (3–5 km depth). However, this sample records a similar time of rapid cooling from 15 to 10 Ma, suggesting that the entire Cookes Range was cooled at this time. Sample 13COOKS2 has an AFT age of 26 ± 3 Ma with an average track length >13 µm, suggesting rapid cooling. Most likely, the rapid cooling observed in the 13COOKS2 sample is the result of isothermal relaxation following the cooling of the Cookes Peak pluton.

Of the four new AFT dates from the Florida Mountains, two are associated with samples that have AHe data (Biddle et al., 2018) and were used to improve existing inverse models. Sample 16FL-04 was modeled with AHe and AFT data and suggests a single period of cooling from ca. 18–12 Ma (Fig. 8), similar to the cooling history present for the Cookes Range. Sample 16FL-05, modeled with AHe, ZHe, and AFT data, is better constrained and records one period of rapid cooling beginning ca. 25 Ma and reaching surface temperatures by ca. 20 Ma.

Six samples from the San Andres Mountains were suitable for thermal history modeling (Fig. 9), and they record similar cooling histories. Modeling of sample 15SA-88 suggests an onset of relatively rapid cooling beginning at 15 ± 5 Ma. Similarly, models for samples 17SAS-40, 18SAS-44, and 18SAS-51 record periods of rapid cooling onset at 17 ± 5 Ma. Rapid cooling (>20 °C/m.y.) in each sample continued until ca. 10 Ma; in the context of their locations (see Fig. 3), this suggests a consistent uplift history for the entire San Andres fault block.

Thermal history modeling of three samples (06CM-01, 06CM-02, and 06CM-04) from the Caballo Mountains suggests rapid cooling beginning ca. 17 ± 2 Ma (Fig. 9). However, sample 06CM-02 may also show an earlier, less defined pulse of cooling from ca. 35–30 Ma, which may be the result of heating from nearby Paleogene dikes. All of the thermal history models from the Fra Cristobal range yield similar results that suggest rapid cooling from above 300 °C beginning at 15 ± 2 Ma and reaching surface temperatures by ca. 12 Ma.

Thermal history models produced from our AHe, AFT, and ZHe data constrain times of cooling across the southern Basin and Range–Rio Grande rift transition. In order to more carefully compare thermal models from different locations, times of rapid cooling were estimated for each model. Comparison between monotonic cooling models and reheating models from the Basin and Range samples (Fig. 7) shows very similar times of cooling to near-surface temperatures. In the reheating models, samples reach high enough temperatures to almost completely remove older radiation damage and result in final cooling that agrees with the monotonic models. Most monotonic cooling thermal models produced in this study show sharp inflection points at the transition from slow cooling to rapid cooling (Figs. 7 and 8). These inflections are used to estimate the onset of cooling for each sample. Similarly, most models preserve a period of rapid cooling and a second inflection where cooling dramatically slows at near-surface temperatures. These two inflections bracket the time period of rapid cooling for each sample. For each sample, we plot the onset of cooling, the end of cooling, and the average of the two, and use this to compare and contrast times of cooling along an east-west transect across eastern Arizona and southern New Mexico (Fig. 10). We also include thermal models from the southern Rio Grande rift from Ricketts et al. (2016) in our discussion.

Collectively, thermal models for Basin and Range samples indicate an onset of cooling that ranges from 34.8 to 24.6 Ma (Fig. 10). Initial cooling older than 30 Ma is observed in the Chiricahua Mountains, the Peloncillo Mountains, and the Burro Mountains, whereas the Burro Mountains and Little Hatchet Mountains also show slightly younger cooling. The total time of rapid cooling observed in each sample varied, ranging from ca. 5–20 Ma. The earliest end of cooling occurred in the Peloncillo Mountains at 29 Ma, whereas all other samples indicate an end of cooling from 20.4 to 13.8 Ma (Fig. 10).

The Cookes Range and Florida Mountains are located in an area that has been described as the transition zone between the Basin and Range and RGR (Fig. 3; Mack, 2004). No thermochronological data existed for either range prior to this work, and so previous work relied on fault relationships to constrain times of uplift (Clemons, 1982; Clemons and Mack, 1988). The Cookes Range and Florida Mountains record slightly different cooling histories despite being roughly along strike of one another, in contrast to other along-strike ranges (Caballo, Fra Cristobal, and San Andres) with more uniform cooling histories. Observed onset of cooling in one sample from the Florida Mountains occurred at 28.8 Ma, similar to the onset of cooling in the southeastern Basin and Range (Fig. 10). All three other samples, however, indicate a younger onset of cooling beginning 17.7–10.9 Ma. The end of extension in the Cookes Range was from 7.5 to 5.1 Ma, whereas end of cooling was slightly earlier in the Florida Mountains at 16.1–11.5 Ma.

Thermal models from the southern Rio Grande rift from this study are combined with additional thermal models from the Organ Mountains and Sierra Blanca from Ricketts et al. (2016). Together, these models indicate initial cooling that ranges from 24.9 to 11.7 Ma (Fig. 10). End of cooling in the southern Rio Grande rift was from 13.6 to 6.7 Ma, which is younger than any samples from the southeastern Basin and Range.

Kernel density estimation (KDE) plots (Vermeesch, 2012) of individual AHe and ZHe dates from each of the two areas (B&R and RGR) were created to evaluate the broad differences in cooling between the areas (e.g., Carrapa et al., 2019). These plots omit details and issues related to the age versus eU relationship, but they can be helpful to see the general trends when large numbers of samples are plotted together, despite the realization that the cooling ages are individually controlled by numerous variables (e.g., Flowers et al., 2015). The higher-temperature (<210 °C) cooling dates revealed by ZHe analyses show a pronounced older range of dates in the B&R, with an age peak ca. 27 Ma, whereas in the RGR, the peak age is at 12 Ma (Fig. 11A). The AHe KDE curves are more similar in the two areas, with the B&R peak at 16 Ma and the RGR peak at 12 Ma (Fig. 11B). AFT ages (Fig. 11C) have more scatter because of low n values, but also indicate older cooling in the B&R relative to the RGR. These data are consistent with the hypothesis that the RGR represents a distinctly younger phase of extension than the southeastern B&R.

Previous workers have noted that in some samples, AHe ages are older than AFT ages, despite evidence that AHe has a lower closure temperature (e.g., Flowers and Kelley, 2011; Green and Duddy, 2018). In a comparison of our AFT versus AHe dates (Fig. 12), most samples lie either along the 1:1 line or below it, indicating that for most samples, data dispersion and “inverted” apatite dates are not observed, and these samples did not likely reside within the AHe PRZ for prolonged periods of time (Flowers and Kelley, 2011).

The observed thermochronological dates obtained in this study and cooling histories interpreted from inverse models can be attributed to multiple processes, including reheating during magmatism and/or volcanism, monotonic exhumation related to faulting, or a combination of the two. Inverse models from the Basin and Range samples are compatible with both scenarios. While sample collection was designed to target footwalls of normal faults to document times of extensional exhumation, multiple sample sites in the southeastern Basin and Range preserve evidence of magmatism that likely reset thermochronological dates. ZHe dates ranging between 29.5 and 24.4 Ma from the Chiricahua Mountains overlap with ages of Oligocene volcanic rocks in the region, including the 29–25 Ma Turkey Creek caldera (du Bray et al., 1997), and similar ZHe dates, despite a wide range in eU concentration, suggest fast cooling through the ZHe PRZ (Fig. 4). Volcanism occurred in the Schoolhouse Mountain caldera, which is adjacent to our Burro Mountains samples (Fig. 3) from ca. 35–33 Ma (Swenton, 2017), and volcanism throughout the region occurred from 36 to 23 Ma (McIntosh et al., 1992; McIntosh and Bryan, 2000; Chapin et al., 2004). We compiled the volcanic rock ages in the region (from the North American Volcanic and Intrusive Rock Database [NAVDAT]; Carlson et al., 2001) included within Figure 3, splitting them into B&R versus RGR regions with a division at 107.5 °W, which coincides with the RGR boundary of Mack (2004). The resulting B&R compilation includes 317 entries from rocks with ages ranging from 40 to 0 Ma (Fig. 10). A KDE (Vermeesch, 2012) of the data shows that the vast majority of rock ages are older than 25 Ma, including multiple peaks at 37, 33, 30, and 27 Ma (Fig. 10). These peaks overlap with the times of cooling in the Chiricahua Mountains, Peloncillo Mountains, and with cooling observed in all but one of the Burro Mountains samples.

These observations suggest it is likely that thermochronological ages from the southeastern Basin and Range were heavily influenced by magmatism. Oligocene calderas erupting within the Mogollon-Datil and Boot Heel volcanic fields (Fig. 3) blanketed much of the region in ignimbrite outflow sheets (McIntosh and Bryan, 2000). In addition, modeling of the influence of post-magmatic cooling on low-temperature thermochronometers suggests that plutons up to 10–15 km depth can reset low-temperature dates in the upper crust, depending on the size and the emplacement temperature of the pluton (Murray et al., 2018). However, Murray et al. (2018) also note that resetting of different thermochronometric systems during magmatism would result in different dates for each system. These relationships indicate that samples collected from the southeastern Basin and Range were at near-surface temperatures by the end of this intense period of magmatism, and AHe, AFT, and ZHe dates were at least partially reset during this event. While these data do not directly capture the time of extension in this region, they are compatible with three possible models for exhumation: (1) Exhumation preceded volcanism, and Proterozoic rocks had already been exhumed to the surface by the time igneous rocks were emplaced in the crust; (2) exhumation was synchronous with magmatism such that the cooling observed in inverse models is driven by a combination of post-relaxation cooling of the crust and exhumational cooling during extension; or (3) exhumation initiated after volcanism ceased and was of limited magnitude such that thermochronological dates do not record unroofing during this event.

Although disentangling the effects of magmatic reheating versus extensional exhumation is difficult, comparison with other regions in the rift offer some insight. Total extension in the rift increases to the south, from 8%–12% in the San Luis Basin in Colorado (Kluth and Schaftenaar, 1994) to 50%–100% in southern New Mexico (Morgan et al., 1986). We note that times of rapid cooling in the low-extension Sawatch Range (north of the San Luis Basin) and the high-extension Organ Mountains (southern New Mexico) both are distinctly younger than Oligocene magmatism (Ricketts et al., 2016). Since total extension in southern New Mexico is much higher than in the northern Rio Grande rift, then by analogy if extensional exhumation had initiated after the main phase of volcanism in the Boot Heel volcanic field, it would have imparted a distinctly younger pulse of cooling observed in the thermal models. For these reasons, we discount this as a possibility, and instead propose that extension occurred either prior to or synchronous with magmatism.

Exposures of Miocene and possibly older basin fill in the southeastern Basin and Range are limited such that the early history of extension in this region is not well constrained (Mack, 2004), but available data are consistent with extensional cooling older than 25 Ma inferred from thermochronological data. In the Boot Heel region of southwestern New Mexico, thick accumulations (250–1000 m) of rift fill suggest that much of it may be Miocene in age (Thompson et al., 1978; Heywood, 1992; Klein, 1995; Chang et al., 1999; Hawley et al., 2000; Shearer and Miller, 2000). At other locations, the presence of Miocene or possible older extensional basins is suggested by thick accumulations of strata in the San Vicente and Mimbres basins (Seager, 1995). Isolated, scattered outcrops of presumably Miocene or older sediments are also located along the western edge of the Playas basin (e.g., Elston et al., 1983) and southwest of the Burro Mountains (Hawley et al., 2000; Hibbs et al., 2000). These observations suggest that the cooling observed in the thermal models is likely a combination of thermal relaxation following Boot Heel magmatism and initial extensional exhumation. This model would amplify the rate of cooling, resulting in more uniform thermochronological dates from different systems.

Similarly, volcanic rock ages are compiled from the southeastern Rio Grande rift from the area included in Figure 3, but east of 107.5 °W longitude (Fig. 10), where widespread ignimbrite volcanism did not play a major role other than in the Organ Mountains (Fig. 3). The igneous rock compilation includes a total of 177 individual rock ages, and the resulting KDE shows only two main peaks at 36 Ma and <1 Ma, though small-volume RGR rift basalts are present with ages of ca. 28 Ma (Clemons, 1979), 10 Ma (Seager et al., 1984), ca. 8–5 Ma (Richard and Amato, 2018), and ca. 2 Ma (Amato et al., 2012). The younger peak, which is absent in the Basin and Range KDE, reflects the presence of Quaternary volcanic rocks.

There are two notable characteristics of thermochronological dates from the southern Rio Grande rift that differ from the southeastern Basin and Range data set. First, many locations contain ZHe dates from this study that are much older than the age of the rift, including the Caballo Mountains, Mud Mountain, San Andres Mountains, and East Potrillo Mountains. This was explored more thoroughly for samples in the Cookes Range, New Mexico, and Franklin and Carrizo Mountains of western Texas, where single-sample ZHe data sets span hundreds of millions of years and were used to reconstruct billion-year thermal histories for these areas (Reade et al., 2020). The presence of older ZHe dates requires that these grains were not completely reset during Oligocene magmatism. Second, many AHe, AFT, and ZHe dates from this region do not coincide with main peaks in the volcanic rock KDE (Fig. 10). This is also reflected in the resulting thermal models, where times of cooling in the southern Rio Grande rift also do not coincide with main periods of volcanism. Rather, cooling occurred during a time of volcanic quiescence and follows the 36 Ma peak by 10–20 m.y. This relationship has also been observed in other thermochronology studies of the Rio Grande rift (Ricketts et al., 2016; Abbey and Niemi, 2020) and suggests that cooling observed in all of the southern Rio Grande rift thermal models is driven by extensional faulting rather than magmatism.

The most commonly cited evidence for the onset of rifting in the southern RGR comes from a study of an inferred half-graben in the Cedar Hills Vent Zone, a linear belt of ca. 35 Ma rhyolite domes and flows concentrated along incipient faults west of Las Cruces, New Mexico (Seager and Clemons, 1975; Mack, 2004). Seager (1973) interpreted the Goodsight–Cedar Hills basin as being related to volcanic subsidence, but Mack et al. (1994a) postulated that because the tuffs were distally derived, the feature was unlikely to represent subsidence related to local eruptions. The Goodsight Mountains lie in between two known calderas in the Mogollon-Datil volcanic field (Fig. 3), the Emory Caldera (ca. 35 Ma; McIntosh et al., 1992) and the Organ Mountains Caldera (ca. 36 Ma; Zimmerer and McIntosh, 2013; Rioux et al., 2016). In addition, half-grabens are common features above collapsing calderas (e.g., Moore and Kokelaar, 1998), particularly in the early stages of caldera collapse, when subsidence is incremental and collapse occurs in several stages instead of catastrophically (e.g., Walter and Troll, 2001). We suggest that these features could reflect early subsidence over a magma chamber and that these observations do not record strong evidence for the initiation of the Rio Grande rift as is commonly cited (e.g., van Wijk et al., 2018).

The earliest voluminous rift-related mafic magmatism comes from the ca. 26–27 Ma Uvas basalt (see age summary in Clemons, 1979), though small-volume mafic rocks are locally present in the region (e.g., McMillan et al., 2000; Vermillion et al., 2020). The Uvas basalt was inferred, based on geochemistry, to have resulted from lithospheric thinning and asthenospheric upwelling (McMillan et al., 2000) and thus represents strong evidence for extension in the rift at this time. We agree with Mack (2004) that the rapid accumulation of the Thurman Formation, which is interfingered with the Uvas basalts and has detrital sanidine ages as young as ca. 27 Ma (Boryta, 1994), could be additional evidence of rapid subsidence related to extension in the RGR at this time. By the Miocene epoch, several major uplifts and basins had begun to develop. Based on the age of sedimentary fill in the hanging wall of normal faults bordering the Caballo Mountains, this range began to form during the late Oligocene through slip along border faults (Mack et al., 1994b), although most of the sedimentary rocks there are most likely early Miocene in age (Seager et al., 1984). In the Winston graben along the western edge of the southern Rio Grande rift, basin-fill strata were deposited at least from 18 to 4.8 Ma based on the ages of two interbedded volcanic flows (Seager et al., 1984; Harrison, 1994). At other locations, subsurface and well data indicate Miocene subsidence and sedimentation in the modern-day northern Mimbres, Mesilla, Malpais, and Tularosa basins (e.g., Mack, 2004). These times are coeval with initial deposition of basin fill in west Texas at 23–22 Ma (Stevens and Stevens, 1990). Thus, sedimentary basin fill and thermochronology both suggest that the majority of basins in the southern Rio Grande rift began to develop sometime during the late Oligocene and continued throughout the Miocene.

The strong geologic evidence for the onset of RGR extension at ca. 27 Ma is entirely consistent with our thermochronology data. The oldest cooling ages of 25 Ma in the RGR suggest that the main phase of Rio Grande rift extension narrowly postdated the onset of mafic magmatism in the rift. If RGR rifting initiated at 35 Ma, as was previously postulated, it should be documented in at least some of our thermochronological thermal histories. Thus, our preferred hypothesis is that the Rio Grande rift in southern New Mexico initiated ca. 27 Ma based on the eruptions of large volumes of basalt or by 25 Ma based on our thermochronological evidence, and that evidence for extension in the rift at 35 Ma is more easily explained through volcanic processes that were clearly ongoing elsewhere in the Mogollon-Datil volcanic field.

RGR formation models include (1) backarc extension related to the Farallon plate followed by extension synchronous with rapid growth of the San Andres transform; (2) oblique extension followed by rotation of the Colorado Plateau; (3) collapse of over-thickened crust following the Laramide orogeny; and (4) passage of a weak mantle plume beneath the RGR (e.g., Hamilton, 1981; Humphreys, 1995; Landman and Flowers, 2013; van Wijk et al., 2018; Abbey and Niemi, 2020). We do not address (3) and (4) because they do not predict specific times of extension that can be tested with our thermochronological data. Instead, we focus our discussion on (1) and (2) and use available age constraints and comparison with times of magmatism. Dickinson (2002) highlights the compound evolutionary history of the Basin and Range taphrogen (a region that has experienced multiple phases of extensional deformation), and he interprets the patterns of extension in space and time to have been driven by separate geodynamic forces. Within our study area, Dickinson (2002) suggests that a ca. 28.5–17.5 Ma phase of extension in the Rio Grande rift, the Texas-Chihuahua border belt, and along metamorphic core complexes in eastern Arizona were driven by backarc extension during slab rollback of the underlying Farallon plate. This was followed by a second phase of extension from 12.5 to 5 Ma that was driven by rapid southern migration of the Rivera Triple Junction and lengthening of the San Andreas fault system, causing renewed, block-style extension in the Rio Grande rift.

The timing of initial Basin and Range extension in southwestern New Mexico from thermochronology data presented in this study is constrained to older than 25 Ma based on comparison of rock cooling with ages of volcanic rocks (Fig. 10). This timeframe overlaps with extension within the Pinaleño-Jackson metamorphic core complex of southeastern Arizona, where ⁴⁰Ar/³⁹Ar biotite data from the mylonite zone document 29–19 Ma extension (Long et al., 1995). However, these ranges lack characteristics of metamorphic core complexes, such as a low-angle master detachment fault, upper-plate normal faults, a thick (0.5–3 km) ductile shear zone, and accumulation of tens of km of slip (e.g., Axen, 2004), suggesting lower magnitudes of extension. The main phase of southern Rio Grande rift extension from modeling of thermochronology data lasted roughly from 25 to 6 Ma, and all but three samples preserve a shorter period of extension from ca. 18–6 Ma (Fig. 10). Most of these samples also only indicate a single pulse of extension, as evidenced by long track lengths and similarity between AHe, AFT, and ZHe dates. These observations are at odds with models invoking two-stage uplift from 28.5 to 17.5 Ma and 12.5–5 Ma, separated by a period of tectonic quiescence (Dickinson, 2002). Rather, they are similar to other thermochronology studies of extension within the Rio Grande rift. Ricketts et al. (2016) modeled AHe and AFT data from Rio Grande rift flank uplifts and noted a similar pulse of extension from 25 to 10 Ma, supporting a model of synchronous extension along the length of the rift (Chapin and Cather, 1994; Landman and Flowers, 2013). This was supported by subsequent work where densely spaced vertical transects from the northern segment of rift indicate fault initiation at 22 Ma (Abbey and Niemi, 2018) and a more regional perspective where thermal modeling indicates 25 Ma fault initiation in both the northern and southern segments of the rift (Abbey and Niemi, 2020).

Analysis of spatial and temporal patterns of magmatism and thermochronological data from the Rio Grande rift suggest that the northern and southern segments of the rift may have initially evolved through a model of oblique strain (Abbey and Niemi, 2020), and later extension was facilitated by clockwise rotation of the Colorado Plateau in the middle to late Miocene (e.g., Chapin and Cather, 1994; Kreemer et al., 2010). Abbey and Niemi (2020) examined vertical transects in central New Mexico and found that at these locations, exhumation occurred from 25 to 8 Ma, similar to 24.9–6.7 Ma extensional exhumation from this study (Fig. 10). In contrast, Abbey and Niemi (2020) note that fault slip in the central Rio Grande rift during the middle Miocene was minimal, and extension was assisted by magmatism associated with the Jemez lineament, which helped to integrate the entire rift system. Similarly, patterns of magmatism and extension in the Basin and Range/Rio Grande rift transition may indicate separate models of extension. In the southeastern Basin and Range, extension either predated or was contemporaneous with magmatism in the Mogollon-Datil and Boot Heel volcanic fields (Fig. 10). Given the similarities between cooling in our data set and the time of extension in the Pinaleño-Jackson metamorphic core complex, we favor a model where extension was contemporaneous with magmatism. If so, magmatism may have facilitated extension in this region similar to the central Rio Grande rift where it crosses the Jemez lineament. In contrast, early Miocene magmatism is virtually absent in the southern Rio Grande rift, and extension may have formed through oblique strain. Differences in the timing and driving forces of extension support a model where the southern RGR is a distinct structural entity from the adjacent Basin and Range (Fig. 2A). These results will need to be integrated with diverse geologic and geophysical data sets to better understand the nature of the boundary between the Basin and Range and southern RGR.

A suite of low-temperature thermochronology data composed of 128 apatite and 63 zircon (U-Th)/He dates and 27 apatite fission-track dates was generated from samples collected in mountain ranges in the Rio Grande rift/Basin and Range transition zone in southern New Mexico. These dates were used to create time-temperature inverse models using the program HeFTy to investigate the onset of rapid exhumation and cooling history of footwall uplifts associated with Neogene extension.

Throughout southern New Mexico, rapid cooling occurred beginning ca. 30–25 Ma in the easternmost Basin and Range, 5–7 m.y. before the main phase of faulting in the Rio Grande rift. ZHe and AHe ages from Basin and Range locations overlap with ages of regional caldera volcanism and are likely reset but suggest that samples were either transported through the AHe and ZHe PRZs either prior to or synchronously with magmatism.

In the southern Rio Grande rift, thermal history modeling and low-temperature themochronology data support contemporaneous onset of rapid extension throughout the entire 1000 km length of the RGR beginning ca. 25–20 Ma. Old ZHe ages are only present east of the Cookes Range, and thermochronology ages do not overlap with rift magmatic ages, suggesting that cooling was driven by exhumation along normal faults in the southern RGR.

The main control on differences between the Basin and Range and RGR remains incompletely understood, but the results of this project show that the cooling histories of the two provinces in southern New Mexico are markedly different. However, the nature of the relationship between the RGR and Basin and Range would benefit from further characterization of the lithosphere and asthenosphere of southern New Mexico through additional geophysical research, such as an E-W seismic-reflection line along the transect of this study, and detailed thermochronologic work focused on the southeastern Basin and Range in areas less influenced by the ignimbrite flare-up.

This work was conducted in part for the M.S. degree research of Michelle Gavel at New Mexico State University (NMSU) and Julian Biddle at University of Texas at El Paso (UTEP). Amato was supported by National Science Foundation (NSF) grant EAR 1624575. Ricketts was supported by NSF grant EAR1624538 and a University of Texas at El Paso Research Institute grant. Gavel received funding from the NMSU Russell Clemons Field Scholarship and the Albuquerque Gem and Mineral Club Student Scholarship. Some samples were collected by NMSU students Sean Gaynor, Colby Howland, Trey Becker, and Chelsea Ottenfeld and by UTEP students Jose Adrian Rubio and Nathan Reade. Becky Flowers and Jim Metcalf at the University of Colorado TRaIL lab in Boulder, Colorado, helped acquire (U-Th)/He data. Pieter Vermeesch created the Density Plotter program to generate the KDE plots. Reviews by Nathan Niemi, Lisa Tranel, and Associate Editor James Spotila improved the manuscript and are much appreciated.