The impact of bedrock meander cutoffs on 50 kyr scale incision rates, San Juan River, Utah

Incision rates of major rivers may reflect the effects of drainage reorganization, hillslope processes, tectonic uplift, climate, the properties of rocks into which rivers incise, and other autogenic processes. On the Colorado Plateau, incision rates along the Colorado River have been interpreted as resulting from abrupt base-level changes produced by the integration of the Colorado River system. Specifically, the integration of the Colorado River in the location of Grand Canyon is thought to have created a knickpoint, enhanced by lithologic contrasts, which is propagating upstream. While evidence exists for a < 1 Ma acceleration of incision on parts of the Colorado River, uncertainty about the processes reflected in shorter-term incision rates muddies comparison with longer-term averages. In this work, we combine a cosmogenic radionuclide depth profile exposure age and post-infrared infrared stimulated luminescence (p-IR IRSL) to date fluvial deposits adjacent to the San Juan River, a major tributary of the Colorado River, near Mexican Hat, Utah. The deposits, resting on a 27 m strath surface, are constrained to be ∼ 28–40 ka, suggesting an incision rate of 678–968 m Myr⁻¹, which is a factor of 7 higher than the long-term rate of ∼ 125 m Myr⁻¹ over the past ∼ 1.2 Myr. We observe fluvial deposits that were abandoned due to a bedrock meander cutoff, which partially explains our accelerated incision rate. We use a simple geometric model, informed by our field data, to demonstrate how planform river evolution may, in some circumstances, increase short-term incision rates relative to long-term averages. These short-term rates may reflect a combination of autocyclic and climatic processes, which limits their ability to resolve longer-term changes in incision rates that may be related to changes in base level or tectonics.

1 Introduction

The landscape of the Colorado Plateau in the American Southwest is known for deep canyons, broad plateaus separated by steep escarpments, and mesas that protrude from otherwise low-relief surfaces. Regional denudation and the resulting landforms are closely tied to the rate of incision of large rivers, which serve as major conduits for mobilized sediment. Despite the centrality of the Colorado River and its tributaries to understanding the complex geologic, tectonic, and climatic history that shaped the Colorado Plateau, significant uncertainty still exists around the fluvial evolution of the region and how the major rivers came to their present configuration.

Integration of the Colorado River and its tributaries is generally believed to have occurred 5 to 6 million years ago when a combination of integration through Grand Canyon incision (Karlstrom et al., 2014), basin spill-over (Blackwelder, 1934; Pearthree and House, 2014), and groundwater sapping (Crossey et al., 2015) connected the Colorado Plateau to the Gulf of California. Reconstruction of how tributaries became integrated and responded to downstream changes in base level relies on analysis of dateable geologic material. These alluvial and fluvial deposits have been used to constrain the incisional history and evolution of the Green (Aslan et al., 2010; Darling et al., 2012), Upper Colorado (Aslan et al., 2019), and Little Colorado rivers (Karlstrom et al., 2017) over the last 1.5 Myr. They have been used to infer the regional tectonic evolution of the Colorado Plateau and the Rocky Mountains (Karlstrom et al., 2012; Pederson et al., 2013; Rosenberg et al., 2014; Aslan et al., 2019). For example, Darling et al. (2012) synthesized multiple decades of studies documenting incision on the Colorado River and its tributaries. Over the > 1 Myr timescale, incision rates below a knickpoint at Lees Ferry, Arizona, are ∼ 170–230 m Myr⁻¹ (Pederson et al., 2002; Karlstrom et al., 2008; Polyak et al., 2008), which is higher than the ∼ 110–130 m Myr⁻¹ rate above the knickpoint (Darling et al., 2012). Cook et al. (2009) found rates on a tributary of the Colorado River as high as 350–600 m Myr⁻¹. These spatial patterns have been interpreted to represent a change in base-level conditions following integration through Grand Canyon (Wolkowinsky and Granger, 2004; Cook et al., 2009), regional tectonism (Karlstrom et al., 2008, 2012), or contrasts in rock properties (Cook et al., 2009; Bursztyn et al., 2015).

Over shorter timescales (< 0.5 Ma), incision rates at Lees Ferry are estimated to be ∼ 350 m Myr⁻¹ over 150 kyr. This higher incision rate, compared to the million-year-scale rate, has been attributed to isostatic feedback from drainage integration (Pederson et al., 2013), although this is contested by Karlstrom et al. (2013). Short-term incision rates along the Upper Colorado River have also yielded higher rates than those averaged over the > 1 Myr timescale, with some as high as > 3700 m Myr⁻¹ where drainage reorganization appears to play a role (Aslan et al., 2019). Others have noted that short-term incision rates can be particularly sensitive to glacial oscillations, which contribute to fluctuating sediment supply, discharge, and isostatic effects in rivers with montane or alpine headwaters (Hancock and Anderson, 2002) as well as the cyclical nature of bedrock incision and aggradation in rivers (Finnegan et al., 2014).

This work focuses on the incision history of the San Juan River, which is a major tributary in the Colorado River system. Lucchitta and Holm (2020) proposed that this tributary reflects an upstream expansion of the Colorado River system following its integration. In this view, an early course of the lower San Juan River formed by exploiting a posited Oligocene to late Pleistocene “Crooked Ridge River” that flowed through the Kletha Valley (Cooley et al., 1969; Hunt, 1969), although this has been rebutted by Heizler et al. (2021) based on ages of deposits that are inconsistent with that location for a paleo-San Juan River.

Discerning the recent incision history of the San Juan River is challenging due to a dearth of dated fluvial or alluvial material in the region. Two notable exceptions to this are the cosmogenic ²⁶Al–¹⁰Be burial dating of sediments near Bluff, Utah (Wolkowinsky and Granger, 2004), which were also recently dated using detrital sanidine (Heizler et al., 2021) and the exposure age dating of surfaces grading to the present location of the San Juan River along the flanks of Navajo Mountain (Garvin et al., 2005). These incision rate estimates at Bluff (> 110 ± 14 and ≥ 125 m Myr⁻¹ for cosmogenic isotopes and detrital sanidine, respectively) and at Navajo Mountain (∼ 400–825 m Myr⁻¹) are separated by ∼ 120 km, employ different surface exposure dating methodologies, and sample different time intervals of the incision history of the San Juan River. This has led to some uncertainty as to whether the differences in incision rates between the two sites represent a methodological bias, a regional spatial variation in incision rate, a regional temporal fluctuation in erosion rate (Darling et al., 2012), or a bias toward higher incision rates over shorter timescales (Finnegan et al., 2014).

Interestingly, Pleistocene river terrace surfaces suitable for constraining the short-term incision rate of the San Juan River downstream of the Bluff site exist, but few have been characterized or dated in the scientific literature. This presents an opportunity to constrain the incision history of the San Juan River over different periods of time at two proximal sites, allowing us to study (1) if the San Juan River short-term (< 500 kyr timescale) and long-term (> 500 kyr timescale) rates are in agreement and (2) what mechanisms may produce disagreement, should it exist. We use geomorphic mapping, in situ cosmogenic ¹⁰Be depth profile analysis, and post-infrared infrared stimulated luminescence (p-IR IRSL) dating methods to constrain the rates and mechanisms that produce the short-term incision history along the San Juan River near Mexican Hat, Utah, ∼ 45 km downstream from Bluff. We find that deposition of the terrace materials occurred during or after 39.8 ± 3.1 ka, which is consistent with the ¹⁰Be depth profile exposure age of 27.9 ± 0.8 ka. The implied 10 kyr incision rate measured in this study (∼ 678–968 m Myr⁻¹) greatly exceeds those measured over the ∼ 1.20–1.63 Myr timescale (84–125 m Myr⁻¹). The modern San Juan River bypasses the dated material, which lies in a cutoff meander entrenched into resistant bedrock. We establish a conceptual model for how meander cutoff can result in a discrepancy in the apparent long-term and short-term incision rates and explore the effect of lithology on formation of planform erosional features, which can preserve datable geologic material, with implications for biases when measuring vertical incision.

2 Geologic setting

2.1 Regional geology

The San Juan River is a major tributary of the Colorado River and drains ∼ 64 000 km² in southeastern Utah, northeastern Arizona, northwestern New Mexico, and southwestern Colorado (Fig. 1). The headwaters of the San Juan River are in the formerly glaciated San Juan Mountains, from which the river flows ∼ 600 km to its confluence with the Colorado River at Glen Canyon. The upper San Juan River occupies a broad valley attributed to Miocene incision (Cooley et al., 1969; Hunt, 1969), though recent controversy about the evolution of the San Juan River has brought into question how and when the San Juan River established its modern course (Hereford et al., 2016; Lucchitta and Holm, 2020; Heizler et al., 2021).

Figure 1 Shaded relief map of the Four Corners region (southwestern Utah). Major physiographic, geographic, and fluvial features are labeled. The red dot indicates the location of sampling sites for Wolkowinsky and Granger (2004) and Heizler et al. (2021). Inset shows San Juan drainage in a yellow polygon with the Colorado River and Green River shown as blue lines. The hollow red box in the detailed map shows the area of Fig. 2. AZ – Arizona, CO – Colorado, NM – New Mexico, UT – Utah, WY – Wyoming.

This broad valley becomes inset by a narrow gorge just downstream of Chinle Creek, Utah, that extends to the confluence with the Colorado River at Glen Canyon, Arizona (Fig. 1). This gorge runs through the Monument Upwarp, a large structural fold dominantly comprised of Pennsylvanian and Permian marine sedimentary rocks with discrete zones of Laramide-style, fault-cored monoclines and anticlines with corresponding synclines (Ziony, 1966). The San Juan River incises through a number of these features, including Comb Ridge monocline, which bounds the Monument Upwarp on the east, and Raplee Ridge monocline (Fig.  1). Downstream of Raplee Ridge, the deeply incised meanders of the Goosenecks led early researchers to infer that the San Juan River predated the Laramide Monument Upwarp (Wengerd, 1950), though Moore (1926) identified some other meanders on the San Juan River that may have formed in place where bounded by more erodible bedrock. Subsequent work has corroborated the fact that incising river meanders are not exclusively inherited from preceding alluvial rivers but may form and change shape in actively incising bedrock settings (Tinkler, 1971; Shepherd, 1972; Barbour, 2008; Finnegan and Dietrich, 2011; Limaye and Lamb, 2014).

We focus on the San Juan Canyon between Raplee Ridge monocline and Gooseneck State Park, near Mexican Hat, Utah (Fig. 2). Raplee Ridge is a ∼ 14 km long monocline with a north–south-oriented, doubly plunging fold axis. The fold forming Raplee Ridge has been variably described in the literature as a monocline and anticline, though here we opt for Raplee Ridge monocline as the topographic expression and fault kinematics (Hilley et al., 2010) are similar to other Laramide monoclines in the area. The core of the monocline is comprised of interbedded erodible shale and resistant limestone layers of the Pennsylvanian Rico and Paradox Formations (Fig. 2), which are exposed in gullies that dissect the monocline surface (Steelquist et al., 2023) as well as where the monocline is bisected by the San Juan River. Downstream of Raplee Ridge, the San Juan River is bounded by the Permian Halgaito Tongue Formation mudstone for ∼ 11 km before again incising into the Rico and Paradox Formations as it enters the Goosenecks west of Mexican Hat, Utah.

Figure 2 Geologic map of the Mexican Hat (Utah) region, modified from O'Sullivan (1965). Boundaries of structure-from-motion surveys and sampling locations are in the center of the map.

2.2 Sample site – Mexican Hat quarry

We focused our study on the area surrounding a gravel quarry 1.6 km northeast of Mexican Hat, Utah, north of the San Juan River (Figs. 2 and 3a). This quarry was active at the time of sampling but is now inactive. O'Sullivan (1965) briefly described these fluvial deposits, saying “… there are at least four levels of terrace gravels, but all cover small areas…. Parts of the town of Mexican Hat are built on the most extensive terrace in this area, which is about 125 feet above river level.” Figure 3d shows an example of the context of these deposits in the San Juan Canyon.

Figure 3 (a) Orthographic images created using structure-from-motion techniques, including approximate image locations (black dots) and look directions (black arrow) of subsequent panels. (b) Fluvial gravel sampled for ¹⁰Be analysis. The outcrop height in the image is approximately 2 m. (c) Poorly consolidated fluvial sand outcrop sampled for IRSL dating, resting on fluvial gravel; the sand height in the image is approximately 3 m. (d) View looking approximately east from the gravel quarry across the San Juan River. Note the multiple bedrock straths (red-brown) mantled by fluvial gravel (gray). Fluvial sands overlying the higher gravel appear as light tan patches where not hidden by colluvium. Raplee Ridge is in the background.

The fluvial deposits at the Mexican Hat quarry are part of a fill terrace resting on a strath surface 27 m above the modern San Juan River surface. The lower part of the deposit is composed primarily of medium to coarse, clast-supported cobbles in a sand matrix (Fig. 3b). Gravel deposits are increasingly undisturbed by mining activity towards the San Juan River, where the exhumed top forms a cobble-dominated gravel pavement. Gravel deposits are ∼ 3 m thick but thin towards a small gully that bisects the deposit.

The gravel quarry is bounded on the northeast by a mesa of Halgaito Tongue Formation and Cedar Mesa sandstone. At the base of this mesa are eroded remnants of sand deposits ∼ 4 m thick (Fig. 3c). Sand deposits throughout the region are largely attributed to aeolian processes (O'Sullivan, 1965); however, this sand variably overlies the sandy gravel and bedrock strath surface. Exposures contain subhorizontal beds with plane bedding, relatively massive upward-fining units, ripples, and silt drapes, which suggest a fluvial origin. This is further bolstered by the absence of high-angle cross-bedding that would indicate aeolian deposition. The sand appears inset into the gravel at some locations, plausibly filling a paleo-channel in the gravel top; however, disruption by quarrying operations means this relationship cannot be traced throughout the area. Similar sandy deposits are more widespread nearby, suggesting sand deposits originally covered the sandy gravel at the quarry, making the gravel pavement an eroded rather than depositional surface. As cosmogenic radionuclide depth profile dating is sensitive to erosion of the deposit surface, which is possible given our observations, we used multiple dating techniques to determine the deposit age.

3 Methods

3.1 Structure-from-motion (SfM) surveying and geomorphic mapping

We completed three aerial surveys of the San Juan River directly downstream of Raplee Ridge using an uncrewed aerial vehicle (UAV). Each survey was processed in Agisoft Metashape (Version 1.6.5) to form three-dimensional red–green–blue-labeled point clouds. To ensure consistency across surveys and proper georeferencing, we established ground control points using a high-precision GPS (lateral precision ∼ 5 cm, vertical precision ∼ 20 cm) for each survey region and used these points to reference image locations during Metashape processing.

We used the SfM point clouds to create paired digital elevation model (DEM) and orthographic photos for each survey region. Using these data in combination with field observations, we mapped fluvial deposits and fill terraces within the San Juan Canyon between Raplee Ridge and Mexican Hat. All fluvial deposits include sandy pebble to cobble gravel, often with scattered boulders, that rests directly on bedrock strath surfaces. Some of this gravel is overlain by several meters of fluvial sands and slightly gravelly sands, resembling those at the quarry. We separated the relative age of fluvial deposits based on the height of the underlying strath surfaces above the modern San Juan River in the SfM datasets. We also mapped gravel deposits in the San Juan Canyon between Raplee Ridge and Mexican Hat using satellite imagery outside of the SfM survey regions. However, the low resolution of topographic data for the region (∼ 10 m per pixel and > 1 m vertical accuracy) and human disturbance at Mexican Hat allow only speculative correlations to be made between these terrace levels and those captured by SfM.

3.2 Cosmogenic radionuclide dating – reassessing Bluff terrace ages

¹⁰Be and ²⁶Al are primarily produced through spallation of Si and O when material is exposed to cosmic-ray bombardment at (or near) Earth’s surface (Lal and Peters, 1967; Nishiizumi et al., 1991), though some contribution from negative-muon capture and fast-muon reactions also occurs (Heisinger et al., 2002a, b). As neutrons rapidly attenuate with depth, the concentration of ²⁶Al–¹⁰Be in rock material can be used to determine the amount of time that material has been at the surface (Lal, 1991). For fluvially transported material, this is typically assumed to be the combination of ²⁶Al–¹⁰Be produced during erosion from its source region (typically referred to as the “inherited concentration”) and ²⁶Al–¹⁰Be produced after terrace abandonment (Granger and Muzikar, 2001). However, if buried deeper than this attenuation depth, the relative abundance of ²⁶Al in previously dosed sediment will decrease relative to ¹⁰Be because of the shorter half-life of the former relative to the latter (Lal, 1991). While individual clasts have different erosional histories and thus different inherited ²⁶Al–¹⁰Be concentrations, well-mixed fluvial sands represent catchment average denudation rates (Granger et al., 1996). Given this averaging effect, in assessing in situ ²⁶Al and ¹⁰Be, we assume that sand sampled from multiple depths shares a single inherited concentration, simplifying the process of dating fluvial material.

Previous work at the Bluff site (Fig. 1) exploited the prior accumulation of ¹⁰Be and ²⁶Al by clasts during exhumation and transport, the near-surface accumulation of ¹⁰Be and ²⁶Al after deposition, and the more rapid decay of ²⁶Al at depth to jointly solve for the exposure age, erosion, and inheritance of clasts (Wolkowinsky and Granger, 2004). These researchers collected seven clast samples from terrace fill gravel > 13 m thick, with samples taken between the surface and 11.7 m depth. In some cases, these samples were amalgamations of 20–30 quartzite clasts. However, since the time of that initial study, the production and decay rates for these cosmogenic isotopes have been revised (e.g., Lifton et al., 2014); thus, we first updated the results from Wolkowinsky and Granger (2004) using modern parameter values to ensure consistency between the estimation of a long-term incision rate at Bluff and our estimation of a short-term incision rate at nearby Mexican Hat (see Sect. 3.3).

The procedure for computing ages and incision rates is documented in Wolkowinsky and Granger (2004). Briefly, we calculated expected concentrations of ¹⁰Be and ²⁶Al given a terrace age, terrace erosion rate, individual clast inheritance (assuming an ²⁶Al to ¹⁰Be production ratio of 7.112), bulk density, and loss of material from the terrace surface as follows (Wolkowinsky and Granger, 2004, their Equation DR1):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& N=N_{\mathrm{inh}}{e}^{\frac{-t}{\mathit{\tau }}}+\left[{\displaystyle \frac{P_{\mathrm{n}}{e}^{-\frac{\mathit{\rho }z}{\mathrm{\Lambda }}}}{(\frac{\mathrm{1}}{\mathit{\tau }}+\frac{\mathit{\rho }\mathit{\epsilon }}{\mathrm{\Lambda }})}}\right]\left[\mathrm{1}-e^{-t(\frac{\mathrm{1}}{\mathit{\tau }}+\frac{\mathit{\rho }\mathit{\epsilon }}{\mathrm{\Lambda }})}\right]\\ & +\left[{\displaystyle \frac{Y{A}_{\mathrm{1}}{e}^{-\frac{\mathit{\rho }z}{L_{\mathrm{1}}}}}{(\frac{\mathrm{1}}{\mathit{\tau }}+\frac{\mathit{\rho }\mathit{\epsilon }}{L_{\mathrm{1}}})}}\right]\left[\mathrm{1}-e^{-t(\frac{\mathrm{1}}{\mathit{\tau }}+\frac{\mathit{\rho }\mathit{\epsilon }}{L_{\mathrm{1}}})}\right]\\ & +\left[{\displaystyle \frac{Y{A}_{\mathrm{2}}{e}^{-\frac{\mathit{\rho }z}{L_{\mathrm{2}}}}}{(\frac{\mathrm{1}}{\mathit{\tau }}+\frac{\mathit{\rho }\mathit{\epsilon }}{L_{\mathrm{2}}})}}\right]\left[\mathrm{1}-e^{-t(\frac{\mathrm{1}}{\mathit{\tau }}+\frac{\mathit{\rho }\mathit{\epsilon }}{L_{\mathrm{2}}})}\right]\\ & +\left[{\displaystyle \frac{B{e}^{-\frac{\mathit{\rho }z}{L_{\mathrm{3}}}}}{(\frac{\mathrm{1}}{\mathit{\tau }}+\frac{\mathit{\rho }\mathit{\epsilon }}{L_{\mathrm{3}}})}}\right]\left[\mathrm{1}-e^{-t(\frac{\mathrm{1}}{\mathit{\tau }}+\frac{\mathit{\rho }\mathit{\epsilon }}{L_{\mathrm{3}}})}\right],\end{array}\end{array}$$

where N (atm g⁻¹) is the measured concentration of ²⁶Al or ¹⁰Be, N_inh is the inherited concentration for ²⁶Al or ¹⁰Be, t (years) is the exposure age, z (cm) is sample depth (allowing for addition material to all samples), τ (yr⁻¹) is the radioactive mean life of ²⁶Al or ¹⁰Be, P_n (atm g⁻¹ yr⁻¹) is the production rate of ²⁶Al or ¹⁰Be by nucleon spallation, ρ (g cm⁻³) is the density of the overburden, and ε (cm yr⁻¹) is the erosion rate of the terrace surface. Following Wolkowinsky and Granger (2004), we used Λ = 160 g cm⁻², A₁ = 170.6, and A₂ = 36.75 at sea level and high latitude; L₁ = 738.6 g cm⁻², L₂ = 2688 g cm⁻², L₃ = 4360 g cm⁻², Y_Al = 4.24 × 10⁻³, and B_Al = 0.192 for ²⁶Al; and Y_Be = 4.91 × 10⁻⁴ and B_Be = 0.023 for ¹⁰Be. We updated production rates for the sample location and elevation from the time-dependent, nuclide-dependent scaling scheme of Lifton et al. (2014) estimated using the CRONUScalc online interface (Marrero et al., 2016a). A₁ and A₂ were scaled for latitude and elevation using a factor of 1.5231 (derived from CRONUScalc). We also adjusted the measured ¹⁰Be concentrations by a factor 1.0351 given the effect of the revised ¹⁰Be concentration in NIST SRM4325 from Nishiizumi et al. (2007), as this standard was used to correct the original Bluff ¹⁰Be measurements.

We then used the reported measured abundances of ¹⁰Be and ²⁶Al to define a misfit between calculated ¹⁰Be and ²⁶Al concentrations as follows:

$$\begin{array}{}\text{(2)}& \mathit{\varphi }={\displaystyle \frac{\mathrm{1}}{v}}\left[\sum {\left[{\displaystyle \frac{N_{\mathrm{Al}}^{\mathrm{pred}}-N_{\mathrm{Al}}^{\mathrm{meas}}}{{\mathit{\sigma }}_{\mathrm{Al}}^{\mathrm{meas}}}}\right]}^{\mathrm{2}}+\sum {\left[{\displaystyle \frac{N_{\mathrm{Be}}^{\mathrm{pred}}-N_{\mathrm{Be}}^{\mathrm{meas}}}{{\mathit{\sigma }}_{\mathrm{Be}}^{\mathrm{meas}}}}\right]}^{\mathrm{2}}\right],\end{array}$$

where ϕ is the misfit, v is the degrees of freedom (6 for these data), and σ is the concentration measurement error for each element. We numerically minimized this misfit using a steepest-descent method to solve for the model parameters that best fit the observed data. As with Wolkowinsky and Granger (2004), we then set all model parameters equal to their best-fitting values, while varying terrace erosion rate and surface exposure age to provide an estimate of the uncertainties in these parameters. The uncertainties are equal to the minimum and maximum values of erosion rate and exposure age for which ϕ is within 1 standard error of the best-fitting values.

3.3 Cosmogenic radionuclide dating – ¹⁰Be depth profile

We used in situ produced ¹⁰Be in quartz grains to determine the age of the fluvial deposits at the Mexican Hat gravel pit. We sampled quartz sand at five depths ranging from 5 cm to 2 m below the present top of the gravel (each sample was taken from a ∼ 2 cm horizon) and processed the samples at Stanford University using a modified procedure from Kohl and Nishiizumi (1992). Full sample preparation description can be found in the Supplement. All samples were analyzed using the accelerated mass spectrometer (AMS) at Lawrence Livermore National Laboratory in Livermore, California, using the 07KNSTD3110 standard with assumed ¹⁰Be $/$ ⁹Be = 2.85 × 10⁻¹².

We used the resulting ¹⁰Be $/$ ⁹Be ratio from the AMS analyzes to determine the concentration of ¹⁰Be at each sample depth and fit a line with the form of Eq. (1) using nonlinear least-squares minimization, solving for exposure age (t) and ¹⁰Be inheritance (N_inh). We assumed no terrace erosion (addressed in Sect. 4.3) and used the bulk density calculated from our reassessment of the Bluff terrace data.

3.4 Infrared stimulated luminescence (IRSL) dating

Minerals such as quartz and feldspar contain luminescence signals that result from radioactive decay in naturally occurring radioisotopes in the surrounding material, with signal accumulating during burial (Huntley et al., 1985; Rhodes, 2011). As sediment grains are transported by geomorphic processes, exposure to sunlight bleaches the luminescence signal. Thus, we can use luminescence dating of buried sediments to infer the time since the last exposure to sunlight (Rhodes, 2011; Buylaert et al., 2012). While previously applied to amalgamated fluvial sediment, recent advances in measuring single-grain luminescence signals in many individual grains taken from a single deposit yield a distribution of luminescence ages (Arnold et al., 2007; Rhodes, 2015; Gray et al., 2018). The single-grain age distributions provide the last depositional age and may also provide information about the depositional setting (e.g., Saha et al., 2021). For example, aeolian transport will likely expose all grains to sufficient sunlight to bleach all luminescence signals, resulting in a unimodal population of single-grain ages. In contrast, transport in a silty fluvial setting like that in the San Juan River may leave a population of grains unbleached, leading to a multi-modal age distribution due to residual luminescence from a more complex bleaching history (e.g., Gliganic et al., 2017). Using single-grain IRSL dating techniques provides insight into the depositional character and timing for the sands at the quarry. In combination with the cosmogenic isotope depth profile, which dates the exposure of the gravel layer after abandonment, we can better constrain stratigraphic relationships and timings of terrace fill deposition and channel abandonment (see Sect. 4.3).

We used post-infrared infrared stimulated luminescence (p-IR IRSL; Reimann et al., 2012; Rhodes, 2015) analysis in single-grain K-feldspars to determine the distribution of last sunlight exposure in grains within the poorly consolidated, fluvial sand at the Mexican Hat quarry. We sampled sand ∼ 1 m above the gravel deposit by inserting an opaque 5 cm diameter metal tube into a freshly cleaned face of the sand deposit and immediately capping it. K-feldspar grains of 175–200 µm diameter from the sedimentary samples were isolated at the UCLA Luminescence Laboratory using the procedures in Rhodes (2015) under dim amber LED light conditions. We then processed the samples through wet sieving and utilized lithium metatungstate to separate the feldspar grains based on density ρ < 2.565 g cm⁻³ (Rhodes, 2015). Luminescence measurements were measured using a TL-DA-20 Risø automated reader equipped with a single-grain IR laser (830 nm, at 90 % of 150 mW) and ⁹⁰Sr $/$ ⁹⁰Y beta source (Bøtter-Jensen et al., 2003). Emissions were detected using a photomultiplier tube with the IRSL signal passing through a Schott BG3–BG39 filter combination. Full descriptions of dose-rate calculation, post-IR IRSL protocol, and anomalous fading test can be found in the Supplement.

We estimated equivalent doses for 158 K-feldspar grains (Table 2) and determined single-grain luminescence age populations following the methods detailed in Saha et al. (2021). The samples showed a high degree of overdispersion (abbreviated OD, 42 %) between single-grain De. Such high OD values may indicate incomplete bleaching during their last transport and depositional event (Fig. 7a, Table 2, Rhodes, 2011). Following the recommendation of Rhodes (2011, 2015), a conservative OD of 15 % was used as a cutoff between well-bleached and partially bleached samples. We therefore used a three-parameter minimum age model (MAM) after Galbraith et al. (1999) assuming a standard OD (i.e., sigmab) of 15 % in the R statistical package for luminescence dating (Kreutzer et al., 2012). The post-IR IRSL ages were calculated in the DRAC 1.2 online calculator (Durcan et al., 2015a). We assumed that the MAM represents the most recent depositional age of the sand deposits at the gravel quarry (Fig. 7).

4 Results

4.1 Geomorphic mapping

Quaternary geologic and geomorphic mapping resolves multiple stages of incision in the San Juan Canyon between Raplee Ridge and Mexican Hat. We identified multiple eroded fill terraces higher than the terrace sampled for ¹⁰Be cosmogenic radionuclide and p-IR IRSL analyses, as well as multiple lower strath terraces. Figure 4a–c show the preliminary results of both SfM orthographic imagery and satellite imagery geomorphic mapping. Fluvial deposits can be found on both banks of the San Juan River; however, straths are typically unpaired in any single reach.

Figure 4 (a) Satellite imagery of the San Juan Canyon near Mexican Hat with both SfM and satellite-derived geomorphic map overlays. Satellite imagery collected 4 October 2017, © Maxar Vivid Imagery, accessed through ESRI World Imagery. (b) Detailed geomorphic map of the upper canyon created using both field observation and SfM-derived orthographic imagery. (c) Detailed geomorphic map of the lower canyon and gravel quarry created using both field observation and SfM-derived orthographic imagery. ¹⁰Be and IRSL dating was performed on fluvial terrace Q3. (d) Strath surface elevations are relative to the modern San Juan River in SfM-derived DEM and orthographic imagery mapping.

The paired orthographic imagery and DEMs from the SfM surveys allow us to classify fluvial deposits based on the elevation of the strath surfaces on which they rest. We identify a series of at least six terraces with elevations ranging from 10–60 m above the modern San Juan River (Fig. 4d). All strath heights are the difference between the elevation of the San Juan River and the elevation of exposed bedrock strath surfaces in the SfM models. Q1 is found in the upstream portion of the surveyed area and is ∼ 13 m above the modern San Juan River. Q2 (∼ 20 m above the modern San Juan River) is found on the southeast side of the San Juan River across from the Mexican Hat quarry and parallels the modern San Juan River. Q3 (∼ 27 m above the river) and Q4 (∼ 41 m above the river) are the most extensive, paralleling the San Juan River in some locations and occupying the ends of minor valleys orthogonal to the modern San Juan River (Fig. 4a). The Q3 fluvial deposits were dated using p-IR IRSL and cosmogenic radionuclide techniques, where they overlie a strath surface 27 m above the modern San Juan River. Two higher but deeply eroded deposits, Q5 (∼ 55 m above the river) and Q6 (∼ 60 m above), mantle isolated hills in the upstream reach of the survey regions and may be either lag gravel remnants of fluvial deposits or colluvium. Q3 is the thickest fluvial deposit at ∼ 3–7 m, while the others are typically several meters but locally can be <1 m. Additionally, we used a combination of field observations and Google Earth to map fluvial material around Mexican Hat. These units are denoted with “Q?” and use ordered letters moving up from the river (“A” is lowest elevation, “B” is next lowest elevation, etc.) as we lack sufficient data to determine exact elevations or correlate them with our deposits in the SfM area. Some of these deposits are older than Q6, with strath heights 100 m above the river.

The location and spatial pattern of fluvial deposits within the study area resolve two cutoff bedrock meanders in direct proximity to our sampling location (Fig. 5a). Eroded fluvial deposits are exposed throughout the smaller meander (southeast of the modern San Juan River) and in parts of the larger one (northwest of the modern San Juan River) where not buried by younger, locally derived alluvium (Fig. 5a). Near the Mexican Hat Rock quarry these cutoff meander features contain Q3, suggesting Q3 predates the meander cutoff. Within this cutoff meander there is no inset strath terrace, and its southeast edge has excellent exposure of the contact between the Q3 strath and its fill. The strath is subhorizontal and has some minor undulations, but we did not observe any evidence of subsequent incision after deposition of the fill, though gullying crosscuts both fill and strath on the southwest edge of the gravel pit. The major strath surface below is Q2, which parallels the modern San Juan River in this area and predominantly appears on the opposite bank of the river from the gravel quarry (Fig. 5a). We interpret this to be the first strath formed after the meander cutoff. Given this, we interpret the meander cutoff to have resulted in the abandonment of Q3. The modern river in the bypass reach flows south, whereas before cutoff, flow in this reach was toward the north. This is the only documented location along the San Juan River where such a local flow reversal has occurred.

Figure 5 (a) Shaded relief map of San Juan Canyon near Mexican Hat, Utah. Mapped fluvial deposits are shown with proposed migrating channel locations (blue lines), the San Juan River prior to meander cutoffs in this study (black dashed lines), and approximate meander cutoff locations (yellow stars). Digital elevation model and derived products created using the $\mathrm{1}/\mathrm{3}$ arcsec US Geological Survey National Map 3D Elevation Project (3DEP, https://www.usgs.gov/core-science-systems/ngp/3dep, last access: 22 September 2021) (b) Shaded relief map of the Mexican Hat region. Cutoff bedrock meanders are shown as dashed blue lines. The red box denotes the extent of Fig. 5a. Digital elevation model and derived hillshade created using the 1 arcsec US Geological Survey National Map 3D Elevation Project (3DEP, https://www.usgs.gov/core-science-systems/ngp/3dep, last access: 22 September 2021).

Between Raplee Ridge and the Mexican Hat Rock quarry, the SfM and satellite imagery mapping suggest bedrock channel migration to the west–northwest, outer bank of the meander, and down-stepping bedrock “scroll bars” along the inside of the first meander bend in the San Juan River downstream of Raplee Ridge (Fig. 5a). There is also a northeast-migrating down-stepped series of strath surfaces on the opposite bank downstream of this meander.

Farther upstream, there are multiple bedrock meander cutoffs (Fig. 5b), particularly where the San Juan River incises into the resistant Pennsylvanian strata (Leopold and Bull, 1979; Harden, 1990). These meanders display far greater ridge-to-channel relief than those downstream, which likely reflects the difference in erodibility between the resistant Pennsylvanian limestones upstream of Raplee Ridge monocline and the mudstones of Halgaito Tongue Formation downstream of the monocline. The sinuosity of the deeply entrenched meander cutoffs is similar to those found in the Goosenecks downstream of Mexican Hat, though no abandoned meanders are observed in that stretch of the river.

4.2 Reassessing Bluff terrace dating

Our reassessment of the Bluff terrace data yielded an exposure age of $t=\mathrm{1.63}\frac{+\mathrm{0.08}}{-\mathrm{0.07}}$ Ma, a terrace erosion rate (i.e., the rate of erosion of the terrace surface itself by processes such as rain splash, wind transport, and surface deflation) of $\mathit{\epsilon }=\mathrm{9.8}\frac{+\mathrm{0.8}}{-\mathrm{0.6}}$ m Myr⁻¹, a deposit density of ρ = 1.77 g cm⁻³, and the addition of no material to the terrace surface (Fig. 6a). Clast inheritance values ranged from 0.112–1.622 × 10⁶ atm g⁻¹ (Fig. 6a inset). This increase in age suggests a long-term incision rate of 84–101 m Myr⁻¹, a slightly lower incision rate than the original estimate for this location of 110 m Myr⁻¹. Our best-fit terrace erosion rate (not to be confused with the river incision rate) overlaps within error with the 14 ± 4 m Myr⁻¹ best-fit rate from Wolkowinsky and Granger (2004) but is slight lower than previously estimated. An estimated density of 1.77 g cm⁻³ and the addition of no eroded material are both reasonable. We discuss the accuracy of this incision rate in more detail in Sect. 5.2.

Figure 6 (a) Recreation of best-fit terrace erosion rate and age from Bluff data (Wolkowinsky and Granger, 2004). Contours are 1 standard deviation. The star represents the best-fit value. Sample location details in Wolkowinsky and Granger (2004). No shielding factor was applied. (b) Depth profile of ¹⁰Be concentrations with error bars. Dotted lines represent 1σ error of best fit as determined by nonlinear least squares. Sample taken at 37.16205° N, 109.84611° W and 1278 m elevation. No surface erosion or shielding factors were applied.

4.3 Mexican Hat terrace age dating

Figure 6b shows ¹⁰Be concentrations as a function of depth at the gravel quarry near Mexican Hat, Utah (Table 1), which yielded an exposure age of 27.9 ± 0.8 ka (R² = 0.996) with an inherited ¹⁰Be concentration of 0.712 ± 0.032 × 10⁵ atm g⁻¹, assuming no erosion has taken place since the exposure of the top of the observed deposit (an assumption addressed below). The catchment-averaged inheritance value is consistent with the 0.27–5.11 × 10⁵ atm g⁻¹ range of ¹⁰Be inheritance values from Wolkowinsky and Granger (2004).

Table 1 Cosmogenic radionuclide depth profile data.

Sampling location details and common values between samples: location (° N/° W) = 37.16205/109.84611, elevation (m a.s.l.) = 1278, thickness (cm) = 2, shielding factor = 1, denudation (mm yr⁻¹) = 0, density^f (g cm⁻³) = 1.77, surface production rate (atm g⁻¹ yr⁻¹) for spallation = 10.1979 and muon = 0.1781, LSDn scaling scheme, taken from CRONUScalc (Marrero et al., 2016b)
^a Depths measured from the top of gravel pavement surface. ^b Isotope ratios were standardized to ¹⁰Be standards prepared by Nishiizumi et al. (2007), 07KNSTD3310, with a value of 2.85 × 10⁻¹² and ¹⁰Be half-life of 1.39 × 10⁶ years. ^c Uncertainties reported at 1σ level. ^d Propagated uncertainties include error in the blank. ^e Blank correction was performed by subtracting ¹⁰Be atoms in procedural blank from all samples. ^f Estimated density for unconsolidated fluvial sediments.

The MAM for the sample yields an equivalent dose (De) of 141.8 ± 9.9 Gy and a p-IR IRSL age of 39.8 ± 3.1 ka (Table 2; Fig. 7a). High overdispersion (42 %) in our sample likely indicates incomplete bleaching from sunlight exposure, a result common from fluvial settings (Fig. 7a, e.g., Gliganic et al., 2017). Our single-grain analysis suggests that the sands are likely fluvial in nature rather than aeolian as originally assumed in O'Sullivan (1965), though other sands that do appear reworked by wind can be found in the region. This is further bolstered by the dominantly subhorizontal < 1 mm bedding in the sands (Fig. 7b) and absence of high-angle cross-bedding that would indicate aeolian deposition. We assume that the 40 ka from MAM represents the most recent depositional age of the fluvial sands, while older grains represent remobilized fluvial material whose luminescence signal was not reset during transport.

Figure 7 (a) Radial plot showing single-grain equivalent. Circles represent individual grains. Filled circles are analyses that fall within the minimum age model results (shaded bar). (b) Horizontal, < 1 mm scale bedding in sand deposit.

Table 2 Dose-rate information and post-IR IRSL ages based on the minimum age model (MAM).

Sampling location details: location (° N/° W) = 37.1636/109.8450, elevation (m a.s.l.) = 1278, depth (m) = 3.2. Grain size used: 175–200 µm. Radionuclide conversion factor after Liritzis et al. (2013); α attenuation factor after Brennan et al. (1991); β attenuation factor after Guérin et al. (2012). Internal K contents were 12.5 ± 0.5 % after Huntley and Baril (1997). Cosmic dose rates following Prescott and Hutton (1994). ^a U, Th, and K contents derived via ICP-OES with relative uncertainties of 5 %. ^b Ages are calculated in DRAC calculator (Durcan et al., 2015b). ^c No fading.

The ages of sand deposition, as measured by p-IR IRSL, and gravel exposure, as measured by ¹⁰Be depth profile, are coherent, with deposition preceding exposure. This finding is somewhat counterintuitive, however, given that the sands overlie the gravel. As the top of the sand deposit is ∼ 4 m higher than the top of the gravel pavement surface, we believe the gravel was largely shielded from cosmogenic rays by the sand deposit. After meander cutoff and channel abandonment, erosion of the Q3 fill would eventually remove the sand and expose the gravel, leading to the observed gravel pavement surface that serves as the top of the cosmogenic isotope depth profile. We view this scenario as likely and believe the ¹⁰Be depth profile exposure age likely underestimates the time of the meander cutoff. Nonetheless, the meander cutoff could have occurred no earlier than the age of deposition of the units that were bypassed (∼ 40 ka) and no later than the exposure age of the gravel package (∼ 28 ka), thus bracketing the age of the meander cutoff.

The 27 m strath elevation above the modern San Juan River suggests a short-term bedrock incision rate between 678 and 968 m Myr⁻¹ over the last ∼ 28–40 kyr, though the higher bound likely overestimates incision as it relies on the gravel exposure age. As our vertical datum is the water surface rather than the riverbed, these incision rates could increase by 5 %–10 %. Additionally, there are undulations in the Q3 strath surface moving away from the sampling locations that could decrease the incision rate by 5 %–10 %. Regardless of these uncertainties, this range of incision rates is higher than the that inferred for the San Juan River ∼ 45 km upstream in Bluff, Utah, particularly given that our revision age of the Bluff deposit from ∼ 1.36 to ∼ 1.63 Ma adjusts the range of possible incision rates to 84–125 m Myr⁻¹. Our calculated incision rate is also on the high end compared to rates measured for the San Juan River near Navajo Mountain (425–800 m Myr⁻¹; Garvin et al., 2005).

5 Discussion

5.1 Paleo-channel geometry, preservation of meander cutoffs, and short-term incision rate

Although meanders are typically associated with alluvial rivers, active bedrock meanders are well-documented in studies of bedrock channel strath formation (Hovius and Stark, 2001; Finnegan and Dietrich, 2011). Multiple locations within our field area (Fig. 5a) show evidence of down-stepping, unpaired strath terraces on the inside of bends, typically a sign of bedrock meandering (Personius et al., 1993; Finnegan and Dietrich, 2011). As noted by Finnegan and Dietrich (2011), bedrock meander cutoffs result in knickpoints that can propagate upstream, which would enhance local vertical incision.

Variations in bedrock erodibility have affected the relief between the channel and surrounding uplands along the river. Cutoff bedrock meanders with floors in the erodible Halgaito Tongue are typically entrenched ∼ 40 m (but up to 100 m) below surrounding mesas in the Cedar Mesa sandstone, the cutoff meanders of this study being a primary example. In contrast, those with floors in the Hermosa Formation (see cutoff bedrock meanders upstream of Raplee Ridge in Fig. 5b) are typically ∼ 200 m below plateaus formed by the resistant Rico Formation. Also, the active channel and better-preserved cutoff meanders in the Rico–Hermosa area have rather symmetric cross-valley profiles, whereas some in the Halgaito Tongue–Cedar Mesa area are markedly asymmetric. Given that fluvial bedrock erosion consists of both vertical and lateral erosive processes, these topographic differences imply a higher rate of vertical relative to lateral erosion in resistant layers, allowing considerable deepening prior to meander cutoff. Conversely, meanders in erosive lithologies appear to incise less before lateral migration causes a cutoff in the meander. This implies a higher frequency of meander cutoffs in erodible lithologies compared to resistant lithologies, resulting in higher variability in measured incision rates depending on the spatial and temporal proximity to cutoff events. It is plausible that bedrock meanders upstream from Raplee Ridge deepened primarily through vertical incision, while in resistant limestones of the Rico Formation lateral migration and cutoff occurred when the river encountered more erodible shale layers in the underlying Hermosa Formation. This would require higher-resolution elevation data and field observation to confirm.

To determine whether the process of meander cutoff can explain the observed discrepancy in short- versus long-term incision rates observed along the San Juan River, we constructed a simple geometric model in which the instantaneous cutoff of a meander strands sediments residing at the upstream end of the meander. In this construction, the preferred wavelength of an initially straight channel is denoted as L and evolves to a particular sinuosity, ζ, at which time the channel bypasses the meander (Fig. 8a). We further assume that the slope of the channel outside the meander (S) remains approximately constant during meander formation (Sieh and Jahns, 1984). Finally, we assume that the length of the bypassing segment is small relative to L.

Figure 8 (a) Schematic diagram of meander cutoff and evolution through time with length of straight channel (L) and channel segment upstream of cutoff (l_d). (b) Plot of time since meander cutoff (t_co) versus the difference between long- and short-term incision rates (ϵ and ϵ_m, respectively) for various channel slope (S) and sinuosity (ζ) combinations.

With these assumptions, the stream-wise length of the channel segment above the meander cutoff (l_d) is

$$\begin{array}{}\text{(3)}& l_{\mathrm{d}}=L\left(\mathit{\zeta }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right),\end{array}$$

which implies an elevation difference across the meander at the moment prior to cutoff of

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }h=SL\left(\mathit{\zeta }-\mathrm{1}\right).\end{array}$$

In this case, Δh represents the elevation difference that results solely from the bypass of the meander bend. From the time of the cutoff until the present (t_co), the total accumulated elevation between the channel bottom and the top of the deposit is

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }h+\mathit{ϵ}{t}_{\mathrm{co}},\end{array}$$

where ϵ is the long-term incision rate. Thus, the measured short-term incision rate (ϵ_m) inferred from the age of the deposit abandoned by cutoff and the current river elevation is

$$\begin{array}{}\text{(6)}& {\mathit{ϵ}}_{\mathrm{m}}={\displaystyle \frac{\mathrm{\Delta }h}{t_{\mathrm{co}}}}+\mathit{ϵ}\end{array}$$

such that the difference between the short- and long-term incision rates will be

$$\begin{array}{}\text{(7)}& {\mathit{ϵ}}_{\mathrm{m}}-\mathit{ϵ}={\displaystyle \frac{\mathrm{\Delta }h}{t_{\mathrm{co}}}}={\displaystyle \frac{SL(\mathit{\zeta }-\mathrm{1})}{t_{\mathrm{co}}}}.\end{array}$$

When the cutoff has happened relatively recently (t_co is small), the short-term incision rate will be strongly influenced by the cutoff process, whereas as time passes (t_co is large), ϵ_m will approach ϵ. The magnitude of the discrepancy at a given value of t_co will vary in proportion to the amount of relief created by the cutoff event (SL(ζ−1); Fig. 8b). For small cutoff meanders with low relief from the upstream to downstream end, the discrepancy in erosion rates will only be appreciable for a short time period, while large cutoff meanders with significant elevation difference across the ends may impact erosion rate discrepancies for 10–100 kyr.

Importantly, each of the terms in Eq. (7) can be coarsely measured or inferred based on field observations, which allows evaluation of the effect of meander growth and cutoff at particular sites. The current value of S between Comb Ridge and Mexican Hat is 1.2 × 10⁻³ based on measurements from the Cortez 1:24 000 Quadrangle (the San Juan River undergoes an elevation drop of 48.8 m over 40.5 km of river profile). Likewise, we approximate the initial meander wavelength of the San Juan River in this area to be around 1.2 km based on the current channel geometry and the cutoff length observed near Mexican Hat. This being the case, the length of the two mapped cutoff meanders in Fig. 5a (∼ 9.1 km) implies a local sinuosity of around 7.6. Finally, we use the ¹⁰Be exposure age dating (∼ 28 ka) to approximate the timing of cutoff of the meander (t_co). These values would suggest a difference between the incision rate since the meander was cut off and a long-term incision rate of ∼ 342 m Myr⁻¹. Using the measured short-term incision rate (678–968 m Myr⁻¹) and the range of measured long-term incision rates (84–125 m Myr⁻¹, depending on dating method), our calculations yield only ∼ 39 %–65 % of the observed discrepancy between short- and long-term incision rates. While appreciable, this effect cannot entirely reconcile the observed divergence in incision rates.

Late Cenozoic climatic cycles are known to result in fluvial terraces in both glaciated and unglaciated catchments worldwide (Westaway et al., 2009). The cycles of incision and aggradation observed both prior to and after meander cutoff near Mexican Hat suggest that other exogenous effects, plausibly glacial–interglacial cycling, modulate the patterns of alluviation and thus incision rate over time. This phenomenon has been specifically invoked to explain the long-term terrace record in the broader Colorado Plateau (e.g., Bridgland and Westaway, 2008). The San Juan Mountains have an extensive history of glaciation (Atwood and Mather, 1912; Johnson et al., 2017), so sediment supply and alluviation in the San Juan River must have been impacted by glacial advance and retreat in the headwaters. The Rocky Mountains lack absolute age constraints for the onset of the last glaciation, making it difficult to directly relate the Mexican Hat fluvial deposit to a particular glaciation, and our ages do not clearly align with any marine isotope stage typically associated with glaciation. Nonetheless, there has been a series of documented landscape adjustments in the San Juan Mountains (Johnson et al., 2011) that occurred after the documented meander cutoff. Such adjustments may reflect variations in runoff, sediment supply, and sediment caliber that might also enhance incision rates downstream. Thus, a combination of local autocyclic processes and regional changes in sediment delivery and runoff to the area together might explain the difference between short- and long-term incision rates along the San Juan River at the Mexican Hat site.

5.2 Long- and short-term San Juan River incision rates relative to the Colorado River and tributaries

The short-term incision rate near Mexican Hat is significantly higher than estimates for the San Juan River near Bluff and is closer to the range observed near Navajo Mountain at the confluence of the San Juan and Colorado rivers (400–825 m Myr⁻¹). By combining dating methods that record deposition (p-IR IRSL) and exposure (in situ cosmogenic ¹⁰Be depth profile), the timing of channel abandonment at the quarry must be between ∼ 27.9 and 39.8 ka, representing an incision rate of 678–968 m Myr⁻¹. This is nearly an order of magnitude higher than the cosmogenic burial isochron ages at the Bluff terrace, originally calculated as 110 m Myr⁻¹ and reassessed here to be 84–101 m Myr⁻¹. Detrital sanidine dating by Heizler et al. (2021) shows a maximum depositional age for the Bluff terrace fill of 1.208 Ma. As the deposit cannot be younger than the grains that comprise it, we take this to be the most robust age. We suggest that the older cosmogenic burial isochron age may be the result of more complex exposure and burial history (as suggested by Hanks and Finkel, 2005) or may represent the onset of deposition rather than terrace abandonment (as noted by the original authors in their reply to Hanks and Finkel). We take the maximum depositional age from detrital sanidine to calculate the minimum long-term incision rate on the San Juan River in this area as 125 m Myr⁻¹, which is closer to our short-term incision rate but still slower by at least a factor of 5.

Our incision rate estimate falls on the upper end of 250–725 m Myr⁻¹ short-term incision rates for the upper Colorado River and its tributaries (Darling et al., 2012; Aslan et al., 2019). As noted above, there are multiple potential explanations for the observed differences in incision rate: methodological biases from exposure age dating, spatial or temporal fluctuations in erosion rate due to drainage integration or tectonism, and biases toward higher incision rates when measuring over short timescales (Finnegan et al., 2014). Given the combination of dating techniques used in this work, we view the first of these (methodological biases) as an unlikely explanation for the discrepancy, at least at the Mexican Hat site.

The greater Colorado River system appears to have had changes in incision rates that may be related to tectonism or drainage reorganization. Short-term incision rates in the Glen Canyon area show an acceleration of incision rate from ∼ 120 to ∼ 300–900 m Myr⁻¹ over the past ∼ 250–500 kyr (Garvin et al., 2005; Cook et al., 2009; Darling et al., 2012), with this acceleration being attributed to the combination of an upstream-propagating incisional wave resulting from Grand Canyon integration and the increased erodibility of the Glen Canyon Group sandstone immediately downstream of Lees Ferry. Darling et al. (2012) specifically noted the uncertainty of the state of such an incisional wave for the San Juan River. They highlighted that while the Bluff terrace (150 m above the modern channel) measured by Wolkowinsky and Granger (2004) corresponds in elevation to those of terraces at Bullfrog near Glen Canyon, long-term incision rates at Bluff were lower than at Bullfrog (116–138 m Myr⁻¹; Darling et al., 2012). However, the detrital sanidine dating at Bluff brings these incision rates into alignment, suggesting an incisional wave may have reached this stretch of the San Juan River and be contributing to elevated short-term incision rates.

Simple detachment-limited bedrock channel incision modeling has shown that lithologic contrasts may enhance upstream propagation of such an incisional wave once it encounters more erodible strata, increasing short-term incision rates (Cook et al., 2009; Darling et al., 2012). Cook et al. (2009) documented such accelerated incision rates in Trachyte Creek, a minor tributary of the Colorado River that enters Glen Canyon roughly ∼ 135 km upstream of the Colorado’s confluence with the San Juan River. Given the similarity of soft-over-hard lithologic contrast between Paleozoic limestones and Mesozoic sandstones along the San Juan River upstream of the confluence with the Colorado River, it is reasonable to assume a knickpoint would propagate similarly due to equivalent lithologic contrasts.

Differences in short-term and long-term incision rates have been observed elsewhere on the Colorado Plateau (Aslan et al., 2019). In the upper Colorado River basin, short-term incision rates show significant variation, even within individual regions, while long-term incision rates are both lower and more consistent with one another (Aslan et al., 2019). Aslan et al. (2019) invoked tectonic and isostatic uplift of the flanks of the Rocky Mountains to explain the spatial and temporal incision rate trends observed. However, in the absence of other evidence for such uplift in the Mexican Hat area, we conclude that this explanation would be insufficient to explain our increased incision rates.

This leaves other factors to produce higher incision rates measured over shorter timescales. Finnegan et al. (2014) noted the sensitivity of incision rates to measurement interval, citing the non-steady-state nature of rivers experiencing alluvial aggradation resulting from factors including climate and sediment supply. Additionally, Aslan et al. (2019) pointed to local factors such as bedrock lithology, glacial history, and drainage reorganization as modulating incision rates, with the largest incision rates exceeding > 3000 m Myr⁻¹ on the Gunnison River. Meander cutoffs provide a novel mechanism by which this “Sadler effect” may occur, complementing previous research on biases in incision rates from strath terraces (Gallen et al., 2015; Finnegan et al., 2014) and on stochastic events that alter channel geometry (Ouimet et al., 2008). Given that the sampled deposits are within a cutoff bedrock meander and proximal to a second such meander, enhanced local incision is probably partially responsible for the almost order-of-magnitude increase we observed at Mexican Hat relative to the long-term rate at the Bluff site. As noted above, the process of meander growth and eventual cutoff may increase short-term incision rates relative to their long-term average, as the shortened channel length after the cutoff may juxtapose sediments deposited on the upstream end of a meander with the bypassed segment of the meander. Our work supports meander cutoff effects accounting for up to 50 % of incision rate discrepancy near Mexican Hat; however, the remainder may be attributable to combinations of top-down (climatic–sedimentological) and bottom-up (drainage integration and tectonic) processes.

5.3 Broader implications

Gilbert (1877) first described the connection between the amount of alluviation and a river’s ability to vertically incise, coining the term “planation” to describe the lateral carving of broad plains when vertical incision is slow. Specifically, he noted that “… downward wear ceases when the load equals the capacity for transportation. Whenever the load reduces the downward corrasion to little or nothing, lateral corrasion becomes relatively and actually of importance” (p. 126). Researchers have documented evidence for lateral erosion in bedrock systems, including their tendency to form meanders, for over a century (e.g., Winslow, 1893; Tinkler, 1971; Stark et al., 2010). The fluid mechanics of meander formation have predominantly focused on alluvial systems (summarized by Seminara, 2006), but similar models have been applied to bedrock rivers as well (Turowski, 2018).

Experiments by Shepherd (1972) demonstrated the effect of alluvial cover on bedrock meanders, showing that alluviation on the inner edge of the meander resulted in a shift towards higher erosion rates on the outer edge. More recent mechanistic models of mixed bedrock–alluvial systems have demonstrated that rates of bedrock meander development (in addition to other important aspects of channel geometry) are significantly affected by patterns of alluviation (Inoue et al., 2017; Turowski, 2018; Inoue et al., 2021). These studies found that increases to initial alluvium thickness decreased vertical erosion while increasing the rates of lateral erosion. Rock tensile strength (one potential measure of erodibility) impacts rates of both vertical and lateral erosion in the models of Inoue et al. (2021). However, vertical erosion rates are much more sensitive to the areal fraction of alluvial cover on the channel bed than lateral erosion rates (Shepherd, 1972). Under alluvial deposition, lateral migration would be expected to continue even if vertical incision is significantly impeded by alluviation; however, lateral erosion rates would vary significantly depending on the lithologic susceptibility of the bank to erosion (Stark et al., 2010). Applying this to a river such as the San Juan with headwaters that have experienced glaciation suggests that the availability of alluvium due to glacial–interglacial cycling should impact meander migration rates and sinuosity (Turowski, 2018). Gilbert (1877) specifically noted the patterns of alluviation and planation in soft rocks along Trachyte Creek, showing that maintenance of a constant river profile resulted in transport-limited conditions in erodible lithologies, creating higher lateral channel migration. This also agrees with the observation of Moore (1926), who noted that entrenched meanders along the San Juan River are best preserved in resistant lithologies but are less commonly observed in erodible lithologies where the “… meandering pattern [has] been effectively obliterated.”

Notably, the high rates of channel migration in an erodible lithology are more likely to result in stranded, dateable geologic material than in a resistant lithology. Increased alluviation in a meandering reach increases the frequency of cutoffs, which strands material suitable for cosmogenic or luminescence dating. This biases these dating techniques towards periods of punctuated changes to local vertical incision rates. Given that alluviation impedes vertical incision rates and increases channel migration (leading to higher rates of channel abandonment that preserves datable material) lithologic context is important for interpreting vertical incision rates. This is particularly important over timescales that are similar to or less than the mean recurrence of punctuated events like meander cutoffs.

When locations with preserved fluvial material are sparse, interpretations of regional incision rates necessarily rely heavily on a handful of sampling localities, some of which may be particularly sensitive to local effects resulting from variations in lithology and the planform evolution of channels. This becomes increasingly true in estimates of short-term incision rate, where the temporal scale is insufficient to average over repeated punctuated episodes of incision resulting from autogenic processes and channel dynamics. When attempting to connect widely separated remnants to understand the current state of a fluvial system, interpretation becomes even more fraught. While far-felt effects of drainage integration will impact incision rates, this signal may be obfuscated by these shorter-timescale and local effects.

6 Conclusions

Fluvial deposits in the San Juan Canyon near Mexican Hat, Utah, dated using p-IR IRSL and cosmogenic radionuclide ¹⁰Be depth profile methods suggest that the San Juan River has experienced a short-term bedrock incision rate of 678–968 m Myr⁻¹ over the past ∼ 28–40 kyr. This rate is within the range of other short-term incision rates throughout the Colorado Plateau. While increases in incision rates elsewhere are attributed to the upstream propagation of knickpoints due to base-level fall as the Colorado River integrated through Grand Canyon and downstream basins, high short-term incision rates in the San Juan Canyon likely reflect, at least in part, cutoff bedrock meanders in the area surrounding the sample site. Evidence for additional incision and aggradation cycles following meander cutoffs in the San Juan Canyon may be correlated with the changing sediment supply during glaciation in the mountainous headwaters of the San Juan River. Preservation of fluvial sediment may be biased by planform rearrangements that also impact short-term vertical incision rates, requiring enhanced scrutiny of the geologic and geomorphic context of the few deposits remaining in a largely erosional landscape, such as the Colorado Plateau.