Zircon (U-Th)/He Thermochronologic Constraints on the Long-Term Thermal Evolution of Southern New Mexico and Western Texas

Thermochronology is a powerful tool to constrain the ages and durations of past geologic events because exhumation leads to cooling, the timing of which is recorded by different isotopic systems. Cooling ages are interpreted in context of closure temperatures of different minerals (e.g., [1]) in which diffusion of radiogenic daughter isotopes slows significantly below a known temperature range. However, there can be a gap in the dates obtained from high-temperature thermochronologic systems such as titanite U-Pb, hornblende ⁴⁰Ar/³⁹Ar, and mica ⁴⁰Ar/³⁹Ar, which record the timing of cooling from ~600 to 300°C [2], versus low-temperature systems such as apatite fission-track and (U-Th)/He, which record the timing of cooling below 120–60°C and 90–30°C, respectively [3, 4]. Multidiffusion domain (MDD) analysis of K-feldspar using the ⁴⁰Ar/³⁹Ar system [5] could fill in the temperature gap between biotite ⁴⁰Ar/³⁹Ar and apatite (U-Th)/He, although it has been shown to be problematic in some instances (e.g., [6–8]), whereas MDD analysis of muscovite shows promise [9, 10]. Zircon fission-track, sensitive to temperatures of ~270–210°C [11], is a widely used thermochronometer to partially fill this temperature range, although there can still be a gap in the thermochronologic history between ⁴⁰Ar/³⁹Ar methods and apatite (U-Th)/He because of variations in closure temperature related to cooling rate.

Limitations in deciphering thermochronologic histories related to this gap may be overcome by the work of Guenthner et al. [12], who describe a helium diffusion model that incorporates an important relationship between zircon (U-Th)/He (ZHe) dates and radiation damage accumulation and annealing in zircon crystals. This relationship is typically expressed by intrasample ZHe dates that sometimes span hundreds of millions of years. For grains that have experienced identical thermal histories, differences in ZHe dates are governed by differences in the effective uranium (⁠$eU=U+0.235Th$⁠) that are reflected as either positive or negative ZHe date-eU relationships. These relationships result from differences in the ability of each grain to retain helium that are dependent on eU values. Positive and negative ZHe date-eU trends are the cumulative result of the overall thermal history and the effects of radiation damage on helium diffusion, and these allow for continuous thermal histories to be reconstructed. Helium diffusion models in zircon suggest a temperature sensitivity window ranging from approximately 210 to 50°C [12, 13], which partially overlaps with zircon fission-track and K-feldspar ⁴⁰Ar/³⁹Ar data at the higher range and with apatite fission-track and (U-Th)/He techniques at the lower end. The large ZHe temperature sensitivity window makes it possible to constrain more complete and continuous thermal histories, offering an opportunity to bridge the gap between higher and lower temperature methods. Previous studies have relied on this intrasample spread in ZHe dates to investigate a wide range of geologic processes, such as the Proterozoic thermal history of rocks (e.g., [14]), timing of Laramide exhumation [15], and development of the South American passive margin [16]. However, additional studies are needed to further assess the efficacy of this method.

Proterozoic crystalline rocks of southwestern North America have experienced a complex tectonic history that spans over a billion years (e.g., [17]). In southwestern New Mexico and western Texas, this history includes, but is not limited to, Mesoproterozoic and Neoproterozoic assembly and break-up of the supercontinent Rodinia, Pennsylvanian-Permian Ancestral Rocky Mountain deformation, widespread Paleozoic and Mesozoic sedimentation, latest Cretaceous-Eocene Laramide compression, and culminating with Neogene extension related to development of the Rio Grande rift (e.g., [18–22]).

⁴⁰Ar/³⁹Ar cooling ages from Proterozoic crystalline rocks in New Mexico typically range from approximately 1600 to 1000 Ma and reflect thermal pulses during intracontinental tectonism and plutonism (e.g., [23, 24]). However, apatite fission-track and (U-Th)/He studies from the same region in New Mexico yield dates that are typically younger than about 100 Ma, reflective of cooling during exhumation associated with the Laramide orogeny and Rio Grande rift extension [25–28]. These two disparate datasets each provide brief snapshots into a more complex thermal history, but a continuous thermal record remains lacking, particularly in the range from about 1000 Ma to younger than 100 Ma. Thus, this region of the southwest U.S. is an exceptional natural laboratory for investigating long-term thermal histories using ZHe thermochronologic methods.

Here, we present the first ZHe dates from three ranges in southern New Mexico and west Texas: the Cookes Range in New Mexico and the Franklin Mountains and Carrizo Mountains of west Texas (Figure 1). ZHe dates for each sample have significant intrasample ZHe date variability which we use to produce a series of forward and inverse thermal history models. We also present ZHe data from several detrital samples and discuss the ways in which these data compliment ZHe dates from crystalline samples. Our new thermal history models capture multiple pulses of cooling and reheating that can be directly tied to known geologic events and provide an important link between disparate higher-temperature ⁴⁰Ar/³⁹Ar cooling ages and lower-temperature apatite fission-track and (U-Th)/He datasets. Finally, we evaluate the use of both inverse and forward models to compare the long-term thermal history of rocks in this region and determine an improved application of the ZHe thermochronology method.

From 1.8–1.6 Ga, Paleoproterozoic growth of southwest Laurentia occurred by progressive accretion of arc terranes, including the Yavapai and Mazatzal Provinces [29–31], where final assembly of the supercontinent Rodinia occurred during Grenville tectonism in the late Mesoproterozoic (Figure 1) [32]. The Grenville orogeny in southwestern Laurentia records arc accretion and continent-continent collision between 1350 and 980 Ma [33]. Some of the Grenville orogenic foreland is exposed in western Texas in the Carrizo Mountains [33, 34]. The ca. 1380–1327 Ma Carrizo Mountain Group (Figure 1(b)), composed of immature clastic rocks, within-plate rhyolitic volcanic rocks, and minor carbonates, likely records rifting of continental crust within a back-arc basin during overall convergence [35–37]. The Franklin Mountains (Figure 1(b)) record continued sedimentation within this back-arc basin from ca. 1260–1240 Ma [36] until continent-continent collision ca. 1150–1120 Ma [36]. North of the Grenville deformation front (Figure 1), magmatism in the Franklin Mountains is represented by ~1.1 Ga plutonic and volcanic rocks, including the Red Bluff granite which was targeted in this study for ZHe dating [38, 39]. Petrological and geochemical studies classify the Red Bluff granite as a within-plate, A-type granite and support a model where these rocks were emplaced within a more regional strain field dominated by NW-SE shortening and orthogonal NE-SW extension related to Grenville convergence [39, 40]. Final shortening along the southern margin of Laurentia occurred between ~1035 and 980 Ma [33, 41].

Models for the breakup of Rodinia suggest diachronous disassembly with early rifting occurring between 780 and 680 Ma, followed by the main rifting phase between 620 and 550 Ma [18, 21, 42, 43]. Supercontinent breakup may have also been accompanied by a major pulse of denudation recorded in ZHe datasets [14]. Breakup was followed by deposition of Paleozoic passive margin sediment on Precambrian basement (e.g., [44]) to create the Great Unconformity, a globally significant feature in the rock record. In southern New Mexico and western Texas, Proterozoic crystalline rocks beneath the Great Unconformity typically consist of 1.45-1.35 and some 1.1 Ga granites intruded into 1.68-1.65 supracrustal rocks of the Mazatzal Province [17]. These rocks are generally overlain by the Cambrian-Ordovician Bliss Sandstone, which define a Great Unconformity that encompasses 600 m.y. to more than 1,000 m.y. of missing time.

Late Paleozoic deformation in western North America resulted in the development of the Ancestral Rocky Mountains (e.g., [19, 45]). In the Cookes Range (Figure 1), Ancestral Rocky Mountain deformation is evident by a near absence of sedimentation from approximately 320 to 290 Ma [19]. In contrast, the Franklin Mountains record near continuous sedimentation during this timeframe [46], suggesting minimal exhumation. In west Texas, deformation locally stripped early Paleozoic sedimentary rocks and exposed Proterozoic basement which was subsequently covered by early Permian rocks [47].

Starting in the Late Jurassic, a continental rift developed in the region of northern Sonora, Mexico, and southern Arizona (e.g., [48]). The resulting rift basin was originally termed the Bisbee basin [49] and now referred to as the Mexican Border Rift (Figure 1) [50]. Marine strata have Tethyan fossils indicating a connection to the Gulf of Mexico (e.g., [51]). Volcanogenic material such as mafic volcanic rocks and pillow basalts were deposited in this rift [52–54]. Several of the samples from this study were collected from the northern rift shoulder in southwestern New Mexico. Rifting in this region ended by the middle Cretaceous based on subsidence studies indicating a transition from rift-related thermal subsidence to the formation of a foreland basin beginning in the Albian [50, 55].

Late Mesozoic to Eocene deformation during the Laramide orogeny involved northeast-directed crustal shortening that resulted in a series of uplifts and basins [22]. These basins and uplifts trend northwest in southern New Mexico and extend into northern Mexico and western Texas; our Carrizo Mountains study site lies at the easternmost extent of Laramide deformation (Figure 1). These structures were extensively dissected by younger continental extension related to development of the Rio Grande rift [22]. Initial rifting in southern New Mexico began during the Oligocene to produce a series of fault-block uplifts and basins that create the modern topographic grain of the region [20]. The Cookes Range is bounded on three sides by normal faults related to development of the Rio Grande rift, and the Franklin Mountains and Carrizo Mountains are both associated with active normal faults that document continued deformation in this region.

Uplifts within the southern Rio Grande rift, such as the Carrizo Mountains, Franklin Mountains, and Cookes Range, expose Proterozoic crystalline basement that has been exhumed to the surface in the footwalls of Cenozoic normal faults. Samples were collected from crystalline basement rocks and clastic coarse-grained sandstone or granule conglomerate sedimentary rocks because they likely contain abundant zircon crystals. In addition, sample locations were targeted to avoid possible effects of widespread Oligocene plutonism. Sample collection along an E-W transect across southern New Mexico suggests that many samples farther to the west have been completely reset, and all ZHe dates are Oligocene-Miocene [25, 56, 57]. These studies found that along this EW transect, the Cookes Range was the westernmost location to preserve older ZHe dates, although it does contain Oligocene intrusive and extrusive rocks (Figure 1(b)) that may have partially affected the ZHe data (as described in more detail below). Therefore, mountain ranges east of and including the Cookes Range were selected because they are removed from the effects of Oligocene magmatism and likely retain a record of the older, long-term thermal history.

Whole-rock samples were processed using standard mineral separation techniques to isolate the zircon fraction. Zircon separates were inspected with a petrographic microscope equipped with a digital camera used for measuring grain dimensions, including length, width, and depth. Ideal zircon crystals for ZHe thermochronology are euhedral with no inclusions and a minimum diameter of 70 μm. Five to ten of the best zircons from each sample were hand-picked and loaded into Nb tubes for analysis. All samples were analyzed at the Thermochronology Research and Instrumentation Laboratory at CU Boulder (CU TRaIL).

Since all samples yield a range of ZHe dates and eU values, a forward modeling approach was first used to constrain the long-term thermal history of each mountain range since crystallization of the sample. Forward modeling allows for the calculation of a thermochronometric age from a determined time-temperature path, using helium diffusion or annealing kinetics [58]. Multiple hypothetical time-temperature paths are investigated that are based on known geologic constraints, including formation of the Great Unconformity, periods of Paleozoic and Mesozoic sedimentation, and possible exhumation during Ancestral Rocky Mountain, Laramide, and Rio Grande rift deformation events (Table S1, Fig. S1). For each of these paths, we vary maximum or minimum temperatures achieved during different segments of the sample’s thermal history. ZHe date-eU curves are calculated from each t-T path using a Matlab script (Guenthner, pers. comm., 2018) that incorporates the zircon radiation damage accumulation and annealing model (ZRDAAM) from Guenthner et al. [12]. Inputs include a specific time-temperature path, zircon eU values, and zircon grain size. These date-eU paths are then compared to ZHe data. Model date-eU paths that are drastically different than the observed data are not viable thermal histories. By varying different segments of the thermal history, the model date-eU paths from each sample are incrementally adjusted to provide a range of permissible thermal histories. This approach assumes that the only control on ZHe dates is helium diffusion related to crystal damage and annealing that evolves through time. However, the spread in ZHe dates could also be influenced by other factors such as crystal size, U and Th zoning in the zircon crystal, and implantation of helium from neighboring grains (e.g., [59–61]). Although possible He implantation is generally not known, and zoning information is not typically obtained in ZHe studies, the effects of crystal size on resulting date-eU curves can be investigated. To do this, we use a method similar to other workers (e.g., [14, 62, 63]) in which three separate date-eU curves are modeled to create an envelope within which the ZHe dates should lie. To create the three curves, we use the $meangrainsize±2standarddeviations$⁠.

We also use inverse modeling to further refine possible thermal histories. Whereas forward modeling was used to test tens of possible paths, an inverse modeling approach can test tens of thousands of possible paths and explore additional possible t-T paths that are permissible by the ZHe data. Inverse modeling allows for the calculation of t-T paths that match measured thermochronometric ages to within a specified amount of statistical error [58]. We use HeFTy v. 1.9.3 [4], which uses a Monte Carlo method to plot possible time-temperatures paths based on entered data. Input data into HeFTy includes U, Th, Sm concentrations, grain radius, measured age, and age uncertainty. HeFTy requires a set of temperature-time parameters for modeling, such as start times, end times, and temperatures. All inverse models explore a total of 50,000 paths, where resulting “acceptable” paths have goodness-of-fit parameters >0.05, and “good” paths have goodness-of-fit criteria >0.5. Complete model constraints and inputs are available in Table S1 and Figure S1 and follow methods outlined in Flowers et al. [64].

Finally, we explore the significance of detrital ZHe dates from samples collected from the Cookes Range and Franklin Mountains. In contrast to igneous rock samples, partially reset detrital zircon grains only share a common post-depositional thermal history, but each individual grain may still retain an older, unique thermal imprint. As a result, they should not necessarily lie along a single date-eU curve, but will instead have potentially highly variable ZHe dates [15]. Analysis of detrital grains entails constructing inheritance envelopes for each location, following methods outlined in Reiners et al. [65] and Guenthner et al. [15]. The inheritance envelope is created by combining a zero-inheritance date-eU curve with multiple maximum inheritance curves that assume different zircon crystallization ages. The zero-inheritance date-eU curve is constructed from the thermal history that begins at the depositional age of the detrital sample and assumes no inherited radiation damage at that time, a situation that can be achieved by either complete resetting of detrital grains at the time of deposition or zero U-Pb ages at the time of deposition. Maximum inheritance curves incorporate the effects of inherited radiation damage by assuming zircon crystallization at surface temperatures prior to the depositional age. These grains then share an identical thermal history as the zero-inheritance curve after deposition. When combined, these different curves yield an inheritance envelope that should contain the observed ZHe dates if the t-T path is a plausible result of the sample’s thermal history.

A total of 37 new ZHe dates are presented, including 10 dates from the Carrizo Mountains, 20 dates from the Franklin Mountains, and seven dates from the Cookes Range. These are combined with 18 ZHe dates previously reported [25, 56], including nine dates from the Franklin Mountains and nine dates from the Cookes Range (Figure 2). In total, this project incorporates 55 individual grain ZHe dates from eight samples to investigate patterns in long-term thermal histories spanning hundreds of millions of years (Table 1).

Ten ZHe dates from two samples were obtained from the Mesoproterozoic Hackett Peak Formation of the Carrizo Mountain Group in the northern Carrizo Mountains (Figures 1(b) and 2, Table 1). These data show a wide range in both ZHe date (6–731 Ma) and eU concentration (52–1729 ppm). There is also a well-defined negative ZHe date-eU correlation that shows a steep trend until eU values of approximately 400 ppm. ZHe dates show a slight positive relationship with grain radius. The oldest ZHe dates correspond to the largest grains, suggesting a possible grain size control on ZHe dates.

Twenty-nine ZHe dates from four samples were obtained from the Mesoproterozoic Red Bluff granite and the Cambrian-Ordovician Bliss Sandstone in the Franklin Mountains (Figures 1(b) and 2). ZHe dates from the Red Bluff granite have a large range from 19 to 401 Ma. Their eU concentrations show a large spread as well, ranging from 63 to 828 ppm. Together, these data have a slightly negative, although poorly defined, ZHe date-eU correlation, where older ZHe dates typically correspond to lower eU values. However, individual samples from the Franklin Mountains do not always show a negative trend, and some (17FR03) show more of a flat trend. The majority of crystalline zircon grains from the Franklin Mountains cluster within a narrow size range of 35-55 μm such that a ZHe date-grain size trend is not observed. ZHe dates from the Bliss Sandstone show a ZHe date range of 62–649 Ma, with eU concentrations ranging from 90 to 374 ppm. Many of these ZHe dates are older than ZHe dates from crystalline basement samples, and together, they define a negative trend date-eU trend. This sample displays a wider range in crystal sizes than the crystalline grains, where most of the grains are older and larger. However, ZHe dates do not appear to increase with larger crystals.

A total of 16 ZHe dates are presented from the Cookes Range (Figures 1(b) and 2). Ten ZHe dates are from a Proterozoic granite exposed in Rattlesnake Ridge along the southern edge of Cookes Range. These ZHe dates range from 63 to 446 Ma with an eU range of 232–945 ppm. These data show a steep negative ZHe date-eU relationship at lower eU values, which transitions to a flatter trend at higher eU values.

An additional six ZHe dates were obtained from the Permian Abo Formation from the Cookes Range. ZHe dates range from 44 to 130 Ma, with corresponding eU values from 56 to 351 ppm. Although this sample does not yield as large of an eU range, these ZHe dates have a slight negative trend where older ZHe dates correspond to lower eU values. In contrast to Franklin Mountains samples, detrital grains from the Cookes Range are mostly younger than crystalline zircon grains. Individual ZHe dates do not seem to show a well-defined trend with grain size.

Hypothetical time-temperature paths were constructed for crystalline samples at each location for two purposes (Figures 3–5). First, these paths are used to eliminate thermal histories that are incompatible with the observed ZHe dates and begin to refine possible thermal histories. Second, these paths are used to test which periods of burial or uplift have the largest effect on the resulting ZHe date-eU pattern and which events have relatively minor effects. To do this, we used a Matlab routine that incorporates the zircon radiation damage accumulation and annealing model (ZRDAAM) from Guenthner et al. [12]. Forward model outputs are calculated ZHe date-eU curves that are predicted from the input data. Hypothetical time-temperature paths were constructed for each mountain range by varying the amount of burial or exhumation from different deformational events in that area. Specific geologic constraints for each location are provided in Table S1 and Figure S1. Each ZHe date-eU curve calculated from an input time-temperature is then compared to the ZHe data collected from each location. Below, we first present preliminary forward models for each of the three study locations. Then, we refine these results and produce best-fit models that represent the closest match to the observed data acquired through this approach.

Four sets of representative t-T paths were constructed for the Carrizo Mountains that each varies a different segment of the history. The first set of representative paths varies the timing of Proterozoic exhumation prior to deposition of Paleozoic strata (Figure 3(a)). Three possible thermal histories are examined to test how sensitive ZHe data are to Proterozoic exhumation. These three paths all cool through the range of 500 to 350°C at 1035 Ma, consistent with ⁴⁰Ar/³⁹Ar hornblende and muscovite data [41]. From there, they diverge to include early exhumation to 15°C at 1000 Ma (blue path), intermediate exhumation from 300 to 15°C from 800 to 700 Ma (teal path), and late exhumation from 300 to 15°C from 600 to 500 Ma (yellow path) (Figure 3(a)). The post-500 Ma segments of these paths are identical. These three paths each yields a date-eU curve that is notably different than the others, suggesting that the Proterozoic cooling history imparts a significant influence on the resulting date-eU relationship. Of the three proposed cooling scenarios, rapid exhumation to the surface from 1050 to 1000 Ma predicts maximum ZHe dates that are similar to the observed dates, although it does not match the negative slop of the data. However, this result supports rapid cooling followed by prolonged residence at the Earth’s surface.

The second set of representative paths builds off the first iteration by using 1050–1000 Ma exhumation from Figure 3(a). Here, we vary the depth of maximum Paleozoic burial after formation of the Great Unconformity (Table S1, Figure 3(b), and Figure S1). Eight thermal histories are examined that vary the maximum burial temperature from 110 to 180°C at 325 Ma, signifying maximum burial before Ancestral Rocky Mountain uplift. The resulting ZHe date-eU patterns vary slightly in width, but there is a larger variation in the maximum ZHe age. However, maximum burial depths of less than 140°C all yield almost identical date-eU paths, suggesting that this dataset is insensitive to lower temperatures. Modeled date-eU paths of 110–140°C predict maximum ZHe dates that are most similar to the observed data.

The third set of representative paths incorporates 1050–1000 Ma Proterozoic cooling and uses a maximum Paleozoic burial temperature of 120°C. This iteration varies depths of Mesozoic burial after Ancestral Rocky Mountain uplift (Table S1; Figure S1). Eight thermal histories are examined that vary maximum depth from 130 to 200°C at 80 Ma (Figure 3(c)). These t-T paths yield date-eU curves with a wide range in maximum ZHe dates and widths. When compared to the observed data, there is a trade-off between matching the oldest ZHe dates and matching the negative date-eU slope of the data. For example, a temperature of 160°C predicts oldest ZHe dates that match the observed data, while a temperature of 170-180°C predicts a negative slope that closely follows the data. Of the calculated date-eU curves, a maximum Mesozoic burial depth of 170°C gives a predicted curve that somewhat matches the steep negative date-eU slope in the observed data, as well as the kink where the slope changes to a much shallower trend.

The final set of representative paths was constructed to vary the transition from Laramide shortening-related uplift to Rio Grande rift extensional uplift. To do this, we vary the temperature at which the sample resided after Laramide uplift and prior to Rio Grande rift uplift (Figure 3(d)). Eight thermal histories are examined that vary the temperature between 60 and 130°C from 40 to 30 Ma. All of the paths are nearly identical, suggesting that the overall ZHe date-eU pattern is not as sensitive to this younger deformational event in the Carrizo Mountains.

Neogene faulting and tilting of the Franklin Mountains expose Proterozoic granite along the eastern base of the range (Figure 1(b)). A total of three Proterozoic samples from the Franklin Mountains were collected from various locations along the range. Although they were all collected from similar elevations, calculated depths for each sample beneath the Great Unconformity varied. For example, although samples 15FR03 and 17FR03 resided at paleodepths of 361 m and 229 m, respectively, sample 17FR04 was located at a significantly greater depth of approximately 1587 m. For modeling purposes, samples 15FR03 and 17FR03 were combined into a single sample because they likely experienced near-identical thermal histories, and together, these grains show a larger range in eU values that allows for tighter constraints on the thermal history. ZHe dates from sample 17FR04 were excluded, and no attempt was made to model this sample, in part, because only four ZHe dates are available for this sample with limited spread in eU values (Figure 2, Table 1).

Four sets of representative paths were constructed for the Franklin Mountains (Figure 4). The four sets of paths are presented in the order of decreasing date-eU curve variability. For example, Mesozoic burial yields date-eU curves that vary widely in terms of maximum ZHe date and the location of the steep negative trend, whereas different Proterozoic paths have little effect on date-eU curves. The first set of representative paths was constructed to vary the amount of Mesozoic burial after regional uplift ceased in the Cretaceous (Figure 4(a)). Eight paths are examined that vary the maximum burial temperature between 80 and 220°C at 80 Ma, signifying maximum burial before Laramide uplift. These paths produce date-eU curves with marked differences. The most notable differences between these curves are the maximum ZHe date and the position of the negative slope, suggesting that this segment of the path has a significant effect on the final ZHe date-eU curve. Of the eight proposed t-T paths, burial to 180°C yields a date-eU curve that closely aligns with the observed data.

The second set of representative paths was constructed to vary the amount of Paleozoic burial after the Great Unconformity. This path includes a period of constant temperature from 250 to 150 Ma because in the Franklin Mountains, lower Cretaceous rocks unconformably overlie Permian rocks, suggesting a period of no deposition [66]. Eight thermal histories are examined that vary the maximum burial temperature between 40 and 180°C from 250 to 150 Ma (Figure 4(b)). A minimum temperature of 40°C was chosen to represent minimal Paleozoic burial, and a maximum temperature of 180°C was chosen because it is near the higher end of the ZHe temperature sensitivity window. The resulting date-eU curves are similar at eU values greater than 400 ppm. At lower values, the curves diverge to form peaks at different ZHe dates. Although multiple resultant date-eU curves roughly match the observed ZHe dates, an intermediate burial temperature of 120°C produces a curve that bisects the data.

Next, the transition from Laramide shortening-related uplift to Rio Grande rift extensional uplift is varied (Figure 4(c)). Eight possible thermal histories are examined that vary post-Laramide temperatures between 40 and 180°C at 40 Ma. Higher burial temperatures yield date-eU curves that are nearly flat. Burial temperatures less than about 120°C produce curves with little variation, yet these curves predict oldest ZHe dates that more closely align with the observed data. A temperature of 60°C is selected because it produces a curve with the highest ZHe peak.

In the final iteration, the Proterozoic thermal history is examined (Figure 4(d)). All of the possible paths are forced to near-surface temperatures shortly after crystallization because the Red Bluff granite intrudes the overlying Thunderbird Rhyolite in some locations, suggesting it was emplaced at shallow levels in the crust. From here, eight different paths are examined that postulate that the Red Bluff granite may have been buried and then exhumed back to the surface prior to deposition of Paleozoic strata. Maximum burial temperatures vary from 15 to 200°C. In addition, two different times of exhumation are investigated, including exhumation from 800 to 750 Ma and exhumation from 650 to 550 Ma. These different t-T paths give predicted date-eU curves with little variability, suggesting that the observed ZHe date-eU values are not strongly affected by the Proterozoic thermal history of these samples. This is consistent with the observed ZHe dates, which are mostly younger than 250 Ma.

Four sets of representative paths were constructed for Cookes Range. The first set of paths was constructed to vary the amount of Paleozoic burial after the Great Unconformity. Eight thermal histories are examined that vary the maximum burial depth to temperatures that range from 60 to 200°C at 325 Ma, signifying maximum burial before Ancestral Rocky Mountain uplift (Figure 5(a)). The higher temperature t-T paths yield date-eU curves with prominent flat pediments at intermediate eU values. This pediment disappears with lower burial temperatures. Although the observed ZHe date-eU relationships do not show this pediment, the higher burial temperatures provide a better match because the maximum ZHe dates are similar to the observed dates, whereas lower burial temperatures predict ZHe dates older than 1000 Ma.

The second set of paths was constructed to vary the amount of Mesozoic burial after Ancestral Rocky Mountain uplift. These paths use a Paleozoic burial temperature constraint of 160°C (Figure 5(a)). Eight thermal histories are examined that vary the temperature from 60 to 200°C (Figure 5(b)). The resulting date-eU curves show dramatic variability. At lower burial temperatures, a prominent pediment is observed, although this pediment gradually diminishes at higher temperatures. A Mesozoic burial temperature of 160°C produces a date-eU curve that lies along the path of observed ZHe dates and captures the broad date-eU pattern at intermediate and high eU values (200–1000 ppm).

In the third iteration, Paleozoic and Mesozoic burial temperatures are held at 160°C, and the transition from Laramide shortening to Rio Grande rift extensional uplift is varied from 20 to 160°C (Figure 5(c)). Lower temperatures of 20–120°C produce nearly identical date-eU curves. The highest temperature of 160°C predicts oldest ZHe dates that are too young to match the observed data. Because lower temperatures appear to be insensitive to the resulting date-eU pattern, a temperature of 80°C closely follows the observed dates and is used in the final iteration.

Finally, the three previous iterations are used to constrain the Phanerozoic segments of the t-T history, and the Proterozoic uplift history is varied (Figure 5(d)). Five possible histories are evaluated that vary the timing of exhumation to the surface. Three paths of either early cooling or constant cooling are nearly identical at eU values greater than 250 ppm. At lower values, they diverge. t-T paths involving rapid cooling younger than 1000 Ma produce distinctively broad date-eU patterns. The three early cooling scenarios yield date-eU curves that are similar to the observed ZHe date-eU trends. The only significant difference between these three date-eU curves is the predicted maximum ZHe date. However, no ZHe grains are available with eU values lower than 200 ppm where these paths diverge.

The above analysis serves as a useful guide for interpreting ZHe datasets. For all three locations, the final forward models are nonunique such that our approach to refine the thermal history only manages to produce a single t-T path out of many. However, this approach does highlight which segments of the t-T history the resulting date-eU curve is most sensitive to. For example, at all three locations varying the magnitude of Mesozoic burial has a drastic effect on the resulting date-eU curves. In the Carrizo Mountains and Cookes Range, the timing of Proterozoic exhumation also produces date-eU curves that show significant differences. In contrast, at all three locations, the post-Laramide/pre-Rio Grande rift temperature does not appear to significantly affect the date-eU curves. This analysis highlights some of the uses of a forward modeling technique for understanding which events may be preserved in the thermal record and which events the data are not sensitive to. Below, we use these results and incorporate the possible effects of grain size to produce our best-fit forward models.

The preliminary forward modeling efforts above fail to yield a single t-T path that adequately aligns with the observed date-eU patterns (Figures 3–5). These results suggest that the effects of radiation damage and annealing on helium diffusion through the crystal lattice [12] may not be the only control on the observed ZHe dates. As described above, other factors that could possibly influence ZHe dates are crystal size, U and Th zoning in the zircon crystal, and implantation of helium from neighboring grains. Although the effects of zoning and implantation are unknown, for some samples, grain size does seem to be related to ZHe date (Figure 2). Here, we produce refined forward model t-T paths for each location that incorporate possible ZHe date effects related to differences in grain size. For each t-T path, we create three separate date-eU curves using the $averagegrainsize±2standarddeviations$⁠, based on similar published methods (e.g., [14, 62, 63]). These curves define a date-eU envelope that shows the combined effects of radiation damage and grain size.

For the Carrizo Mountains, we use the average grain size of $53±17μm$ (±2 standard deviations) to create a date-eU envelope of the best-fit preliminary curve (Figure 6(a)). This t-T path includes cooling to surface temperatures from 1050 to 1000 Ma, early Paleozoic burial to 120°C, exhumation to surface temperatures by 285 Ma, Paleozoic through Mesozoic reheating to a maximum temperature of 170°C at 80 Ma, and a post-Laramide temperature of 60°C before final cooling to surface temperatures. This best-fit t-T path yields a date-eU envelope that encompasses all but the oldest ZHe date.

The best-fit forward model for the Franklin Mountains samples uses an average zircon grain size of $43±10μm$ to produce a date-eU envelope (Figure 6(b)). This t-T path includes no Proterozoic burial, Paleozoic burial depths to a temperature of 120°C, renewed Mesozoic burial to a maximum temperature of 180°C, post-Laramide residence at a temperature of 60°C, and final uplift to the surface during Rio Grande rift extension. The best-fit t-T path yields a date-eU envelope that encompasses nine of the 18 total ZHe dates. This model fails to account for the oldest as well as the youngest ZHe dates at low eU values. The large spread in ZHe dates from 19 to 402 Ma at eU values <200 ppm suggests that other factors, such as crystal zoning, may also significantly affect ZHe dates in the Franklin Mountains.

The best-fit forward model for the Cookes Range uses zircon grain sizes of $46±9μm$ (Figure 6(c)). Unlike the other two locations, here, we show the results of three separate best-fit models. Two models involve early cooling, whereas the third model is constant cooling from the time of crystallization. Each of the resulting date-eU envelopes are nearly identical at intermediate and high eU values, but diverge at low eU values where no ZHe dates are available. The resulting best-fit t-T paths encompass nearly all of the observed ZHe dates. However, the yellow and green paths are slightly skewed towards lower eU grains and do not encompass the oldest ZHe date or the three grains at eU values of ~300 ppm. The blue path appears to yield a slightly better fit with the overall ZHe date-eU trend.

The forward modelling method described above is successful at finding possible t-T paths that yield date-eU curves that are consistent with the observed data. However, with that approach, it is only feasible to investigate a limited number of possible t-T paths. Here, an inverse modeling approach is used to continue to refine the thermal histories and explore additional possible t-T paths that may have been overlooked during the forward modeling analysis. In this section, we present a single inverse model for each of the three study locations. The inverse models use synthetic grains that represent binned averages of the observed ZHe data. This is necessary because of the large number of individual ZHe dates for each location and the limited number of inputs available in HeFTy. Complete inverse modeling methods and data inputs are provided in Table S1 and Figure S1.

The final inverse model for the Carrizo Mountains produced nine paths with a good fit to the input synthetic grains (Figure 7(a)). These paths suggest that after metamorphism, the sample remained at elevated temperatures until rapid uplift at approximately 1035 Ma based on ⁴⁰Ar/³⁹Ar hornblende and muscovite ages [41]. ZHe data also document the timing of this pulse of exhumation and suggest that the sample was cooled to near-surface temperatures at this time (Figure 7(a)). However, some t-T paths also suggest that the sample could have remained at shallow levels in the crust and was further exhumed to the surface during a second event at approximately 600–500 Ma. Increasing temperatures during Paleozoic sedimentation are not well constrained by the model, as maximum burial temperatures from the nine t-T paths range from 50 to 180°C prior to Ancestral Rocky Mountain uplift which exposed Proterozoic basement to the surface by around 285 Ma. The Mesozoic segment of the inverse model is tightly constrained and suggests that the sample reached maximum temperatures of 150–160°C. Inverse modeling shows minimal cooling associated with the Laramide orogeny, and most paths remain at temperatures >130°C until approximately 30 Ma, when the sample was rapidly cooled to near-surface temperatures.

A total of 13 good paths were obtained from inverse modeling of samples collected from the Red Bluff granite (Figure 7(b)). Field relationships suggest that this granite was emplaced at relatively shallow crustal levels, as it intrudes into and cuts an ignimbrite within the Thunderbird Group. These volcanic rocks are interpreted to be the erupted equivalent of the Red Bluff granite, which is supported by geochemical evidence [67] and overlapping geochronologic ages [38]. Inverse models suggest that the Red Bluff granite was then reheated to temperatures ranging from 50 to 250°C. While this large uncertainty likely suggests the ZHe data from the Franklin Mountains may be largely insensitive to the Proterozoic history, the observation that all good paths require some amount of heating is consistent with a model of possible Proterozoic reburial following crystallization. All 13 paths show relatively monotonic reheating throughout the Paleozoic and Mesozoic, and converge within a temperature range of 155–190°C by the time of incipient Laramide deformation. These paths then cool rapidly from approximately 80 to 50 Ma to temperatures less than 100°C. A period of slower cooling is observed beginning at 50 Ma, and by 25–15 Ma, all paths show a final rapid pulse of cooling to surface temperatures.

The inverse model for the Cookes Range yielded a total of nine good paths that are consistent with the available geologic constrains and synthetic grains (Figure 7(c)). Following crystallization of this sample, all t-T paths show rather monotonic Proterozoic cooling towards surface temperatures. However, all t-T paths preserve a pulse of more rapid cooling between 800 and 500 Ma that terminated at near-surface temperatures. As with previous inverse models from the Carrizo and Franklin Mountains, the Paleozoic thermal history of this sample is relatively unconstrained. However, the nine paths suggest the sample reached maximum burial temperatures that range from 85 to 180°C prior to Ancestral Rocky Mountain uplift. The uplift history during the Ancestral Rocky Mountain deformational event is also relatively unconstrained, with uplift temperatures ranging from 40 to 165°C. However, during renewed heating during Mesozoic burial, the paths converge to a maximum temperature of 140–180°C. This sample shows a strong component of cooling that is coeval with Laramide deformation. All paths display rapid cooling to lower temperatures of 50–70°C beginning 60–70 Ma and ending by approximately 45 Ma. The sample remained within this temperature range until approximately 25 Ma, when it was cooled to near-surface temperatures.

Two Paleozoic detrital samples were collected for ZHe thermochronology from the Cookes Range and the Franklin Mountains. ZHe data from these locations were obtained in order to compare to ZHe data from crystalline basement samples, test whether zircon grains had been reset after deposition, and test for possible influences of magmatism. A total of six ZHe dates were obtained from the Permian Abo Formation from the Cookes Range [25]. Here, the Abo Formation consists of fine- to medium-grained sandstone. An additional seven ZHe dates were obtained from the Cambrian-Ordovician Bliss Sandstone in the Franklin Mountains, including three presented by Biddle et al. [25]. Ideally, more detrital zircon grains would have been selected, but these samples did not yield sufficient zircon grains that were euhedral and of sufficient size.

For both samples, we rely on the inverse modeling constraints described above to aid with the construction of inheritance envelopes. To do this, we use the best-fit path from the inverse models as the assumed post-depositional thermal history of the detrital samples (Figure 8). The best-fit path highlighted in red in Figure 8 represents the thermal history used to create the zero inheritance curves. The maximum inheritance curves are created using published detrital zircon U-Pb ages for the Permian Abo Formation and Cambrian-Ordovician Bliss Sandstone in southern New Mexico. Below, we describe each sample in more detail.

The Abo Formation in the Cookes Range lies approximately 600 m above the Proterozoic granite based on thicknesses of older Paleozoic units [68]. Assuming a geothermal gradient of 25°C/km, this is equal to a temperature difference of 16°C. This difference is likely minimal in the resulting date-eU curves such that the thermal history of the Proterozoic granite is likely similar to the post-depositional thermal history of the overlying detrital sample. The zero-inheritance date-eU curve that is constructed using the best-fit path from the inverse model is combined with three other maximum inheritance curves to create an inheritance envelope (Figure 8(a)). The maximum inheritance curves use zircon crystallization ages of 1250 Ma, 1465 Ma, and 1785 Ma, based on peaks in U-Pb detrital zircon ages from the Abo Formation in southern New Mexico [69]. When compared with the inheritance envelopes, detrital ZHe dates all fall along the zero-inheritance curve or below the inheritance envelope (Figure 8(a)). If the post-depositional thermal history is correct, then these observations would suggest that the three grains that lie along the zero-inheritance curve could be derived from a ca. 280–300 Ma source. However, probability density plots do not show abundant Paleozoic detrital zircon U-Pb ages [69]. This would also not account for the three ZHe dates that fall below the inheritance envelope. Alternatively, all ZHe detrital grains could be derived from Proterozoic sources if the post-depositional history is modified to include a thermal pulse to partially reset ZHe dates. A granodioritic stock and several basaltic and dacitic dikes and sills with reported K-Ar and ⁴⁰Ar/³⁹Ar ages that range from 38 to 45 Ma [68, 70, 71] intrude Paleozoic and Mesozoic rocks in the Cookes Range [68], possibly resulting in partial resetting of zircon grains. Partial resetting of detrital ZHe dates from the Abo Formation is also supported by the date-eU plots in Figure 2, where Abo Formation ZHe dates are all younger than the Proterozoic ZHe dates at low eU values.

The Bliss Sandstone sample was collected approximately 300 m stratigraphically above the Proterozoic samples used for inverse modeling, suggesting a temperature difference of only 7–8°C, assuming a geothermal gradient of 25°C/km. The best-fit inverse path should therefore be a good approximation of the postdepositional thermal history of the Bliss Sandstone (Figure 8(b)). Maximum inheritance curves were constructed using zircon crystallization ages of 1250 Ma, 1465 Ma, 1625 Ma, 1710 Ma, and 1850 Ma based on U-Pb detrital zircon age probability plots from the Bliss Sandstone in southern New Mexico [72]. The resulting inheritance envelope encompasses all but one ZHe date from the Bliss Sandstone (Figure 8(b)). These observations are consistent with the zircon grains in the Bliss Sandstone in the Franklin Mountains being derived from Proterozoic sources, similar to other regions of southern New Mexico. These results also suggest that the post-depositional t-T path for the Bliss Sandstone (red curve in Figure 8(b)) is a plausible solution for its thermal history.

The analysis above incorporates multiple methods to constrain a ZHe date-eU dataset. Southern New Mexico and western Texas have experienced a prolonged and complex tectonic history and, while the ZHe dataset may not be sensitive to all tectonic and sedimentation events, the datasets do constrain parts of the overall history with varying confidence. Below, we first comment on similarities and differences between the forward and inverse modeling methods. We then use the inverse modeling and detrital inheritance curve results to explore the more regional tectonic significance. These results are used to comment on multiple tectonic events in the study area spanning more than a billion years of Earth history.

Our modeling approach uses two methods for investigating thermal histories using ZHe date-eU relationships. The forward modeling approach was designed to incrementally narrow the possible thermal history by varying separate segments of the t-T path. This approach was successful at producing best-fit paths for each region that yield date-eU curves that are consistent with the observed ZHe data (Figure 6). However, because the entire thermal history of each region is considered (over a billion years of time), the best-fit paths are inevitably simplistic, allowing for the possibility that additional t-T paths might exist that are also consistent with the data. Our next attempt to constrain each thermal history was to use an inverse model approach to test for this possibility by investigating a total of 50,000 paths per sample, rather than tens of paths with the forward model approach.

In the Carrizo Mountains, the best-fit forward model and the inverse modeling paths are similar (Figure 7(a)), although the inverse modeling results allow for more complexity in the sample’s thermal history. For example, the inverse modeling paths suggest that the sample may have remained at shallow depths in the crust from 1000 to 600 Ma, followed by a second pulse of uplift. The inverse modeling paths also allow for a wider range of Paleozoic burial temperatures prior to Ancestral Rocky Mountain uplift. The largest difference between the two modeling results is during 80–40 Ma Laramide deformation. The forward model predicts significant Laramide uplift, followed by a smaller pulse of Rio Grande rift exhumation. The inverse model results, however, suggest that Laramide exhumation was minimal, and much of the final cooling was accomplished within the last 30 Ma (Figure 7(a)).

In the Franklin Mountains, the best-fit forward model and the inverse model paths are also similar for most of the thermal history (Figure 7(b)). The biggest inconsistency between the two is from 1000 to 500 Ma, where the best-fit forward model remains at surface temperatures, whereas the inverse model paths all show some amount of reheating prior to Paleozoic deposition. However, both approaches predict similar Mesozoic burial temperatures followed by more recent cooling to surface temperatures.

Similar to the other two locations, comparison of forward and inverse modeling results from the Cookes Range also shows the widest discrepancy during the Proterozoic (Figure 7(c)). Whereas our forward modeling analysis yielded multiple Proterozoic cooling paths that are consistent with the observed data, the inverse modeling results are most consistent with steady cooling from 1600 to 500 Ma, rather than earlier pulses of rapid cooling. The inverse models also suggest that the Proterozoic history of this region is more complex than the simplistic forward models, possibly involving multiple pulses of cooling separated by hundreds of millions of years. The Phanerozoic segment of the forward model is consistent with the inverse model results, although the inverse model paths allow for more variability.

Overall, the comparison between forward and inverse modeling results is similar at all three locations. However, in instances such as this where the goal is to constrain the entire thermal history from ZHe data, the forward modeling approach is overly simplistic and is unable to explore the full range of possible thermal histories. Thermal histories from these regions are complex and a thorough investigation is not plausible using a forward model approach where tens of paths are investigated. An inverse modeling approach that tests thousands of paths is better equipped to fully explore possible thermal histories in order to more confidently interpret the results. As a result, we will refer to the inverse model results during the remainder of the discussion below.

In the Carrizo Mountains, inverse modeling results combined with existing ⁴⁰Ar/³⁹Ar data [41] suggest rapid cooling to temperatures <150°C at 1030–1000 Ma (Figure 7(a)). This period of cooling is similar to ⁴⁰Ar/³⁹Ar cooling ages and likely reflects movement along the Streeruwitz thrust during continent-continent collision (Figure 1(b)). All but three of the paths remain at temperatures higher than approximately 75°C until 600 Ma. These results are consistent with sedimentological analyses of the synorogenic Hazel Formation in the footwall of the Streeruwitz thrust. Metavolcanic and metasedimentary clasts derived from the Carrizo Mountain Group are notably absent in the Hazel conglomerate, suggesting that these rocks were most likely not exhumed to the Earth’s surface at this time [73]. Instead, most inverse paths preserve a second pulse of cooling to near-surface temperatures from 600 to 500 Ma (Figure 7), which coincides with the final breakup of Rodinia (Figure 9). Although slightly less constrained, final exhumation in the Franklin Mountains and Cookes Range occurred within a similar timeframe. In the Franklin Mountains, the Red Bluff granite was emplaced at shallow depths and subsequently buried, consistent with geochemical studies of the Red Bluff granite that suggest it formed within an extensional stress field [39]. Good paths for the Franklin Mountains and Cookes Range both preserve a pulse of cooling to near-surface temperatures between approximately 800 and 500 Ma (Figures 7(b) and 7(c)). These results are similar to previous ZHe thermochronologic investigations. DeLucia et al. [14] presented ZHe data and inverse modeling results from the Ozark Plateau of Missouri that document pulses in exhumation from 850 to 680 Ma and 225 to 150 Ma, coinciding with the breakup of Rodinia and the breakup of Pangaea. They suggest a genetic relationship between continental uplift/denudation and the supercontinent cycle that is also observed in rocks from southern New Mexico and western Texas.

In southern New Mexico and western Texas, Ancestral Rocky Mountain intracontinental deformation resulted in a series of uplifts and sedimentary basins (e.g., [19]). The sedimentary record preserved in Ancestral Rocky Mountain basins suggests diachronous subsidence histories that young to the west, a pattern that is interpreted to reflect diachronous closure of the Ouachita suture [74]. Within the study area, the Pedregosa and Orogrande basins experienced similar times of peak subsidence rates from approximately 300 to 285 Ma [74]. Inverse models for the Carrizo Mountains and Cookes Range include Paleozoic constraint boxes that allow for a wide range of burial temperatures (50–300°C) from 450 to 290 Ma (Figure 7). Good paths are not sensitive to peak burial temperatures in either model other than they could not have been heated to temperatures >180°C. These paths are also not sensitive to the timing of reheating, with paths reaching peak temperatures from 450 to 300 Ma.

These observations suggest that, at least in this study area, ZHe datasets are not sensitive to either the timing or magnitude of Ancestral Rocky Mountain deformation (Figure 9). Instead, ZHe date-eU patterns are largely controlled by the timing and magnitude of both older and younger events that affected the thermal history of these rocks. These observations are useful for possible future studies. In order to investigate burial and exhumation related to the Ancestral Rocky Mountain event using ZHe thermochronology, we suggest either focusing on a location with a less complicated tectonic history so that the date-eU pattern is largely controlled by Ancestral Rocky Mountain deformation or combining ZHe data with additional datasets that constrain either the older or younger tectonic history.

Southwestern New Mexico, southeastern Arizona, and northern Sonora, Mexico experienced rifting starting in the Late Jurassic which created the Bisbee basin [48, 49]. The rift basin strata in southern New Mexico reached a maximum thickness of about 2800 m by Albian time [48]. Near the base of the section (TSA1 of Lawton et al., in press), these strata are interlayered with asthenospherically derived mafic volcanic rocks [52]. The strata on the flanks of the basin are much thinner or nonexistent, as the rift shoulder was a topographic high during basin development. The samples from this study were all along the inferred edge of the rift (Lawton et al., in press), but they still could have experienced the elevated geothermal gradients that are common in continental rifts (e.g., [75, 76]). In addition, the injection of mafic magmas into the upper crust in the areas near our study area could have been a factor in conducting heat from the mantle. Gradients in rifts can exceed 60°C/km [77] but can decrease rapidly after extension ends. These elevated geothermal gradients are a likely cause of the maximum temperatures of all three sample locations. In the Bisbee basin, rifting was followed by shortening resulting in the creation of a basin in Albian time dominated by dynamic subsidence in response to arc-continent collision in Mexico (Lawton et al., in press). Thus, geothermal gradients likely reverted to normal values by 100 Ma in the study area.

Southern New Mexico and western Texas are at the eastern limit of Laramide deformation. In southern New Mexico, NE-directed shortening produced a series of NW-trending basins and uplifts that range in age from Late Cretaceous to Eocene [22]. Most, if not all, Laramide faults are reactivated normal faults that developed during the formation of the Mexican Border rift [22, 78], and many thrust faults have fault throws that exceed several kilometers (e.g., [79]). Faults and folds related to Laramide shortening follow the Texas-Mexico border towards Big Bend National Park [80]. However, these structures do not extend to the Carrizo Mountains, where Laramide deformation consists of small-amplitude folds and minor thrust faults [81].

Laramide magmatism in southwestern New Mexico occurred in several pulses between about 75–70 Ma and 60–55 Ma, and 45–40 Ma [82, 83]. Evidence includes tuffs [82], andesite boulders in basins [84], plutonic rocks such as the Sylvanite complex [79], and various intrusions associated with ore deposits [85]. Our inverse modeling results are designed to include only exhumation and cooling during Laramide uplift, although the nonuniqueness of these models does not preclude some amount of Laramide reheating. In addition to the exposed volcanic and plutonic rocks, deeper Laramide plutons could also have affected ZHe dates. For example, Murray et al. [86] provide numerical models that suggest midcrustal plutons can reset low-temperature thermochronologic ages in upper crustal rocks. This is a result that is not explored here, but could also be a contributing factor in the observed ZHe data.

Cooling that is coeval with the Laramide orogeny is preserved in inverse models from the Cookes Range and the Franklin Mountains. In the Cookes Range, paths show a pulse of cooling beginning 60–70 Ma and ending by 45 Ma (Figure 9). The Franklins preserve a similar time of cooling from approximately 80 to 50 Ma. Although scattered Laramide basins in southern New Mexico were present during the Late Cretaceous, the majority of basins and uplifts developed largely during the Paleocene–Eocene [22], overlapping with the inverse models from the Cookes Range and Franklin Mountains. In contrast, the Carrizo Mountains show almost no cooling during this timeframe, and Cenozoic cooling to near-surface temperatures was delayed until after 30 Ma (Figure 9). These observations are consistent with the location of the Carrizo Mountains at the edge of the belt of Laramide shortening where total deformation was minimal.

The Neogene Rio Grande rift of southern New Mexico disrupted all previous tectonic elements and imparted the present topographic grain in the landscape. At the latitude of El Paso, Texas, the rift makes a sudden bend and structures south of here trend NW-SE instead of N-S as they do in the central and northern segments of the rift. In Colorado and New Mexico, combined apatite fission-track and apatite (U-Th)/He data suggest that extension was largely synchronous across this region, and these data record a main pulse of cooling from 25 to 10 Ma that is interpreted to reflect extensional exhumation (e.g., [28, 87–89]).

Inverse models from all three study locations preserve a pulse of cooling that is coeval with the development of the Rio Grande rift. Rapid cooling began at approximately 25 Ma in the Cookes Range, 25–15 Ma in the Franklin Mountains, and 30 Ma in the Carrizo Mountains (Figure 9). Although this timeframe is only a fraction of the entire thermal history investigated (over a billion years), it seems to have a significant effect on the resulting ZHe date-eU curve at each location. These results are consistent with previously published low-temperature thermochronologic data and modeling from the southern rift that suggest that the southern segment of the rift was active at the same time as the central and northern segments and that cooling of rocks was accomplished through fault exhumation rather than through magmatic injection [28, 87].

A total of 55 individual grain ZHe dates are presented from three ranges in southern New Mexico and western Texas. ZHe dates span hundreds of millions of years and record long-term thermal histories of Proterozoic crystalline rocks. These rocks have experienced a prolonged tectonic history involving multiple periods of cooling and reheating. A combination of forward and inverse modeling techniques was used to investigate plausible thermal histories of these samples. We find that although a forward modeling approach is advantageous for quickly comparing several dissimilar thermal histories, an inverse modeling approach can further refine these results and can test tens or hundreds of thousands of possible t-T paths.

Our thermal history inverse modeling results lead to the following conclusions: (1) the Red Bluff granite in the Franklin Mountains was likely buried and exhumed prior to deposition of Paleozoic sediments, and the unconformity at this location is therefore a compound erosional surface; (2) Proterozoic exhumation at all three locations occurred between 800 and 500 Ma, coinciding with the break-up of Rodinia; (3) all three study sites record elevated temperatures at approximately 100 Ma that likely reflects elevated heat flow during continental rifting to form the Mesozoic Bisbee basin. These results suggest that elevated heat flow within the rift extended into the rift shoulder, rather than being confined to the rift axis; (4) Laramide cooling occurred from 70 to 45 Ma in the Cookes Range and 80 to 50 Ma in the Franklin Mountains. In contrast, the Carrizo Mountains shows no cooling during Laramide shortening, which likely reflects its location at the edge of observed Laramide deformation; and (5) all three study sites preserve a pulse of cooling that begins at 30–25 Ma, simultaneous with observed cooling in the central and northern segments of the rift.

In contrast, these data provide almost no information on the timing or magnitude of Ancestral Rocky Mountain deformation. These data provide important information on the broad thermal history of southern New Mexico and western Texas and highlight the effectiveness of using the ZHe method in tectonically complex regions. These data also help to bridge the gap between lower temperature systems such as apatite (U-Th)/He and fission-track and higher temperature methods of thermochronology such as ⁴⁰Ar/³⁹Ar.

The data supporting the results of this study include zircon (U-Th)/He dates, grain sizes, and U, Th, Sm, and He concentrations. This information, along with sample locations, is included in Table 1.

The authors declare that they have no conflicts of interest.

JWR was supported by NSF-EAR 1624538 and JMA was supported by NSF-EAR 1624575. ZHe data were obtained from the (U-Th)/He thermochronology lab at CU Boulder. We thank Rebecca Flowers and James Metcalf for their help with data acquisition and discussions and Michelle Gavel for assistance with mineral preparation. Two anonymous reviewers provided very useful feedback that helped improve the quality of this manuscript.