Pleistocene to recent geomorphic and incision history of the northern Rio Grande gorge, New Mexico: Constraints from field mapping and cosmogenic ³He surface exposure dating

The Rio Grande of the southwestern United States has long been used as a natural laboratory for understanding fluvial processes, including the mechanisms and timescales of basin capture and drainage integration (e.g., Dethier et al., 1988; Connell et al., 2005; Mack et al., 2006), feedbacks between tectonic processes, river morphology and incision (e.g., Chapin and Cather, 1994; Mack et al., 2011; Repasch et al., 2017), and the role of local and global climatic events on river incision and aggradation (e.g., Dethier, 2001; Pazzaglia, 2005). The Rio Grande is the axial river drainage of the Rio Grande rift (Fig. 1) and flows through structural basins from its headwaters in the San Juan Mountains to its mouth at the Gulf of Mexico. Due to its considerable length, proximity to active structures, and sensitivity to climatic fluctuations, among other factors, the geomorphology of the river changes dramatically along its length, from occupying deeply incised canyons in northern New Mexico to more broad and shallow floodplains farther south.

Within the larger Rio Grande system, the Rio Grande in northern New Mexico is especially unique and enigmatic for a number of reasons. Firstly, the river valley along this reach is narrow and steep walled—in several places a gorge—and it hosts the largest knickpoint of the river’s longitudinal profile, located near the Red River fault zone (RRFZ) (Repasch et al., 2017; Fig. 1D). These morphologic features suggest that this part of the river is in disequilibrium and is actively responding to some form of recent perturbation (Pazzaglia et al., 1998; Whipple, 2001; Bierman and Montgomery, 2013). Secondly, the northern Rio Grande is proximal to several active rift-related structures, including basin-bounding normal faults, and is within a zone of long-wavelength dynamic uplift associated with rift-related asthenospheric upwelling (Moucha et al., 2008; Karlstrom et al., 2012). This active tectonic environment has the potential to influence incision. Third, the northern Rio Grande is the focus of a debate about the age of integration of the northern and southern San Luis Basin in mid-late Pleistocene time (Machette et al., 2013; Repasch et al., 2017). Integration expanded the Rio Grande watershed by ∼22,000 km² (Wells et al., 1987), creating a potential driver for increased incision. Finally, this section of the Rio Grande is proximal to the river’s headwaters and drainage networks in the San Juan and Sangre de Cristo mountain ranges, each relevant sites to understand the effects of regional glacial-interglacial climate change over timescales of 10³–10⁴ years (e.g., Armour, 2002; Pierce, 2003; Pazzaglia, 2005; Guido et al., 2007; Thackray, 2008). These climatic fluctuations have recently received considerable attention as a primary driver of fluvial incision in western U.S. river systems (e.g., Dethier, 2001; Connell et al., 2005; Mack et al., 2011).

These unique aspects of the northern Rio Grande make it a critical place large enough to examine the geomorphic responses of a river system to a remarkable range of external (i.e., climate and large-scale tectonics) and internal (i.e., drainage integration, local tectonics, local sedimentation, and landslide damming) driving forces, yet smaller and with less inherent complexity than, for example, the Colorado River. Here we use detailed field mapping and cosmogenic ³He surface-exposure geochronology of fluvial terraces to examine the incision history of this ∼5 km reach of the northern Rio Grande gorge in New Mexico. We use these observations to determine the timescales, rates, and potential drivers of incision in this section of the gorge, as well as how the proposed incision history here compares to other major river systems in the west-central United States.

This study focuses on a section of the Rio Grande gorge located within the tectonically active southern San Luis Basin of the northern Rio Grande rift in New Mexico (Fig. 1). The San Luis Basin is a 10–100-km-wide, east-tilted, asymmetric half-graben that has experienced left-oblique extension since the late Oligocene, with the most rapid phase of regional extension occurring during the middle to late Miocene (Chapin and Cather, 1994). The basin is bounded on the east by the Sangre de Cristo normal fault and is transected by left-stepping en-echelon normal fault strands of the Red River fault zone (RRFZ) (Fig. 1). The basin is occupied primarily by basaltic lavas of the Taos Plateau volcanic field; these lavas are interlayered with and underlain by Pliocene and late Miocene alluvial and fluvial sediments of the Santa Fe Group. The youngest basaltic lavas, with eruptive ages between 5.2 and 2.8 Ma, the Servilleta Basalt (Lipman and Mehnert, 1979), compose much of the upper portions of the Rio Grande gorge in the study area, with thin layers of unconsolidated to weakly consolidated Neogene and Quaternary sediments locally preserved between and around them (Baldridge et al., 1980; Dungan et al., 1984; Bauer et al., 2015; Repasch et al., 2017; Fig. 1).

A unique aspect of the geomorphic and incision history of the northern Rio Grande River is its hypothesized youthful history, with incipient incision and gorge formation occurring after a “final,” large-scale drainage basin integration event in the northern San Luis Basin sometime between 440 and 800 ka (Machette et al., 2013; Repasch et al., 2017). The northern San Luis Basin was the site of an ancient high-altitude pluvial lake—Lake Alamosa—which persisted from the Pliocene to middle Pleistocene (Fig. 1). Machette et al. (2013) posit that the overtopping of the lake’s hydrologic sill and the southward draining of the lake were likely driven by increased alluvial and lacustrine basin fill from extensive deglaciation during marine isotope stage 11 (MIS 11; 424–374 ka), following what is believed to be one of the most extensive Pleistocene glacial periods recorded (MIS 12; Lisiecki and Raymo, 2005). A different interpretation by Repasch et al. (2017) invokes continuous uplift of the southern Rocky Mountains as the cause for this spillover event and indirectly brackets the timing of spillover to ca. 800 ka based on the presence of San Juan Mountains–derived sediments recovered in the Gulf of Mexico (Galloway et al., 2011; Repasch et al., 2017). Earlier workers postulated that this spillover event occurred between 690 and 300 ka, likely 670 ka, based on the argument that the San Luis Valley did not drain southward into the Rio Grande until deposition of the youngest Alamosa Formation sediments (Wells et al., 1987; Rogers et al., 1992).

Regardless of timing, the overflow of ancient Lake Alamosa likely fully integrated the Rio Grande River from its current upper reaches in Colorado to its earlier headwaters in New Mexico. This event added recharge areas in the San Juan Mountains and northern Sangre de Cristo Mountains, both of which show evidence for Pleistocene glaciation, particularly during the Last Glacial Maximum (MIS 2) (Armour, 2002; Guido et al., 2007). The integration event had profound effects on fluvial dynamics, including but not limited to substantial base-level lowering, pronounced incision, and in particular, knickpoint formation and migration (Machette et al., 2013). These effects were most likely superimposed on an underlying signal of northward headward erosion and stream capture of the Rio Grande into the southern San Luis Basin (Wells et al., 1987).

This study specifically focuses on an ∼5 km portion of the Rio Grande gorge near La Junta Point (Figs. 1 and 2). The gorge here ranges from 230 to 250 m deep (averaging 240 m) and 780–2000 m in width. Gorge width generally increases southward. At and directly (20–80 m) below the gorge rim, cliff faces expose bedrock, whereas landslides, colluvium, and planar geomorphic surfaces are the primary surficial features in the lower ∼180 m of the gorge (Figs. 2 and 3). Bedrock exposures on both the eastern and western gorge rim include multiple flows of Servilleta Basalt, as well as more localized lobate dacitic flows (Figs. 2 and 3; Fig. S1¹). Exposed thicknesses of Servilleta Basalt within the study area range from ∼20–50 m, though the maximum thickness is on the order of 200 m immediately to the south (Ozima et al., 1967; Bauer et al., 2015). Between upper and lower Servilleta Basalt, as well as intercalated within the lower Servilleta, are sedimentary packages of Santa Fe Group silts, sands, and gravels of various thicknesses (Fig. S1 [see ^{footnote 1}]). Thicknesses of sediments range from ∼5–10 m to 70 m in the southern end of the study area (Bauer et al., 2015).

Detailed geologic mapping of gorge surficial deposits was facilitated by a U.S. Geological Survey, airborne light detection and ranging (LiDAR)–derived, 1 m horizontal × 0.12 m vertical resolution digital elevation model, derivative hill shade and slope maps, submeter-resolution aerial photographs, and high-resolution photo panoramas of the canyon walls acquired at five sites on the gorge rims (e.g., Figs. 2, 3, and S1 [see ^{footnote 1}]). These products allowed us to recognize previously undocumented subhorizontal planar or near-planar geomorphic surfaces that we then examined and mapped in detail. Among these surfaces, fluvial terraces were distinguished on the basis of well-defined terrace tread and riser couplets and the presence of associated fluvially modified clasts. Landslide deposits were differentiated on the basis of their spatial extent, morphology, and degree of coherence. Recognition of canyon wall and rim bedrock units relied upon previously published mapping (Kelson et al., 2008; Bauer et al., 2015), augmented by LiDAR hill shades and the photo panoramas.

Cosmogenic surface-exposure dating provides an estimate of the length of time that a rock has been exposed within ∼2 m of Earth’s surface. In situ terrestrial cosmogenic nuclides (TCNs) (e.g., ³He, ¹⁰Be, ²⁶Al, among others) are produced by the interaction of secondary cosmic radiation (primarily neutrons) of different energies with target minerals in exposed rocks (Gosse and Phillips, 2001). Cosmogenic ³He is particularly useful in geologic studies (e.g., Marchetti and Cerling, 2005; Foeken et al., 2009) because it is a stable nuclide that has the highest production rate of all TCNs, as well as a low detection limit on a noble gas mass spectrometer. ³He is produced primarily via spallation reactions on O, Mg, Si, Ca, Fe, and Al within olivine, pyroxene, hornblende, and garnet crystals. Global production rates are internally consistent between olivine and pyroxene and are among the most well constrained of all TCNs (Goehring et al., 2010). We targeted pyroxene within fluvially modified Servilleta Basalt boulders on terrace riser levees, primarily, and terrace treads to estimate the time of terrace abandonment.

Servilleta Basalt clasts were sampled from terrace levees and proximal treads following recommendations outlined in Gosse and Phillips (2001). Levees are commonly constructed during a single flood event during the end of terrace planation (Stewart and Lamarche, 1967), thus emphasis was placed on sampling boulders near or on levees in order to most closely date the timing of abandonment for each terrace, as opposed to the timing of terrace assembly. Samples came from only the upper 10 cm of the most horizontal exposed boulder surfaces (no surfaces sloped more than ∼25°). Only boulders that showed clear evidence for fluvial modification (e.g., sculpting, fluting, smooth, and polished surfaces) were sampled. Boulder morphologies and sizes are provided in the Supplemental Files (see ^{footnote 1}). Four terraces on the more remote western side of the river and one northern terrace were not sampled; at least two samples (more commonly three) from different boulders were collected from all other terrace treads or levees (Figs. 2 and 4). Ten large smooth boulders (>1 m in height) from the modern river plain were also sampled and combined to quantify average inheritance. All dated samples contained variable amounts of pyroxene and olivine (micro)phenocrysts. Pyroxene was targeted for mineral separation and gas analysis due to its comparatively lower degassing temperature (∼1150 °C).

Samples were crushed, milled, and sieved to a 125–425 µm size, followed by removal of ferromagnetic material with a hand magnet and adhering rock dust by rinsing with deionized water. The vast majority of pyroxene microphenocrysts still had volcanic glass attached after these steps. Purification was achieved with two 24 h hydrofluoric-acid leaching steps using a process modified from Bromley et al. (2014). The leaching process has the added benefit of producing more accurate measurements of cosmogenic helium, because leaching removes any non-helium-bearing minerals and groundmass that can be aggregated with the crystals of interest (Bromley et al., 2014). Purity was confirmed via microscope observations. Eighty percent of the samples contained enough pyroxene for analysis.

³He concentrations were determined at Berkeley Geochronology Center (BGC) and the Jackson School of Geosciences (JSG). At BGC, samples ranging from 78 to 150 mg were encapsulated in a tantalum packet and heated under vacuum via a 150 W, 810 nm diode laser. Temperature of the packet was monitored with an optically coaxial pyrometer; all samples reached 1200 °C during a single heating step of 15 min. Re-extraction attempts with a subsequent heating step at 1200 °C for all aliquots revealed an He signal indistinguishable from blank; thus, all gas was assumed to be released in the initial heating step. Extracted gas was then purified by reaction with a SAES getter, frozen with activated charcoal to 11.5 K, and released into a MAP-215 noble gas mass spectrometer at 33 K. The ⁴He signal was measured on a Faraday cup, and the ³He signal was measured on a continuous dynode electron multiplier. Gas amounts were quantified and calibrated with a custom-mixed standard having an absolute ³He/⁴He ratio of 6.116 × 10⁻⁴. Reported uncertainties for the ³He signal represent the reproducibility of multiple mineral aliquots of the same sample and include: the internal measurement precision of each analysis, uncertainty in nonlinearity, uncertainty in blank subtraction (assumed to be negligible), and uncertainty in reproducibility of the custom He standard. Multiple aliquots of individual samples (indicated by a dashed numbering scheme in Table 1) typically reproduced to better than 5%. Samples analyzed at the JSG employed a similar laser system and analytical setup with some exceptions (see Supplemental Files [^{footnote 1}]). Comparisons of interlaboratory results indicate that JSG’s normalized gas amounts are systematically lower than BGC’s, typically by ∼25%. This difference was resolved by application of an empirical correction factor derived from analyses of replicate sample aliquots at BGC (Fig. S2).

³He surface-exposure ages were calculated by applying the empirically derived 1991–2000 Lal-Stone time-independent scaling model to the total ³He concentrations using the CRONUS-Earth ³He exposure age calculator (Lal, 1991; Stone, 2000; Phillips et al., 2016). The ³He nuclide production rate assumed in this model is 118 ± 18 atoms/g/yr and is subsequently scaled to each sample location and elevation. Scaling factors for local topographic shielding are applied in accordance with the data reported in Table 1. Corrections for snow and loess cover are unwarranted because of the relatively high average height of the sampled boulders (>0.5 m). Average winter snow depth within the gorge is difficult to accurately determine; however, at this elevation, it is highly unlikely that snow would have covered the samples deeply or long enough to have an appreciable effect on the exposure age (Gosse and Phillips, 2001). All samples are assumed to have negligible (<1%) magmatic ³He because of the youthful eruptive ages (<5 Ma) of Servilleta Basalt flows. Exposure ages calculated from multiple aliquots of the same samples at both laboratories overlap within error after a correction factor was applied to the JSG data (see Supplemental Files [^{footnote 1}] for details). For samples with multiple aliquots, the sample exposure date was assigned using a weighted mean of all aliquots.

Landslide deposits and colluvium fill most of the lower portion of the gorge (Fig. 2). The gorge in this area is both wider and deeper than upstream and downstream (Fig. 1E), and landslide deposits are concomitantly more prevalent. With rare exceptions, landslide head scarps in the gorge walls cannot be recognized because of mantling talus, but cuspate sections of the gorge rim provide evidence for the location of some.

Landslide deposits span a range of characteristics not individually mapped. Blocky, discontinuous, arcuate to linear ridges and piles of angular basalt and dacite boulders compose much of the hummocky terrain on the west side of the gorge (Fig. 2). Basalt boulders of some deposits have a black shiny varnish that distinguishes them from colluvium or rock falls of more recent origin. The west side hummocky deposits were likely derived from canyon wall slumps that turned into debris flows upon losing strength during descent. At the other end of the spectrum are large, coherent Toreva blocks (Reiche, 1937) on the east side of the gorge that partially retain the stratal layering of the upper gorge walls (Fig. 3). The distinctly planar tops of some of these blocks are capped by columnar jointed and varnished Servilleta basalt that is inclined 30°–50° toward the canyon walls, clearly demonstrating that block displacement was along listric surfaces (Fig. 5). Slip surfaces may have soled and flattened into weaker interlayers of unconsolidated Santa Fe Group sediment and paleosols (i.e., Lama Formation of Bauer et al., 2015), which, except near the gorge rims (Fig. 2), are now entirely covered. V-shaped depressions originally formed between the back-tilted tops of some Toreva blocks, and the gorge wall subsequently served as loci for sediment accumulation. Resultant fill forms near-planar surfaces, some high on the east gorge walls (Fig. 2), that superficially resemble terrace treads. Deposits along the axis of these surfaces are, however, dominantly unstratified aeolian silt, locally to depths of greater than 50 cm, though pebble-sized clasts and sand are also common, with a variable mix of rounded and angular clasts associated with small alluvial fans (Fig. 2) sourced from gorge walls above. Some of these surfaces contain well-varnished and rounded boulders of Servilleta Basalt, but none show evidence of fluvial modification (e.g., river polishing, sculpting, and fluting), having been rounded instead by spalling. We did not observe soils on landslide surfaces or within depressions associated with them.

Seventeen fill or fill-cut terraces were identified within the La Junta Point area, with tread surface areas spanning 1498–43,443 m². Terraces were assigned to six numbered terrace levels, with average elevations above the modern river grade of 7 m (Qt1), 10 m (Qt2), 21 m (Qt3), 28 m (Qt4), 40 m (Qt5), and 50 m (Qt6) (Figs. 2 and 6). Terrace levels Qt5 and Qt6 are the least represented, with only one instance of each along the central-eastern portion of the gorge. Terrace level Qt4 is the most extensive, with multiple instances up, down, and across river, including paired sets in the central portion of the study area (Figs. 2, 3, and 6).

Terraces are built upon risers of alluvium composed primarily of subrounded to well-rounded and water-polished Servilleta Basalt boulders and cobbles typically ranging in size from 0.01 to 1 m³ (Fig. 4), but locally reaching up to 3 m³ in the riser of Qt5, with rounded pebble- to sand-sized sediment filling interstitially. Heights of terrace risers typically range from 4 to 10 m. Terrace risers are nowhere exposed in cross section; so they cannot be described or diagrammed in detail. Where exposed in plan view, risers comprise steep, buttress-like faces of stacked, rounded boulders and cobbles of Servilleta Basalt. Terrace treads are planar surfaces capped by aeolian silt above moderately to very well-rounded, water-polished cobbles and pebbles of Servilleta Basalt, intermediate to silicic volcanic rock, quartzite, and vein quartz in that relative order of abundance. Like terrace risers, tread surfaces are not exposed in cross section, precluding observation of tread profiles or soil development. Isolated, partially buried, well-varnished and rounded basalt boulders are also sparsely distributed on these surfaces but are most concentrated along the margin of the terraces within riser-parallel levees that can exceed a meter in height and width (Fig. 7). The majority of boulders found on each terrace tread, levee, and riser show one or more of the following features as evidence of predepositional fluvial transport: (1) sculpting of surfaces (Fig. 4B); (2) a smooth, river-polished surface on all exposed areas; (3) smooth-walled grooves, pits, and flutes, typically at a small scale (Fig. 4C); and (4) potholes in either upright or random orientations, the former indicating formation in situ.

Terrace dates range from 69.0 +8.4/−9.2 ka to 24.3 +7.6/−6.7 ka, with some terraces showing scatter and others exhibiting clusters and distinctive outliers (Fig. 8). Scatter in cosmogenic data sets of fluvial terraces can be caused by a range of processes, including predepositional exposure history of clasts (i.e., inheritance); prolonged periods of terrace construction and/or abandonment; postdepositional mixing; and/or postdepositional clast erosion (e.g., Repka et al., 1997; Gosse and Phillips, 2001; Dühnforth et al., 2012). In particular, recent studies (Dühnforth et al., 2012; Foster et al., 2017) have documented repeated depositional and erosional cycles within single strath terraces on the Colorado high plains and the consequent difficulties in interpreting terrace tread boulder exposure ages. The strath terraces analyzed in those studies are situated on virtually unconfined plains, where rivers emerging from the mountain front had abundant space to laterally planate and avulse. In contrast, the fill and fill-cut terraces of this study occupy the lower quarter of a nearly linear gorge no more than 2 km wide at its rim (Fig. 1E), leaving little room for significant Rio Grande avulsion and lateral planation. For this reason, and because we sampled only large (typically 0.5 m³), in situ boulders on terrace levees, we consider complex histories of terrace occupation and abandonment, as well as postdepositional mixing as unlikely contributors to the observed scatter, leaving inheritance and boulder spallation and/or erosion as remaining sources.

Boulders from terraces Qt3, Qt4, and Qt6 exhibit dates that appear to be statistical outliers (greater than 2-sigma outside the mean). The outlier on Qt6 is >20 k.y. younger than the other cluster of dates from that terrace, likely reflecting anomalous erosion and/or spallation of this particular sample (Fig. S3 [^{footnote 1}]). By contrast, Qt3 and Qt4 outliers are >20 k.y. older than the clustered means, likely reflecting inheritance. The presence of statistical outliers separated from clustered data differs from some cosmogenic data sets that show log-normal distributions with long tails (e.g., Prush and Oskin, 2016); this suggests that the terraces may record inheritance in erratic boulders (the outliers) as opposed to varying amounts of inheritance affecting all samples.

That inheritance is a factor that in some samples is supported by the amalgamated sample of boulders from the modern river plain, which yields an exposure date of 17.5 ± 2.8 ka. This inheritance value is similar in magnitude to the outlier boulders on Qt3 and Qt4. A simple interpretation would be that the modern river plain records a similar distribution of ages as preserved on the terraces (i.e., a cluster of similar ages—at or close to zero for the modern river), with an outlier boulder or two with substantial inheritance. If this is the case, an anomalously older exposure age is expected after amalgamation. Some studies have suggested using active channel inheritance values to correct for inheritance in terrace ages from the same catchment (e.g., Armstrong et al., 2010; Le Dortz et al., 2011), but more commonly, scatter in the active channel ages is far too large to justify making this correction (e.g., Matmon et al., 2005; Mériaux et al., 2005; Owen et al., 2011; Chevalier et al., 2012; Cording et al., 2014; Gray et al., 2014). Additionally, studies in which other independent geochronometers have been used indicate that this correction substantially overestimates the inheritance value; Blisniuk et al. (2012), for example, demonstrated that correcting cosmogenic ¹⁰Be alluvial fan ages by the average age of active channel samples yielded ages younger than (and therefore incompatible with) the minimum depositional ages constrained by U-series on pedogenic carbonate. Based on this published body of work, we do not consider a correction for inheritance based on the active channel samples to be warranted here.

After removing statistical outliers, the remaining dates were used to estimate abandonment ages for each terrace using probability distribution functions, assuming a Gaussian uncertainty (Zechar and Frankel, 2009; Fig. S4 [see ^{footnote 1}]). Physically, a Gaussian uncertainty model implies that boulder mean ages reflect a balance between clast erosion and inheritance. An alternative interpretation of the data is that even the well-clustered samples may be dominated by small (e.g., ∼5–10 k.y.) amounts of inheritance (with negligible boulder erosion), in which case the youngest exposure age would provide a closer approximation to the true abandonment age. This interpretation would produce ages that fall within error of the Gaussian distributions and therefore would not significantly affect our interpretations of terrace abandonment timing (Fig. 8). We favor the Gaussian distribution interpretation here because overall it is more conservative in estimating errors. Qt6 is thus interpreted to have been abandoned at ca. 69.0 ka, Qt5 at ca. 36.7 ka, Qt4 at ca. 26.9 ka, Qt3 at ca. 25.3 ka, and Qt2 at ca. 24.3 ka. Although dates for Qt4–Qt2 overlap within uncertainty, the progressive younging of ages with decreasing elevation above the modern river lends support to our abandonment age interpretations.

The spatial relationships between large-scale landslide deposits and fill terraces illuminate the relative sequence of events since the onset of fluvial incision in this area. The relative timing of landsliding can be established in a few instances, from which it is clear that landslides both predate and postdate terrace formation. The best examples are near the river, where the youngest slides displace earlier slide deposits and partially bury or remove sections of terraces, some themselves inset into earlier slide deposits (see below and Fig. 5), and the older river alluvium into which terraces were cut (e.g., Fig. 2D). However, the planarity, lack of tilting, and preservation of paleoriver elevation gradients on all terrace treads (Fig. 6) suggest that the preserved terraces themselves have not been displaced by landslides subsequent to their formation. There is also evidence that several major Toreva blocks are located above the sequence of fill terraces; in many places, the projected landslide toes (along their maximum listric slip surface) appear to be truncated by the terraces (Fig. 5). This relationship suggests that the terraces formed and were abandoned after an early, prolonged, but subsequently less active phase of large-scale landsliding and gorge widening that diminished or ceased prior to the abandonment of Qt6 at ca. 70 ka. This phase, accompanied by cliff retreat and progressive gorge widening, has removed any earlier fluvial record.

The size and continuity of the largest landslide deposits on both sides of the river suggest that landslides would have intermittently dammed the river. Direct evidence for this, such as described and documented ∼110 km down river in White Rock Canyon (Reneau and Dethier, 1996a, 1996b) is, however, lacking. Although searched for, lacustrine sediments or other aggradational deposits associated with damming or outburst events are not evident, and fill and fill-cut terraces with levees, as described below, appear inconsistent with such an origin (cf. Reneau and Dethier, 1996a, 1996b; Dethier and Reneau, 1996). We thus suggest that landslide dams either did not persist for sufficient periods of time to leave a record or that such a record was subsequently buried or removed by further landslides prior to terrace formation. We have not identified any evidence that the older alluvial fill and terraces developed upon it are associated with base-level fluctuations resulting from landslide dams.

Based on our interpreted terrace ages, at least two, but likely three periods of aggradation were necessary to generate the ∼50 m of fill that composes the terrace sequence. Qt6 and Qt5 are each relatively limited in spatial extent and have older exposure ages than the lower four terrace levels; thus it is highly unlikely that all terrace fill was accumulated during a single period of aggradation. Because of the substantial difference in age and height between Qt6 and Qt5, these terraces likely represent individual aggradation and incision events. Qt4 is also likely to represent an individual aggradational event based on its extensive preservation and width. The short duration over which Qt4 through Qt2 were abandoned (∼1–5 k.y.) suggests they may be fill-cut terraces incised into a single extensive fill deposit represented by Qt4. While these ages indeed overlap within uncertainty, the rapid incision suggested by these ages is not surprising given the relatively small differences in elevation between them and that incision was through alluvial fill rather than bedrock. The alluvial fill required to generate Qt4–Qt2, if confined to this widest ∼5 km reach of the gorge, comprises a volume of ∼0.01 km³ deposited within 4–5 k.y.

To examine whether the terraces described here record a response to local and/or global climate change, we compare their abandonment ages to both the global marine oxygen-isotope stage (MIS) glacial-interglacial climate record (Lisiecki and Raymo, 2005) and to more proximal indicators of western U.S. climate change (e.g., Gosse et al., 1995b; Benson et al., 2005; Bhattacharya et al., 2018). In considering the global record, fill terrace formation and abandonment have often been attributed to climate-driven fluctuations in hillslope sediment supply and river transport capacity during glacial to interglacial transitions (Bull, 1991, 2000; Vandenberghe, 2003). In particular, several workers have suggested that glacial and/or pluvial periods should correspond to a fluvial regime in which sediment supply overcomes stream power, leading to net aggradation, whereas the opposite regime dominates during interglacial and/or interpluvial periods, leading to net incision (Bull, 1991; Hancock and Anderson, 2002), even in the absence of glaciers within the headwaters (Pierce and Scott, 1983). However, this relationship is complex and is intimately related to a river system’s sediment budget and supply, which can be influenced or even dominated by non-climatic factors such as base-level change (reviewed in Blum and Törnqvist, 2000), hillslope strength (e.g., Carroll et al., 2006), river avulsion (e.g., Dühnforth et al., 2012; Foster et al., 2017), or tectonic tilting (e.g., McMillan et al., 2002; Riihimaki et al., 2007). Dühnforth et al. (2012) and Foster et al. (2017) provide particularly well-documented examples of the decoupling of glacial-interglacial cycles and aggradation and/or incision regimes resulting from relatively rapid incision events following prolonged periods of terrace development. As stated above, because river avulsion at the time of terrace formation must have been highly restricted within this reach of the Rio Grande, it is unlikely that is a complicating factor in this instance.

Of relevance to the terrace abandonment ages in this study are MIS 4 through 1, with start dates at 71 ka, 57 ka, 29 ka, and 13 ka, respectively (Fig. 9). As demonstrated in Figure 9, four out of the five terraces dated have abandonment ages that overlap with glacial periods MIS 4 or MIS 2. If we assume that aggradation prior to incision through fill was relatively rapid, terrace formation in our study area is generally consistent with enhanced aggradation during glacial periods, with the exception of Qt5, which appears to represent an aggradation and incision event within interglacial MIS 3 (although its age uncertainties are the largest of all the terraces, and it has the lowest sample size).

In Figure 10, we compare our terrace sequence ages with a compilation of other quantitatively dated geomorphic features (e.g., terminal moraines, recessional moraines, striated bedrock, and fill terraces) in the greater Rocky Mountain region. Several of the ¹⁰Be dates in this compilation were originally calculated using higher production rates than are in use presently; therefore, we recalculated them using a ¹⁰Be production rate of 4.01 atoms/g/yr (Borchers et al., 2016), with resulting age increases of ∼12%–18%. The compilation demonstrates that, as discussed by previous authors, evidence for regional glaciation during MIS 4 is sparse, and few glacial deposits have been identified. Nonetheless, some lake highstands and vegetation records (e.g., Jimenez-Moreno et al., 2007) and outwash terraces (e.g., Sharp et al., 2003) at higher latitudes, along with fill terraces (including our Qt6) at lower latitudes, overlap within error with this period and may record glaciation between the better preserved Pinedale and Bull Lake glaciations. The ∼36 k.y. age for the Qt5 terrace preserved in our study area is anomalous relative to other deposits in the region, with the exception of a fill-cut terrace preserved along a tributary of the Colorado River (Nankoweap) (Anders et al., 2005) and a strath terrace preserved along the Colorado Front Range (Foster et al., 2017). Terraces Qt4–Qt2 all appear to be closely associated with the LGM and/or Pinedale glaciation and associated increase in runoff, discharge, and sediment flux. They are coeval with Pinedale fill terraces dated in Boulder Canyon along the Colorado Front Range (Schildgen et al., 2002) and with terminal moraines in the Wind River Range (Gosse et al., 1995b) and the Yellowstone area farther north (Licciardi et al., 2001). The Qt4 terrace overlaps within error but may slightly predate a recessional moraine and outwash terrace sequence preserved in the Animas Valley to the northwest in the San Juan Mountains (Fig. 1; Guido et al., 2007). Modeling of moisture variability in the southwestern United States during the LGM based on a robust leaf-wax biomarker proxy record from cores retrieved from the Gulf of California indicates that North American monsoon-derived rainfall decreased during this time and is largely tied to glacial-interglacial cycles of North American ice cover extent and associated atmospheric patterns (Bhattacharya et al., 2018), thus complementing the regional record and trends observed farther north (Fig. 10). These timing relationships and the extensive development of the Qt4 terrace in our study area suggest it reflects punctuated incision in response to the onset of deglaciation and the early stages of Pinedale ice retreat.

We first consider constraints on early time-averaged bedrock incision over the first ∼190 m of gorge incision using two brackets: the onset of gorge incision and the timing of development of the highest terrace (Qt6) riser. The onset of incision is bracketed by either (1) the age of Lake Alamosa overflow at 440 ka (Machette et al., 2013) or 670 ka (Wells et al., 1987; Rogers et al., 1992); or (2) the 800 ka oldest age of northernmost Rio Grande sediments recovered in the Gulf of Mexico (Galloway et al., 2011; Repasch et al., 2017). The younger age bracket is approximated by the timing of abandonment of the highest terrace in the gorge (Qt6 at ca. 69 ka). The estimated incision rates using these brackets are ∼512 m/m.y., 316 m/m.y., and 260 m/m.y. from 800 to 69 ka, 670–69 ka, and 440–69 ka, respectively (Fig. 11A). Independent of the age bracket chosen, these minimum incision rates for the northern Rio Grande over this time interval are significantly higher than those estimated for nearby river systems, including the Rio Chama, the Colorado River, and the Wind River (Fig. 11A).

In Figure 11B, we examine the late Pleistocene to Holocene (<126 k.y.) incision rates as constrained by fill terraces for the Rio Grande and nearby river systems. Consistent with an earlier compilation (Dethier, 2001), rates computed over these timescales are substantially faster than those estimated over the middle Pleistocene. This observation was interpreted by Dethier (2001) to reflect increased amplitude of pluvial cycles during late Pleistocene time but may in part also be due to a “Sadler effect” (Sadler, 1981; Pederson et al., 2006; Gallen et al., 2015). Nonetheless, the average incision rate estimated for the Rio Grande using all available data is nearly twice the incision rates calculated for the Colorado and Wind Rivers (Fig. 11B), thus mimicking the trend observed for longer timescales shown in Figure 11A. This average rate is also very similar to the rate computed using the age and distance from the Qt6 terrace to the modern river plane in our study area (∼724 m/m.y.). A more conservative, but still high, rate of ∼515 m/m.y. is estimated for the time period between Qt6 and Qt4, which, as discussed in the Relationships between Terraces and Global and Regional Climate Proxies section, likely represents climate peaks during MIS 4 and MIS 2, respectively. Independent of the methods used, these calculations suggest that Rio Grande incision rates have been significantly faster than adjacent river systems from the middle Pleistocene to the present day.

The fast incision rates calculated for the Rio Grande since the middle Pleistocene suggest an underlying driver for incision that is unique to this river system. Potential sources include proximal tectonic activity, dynamic topography, drainage basin integration, transient knickpoint migration, and regional climate oscillations. Active slip along normal fault systems oriented at high angles to the Rio Grande, such as the Red River fault system, has the potential to locally increase stream gradient. However, there is no direct evidence of late Pleistocene slip on this structure (Ruleman et al., 2013). Joint seismic-geodynamic modeling of density anomalies beneath the Colorado Plateau and surroundings does predict concentrated dynamic uplift beneath the eastern flank of the Rio Grande rift (Moucha et al., 2008), but the estimated uplift rates of ∼35 m/m.y. are too slow to account for the much faster incision estimated for the Rio Grande. Additionally, these same uplift rates should also affect the Rio Chama, yet much slower middle Pleistocene incision rates are estimated there (Fig. 11A; Dethier and McCoy, 1993). As discussed above, planation of, and incision through, the Qt4 terrace and/or fill appears to closely track the start of Pinedale ice retreat during MIS 2, pointing to climate oscillations as a driver for incision in our study area. However, although the amplitude of these oscillations is increasing through time and may be inducing overall faster incision for all river systems (e.g., Dethier, 2001), there is no reason this effect should differentially influence the Rio Grande. This analysis leaves capture of the San Luis basin during drainage basin integration and concomitant transient knickpoint migration as the most likely explanation for the fast incision rates that characterize the Rio Grande. Integration of the northern San Luis Basin, independent of whether it represents an integration or re-integration event, and independent of the precise timing, likely had a substantial impact on incision rates through a significant increase in stream power (Wells et al., 1987; Machette et al., 2013; Repasch et al., 2017), despite the likelihood of an initial release of stored sediment post-overflow. Dynamic response to this process is likely still ongoing as evidenced by the ∼135 m knickpoint located ∼8 km north of the Rio Grande–Red River confluence (Machette et al., 2013; Fig. 1D). The current knickpoint nearly coincides with where the gorge abruptly widens southward, from ∼500 m to more than 1000 m in width (Fig. 1E). This widest gorge reach, which is not associated with any marked changes in the lithology of the gorge walls, persists southward to the Red River confluence, where it begins to gradually narrow, diminishing to less than 500 m ∼6 km south of the confluence. A migrating knickpoint, initiated near the present Red River confluence, with an associated steeper river gradient and greater stream power, must have played a role in increased transient rates of incision and gorge widening, which in turn provided space for large landslides and later terrace development, both hallmarks of this widest section of the gorge.

Cosmogenic ³He surface exposure dating and LiDAR-based geomorphic mapping of fill and fill-cut terraces in the Rio Grande gorge near Questa, New Mexico, provide constraints on the timing and rates of river incision and gorge development. Early gorge incision was punctuated by large landslides that generated coherent Toreva blocks and less coherent mass wasting deposits. The narrower, lower ∼50 m of the gorge preserve fill and fill-cut terraces interpreted to record at least three aggradation-incision geomorphic cycles over the past ∼70 k.y. The formation and abandonment timing of four of the five terrace levels overlap with global marine isotope stages 4 and 2. The Qt6 terrace adds to other cryptic evidence from the southern Rockies for regional glaciation during MIS 4. The formation and incision of terraces Qt4–Qt2 appear coincident with the early stages of Pinedale ice retreat in the Rocky Mountain region. Mid-Pleistocene to present bedrock and sedimentary and/or bedrock incision rates inferred from independent estimates of the timing of initial gorge incision are >260 m/m.y. during the middle Pleistocene and >515 m/m.y. from the late Pleistocene to the present. These rates are approximately twice as fast as rates estimated over the same timescales for the Colorado River, Wind River, and Rio Chama, suggesting an incision driver that is unique to the Rio Grande system. The most plausible driver for ongoing fast incision is increased stream power following capture of the northern San Luis Basin between 440 and 800 ka, along with knickpoint migration through this reach of the river as a consequence of this drainage basin integration event, with dynamic uplift and climate oscillations as superimposed but lower magnitude signals.

We thank Desmond Patterson, Greg Balco, Marissa Tremblay, and Daniel Stockli for assistance with ³He cosmogenic measurements. David P. Dethier and an anonymous reviewer provided extensive comments and suggestions that substantially improved the manuscript. Ren Thompson, U.S. Geological Survey, provided airborne LiDAR data prior to public release. We also are very grateful for several discussions in the field with Ren, Cal Ruleman, Peter Gold, and Paul Bauer.