Cenozoic incision history of the Little Colorado River: Its role in carving Grand Canyon and onset of rapid incision in the past ca. 2 Ma in the Colorado River System Open Access

Little Colorado River (LCR) is one of the largest drainage basin area tributaries to Colorado River (CR). As shown in Figure 1, its headwaters are in the White Mountains and Springerville volcanic field at the southern edge of the Colorado Plateau. It flows in a broad valley across highly erodible Mesozoic strata of the Holbrook and Winslow areas, enters a deep bedrock gorge in Paleozoic rocks near Cameron, and has its confluence with the CR in Grand Canyon 62 river miles below Lees Ferry (Stevens, 1983). Its modern river profile is plotted along with the CR profile in Figure 2. The LCR profile departs from a concave-up “equilibrium” profile in its lower and upper reaches; the lower 75 km has a convex-up steepened reach in resistant Paleozoic rocks below a knickpoint (LCR knickpoint). It also has sharp bedrock knickpoints including a ca. 20 ka basalt flow at Grand Falls (Duffield et al., 2006) and several knickpoints in basalts in the headwater regions.

The overall goal of this study is to understand linkages between the Little Colorado River and Colorado River systems and their influence on carving Grand Canyon and to test models for landscape evolution in the region. We synthesize new and published data on the incisional history of the LCR region over the past 70 Ma and present the following data sets. (1) HeFTy modeling of low-temperature thermochronology data (from Flowers et al., 2008) provides constraints on the 70–15 Ma denudation episodes. (2) Data on the Bidahochi Formation and associated volcanic rocks (Dallegge et al., 2001; Dickinson, 2013) provide a record of the 16–6 Ma time interval. (3) Incision rate data over the past 7 Ma are derived from ⁴⁰Ar/³⁹Ar plus older K-Ar ages of basalt surfaces that form inverted topography and/or overlie gravel terraces in the LCR and its tributaries (Damon and Spencer, 2001; Holm, 2001a). (4) We also add ⁴⁰Ar/³⁹Ar detrital-sanidine dates on elevated surfaces (Cooley et al., 1969) and paleoriver deposits (Hereford et al., 2016). (5) For the past 1 Ma, we report new ⁴⁰Ar/³⁹Ar dates on basalts that entered the paleo-LCR channel and U-series dates on travertine-cemented CR and LCR gravel terraces (Embid, 2009; Crow et al., 2014).

Modern elevations are used to reference the relative landscape position of constraint points, and the modern LCR river profile is used for estimating incision magnitudes. However, it is probable that both relative and absolute elevations have changed due to differential surface uplift (Hunt, 1969, p. 113) resulting from a combination of changes in mantle buoyancy (Karlstrom et al., 2008, 2011; Crow et al., 2014), fault dampened incision (Pederson et al., 2002; Karlstrom et al., 2007), and isostatic rebound due to differential erosion (Lazear et al., 2013).

Powell (1879) had little scientific comment on the LCR. Dutton (1882, p. 204) thought it was as “old as the Colorado River itself” and viewed it as early Tertiary.

Early workers interpreted the erosional history of the SW Colorado Plateau to have involved Colorado Plateau planation followed by entrenchment of the Colorado River and Grand Canyon incision (Dutton, 1882; Davis, 1901; see summary in Ranney, 2014). However, these various erosional styles and processes (regional denudation, cliff retreat, and canyon incision) are all entwined and have been driven by major tectonic base-level fall and/or headwater uplift episodes that reorganized drainages. These episodes were the 70–50 Ma Laramide orogeny, 35–25 Ma ignimbrite flare-up and related caldera eruptions, 25–15 Ma extension of the Basin and Range, and 6–5 Ma integration of the Colorado River (Karlstrom et al., 2011).

Cooley et al. (1969) applied and supported the peneplain concept that stepped surfaces had formed during periods of erosion that beveled the landscape (Great Denudation of Dutton, 1882; Plateau Cycle of Davis, 1901; Robinson, 1907; Gregory, 1947; Childs, 1948). These workers described numerous regional surfaces, some gravel or basalt capped. Highest surfaces that were interpreted to predate Grand Canyon incision were named the Valencia surface (pre-lower Bidahochi now known to be pre–16 Ma) and Hopi Buttes–Zuni surface (pre-upper Bidahochi, hence pre–6–8 Ma). Surfaces interpreted to have developed during Grand Canyon incision include the Black Point surfaces (terraces averaging 137 and 213 m above the modern LCR) and Wupatki surfaces (terraces averaging ∼9, 15, 26, 54, and 76 m above the LCR). Cooley et al.’s (1969) surfaces and terraces are listed in Table 1 with new age control from this paper.

Wilson (in McKee et al., 1967), writing the consensus hypothesis of the 1964 Flagstaff meeting at the Museum of Northern Arizona, implicated the LCR to reconcile two ideas. (1) Most workers envisioned a Miocene ancestral Colorado River that originated in the Rocky Mountains and flowed to southern Utah lakes and perhaps into northern Arizona (Hunt, 1969, his fig. 69). But, (2) the Muddy Creek Formation constraint indicated that the Colorado River could not have flowed across Grand Wash Cliffs in the location of the modern Grand Canyon until after 5–6 Ma (Blackwelder, 1934, p. 557; Longwell, 1946, p. 829; Lucchitta, 1966). The Muddy Creek constraint has been questioned (Wernicke, 2011, p. 1306; Flowers et al., 2012; Flowers and Farley, 2013) but remains supported by most workers (Lucchitta, 2013a; Karlstrom et al., 2014). To reconcile these apparent conflicts, McKee et al. (1967) routed an ancestral Colorado River south from Utah, along the east side of the Kaibab uplift, to flow SE out the LCR valley, possibly to the Rio Grande drainage (their fig. 18). However, Rio Grande connections have not been supported by subsequent work nor have ancestral CR gravels been recognized in the Bidahochi Formation (Dallegge et al., 2001), nor in the 4 Ma to 55 Ma sedimentary record of the Springerville area (Cather et al., 1994). The concept of a Miocene southerly exit of an ancestral CR to the Salt River was revived by Potochnik (2011) but was not supported by Dickinson (2013) nor Karlstrom et al. (2014) based on a lack of recognized CR deposits and on the evidence for a long-lived NE paleoslope for the rim gravels (Cooley and Davidson, 1968). This Mogollon slope persists for tributaries to the LCR today (Fig. 1; Condit, 1991; Holm, 2001a).

During lower Bidahochi time (now known to be 16–13 Ma; Dallegge et al., 2001), the McKee et al. (1967) synthesis proposed that drainage from the S-flowing ancestral upper CR terminated in an internally drained Bidahochi basin (their fig. 21). The term Hopi Lake was used for this lake by Williams (1936), Repenning and Page (1956), Shoemaker et al. (1961), Sutton (1974), Scarborough (2001), and Dickinson (2013). West of the Kaibab uplift, McKee et al. (1967) proposed that an independent W-flowing “Hualapai drainage system” was established on the Hualapai plateau in pre-Miocene time (their fig. 20) with flow paths that were controlled by active faults and monoclines and potential overall northward flow toward Claron and Uinta basins (Davis et al., 2010; Dickinson et al., 2012; Karlstrom et al., 2014).

Integration of the modern Colorado River system was proposed by McKee et al. (1967) to have involved post-early Pliocene drainage reversal of the previously S-flowing LCR as it was captured by a headwardly eroding Hualapai drainage system across the East Kaibab uplift (their fig. 22; but note that, in 1967, the Miocene–Pliocene boundary was thought to be ca. 12 Ma). Cooley and Akers (1961) and Hunt (1969) proposed a pre-Bidahochi LCR that crossed the Kaibab uplift south of Grand Canyon (as also supported by Karlstrom et al., 2014 and in this paper), but Babenroth and Strahler (1945), Strahler (1948), McKee et al. (1967), and Lucchitta (1988) envisioned breeching of the uplift in Pliocene time by an ancestral CR, not a paleo-LCR. Proposed CR integration mechanisms have been diverse: headward erosion (McKee et al., 1967), lake spillover from upstream lakes (Blackwelder, 1934; Scarborough, 2001), karst piping (Hunt, 1969; Hill et al., 2008; Pederson, 2008; Hill and Polyak, 2010, 2014), groundwater sapping (Crossey et al., 2015), and linkage of older paleocanyons (Karlstrom et al., 2014). Aspects of all of these mechanisms probably operated during 6–5 Ma integration of the CR through the Grand Canyon region.

Wernicke (2011, his fig. 9A) depicted the LCR as part of an E-flowing 70 Ma “California” river system that originated in the Sierra Nevada arc and flowed through Grand Canyon, then southeast out the LCR to the San Juan Basin. At 55 Ma, in his model, a divide formed in the middle of the modern LCR valley, and rivers flowed west in an “Arizona River” through a nearly completely carved Grand Canyon and east to the Eocene Baca basin and the Gulf of Mexico. This model was supported by Flowers and Farley (2012, p. 3; 2013, p. 143c). However, there is no evidence for a paleo-LCR that had a role in either of these proposed pre–50 Ma river pathways. Instead, as documented here and by Flowers et al. (2008), thermochronology suggests that Moenkopi Formation rocks currently exposed in the LCR valley were buried beneath 2–3 km of Mesozoic strata until 25–15 Ma. Also, Marble Canyon was beneath 2–3 km of strata until 6–4 Ma integration of the CR (Flowers et al., 2008; Lee et al., 2013; Karlstrom et al., 2014).

In contrast to models for 70–30 Ma paleocanyons that followed the same path as modern Grand Canyon (Wernicke, 2011; Flowers and Farley, 2012, 2013), others suggest that the Laramide drainage system involved rivers that flowed generally north across the Grand Canyon–Colorado Plateau region (Fig. 1, McKee et al., 1967; Young, 2001; Karlstrom et al., 2013a, 2014). A fragmentary record of the 70–30 Ma landscape components and paleoriver deposits is preserved as coarse gravels that rim the southern Colorado Plateau and contain unroofed basement lithologies from the Kingman and Mogollon uplifts. These gravels are referred to as the Rim gravels (Cooley and Davidson, 1963; Peirce, 1984), Music Mountain Formation (Young, 2001), and Mogollon Rim Formation (Potochnik, 1989, 2001; Potochnik and Faulds, 1998). They are now known to be of disparate age. In the western Grand Canyon region, the Music Mountain Formation (Young, 2008) is 65–55 Ma (Young and Hartman, 2014; Hill et al., 2016) and is overlain by the Buck and Doe conglomerate that has a 24 Ma ash in its upper parts (Young and Crow, 2014). In eastern Arizona and western New Mexico, the Eagar Formation and the correlative Baca Formation are 50–40 Ma (Cather, 2004). These are overlain by the 37–34 Ma Mogollon Rim Formation on the Colorado Plateau, which is correlated with the 38–20 Ma Whitetail Formation in the Salt River Canyon (Potochnik, 2001). For Colorado Plateau Rim gravels, local base level remained at similar elevation from 65 to 25 Ma, rivers were sourced in the contemporaneously uplifting Mogollon highlands, regional paleoslopes were generally toward the north and east, and flow was initially toward the retreating Cretaceous seaway (Potochnik, 1989, 2001).

The LCR valley has not been proposed as an important pathway for Laramide and post-Laramide drainage (except Wernicke, 2011); instead most workers envision a broad alluvial plain with major rivers potentially reaching San Juan and Baca basins (Potochnik, 2001; Cather, 2004) and/or Uinta Basin (Davis et al., 2010; Dickinson et al., 2012). Paleocanyons developed 70–30 Ma such as a ∼1-km-deep Music Mountain paleovalley (Young, 2008) and 1000–1400-m-deep Laramide Salt River paleocanyon (Potochnik and Faulds, 1998). Additional thermochronologic evidence for km-scale of Laramide paleotopography is discussed below but is related to retreating cliffs of Mesozoic strata rather than to paleocanyons. Paleocene to Oligocene times involved aggradation in the southern Colorado Plateau. Similar to the “turnaround” from an aggradational to an erosional landscape in the Rocky Mountains (McMillan et al., 2006), this transition took place in some locations after 25 Ma and after 10–6 Ma in the southern Colorado Plateau (Karlstrom et al., 2011).

The region of the confluence of the CR and LCR is of key importance for river evolution models. Lucchitta et al. (2011, 2013b) proposed that a Miocene (?) paleoriver that they named the Crooked Ridge paleoriver carved across the East Kaibab uplift, a model now refuted by ⁴⁰Ar/³⁹Ar sanidine ages of 2 Ma on these paleoriver deposits (Hereford et al., 2016). Figure 3 is a photograph of the CR-LCR confluence looking east. The CR (green water) is dominantly snowmelt from the Rockies; the LCR (turquoise blue water) is carbonic spring water from Blue Springs (Crossey et al., 2009, 2015). To visualize the problem, the white dashed lines that bottom at 1600 and 1200 m elevation show speculative depths of the East Kaibab paleovalley at the ca. 6 Ma time of CR integration (Lee et al., 2013; Karlstrom et al., 2014). Alternatively, Cape Solitude (CS, 1830 m) has been proposed as point of spillover of Hopi Lake to drive CR integration (Scarborough, 2001). The 1145 m and 985 m lines show the elevation of the 2 and 1 Ma CR-LCR confluence based on steady incision rate data (Crow et al., 2014). The red line in the distance shows the Gap, a paleocanyon-shaped “wind gap” in the Jurassic Navajo Sandstone along the Echo Cliffs monocline, and the dashed red lines show the farthest west outcrops of the ca. 2 Ma White Mesa alluvium (Hereford et al., 2016; also referred to as the Crooked Ridge paleoriver deposits by Lucchitta et al., 2011, 2013b). Big Creek is a dry tributary to the LCR.

This paper integrates a variety of methodologies to understand the Cenozoic landscape evolution of the LCR valley.

Bedrock incision estimates were made using strath heights of terrace gravels or base of basalt flows that had channel-like shape suggesting they flowed in the LCR or a tributary paleochannel. Incision rate data points are listed in Table 1. Description of bedrock incision rate calculation methods, individual incision points of this study, and their uncertainties are in section DR-1 of the Supplemental Materials (see ^{footnote 1}). For samples distant from the LCR mainstem, we projected the height of the strath or flow base to the closest tributary and reported them as tributary incision rates. In the headwater regions, we digitized the elevation of bases and tops of basalt flows, gravel terraces mapped by Sirrine (1958), and travertine straths and projected them perpendicular to the river to estimate incision rates (Embid, 2009).

This paper presents thermal history models of apatite (U-Th)/He (AHe) thermochronologic data from Flowers et al. (2008) using the program HeFTy 1.7.4 (Ketcham , 2005). Constrained thermal paths are generated by inverse modeling using the radiation damage accumulation and annealing model diffusion model (RDAAM) (Shuster et al., 2006; Flowers et al., 2009). Section DR-2 in the Supplemental Materials (see ^{footnote 1}) lists sample locations and information. HeFTy takes predicted results from time-temperature (t-T) paths, compares them to observed data, and calculates a goodness of fit (GOF). Goodness of fit indicates the probability of failing the null hypothesis that the modeled data and the observed data are statistically different. High values of GOF indicate a high probability that the null hypothesis can be rejected and there is hence a good match to the measured data. Goodness of fit values >0.05 are defined as acceptable agreement between modeled and observed data; values >0.5 are regarded as good fits (Ketcham, 2005). A temperature of 15 ± 10 °C is used to represent near-surface conditions, which encompass the 25 °C maximum surface temperature used for conversions of temperature to depth in this paper as well as probably more realistic 10 °C surface temperatures.

⁴⁰Ar/³⁹Ar dating was conducted by the New Mexico Geochronology Center. Analytical methods and results are reported in section DR-3 of the Supplemental Materials (see ^{footnote 1}). These build on the data presented in Hereford et al. (2016).

Uranium-series analysis of travertine was done on a Neptune inductively coupled plasma mass spectrometer (ICPMS) in the Radiogenic Isotope Laboratory at the University of New Mexico to establish absolute, high-resolution ages of periods of travertine deposition. Calculations used the following decay constants: λ²³⁰ of 9.1577 ± 0.0278 × 10^–6 per year and λ²³⁴ of 2.8263 ± 0.0057 × 10^–6 per year (Cheng et al., 2000; note dates were not recalculated according to Cheng et al., 2013). Errors reported are 2σ confidence intervals. Two laboratory blanks, analyzed for quality assurance, had measured values of <30 pg of ²³²Th and ²³⁸U. Ages were corrected using an initial ²³⁰Th/²³²Th activity ratio of 0.8 ± 0.4. Section DR-4 in the Supplemental Materials (see ^{footnote 1}) reports the analytical data.

A geologic map of the LCR valley is shown in Figure 4. The LCR region is a broad strike valley cut mostly on the Triassic Moenkopi Formation. Precambrian basement is exposed in the Mogollon highlands to the southwest, and progressively younger Phanerozoic strata are exposed to the north and northeast. The now-eroded thickness of Jurassic and Cretaceous strata that once covered the LCR region is not well known, but estimates can be made based on nearby sections. As shown in Figure 5, Upper Cretaceous strata in southern Utah, just north of Lees Ferry (Doelling and Willis, 2006), have a mean thickness of 1667 m (range 1194–2142 m). This is comparable to the thickness of 1785 m (1325–2245 m) for the Kaiparowits Plateau (Tindall et al., 2010). An incomplete section of 511 m (312–710 m) is preserved on Black Mesa (Nations et al., 1995). Jurassic strata are estimated as 1107 m thick (619–1492 m) in southern Utah (Doelling and Willis, 2006) to 1680 m (960–2400 m) exposed in the Vermillion cliffs near Lees Ferry (Billingsley, 1989). Thus, based on stratigraphic thickness to the north, presently exposed lower Triassic rocks of the LCR valley may have been buried by up to 3–3.5 km of Mesozoic strata consisting of 1.7 km of upper Cretaceous rock, 1.1–1.7 km of Jurassic rock, and several hundred meters of Triassic strata. Assuming a geothermal gradient of 25 °C/km and a surface temperature range of 5–25 °C, these stratigraphic data suggest that lower Triassic rocks of the LCR valley would have resided at temperatures of 80–110 °C at 85–70 Ma. This is in agreement with the ≥80 °C temperatures suggested by Flowers et al. (2008) and with the 80–120 °C temperatures that best fit the thermochronologic data using the RDAAM model (this paper).

However, this 3–3.5 km pre-Laramide depth-of-burial for the LCR samples indicated by the thermochronologic data contradicts Molenaar (1983), who suggested that Cretaceous rocks tapered to a feather edge along the southern Colorado Plateau margin. The original thickness of Jurassic strata in the study area is unknown and is complicated by a sub–Dakota Formation unconformity with progressively thinner earlier Mesozoic rocks to the south that indicates that most or all of the Triassic and Jurassic section may have been denuded prior to ca. 94 Ma (Dickinson et al., 1989; Potochnik, 2001, his fig. 4), at least in areas near the New Mexico border.

The distribution of AHe thermochronology samples shown in Figure 1 includes 14 published AHe samples (Flowers et al., 2008). All but one sample is from the lower Moenkopi Formation of Triassic age, which has a depositional age of 247–241 Ma (Dickinson and Gehrels, 2009). Sample M is from Permian strata (Fig. 4). The samples can be grouped into three cross-valley NE-SW transects: Cameron transect closest to Grand Canyon (samples A, E, F, G, J, and N), Winslow transect (samples B, I, and L), and Holbrook transect (samples C, D, H, K, and M). All three transects have samples that extend SW from near the LCR level up the ∼0.5°–2°, NE-dipping Mogollon topographic slope (Fig. 1). The present LCR valley is broad and low relief, with total vertical relief between samples of 700 m over a lateral distance of 100 km (Fig. 1). The apatite grains were deposited as detrital sedimentary grains in the Mesozoic sandstones and hence have crystallization ages older than the depositional age of the Moenkopi Formation. Based on analyzed detrital zircons from this formation (Anderson, 2006), >50% of the apatites are also likely to have formed in the Precambrian. AHe ages are all <0–70 Ma, suggesting that cooling took place due to denudation following the 80–75 Ma onset of the Laramide orogeny (Goldstrand, 1994; Tindall et al., 2010).

Between four and ten apatite grains were analyzed from each sample and yielded a wide range of AHe ages within individual samples (Table DR-2.1 in the Supplemental Materials [see ^{footnote 1}]). This can explained by the RDAAM model (Shuster et al., 2006; Flowers et al., 2009), where crystals develop varying amounts of lattice damage depending on their effective uranium concentration (eU = [U] + 0.235[Th]). Higher eU causes increased lattice damage, higher closure temperatures, and hence older ages. In this data set, the correlation between age and eU for the grains as a whole is weak (R² = 0.13; Fig. DR-2.2B [see ^{footnote 1}]) but has a positive slope as predicted by the RDAAM model.

Flowers et al. (2008, their fig. 7) interpreted the data from these samples in terms of zones of diachronous cooling of the Kaibab surface below 45 °C, with 60–50 Ma for samples I–M high on the Mogollon slope above 1900 m elevation; 50–35 Ma for sample H; 28–18 Ma for samples A–G near the level of the LCR; and <6 Ma for the Marble Canyon region. These zones agree with the distribution of the Rim gravels and the Bidahochi Formation.

To better evaluate the cooling history of each sample and implications for denudation, we applied the RDAAM model using HeFTy to generate cooling paths for each sample. Our modeling approach involved a number of steps that were applied systematically to all samples. Geologic constraint boxes were imposed on the models as summarized in section DR-2.4 and Figure DR-2.6 (see ^{footnote 1}). These included a Precambrian crystallization age of most of the apatite grains (Anderson, 2006), early Triassic deposition (250–230 Ma), heating to a pre-Laramide (90–80 Ma) temperature range of 40–140 °C due to burial by an unknown thickness of Mesozoic strata, then cooling to modern surface temperatures of 10–25 °C. The pre-Laramide portion of the path allows for buildup of radiation damage in crystals that reside above or within the AHe partial retention zone for extended time periods (Fox and Shuster, 2014). Using these geologic constraints, and with a broad assigned error of ±20% for corrected AHe ages, HeFTy time-temperature paths were generated that predict most of the dated grains constraints according to the RDAAM model (Table DR-2.1 [see ^{footnote 1}]) with a statistical GOF >50% (Ketcham, 2005). Samples in proximity to Rim gravels must have been at near-surface temperatures in the mid-Tertiary, and these were run with and without an appropriate constraint box (Fig. DR-2.6 in the Supplemental Materials [see ^{footnote 1}]). Similarly, samples in proximity to Lower Bidahochi Formation must have been at near-surface temperatures by 16 Ma, and these were run with and without this constraint box (Fig. DR-2.7 [see ^{footnote 1}]).

Details of the HeFTy modeling are described in section DR-2 with a page per sample in section DR-2.5 that shows the full time-temperature path (from the Precambrian), the post–100 Ma blow-up of the cooling history, the number of path attempts needed to generate 100 good paths (GOF = 0.5), the number of acceptable paths (GOF = 0.05), and the eU-age plot identifying the grains that were modeled together. To evaluate how well the best-fit time-temperature path predicts the AHe data, we show its GOF for each grain constraint. All time-temperature models are based on four to seven analyzed apatite grains; half of the samples generated successful HeFTy models using all analyzed grains. For the others, one to three grains were excluded empirically in the process of trying to generate acceptable and good-fit time-temperature paths (Table DR-2.1 in the Supplemental Materials [see ^{footnote 1}]).

Figure 6 shows results of models for the southeastern (Holbrook) transect. Samples C, D, H, K, and M extend from 1668 to 2025 m elevation, and the different cooling paths are summarized in the bottom panel by the HeFTy-generated weighted mean paths. Note that we show these smoother weighted mean paths rather than the more jagged best-fit paths for comparison, but the weighted mean path generally predicts the data less well than the best-fit t-T path (also shown in each model). All of these samples yielded multigrain “good” (GOF >0.5) time-temperature paths using a majority of apatite grains (C = 5/5, D = 3/4, H = 5/7, K = 5/5, M = 4/4; Table DR-2.1 [see ^{footnote 1}]). We show our preferred t-T models: with the Bidahochi constraint imposed for C and D, and the Rim gravel constraint imposed for H, but Figures DR-2.6 and -2.7 (see ^{footnote 1}) show that HeFTy models run with and without these geologic constraint boxes are very similar. In this transect, lower elevation samples (also the samples farther northeast) remained hotter than higher elevation (southwest) samples until all paths converge at ca. 18 Ma. The models show that near–river-level samples (C and D) resided at post-Laramide temperatures of 80–120 °C from 70 to 30 Ma. The mean-path temperature of river-level versus Mogollon rim samples between 50 and 30 Ma varies by as much as 80 °C (Fig. 6, lower right), which far exceeds the expected ∼17 °C difference predicted by relative elevations. This suggests that cooling took place from SW to NE (Flowers et al., 2008) and that the relatively sharp but diachronous cooling (M, K, and H versus C and D) suggests cooling was influenced strongly by cliff retreat. River-level samples (C and D) cooled from ∼80–30 °C between 25 and 15 Ma, and converged with modeled temperatures of higher elevation samples by 16 Ma, in agreement with the Bidahochi geologic constraint that a broad paleovalley had formed by 15 Ma.

Figure 7 shows HeFTy-generated, time-temperature paths for the Winslow transect. Samples B, I, and L extend from 1519 m to 2122 m. Paths shown for these samples were generated using most grains that modeled successfully together in HeFTy (B = 7/8, I = 7/7, L = 4/4; Table DR-2.1 [see ^{footnote 1}]). We use the Bidahochi constraint for B and the Rim gravel constraint for I, but Figures DR-2.6 and -2.7 (see ^{footnote 1}) show that these geologic constraint boxes do not markedly change the overall t-T paths. Similar to the Holbrook transect, these samples show that the lower elevation (NE) river-level sample (sample B; Fig. 7) resided at higher temperature than the higher elevation (SW) samples. The best-fit path shows rapid cooling of the river-level sample at ca. 20 Ma, whereas the highest elevation sample had cooled through ∼40 °C before ca. 50 Ma. The intermediate elevation sample (I) shows an intermediate post-Laramide temperature of residence consistent with NE cliff retreat.

Figure 8 shows t-T paths for the northernmost (Cameron) transect, nearest Grand Canyon. Cameron transect samples A, E, F, G, and J (plus N in section DR-2.5N) extend from 1417 m elevation (sample A) to part way up the north side of the Gray Mountain uplift (sample E, 1910 m), which is a southern continuation of the East Kaibab uplift. Samples J and N are located on the south Kaibab Plateau north of Flagstaff (Fig. 1). All of the samples in this transect were modeled in HeFTy with greater than half of the dated grains (A = 4/6, E = 5/8, F = 5/8, G = 5/9, J = 6/10; Table DR-2.1 [see ^{footnote 1}]). We use the Bidahochi constraint for A, but the rest are shown without constraint boxes (Figs. DR-2.6 and DR-2.7 [see ^{footnote 1}] show that these geologic constraint boxes do not markedly change the overall t-T paths). This transect shows the same pattern as the other transects with samples J and N cooling in the Laramide and the rest cooling rapidly from 70 °C to near-surface temperatures after 30–20 Ma.

AHe thermochronology constraints on the denudational history of the southwest side of the LCR valley for the area of the Winslow and Holbrook transects are shown on the left side of Figure 9. Geologic constraints on evolution of 16–6 Ma paleosurfaces preserved by the Bidahochi Formation on the northeast side of the LCR valley are shown on the right side. Each data set can be compared and perhaps calibrated by the other. Shoemaker et al. (1961), Cooley et al. (1969), and Cather et al. (2008) described a major denudational event that postdated the 35–27 Ma aggradation of the Chuska Sandstone in the Four Corners region and predated post–6 Ma CR integration (red arrow of Fig. 9). Cather et al. (2008) estimated average rates of 119 m/Ma for this denudational event. This denudation lasted until the 15.8–13.7 Ma deposition of the lacustrine lower Bidahochi Formation at modern elevations of 1770–1950 m (orange in Fig. 9). In the headwater region of the LCR, this denudation produced a Late Oligocene to Early Miocene regional angular unconformity with the exhumation of at least 1230 m of section by ca. 16 Ma (Cather et al., 2008). This was later followed by a hiatus in the record (non-deposition?) from 13.7 to 8 Ma, then by 8–6 Ma aggradation of the upper (fluvial) Bidahochi Formation (blue arrow in Fig. 9), then renewed incision to form the modern LCR valley in the past 6 Ma.

Thermochronology samples from the Winslow and Holbrook transects near the cross section line are projected into a regional NE-SW cross section (30 times vertical exaggeration) in Figure 9. Stratigraphic thicknesses (from Fig. 5) are shown along the left axis. Note these are maximum thicknesses as Jurassic rocks may have been stripped prior to Late Cretaceous deposition. Inferred paleosurfaces through time at 10 Ma intervals were constructed by taking the weighted mean temperatures from the cooling paths (Figs. 6–8) and assuming a 25 °C surface temperature and a 25 °C/km geothermal gradient; these assumptions are inferred to produce minimal paleosurface elevations.

Thermochronology time-temperature paths suggest that cooling at the level of the currently exposed Triassic strata took place progressively from SW to NE and that most samples underwent a period of relatively rapid cooling of 40–50 °C but at different times. The inferred paleosurfaces in Figure 9 suggest that this cooling may have been due to northward cliff retreat of ∼1.7 km of upper Cretaceous strata (Gray Cliffs of the Grand Staircase) and that Cretaceous strata remained above samples B and C at the LCR river level until after 20 Ma. Cooling due to retreat of Jurassic and Triassic cliffs (>1 km) would have been another ∼35 °C in areas such as Lees Ferry where the Jurassic Glen Canyon and San Rafael groups are present but would have been less (perhaps <25 °C) in the Cameron transect where the Glen Canyon Group was partially present, and near zero to the Winslow and Holbrook transects and St. Johns area, where upper Cretaceous strata rested on Triassic rocks at 90 Ma (Dickinson, 2013).

Figure 10 shows our interpretation of the cliff retreat progression. For the Laramide time frame, all samples show a similar 90–80 Ma pre-Laramide peak temperature of >80–100 °C that corresponds to burial by ∼3–3.5 km of Mesozoic strata throughout the LCR valley and the Mogollon Rim. Similar thicknesses of upper Cretaceous strata may have extended SW across parts of the Arizona transition zone as shown by the K(80) line in Figure 10 that was drawn southwest of samples with apatite fission-track ages of 60–90 Ma in basement rocks (unlabeled red circles in Fig. 10). These samples suggest that 3–4 km of combined Phanerozoic cover and eroded basement were removed in these locations by 90–60 Ma. The highest elevation Moenkopi samples (J, K, L, M, and N) cooled rapidly between 70 and 50 Ma and are interpreted to record NE retreat of all of the Cretaceous plus any Jurassic section. We envision that 70–50 Ma paleorivers were bounded to the northeast by cliffs (K(50) line of Fig. 10) somewhat analogous to the Gray Cliffs section of the modern Grand Staircase. Samples H, I, F, and G remained at temperatures of 60–80 °C after 50 Ma. The resulting post-Laramide (ca. 50 Ma) topography based on thermochronology is compatible with outcrops of 65–55 Ma Music Mountain Formation gravels near sample N and 50–35 Ma Rim Gravels near samples J, H, I, K, and L.

In addition to the progressive cooling by NE cliff retreat, a pulse of denudation is recorded by samples A, B, C, and D near river level in the LCR valley; these remained >60 °C from 50 to 25 Ma, corresponding to preservation of 1–2 of Mesozoic strata that still covered these rocks until ca. 25 Ma. All samples in the LCR valley cooled to near-surface temperatures by 15 Ma, in agreement with geologic constraints for a 16 Ma Lake Bidahochi basal surface. A broad paleovalley shape at 15 Ma is suggested by subequal elevations of the reconstructed 15 Ma paleosurface above the Mogollon slope to the southwest, with highlands of the Chuska-Defiance uplift to the northeast (Fig. 9). Hence we infer that the NE rim of the LCR paleovalley at 25–15 Ma mimicked the modern Black Mesa topography.

Additional thermochronologic evidence for such a paleovalley is seen in Grand Canyon (Lee et al., 2013; Karlstrom et al., 2014) in rim and river rocks (samples 1–5 of Fig. 1). These samples are separated vertically by 1.5 km and resided at different post-Laramide temperatures compatible with a ∼20 °C/km geothermal gradient, followed by convergence of cooling paths 25–15 Ma attributed to carving of the East Kaibab paleocanyon. The north rim of the East Kaibab paleocanyon is inferred to have had cliffs of Jurassic and some Cretaceous rock at 15 Ma as shown by line K(15) in Figure 10.

⁴⁰Ar/³⁹Ar spectra for four dated basalts are shown in Figure 11, and analytical data are presented in section DR-3 (see ^{footnote 1}). All four of the dated samples yielded reliable eruption ages between 9.17 Ma and 0.339 Ma.

Red Butte basalt is an isolated basalt-capped butte near the south rim of Grand Canyon that overlies Triassic rocks (RB of Figs. 1 and 2). A new ⁴⁰Ar/³⁹Ar plateau age on this basalt is 9.17 ± 0.04 (Fig. 11), a refinement of the older K-Ar date of 9.73 ± 0.91 Ma (Billingsley, 2001). No gravels have been found beneath this basalt.

A basalt at Grassy Mountain in the Shivwits volcanic field (SV of Fig. 1) was reported by Lucchitta and Jeanne (2001) to overly Precambrian clast conglomerates on the west side of the mesa that were interpreted to be derived from south of Grand Canyon. However, no far-traveled gravels were observed by us, and instead, the basalt overlies locally derived gravels. ⁴⁰Ar/³⁹Ar dating yielded a 5.47 ± 0.05 plateau age.

A new ⁴⁰Ar/³⁹Ar date of 891 ± 13 ka (Fig. 11) was derived for the Black Point flow (location 5 of Fig. 1); this agrees with an independent ⁴⁰Ar/³⁹Ar date from Hanson (2007) of 873 ± 8 ka. This age corrects an earlier inaccurate K-Ar age of 2.43 Ma that was applied to the Black Point flow (Cooley et al., 1969; late Black Point surface) by several workers (e.g., Billingsley, 2001; Holm, 2001a). The Black Point flow originated in the San Francisco volcanic field and flowed into the LCR channel as evidenced by thick alluvial sediments encountered in a drill hole beneath the flow (section DR-1 [see ^{footnote 1}]; Haines and Bowles, 1976).

Tappan flow (location 7 of Fig. 1) was recognized by Colton (1937) to have flowed >50 km north down Tappan Wash from the San Francisco volcanic field and filled an earlier LCR river course (Rice, 1977). Our sample directly overlies LCR gravels. A new ⁴⁰Ar/³⁹Ar date on basalt from within a meter of the base of the flow is 339 ± 5 ka (Fig. 11).

The White Mesa alluvium was considered to be Early Miocene (Crooked Ridge paleoriver of Lucchitta et al., 2011, 2013b). Detrital-zircon (DZ) analysis challenged this assertion because a single zircon age of 15 ± 4 Ma indicated a much younger maximum depositional age (Fig. 12). ⁴⁰Ar/³⁹Ar dating of sanidine was applied to an ash bed that is interbedded with White Mesa alluvium on the Moenkopi Plateau (location 3 of Fig. 1). This age was 1.99 ± 0.002 Ma (Hereford et al., 2016). Similar ages on detrital sanidine (DS) of 2.02 ± 0.02 Ma at the base to 1.84 ± 0.05 Ma at the top of the section (Fig. 12) were found in the White Mesa alluvium of Crooked Ridge and White Mesa paleovalley, all at similar landscape positions. These dates indicate a ca. 2 Ma direct age for deposition of the ash and for the White Mesa paleoriver (Hereford et al., 2016). The full ⁴⁰Ar/³⁹Ar K-feldspar age distribution (Fig. 12A) ranges between 1.8 and 1600 Ma and is dominated by grains younger than 40 Ma (n = 105). These grains are most certainly sanidine. Detrital grains between 40 and 250 Ma are most likely sanidine, and grains older than 500 Ma represent K-feldspar (microcline and orthoclase) from Precambrian basement. The Precambrian K-feldspar ages reflect the time that the basement source cooled through ∼250–175 °C (Timmons et al., 2005). The Mesozoic grain ages suggest reworking of the Dakota Sandstone (Hereford et al., 2016).

In addition to the breakthrough about the young age for the White Mesa alluvium, the DS data provide new insights about provenance. The detrital-sanidine (DS) data for Cenozoic crystals show a 100-fold improvement in precision of DS single-crystal ages compared to those of DZ (Fig. 12B). This allows very high precision correlation of individual grains to specific ignimbrite eruptions. Many DS nodes are indistinguishable in age and K-Ca ratios from sanidine from San Juan volcanic field (SJVF) ignimbrites (Fig. 12C). For example, one population of White Mesa alluvium grains yields a weighted mean age of 28.209 ± 0.013 Ma and corresponding K/Ca values near 70; these values are distinctive of the 5000 km³ Fish Canyon Tuff sanidine, one of the world’s largest single eruptions, sourced from the La Garita caldera of the San Juan Mountains (Lipman and McIntosh, 2008). The next oldest population of DS grains is 28.685 ± 0.021 Ma with a K/Ca of ∼20 and correlates to Blue Mesa Tuff also from SJVF. None of the dated grains yielded ages between 28.3 and 28.6 Ma, a major pulse of volcanism in New Mexico’s Mogollon Datil volcanic field (McIntosh et al., 1992). Correlating prominent age peaks as well as temporal gaps in the White Mesa DS data indicate that the ultimate primary source for Oligocene detritus in the White Mesa alluvium was the San Juan Mountains of Colorado rather than volcanic fields in New Mexico.

Many of the White Mesa DS age peaks are currently uncorrelated. For instance, the young grains near 2 Ma were first thought to correlate to the large Huckleberry Ridge eruption of Yellowstone. However, upon closer comparison to manuscripts by Ellis et al. (2012) and Singer et al. (2014), the DS peak at 2.0 ± 0.02 Ma (normalized to the same standard ages) is distinguishably younger than the published 2.15–2.10 Ma Huckleberry Ridge ages. The DS data are chemically distinct as well: K/Ca for this DS peak is ∼100, whereas Huckleberry Ridge is ∼10. Similarly, the 10–5 Ma grains are currently uncorrelated ages but would be compatible with northern Great Basin sources (Perkins et al., 1998; Anders et al., 2009; Best et al., 2013).

Precise DS data near 18.8 Ma are shown in Figure 12D. Hereford et al. (2016) suggested that these grains could have originated either from the Peach Springs tuff near Kingman, Arizona, or the Apache Leap tuff near Phoenix, Arizona. Hereford et al. (2016) could not distinguish a Peach Springs versus Apache Leap source because the published ⁴⁰Ar/³⁹Ar geochronology of these tuffs was too imprecise (i.e., existing geochronology data; McIntosh and Ferguson, 1998; Ferguson et al., 2013) could not decipher an age difference between the tuffs). Also, the scatter in the published K/Ca data prevented chemical diagnosis (Fig. 12D). We have re-dated both tuffs with the ultra-precise ARGUS VI mass spectrometer and find that, even at ca. 5 ka precision, we cannot determine an age difference between these tuffs (Fig. 12D). However, the K/Ca data for White Mesa DS grains (K/Ca = 45–55) are a better match to the higher values from the Apache Leap tuff (K/Ca = 30–55) indicating Apache Leap tuff was a source for White Mesa sanidine. This might support the conclusions of Potochnik (2011) for northerly flow of a paleo–Salt River into the proto–Little Colorado River basin. However, wind transport of sanidine detritus is likely and would preclude firm conclusions about fluvial connections between White Mesa alluvium and the areas south of the Mogollon Rim. Thus, the White Mesa alluvium contains detritus ultimately derived from both the San Juan volcanic field and the Superstitions Mountain ignimbrite eruptions.

Reversal of the NE- and E-flowing post-Laramide rivers to SW-flowing streams along the Mogollon Rim in the LCR headwater region took place after 18.6 (and perhaps after 14.5 Ma) for the Salt River as a result of Basin and Range extension and collapse of the Mogollon highlands (Faulds; 1986; Potochnik and Faulds, 1998; Potochnik, 2001). This drainage reversal may have led to the 16–6 Ma period of internal drainage in the LCR valley (Bidahochi Formation). From 16 to 6 Ma, the LCR paleovalley was a shallow Bidahochi depositional basin (Fig. 1; McIntosh and Cather, 1994; Scarborough, 2001) with a relatively stable local base level. The lower Bidahochi lacustrine facies was deposited from 15.84 to 13.71 Ma as constrained by interbedded ash layers; it reached a total stratigraphic thickness of ∼100 m (Dallegge et al., 2001). The depositional environment is interpreted to be a shallow playa or marsh, with no evidence that Hopi paleolake was ever a deep lake (Dickinson, 2013). Figure 1 shows structure contours of the base of the Bidahochi Formation (based on outcrops and drill-hole elevations) that outline the depocenter (Fig. 1; Cooley and Akers, 1961).

Using the landscape position of dated ash beds as long-term incision and denudation constraints, the post–16 Ma long-term average bedrock incision rate for the LCR valley has been very slow (17–31 m/Ma for points 9A and 9B of Table 1). These slow, post–16 Ma averaged rates indicate that 16–6 Ma aggradation has been only slightly exceeded by post–6 Ma incision. From 16 to 6 Ma, aggradation of sediments in the central reaches filled paleovalleys with fine-grained lacustrine deposits in the lower Bidahochi Formation and with coarser fluvial upper Bidahochi Formation (Dallegge et al., 2001; Dickinson, 2013). In the headwater region, gravels were deposited until ca. 6–4 Ma (McIntosh and Cather, 1994), suggesting a generally aggradational setting from 16 to 6 Ma.

Figure 13 shows a series of erosional surfaces that record denudation of the LCR valley generalized from mapping of Cooley et al. (1969), with new geochronologic control from this paper. Surfaces high in the landscape in the Chuska Mountains and on Black Mesa were named the Valencia surface and interpreted to be the result of 26–16 Ma (post–Chuska Sandstone) denudation (Cooley et al., 1969). Cooley et al. (1969) also “correlated” the sub-basalt “surface” at Red Butte on the south rim of Grand Canyon as Valencia surface (Cooley and Akers, 1961). In addition to the Red Butte basalt, other locations in the region show 10–6 Ma basalts overlying Triassic strata (Fig. 3): Flagstaff area (Holm, 2001a), Red Butte, Anderson, and McMillan mesas on the Mogollon slope (Holm, 2001b), and a 6.41 Ma basalt flow at Chevelon Butte (Cather et al., 2008). At Long Point in Arizona, a 6.7 Ma basalt flow overlies the Music Mountain Formation at ∼1800 m elevation suggesting no significant degradation there between 65 and 7 Ma. Similarly, 7–5 Ma basalts of the northern Shivwits Plateau on the North Rim of Grand Canyon overlie Triassic rocks (Lucchitta and Jeanne, 2001). Although not a single “surface or peneplain,” the regional distribution of 10–6 Ma basalts overlying Triassic strata indicate a relatively widespread, locally aggrading, moderate relief landscape carved into lower Triassic strata. By 15–6 Ma, Jurassic (Vermillion Cliffs–type) and Cretaceous (Gray Cliffs–type) escarpments had retreated north of Grand Canyon and east of Bidahochi Formation outcrops.

This section presents new incision control points for the LCR in the context of a synthesis of all published incision constraints (Table 1 and section DR-1 [see ^{footnote 1}]). Figure 14A shows the longitudinal profile of the LCR and some important tributaries, the bedrock that underlies the LCR profile, and the rim of the deeply incised LCR canyon as it enters Grand Canyon at the left edge of the diagram. Red dots represent incision constraints projected into the profile keyed to Table 1.

Modern surface water of the LCR is from several sources. In the lower reach, Blue Springs is a spring source of perennial flow (∼220 cubic feet per second; Crossey et al., 2009) for the LCR that comes from the Redwall-Muav aquifer. In the central reach, flow is ephemeral and controlled by spring runoff and flash floods. In the upper reaches, surface flow is from spring sources and snowmelt from the White Mountains.

At the CR-LCR confluence, estimating total magnitude of post–6 Ma incision depends on the depth of the postulated 25–15 Ma East Kaibab paleocanyon, which would have formed the initial base level for integration of the CR with the paleo–LCR valley 5–6 Ma. Scarborough (2001) used an elevation of 1830 m and envisioned that spillover of Hopi Lake at Cape Solitude initiated downward CR integration. However, modeling of thermochronologic data from near the confluence suggests that river-level rocks had cooled to temperatures of 20–30 °C by 10 Ma and that East Kaibab paleocanyon may have been cut to below the level of the Kaibab limestone in Eastern Grand Canyon (Lee et al., 2013). Dickinson (2013) rejected the spillover model based on conclusions that Hopi Lake was not at the correct elevation and was never a deep lake, but Crossey et al. (2015) cited groundwater sapping as a possible additional mechanism for downward integration. A speculative range of minimum elevations (maximum paleocanyon depths) for the ca. 6 Ma confluence is 1200–1600 m (Fig. 3). Lee et al. (2013, their fig. 9E) showed a preferred elevation of ∼1600 m (Fig. DR-1.1 [see ^{footnote 1}]), which we use for paleoprofile reconstruction in Figure 14A (point 1A of Table 1).

Positions of the 4–2 Ma paleoprofiles in eastern Grand Canyon are not constrained unless one assumes semi-steady incision (e.g., Karlstrom et al., 2008). Polyak et al. (2008) reported a 2.68 Ma age from a cave speleothem in Kwagunt Creek that is ∼5 river miles upstream from the confluence. The sample is at an elevation of ∼1290 m, 260 m above the Kwagunt Creek tributary and 445 m above the river. They assumed the water table was flat and that the speleothem formed during downward passage of the water table that tracked river level as CR-incised Grand Canyon. If so, this date would indicate the CR level was at 1290 m at 2.68 Ma (point 1B, Fig. 14A), and the 2 Ma paleoprofile was at ∼1200 m. However, this cave is in a tributary; so the incision rate and magnitude are maxima (Table 1; Karlstrom et al., 2008; Crow et al., 2014). Nevertheless, an elevation of ∼1200 m for the 2 Ma confluence is close to the estimate of ∼1145 m obtained by assuming that the measured Quaternary rates of 160 m/Ma (documented over the past 623 ka; Crow et al., 2014) may have been semi-steady back to 2 Ma (Fig. 3).

About 4 km downstream from the confluence along the mainstem of the Colorado River, a bedrock strath 70 m above the river was dated at 623 ka by Crow et al. (2014) yielding a bedrock incision rate of 147 m/Ma (assuming a depth to bedrock of ∼21 m). This rate was refined by strath-to-strath dating from terrace flights between river miles 57 and 65, both above and below the LCR confluence, that indicate a semi-steady average bedrock incision rate of 160 m/Ma over the past 623 ka (Crow et al., 2014) placing the CR-LCR confluence at 985 m elevation at 1 Ma (Fig. 3) and 925 m at 0.623 Ma (point 1C of Fig. 14A and Table 1).

In the central reaches of the LCR, above the LCR knickpoint near Cameron, long-term post–6 Ma bedrock incision rates can be estimated based on the position of the 8–6 Ma eruptive surface for the Hopi Buttes basalts. This was a long-lived, low-relief surface (Hopi Buttes surface of Cooley et al., 1969). Many of the 8–6 Ma eruptions formed maars where lavas interacted with groundwater in what was likely the bottom of the LCR valley. This eruptive surface has been variably dissected; it is essentially intact in the northeastern parts of the Hopi Buttes (elevation ∼1900 m), where it is marked by thin, far-traveled flows that interfinger with sediments of the upper (fluvial) Bidahochi Formation. Farther southwest, toward the modern LCR, the volcanic field offers exposures of maar volcanic features that developed as much as 300 m below the eruptive surface (White, 1991; Dallegge et al., 2001). Calculated rates for individual dated basalts (Table 1) are interpreted as maximum denudation rates relative to the modern LCR as they are at large distances (18–75 km) from the modern river and lack paleo-LCR gravels. However, assuming a laterally extensive low-relief surface, LCR valley denudation rates range from 31 to 81 m/Ma (Table 1). Coliseum maar is within 18 km of the LCR and gives a denudation rate of 50 m/Ma, which is similar to the rate of 65 m/Ma calculated using the average age (7.2 Ma) and height (470 m) of the reconstructed eruptive surface (Table 1; #10).

On the Mogollon slope southwest of the LCR (Fig. 1), Holm (2001a) used basalt-capped surfaces in tributaries to postulate average incision rates of ∼100 m/Ma (33–120 m/Ma) over the past 6 Ma. The best individual rates are from basalts closest to the LCR. These give average rates of 51 m/Ma since 6.41 Ma (Chevelon Butte; #14), 107–115 m/Ma over the past 1.92–1.63 Ma (East Sunset Mountain; #13), and 148 m/Ma over the past 0.81 Ma (Woodhouse Mesa; #6). These data suggest increasing incision rates through time in the past 6 Ma (Fig. 14B).

Along Moenkopi Wash, a tributary north of the LCR, new dating of ashes can be used to calculate tributary incision rates. The 1.993 ± 0.002 Ma ash in a deposit of White Mesa alluvium on top of Moenkopi Plateau (#3 of Table 1; Hereford et al., 2016) yields a tributary incision rate of 120–155 m/Ma when the planar tread of that deposit is projected along its slope to Moenkopi Wash (Fig. DR-1.2 [see ^{footnote 1}]). In Blue Canyon at the edge of Black Mesa (#4), an ash correlated to the ca. 1 Ma Bishop Tuff was found in the L2B deposit. Projecting the tread of that deposit to Moenkopi Wash yields a similar incision rate of 120–125 m/Ma. The similar rate of ∼120 m/Ma over the past 2 Ma for tributaries to the north and south of the LCR implies this rate for the mainstem LCR and was used to reconstruct the ca. 2 Ma paleoprofile in the central reaches of Figure 14A.

The new ⁴⁰Ar/³⁹Ar date of 891 ± 13 ka for the Black Point flow (#5) agrees with an independent ⁴⁰Ar/³⁹Ar date from Hanson (2007) of 873 ± 8 ka. The Black Point flow originated in the San Francisco volcanic field and flowed into the LCR channel (Fig. DR-1.3 [see ^{footnote 1}]). Rice (1980) reported LCR gravels atop the flow, and drill data indicate it overlies 22 m of LCR gravel. Gravels overlie a bedrock strath 151 m above the modern river (Fig. DR-1.4 [see ^{footnote 1}]; Haines and Bowles, 1976) giving a bedrock incision rate of 193 m/Ma, somewhat higher than the 160 m/Ma post–623 ka rate at the CR-LCR confluence. The Black Point flow age constrains the 122–152 m terraces (L-2B of Cooley et al., 1969) to be ca. 1 Ma. The Bishop tuff at Blue Canyon is in gravel of the L-2B surface, which provides another direct date of L-2B. This surface extends downstream as far as Tuba City and is projected to Moenkopi Wash in Figure DR-1.2 (see ^{footnote 1}) to give heights of 118–126 m above Moenkopi Wash and corresponding tributary incision rates of 118–126 m/Ma.

Tappan flow (#7) was mapped in detail by Rice (1977). It flowed >50 km down Tappan Wash, then thickened dramatically (from 15 to 30 m), filled the steep-walled LCR paleochannel, and flowed both upstream (1.6 km) and downstream (∼11 km) within the earlier LCR river course (Fig. DR-1.3 [see ^{footnote 1}]). The new ⁴⁰Ar/³⁹Ar date on basalt from within a meter of the base of the flow is 339 ± 5 ka (Fig. 11D). At this location, Tappan flow overlies LCR gravels that sit on a bedrock strath 57 m above the modern LCR channel giving an incision rate of 167 m/Ma. This sample is from just above the LCR knickpoint and may constrain the 46–61 m terrace (L-3B of Cooley et al., 1969) to be ca. 350 ka and, together with the Black Point incision constraint, suggests that this region, like the CR, has undergone semi-steady average bedrock incision of ∼170 m/Ma over nearly the past one million years.

Shadow Mountain basalt flows give an ⁴⁰Ar/³⁹Ar isochron date of 280 ± 50 ka (Conway et al., 1997). The elevation of the lowest flow base is 50 m above the level of Moenkopi Wash (Fig. DR-1.2 [see ^{footnote 1}]), which gives a maximum tributary incision rate of 179 m/Ma, similar to the mainstem rate of 167 m/Ma from the Tappan flow over the same ∼300 ka time interval. Rice (1980) mapped basalt-bearing gravels and tephra that were derived from Shadow Mountain in Hopi Trail Wash within 2 km of the LCR confluence; these project 70 m above the LCR and would give an incision rate of 250 m/Ma using the newer ca. 0.28 Ma ⁴⁰Ar-³⁹Ar age of Conway et al. (1997), whereas Rice (1980) reported 113 m/Ma using an older 0.62 Ma K-Ar date on Shadow Mountain (Damon et al., 1974; see section DR-1 [see ^{footnote 1}]).

Figure 15 shows a geologic map of the headwater region of the LCR at the eastern edge of the Springerville volcanic field. Basalt flows from volcanic fields have “run-out” geometries suggesting they followed paleodrainages, and they now stand as topographically inverted elongate mesas. Results presented here are preliminary as only a few of these mesas have gravels mapped beneath the flows, and ages are K-Ar ages that are rarely from the optimum locations to calculate incision rates. Figure 16 shows the base of basalts projected into the LCR profile, and red stars show incision rate control points listed in Table 1. A 6.8 Ma basalt caps thin beds of Fence Lake and/or Bidahochi gravels at Flat Top (#16), and a nearby 6.03 Ma basalt overlies river gravels; both indicate long-term average incision rates of 40–45 m/Ma for the headwater region (Embid, 2009). However, another 6.52 Ma flow is at a considerable lower elevation relative to the LCR (#17) and gives a lower incision rate of 31–38 m/Ma relative to the LCR. Incision rates of ∼23 m/Ma are obtained from ca. 6 Ma basalts in eastern tributaries to the LCR. Figure 14B shows considerable post–6 Ma incision rate variation in the headwaters with rates of 10–50 m/Ma.

Post–6 Ma flows closest to the modern LCR show a general relationship where younger flows are inset into older flows (Fig. 16). Based on its paleochannel-like geometry, we infer that the 3.67 ± 0.12 Ma Coyote Wash flow marks the paleochannel of the LCR, although we have not mapped gravels at the base of the flow. The lobe shows an uneven height above the modern river profile, perhaps compatible with a graben between the Coyote Wash fault and inferred faults (Embid, 2009). The height of the base of the flow at Lyman Lake is 1950 m, giving an incision rate of 37 m/Ma (Fig. 16A), similar to average post–6 Ma rates.

A prominent but higher flow lobe is called Black Mesa; it has the 2.46 ± 0.04–2.37 ± 0.03 Ma Mesa Parada flow possibly inset on its northeast edge (Fig. 15). The flow(s) appear to have originated in the Red Hill–Quemado area to the east and filled a paleovalley cut into older basalts (McIntosh and Cather, 1994). Black Mesa flows arc around St. Johns Dome and fill a channel cut in Bidahochi gravels from 2290 to 1980 m elevation at the eastern edge of the study area. The higher projected height of this flow relative to the younger Coyote Wash flow (Fig. 16) is more compatible with using the 6.03 Ma date to the SE (#19 and #37) rather than the 2.37 Ma Mesa Parada date to estimate incision rates. The NW tip of Black Mesa is 250 m above the LCR and would yield a provisional LCR incision rate of 41 m/Ma using the 6.03 Ma date, whereas the local incision rate from sample #37 to Coyote Creek is 35 m/Ma.

A regression of the basalt data along the LCR from south of Springerville to St. Johns (inset to Fig. 16) suggests a semi-steady incision rate of 41 m/Ma from 6 to 1 Ma. Outlier points may reflect complex variation expected in incision rates across paleo-knickpoints (#21) and uncertain dating of flows (#22). A similar incision value of 41 m/Ma is obtained from Largo Creek–Carrizo Wash tributary to the LCR based on a 2.46 Ma basalt flow 100 m above the creek (Love and Connell, 2005). Basalt data from eastern tributaries of the LCR (#35–44 of Fig. 1 and Table 1; McIntosh and Cather, 1994) show somewhat slower but still semi-steady incision rates of 23 m/Ma over the past 6 Ma. Thus, given available dating and projection of basalt surfaces, we interpret the headwater reaches and its tributaries of the LCR to record semi-steady bedrock incision of ∼20–40 m/Ma from 6 to 1 Ma.

Figure 15 shows gravel terraces that were mapped by Sirrine (1958), and Figure 17A projects gravel straths into the LCR profile (Embid, 2009). The data are preliminary and would benefit from additional terrace mapping, correlation, and dating of gravels, for example with detrital-sanidine methods. Nevertheless, younger terraces are shown to be inset into older basalts and older terraces supporting the general notion of progressive incision since 6 Ma. Remnants of the oldest gravel (Tg1) are inset into Eagar Formation along Coyote Creek at elevations of 2165–2195 m. Up to 15 m thicknesses of Qg2 gravel underlie much of the Springerville volcanic field at elevations 60–100 m above the LCR riverbed. Based on their heights and comparison to dated basalts, Tg1 is likely >2 Ma, and Qg2 is ca. 1–2 Ma. Projections of Tg1 and Qg2 suggest convexities that are similar in shape to the modern knickpoints. Qg3 gravel (3–5 m thick) is inset ∼10 m below Qg2 and is preserved in Big Hollow Wash beneath 1.68–1.87 Ma basalts. Undated travertines along Carrizo Creek are deposited on this surface.

Uranium-series dates on travertine platforms, mound spring complexes, and on travertine-cemented gravels in the LCR headwater region near Lyman Lake also provide incision constraints. Figure 17B shows that platforms and mounds extend to variable elevations in the landscape due to artesian spring accumulations (Embid, 2009). However, where travertines cement gravels in strath terraces, travertine ages show semi-steady bedrock incision rates of 38 m/Ma from 300 to 100 ka and increased rates since then (Fig. 17C). The oldest dated travertine sample that overlies a strath is 304 ± 13 ka (#27) and was sampled from the southern end of the Salado travertine platform at an elevation of 1827 m; it overlies 1 m of gravels that in turn rest on a strath 26 m above the LCR giving an average incision rate of 87 m/Ma. A sample (#30) 5.5 m higher in the travertine platform gives an age of 272 ka and an incision estimate (from the same strath to the river) of 109 m/Ma (Fig. DR-1.6 [see ^{footnote 1}]). A nearby mound next to Lyman Lake dam (#28) of 281 ± 15 ka formed 60 m higher atop river gravels on a well-exposed strath; this sample gives an incision rate of 153 m/Ma. Sample #29 is from the northern end of Lyman Lake, and it gives an age of 255 ± 6 ka and an incision rate of 94 m/Ma (Fig. DR-1.5 [see ^{footnote 1}]). Sample # 31 gives an age of 224 ± 3 ka with a strath at 30 m and an incision rate of 134 m/Ma (Fig. DR-1.7 [see ^{footnote 1}]). Samples of ca. 100 ka travertines, #32–34, are all from the bases of individual mounds between Lyman Lake and Salado Springs overlying river gravel and/or strath contacts and are dated 101 ± 1, 97 ± 1, and 78 ± 2 ka, respectively. They show that the average river incision rate in the past 100 ka increased to ∼278 m/Ma (average of three data points).

Figure 14 provides a summary of the differential incision history in both space and time for the LCR. Basalt data from the upper reaches of the LCR and its tributaries suggest semi-steady incision at rates of 30–40 m/Ma from ca. 6 Ma to 300–100 ka. This rate is similar to slow denudation rates well away from the LCR at Anderson Mesa. Central reaches of the LCR give higher average rates from 6 Ma to present. For example, the incision magnitude since the formation of the Hopi Butte eruptive surface at ca. 7.2 Ma is ∼470 m, yielding a long-term average incision rate of 65 m/Ma. However, our new data suggest this long-term rate is best interpreted in terms of slow 30–40 m/Ma rates from 6 to 2 Ma (similar to upper reaches) followed by accelerated incision in the past 2 Ma. Thus, inferred paleoprofiles from 6 to 2 Ma in Figure 14A show relatively slow incision between 6 and 2 Ma.

The new ⁴⁰Ar/³⁹Ar sanidine age of ca. 2 Ma for White Mesa alluvium both on the Moenkopi Plateau (Blue Point tuff) and in the White Mesa paleovalley gives incision rates of 137 m/Ma (range of 119–155 m/Ma) for Moenkopi Wash over the past 2 Ma. This wash has its modern confluence with the LCR just above the LCR knickpoint, and we infer that its paleoconfluence at 2 Ma was also above the knickpoint (Fig. 14A). Tributaries typically incise at the same rate as the mainstem at their confluence, and we note that there were similar average incision rates for tributaries north (Moenkopi Wash at 137 m/Ma over 2 Ma) and south (East Sunset Mountain relative to East Clear Creek at 115–120 m/Ma over 1.92–1.63 Ma) of the LCR. We interpret this to indicate that the average incision rate of the central reaches of the mainstem LCR over the past 2 Ma was also 115–140 m/Ma and hence ∼3 times faster than the 6–2 Ma rates of 30–40 m/Ma in the headwaters.

Post–1 Ma rates were perhaps still higher: 122 m/Ma for the Blue Canyon tuff (#4), 170 m/Ma in the mainstem at Black Point basalt (#5), and 148 m/Ma at Woodhouse Mesa (#6). Youngest incision points give rates of 160 m/Ma over 625 ka in the mainstem CR (#1C) and 167–179 m/Ma over 300 ka in the LCR and Moenkopi Wash tributary (#7A and #7B). The overall data suggest an increase in incision rates from 40 to 170 m/Ma in the lower and central LCR over the past 2 Ma (upper curve in Fig. 14B).

Tributaries are required to keep up with the incision rate of the mainstem (or form a hanging valley) such that incision rates at the confluence of the LCR and CR are assumed to have been the same, and hence incision rates in eastern Grand Canyon may also have accelerated in the past ∼2 Ma. This section discusses possible driving forces that may help explain the post–2 Ma increase in incision rates in the lower LCR and perhaps in eastern Grand Canyon. This discussion includes a summary of recent models for the Grand Canyon region that variably emphasize tectonic, geomorphic, and/or climatic drivers for landscape evolution. Multiple forcings are probable, but their relative importance remains in debate; thus this section also includes ideas about future work to test alternate hypotheses.

Epeirogenic uplift of the southwestern Colorado Plateau relative to the Basin and Range (Karlstrom et al., 2008) and of the Rocky Mountains relative to the Colorado Plateau (Karlstrom et al., 2011) in the past 10–6 Ma may have been driven by changes in mantle buoyancy. Moucha et al. (2008, 2009) invoked large-scale mantle convection. Van Wijk et al. (2010) invoked edge-driven smaller-scale mantle convection. Levander et al. (2011) invoked lithospheric delamination and asthenospheric return flow. Crow et al. (2011, 2014) invoked a migrating zone of mantle modification recorded at Earth’s surface by an inboard sweep of basaltic volcanism. However, the timing, scale, and mechanisms, as well as the existence (Pederson et al., 2013), of post–10 Ma mantle-driven surface uplift remain controversial.

Most mantle uplift models propose regional-scale (western United States) uplift. The Crow et al., (2014) model is of most appropriate time frame, location, and spatial scale to explain the acceleration of incision in the LCR region. As shown in Figure 18, progressive NE-migration of a sharp velocity contrast along a step in the lithosphere-asthenosphere boundary that is now located near Lees Ferry may be manifested by a trailing sweep of basaltic volcanism that is progressively becoming younger and more asthenospheric toward the Plateau center (Crow et al., 2011). In this model, the sweep of volcanism itself can uplift Earth’s surface due to constructional topography (e.g., San Francisco Peaks reach 3.8 km in elevation), and the asthenospheric upwelling that is driving the volcanism may add buoyancy that adds several hundred meters of long-wavelength doming and differential surface uplift. The latter was posited to be driving the observed higher incision rates in the past 0.6 Ma in eastern Grand Canyon relative to central and western Grand Canyon (Crow et al., 2014). Figure 18 expands this concept, adds the lower LCR region of higher river incision (dashed oval), and suggests that the Mesa Butte fault system, which the San Francisco Peaks are in part built along, may be a zone that is focusing NE propagation of large-volume volcanism toward the large-velocity contrast centered under Lake Powell. Increased subsurface magmatic flux and mantle buoyancy may have started passing beneath the CR-LCR confluence region at ca. 2 Ma. This is consistent with NE migration of volcanism in the San Francisco volcanic field at a rate of 30 km/Ma in the past 2 Ma (Tanaka et al., 1986). Surface doming would need to have been on the order of 100–200 m over the past 2 Ma to explain the four-fold difference in incision magnitude (and rates) of 60–80 m (30–40 m/Ma) in the LCR headwaters versus 280–320 (140–160 m/Ma) in the lower reaches. We assume that river incision would have kept pace with uplift (Braun, 2010), likely by steepening of the profile, and that isostatic response (rebound) to the ∼220–240 m of differential denudation since 10 Ma might account for nearly half (∼100 m) of the differential uplift (Lazear et al., 2013).

Headwater uplift and base-level fall are difficult to distinguish geomorphically, and both cause changes in river gradients. The profile of the LCR shows a major knickpoint in its lower reaches that suggests base-level fall, but headwater plateau-like uplift may have caused a downstream knickpoint and higher incision rates to form at the edge of the uplifted area, similar to other examples (Whipple, 2004; Rosenberg et al., 2014).

Several studies have suggested that rivers are responding to surface uplift associated with a NE-trending zone of Quaternary magmatism of the Jemez lineament that extends across New Mexico from the Springerville volcanic field to the Great Plains (Wisniewski and Pazzaglia, 2002; Nereson et al., 2013; Channer et al., 2015). Rivers that cross suborthogonally across the Jemez lineament have convex profiles, convex terrace flights, and tilted basalt paleosurfaces that are interpreted by these studies to reflect several hundred meters of broad epeirogenic doming above the lineament. The Springerville volcanic field has been volcanically active for the past 10–6 Ma, and the Hopi Buttes field was active 8–6 Ma such that the turnaround from aggradation to incision after 6 Ma and the semi-steady 40 m/Ma incision rates since then might reflect magmatically driven uplift. Springerville has also had voluminous basaltic volcanism in the past 2 Ma; but present data suggest that the Springerville headwater reaches have had semi-steady incision from 6 to 0.3 Ma with no marked increase in incision at 2 Ma, whereas this acceleration seems to be localized in lower reaches.

Geomorphic explanations for differential incision invoke transient knickpoints (Darling et al., 2012; Donahue et al., 2013; Aslan et al., 2014). Several studies have proposed passage of a transient wave of incision in the mainstem CR that now forms the regional-scale knickpoint at Lees Ferry (Pelletier, 2010) or perhaps has migrated to Cataract Canyon (Fig. 2; Cook et al., 2009). Cook et al. (2009) noted that the LCR and other CR tributary channel profiles are downturned in their lower reaches with knickpoints 150–200 m above their confluence, consistent with a transient that has passed the CR-LCR confluence in eastern Grand Canyon. Their modeling explored the incision rate changes expected at a given location due to passage of a knickpoint, and they modeled a scenario in which the Colorado River in the Glen Canyon region underwent a pulse of rapid incision in the past ∼500 ka. Darling et al. (2012) suggested incision in the Glen Canyon region has accelerated in the past 300 ka. Abbott et al. (2015) invoked the Cook et al. (2009) model and suggested passage of an incision wave through central Grand Canyon between 500 and 400 ka. However, such recent timing is not compatible with Grand Canyon incision data for semi-steady incision (Crow et al., 2015). The Cook et al. (2009) modeling (their fig. 12) suggests that incision rates would peak and then decline dramatically at any given point following passage of the transient. However, existing data in the CR for semi-steady incision over 625 ka in eastern Grand Canyon, and ∼4 Ma in western Grand Canyon, and a similar steady incision rate of 160–170 m/Ma in the mainstem of the lower LCR based on the 0.34 Ma Tappan and 0.89 Ma Black Point basalt flows do not support the modeled changes in incision rates expected from passage of a transient knickpoint in the CR or LCR within the past ∼1 Ma. Instead, incision data favor a longer-term, steady-state process such as tectonic uplift.

The Lees Ferry and LCR knickpoints are both located near the Kaibab Limestone– Moenkopi Formation contact. The steeper parts of the river profiles are in more resistant carbonate-dominated strata of the Paleozoic section below the knickpoint and shallower gradients are observed in more erodible sandstone-shales of the Triassic section above (Figs. 2 and 14). Cook et al. (2009) proposed that the arrival of an incision transient at the Lees Ferry upstream-dipping lithologic knickpoint may have initiated a pulse of incision to sweep rapidly (400 m/Ma) through the upper weaker rocks until it encountered stronger rocks far upstream in Cataract Canyon. Similar models for rapid knickpoint propagation through weak units were proposed for the Mancos Shale by Darling et al. (2012) and Aslan et al. (2010, 2014). Donahue et al. (2013) posited that highest incision rates within the harder basement rocks of the knickpoint region of the Gunnison River relative to slower rates in upstream and downstream reaches underlain by erodible Mancos Shale were evidence that the Gunnison knickpoint is a young transient feature of the river system.

Darling and Whipple (2015) expanded the discussion of knickpoints, tributary profile patterns, and bedrock strength to explain profile patterns that form when tributary drainages have easily eroded strata overlying stronger strata. They noted that knickpoints form due to more rapid erosion within the canyon than on the surrounding lithologically controlled, low-relief benches. But, similar to their conclusion for the Hualapai Plateau of westernmost Grand Canyon, the region above the LCR knickpoint is an older and long-lived stripped surface (e.g., Childs, 1948), not a lithologically controlled bench. Importantly for understanding the LCR knickpoint, there are similar incision rates above, within, and below the knickpoint: Incision through Mesozoic strata at Black Point above the knickpoint has averaged 170 m/Ma over 890 ka; incision of the Tappan flow at the LCR knickpoint has averaged 167 m/Ma over 342 ka; and incision in the mainstem CR near the confluence region has been semi-steady over the past 625 ka at 160 m/Ma. Thus, we infer that the LCR knickpoint is a semi-stable feature of the LCR profile that has been pinned (or hung up) at the Paleozoic–Mesozoic contact for most of the past 2–1 Ma (Fig. 14A) and reflects the increased gradient needed to carve Paleozoic carbonates at about the same rate as Mesozoic strata. Additional strath-to-strath incision data in multiple locations are needed to test if the same explanation may be true for the Lees Ferry knickpoint. Our present view is that the prominent knickpoints at Lees Ferry and the LCR are bedrock-influenced knickpoints that resulted initially from headwater uplift and/or base-level fall that accompanied 5–6 Ma downward integration of the CR and LCR through Grand Canyon, with amplification of the knickpoint caused by renewed 2 Ma differential uplift.

There are several postulated changes in post-Miocene climate that may have influenced landscape evolution. A mid-Miocene warm period (17–15 Ma—Zachos et al., 2004; or 17–12 Ma—Chapin, 2008) was followed by cooling temperatures since then (Fox and Koch, 2003). Chapin (2008) argued that opening of the Gulf of California by ca. 6.4 Ma resulted in an intensification of the North American Monsoon and increased erosion in the Southwest. This proposed climate shift may have combined with a global δ¹⁸O ∼2‰ negative seawater shift starting ca. 5.3 Ma (Zachos et al., 2004) that may reflect increased erosion. A mid-Pliocene period (ca. 3 Ma) that was about ∼3 °C warmer and more humid than today (Thompson and Fleming, 1996; Haywood and Valdes, 2004) was followed by cooling temperatures since ca. 2.7 Ma. This has been proposed to have led to an intensification of glacial climates (Haug et al., 2005) accompanied by increased seasonality of climate and higher erosivity of geomorphic systems, perhaps globally (Molnar, 2004; Molnar et al., 2006). At shorter time scales, oscillating Quaternary glacial-interglacial conditions have caused cyclic changes in bedrock incision versus aggradation recorded by river terraces (cf. Bull, 1991; Hancock and Anderson, 2002; Pederson et al., 2006). For time periods of <150 ka (e.g., Karlstrom et al., 2013b; Pederson et al., 2013), the terrace record may not have averaged out glacial-interglacial oscillations and hence may give maximum bedrock incision rates.

The effects of climate changes on incision rates and patterns in the Rocky Mountain– Colorado Plateau region are not well documented but would be expected to have regional effects rather than the observed sub-regional incision variations (e.g., Nereson et al., 2013; Rosenberg et al., 2014). The 6–5 Ma climate changes seem likely to have influenced the change from aggradation to erosion ca. 6–5 Ma in the Southwest and the integration of the Colorado River system. Whether climate change was the main driver or interacted with tectonic influences requires additional data on spatial variations and temporal incision variations and more precise dating of terraces and paleoriver deposits of this age. Climate change to more erosive climate at ca. 2.6 Ma is a possible contributor to the acceleration of incision in the past 2 Ma in the LCR and eastern CR regions but is hard to reconcile with continued steady incision in western Grand Canyon and the LCR headwaters since ca. 6–4 Ma and cannot by itself explain the observed differential incision.

The incision and denudation history of the Little Colorado River region over the past 70 Ma was punctuated by three denudation pulses, each linked to river integration and/or canyon- carving episodes and coupled to pulses of cliff retreat. Post-Laramide denudation at 70–50 Ma stripped Paleozoic and Mesozoic strata NE off the Mogollon highlands via cliff retreat. Mid-Tertiary canyon carving of the LCR paleovalley and East Kaibab paleocanyon took place 25–15 Ma and caused the regional post–Chuska sandstone denudation of the southern Colorado Plateau. Post–6–5 Ma denudation was initiated by headwater uplift and/or base-level fall and resulting downward integration of the CR to the Gulf of California. Six to 5 Ma CR integration across the Vermilion Cliffs at Lees Ferry allowed drainage from the central Colorado Plateau to flow south along the east side of the Kaibab monocline to a confluence with the paleo-LCR and its East Kaibab paleocanyon (speculatively at ∼1600 m elevation) near the modern CR-LCR confluence.

Long-term post–6 Ma LCR average bedrock incision measured using dated basalts shows that incision magnitude and average rate decrease upstream on the LCR: ∼800 m (>130 m/Ma) at the confluence, ∼450 m (∼75 m/Ma) in the Hopi Buttes area, ∼240 m (∼40 m/Ma in the LCR headwaters), and ∼100 m 150–200 m (23 m/Ma) in the eastern LCR tributaries east of the Springerville volcanic field. Most of this differential incision is attributed to uplift in the region of the CR-LCR confluence in the past 2 Ma.

Perhaps the most surprising and impactful advance of this study (and in Hereford et al., 2016) is the discovery that the White Mesa alluvium and paleovalley are ca. 2 Ma, not a Miocene paleoriver (Lucchitta et al., 2013b), nor a tributary to the fluvial ca. 6 Ma upper Bidahochi Formation. Instead this alluvial system records a modest-size tributary to the Little Colorado River analogous to modern Moenkopi Wash. Due to its present gradient and high landscape position, this young age requires that there was relatively little incision in the lower LCR valley between 6 and 2 Ma compared to more rapid incision of the lower reaches of the LCR in the past 2 Ma. The acceleration of regional denudation and incision of the Little Colorado River region that started between 2 and 1 Ma has beheaded the NE extension of the drainages on Black Mesa, dramatically deepened the region of the San Juan–Colorado River confluence, and carved much of the deep LCR gorge below Cameron.

New ⁴⁰Ar/³⁹Ar and U-series geochronology provides late Quaternary bedrock incision constraints throughout the LCR system over the past 6 Ma. Incision rate of the lowermost LCR has kept pace with that of the CR at the confluence over the past ∼1 Ma that has averaged >160 m/Ma (Marble Canyon) to 167–193 m/Ma (lower LCR valley). In the middle reaches of the LCR, elevated basalt flows suggest that acceleration of incision started 2.0–1.5 Ma. Basalt constraints from headwater reaches and eastern LCR tributaries suggest that incision rates have been steady at 20–40 m/Ma for the past 6–0.3 Ma. Uranium-series dates on travertine that cements gravels directly above bedrock straths suggest late Quaternary average incision rates of ∼38 m/Ma from 1 Ma to 100 ka and higher incision rates of 278 m/Ma in the past 100 ka.

At regional scale, the ∼40 m/Ma steady incision rate of the upper LCR is attributed to epeirogenic uplift of the entire Colorado Plateau. The 100 m/Ma higher incision rate of the lower reach and LCR-CR confluence region is interpreted to reflect sub-regional epeirogenic uplift of the eastern Grand Canyon relative to western Grand Canyon as proposed in earlier papers (Karlstrom et al., 2008; Crow et al., 2014). The new data for acceleration of incision in the past 2 Ma in the lower LCR implies a similar change in the CR system, and we propose surface uplift of 50–60 m/Ma near the CR-LCR confluence relative to both the LCR headwaters and central Grand Canyon to explain incision variations. A mechanism for this uplift may involve NE- propagation of mantle-derived magmatism from the San Francisco volcanic field in the past 2 Ma.

Bedrock erodibility is an important control for river system evolution. Regional cliff retreat during the 70–50 and 25–15 Ma time periods was accomplished by rapid retreat across thick weak strata such as the Chile Formation and Mancos Shale. The cliff retreat episodes were a response to base-level changes in trunk rivers. The present LCR and Lees Ferry knickpoints are viewed as bedrock-controlled knickzones (rather than rapidly moving incision transients) because of new data for similar incision rates below, within, and above the knickpoint in the LCR over the past ∼1 Ma.

Tectonic influences on river system evolution remain controversial and require additional testing. At the scale of the LCR-CR system, the lowermost CR has arguably remained near sea level for the past 5 Ma (Crossey at el., 2015; but c.f. Spencer et al., 2013). Differential incision data suggest that the Colorado Plateau has gone up hundreds of meters relative to lower reaches during this time (Karlstrom et al., 2007, 2008; Crow et al., 2014). Likewise, the Rocky Mountain region may have been uplifted several hundred meters relative to Colorado Plateau over the past 10 Ma (Karlstrom et al., 2011). Several regions above the Jemez lineament have been proposed to have been broadly uplifted over the past 5 Ma (Nereson et al., 2013; Channer et al., 2015). This paper postulates surface uplift associated with the San Francisco volcanic field. We favor a component of magmatically driven surface uplift as a mechanism to increase river gradients and introduce knickpoints. Climate changes at 6–5 Ma and ca. 2.6 Ma are likely to have interacted with any tectonic changes to influence regional river integration and incision rates but do not explain differential incision at the scales of the above examples.

This study shows the importance of differential incision studies at regional spatial scale and long (several Ma) temporal scales and the power of combining thermochronologic and geologic data sets. Future progress and additional tests of competing models require “differential thermochronology” studies linked to additional strath-to-strath incision data at numerous localities along the major rivers of the CR system. Precise dating of paleoriver deposits using a variety of techniques is needed, and detrital-zircon and sanidine dating complement each other to constrain maximum depositional ages and provenance of paleoriver deposits. Sanidine dating in particular is shown in this paper to have led to a breakthrough in understanding Cenozoic evolution of LCR and CR river systems of the Colorado Plateau–Rocky Mountain region. Another key broader implication of this work involves quantification of differential uplift magnitude and spatial patterns that can provide input for geodynamic modeling of links between dynamic topography and drainage evolution.

We acknowledge the decades-long contributions of Maurice Cooley on the Little Colorado River valley. Our research was supported by National Science Foundation grant EAR-1348007 to Karlstrom and Shuster. The paper benefitted from reviews by Andre Potochnik, Kelin Whipple, and Stuart Thomson.