Zircon (U-Th)/He thermochronology reveals pre-Great Unconformity paleotopography in the Grand Canyon region, USA

The Great Unconformity is exposed along the length of the Grand Canyon in northwestern Arizona, United States (Fig. 1) and separates the Cambrian Tonto Group from the underlying Paleoproterozoic basement or Mesoproterozoic-Neoproterozoic Grand Canyon Supergroup. It represents as much as 1.2 b.y. of missing time (Timmons and Karlstrom, 2012). Recent studies have identified various events potentially associated with the Great Unconformity erosion surface that include >800 Ma Rodinia amalgamation, ca. 800 Ma early Rodinia breakup, 717–635 Ma Cryogenian Snowball glaciations, and ca. 580–500 Ma late Rodinia breakup and the Pan-African Orogeny (e.g., DeLucia et al., 2018; Keller et al., 2019; Flowers et al., 2020). Evidence of erosion during all of these periods is preserved in the Grand Canyon Supergroup of the Upper Granite Gorge (UGG; Fig. 1C); in unconformities within the Unkar Group (>800 Ma), disconformities between the Cardenas Basalt, Nankoweap Formation and the Chuar Group (ca. 800 Ma), and the unconformity separating the Chuar Group and Sixtymile Formation/Tapeats Sandstone (spanning ca. 730–520 Ma; Karlstrom et al., 2020). The Lower Granite Gorge (LGG) does not preserve the Grand Canyon Supergroup, which makes it unclear whether the LGG and UGG share a common Neoproterozoic history. Together, these geologic relationships suggest a multiphase and possibly spatially variable history of Great Unconformity development. Here we present ZHe data to decipher the origin of this feature in its iconic Grand Canyon exposure.

The UGG and LGG of the Grand Canyon expose 1.8–1.4 Ga basement, which remained at depths consistent with temperatures >400 °C (∼12–15 km) until ca. 1.4 Ga (Williams and Karlstrom, 1996; Dumond et al., 2007). In the UGG, the Proterozoic Grand Canyon Supergroup occurs on top of basement, and the full Supergroup and Sixtymile Formation (∼3 km thick in total) are only preserved in the easternmost part of the gorge (Fig. 1). The region is cut by faults that offset the basement and Supergroup (Timmons et al., 2005), but only small offsets are apparent in the Phanerozoic units, which indicates that Precambrian tectonism is responsible for most of the observed displacement. In the LGG, the Great Unconformity is defined by Tonto Group Tapeats Sandstone overlying basement, whereas in the UGG, ca. 1255 Ma, Unkar Group rests on basement. It is unclear whether the Supergroup originally extended over the LGG and was largely removed by the sub-Tapeats unconformity or if the unconformity in the LGG is a composite surface with the Tapeats capping older topography. Previous studies have suggested that the Chuar basin was restricted in mid-Chuar time from the proposed Tonian intracontinental seaway (e.g., Dehler et al., 2017; Rooney et al., 2017). This restriction could have been caused by paleotopography. Throughout the Grand Canyon, the Tapeats is succeeded by Paleozoic strata with an Ordovician-Devonian hiatus. These units were buried by Mesozoic foreland deposits that were later removed (DeCelles, 2004). Previous apatite fission-track and apatite (U-Th)/He data document Phanerozoic burial temperatures >80 °C for river-level samples and help constrain subsequent erosion history (e.g., Dumitru et al., 1994; Flowers et al., 2008; Flowers and Farley, 2012; Lee et al., 2013; Winn et al., 2017).

Rocks cool as they are exhumed, and this cooling history—and by proxy, exhumation history—can be recorded by ZHe thermochronology (e.g., Reiners et al., 2002). This method exploits the radioactive decay of U and Th to He. At temperatures >220 °C, He will diffuse completely out of a zircon crystal; at lower temperatures, the He will be retained. The exact temperature-diffusion relationship varies due to radiation damage, which accumulates and anneals with tim as a function of temperature (Guenthner et al., 2013; Ginster et al., 2019). Damage is proxied by effective uranium concentration (eU) for a zircon suite that underwent the same thermal history, or by α-dose estimates. With increasing eU, or α-dose up to ∼1 × 10¹⁸, zircon becomes more He retentive, but at higher damage the He retentivity decreases. This can cause positive and negative date-eU correlations at low and high damage, respectively. Thermal histories to explain a given ZHe data set can be explored using radiation damage accumulation and annealing models for He diffusion, which can include various damage annealing kinetics (Guenthner, 2021). Other factors can affect the (U-Th)/He date and include α-ejection, He implantation, inclusions, eU zonation, and grain size. With appropriate information, some of these effects can be corrected for or avoided (see the Supplemental Material¹ for more detail). Especially important to this study is eU zonation. (U-Th)/He dates for zoned grains may differ from their unzoned counterparts with the same bulk eU. Variability in zonation patterns between grains can introduce dispersion into date-eU relationships, and these effects are magnified by small grain size (e.g., Hourigan et al., 2005; Farley et al., 2011; Ault and Flowers, 2012).

We acquired ZHe data for four samples each from the LGG and UGG (Tables S1 and S2 in the Supplemental Material). Seven of these samples are Precambrian granitoid basement collected near river level, and one is the 729 ± 0.9 Ma Walcott Member Tuff near the top of the Chuar Group (Fig. 1D). To better understand the effects of eU zonation on ZHe dates and their interpretation, we obtained single U, Th, and Sm concentration profiles for 7–8 zircon grains per basement sample using depth-profiling by laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) (Fig. S2). Zonation data were not acquired for the tuff sample because all zircon of sufficient size were dated before zonation analysis. See the Supplemental Material for details.

The LGG ZHe data fall on a single negative date-eU trend spanning 740 ± 27 Ma to 69 ± 4 Ma (Fig. 2A). There is no correlation between date and grain radius (Fig. S1A). Most zircon zonation profiles for these samples have rims enriched in parent nuclides relative to cores and there is limited intrasample variability in eU zonation patterns (Figs. S2 and S3).

ZHe data patterns vary among the UGG samples (Fig. 2B). Samples CP06–52 and UG90–2 yield low eU zircon with maximum dates >700 Ma and lack obvious date-eU correlations. In contrast, despite zircon with comparably low eU, the other UGG samples (UG96–1 and EGC1) yield ZHe dates all <400 Ma with one exhibiting a negative date-eU trend and the other a positive trend. As with the LGG samples, there is no apparent relationship between ZHe date and grain radius (Fig. S1B), and most zircon have rims higher in eU than cores (Figs. S2 and S3). Sample UG90–2 shows high intrasample variation in its eU zonation pattern (Fig. S2E), which may help explain its substantial ZHe date-eU scatter.

ZHe data for the LGG and UGG document differences in basement thermal history. In the LGG, the Neoproterozoic dates record a portion of the Proterozoic time-temperature (t-T) path, and the consistency in the date-eU pattern across samples suggests a shared thermal history (Fig. 2A). In the UGG, data patterns differ from those in the LGG (Fig. 2B), which implies differences in t-T paths across the ∼100 km separating these sample suites. In addition, inter-sample variability in the UGG data, with spatially alternating basement samples with low eU zircon that yield either Neoproterozoic results or much younger ZHe dates, points to more abrupt differences in t-T paths at the tens of kilometers scale. Moreover, Chuar sample EGC1 is stratigraphically higher and younger (729 Ma) than the other samples in the UGG (1.7 Ga; Fig. 1D) but yields post-729 Ma ZHe dates, which also suggests differing thermal histories across short spatial wavelengths.

Broad uniformity in Phanerozoic sedimentary thickness and resultant burial heating across the region implies that the spatial differences in thermal history recorded by ZHe must date to the Proterozoic. Paleozoic sedimentary rocks across the region thicken slightly westward (Beus and Morales, 2003; Timmons and Karlstrom, 2012), and Mesozoic burial thickened eastward (Robinson Roberts and Kirschbaum, 1995; DeCelles, 2004; Wernicke, 2011) but both over spatial wavelengths too large to explain the variability in ZHe data patterns. Instead, we suggest that variable Neoproterozoic burial and exhumation histories across small spatial scales induced by Neoproterozoic faulting during deposition of the Grand Canyon Supergroup is the most likely explanation for the data set.

In the Neoproterozoic, grabens and half-grabens offsetting the basement and Supergroup (Timmons et al., 2005) created conditions for disparate mid-late Proterozoic burial and exhumation histories across major faults. Fault systems in the UGG (Fig. 1B) were activated multiple times during the Proterozoic and culminated in normal faults during the Neoproterozoic based on observations such as offsets in basement and Supergroup-equivalent units to the north of the Grand Canyon, reverse offsets within the Unkar Group units, and reconstruction of pre-Laramide extensional offsets (Shoemaker et al., 1978; Timmons et al., 2001, 2005; Beus and Morales, 2003). The Chuar Syncline and bounding Butte Fault in the eastern UGG (Fig. 1B) were active during the Tonian as documented by stratigraphic thinning and were reactivated in the late Neoproterozoic to Early Cambrian as indicated by incision below the Cambrian Sixtymile Formation (Elston and McKee, 1982; Timmons et al., 2001; Karlstrom et al., 2020).

In the LGG, the absence of Grand Canyon Supergroup suggests that Proterozoic deposition may have been restricted to the UGG east of the Sinyala Fault System. To test this hypothesis, we carried out inverse thermal history simulations of the LGG data using the HeFTy software package (Ketcham, 2005) and the ZRDAAM model (Guenthner et al., 2013) for two endmember t-T histories (Fig. 3): (1) the Supergroup hypothesis (SG), applying the Supergroup burial and exhumation history as preserved in the eastern UGG (Elston and McKee, 1982; Timmons et al., 2005; Dehler et al., 2017; Rooney et al., 2017; Karlstrom et al., 2018), and (2) the Neoproterozoic exhumation hypothesis (NeoExh), in which the LGG was exhumed synchronously with Supergroup deposition in the UGG. Exhumation begins at 823 ± 26 Ma and represents the likely onset of normal faulting that accommodated the Chuar Group, as dated by K-Ar in the UGG (Elston and McKee, 1982), and is consistent with ca. 782 Ma detrital zircon in the base of the Chuar Group (Dehler et al., 2017). Phanerozoic constraints are the same in both models. LGG samples were modeled together (Table S3), and representative eU zonation profiles for each sample were used (Fig. 3D; Table S4). The HeFTy implementation of the widely used ZRDAAM model with fission-track annealing kinetics enables inclusion of zonation profile inputs, so modeling was done using this approach to honor this complexity. Model details are provided in the Supplemental Material and Tables S3–S7. The NeoExh model yielded t-T paths with better fits to the data than the SG model (Fig. 3A). This remains true when endmember combinations of grain size and observed zonation profile are used (Fig. 3C). These outcomes imply that of the two hypotheses tested, the NeoExh model is most consistent with the LGG ZHe data, compatible with the preserved Supergroup extent.

In the UGG, the spatial heterogeneity in data patterns suggests variability in the timing and/or magnitude of Proterozoic burial and exhumation across normal faults. To test this, we performed t-T forward and inverse models of several UGG samples (see the Supplemental Material text and Tables S8–S12). The outcomes illustrate that differences in the Proterozoic thermal history are required to explain the UGG data if the same Phanerozoic thermal history is assumed (Figs. S4 and S5). This is consistent with Neoproterozoic fault-induced variability in Supergroup burial, as also implied by the ZHe data patterns and preserved geologic constraints.

We interpret the different thermal histories of the LGG and UGG and within the UGG as caused by late Meso-Neoproterozoic faulting that produced paleotopography and syntectonic deposition and erosion. The “Upper” and “Lower” basins were likely separated by a paleo high bounded on either side by fault systems as is suggested by west-dipping normal faults between the LGG and UGG (Fig. 1A) and supported by inverse t-T modeling. Variation in thermal history among UGG samples can be explained by relationships to paleotopographical features inferred from preserved geology (Fig. 1): UG96–1 was in a paleo low in the hanging wall of the Crystal Fault (Timmons et al., 2001), where it underwent greater Neoproterozoic burial and associated He loss and now yields younger ZHe dates than sample CP06–52, which was located on a Neoproterozoic paleo high on the footwall of a normal fault (Timmons et al., 2001).

Our study outcomes are consistent with multiphase faulting and erosion in the Grand Canyon region contributing to Great Unconformity development over a protracted Proterozoic interval. Figure 4 shows our schematic reconstruction of the deposition, erosion, faulting, and paleotopographic history, with each time slice corresponding to a known faulting period. The Unkar Group was deposited in a fault-bounded basin at ca. 1255 Ma (Fig. 4A), and syn-depositional tectonic activity continued through ca. 1100 Ma (Fig. 4B). After Unkar deposition, normal faulting and erosion occurred at ca. 830–800 Ma with the onset of Chuar Group deposition at ca. 780 Ma (Elston and McKee, 1982; Dehler et al., 2017), while normal-fault exhumation of the present-day LGG began simultaneously (Fig. 4C). This geometry may have isolated deposition of the Grand Canyon Supergroup from areas farther west and thus from the proposed Tonian intracontinental seaway (Dehler et al., 2017). Tonian sedimentary rocks were deposited syn-tectonically in the deepening Chuar Syncline with shallower burial elsewhere in the UGG. The final pre-Tonto Group tectonic and erosion event occurred at ca. 520–510 Ma (Fig. 4D). This model proposes that the Great Unconformity in the Grand Canyon developed via multiple erosional events driven by tectonism differing on the scale of tens of kilometers, which indicates that small scale topography played an important role in erosion and deposition during the protracted breakup of the Rodinian supercontinent.

This work was supported by U.S. National Science Foundation grants EAR-1822119 and EAR-1916698 to R. Flowers and F. Macdonald and a University of Colorado–Boulder Chancellor's Fellowship to B. Peak. We thank Karl Karlstrom for organizing Grand Canyon trips resulting in sample archives. We thank Jim Metcalf for help with ZHe data acquisition, Emmy Smith for locating archived separates, and Mark Pecha at Arizona LaserChron (Tucson, Arizona, USA) for providing reconnaissance data. We thank David Foster and two anonymous reviewers for feedback that improved this manuscript.