New thermochronometric constraints on the Tertiary landscape evolution of the central and eastern Grand Canyon, Arizona

Near the end of Cretaceous time, the area of the eastern Grand Canyon was buried under 2–3 km of sediment shed from the Sevier and Mogollon highlands to the west and south (Elston and Young, 1991; Burchfiel et al., 1992; Dumitru et al., 1994; Flowers et al., 2008). Today, most of the Mesozoic and Cenozoic sediment has been removed and the Colorado River has incised the 1.6-km-deep Grand Canyon into the Paleozoic and Precambrian rocks of the southwestern Colorado Plateau. The earliest studies by J.W. Powell and C.E. Dutton noted the curious course of the Colorado River as it bisects topographic highs and largely ignores structural relief, yet appears to have integrated along its present course since 6 Ma (Hunt, 1956, 1969; McKee et al., 1967). Because of erosion, the lack of a Cenozoic rock record has hampered the descriptions of fluvial systems and landscape morphology that are critical to identifying the scenarios and processes that carved the Grand Canyon; for this reason the canyon has remained an enigmatic and controversial landscape feature since the first scientific description more than 100 years ago.

Given that numerical constraints on landscape evolution are essential for resolving the controversial origin of the Grand Canyon, the goal of this study is to present new low-temperature thermochronometric constraints on landscape development utilizing both apatite fission track (AFT) and (U-Th)/He (AHe) data. This new synthesis helps provide substantial resolution to recent debates about the age of the Grand Canyon (Pederson, 2008; Polyak et al., 2008; Karlstrom et al., 2008; Flowers et al., 2008; Wernicke, 2011; Luchitta et al., 2011).

The modern Colorado River flows westward across the southwestern margin of the Colorado Plateau. However, paleocurrent direction indicators and provenance studies in the early to middle Cenozoic gravels that fill the older paleocanyons indicate a regional northeast-flowing drainage system (Potochnik, 2001; Young, 2001), as shown in Figure 1. Rim gravels from the surrounding plateau margin contain Precambrian basement clasts shed from the eroding source terrain of the Cretaceous Mogollon highlands to the west and south (Bilodeau, 1986; Elston and Young, 1991; Burchfiel et al., 1992). These rim gravels are the few remaining sedimentary relicts of a denudation event that is well documented from earlier thermochronometer studies (Dumitru et al., 1994; Kelley et al., 2001; Flowers et al., 2008) and supported by our results here. Flowers et al. (2008) presented extensive thermochronometric data from the region that indicate late Mesozoic to middle Cenozoic denudation and outline the broad southwest to northeast erosional patterns of the southwestern Colorado Plateau. However, the lack of high-density data sets that target areas of high topographic relief results in limited insight when deciphering the distinct exhumation histories across the length of the Grand Canyon area.

Most workers agree that the Colorado River system as it is today did not exit the Colorado Plateau until slightly before 5.5 Ma (Lucchitta, 2003; Lucchitta and Richard, 2001; Pederson, 2008, Karlstrom et al., 2012). The western Grand Canyon took its shape after emplacement of the Late Miocene Shivwits basalt (Lucchitta and Richard, 2001) and after deposition of the 12–6 Ma Muddy Creek Formation and Hualapai Limestone (Castor and Faulds, 2001; Faulds et al., 2001; Spencer et al., 2001). Following introduction of the Colorado River to the Grand Wash trough, eventual integration quickly followed multiple spillover events in the Basin and Range to a lowered base level in the Gulf of California that had been opening since the latest Miocene (McDougall et al., 1999; Oskin and Stock, 2003). Sediments from the Colorado Plateau first reached the Gulf of California ca. 5.4 Ma (Dorsey et al., 2005), indicating that the Colorado River system as we define it today had finally achieved its modern course by that time.

However, a few studies have suggested the Colorado River (or a precursor river) may have exited the Colorado Plateau at the present location much earlier (Hill and Ranney, 2008; Polyak et al., 2008); Wernicke (2011) hypothesized that the Grand Canyon is actually Late Cretaceous in age and was carved into the now-absent Mesozoic strata by one of the northeast-flowing drainages that are well documented on the southwest Colorado Plateau boundary. The data used to provide evidence for the Late Cretaceous paleocanyon of the California River hinges on the thermal modeling of a few AHe data samples. Luchitta et al. (2011) presented evidence for Middle Miocene drainage to the east of the Grand Canyon. This proposed drainage system flowed to the southwest, possibly crossing the Kaibab Uplift through a paleocanyon. This is consistent with other work (Flowers et al., 2008; Pederson, 2008) and new evidence (in the following) for broad middle Cenozoic planation of the region and a more modest and somewhat younger paleocanyon dissecting the Kaibab Uplift.

AHe and AFT dating are both well-established thermochronological techniques and are widely applied in geological, tectonic, and geomorphologic studies (e.g., Ehlers and Farley, 2003; Farley, 2000; Farley and Stockli, 2002; Stockli et al., 2000). AHe dating of apatite is based on the radiogenic production of helium from the alpha decay of ²³⁵U, ²³⁸U, and ²³²Th. Helium is completely diffused from apatite at temperatures above ∼80 °C and almost totally retained below ∼40 °C, depending on grain size, cooling rate, and effective uranium concentration (eU). The effect of eU on ages is particularly important and is addressed in more detail herein. AFT dating is based on the decay of trace ²³⁸U by spontaneous nuclear fission (e.g., Fleischer et al., 1975; Dumitru, 2000). The resulting fission tracks are partially or entirely annealed by thermally induced recrystallization at elevated temperatures, causing reductions in both the lengths of individual tracks and the fission track ages. Temperature and chemical composition appear to be the predominant factors controlling fission track annealing in apatite. The quantitative understanding of fission track annealing kinetics can be used to constrain the thermal evolution of samples between ∼60 and 110 °C.

Samples were collected from four age-elevation transects, four boreholes, and a lateral transect along the axis of the Colorado River. In total, 69 samples were originally collected from the vertical transects and the boreholes; 45 of these samples provided adequate apatite for single-grain AHe age dating. The 13 river samples are a subset of a larger sample set collected for fission track analysis and all but 5 have been previously published (see Kelley et al., 2001; Kelley and Karlstrom, 2012).

The four vertical transects are shown from a pseudo-aerial perspective in Figure 2. Two of these transects were collected from the area of highest relief in the Grand Canyon, along the North Kaibab Trail and South Kaibab Trail. The Kaibab transects span a vertical distance of 1573 m. On the western flank of the Kaibab Uplift, an additional vertical transect was collected on the Bob Hall and Thunder River Trails. The vertical relief of this transect was 1478 m. A final vertical transect was collected in Saddle Horse Canyon near Toroweap Overlook in the central Grand Canyon. However, the imposing topographic relief of the region made it impossible to gather samples below the Mississippian Redwall Limestone cliffs. Borehole samples from uranium exploration wells on the Coconino Plateau to the southwest of the eastern Grand Canyon were analyzed to investigate the broader exhumation patterns of the area. The shallowest and stratigraphically highest sample (BM-1; see Table 1) was at an elevation of 1619 m at a depth of 31 m beneath the Kaibab Limestone. Conversely, the deepest and lowest elevation samples come from the SBF borehole (SBF8) at a depth of 552 m and an elevation of 1258 m. The river samples were collected along the banks of the Colorado River from Lee’s Ferry in the east to Quartermaster Canyon in the western Grand Canyon. The 13 samples used in this study were collected for AFT analysis and most were published previously (Kelley et al., 2001; Kelley and Karlstrom, 2012); they have been reanalyzed for AHe ages and are presented here (Fig. 3; Table 2).

Inverse thermal modeling was used to identify cooling histories that could produce the observed thermochronometer ages, track lengths, and age-eU relationships. These thermal histories were calculated using the software package HeFTy version 1.7.4 (Ketcham, 2005). To run a thermal simulation in HeFTy, the user must input between one and six thermochronometric models. Each of these models can be an AFT model or AHe model. The simulation can use any combination of these models. In the inverse-model simulation mode, a set of ages and track length distributions are calculated for a randomly generated thermal history for each of the models provided for that simulation. Each of these modeled age sets are compared to the observed ages and lengths and a goodness of fit (GOF) is calculated. The GOF indicates the probability of failing the null hypothesis that the modeled data and the observed data are different. GOF values >0.05 are defined as acceptable agreement between modeled and observed data; values >0.5 are regarded as best fits indicating that the observed data support the modeled data. For the samples for which we have AFT data (i.e., river-level samples), all simulations include an AFT model and an AHe model. For detrital samples, each simulation includes a number of AHe models, from one to three depending on the spread in eU.

The annealing model used to model AFT ages is that of Ketcham et al. (2007). Additional input parameters for the AFT models include a chlorine weight percent of 0.10%, a default initial mean track length based on the chlorine weight percent of 16.17 um, and a track length reduction standard of 0.893. The GOF calculation method is Kuipers statistic test. The zeta mode is traditional. Table 2 outlines the observed parameters that were used to model AFT ages.

Parameters for the apatite AHe models include the following. (1) The equivalent spherical radius of the grain (or grains) calculated from measured dimensions observed under a stereomicroscope and listed in Supplemental File 1¹. (2) The alpha stopping distances defined by Farley et al. (1996). (3) The alpha calculation is set to “ejection only.” (4) An r_mr0 value of 0.83 is used as is typical for apatite-CaF. (5) The uncorrected age and associated error input is based on the averaged uncorrected age for the specific model group. For most detrital samples, this averaged age does not correspond to an average of all analysis but to the average age of the eU group (as outlined in Supplemental File 1 [see ^{footnote 1}]). The error is the standard deviation of the age group. (6) The observed group average uranium, thorium, and samarium concentrations as measured by inductively coupled plasma–mass spectrometry are given in parts per million. (7) The most recent diffusion model is used (i.e., RDAAM [radiation damage accumulation and annealing model]; Flowers et al., 2009). As a result, the AHe thermal modeling techniques used in this study account for the variability of closure temperatures caused by alpha damage.

The correlation between age distributions and uranium concentrations is well documented and was discussed in previous studies of the Grand Canyon region (Shuster et al., 2006; Flowers et al., 2009; Wernicke, 2011). This correlation is also observed in the data presented here. Because uranium and thorium concentrations are pivotal parameters in determining the diffusion kinetics of apatite, they are critically important when interpreting cooling ages. The result of a sample having a large range of eU values is that of a poorly reproducible single-grain age distribution. However, there is a large amount of information stored in the age-eU relationship. This is particularly true for samples that undergo protracted cooling through the partial retention zone (i.e., 40–80 °C). In this way, grains of different eU effectively have different closure temperatures; in these certain scenarios, modeled thermal histories can therefore be much more constrained than a sample with homogeneous eU. In this study, each sample’s age-eU relationship was evaluated to identify samples that might allow the modeling of constrained thermal histories. Those relationships are outlined in Supplemental File 1 (see ^{footnote 1}). Theoretically, it would be advantageous to model each grain individually when modeling the thermal history of a particular sample. However, the inherent variability in individual grain analysis often makes it impossible to model all grains individually and to produce a common and/or constrained thermal history. For this reason, individual analyses are grouped by natural breaks in the age-eU relationship to provide more characteristic modeling parameters. Supplemental File 1 (see ^{footnote 1}) outlines the natural breaks and groups used in this study. Each of these groups is treated as an individual thermochronometric model when generating thermal histories in HeFTy. Because the spread of eU in crystalline basement samples is often too narrow to allow age-eU grouping, the modeled thermal histories from AHe data alone are relatively unconstrained. Fortunately, these samples were initially collected for fission track analysis and we have AFT data to model with the AHe data and further constrain the possible thermal histories, particularly at temperatures between 90 and 60 °C.

For some AHe analyses, apatite grains had very high eU values relative to the other analyses in that sample. Some of these grains exhibit ages that are younger than equivalently sized lower eU grains from the sample; this is the converse of the usual association of high eU and older ages. Analyses SBF6–6, SBF7–6, Wate1810–6, and GCNK2–3 all showed this behavior (illustrated in Supplemental File 1 [see ^{footnote 1}]). We suggest that the extreme degree of alpha damage in these grains leads to the creation of fast diffusion pathways. The diffusion kinetics of these grains are poorly understood and difficult to quantify. As a result, they are eliminated from age-eU grouping and not included in modeling efforts. Another complicating scenario is one in which apatite grains are not heated enough to diffuse all helium from the grain and are thus only partially reset. The result would be grains that have an inherited helium constituent that leads to older ages, particularly so for larger grains or those with higher eU values.

Because AHe ages are a function of cooling history as well as grain size and effective uranium concentration (Shuster et al., 2006), we do not report average ages and their associated errors. Instead, ages here are generally discussed in terms of the grain-to-grain age span of a sample.

The north and south Kaibab transects were collected from the area of greatest relief in the Grand Canyon. The individual grain ages from the highest samples on the North Rim were 40.6–77.8 Ma, and from the highest samples on the South Rim were 21.1–66.7 Ma (see Supplemental File 2²). In contrast, the ages for the canyon-bottom samples of each transect had an age span of 35.7–52.2 Ma for the north and 17.8–40.5 Ma for the south. The Bill Hall–Thunder River transect is from the western flank of the Kaibab Uplift. Similar to the Kaibab transects, the sample ages were 53.0–87.1 Ma near the rim and 22.9–45.5 Ma at the river. For the Saddle Horse Canyon transect, ages from the highest sample came from an elevation of 1669 m and spanned 17.8–56.8 Ma. The lowermost sample came from an elevation of 1304 m and yielded ages of 10.5–90.6 Ma.

Borehole samples proved to be difficult to produce AHe ages from because the limited stratigraphic range in the boreholes led to a small amount of usable apatite yield. However, 15 samples from 4 boreholes yielded sufficient apatite to produce single-grain AHe ages. Individual grain ages varied widely but correspond roughly with depth and ranged from 147.8 Ma (SBF2) to 12.3 Ma (Wate1255).

AFT ages range from 73.5 ± 8.9 Ma (02GC128) west of the Hurricane fault to 28.4 ± 2.9 Ma (01GC92) at Lee’s Ferry. AFT ages exhibit a general increase toward the west (see Fig. 3). However, abrupt changes in AFT ages often correlate with structural offset of Laramide reverse faults, as discussed in Kelley et al. (2001). The AHe ages from river-level samples often exhibit less age spread than those from the vertical transects; we suspect that this is due to higher ultimate burial temperatures and the more homogeneous uranium concentrations of the crystalline rocks from the river-level samples. AHe ages from the core of the Kaibab arch ranged from 21.7 Ma (98GC16) to 55.9 Ma (98GC11). Samples from the Lee’s Ferry region yielded AHe ages that ranged from 5.6 Ma (01GC90) to 20.3 Ma (01GC93). Laramide and Sevier cooling is the main event recorded in three samples located west of the Hurricane fault in the western Grand Canyon. A robust AFT age of 62.8 ± 4.0 Ma for sample 01GC86 (river mile 243) with relatively long mean track lengths of 13.0 ± 0.4 μm coupled with 29.2–72.3 Ma AHe ages on low eU (11–17) grains constrain an episode of rapid cooling ca. 65–75 Ma, and residence at temperatures of ∼50 °C through most of Cenozoic time. The AFT age of 68.7 ± 3.8 Ma for sample 01GC87 from the Meriwitica monocline is much younger than the AHe ages of 69.5–91.9 Ma for this sample. The high uranium content of the apatite is responsible for the age inversion, as is seen in samples that reside in the partial retention zone for a long time (Flowers and Kelley, 2011). The mean track length of 12.1 ± 0.4 μm and the scatter in the AHe ages indicate that the Meriwitica area had a more protracted cooling history compared to that at river mile 243. Sample 01GC89 contained only 4 grains in the AFT sample, so this age of 63.2 ± 20.3 Ma is not robust.

Dumitru et al. (1994) presented AFT data that indicate that the Kaibab Limestone surface atop the Kaibab Uplift was heated to temperatures of ∼90–110 °C by burial and started cooling in the Late Cretaceous. In addition, AFT data from the Supai Formation on the crest of the Kaibab Uplift indicate that cooling through ∼110 °C was in progress by 90 Ma (Kelley and Karlstrom, 2012). Correspondingly, the rim samples near the Kaibab Uplift were modeled with a beginning constraint of 90–130 °C at 80–130 Ma. Because the Toroweap vertical transect is located west of the Kaibab Uplift, and because the maximum burial temperatures decreased to the west, the Saddle Horse Canyon samples have a beginning constraint of 70–120 °C at 80–130 Ma. The river-level samples are modeled with a beginning constraint of 100–140 °C at 75–80 Ma where the higher temperatures relate to deeper burial. Borehole data are modeled with a relatively wider set of constraints because no AFT data are available to help constrain maximum burial temperatures. However, because the samples come from the west of the Kaibab Uplift, and because the maximum burial temperatures decreased to the west, we expect the maximum burial temperatures of these samples to be slightly less than the samples of equivalent stratigraphic depth from the Kaibab Uplift. Beginning constraints of 70–140 °C at 100–140 Ma were used here and vary according to the depth of the sample beneath the modern surface (illustrated by the red boxes in Supplemental File 3³).

In total, 27 samples yielded either sufficient variations in eU or additional fission track data to allow rigorous modeling of sample thermal histories. These samples come from all 4 age-elevation transects, all 4 boreholes, and 11 of the river-level samples. The degree of constraint on viable thermal histories is directly proportional to the number of analysis, spread in eU, and the availability of fission track data. Figure 4 shows the best-fit (i.e., 95% confidence interval) modeled thermal histories for 16 representative and well-constrained samples from the entire study area. Major cooling episodes can be identified in the thermal histories shown in this plot. For a detailed description of modeling results, see Supplemental File 3 (see ^{footnote 3}), where the detailed modeling results and constraints used for each sample are individually listed.

The broadest and least constrained cooling event observed in almost all thermal histories is a Late Cretaceous cooling episode, which is either observed or allowable in all samples excepting those east of the Kaibab Uplift. The least constrained thermal histories come from the vertical transects and borehole samples because of the lack of fission track data that would have helped define the thermal history at temperatures above ∼70 °C. As a result, data here help resolve little in terms of the Late Cretaceous unroofing beyond that described previously (see Dumitru et al., 1994; Kelley et al., 2001; Flowers et al., 2008). However, one notable observation is that the samples near Lee’s Ferry fail to show any indication of Laramide cooling, while the westernmost sample (01GC86) indicates cooling to ∼50–60 °C by the early Cenozoic. We suggest that this observation supports a trend in decreasing magnitude of Late Cretaceous unroofing to the east, as suggested previously (Flowers et al., 2008).

A period of slow cooling rates, typically <5 °C/m.y., occurred from the early Cenozoic to the middle Cenozoic. This period of slow cooling is particularly evident in river-level samples, except those near Lee’s Ferry. River-level samples in the core of the Kaibab Uplift are between 70 and 90 °C while rim samples from the Kaibab Uplift are between 40 and 60 °C. As discussed in the following, this early to middle Cenozoic temperature gradient has important implications for the existence or absence of paleorelief in the early to middle Cenozoic time.

Following the early to middle Cenozoic period of slow cooling rates, a widespread cooling event occurred in the latest Oligocene to Early Miocene (28–20 Ma). This cooling event, as much as 75 °C in a period of 6 m.y., is particularly evident in samples from river level in the Kaibab Uplift (i.e., 98GC11, 98GC16, 98GC20, 01GC103) and occurs in virtually all samples west of Lee’s Ferry except those farthest from the modern canyon (i.e., the Wate and SBF borehole samples). During this cooling event, river-level samples from the core of the Kaibab Uplift cooled 30–50 °C, bringing them to <∼30 °C, approximately the same temperatures as the rim samples. The relative temperatures of the rim and river-level samples have important implications for the existence of paleotopography, as discussed in the following. A final cooling event is observed post–10 Ma in the borehole samples and Lee’s Ferry samples; samples cooled from 40–60 °C to near surface temperatures. This cooling event is significant because it indicates that the Lee’s Ferry samples were 20–50 °C warmer than samples of nearly equivalent elevation to the west throughout most of the Miocene.

Although most modeled samples yielded comparatively compatible results within the confines of other geologic constraints, a few samples yielded anomalous results, including some of the borehole samples. In an effort to accommodate the likelihood that samples to the west of the Kaibab Uplift underwent lower maximum burial temperatures than those from the Kaibab Uplift, the borehole samples were modeled with wider beginning constraints. The modeling results indicate that cooling from maximum temperatures commenced in the Early to middle Cretaceous, not the Late Cretaceous as in most other samples. Although it is certainly a possibility that the Coconino Plateau underwent much earlier unroofing than the nearby Kaibab Uplift, we offer the alternate possibility that some of these samples may have inherited helium due to incomplete thermal resetting, i.e., SBF6, SBF7, Wate1255, and Wate1810. Their stratigraphically high and relatively westward location makes them more likely candidates for such a scenario. We observed AHe-AFT age inversions in the two westernmost river-level samples; due to decreasing maximum burial temperatures to the west, we suspect that even these stratigraphically lower samples might be retaining inherited helium. Because of the obvious complications with these two samples, they were not modeled as the rest of the river-level samples were. Unfortunately, quantification of the amount of inherited helium is difficult, if possible, to accomplish. We strongly suspect that inherited helium is not an issue with other samples in this study because all other sample locations underwent temperatures sufficient to diffuse helium completely from the apatite grains (i.e., >90 °C).

To obtain characteristic thermal histories for particular regions of interest (e.g., river level at the core of the Kaibab Uplift), we group samples in such regions, and plot their time-temperature paths atop each other to identify common portions of their paths. This technique and its results are outlined in Figure 5A. Upon creating characteristic thermal histories for each region, we are able to easily compare thermal histories between regions of interest. As identified in Figure 5B, four important observations are outlined here.

1. The elevation difference between rim and river samples from the western flank of the Kaibab Uplift (brown and red lines in Fig. 5) is 1.5 km and the observed temperature difference is ∼38 °C for a calculated thermal gradient of 25 °C/km in the middle Cenozoic. This value closely matches the modern geothermal gradient and supports the validity of the modeled thermal histories during that time. The observed thermal offset of the rim and river-level samples is evidence that argues against the existence of large topographic relief above these samples throughout the middle Cenozoic (discussed in more detail in the following).

2. A 28–20 Ma cooling event is observed in all sample groups except those few samples farthest from the modern canyon (see Fig. 5). In this cooling event, the deepest river-level samples are rapidly cooled to near-surface temperatures and the temperature offset between rim and river-level samples of ∼38 °C is lost. We interpret this to identify the creation of topographic relief through partial incision of a paleocanyon in the same location as the eastern Grand Canyon.

3. Figure 5 illustrates that the only sample group that does not exhibit rapid cooling during the Early Miocene comprises the deeper borehole samples, which are farther away from the modern canyon. Because the effect of relief generation on isotherm deflection decreases with distance from the relief and increases with the magnitude of relief, we suggest that these samples were not affected because the proposed Miocene canyon was of smaller scale than the modern canyon and perhaps did not exist at all on the Coconino Plateau. The eventual rapid cooling of the deeper borehole samples occurred after 10 Ma. For the reasons just stated, this observation is attributed to final integration of the Colorado River ca. 5.5 Ma and reinitiated incision leading to increased relief.

An additional way to synthesize the numerous modeled thermal histories is to extract the modeled temperature range for each sample at specific time intervals. It is then possible to plot the elevation of each sample against the extracted thermal range. These temperature-elevation plots allow insight into the structural and topographic evolution of the entire modeled sample set through time. For example, in a hypothetical scenario A (Fig. 6A), where samples underwent cooling in a setting of low-relief erosional exhumation (i.e., horizontal isotherms), one would expect all samples to be on a linear trend of increasing temperatures with depth. The slope of this line would approximate the geothermal gradient during the time interval of sample cooling (paleogeothermal gradients). In contrast, if the samples were hypothesized to cool in an environment of dissection and high canyon relief (Fig. 6B), then the samples close to the canyon would plot off of the linear trend toward older ages because the downward deflected isotherms in the area of the canyon would cause samples of equivalent elevation to cool sooner. Accordingly, if samples underwent differential structural uplift during cooling, then the linear trend would be offset by a calculable amount (Fig. 6C). With these hypothetical scenarios, we are able to test for significant paleorelief and/or large magnitude structural deformation through the Cenozoic.

Time-temperature plots were produced for the modeled samples from the Grand Canyon; Figure 7 shows selected time intervals that are useful in constraining the landscape evolution of the region. The earliest time slice (at 70 Ma) illustrates interesting departures from an overall linear trend of the data set that occur in the samples from the Kaibab Uplift. These samples resided at relatively higher temperatures until 60–50 Ma, at which time they cooled to join the average linear trend shared by all other samples of the central Grand Canyon region (Fig. 7). As predicted in Figure 6C, this observation is consistent with uplift and subsequent unroofing of the Kaibab Uplift in the late Laramide. Previous to formation of the monocline, which currently exhibits structural relief of ∼1 km, the samples would have resided ∼1 km deeper (and therefore 20–30 °C warmer) prior to uplift and unroofing of the monocline. We find that this estimation is in good agreement with the observed thermal histories and therefore loosely constrains the formation and subsequent unroofing of the Kaibab Uplift through the middle Paleocene to the Early Eocene.

The relations between the thermal histories of the age-elevation transect samples and the borehole samples from the Coconino Plateau are particularly informative. As shown in Figure 7, borehole and vertical transect samples are on the same linear trend from 50 Ma until 28–20 Ma. As outlined in Figure 6A, this observation is additional evidence that low-relief planation, rather than significant paleocanyon relief (cf. Wernicke, 2011), characterized this landscape from 80 Ma until 28–20 Ma. This is supported by the additional observation that the calculated paleogeothermal gradient (i.e., slope of the linear trend) closely matches the modern gradient of 25 °C/km.

Samples from the Lee’s Ferry to the east of the Grand Canyon also plot off of the average linear trend observed in most samples throughout the Cenozoic (Fig. 7). Consistent with previous thermochronologic works (Kelley et al., 2001; Flowers et al., 2008), Lee’s Ferry samples were at higher temperatures throughout the Tertiary as a result of wedge-like beveling of the Mesozoic strata in the early and middle Cenozoic and late unroofing late in the Cenozoic.

A final observation from the temperature-elevation plots is the post–20 Ma deviation toward higher temperatures indicated by the deeper borehole samples (Fig. 7). The thermal modeling indicates that these samples resided at a higher temperature than samples of nearly equivalent elevation nearer the Grand Canyon. For the same reason that the linearity of ages observed in the early Cenozoic argue against paleorelief, the large discrepancy in modeled temperatures after 20 Ma and departure from linearity argue for the creation of paleorelief in the form of an ancestral eastern Grand Canyon between 28 and 20 Ma.

In order to produce graphical unroofing cross sections along the course of the modern Colorado River, we calculated the depth of each river-level sample at selected times. Using the extracted thermal ranges from the temperature-elevation plots of Figure 7 (and Supplemental File 3 [see ^{footnote 3}]), one can calculate the amount of overburden above a particular sample when a geothermal gradient and average surface temperature are assumed. A geothermal gradient of 25 °C/km and average surface temperature of 10 °C are used here, as in similar studies and supported by results here (Flowers et al., 2008; Wernicke, 2011). To avoid elevation-related complications due to the effects of regional and local uplift and subsidence, absolute sample elevations are avoided in overburden calculations. Rather, stratigraphic depths are normalized to the Kaibab Limestone (KND, Kaibab normalized depth) and strata are reconstructed where the Kaibab Limestone is absent above the sample, as shown in Figure 8.

The result of this technique is illustrated in Figure 9. Figures 9A–9C illustrate how the region of the central and eastern Grand Canyon underwent low-relief planation from early in the Cenozoic until 28–20 Ma, despite the late Laramide monoclinal flexure of the Kaibab Uplift. At 28–20 Ma, considerable relief formed in the eastern Grand Canyon, specifically across the Kaibab Uplift (Fig. 9D). This paleocanyon is required because the calculated burial depth for the samples is less than the stratigraphic depth below the Kaibab Limestone (the modern rim surface). Although modeling efforts here are unable to define precisely the depth of this paleocanyon, the errors used here constrain the scale of the paleocanyon to 500–1000 m, about half the relief of the modern canyon here and equivalent to the structural offset of the Kaibab Uplift. However, samples to the east (01GC-90, 01GC-92, 01GC-93) of the Kaibab Uplift indicate that 1.7–2.2 km of Mesozoic strata persisted in the area of Lee’s Ferry. As late as 10 Ma, these strata remained in the area of Lee’s Ferry (Fig. 9E).

In summary, these results indicate that the middle Cenozoic Grand Canyon landscape, prior to the ca. 6 Ma integration of the Colorado River, was characterized by (1) the removal of all discernible Mesozoic strata along the axis of the Colorado River; (2) exposure of the Kaibab Limestone along the Kaibab Uplift as well as a localized, 0.5–1.0-km-deep canyon cut across it; (3) escarpment retreat proceeded to the northeast, but a 1.7–2.2 km stack of Mesozoic strata remained in the area of Lee’s Ferry.

The unroofing history outlined in this study offers a new level of detailed constraints for the landscape evolution of the eastern Grand Canyon region. Thermochronologic ages and their modeled thermal histories reveal spatial and temporal variability in erosion through time, characterized by three major erosional episodes. (1) Late Cretaceous to early Cenozoic unroofing removed a majority of the Mesozoic strata on the southwestern margin of the Colorado Plateau during uplift and subsequent beveling and northeast-directed paleocanyon formation, as established previously (Young, 2001; Dumitru et al., 1994; Flowers et al., 2008). In contrast, the topography in the neighboring central and eastern Grand Canyon region was one of one of relatively low relief throughout the early to middle Cenozoic, and evidence argues against the presence of an early Cenozoic eastern Grand Canyon. (2) A major exhumation episode is evident in the latest Oligocene to Early Miocene (28–20 Ma), as clearly recorded in samples from the core of the Kaibab Uplift. In order to explain this modeled cooling history, erosion of a canyon through the Kaibab Uplift is required. Other paleogeographic constraints strongly suggest that the uplift was cut by a west-flowing drainage (Pederson, 2008; Lucchitta et al., 2011). However, as escarpment retreat proceeded to the northeast across the area of Marble Canyon, a package of Mesozoic strata 1.7–2.2 km thick remained around the area of Lee’s Ferry until as late as 10 Ma. (3) After 10 Ma, Marble Canyon and Lee’s Ferry samples show more than 1 km of unroofing, probably shaped by the late integration of the Colorado River. This follows the general trend of erosion and landscape evolution being increasingly youthful and active across the region toward the central Colorado Plateau.

Funding for this project was provided by the University of Kansas Department of Geology and National Science Foundation grants EAR-0408787 and EAR-0414817. We thank Eugene Szymanski, Alec Waggoner, Kwan Ye Chen, Steve Sloan, and Randy Ackerman for helping to haul hundreds of pounds of rocks many miles through the Grand Canyon, and Becky Flowers and an anonymous reviewer for their thoughtful reviews of this manuscript.