We utilized field measurements of erosion rates and topographic analyses to constrain the timing and magnitude of landscape rejuvenation on the western flank of the Rocky Mountains in central Idaho, United States. Deeply incised canyons of the Clearwater, Salmon, and Snake Rivers dissect a broad region of roughly 8 × 104 km2. Along the Salmon River, an observable break in slope separates relict landscapes of low relief (<400 m valley depth) from high-relief landscapes (1200–1600 m valley depth) adjusting to base-level forcing. The 10Be cosmogenic radionuclide concentrations in river sediment record basinwide erosion rates that increase from 0.05 mm/yr ± 0.008 in the low-relief topography to 0.12 mm/yr ± 0.016 in the adjusting, high-relief landscapes over the last 103–104 yr and are consistent with longer-term estimates of erosion. Using the covariance of erosion rates and channel morphology, we calibrated a 1-dimensional river incision model to constrain the dynamics of incision along the Salmon River. More than 105 model runs explored uncertainty and assumptions and found that increased incision initiated roughly 9.5 ± 2 Ma and persists to the present. New constraints on the distribution of erosion processes at locations within a 400 km transect across central Idaho suggest northward surface tilting. In light of these data, we offer a new hypothesis that attributes late Miocene landscape rejuvenation of central Idaho to surface uplift driven by density changes in the mantle-lithosphere precipitated by the Yellowstone plume. We demonstrated the hypothesis through a simple model of flexure of an elastic plate subject to a buried buoyant load, and we found that density changes extending 200 km north of the Snake River Plain can reproduce the south-north distribution of uplift with reasonable values for elastic thickness and anomalous density.

The topographic evolution of a region is a consequence of the interplay between constructive processes that uplift the landscape relative to base level and destructive processes that denude the resulting elevation gradients. Geologic events that increase surface elevations drive a transient incisional response in the rivers flowing across the uplifted region, leaving deeply incised canyons in their wake. The transient landscape form offers an opportunity to calibrate and test incision rules and provide insight into the spatial and temporal dynamics of topographic evolution. The incision history throughout the Northern Rocky Mountains and the Inland Northwest, USA, provides a record of landscape evolution that offers insight into the recent tectonic and geologic events of the Northern Rocky Mountains that have shaped the topography (Beranek et al., 2006; Reidel and Tolan, 2013). For example, in central Idaho, the Salmon, Snake, and Clearwater Rivers have incised 1-km-deep canyons into a low-relief surface (Fig. 1) suggesting a recent rejuvenation of the landscape. Early investigations identified areas within the Inland Northwest as part of a deeply dissected plateau (inset Fig. 1B) where the elevated “peneplain” is the remnant of an old base level or erosion surface (Lindgren, 1904; Umpleby, 1912; Blackwelder, 1912), with some controversy over the age of incision as being either Eocene (Umpleby, 1912; although it is unclear why) or late Miocene to Pliocene (Blackwelder, 1912), based the unwarped condition of the peneplain relative to the pronounced deformation of Eocene sediments and lavas in adjacent regions. In this study, we investigated the enigmatic timing and cause of canyon formation in central Idaho through quantitative study of the transient topographic response to base-level change.

Significant advances made on the timing and processes of major Cenozoic geologic events in the region offer potential mechanisms that may have influenced landscape evolution and generated these high-elevation, low-relief surfaces. These mechanisms include the drainage capture of Lake Idaho by the proto–Snake River (Wood and Clemens, 2002), Columbia River Basalt eruptions (Reidel et al., 2013), uplift associated with the Yellowstone hotspot (Suppe et al., 1975; Smith, 1978; Pierce and Morgan, 1992; Waschbusch and McNutt, 1994), and surface uplift of the Wallowa Mountains (Fig. 1A) following convective removal of a dense lithospheric root (Hales et al., 2005). Studies on the erosion history of the region provide spatial and temporal information on the surface effects driven by these potential geologic processes. The long-term exhumation rate estimated from 40Ar/39Ar cooling ages is 0.05–0.06 mm/yr since 78 Ma near Elk City, Idaho (Fig. 1A; Lund et al., 1986), and fission-track dating measures 0.03–0.1 mm/yr since 50 Ma in the Boise Mountains (Sweetkind and Blackwell, 1989). Exhumation measured in the Boise Mountains (Sweetkind and Blackwell, 1989) agrees with (U-Th)/He dating in the Pioneer-Boulder Mountains (Fig. 1A; Vogl et al., 2014), and both suggest a 3-fold increase in exhumation beginning 11–8 Ma following a long period of little to no exhumation. We add to these data new constraints on Miocene incision at locations between 200 and 350 km north of the Boise Mountains and Pioneer-Boulder Mountains (Fig. 1A) based on quantitative analyses of transient topographic response to base-level change. In tributaries of the Salmon River (white box, Fig. 1A), a river incision model calibrated with field-measured erosion rates was used to estimate the timing and rate of rejuvenated incision. In tributaries of the Clearwater River (white polygons, Fig. 1A), topographic analyses compared base-level change relative to the Salmon River. We propose a hypothesis where a density change in the mantle-lithosphere beneath central Idaho initiated surface uplift that can explain the regional pattern of surface processes. Finally, we modeled the flexure of an elastic plate subject to a buoyant buried load to consider the lithospheric rigidity, load size, and density anomaly required to result in observed values.

Transient Incision of Central Idaho

The Salmon, Snake, and Clearwater basins drain much of the Northern Rockies of Idaho into the Columbia River (Fig. 1A). Canyons in these river systems often exceed 1 km depth and cut through hard igneous and metamorphic rocks. Throughout these river systems, particularly in the Salmon and Clearwater basins, the landscape is dissected by an obvious break in slope that separates elevated relict landscapes of low relief and adjusting landscapes of high relief. The headward advance of the landscape adjustment to the new base-level boundary condition is depicted in the contrast of surface slopes (Fig. 1B), which manifests as broad patches of low-slope terrain drained by steep river gorges dissecting an area of roughly 8 × 104 km2. Tributaries of the Salmon River (white box, Fig. 1A; Fig. 2) offer an appropriate landscape for study because they cut through relict and adjusting surfaces bearing quartz-rich Idaho Batholith rock, which is well suited for measuring cosmogenic 10Be–derived erosion rates. We can confidently attribute formation of knick zones in the Salmon River tributaries to a change in base level because their flow paths do not cross lithologic boundaries, which could introduce a slope bias related to resiliency of different rock types. Furthermore, the consistency in knick-zone elevations in Salmon River tributaries (Fig. 2) is evidence of uniform vertical celerity from a common base-level fall (Niemann et al., 2001; Gallen et al., 2013). While many researchers have recognized widespread landscape transience in central Idaho, the driving mechanism remains unclear.

New constraints on the spatial and temporal details of the transient topography will shed new light on the driving mechanism(s) of landscape rejuvenation in central Idaho. Catchment-averaged erosion rates were used to estimate the rates of Salmon River incision pre- and post-base-level change. Channel profile analyses identified the location of tributary knickpoints, indicated apparent differences in incision between the Salmon and Clearwater Rivers, and helped to quantify the morphological response to base-level forcing. The covariance of erosion rates and channel slope between the adjusting and relict landscapes allowed us to calibrate a 1-dimensional (1-D) numerical model that simulates the incision history of the Salmon River (e.g., Harkins et al., 2007; DiBiase et al., 2010; Legrain et al., 2015).

Cosmogenic Radionuclide–Derived Erosion Rates

Measurements of the 10Be concentrations in quartz-rich stream sediment have provided catchment-averaged erosion rates (Brown et al., 1995; Granger et al., 1996; Bierman and Steig, 1996) that are insensitive to short-term perturbations in erosion and capture the long-term signal pertinent to the geomorphology of a landscape (Kirchner et al., 2001; Bierman and Nichols, 2004; von Blanckenburg, 2005). Our sampled catchments range from 1 to 44 km2 and contain no variations in lithology or signs of landslides that could bias results due to uneven quartz distributions or sediment supply spiked with low 10Be concentration from deep-seated quartz, respectively.

Quartz separation was carried out according to the methods of Kohl and Nishiizumi, 1992), and 10Be isolation and concentration measurements were done at the PRIME Laboratory (Purdue University). Cosmogenic 10Be production rates for each basin were calculated pixel by pixel from a 10 m digital elevation model (DEM; USGS, 2014) while correcting for latitude, altitude (Dunai, 2000), and snow cover (methods described below). The resulting erosion rate was spatially and temporally averaged through the basin over a time scale equal to the penetration depth (0.6 m in rock) divided by the calibrated erosion rate, in this case, on the order of 4000–13,000 yr. As shown in other parts of the Idaho Batholith, cosmogenic radionuclide (CRN) concentrations yielded erosion rates similar to estimates over 106 yr time scales (Sweetkind and Blackwell, 1989; Kirchner et al., 2001). Thus, millennial-scale erosion rates are assumed to reflect long-term erosion rates in the region, though our modeling approach also assessed the implications of this assumption by testing the full range of error on erosion rates.

Snow Cover Corrections

Snow cover attenuates cosmic rays, thereby preventing 10Be accumulation in the underlying bedrock, which leads to an overestimation of erosion rate. Shielding corrections for snow cover have been approached in a variety of ways (Gosse et al., 1995; Ivy-Ochs et al., 1999; Benson et al., 2004; Schildgen et al., 2005). In the nearby Cascade Mountains of Oregon, corrections for snow cover reduced cosmogenic 3He production rates by 4%–19% (Licciardi et al., 1999). Our method corrects for snow shielding based on historic climate data of snow pack density with elevation and attempts to account for warming temperatures over the residence time of rock in the cosmic-ray attenuation depth, which extends into the Last Glacial Maximum (LGM).

The modern relationship of elevation to “fraction of shielded cosmic rays” (F) was determined from snow water equivalency data from four nearby SNOTEL sites (Banner Summit, Bear Basin, Secesh Summit, and Long Valley) in the Salmon Basin in Idaho (USDA Natural Resources Conservation Service, 2014). We estimated the elevation-F relationship at 13 ka to have the same slope as the modern relationship, but we adjusted to a larger degree of shielding to reflect colder temperatures at the end of the LGM. In most glaciated terrains, steep valley walls are subject to frequent snow avalanching, exposing rock for much of the year. Near the Salmon River study area, elongated U-shaped valleys indicative of glacial processes during the LGM extend down to ∼2300 m elevation. We made an assumption that one third of this landscape was exposed, while the other two thirds were covered by deep snow and ice, thereby pinning the 2300 m elevation to a snow shielded fraction of 2/3 during the LGM. We note this assumption is ad hoc, but it allows the shielding factor to adjust with warming climate. The correction for CRN production was forward modeled for each pixel by adjusting the elevation-F relationship from the paleovalue to modern value between 13 and 9 ka. Overall, these snow cover corrections reduced nuclide production rates by as much as 16% in relict basins and 5% in adjusting basins.

Channel Profile Analysis

We measured changes in channel slope for a given drainage area through the channel steepness index, ksn (Wobus et al., 2006). Channel slope (S) tends to decrease with contributing drainage area (A) through a power-law relationship (Morisawa, 1962; Hack, 1973; Flint, 1974; Howard and Kerby, 1983):
where ks is referred to as the steepness index, and θ is the concavity index. Regression of log-transformed slope and drainage area along the stream fits single values for ks and θ in relict and adjusting reaches. In order to avoid complications that arise from the sensitivity of ks to the best-fit concavity index, we used a fixed reference concavity (Wobus et al., 2006) equal to the average concavity of relict channels (0.465). Calculating ksn along the tributary profile is useful for quantifying the slope differences between relict and adjusting landscapes, as well as identifying the location of knick zones (Fig. 2).

Elevation data were extracted from a 10 m DEM (USGS, 2014) for each drainage basin delineated by the contributing drainage area upstream of each cosmogenic sample site. For each tributary basin, the elevation profile of the main channel draining >2 km2 was delineated with a flow direction algorithm in ARCGIS software and MATLAB scripts that followed the path of pixels downstream while recording elevation, streamwise distance, and contributing drainage area data. A 55 pixel (550 m) smoothing window was used to remove inherent noise in the river profile (Wobus et al., 2006). We report the average ksn values for relict and adjusting portions of study tributaries, as well as knickpoint elevations and percent of basin area that is relict landscape.

River Incision Model

The river incision model was used to determine the time required for the transient break in channel slope (knick zone) to migrate from the tributary mouths to their current locations, acknowledging uncertainties in model calibration (example in Fig. 3). We calibrated the model using channel profile analysis and erosion rate data to account for landscape response to base-level fall in both rate and topographic form. Model simulations focused on tributaries where the transient signal has not yet propagated all the way through the system. We did this for two reasons: (1) The ability to reproduce the location of the knick zone is important to distinguish between different modeled scenarios, and (2) by projecting the concavity and steepness of the relict portion of the river profile, the paleoriver profile is re-created (the blue line in Fig. 3; e.g., Schoenbohm et al., 2004) and provides the starting point for the river incision model. We modeled the transient evolution of three tributaries (Carey Creek, Fall Creek, and Pony Creek; Fig. 2) that contain a clear transient signal and show no evidence of modification by glaciation or significant landslides.

The change in elevation (z) over time (t) at each node along the river profile was determined by
where U (described in “Modeled Base-Level Scenarios” section below) is spatially uniform rock uplift, and erosion is set by the stream power erosion model (Whipple and Tucker, 1999; Snyder et al., 2000), in which K is rock erodibility, and drainage area (A) is a proxy for stream discharge. The m and n exponents depend on the hydraulic control on bedrock erosion (i.e., unit stream power or shear stress), basin hydrology, and erosion process. For equilibrium channels in which uplift and rock erodibility are uniform along the river profile, m and n strongly influence channel concavity (Moglen and Bras, 1995; Whipple and Tucker, 1999):
CRN-derived erosion rates and channel steepness data were used to calibrate erosion in the incision model by constraining the bedrock erodibility of Idaho Batholith granodiorite in the Salmon River landscape. Combining Equation 1, 2, and 3 reveals that K and n scale the relationship between erosion rate and channel slope,

Under the assumption that K is constant along the stream profile, which is a reasonable assumption for uniform lithology (Snyder et al., 2000), this relationship guides calibration of stream power in Equation 3. To constrain K, we used the standard deviation of CRN-derived relict basin erosion rates for E and the mean ksn for relict reaches in Equation 4. This approach resulted in a range of K values to be tested under different base-level scenarios that account for uncertainty in the way in which the CRN-derived erosion rates reflect longer-term rates.

Modeled Base-Level Scenarios

The transient incision of the river profile can be simulated in two ways. In the first case, the elevation of the mouth of the tributary is held constant, and the transient is due to uplift of the entire tributary profile with respect to that fixed location. In the second case, the elevation of the steady-state tributary profile is held constant while the elevation of the tributary mouth drops with respect to the geoid (i.e., no change in rock uplift, and the drop in base level is due to the downcutting of the Salmon River). Both scenarios produce the same river profile dynamics. The model alone, therefore, cannot distinguish between uplift and downcutting without additional constraints or paleo-elevations.

To constrain the timing of increased incision along the main-stem Salmon River, the model began from the steady-state condition represented by the reconstructed paleoprofile equilibrated to a rock-uplift rate equal to the relict erosion rate. Beginning from the equilibrium state, the downstream end node of the paleoriver profile was lowered to its modern position over a range of time scales (or model run times) up to 20 m.y. This simulated the drop in base level that drives increased incision. We used the model to explore different dynamics of base-level lowering within these time scales. For example, if base-level lowering is monotonic and persists to the present, then the duration of increased incision is spread out over 100% of the run time. To simulate more drastic, short-lived events, the total base-level fall may occur over as short as 10% of the run time, followed by a long period with no change in base level. In addition to these scenarios, which force a discrete perturbation of increased incision, we also tested scenarios of increasing and decreasing base-level lowering rate, which are analogous to a prolonged ramp-up or ramp-down of incision rate.

Model Assessment

The simulations produced a transient river profile for each calibrated parameter set, timing of increased incision, and variation of base-level signal, resulting in 5 × 104 model runs for each transient tributary. A comparison of the real transient profiles to simulation profiles with a Χ2 misfit analysis (Rosenbloom and Anderson, 1994) can determine the quality of each model fit:

The parameters N and υ are the number of nodes along the profile and free parameters in the model run, respectively, simi is the simulated elevation, and obsi is the observed elevation. The tolerance is the standard deviation of elevation differences between the smooth theoretical paleoprofile (Eq. 1) and the observed profile using only the relict domain of the river profile, which we measured to be 30 m for the three transient tributaries. High values of Χ2 suggest a poor fit well beyond the noise of the modern DEM, while a small value suggests a close reproduction of the modern profile. Simulations that resulted in Χ2 ≤ 2 were considered acceptable fits. Χ2 = 2 with a tolerance of 30 m means that the average difference between the simulated and observed elevation profiles is 60 m.

Timing of Salmon River Incision

Results from field work in the Salmon River show that the mean CRN-derived erosion rate (corrected for snow cover) is distinctly different between the steep, adjusted portion of the landscape (mean = 0.12 mm/yr, n = 3) and the relict portion (mean = 0.051 mm/yr, n = 4), rates that are consistent with previous studies over a range of time scales (Lund et al., 1986; Sweetkind and Blackwell, 1989; Kirchner et al., 2001; Meyer and Leidecker, 1999). Additionally, erosion rates in the few basins that contained both relict and adjusted basins are between these end member rates (Fig. 2; Table 1) and are in line with numerical models of transient effects on cosmogenic nuclides (Willenbring et al., 2013). This consistency gives us confidence in our use of cosmogenic nuclides to infer erosion rates. The mean erosion rate increases by a factor of ∼2.4 (from 0.05 mm/yr to 0.12 mm/yr) between the relict and adjusting landscapes, while channel steepness increases by a factor of ∼3.8 (from 80 m0.9 to 300 m0.9; Figs. 2 and 3; Table 1). These results show that erosion increases with channel steepness to the 2/3 power (n = 2/3, Eq. 1) and is consistent with a shear stress model of river incision (Howard and Kerby, 1983). The rock erodibility parameter (K) was calibrated with Equation 4 using mean ksn and a range of E set by the standard deviation of mean CRN-derived erosion rates of relict basins (Table 1). We calibrated the model based on mean values across the landscape rather than performing a calibration for each transient basin (Fig. 3). This prevented a single basin from skewing model results and weakened any bias that may be present in CRN concentrations in detrital samples or variations of erodibility in Idaho Batholith granodiorite.

The Χ2 misfit analyses tested various incision histories following the onset of increased incision. Results show that acceptable fits of transient profiles only occurred when a constant uplift rate (or constant downcutting) persisted over the entire run duration (Fig. 4). This implies that incision along the Salmon has remained high since initiation of gorge formation. Scenarios involving discrete perturbations or gradients in the base-level lowering rate produced radically different profiles that were consistently rejected in the Χ2 misfit analyses. Simulations where n ≥ 1 generated poor Χ2 fits, so we focus on n = 2/3 for all transient profile simulations.

The distribution of acceptable fits (Fig. 5) constrains the upper and lower bounds on the timing of increased incision in the transient tributaries. Various K (erodibility coefficient) values were calibrated for both the raw and snow-corrected CRN production rates. Each set of K values was tested over the range of run times to produce the two distributions of possible incision times shown in Figure 5. A comparison of the two distributions of raw and snow-corrected data shows that greater erodibility (raw CRN data) produced fewer acceptable fits, but no difference in the lower bound. Considering just the model runs from snow-corrected data, the acceptable fits among the different tributaries overlap at 9.5 ± 2 Ma for the onset of Salmon River canyon formation. The overlap of distributions using the raw CRN data suggest a timing of 8.2 ± 0.6 Ma, which is within the bounds of the snow-corrected timing. Conversely, if erosion rates over the last 10 m.y. were even slower than our snow-corrected CRN concentrations suggest, then results would show an older timing of increased incision.

Knickpoint Migration Along the Main-Stem Salmon

The timing of increased incision on the Salmon tributaries should be considered a minimum, because the incision and base-level signal probably initiated near the junction with the Clearwater River, which is where the degree of landscape dissection begins to rapidly decay to the north and west. The horizontal celerity of the incision signal derived from the stream power equation (Niemann et al., 2001),
shows that the upstream migration of knick zones depends on the uplift rate, rock erodibility, and drainage area (proxy for discharge), where D is horizontal distance, t is time, U1 is an initial uplift rate, and U2 is a subsequent uplift rate. Knickpoints initiating at the Snake-Clearwater junction would have propagated up to our study area in 2.4 m.y., using the mean Idaho Batholith erodibility in our Salmon tributaries (K = 2.7e–6 m0.4/yr). Knickpoints in basaltic basins downstream of our Salmon River study area have traveled a greater distance, albeit not exceptionally farther, suggesting the erosional resistance is similar or less than the granitic basins. The lack of knickpoints in tributaries located downstream of our study area that are composed of metasedimentary rocks suggests much less resistance, which is conceptually consistent with the highly fractured rocks in this zone. Thus, our timing of 9.5 ± 2 Ma should be considered a minimum, although it is possible that a knick zone originating near Lewiston, Idaho (Fig. 1A), formed as early as 13.9 Ma, if we consider the maximum time of knickpoint migration from the Snake-Salmon confluence (2.4 Ma) and the maximum timing of incision at the Salmon study area (11.5 Ma).

North-South Constraints on Late Miocene Landscape Evolution of Central Idaho

The results presented here provide the first quantitative constraints on Salmon River gorge formation. To develop a regional pattern of late Miocene landscape evolution in central Idaho, we included areas south and north of the Salmon River gorge. Constraints on the magnitude of incision at other locations will provide further insight on the possible mechanisms driving regional incision. In the Pioneer-Boulder Mountains, roughly 170 km south of our field area, (U-Th)/He dating suggests exhumation rates significantly increased (up to 4000 m total incision in valleys) beginning between 11 and 8 Ma, following a long period of little to no exhumation (Vogl et al., 2014). Compared to the Salmon River study area, the Pioneer-Boulder Mountains experienced greater than twice the magnitude of incision. North of the Salmon River study area, topographic analyses reveal a northward decay in the vertical migration of base level–driven transient knickpoints (Fig. 6). The three tributaries chosen for these analyses are predominately underlain by one rock type to ensure that changes in channel steepness are not mistaken for changes in rock strength. The Lochsa River shows strong evidence of a transient knickpoint by a spike in steepness (Fig. 6A). Farther north, in the Weitas River of the Clearwater and the Little North Fork of the Clearwater, the signal is weaker, perhaps due to the erodibility of the underlying rock (Figs. 6B and 6C). The anomalous deviations in channel steepness indices show a northward decline in knickpoint elevations, suggesting a lesser magnitude of base-level fall, or rather, north-facing tilt of the regional landscape. Roughly 100 km farther north, cosmogenic 10Be concentrations in the Fernan Creek watershed (Fig. 1A) constrain erosion rates to 0.049 ± 0.017 mm/yr (Parker, 2016), which are roughly equivalent to the erosion rates we measured in the Salmon River relict landscape and the long-term Eocene exhumation rates of the Idaho Batholith (Lund et al., 1986; Sweetkind and Blackwell, 1989). In addition to the decay of knickpoint elevations and slow millennial erosion rates in Fernan Creek, topographic relief decreases north of the Clearwater River (Fig. 1A), suggesting that this area experienced little or no late Miocene increased incision, while areas to the south did. The additional topographic analyses combined with CRN-derived erosion rates in Fernan Creek, Salmon River erosion modeling, and exhumation of the Pioneer-Boulder Mountains suggest that an area of the landscape roughly the extent of the Idaho Batholith (Fig. 1C) was tilted to the north.

Recent Drivers of Landscape Evolution in the Inland Northwest

Our timing of Salmon River incision and the rapid exhumation of the Boise Mountains (Sweetkind and Blackwell, 189) and Pioneer-Boulder Mountains (Vogl et al., 2014) suggest a late Miocene age of onset of landscape rejuvenation. Figure 7 places these results temporally with other geologic events of the late Miocene, including the integration of the Snake River Plain into the Columbia Basin, the evolution of the Yellowstone plume path, and the emplacement of Columbia River Basalt. These geologic events provide the motivation for hypotheses on the origin of landscape rejuvenation in central Idaho. Next, we discuss how processes related to these events may be responsible for the landscape rejuvenation.

The spatial and temporal extents of lacustrine deposits in southwestern to south-central Idaho suggest the prior existence of a large lake that covered much of the modern western Snake River Plain (Wood and Clemens, 2002). The capture of Lake Idaho and the Snake River Plain by the proto–Snake River occurred sometime between 6.4–1.7 Ma (Wood and Clemens, 2002), thereby tripling the drainage area of the Snake River to 179,000 km2 downstream of the spillover point near Weiser, Idaho (Fig. 1A). This integration would have increased the stream power at the confluence with the Salmon and Clearwater and potentially forced increased incision along these rivers. The emplacement of Columbia River Basalts from 17.5 to 6 Ma (Reidel et al., 2013) also has the potential to have influenced the topographic evolution of the region. Basalt flows emerging from dike swarms along the western Idaho border (Fig. 1C), near Lewiston, Idaho, intermittently filled in low topography and valleys of the Columbia River drainage. This could potentially have driven a base-level rise, temporarily slowing erosion and creating low-relief relict topography upstream of basalt emplacement. After the cessation of eruptions, a return to rapid incision would have been driven by the increase in surface elevation relative to base level. The Yellowstone plume has the ability to impact landscape evolution across central Idaho in different ways. Several researchers have suggested that the high topography of the Yellowstone hotspot is supported by convective thermal buoyancy of a mantle plume (Suppe et al., 1975; Smith, 1978; Pierce and Morgan, 1992; Waschbusch and McNutt, 1994) at ∼150 km wavelength (Wegmann et al., 2007). Additionally, the Yellowstone plume could cause a change in lithospheric density that would drive uplift due to isostasy. Interaction between a spreading plume head and the edge of older cold lithosphere would lead to decoupling of the crust from the mantle-lithosphere and delamination at the lithosphere boundary (Burov et al., 2007), where less dense material flowing in to replace the delaminated lithosphere would result in a density change. Alternatively, tomography shows that the modern Yellowstone plume stem dips 65° WNW beneath central Idaho (Yuan and Dueker, 2005). If plume geometry or plume flux differed in the past where the plume stem had more direct interaction with the mantle-lithosphere beneath an aging Idaho Batholith, it may have facilitated a change in density through thermal alteration.

The timing of both the Snake River Plain integrating into the Columbia Basin and the emplacement of the Columbia River Basalt does not correlate well with the timing of incision along the Salmon River or with increased exhumation of the Pioneer-Boulder Mountains (Fig. 7). Notwithstanding potential variability in our predicted timing of events, particularly the lag between peak basalt emplacement and incisional response, the constrained timing of Salmon River incision alone does not answer the origin of landscape rejuvenation without uncertainty. However, the apparent disparity in incision magnitudes offers an additional criterion. If incision were related to drainage capture or basalt emplacement, we would expect increased incision to occur with a similar timing and magnitude in all of the subbasins. The difference in knickpoint elevations among the Salmon and Clearwater basins and the different magnitudes of exhumation among Fernan Creek, the Salmon River, and the Pioneer-Boulder Mountains suggest otherwise. Landscape transience does not appear to have been initiated by either Snake River Plain integration or Columbia River Basalt emplacement, but Snake River Plain integration may have played a part in continuing increased incision to the present. Diachronous base-level changes might be distinguishable as separate anomalies in transient river elevation profiles. While we did not observe such evidence in our studied tributaries, we cannot rule out that perhaps increased stream power from the subsequent integration with the Snake River Plain drove some amount of incision post–6.4 Ma.

Our preferred explanation for the regional pattern of incision is related to density changes in the mantle-lithosphere caused by the Yellowstone plume interacting with an intracontinental ocean-continent boundary, because it provides the most agreeable timing and explanation for northward down-tilting of the landscape. The western edge of the Idaho Batholith (Fig. 1C) roughly marks the north-striking subvertical tectonic boundary where island-arc terranes to the west were sutured to the North American craton to the east during the mid-Cretaceous (Lund and Snee, 1988). Due to plate motions, the Yellowstone plume passed from thin oceanic lithosphere into thicker cratonic lithosphere during the late Miocene. The onset of Salmon River canyon formation (9.5 ± 2 Ma) and rapid exhumation of the Pioneer-Boulder Mountains (11–8 Ma) coincide with the ca. 12–9 Ma intersection of the Yellowstone plume stem with the intracontinental ocean-continent tectonic boundary, the NW to SE rhyolitic volcanism and formation of the western Snake River Plain graben along the edge of the Idaho Batholith (Fig. 1C), and the beginning of hotspot-related volcanism in southern Idaho (Fig. 1A; Wood and Clemens, 2002; Shervais et al., 2002). In our hypothesis, the Yellowstone plume initiated a density change in the mantle-lithosphere beneath the southern end of the Idaho Batholith, either through delamination or thermal processes. Northward tilting of the central Idaho surface is driven by a source of uplift generated at the southern end of the Idaho Batholith that decays to the north. The plume path model of Shervais and Hanan (2008) provides a reasonable premise. They propose that a spreading plume head beneath thin, accreted oceanic lithosphere impinged against the thicker North American craton ca. 16.6 Ma. Subsequently, the plume stem was deflected to the south against the edge of the craton, culminating in the Breneau-Jarbidge eruptions (Fig. 1A) at ca. 12 Ma. In this scenario, the spreading plume head could have initiated delamination (e.g., Burov et al., 2007), with the deflecting plume stem establishing a channel that continues to remove material. Seismic tomographs of the modern Yellowstone plume and of central Idaho support this hypothesis. Tomography of the modern Yellowstone plume shows a low-velocity structure adjacent to a high-velocity structure interpreted as possible lithospheric downwelling next to the plume (Yuan and Dueker, 2005). Tomography of central Idaho along the 44th parallel shows significantly slower mantle seismic velocities, suggesting lithosphere may have been removed, and warmer mantle has filled in this region (Shen et al., 2013). However, if the geometry of the plume stem differed in the past due to pulses of high flux, the plume may have thermally reduced the density of the base of the lithosphere farther north of the hotspot track. Seismic tomography data in this region are subject to a range of interpretations (Hales et al., 2005; Yuan and Dueker, 2005; Schmandt and Humphreys, 2010; Shen et al., 2013), but none of them discounts a change in the lithospheric density structure between south and central Idaho.

There are other mechanisms associated with the plume that may have driven uplift in some part. While dynamic elevation of the modern Yellowstone hotspot is supported at a wavelength of ∼150 km (Wegmann et al., 2007), which is narrower than the >350 km north-south swath of landscape transience observed west of Yellowstone, we acknowledge the possibility that dynamic topography could extend farther north if the lithospheric rigidity is stiffer than assumed. Also, pressure-driven asthenospheric flow caused by the spreading plume head could drive dynamic topography (e.g., Colli et al., 2014). However, both of these effects should decay as time passes, and our results suggest persistent incision without decay. Ultimately, we prefer a hypothesis of a density change in the mantle-lithosphere beneath central Idaho caused by interaction with the plume stem for three reasons. (1) Uplift from the convective buoyancy of the plume should decay over time as the hotspot migrated east, and our results do not support this. (2) We have no reason to think the lithosphere beneath our study area is stiff enough to support dynamic topography from convective buoyancy at a 350 km wavelength. (3) Tomography of the modern Yellowstone plume (Yuan and Dueker, 2005) and a transect along the 44th parallel in central Idaho (Shen et al., 2013) encourage the interpretation of removal of dense mantle-lithosphere or anomalously buoyant density. In this framework, the southern end of the Idaho Batholith uplifts due to isostasy, while the rigidity of the lithosphere supports uplift as far north as the Clearwater Basin, and long-term uplift is sustained by either long-lasting density change or less dense material flowing in to replace the root. It is beyond the limits of this study to further differentiate the landscape response to either delamination or thermal change, and therefore, we describe the preferred driving mechanism as a nonspecific density change in the mantle-lithosphere.

Miocene Surface Uplift and Crustal Flexure Driven by a Buoyant Load

In this section, we explore the ability for a buoyant density anomaly to cause a northward-decreasing trend in surface uplift. We used a 1-D analytical model of elastic flexure of a thin plate subject to a load of uniform height and given width (Watts, 2001) to simulate the topographic response to a density anomaly. The model introduces a buoyant rectangular load to the mantle with a fixed load width of 200 km, roughly the distance of the plume stem deflection against cratonic lithosphere (Fig. 1A) proposed by Shervais and Hanan (2008). Model runs differed by varying the effective elastic thickness (Te) of the deforming plate between two end members, 10 and 40 km (typical of thinned continental lithosphere with high heat flux). Solutions to the thickness and density anomaly necessary to deform the surface to within 200 m of three uplift constraints are reported in Table 2. We used constraints of surface uplift and rock uplift (uplift without the effects of erosion) since 10 Ma to account for uncertainty in the rate of landscape response, particularly at the Pioneer-Boulder Mountains. For example, the landscape could have responded rapidly to root detachment but maintained extended uplift through feedback to isostatic uplift or continued dynamic support from mantle flow around the detached root. The three constraints of surface uplift involve: (1) the lack of evidence for increased incision north of the Clearwater River, (2) the elevation of relict surfaces in the Salmon River, and (3) the magnitude of exhumation in the Pioneer-Boulder Mountains (Fig. 8). The constraints of rock uplift do not account for elevation decrease from surface erosion (∼500 m over 10 m.y. in the relict landscape along the Salmon River). The first constraint is the proposed hinge zone where there is a transition to topography that appears to be unaffected by base-level fall. River valleys of this landscape share a common elevation with the Clearwater and Salmon valleys (Fig. 8), which provides a reasonable baseline elevation from which to begin flexure simulations.

Accurate fits of surface uplift were possible for both end-member values of Te, with distinct differences in the shape of the deformed surface. Based upon these results, a stiffer lithosphere is more capable of reproducing the observed landscape transience. Better fits with thin Te values may require wedge-shaped rather than rectangular loads, or remnant hotpot topography (supported at 150 km wavelengths) to create the elevation rise between the Salmon study area and the Pioneer-Boulder Mountains, which is not accounted for in this model. Accurate fits of rock uplift are only possible with Te = 40 km. In all three solutions, the resulting density anomaly ranged between 200 and 300 kg/m3 less than mantle density, which is typical of the difference between eclogite and peridotite. Regardless of how realistic the pure rock-uplift scenario is, it accounts for the maximum landscape response to an upper-mantle density anomaly, regardless of how long lasting it may have been. Overall, this analysis demonstrates that a change in density over a scale of 15–40 km thick and 200 km north of the Snake River Plain has the potential to drive an isostatic response and landscape rejuvenation as far north as the Clearwater Basin.

Assessment of Assumptions and Uncertainties

As with any modeling study, we made assumptions that could impact our estimates of canyon incision if the model does not capture the essence of river dynamics appropriately. Here, we discuss the implications of errors in these assumptions. First, our analysis included only a few basins per topography type. More basins for each topography type (e.g., relict and adjusted) would improve the robustness with which we constrained erosion rates within the relict and adjusting landscapes. Also, the snow shielding corrections for CRN production rates required an estimate of its variation with elevation, and it could therefore impart a false difference in erosion rate between relict and adjusting landscapes if this value is incorrect. It is not known how much snow covered the landscape over the last 13,000 yr, limiting our ability to formally assess this uncertainty. Local factors also could impact snow cover, such as wind direction, hillslope aspect, and changes in mean annual precipitation. Nevertheless, snow cover impacts cosmic-ray attenuation, and the effect was likely greater during the LGM. We acknowledge that if the impact was greater than our snow cover model predicts, then we are overestimating the erosion rate. Slower CRN erosion rates would ultimately push the onset of base-level fall farther in the past due to less erosional efficiency in our river incision model. Also, our assumption that the Idaho Batholith has a consistent erodibility throughout does not consider the effects that variation in fracture density may have on erosional efficiency. While, the vertical migration rate of knickpoints is not sensitive to substrate erodibility (Niemann et al., 2001; Gallen et al., 2013), it can be sensitive to erosion process, which could vary depending on bedrock fracture spacing. This may be responsible for relatively small difference in observed knickpoint elevations. However, the time it would take to incise the various observed knickpoint elevations (1 m.y. to incise 120 m, assuming a 0.12 mm/yr rate of erosion) is within error of our analysis. Although we cannot formally assess all the uncertainty in this analysis, our modeling approach allowed for variation in erosion rates and erodibility to attempt to account for the range of possible outcomes. This led to the spread of acceptable X2 fits seen in Figure 5. We also note the ability of our incision model calibration to accurately reproduce knickpoint profiles. If the erosional efficiency parameters (K and n from Eq. 4) were much different than what was calibrated, the shape of the simulated profiles would be fundamentally different from the observed profiles. Finally, to place the results in a regional context in support of quantitative evidence of northward tilting, we acknowledge that this assessment is limited to four locations over a 350 km distance. However, the northward-decreasing topography and relief (Figs. 1A and 8) further support this conclusion. An exhaustive search of base level–driven transient knickpoint elevations throughout the Inland Northwest could provide more insight into the spatial extent of surface tilting. While other geologic processes may have contributed to sustaining incision to the present, the degree to which they have contributed remains poorly understood. For instance, we acknowledge the possibility that the original uplift event may have ceased, and an isostatic response to erosion and transport of sediment out of the landscape may be maintaining uplift to the present.

This study demonstrated that mantle processes most likely are driving the observed landscape transience in central Idaho, but robust constraints on the mechanism causing landscape transience are lacking. A simple model of crustal flexure shows that a 6%–9% density reduction along the Yellowstone plume path is capable of driving landscape rejuvenation over an area of roughly 80,000 km2 in central Idaho. Exploration with different models of dynamic topography resulting from plume and lithosphere interaction at an inactive ocean-continent plate boundary could allow more precise predictions of the temporal and spatial evolution of transient uplift. These predictions could be tested with surface process analyses over a broader extent than shown here and perhaps provide greater insight on the coupling between mantle and surface processes.

The landscape of central Idaho, USA, is deeply dissected by steep canyons that were set off by rapid incision. River profile analyses confirm that transient knickpoints separate relict and adjusting landscapes along Salmon River tributaries. Erosion rates determined from 10Be concentrations show an increase in erosion rate by a factor of ∼2.4 between the relict and adjusting landscapes. River steepness increases by ∼3.8 between the relict and adjusting landscapes. Our river incision model calibrated with field-measured erosion rates suggests that incision rates along the Salmon River increased ∼3-fold beginning 9.5 ± 2 Ma and persisted at that rate to the present. The synthesis of Fernan Creek erosion rates (Parker, 2016), knickpoint elevations in the Clearwater Basin, Salmon River incision rates, and exhumation rates of the Pioneer-Boulder Mountains from Vogl et al. (2014) suggests a regional base-level change beginning 13–8 Ma and a broad region (roughly the size of the exposed Idaho Batholith) of northward-tilting uplift north of the Yellowstone hotspot track. These results support a scenario in which most of the observed landscape rejuvenation in central Idaho is driven by mantle processes associated with plume and lithosphere interaction at an intracontinental boundary between oceanic and cratonic lithosphere. Crustal flexure modeling shows that a buoyant density anomaly beneath the Idaho Batholith reaching 200 km north of the Snake River Plain could drive tilting of the landscape and ultimately transient incision >350 km north of the hotspot track. Our hypothesis of a density change in the mantle-lithosphere beneath the Idaho Batholith provides a framework and an explanation for the enigmatic high topography of central Idaho and long-wavelength increase in erosion across the region. However, the degree to which other ancillary processes may have contributed to landscape evolution, and the mechanism of density change remain poorly understood.

The University of Idaho Office of Research and Economic Development supported this study. Cosmogenic nuclide concentrations were measured at PRIME Laboratory at Purdue University. Yanites was supported by National Science Foundation grant EAR- 1727139. This work benefited from thoughtful reviews by Reed Lewis, Karl Wegmann, Anke Friedrich, and anonymous reviewers; however, any errors or inaccuracies presented here are solely the responsibility of the authors.

Gold Open Access: This paper is published under the terms of the CC-BY-NC license.