Zircon (U-Th)/He thermochronology of Grand Canyon resolves 1250 Ma unroofing at the Great Unconformity and <20 Ma canyon carving

The Precambrian-Cambrian crystalline basement–sedimentary contact, or the Great Unconformity (GU), represents a prominent erosion surface observable across North America (Marshak et al., 2017) and globally (Ronov et al., 1980). Researchers have debated whether the erosion processes that created the GU are related to the Rodinia supercontinent cycle (Timmons et al., 2005; DeLucia et al., 2018; Flowers et al., 2020) or Cryogenian glaciations (Keller et al., 2019). The timing of these weathering processes has broad implications for resolving the GU's association with key changes in oceanic and atmospheric chemistry (Husson and Peters, 2017), and subsequent biomineralization during the Cambrian explosion (Peters and Gaines, 2012). The current debate as to the formation of the GU therefore centers on the timing of basement erosion, whether erosion was globally synchronous, and the composite or singular nature of the GU erosion surface.

The Grand Canyon (southwestern United States; Fig. 1) represents an important geologic laboratory where the GU was first defined (Dutton, 1882) and is exceptionally well exposed. Fault blocks of the Grand Canyon Supergroup show that the GU is a composite erosion surface that records exhumation of Proterozoic crystalline basement during several pulses of intracratonic tectonism (Karlstrom and Timmons, 2012). The amount of time “missing”, represented by each of the unconformities, is shown in the inset in Figure 1 (Karlstrom et al., 2020a; Table S1 in the Supplemental Material¹). This highlights that the duration of the GU varies depending on the ages of local strata (e.g., the sub-Tonto Group unconformity has between 200 and 1300 Ma of time “missing”).

We applied the zircon (U-Th)/He (ZHe) thermochronometer to model the thermal histories of basement rocks of the eastern Grand Canyon. Our models incorporate constraints from previously published apatite (U-Th)/He (AHe) and fission-track (AFT) data (see a summary in Karlstrom et al. [2020b]), and hornblende, mica, and potassium feldspar (K-spar) ⁴⁰Ar/³⁹Ar results (Timmons et al., 2005). Our goals were to (1) test congruence among ⁴⁰Ar/³⁹Ar, ZHe, AHe, and AFT data; (2) use multiple inverse modeling approaches to assess the resolving power of ZHe data as applied to deep-time thermal histories; and (3) examine implications at both ends of the timescale to address debates about the GU's singular versus composite nature, as well as the age of carving of the Grand Canyon (c.f. Flowers and Farley, 2012, 2013; Karlstrom et al. 2014, 2020b; Winn et al., 2017).

The ZHe data set is composed of 37 single zircon grain aliquots from 7 sample locations (∼5 grains from each sample) that span river miles (RM) 80–236 (Fig. 1; measured downstream from Lees Ferry; see Table S1 and the Supplemental Material text for analytical methods). Six samples from eastern Grand Canyon basement rocks ranging in age from 1.71 to 1.84 Ga (zircon U/Pb ages) were collected between RM 80 and 115 (Fig. 1). These samples were collected from approximately the same depths below the Unkar and Tonto Groups (Fig. 1 inset), had a demonstrably similar post–1.4 Ga thermal history, and are the focus of our modeling. They are compared with a single 1.1 Ga (zircon U/Pb ages) diabase sample from the Grand Canyon Supergroup that had a different early cooling history, and one 1.68 Ga (zircon U/Pb age) basement sample from the western Grand Canyon (sample K12-236-6).

ZHe dates show large intra- and inter-sample age variations that ranges from ca. 3 Ma to ca. 543 Ma (Table S1). Figure 2 shows the ZHe data set plotted in date-effective eU (eU = U + 0.235 × Th) space. Differences in eU concentration among zircon crystals with the same time-temperature (t-T) history are proportional to differences in radiation damage, which results in a range of crystal-specific He diffusivities due to the damage-diffusivity relationship (Guenthner et al., 2013). Individual samples show modest variation in eU concentration, but large variation (108–2131 ppm) for the whole data set. Older ages correspond to lower eU (e.g., 542 Ma at 108 ppm) and younger ages to higher eU (e.g., dates of 3–23 Ma at 1282–1477 ppm; Fig. 2). An inflection in the date-eU trend of eastern Grand Canyon ZHe dates occurs at ca. 20 Ma and all points between ∼725 and 2200 ppm eU are younger than 23 Ma, with 7 out of 10 grains in the range of 3–9 Ma (Fig. 2).

We applied several modeling approaches that model the date-eU response to different thermal histories in the context of the damage accumulation and/or annealing and He diffusion relationship in zircon (Guenthner et al., 2013). Our focus is on the basement rocks in the eastern Grand Canyon, RM 80–115, as they probably share a similar low-temperature (<350 °C) thermal history post 1.4 Ga (Timmons et al., 2005; Karlstrom et al., 2020b). Forward thermal history modeling is discussed in the Supplemental Material, as those results are complementary to our inverse modeling. We applied inverse modeling using two different software packages: HeFTy (version 1.9.3; Ketcham, 2005) and QTQt (Gallagher, 2012; http://www.iearth.org.au/codes/QTQt/). HeFTy was used to assess compatibility among ZHe and other available thermal history constraints: multi-mineral ⁴⁰Ar/³⁹Ar data, AHe and AFT data, and the stratigraphic constraints of Unkar and Chuar Group deposition (Fig. 3; see Table S2 for constraint box details). We used QTQt to explore the key regions of t-T space that the ZHe data alone are most likely to predict when left unconstrained.

HeFTy damage models are limited to testing dispersion that can be quantified by the damage-diffusivity relationship alone. Our data set shows date dispersion that exceeds the first-order trend between date and eU, and could result from the influence of U and Th zonation, which can be large and lacks a routine approach for quantification (Guenthner et al., 2013). Previous researchers have modeled ZHe and AHe data with similar second-order dispersion through a synthetic grain approach (Flowers and Kelley, 2011; Murray et al., 2016; DeLucia et al., 2018; Flowers et al., 2020), with dates binned based upon eU concentrations, and mean values used for the date (±1σ standard deviation), eU concentration, and grain size. For our ZHe inputs, we combined the eastern samples (the 1.1 Ga diabase sample excluded) because these zircon grains share a common post–1700 Ma thermal history, and we split our grains into three bins: 200–400 ppm, 400–800 ppm, and >1200 ppm (there were no grains between 800 and 1200 ppm). For AHe and AFT inputs, we selected previously published samples that spanned the same river miles as our eastern Grand Canyon samples, and that had results for both thermochronometers: 98GC20 and 98GC11 (Lee et al., 2013). These inputs are provided in Tables S2 and S3.

Figure 3 shows HeFTY model output for our ZHe inputs, AHe and AFT inputs from thermochronometer 98GC11, and constraint boxes as annotated in the figure and detailed in Table S1. Outcomes for both 98GC11 and 98GC20 (Fig. S2) are similar. From 200,000 attempted paths, the model produced 5 “good” (goodness of fit [GOF] > 0.5) and 3402 “acceptable” (GOF > 0.05) t-T paths, which shows to first order that the ZHe data are compatible with known constraints. The number of good paths is small due to the number and size of constraint boxes, and the demands of using three different thermochronologic systems. These paths support rapid cooling prior to 1255 Ma in agreement with the geology and K-spar multi-diffusion domain (MDD) models (Timmons et al., 2005). Good paths suggest that maximum Laramide reheating as high as 140 °C was achieved by ca. 80 Ma, followed by rapid cooling to ∼90 °C by 55 Ma, slow cooling from 55 to 20 Ma, and rapid cooling to near surface temperatures after 20–10 Ma, which agrees with previous AHe and AFT thermal history modeling (Lee et al., 2013).

Our HeFTy models are purposely constrained to test compatibility of the ZHe data with all known geologic and thermochronologic Grand Canyon constraints. In contrast, we used QTQt to explore the sensitivity of the ZHe data alone with a single constraint box (1750 ± 100 Ma, 750 ± 50 °C) that represents the range of U/Pb dates of the basement samples. All single-grain ZHe data were used, and, due to the lack of constraint boxes, it is expected that the QTQt-generated t-T paths will traverse regions of t-T space not in agreement with geological or other thermochronologic constraints. The goal was therefore to assess the QTQt results in terms of sensitivity of the ZHe data to the most likely thermal scenarios, and explore the key t-T segments that have the most leverage for a date-eU curve. The post-crystallization thermal history produced by the QTQt model shows that the ZHe data is sensitive to cooling from ∼750–200 °C from ca. 1750 to 1200 Ma (Fig. 4). Obvious disagreements with the observed geologic record are apparent for the post–1200 Ma thermal history, which we attribute to a lack of sensitivity of the ZHe data to identify Neoproterozoic to Devonian unconformities. Importantly, the late Phanerozoic portion (Fig. 4 inset) shows Cenozoic reheating much later than known maximum Laramide burial conditions at ca. 80 Ma but then shows rapid cooling at ca. 7 Ma.

The contrast between HeFTy and QTQt models highlights the promise of ZHe data alone, but also the necessity of combining multiple geologic and thermochronologic observations. QTQt paths are most likely to have early (1800–1200 Ma) and late (7–0 Ma) cooling in broad agreement with other thermochronometers. But limitations of deep-time thermochronology based solely on the ZHe method are apparent for the time period from 1200 Ma to 300 Ma, when geologic constraints indicate basement rocks were within a few kilometers of the surface during Unkar and Chuar Group deposition but ZHe data alone provide no constraints. HeFTy models show that even with additional stratigraphic constraints, a wide range of t-T paths are plausible, and hence ZHe is relatively insensitive to deposition and unroofing of the Unkar, Chuar, and Tonto Groups (see also forward models in the Supplemental Material). Insets for both models (Figs. 3 and 4) show that the ZHe data are sensitive to late cooling caused by canyon carving; however, the trajectory of this cooling depends upon the inputs used. HeFTy produces a multi-stage cooling path with a second pulse of major cooling from ∼80 °C to near-surface temperatures beginning at 20 Ma, whereas QTQt shows a single pulse of cooling at ca. 7 Ma, matching the numerous young dates (3–9 Ma) over a broad range of eU concentrations (∼700 to ∼1500 ppm).

Deep-time ZHe thermal histories from basement rocks within a few hundred meters below the GU surface in the Grand Canyon show a dominant unroofing event in the eastern Grand Canyon at 1350–1250 Ma. ZHe inverse models for the eastern Grand Canyon all require a Mesoproterozoic phase of cooling to predict the first-order trend between older ages and lower eU concentrations between ∼50 ppm and 500 ppm. Peak et al. (2021) also reported ZHe data from basement rocks of the Grand Canyon that show overlap and similar scatter to our data in both the western and eastern Grand Canyon. Our data differ by having more <25 Ma ages at 700–2000 ppm eU, which we purposefully picked (i.e., Ault et al., 2018). These young ages cause our thermal models to be sensitive to both Precambrian and Cenozoic t-T segments, and this may explain the different interpretations from similar data. In addition, Peak et al. (2021) interpreted differences in supergroup deposition between the eastern and western Grand Canyon; however, given our model results, we argue that ZHe data (theirs and ours) alone are relatively insensitive to 1–2-km-scale Mesoproterozoic and Neoproterozoic faulting (Karlstrom and Timmons, 2012) and do not provide evidence for an absence of Grand Canyon Supergroup deposition in the western Grand Canyon region as proposed by Peak et al. (2021).

A major pulse of Mesoproterozoic basement cooling and unroofing, as shown by the ⁴⁰Ar-³⁹Ar data in the western Grand Canyon region (Karlstrom, et al., 2010) and our combined data in the eastern Grand Canyon, is consistent with regional uplift due to far-field stresses related to Grenville-age orogenesis (Karlstrom et al., 2001; Timmons et al., 2005, Mulder et al., 2017). At the local scale of the eastern Grand Canyon, these stresses produced rock exhumation (i.e., cooling) in monoclines and up-thrown footwall blocks associated with Proterozoic normal faults (Timmons et al. 2005). Such events are broadly related to Rodinia assembly (ca. 1300–900 Ma), but notably predate Rodinia rifting (ca. 800–550 Ma), a mechanism suggested for exhumation of the GU in other North American locations constrained by ZHe data (DeLucia et al., 2018; Flowers et al., 2020). Our data are not compatible with a major pulse of snowball Earth–related erosion at this location (Keller et al., 2019). However, given the contrast between modeling results presented here and by Ricketts et al. (2021) that show a Mesoproterozoic cooling pulse (1350–1250 Ma) in the Rocky Mountains, and the Neoproterozoic signal of cooling (850–680 Ma) from the cratonic interior (DeLucia et al., 2018), a broader deep-time thermochronologic transect across Laurentia is needed to fully understand the multiple mechanisms for the profound multi-stage basement exhumation that created the GU composite erosion surface.

A key methodological finding of our study is that ZHe modeling is highly sensitive to both early and late components of the thermal history and can resolve aspects of the Cenozoic canyon evolution. Our joint inversion of the ZHe, AFT, and AHe data shows peak pre-Laramide temperatures of 140 °C, and our new ZHe results agree with the previously reported eastern Grand Canyon AHe and AFT data (Flowers et al., 2008; Lee et al., 2013) as interpreted by Karlstrom et al., (2014, 2020b). Our models and data support multi-stage incision of the eastern Grand Canyon at 20–10 Ma, and after 7 Ma and do not support rapid cooling to surface temperatures at 70—50 Ma as proposed in “old” Grand Canyon models (Wernicke, 2011; Flowers and Farley, 2012, 2013). This result is seen in two different types of inverse models, which show that the ZHe method alone can help resolve multiple cycles of cooling and reheating across the geologic timescale, with improved precision coming from the integration of ZHe with multiple thermochronometers.

W. Guenthner acknowledges support from U.S. National Science Foundation grants NSF-EAR1848013 and 1735788. K. Karlstrom and M. Heizler acknowledge NSF-EAR grant 9902955. K. Karlstrom acknowledges NSF-EAR grants 1348007 and 1955078. We thank four anonymous reviewers for extensive comments that improved this paper.