Thermochronologic constraints on the origin of the Great Unconformity.

The origin of the phenomenon known as the Great Unconformity has been a fundamental yet unresolved problem in the geosciences for over a century. Recent hypotheses advocate either global continental exhumation averaging 3 to 5 km during Cryogenian (717 to 635 Ma) snowball Earth glaciations or, alternatively, diachronous episodic exhumation throughout the Neoproterozoic (1,000 to 540 Ma) due to plate tectonic reorganization from supercontinent assembly and breakup. To test these hypotheses, the temporal patterns of Neoproterozoic thermal histories were evaluated for four North American locations using previously published medium- to low-temperature thermochronology and geologic information. We present inverse time-temperature simulations within a Bayesian modeling framework that record a consistent signal of relatively rapid, high-magnitude cooling of ∼120 to 200 °C interpreted as erosional exhumation of upper crustal basement during the Cryogenian. These models imply widespread, synchronous cooling consistent with at least ∼3 to 5 km of unroofing during snowball Earth glaciations, but also demonstrate that plate tectonic drivers, with the potential to cause both exhumation and burial, may have significantly influenced the thermal history in regions that were undergoing deformation concomitant with glaciation. In the cratonic interior, however, glaciation remains the only plausible mechanism that satisfies the required timing, magnitude, and broad spatial pattern of continental erosion revealed by our thermochronological inversions. To obtain a full picture of the extent and synchroneity of such erosional exhumation, studies on stable cratonic crust below the Great Unconformity must be repeated on all continents.

One of the most profound divides in Earth history may be the one that separates rocks containing abundant macroscopic fossils from those that do not, a dividing line that is implicit in the name of Earth’s current geological eon—the eon of visible life, the Phanerozoic. For nearly as long as the significance of this dividing line has been appreciated, and before the name Phanerozoic was yet coined (1), it has been associated with another phenomenon—the frequent occurrence of one or more significant unconformities below the oldest rocks containing abundant macroscopic fossils (2). This phenomenon, taking its name from a particularly charismatic occurrence at the Grand Canyon (3), has subsequently been referred to by some authors as the Great Unconformity (4, 5). While lacunae in the geologic record are common (6), those below the oldest rocks of the Phanerozoic are frequently large—in many cases even juxtaposing undeformed sedimentary rocks above, with crystalline igneous or metamorphic basement below (4). The presence of the Great Unconformity in the rock record is significant because the erosion required to create the unconformity and the widespread burial that preserved it are both equally important. The crucial defining feature of the Great Unconformity is that erosion occurred across a vast area, especially the cratonic interior. The most quantitative reflection of this feature is arguably provided by the coeval stepwise increase in preserved sediment abundance per unit time across the unconformity—a step change first accurately quantified by Ronov et al. (7, 8) and Ronov (9) and observed on every continent with the possible exception of Africa. This fivefold discontinuity in global preserved sediment abundance (8) suggests profound changes in both erosional and depositional processes (5) and in any event provides a quantitative metric for the significance of the Great Unconformity as a global feature.

The Great Unconformity is, however, far from the only significant phenomenon associated with the emergence of the Phanerozoic world. The transitional Neoproterozoic era showed several significant changes in the Earth system, including the gradual breakup of the supercontinent Rodinia from ∼825 Ma to ∼570 Ma (10–12), possibly significant fluctuations in atmospheric oxygen (13), and two severe failures of Earth’s silicate weathering feedback (14) within the Cryogenian Period (717 to 635 Ma) that glaciated the continents down to the equator (15, 16). This interval culminated in the Ediacaran Period (635 to 540 Ma) when the appearance of a more diverse biosphere (17), especially macroscopic multicellular organisms, set the stage for the dramatic diversification of visible metazoan life in the earliest Cambrian (18–20). Perhaps the most marked and nonuniformitarian of these events were the hypothesized low-latitude glaciations (15, 16, 21). Maximization of silicate weathering sensitivity due to concentration of Rodinian continents near the equator favored ice-house conditions, and glaciation is thought to have initiated when sea ice advanced within ∼40 of the equator (22). The sea-ice/albedo positive feedback overwhelmed the silicate weathering negative feedback and continental glaciations extended to low latitudes in three episodes, the Sturtian (717 to 659 Ma), Marinoan (641 to 635), and Gaskiers (∼580 Ma)—of which the Sturtian and Marinoan were global “snowball Earth” events (15, 16, 22–24). The proximal trigger for all three glaciations, however, remains a matter of debate (25–27).

Recently, Keller et al. (5) proposed that widespread erosion by continental ice sheets during these Neoproterozoic glacial intervals may be responsible for the anomalous concentration of unconformities at the end of the Precambrian. If correct, a link between continental glaciation and kilometer-scale cratonic exhumation would have dramatic implications for our understanding of the long-term preservation, composition (via increased sediment subduction and relamination), and freeboard of continental lithosphere and could help explain a wide set of puzzling observations across several related disciplines (5)—including the much-discussed increase in apparent high-latitude terrigenous sediment flux coincident with Laurentide glaciation (28). However, this proposal has not been without controversy (29–31). While some of this controversy may be attributable to differences in terminology, significant points of contention remain—primarily, whether Neoproterozoic glaciation did or did not cause significant upper crustal exhumation. Resolving these differences is critical to our understanding of the Neoproterozoic Earth system and the couplings and feedbacks between tectonic, climatic, and biogeochemical processes therein.

Over the past century, the term “Great Unconformity” has acquired multiple overloaded meanings. Historically, the term was first applied by Clarence Dutton (3) from the rim of the Grand Canyon (United States) to the unconformity at the base of the flat-lying Phanerozoic sedimentary succession within (in some regions a disconformity and in other regions a nonconformity)—although, at the time, Dutton did not yet know the true age of the rocks involved. Subsequently, it has been variously used to denote one or more of the following:

An unconformity at or near the base of the Phanerozoic that is separating rocks that contain visible fossils from those that do not (2), either in general or at a specific locality;

A basement nonconformity, either in general or at a specific locality, often with the (perhaps implicit) additional requirement that the involved basement be Precambrian in age (29, 32); or

A broader phenomenon evidenced qualitatively by the observation (2, 4) that unconformity 1 and nonconformity 2 frequently coincide (especially relative to what one might expect by chance), suggesting the existence of a globally widespread exposure surface (4, 5)—an inference quantitatively confirmed by the global step in preserved sediment abundance first observed by Ronov et al. (8).

This variation in meaning invites confusion and controversy as to the synchroneity or diachroneity of the Great Unconformity, depending on which (or which combination) of the above meanings is intended. On the one hand, individual physical unconformity surfaces are ubiquitously composite in origin, with later episodes of erosion capturing and subsuming previous erosional surfaces. On the other hand, the set of unconformities spanning the base of the Phanerozoic (i.e., definition 1) are in a sense synchronous by nature, as is consequently, to some extent, the broader phenomenon implied in definition 3. Thus, we apply the term Great Unconformity to the temporal correlation of unconformities in the late Precambrian (definition 3), whereas, for example, the usage in Flowers et al. (29) is more aligned with definition 2—asserting diachronous worldwide development of many Great Unconformities in the Neoproterozoic.

After accounting for such semantic differences, remaining points of disagreement center on the question of whether or not Neoproterozoic glaciations were significantly erosive. Relatedly, while in no means mutually exclusive with glacial erosion, it also remains entirely worthwhile to quantify the relative contributions to Neoproterozoic crustal exhumation in different regions from such known tectonic events as Rodinia assembly, Rodinia breakup, and Pan-African orogeny. In principle, thermochronology, which allows us to determine time–temperature (and thus exhumation) histories, is well suited to resolve such questions. However, recent attempts (29–31), taken individually, fall short of truly resolving the critical questions.

First, the uncertainty of time–temperature (t–T) paths derived from a single thermochronometer can be large for older rocks—a problem sometimes exacerbated by the use of suboptimal inversion methodologies—making it difficult to discern between glacial and tectonic drivers by timing alone. Second, the magnitudes of both glacial and tectonic erosion are expected to be spatially heterogeneous. Fortunately, however, glacial and tectonic processes predict distinct spatial patterns of exhumation—with tectonic erosion focusing in tectonically active regions near cratonic margins and ice-sheet glacial erosion focusing in regions of wet-based ice—namely, in the models of Donnadieu et al. (33), broad regions of the low-latitude cratonic interiors away from ice divides, narrowing to a more “hit-or-miss” pattern at cratonic margins where basal slip is focused into only a few rapid outlet ice streams, as is observed at modern Greenland and Antarctic ice margins. Consequently, to resolve the relative contributions of all such climatic and tectonic drivers of erosion in the Neoproterozoic, not to mention their potential interactions, we require higher-resolution t–T paths from localities that can address the spatial pattern of Neoproterozoic exhumation at a global scale. Here we report robust Bayesian thermochronological inversions to test these hypotheses and our results show a widespread pattern of nearly synchronous Cryogenian rock cooling across North America that is interpreted as multiple kilometers of erosional exhumation due to ice-sheet glaciation.

Deep-Time Thermochronology

Thermochronology allows us to estimate the temperature that a mineral crystal has experienced over time and its position in the continental crust given a particular thermal structure. It provides a potential test for the contrasting hypotheses regarding the proposed link between widespread glaciation and cratonic exhumation, specifically linking snowball Earth glaciations to the phenomenon of widespread unconformity spanning the late Neoproterozoic. The use of multiple thermochronometers with varying temperature sensitivities is critical for such deep-time applications, because the parameter space of possible t–T paths grows with increasing timescale (34). Moreover, although a multichronometer approach can be time and resource intensive, the improved resolution critically allows for model results to be independently validated by testing against known geologic constraints, rather than merely forcing the model to fit such constraints a priori.

Recent reports in the thermochronologic literature indicate that nearly continuous thermal histories can be constrained using a multimethod approach (400 C and lower) that involves jointly inverting these data to effectively explore Precambrian histories, supplemented by existing high-temperature metamorphic data and stratigraphic constraints (34–37). In this context, the inclusion of medium-temperature (100 to 300 C) thermochronometers such as K-feldspar 40Ar/39Ar and zircon (U–Th)/He are especially important, since low-temperature systems (<100 C) tend to record only the most recent Phanerozoic overprints from burial reheating. A robust multichronometer approach featuring a full range of temperature sensitivities, however, should allow us to see past such overprints and accurately constrain the erosion history of ancient crystalline basement over ∼Ga timescales. To this end, we consider the following range of thermochronometers:

Potassium Feldspar 40Ar/39Ar Dating.

Potassium feldspar is notable for its ubiquity in crustal rocks, for containing appreciable amounts of radiogenic argon, and for containing domains of differing diffusion radius (38). The degassing behavior of domains can be characterized during laboratory 40Ar/39Ar step-heating experiments and mathematically modeled to determine the number of domains, relative size distribution, and kinetic parameters specific to each sample (39). This information can in turn be inverted to yield a continuous thermal history record between ∼350 and 150 C (34, 40) and provides a crucial link between high- and low-temperature thermochronometers.

Zircon (U–Th)/He and Apatite (U–Th)/He Dating.

Helium diffusivity in zircon and apatite is modulated by accrued alpha-radiation damage from radioactive decay in the crystal lattice (41–43). Higher-radiation damage in apatite correlates with higher He retentivity (i.e., lower diffusivity) (41). High-U zircon grains with greater radiation damage experience faster He diffusion rates over geologic time, whereas the opposite is true for low-U grains. Heating of these minerals causes annealing of accumulated radiation damage. Given certain t–T conditions and mineral chemistries, radiation damage effects manifest as large intrasample He date variation. Individual grains accumulate a predictable amount of radiation damage as a function of their U and Th concentration and t–T path, and multiple grains from the same sample with different U and Th concentrations will therefore each have a different respective He diffusivity and behave as an independent thermochronometer. The “effective uranium” of any grain can be represented by the single parameter eU (= [U] + 0.238 × [Th] + 0.0012 × [Sm]) (44), which weights the He contribution from each parent by its alpha-decay productivity. Date-eU trends provide much more powerful and informative thermal history information than any one date (34, 45). The use of many single-crystal dates provides useful information that can be inverted for thermal history, often spanning ∼200 to 40 C over a range of <100 ppm to >2,000 ppm eU for zircon and ∼100 to 50 C over <10 ppm to <200 ppm eU for apatite grains.

Apatite Fission-Track Dating.

The apatite fission-track (AFT) method is sensitive to temperatures between ∼110 and 60 C for most rocks that incorporate common apatite and for this reason is useful for determining upper crustal erosion and burial histories. Fission-track dating is based on quantifying (counting) the damage trails created from the energetic fission of 238U, which happens continuously at a known rate in the mineral crystal lattice (46). These “fission tracks” are then related to the amount of uranium present in a counted grain area to calculate an apparent “age” for an apatite grain or approximate time over which appreciable fission tracks have accumulated in the crystal (47). The production of fission tracks is continuous across a sample’s thermal history. Tracks initially have an etched length of ∼16 to 17 µm and shorten with heating, being totally annealed at >120 C (48, 49); thus, each track has a different age and records a different portion of the thermal history. Annealing resistance is also influenced by apatite chemical composition, notably Cl and other elemental substitutions (50, 51). Track lengths are measured since they can be used to model the style and magnitude of cooling (or partial reheating) experienced during a rock’s thermal history (52).

Evaluating Published Thermochronology Data from North America.

We examined previously published thermochronology data from the North American interior spread across the continent to adequately test models of the first-order spatial and temporal pattern of Neoproterozoic crustal exhumation (Fig. 1). Data were compiled from the East Lake Athabasca region (Saskatchewan, Canada) (53–55), Archean terranes in the Minnesota River Valley (Minnesota) (43, 56), the Ozark Plateau (Missouri) (36), and the Pikes Peak Batholith (Colorado) (29). The cratonic interior of North America provides an ideal locality for testing the various Great Unconformity formation hypotheses (29) when compared to paleo-margin locations because the craton has remained tectonically stable over the last ∼1.8 Ga, which alleviates most concerns about more recent, extensive thermal disturbances. In some situations, this allowed us to jointly model samples collected from a broader area of up to 100 km (i.e., Minnesota), under the assumption that over this scale these cratonic rocks have experienced similar thermal histories.

Fig. 1.

North American location map for previously published thermochronology datasets discussed in this paper. Sample locations (triangles): A, East Lake Athabasca region; M, Archean Minnesota River Valley terranes; O, Ozark Mountains; P, Pikes Peak batholith. Map shows structural and geologic features of the United States and Canada, adapted from Whitmeyer and Karlstrom (68) and Marshak et al. (71). Precambrian exposed outcrop is in pink and Phanerozoic orogens are in orange shading. Red lines are the edge of Cenozoic rifting in the west and the Appalachian front in the east from Marshak et al. (71). Major highlighted rifts are in gray that were active in the mid–late Neoproterozoic. MCR, Midcontinent Rift; OKA, Oklahoma aulacogen; RFR, Reelfoot Rift. Note that regional faults were active in the late Neoproterozoic at the Pikes Peak (69, 75) and Ozarks locations (140), whereas faulting near the Athabasca and Minnesota samples predated Cryogenian time [i.e., ca. 1.9 to 1.65 Ga (61) and ca. 1.9 Ga Penokean orogeny and/or 1.1 Ga MCR, respectively].

The QTQt software package (57) was used for Bayesian t–T inversion. Thermal-history modeling is often conducted using a simple Monte Carlo approach by searching for and selecting a subset of “acceptable” paths from a finite set of randomly generated t–T paths (58). However, the large parameter spaces of deep-time thermochronology are arguably better suited to an adaptive inversion methodology such as the reversible-jump Markov chain Monte Carlo (rjMCMC) approach used by QTQt (34). A key aspect of the rjMCMC method as implemented in QTQt is that the complexity of thermal-history solutions is inferred from the data rather than being defined a priori (57, 59). Beyond this, our approach differs from many routine thermochronometric studies by using Bayesian statistical methods for both the search algorithm and data uncertainty treatment, the generation of many more t–T paths (several orders of magnitude) during the course of modeling, and a distinctly empiricist philosophy regarding geologic “constraint box” implementation (). That is, we greatly minimize the use of such constraint boxes that force the model to take an expected path, allowing us in such cases to instead observe the ability of the model to independently infer geologically plausible paths from the thermochronologic data alone. We present the resulting t–T histories arranged in order from the cratonic interior outward toward the paleo-Laurentian margins (Fig. 1). The more interior locations generally include more thermochronometric systems and longer modeled time intervals and are characterized by less interpretive complexity (Fig. 2).

Fig. 2.

QTQt time–temperature inversions of thermochronology data from North America. Relative probability is proportional to t–T path density, where darker colors (or higher saturation) denote higher relative probability. Unless otherwise noted on the panels, hierarchical Bayes “error resampling” (57, 138) was implemented within QTQt (scaled from 1 to 100 times, with a value of 1 equal to the input uncertainty) and more complex models were accepted for equivalent likelihood (see Materials and Methods for details). Cyan bars are the time intervals for the respective Sturtian, Marinoan, and Gaskiers glaciations and white t–T boxes are geologic constraints. Plots showing observed and predicted data for each simulation are in . (A and B) Inversion results for East Lake Athabasca (Chipman domain) including the K-feldspar 40Ar/39Ar MDD age spectrum reported by McDannell et al. (55) and (U–Th)/He data reported by Flowers et al. (53) and Flowers (54). The modeled ZHe dates are from nearby sample 00-196C a few kilometers away. (A) Model without geologic constraints. (B) Model with a constraint box at 1,650 ± 50 Ma between 25 ± 25 C or the time of required basement exposure prior to Athabasca Basin formation (34) and a box at 545 ± 90 Ma and 20 ± 20 C to include uncertainty in surface exposure prior to Paleozoic burial onset in the adjoining Western Canada Basin. The Sturtian cooling trend is present in both models, with or without boxes. (C) Inversion results without constraint boxes for the Minnesota River Valley terranes data reported by Miltich (56) and Guenthner et al. (43). A separate model is shown in implementing geologic constraints of Sioux Quartzite deposition (1,695 ± 65 Ma) and a Precambrian-Cambrian near-surface constraint (600 ± 100 Ma) prior to Paleozoic burial. The latter box honors the paths at low temperatures in the unconstrained model. The Minnesota ZHe data underwent an empirical form of hierarchical Bayes resampling (Materials and Methods and ). (D) Ozarks model result with enforced geologic constraints as described by DeLucia et al. (36), except with an expanded “Cambrian” box; see text for details. An additional “no constraint” model is shown in . (E and F) Model inversion results for Pikes Peak batholith ZHe data from Flowers et al. (29). (E and F) More complex models were allowed and data underwent error resampling (E), whereas in F dates were randomly sampled within the assigned 10% SD and more complex models were accepted only if they improved model predictions with respect to the input ZHe data. Importantly, the Pikes Peak simulations do not incorporate constraint boxes (Fig. 3). An alternate model for Pikes Peak is shown in . See for linked QTQt files.

Fig. 3.

(A) The Flowers et al. (29) HeFTy (58) time–temperature model for Pikes Peak showing constraints used in their modeling (see text and for discussion of the nature of these constraints). (B) Pure Monte Carlo simulations where a simple script was used to generate random paths to pass through constraint boxes without including thermochronologic data. The simulation in B shows 500 randomly drawn paths (green) and a subset of 30 paths (purple) randomly drawn from those 500 to more clearly show overall path behavior. Path colors are only meant to resemble the default HeFTy scheme (58). All boxes are the same as in the Flowers et al. (29) model. Our model in B forces paths only through boxes and is very similar to the Flowers et al. (29) result (their figure 4 or A above). It is important to note that their best-fitting paths (in magenta, A) are nearly indistinguishable from a random sampling of 30 of our 500 Monte Carlo paths. A separate model in shows the result of utilizing only the Phanerozoic geologic constraints and the removal of the Precambrian interpretive boxes, also without thermochronology data. Results show that either early cooling to near-surface conditions (i.e., a Rodinia tectonic scenario) or late cooling during a Cryogenian glacial cooling scenario is allowed. Fig. 2 shows the results of modeling thermochronology data only (without boxes). The model is truncated at 200 C for plotting. The Monte Carlo script is included as a supplemental file; see .

East Lake Athabasca, Canadian Shield, Saskatchewan, Canada.

The East Lake Athabasca region lies along the Snowbird Tectonic Zone in the western Canadian Shield at the margin of the remnant ca. 1,700 to 1,650 Ma Athabasca Basin (60). High-temperature U–Pb (titanite, apatite, rutile; ∼650 to 400 C) and 40Ar/39Ar (hornblende, muscovite, biotite; ∼550 to 300 C) geochronology constrain episodic, post–1,900 Ma exhumation of granulites from the deep crust to the surface by 1,650 Ma (61). K-feldspar 40Ar/39Ar multidiffusion domain (MDD) data (∼350 to 150 C) seamlessly link published high- and low-temperature data and establish rapid cooling and exhumation to the near surface by 1,650 to 1,600 Ma (34, 55). Low-temperature thermochronological studies utilizing zircon (U–Th)/He (ZHe), AFT, and apatite (U–Th)/He (AHe) data imply thermal resetting and burial heating of the Athabasca region during the late Proterozoic and again during the early Paleozoic (34, 53, 54). This sample suite was also recently remodeled (without AFT data) in McDannell and Flowers (34), providing similar results using a different t–T search algorithm and different explicit model boundary conditions. This dataset is the most robust of all the locations studied due to the greater quantity of high-quality thermochronologic data.

Our QTQt thermal history simulations demonstrate rapid cooling to the surface by 1,600 Ma, reheating to ∼120 to <150 C, followed by cooling to the surface again from 750 to 600 Ma (Fig. 2 ). Minor reheating ensued during Cambro-Ordovician through Devonian time, in agreement with early deposition in the nearby Western Canada Basin. Geologic constraints are enforced in the model at the time of presumed cratonic basement exposure prior to Athabasca Basin formation and at the regional basement nonconformity prior to Paleozoic sedimentation (Fig. 2). A noteworthy outcome is that the integration of multiple Athabasca thermochronometers containing redundant or complementary kinetic information constrains a broad range of t–T space and yields similar model results for both the “unconstrained” (Fig. 2) and “constrained” (Fig. 2) models (i.e., no constraint boxes compared to the model with constraint boxes). These t–T models suggest ∼3 to 4 km of exhumation during the Sturtian and Marinoan Snowball glaciations in this intracratonic setting (assuming a 25 to 35 C/km paleo-geothermal gradient and 20 C surface temperature, used throughout this paper for any exhumation calculations).

Minnesota River Valley Terranes, Southwestern Minnesota.

Minnesota hosts some of the oldest exposed rocks in the United States (Fig. 1), including the 3.5 Ga Morton and Montevideo Gneiss units—both of which are intruded by the 2.6 Ga Sacred Heart Granite (62). Paleoarchean Minnesota River Valley Terranes (MRVT) make up the southernmost extension of the Canadian Superior Province and lie west of the 1,100 Ma Midcontinent Rift (MCR) (Fig. 1). The late Proterozoic surface history of the southwestern Minnesota Archean basement is poorly known; however, the preserved Sioux Quartzite to the south of these samples is a unit that was deposited at ca. 1,760 to 1,630 Ma (63). Regional geologic relationships demonstrate that the Archean crystalline basement was exposed (64) prior to burial during Cambrian through Devonian time, followed by burial again in the Jurassic-Cretaceous (65).

We modeled the ZHe and AHe data reported by Miltich (56) and Guenthner et al. (43). The QTQt model results (Fig. 2) suggest cooling ensued ca. 750 to 650 Ma after maximum heating by ca. 800 Ma that obscures the pre-1,000 Ma history. It is conceivable that the reheating that concluded by 800 Ma was due to burial by erosional detritus shed from the nearby Grenville orogenic belt (66). The Phanerozoic model history is characterized by Cambrian through Devonian reheating, followed by cooling and a second reheating event that peaks in the Cretaceous, both of which agree with the preserved regional geology (). An amphibolite inclusion from the Sacred Heart Granite dated by C. Naeser in 1974 yielded an AFT age of 460 ± 45 Ma (1σ) (67), which is in broad agreement with our Phanerozoic model results showing cooling through the fission-track partial annealing zone (120 to 60 C) after 500 Ma. These samples are deep in the continental interior and there is no evidence for faulting associated with Rodinia breakup in Minnesota. The effects of 1,100 Ma MCR faulting and rifting were localized and would have been followed by thermal subsidence. Regardless, all of the MCR events preceded the Cryogenian, and thus we anticipate all cooling from 200 C at ca. 720 to 650 Ma to be associated with >4 km exhumation resulting from Cryogenian glacial erosion.

Ozark Plateau, St. Francois Mountains, Missouri.

The Ozark thermochronology dataset was published by DeLucia et al. (36), who carried out a combination of forward models to test endmember geologic scenarios under different conditions, as well as inverse t–T models to explain their ZHe data. Individual zircon dates were binned by eU and averaged to create “synthetic” data for use in the HeFTy software (58); also see Pikes Peak Batholith, Colorado below. They interpreted their HeFTy model results as burial due to Rodinia assembly and Grenville orogenesis from 1,200 to 1,000 Ma followed by significant Neoproterozoic cooling of ∼220 to 200 C that they related to the breakup of supercontinent Rodinia. They concluded that exhumation led to increased weathering and CO2 drawdown, triggering snowball glaciation. Sedimentary burial over the course of the Paleozoic-Mesozoic abruptly ceased with rapid cooling from 225 to 150 Ma, interpreted as uplift and exhumation from the breakup of supercontinent Pangaea. The Ozark Plateau, like western Colorado (see below), was near the paleo-cratonic margin in the Neoproterozoic-Paleozoic undergoing normal faulting and regional extension (68–70). This area hosts extensive structural lineament systems and faults, including the Ste. Genevieve, Cottage Grove, and Rough Creek fault zones, the larger Reelfoot Rift (Fig. 1), and the active New Madrid Seismic Zone (71). The ∼7.5 km of structural relief that exists between the Great Unconformity exposed in the St. Francois Mountains and the buried Great Unconformity surface in the adjacent Illinois Basin attests to late Precambrian tectonic deformation (71).

Our QTQt inversions are shown with geologic constraints (Fig. 2) and yield results broadly consistent with those in DeLucia et al. (36), albeit with a greater number of t–T paths generated during the course of modeling and the use of single-grain ZHe data. A model without explicit t–T constraints () clearly shows that the thermal event that set the AFT data also obscures the sensitivity of the ZHe data prior to that time. Therefore, we used the same modeling constraints as DeLucia et al. (36). Following ca. 1,450 Ma granite emplacement from high temperatures, the geologic constraint at 1,365 ± 15 Ma and 50 ± 50 C represents cooling of surficial rhyolite or hypabyssal granites to near-surface temperatures (Fig. 2). The late Cambrian Lamotte sandstone rests unconformably on the Great Unconformity surface and is represented by the constraint box at 560 ± 75 Ma and 20 ± 20 C. Pevehouse et al. (72) suggest weathering and soil formation occurred in the Ozarks after the last major Neoproterozoic glaciation. We have expanded the surface t–T box to include Ediacaran time to account for possible subaerial exposure prior to sandstone deposition (72) and to accommodate elevated Cambrian ocean temperatures (73). The Ozarks ZHe inversion shows reheating between ca. 1,300 and 800 Ma and cooling to surface temperatures by the Cambrian (Fig. 2). The timing of cooling from peak temperatures of ∼250 to 200 C is poorly constrained between ca. 800 and 650 Ma, albeit still consistent with both “Rodinia breakup” exhumation and snowball Earth glaciations.

Pikes Peak Batholith, Colorado.

Flowers et al. (29) published a ZHe dataset from the Pikes Peak batholith in Colorado. They modeled synthetic ZHe data (see for further discussion) collected from samples below the Great Unconformity surface and other fault block locations in their study area. Flowers et al. (29) interpreted their t–T results from this single location as unroofing due to global tectonic activity related to supercontinent Rodinia assembly and/or breakup. While such a model would not be incompatible with a glacial model for the origin of the Great Unconformity, given the tectonic activity of the Pikes Peak region in the Neoproterozoic (as shown by the fault-bounded nature of many Tavakaiv bodies; see below), several aspects of their interpretation warrant a critical reexamination. Their (29) t–T modeling hinges on assuming shallow emplacement of the enigmatic Tavakaiv quartzite injectites* (75–77) near the paleosurface at 676 ± 26 Ma from hematite (U–Th)/He data published by Jensen et al. (78) (Fig. 3). However, the depth of Tavakaiv emplacement is uncertain due to an unknown emplacement mechanism and the hematite He data can be interpreted as either mineralization or cooling ages (78). The cooling-age interpretation (our preferred model) requires Neoproterozoic burial reheating (78), which is anticipated near Pikes Peak given the striking similarities between detrital zircon U–Pb age distributions for the Tavakaiv dikes and regional Neoproterozoic reference ages (75, 76) (see for details). Given the enigmatic nature of Tavakaiv emplacement, their model design could be more accurately described as a compatibility test between the thermochronologic and detrital zircon data; however, the authors presented shallow Tavakaiv emplacement as an a priori constraint and forced their t–T paths to conform to this constraint.

Regardless of the interpretive framework to explain the thermochronology data, the t–T models published by Flowers et al. (29) were largely controlled by their use of constraint boxes (79) in the HeFTy software (58) (Fig. 3). We verified this by generating random Monte Carlo t–T paths using a simple script that incorporated only their constraint boxes without including thermochronologic data (Fig. 3). In our model (Fig. 3), random paths were simply forced through the boxes, yielding the same results as in Flowers et al. (29). The box control on modeling is evident from specific placement of their Great Unconformity “exploration field” (Fig. 3, blue box). This interpretive box (and the other Precambrian boxes) prevents exploration and forced cooling to occur prior to (or by) 720 Ma because paths are required to be between 20 and 0 C from 1,000 to 720 Ma in the model. There is no physical geologic evidence to support a pre-720 Ma surface condition and it is not demanded by the ZHe data (see for details). We ran additional “no data” Monte Carlo simulations without Precambrian surface constraints and there are Neoproterozoic cooling paths that satisfy either the glacial or tectonic exhumation hypotheses when not forced to cool to surface temperatures prior to 720 Ma (). The full range of possible t–T paths is also shown after removal of the nested Paleozoic and Mesozoic boxes derived from their assumptions regarding the Pikes Peak history in the Phanerozoic. The example in Fig. 2 shows the results of using only the thermochronology data to resolve the thermal history without relying on interpretive t–T boxes.

To better understand the thermal history of the Pikes Peak region, we remodeled the Flowers et al. (29) Pikes Peak ZHe dataset using QTQt. Importantly, we applied no constraint boxes; any variations from uniform path density in the results reflect instead the information contained in the 12 measured single-grain ZHe dates from their GU surface samples F1936 and F1937 (29). The resulting t–T history (Fig. 2 ) exhibits Neoproterozoic cooling from ∼220 to 200 C at ∼745 to 700 Ma to near-surface temperatures by ∼660 to 600 Ma. The model in Fig. 2 is an alternate version where t–T points were accepted only if they resulted in better prediction of the observed dates (i.e., model paths are only as complex as necessary to optimize the data fit between the model and the observations). The latter model is shown only to provide a lower limit on the complexity required to reproduce the ZHe data and reduces noise in Fig. 2. It is obvious that the greatest resolution lies near 200 C at ca. 700 to 660 Ma (constrained by the oldest ZHe grains), followed by cooling to surface before 600 Ma, and a late reheating event to <150 C at <100 Ma, presumably due to burial from the Laramide orogeny. Any heating that may have occurred between 600 and 100 Ma must be <150 C and is not necessarily required or well resolved by the Pikes Peak ZHe data (). Mid-Paleozoic burial is also not required and basement rocks are not presently mantled by sedimentary cover in the field (29). The Tavakaiv quartzite injectite emplacement age of 676 ± 26 Ma from Jensen et al. (78) and the geologic constraint of basement being exhumed to the surface prior to Sawatch sandstone deposition in the Cambrian are honored in our simulation without imposing t–T constraint boxes (Fig. 2 ). It is possible that faulting, Tavakaiv emplacement, and basement exhumation were coincident near 700 to 650 Ma due to Rodinia breakup along the cratonic “margin” and Snowball ice-sheet dynamics (see for further discussion). The results of our t–T inversion for Pikes Peak basement ZHe data offer support for this scenario while still honoring the interpretation of coeval hematite resetting/cooling and injectite emplacement from 200 C to near-surface conditions during the Cryogenian (78). The Neoproterozoic cooling segment in our model is consistent with both the Sturtian and Marinoan glaciations and Rodinia breakup resulting in up to ∼5 to 7 km of erosional exhumation.

Reconciling Neoproterozoic Exhumation Trends

Spatial Patterns of Tectonic and Glacial Erosion of Continents.

McDannell et al. (55) and DeLucia et al. (36) came to the conclusion that kilometer-scale Neoproterozoic exhumation occurred after 1 Ga within the North American interior and linked this to formation of the Great Unconformity due to Rodinian geodynamics and/or snowball Earth glaciations. These two hypotheses are not mutually exclusive—it is possible that both tectonics and glaciation contributed to global Earth system disruption (80, 81) during formation of the Great Unconformity. Glaciation would be most effective as a driver of erosion in regions with preexisting topography (be it from rifting or orogeny); therefore, erosional synergy between tectonics and ice sheets is a possibility (82). Ultimately with respect to the Great Unconformity, it may be that the generally accepted reconstruction(s) of more concentrated equatorial packing of the Rodinian continents (11, 83), along with the unique environmental conditions of the Neoproterozoic, proved to be a time of geologic serendipity unlike most any other in Earth history.

Direct and meaningful comparisons between tectonic and glacial unconformity hypotheses are complicated by the fact that there are precise estimates for the timing of Snowball glaciations (23), whereas the timing and duration of Rodinia assembly and breakup remain incompletely understood due to discrepancies between paleomagnetic and geologic data (11, 83, 84). Rodinia assembly and breakup occurred episodically and diachronously over at least 250 Ma for each phase, with timing dependent upon location (10, 11). Invocation of Rodinian tectonics as a primary, global cause of the Great Unconformity partly requires a consensus or at least reconciliation of the myriad configurations of the supercontinent (11, 70, 83, 85–88) to construct valid geodynamic models of uplift during the supercontinent cycle. Otherwise, any thermochronologic cooling signal can simply be attributed to “Rodinian tectonics” in the Neoproterozoic. Notwithstanding Rodinia’s exact arrangement, the majority of rift-related deformation and exhumation would have been confined to cratonic margins or to localized horst–graben systems (89). A question that arises by appealing to “tectonics” as a global cause of the Great Unconformity is, Why do we not observe an equivalent hiatus as a result of the assembly and breakup of other supercontinents such as Pangaea? If supercontinent cyclicity caused global unconformities akin to the Great Unconformity, we anticipate that the North American Sauk Sequence (as currently defined) would instead occur in the late Mesozoic due to capture by Pangaean erosion. The lull in Pangaean sediment volume (8) during supercontinent breakup is apparently instead due to nondeposition during a sea-level low stand—and is not accompanied by the same stepwise difference in sediment volume that occurs prior to the beginning of the Phanerozoic (5).

The dynamics of supercontinent breakup remain poorly understood (90), but remain a focus of discussion here since the timing of rifting in North America closely overlaps with Snowball glaciations and the timing of cooling in our t–T inversions. Mantle-plume push (i.e., “bottom–up” processes) (91) and plate boundary dynamics (i.e., subduction retreat or “top–down” processes) (92) both govern supercontinent breakup (90, 93). Mantle plumes initiate breakup (94), as evidenced by large igneous province eruptions that are either the cause or the manifestation of supercontinent demise (90). Successful rifting results in a passive margin and the high number of passive margins during staged Rodinia disassembly (95) implicates Laurentian margin rifting as the dominant mode and locus of tectonic activity during the Neoproterozoic. Longstanding models suggest supercontinents insulate the mantle, causing upwelling and breakup (96); however, recent work suggests that subduction plays a dominant role in subcontinental mantle upwellings (97). Laurentia may not have had well-established margin subduction zones until ca. 600 to 540 Ma (98), which broadly explain the formation (i.e., subduction-related dynamic topography) of North American cratonic unconformities (99) in the Phanerozoic (100)—leaving early Neoproterozoic continental dynamics an open question.

A dynamic topographic response to mantle convection anomalies can produce low-amplitude surface uplift (101), tilting, and erosion across a continental interior over a few million years (102), although this often involves a complex interplay between plate motions and mantle swell position, topography and drainage network organization, and climate change (103)—which are exceedingly difficult to quantify in the Proterozoic. The erosional response to dynamic uplift is proportional to the upwelling wavelength (104); therefore, dynamic topography would be required at the scale of the North American continent to induce widespread erosion that agrees with our models. Continental erosion would likely be limited within the interior (<1 to 2 km) and occur relatively slowly over many tens of millions of years (102, 105) in the absence of significant modification of the cratonic lithosphere (55). This is considerably less than the amount of unroofing suggested by our t–T models. However speculative, an episode of widespread kilometer-scale epeirogenic uplift associated with a thermally buoyant Rodinia supercontinent (106, 107) may have led to increased continental exposure and the formation of the Great Unconformity on multiple continents. Erosional detritus would have in turn influenced ocean chemistry and atmospheric CO2 concentrations that contributed to snowball Earth glaciations (22, 108–110).

Conversely, Snowball glaciations could have been the main driver of erosion that created the Great Unconformity. Through a combination of wet-based glacial sliding and lowering of erosional base level, global glaciations in the late Neoproterozoic could have removed several kilometers of rock (including cratonic sedimentary rocks) to produce the Great Unconformity surface. Notably, this would not require incision rates any different from those observed in modern ice-sheet environments. A scenario where modest continental ice-sheet incision rates are effectively constant at 0.05 to 0.1 km/My yields 2.9 to 5.8 km of exhumation over the Sturtian glacial interval alone. Large amounts of exhumation could be accomplished either at lower rates for prolonged periods of basal ice sliding or more rapidly over short intervals during deglaciation. For example, Cowton et al. (111) indicated that the modern Greenland ice sheet erosion rate is ∼2.2 to 7.4 km/My (from the ice margin to >50 km inland) during the deglacial phase, which is at least an order of magnitude higher than previously established ice-sheet erosion rate estimates (112) and places incision rates on par with empirical estimates of ∼1 to >10 km/My for temperate glaciers (113). The results of Neoproterozoic ice-sheet simulations demonstrate that only high-latitude Rodinian cratons (i.e., not Laurentia) would have been characterized by cold-based ice, with low-latitude interior basal ice temperatures near 0 C and continental basal sliding displacement rates of ∼1 to >10 m/y (33). Furthermore, glacial incision is expected to increase with decreasing latitude (113) and the low-latitude position of Rodinia during the late Neoproterozoic favored increased continental weatherability and precipitation rates (114), thus creating a relationship where erosion would be maximized with lubricated basal ice increasing sliding—leading to more rapid erosion (115).

Cratonic interiors provide the only location to truly test and differentiate the hypotheses of pre-, syn-, or post-Cryogenian formation of the Great Unconformity. Timing is a key component of this signal, but spatial pattern and magnitude of exhumational rock cooling are also critical. Tectonic rifting and glacial erosion will produce opposing spatial patterns of exhumation and different magnitudes of crustal unroofing across a continent. The majority of exhumation associated with supercontinent assembly and breakup would be limited to compressional orogenic belts and extensional (faulted) rift margins, respectively. Rifting will show large exhumation narrowly restricted to continental margins, where tectonic activity is highest, whereas stable continental interiors will experience little to no erosion or even deposition. In addition to orogenic erosion, intraplate stresses manifest as continental extension (69), causing subsidence and burial across a craton (116–118). This is hypothesized for the Rodinian interior during terminal assembly and incipient breakup (66) and agrees with the consistent heating signal seen in our thermochronological inversions (Fig. 2). We would expect most tectonic uplift and erosion to occur during early supercontinent assembly and orogenesis, rather than breakup. Thus, the rock-cooling signals for Rodinia assembly (ca. 1,300 to 900 Ma) (11), major rift breakup phases (ca. 850 to 680 Ma) (12, 98), and snowball Earth glaciations (ca. 720 to 635 Ma) should be rather distinct in terms of timing and location. As an example, recent work by Ricketts et al. (31) apparently shows exhumation that broadly aligns with exhumation during Rodinia assembly in the southwestern United States. While they did not jointly invert 40Ar/39Ar and zircon (U–Th)/He data, early or episodic Neoproterozoic exhumation may nevertheless be expected locally, since western North America was undergoing active tectonism during that time (69).

In contrast, long-term glacial erosion will produce high-magnitude exhumation over areas of thousands of square kilometers, with ice-sheet margins experiencing either very little or extremely high incision due to fluctuating ice dynamics and runoff (33, 108). The timing of cooling in our models is coincident with both rifting and glaciation in western North America. We would expect the tectonic versus glacial signals of exhumation on a paleo-cratonic margin (or at least areas experiencing Rodinian syn-rift breakup deformation) to be nearly indistinguishable from one another—as observed in the Pikes Peak region. Locations such as Athabasca are too far from continental margins to have experienced >3 km of erosion over a short interval solely due to rifting. Moreover, if there is widespread erosion during a “hard snowball” glaciation, ice would have to be a dominant erosive agent. The only foreseeable way to obtain a consistent, high-magnitude, and synchronous Cryogenian unroofing signal at the continental scale is through ice-sheet glaciation. Our t–T model results demonstrate the viability of such an exhumation pattern across North America.

Deep Continental Ice-Sheet Erosion.

Widespread, deep Neoproterozoic glacial erosion (119) may appear to contradict the oft-held perception that continental ice sheets cannot deeply erode the upper crust (120, 121). Early estimates of physical erosion as a result of the Laurentide glacial episode were that Approximately 120 m of rock was removed over the last 3 My across upper North America (122), a rate that would equate to some ∼2.5 km over the duration of the Cryogenian glaciations. However, Laurentide glaciation is perhaps a poor analog to Cryogenian glaciations, since continental freeboard was fundamentally different (i.e., lower) during the Cryogenian, providing a gravitational potential energy gradient essential for deep glacial incision (5). Net base-level fall during snowball Earth termination (and shortly after glaciation) is predicted to be the greatest (up to –600 m) in continental interiors, decreasing toward margins (22), whereas estimates suggest less than –120 m relative sea-level fall during the Laurentide (123). The Laurentide ice divide was positioned over the Hudson Bay Basin where preservation of sedimentary strata was likely due to low rates of basal sliding (124). However, the simple observation that, beyond this ice divide, the thickest parts of the Laurentide ice sheet (125) match the currently exposed extent of the Canadian Shield implies that any continental ice sheet is capable of denuding the craton.

An underappreciated aspect of the “deep erosion” argument is that continental-scale exhumation need not imply that most of the crust removed was crystalline basement; on the contrary, a substantial portion of the eroded crust may well have been intracratonic sedimentary rocks deposited during the Proterozoic across the continental interior (66, 126). Geology and our inversions directly indicate burial heating of basement was probably due to thick Proterozoic cover for (at least) the Athabasca region and the Ozark Plateau. In support of this, global average zircon Hf and O isotope anomalies were interpreted as old crustal material from the Earth’s surface being subducted and incorporated into new magmas in the Neoproterozoic (5). The Hf and 18O isotopic signatures require only surface exposure and subduction of crust containing ancient zircons—whether that material was directly sourced from Precambrian basement or recycled from Proterozoic basins makes little difference. Ocean basins serve as the main repository for sediments produced during ice-sheet denudation (119, 122) and, due to the shorter oceanic crust lifecycle (compared to continental crust), provide one explanation for the reduced survival rate of Proterozoic detritus that is evident in the Ronov (9) compilations. This conceptually agrees well with the observation that many Archean and Proterozoic terranes have experienced relatively modest amounts of net crustal erosion (127), partially explains the variability and regional lack of evidence for snowball Earth glacial incision (128), and agrees with time-averaged measurements of net continental exhumation rates that approach zero over gigayear timescales (129).

Thermochronologic Support for a Glacial Unconformity

The anomalous abundance of unconformities near the Proterozoic-Phanerozoic boundary—each one different, and frequently composite, but evidently captured by a globally widespread erosive event—is what makes the Great Unconformity unique. Neoproterozoic glacial erosion, which we interpret as the primary cause of the Great Unconformity, is detected in North American thermochronometry without making numerous assumptions about past conditions. We stress that assumptions about past geologic conditions should not be prescribed as evident or imposed in lieu of quantitative thermochronology in thermal-history models. Our thermochronological inversions honor the measured isotopic data and physical geology, while demonstrating that the late Proterozoic basement nonconformity is a feature that 1) manifests as large-magnitude erosion between ca. 700 and 635 Ma, 2) maintains consistency across North America for multiple locations over 1,000 km, and 3) can be interpreted as widespread (albeit likely spatially heterogeneous) erosional unroofing of at least 3 to 5 km. Collectively these features can only be readily satisfied by a Cryogenian glacial model for exhumation of rocks sampled from both proximal and distal reaches of exposed Laurentian cratonic basement. It is important to note that this major denudation event does not preclude later, minor sub–kilometer-scale erosion (or nondeposition) that undoubtedly occurred across the craton prior to Cambrian flooding of the continent. The removal of ∼3 to 5 km of thick Mesoproterozoic basin rocks and upper crust from the craton likely caused a disturbance to the stable crustal thermal structure—leaving it warm and isostatically buoyant and thereby inhibiting extensive deposition until Paleozoic transgressions during Pannotia-Gondwana plate reorganization (130).

Development of the Great Unconformity as a physical surface is constrained in this work only between the Cryogenian erosion pulse observed in our t–T models and the age of the overlying sediments—therefore, we do not rule out a multistage or multiprocess model for the individual unconformity surfaces associated with the Great Unconformity as a broader phenomenon. However, to create and subsequently preserve a widespread unconformity by aggradation, most topographic relief must be removed and the landscape needs to be at (or below) base level (6)—which is difficult to achieve by fluvial or hillslope processes alone. It may be that continental-scale glaciation is the only foreseeable process that can account for both the formation and preservation of the Great Unconformity. Major unconformities, or significant step changes in North American (or global) sediment abundance, are not observed during other times of equatorial continental assembly, potentially invalidating supercontinent tectonic activity as the primary or sole driver of Neoproterozoic exhumation. In our view, it is not a coincidence that the thermochronologic inversions shown here demonstrate nearly synchronous exhumation transpiring across a vast region of North America during a known period of apparent worldwide glaciation. We present a more comprehensive appraisal for the origin of the Great Unconformity within North America that serves as a template for assessing exhumation globally to necessarily test further the hypothesis of a glacial origin due to snowball Earth conditions in the Neoproterozoic.

Materials and Methods

Inverse t–T simulations are presented for samples from the North American interior and were modeled using the QTQt v. 5.8.0 software (57). The QTQt program utilizes Bayesian statistics and a rjMCMC search method. We modeled K-feldspar 40Ar/39Ar, ZHe, AFT, and AHe data, implementing the MDD model of Lovera et al. (39), the zircon radiation damage accumulation and annealing model (ZRDAAM) of Guenthner et al. (43), the AFT multikinetic annealing model of Ketcham et al. (131), and the AHe radiation damage (RDAAM) kinetic model of Flowers et al. (42) for each respective thermochronometer in our surveyed datasets. To encourage thorough exploration of t–T space, more complex models were accepted for equivalent likelihood and proposal jumps were rejected if they were proposed outside of the general prior (t–T model space) in QTQt. A total of 1,000,000 models were completed for each example, with 500,000 burn-in iterations that were discarded and an additional 500,000 iterations retained post burn-in for each simulation. The acceptance rates were within the recommended range of ∼0.2 to 0.5 and the sampling distribution reached stationarity under these conditions, which collectively signify model convergence (57).

Quantification of Data Uncertainties.

Currently, uncertainties related to eU estimation (132), U–Th isotopic zonation (133, 134), and imperfect grain geometries are not easily or routinely characterized; therefore, it is reasonable to assume that single-grain date uncertainties at the 2σ level are underestimated for both zircon and apatite (U–Th)/He thermochronometry. It is customary for analytical errors to be calculated from the propagated uncertainty from U, Th, and He measurements. Uncertainties are on the order of ∼1 to 5% and typically about 2 to 3% (135). However, the uncertainties including the Ft correction for alpha ejection are commonly greater, and the reproducibility of laboratory age standards yields total uncertainties nearer to 8 to 10% for zircon and ∼6 to 7% for apatite (135). These error estimates are more realistic, yet still conservative, and correspond to two SDs typically observed on replicate single-grain laboratory age standard Fish Canyon Tuff zircon and Durango apatite analyses (132, 135, 136). The age reproducibility estimated for large numbers of replicate analyses of natural AHe samples is much worse, on the order of 15 to 20% or more (137). We usually applied 6% uncertainty for AHe dates (typical Durango apatite reproducibility) and 8 to 10% for ZHe dates (135) if reported uncertainties were less than these values before modeling. During modeling, dates were randomly sampled from a normal distribution centered on the reported/assigned error (scaled from 1 to 100 times the input error), which we refer to as “error resampling,” a form of hierarchical Bayes resampling utilized in QTQt where the data are used directly for uncertainty inference and the variance of the data errors is estimated from their most probable value, given the data (57, 138). In scenarios where there are abundant, dispersed data of varying quality (i.e., Minnesota dataset), another type of empirical Bayes resampling was utilized to explore ZHe date uncertainties. The aim was to expand uncertainty accounting where the prior hyperparameters (i.e., observed dates) will have a prior distribution that expresses their initial uncertainty and a posterior distribution that is determined by the data directly (138). The individual date errors were treated as hyperparameters drawn from a probability distribution and the data variance was used to infer date uncertainty. Importantly, observed dates were modeled but the weighted uncertainty was inferred from the scatter of the data as determined by the SD of the data weighted by a Gaussian kernel in eU space (σ eU = 100 ppm). The empirical Bayes resampling code is available as a Jupyter notebook from https://github.com/kmcdannell/helium-empirical-bayes.git.

Athabasca.

We modeled the K-feldspar MDD sample 02-123A from McDannell et al. (55). Refer to McDannell and Flowers (34) for further information on sample data. QTQt modeling information is as follows: general prior (t–T model space) 900 ± 900 Ma and 200 ± 200 C with an imposed 10 C/My maximum heating/cooling rate. Model was truncated at 300 C for plotting purposes.

Minnesota.

We modeled the ZHe and AHe samples contained primarily in the Miltich (56) thesis and Guenthner et al. (43). The Minnesota ZHe samples underwent empirical Bayes resampling due to the greater number of scattered ZHe (n = 22) dates and the extreme timescale involved in modeling (approximately two to three times that of other examples). The majority of reported MRVT (U–Th)/He dates ranged from ca. 925 to 10 Ma (zircon) and ca. 1,725 to 125 Ma (apatite) (56). Extreme age overdispersion of over 1 Ga affected the apatite grains, which were noted as poor quality by Miltich (56). We refrained from modeling the oldest uncorrected dates because they were typically characterized by very small grain sizes (∼30- to 40-µm halfwidths) and were much older than the more numerous ca. 300 to 200 Ma grains. Most raw (no Ft correction) AHe dates ranged from about 270 ± 90 Ma over a range of 37 ± 34 ppm eU. We conservatively applied 10% errors to the MRVT apatites (n = 11 of 16 total analyses) due to the questionable quality of the data—but did not utilize hierarchical error resampling since the dataset likely contains both representative and extreme outlier ages. In this case, error resampling would incorrectly treat all observed dates as reliable, yet more uncertain than initially quantified. The oldest dates were excluded as clear outliers because they were much older than the mean age and during simulation trials they were among the highest misfit grains in the inversions (i.e., grains older than ∼400 Ma were instead always predicted between ca. 200 and 350 Ma). The remaining dates form a positive date–eU trend that “plateaus” at high eU and generally aligns with the RDAAM expectations. QTQt modeling information is as follows: general prior (t–T model space) 1,500 ± 1,500 Ma and 150 ± 150 C with an imposed 5 C/My maximum heating/cooling rate. Constraint boxes represent Sioux Quartzite deposition at 1,695 ± 65 Ma and 40 ± 40 C and late Precambrian basement exposure 25 ± 25 C prior to late Cambrian Mt. Simon sandstone deposition (600 ± 100 Ma; the unconstrained MRVT model shows solutions at near-surface temperatures during this entire interval, supporting box placement).

Ozarks.

We remodeled ZHe (samples 14OZ01 and 14OZ11; n = 10), AFT (sample 14OZ07), and AHe data (sample 14OZ11; n = 6) collected from basement below the Great Unconformity surface in the St. Francois Mountains of Missouri from DeLucia et al. (36). The ZHe samples that provided the broadest range in dates and eU were chosen for modeling (∼1,050 to 180 Ma and ∼400 to 1,800 ppm eU). The dates from the other samples cluster around ∼700 to 600 Ma. The AFT sample central age is 185 ± 16 Ma (n = 20) and mean track length is 13.54 ± 1.23 µm (n = 78) with a mean Dpar (track etch pit diameter) of 1.75 µm. The AHe sample contains six grains (<15 ppm eU) with dates between ∼210 and 150 Ma. This information alone signifies heating to temperatures >100 to 120 C near 200 Ma to cause thermal resetting of the AFT system followed by relatively rapid cooling through ∼110 to 60 C. QTQt modeling information is as follows: general prior (t–T model space) 725 ± 725 Ma and 150 ± 150 C with an imposed 5 C/My maximum heating/cooling rate. Error resampling (1 to 100 times) for ZHe data and complex models allowed for both scenarios.

Pikes Peak.

We remodeled zircon (U–Th)/He data from Pikes Peak samples F1936 and F1937 collected from Great Unconformity surfaces reported by Flowers et al. (29). The 12 single-grain dates span between ∼1,000 to 45 Ma and ∼30 to 2,000 ppm eU. QTQt modeling information is as follows: general prior (t–T model space) 538 ± 538 Ma and 150 ± 150 C with an imposed maximum heating/cooling rate of 5 C/My. Error resampling (1 to 100 times) for ZHe data and more complex models allowed for Fig. 2. The Fig. 2 model did not undergo error resampling and more complex models were rejected for equivalent likelihood values. Therefore, proposed t–T paths were accepted only if they provided a better fit to the data.