Geologists frequently debate the origin of iconic river canyons, as well as the extent to which they record climatic and tectonic signals. Fluvial and hillslope processes work in concert to control canyon evolution; rivers both set the boundary conditions for adjoining hillslopes and respond to delivery of hillslope-derived sediment. But, what happens when canyon walls deliver boulders that are too large for a river to carry? River canyons commonly host large blocks of rock derived from resistant hillslope lithologies. Blocks have recently been shown to control the shapes of hillslopes and channels by inhibiting sediment transport and bedrock erosion. Here, we developed the first process-based model for canyon evolution that incorporates the roles of blocks in both hillslope and channel processes. Our model reveals that two-way negative channel-hillslope feedbacks driven by block delivery to the river result in characteristic plan-view and cross-sectional river canyon forms. Internal negative feedbacks strongly reduce the rate at which erosional signals pass through landscapes, leading to persistent local unsteadiness even under steady tectonic and climatic forcing. Surprisingly, while the presence of blocks in the channel initially slows incision rates, the subsequent removal of blocks from the oversteepened channel substantially increases incision rates. This interplay between channel and hillslope dynamics results in highly variable long-term erosion rates. These autogenic channel-hillslope dynamics can mask external signals, such as changes in rock uplift rate, complicating the interpretation of landscape morphology and erosion histories.

River canyons, steep-walled valleys often developed in bedrock, evolve through a combination of deepening by river incision and widening by hillslope processes. Considerable effort has been expended on establishing the timing and mechanisms of canyon evolution (e.g., Cook et al., 2009; Schildgen et al., 2009; Flowers and Farley, 2012), with a focus on understanding landscape response to climatic and tectonic forcing. The traditional view of canyon erosion holds that river incision, driven by tectonics, climatic perturbations, or changes in substrate erodibility, lowers the canyon bottom. Adjacent hillslopes then respond to river incision by rockfall, landsliding, and/or diffusive sediment transport (e.g., Mudd and Furbish, 2007; Gallen et al., 2011). Under the assumption that sediment delivered to the channel is mobile, patterns and time scales of river incision control hillslope form and dominate canyon evolution. The majority of prior work on canyon development has embraced this view and drawn conclusions about the timing of canyon evolution under the assumption that hillslopes respond passively to river incision. However, canyon-confined rivers do not operate in isolation from their adjacent hillslopes (Egholm et al., 2013; Attal et al., 2015; Shobe et al., 2016, 2018; Bennett et al., 2016; Golly et al., 2017; DiBiase et al., 2018). Hillslopes make up the majority of canyon plan-view area and are often the primary source of sediment to the rivers. Steep canyon walls with substantial bare bedrock exposure and sufficient fracture density commonly release large pieces of rock (with diameters of several meters) into the channel (Howard and Selby, 1994; Glade et al., 2017; DiBiase et al., 2018; Glade and Anderson, 2018). Large grain delivery to rivers can inhibit incision over large spatial and temporal scales (Shobe et al., 2016, 2018), even damming rivers for short periods of time (Korup et al., 2006; Ouimet et al., 2007; Castleton et al., 2016). We propose that slowing or cessation of river incision must then influence the hillslopes by reducing the rate of hillslope steepening. Block delivery therefore acts as a negative feedback on both river and hillslope erosion. To constrain canyon evolution rates and process dynamics, it is critical to understand the interactions between canyon-confined rivers and their adjacent hillslopes.

We developed a numerical model that explicitly treats block dynamics both on hillslopes and in channels. We then examined the influence of these negative channel-hillslope feedbacks on river canyon evolution, with two guiding questions:

  • (1) Are negative channel-hillslope feedbacks necessary and sufficient to explain the cross section and planform shapes of natural canyons?

  • (2) How do these feedbacks affect long-term erosion dynamics in river canyons responding to base-level fall?

Block delivery is likely important in any block-producing landscape. However, as a simplified test case, we focused on river canyons in layered rock (Figs. 1 and 2) in which a resistant cap rock (e.g., sandstone) overlies softer rock (e.g., shale). In this simplified geologic setting, the cap rock acts as a line source of blocks, with block size dictated by cap-rock thickness and joint spacing. The softer, underlying layer produces soil but no blocks. Canyons developed in layered rock often exhibit key morphologic features, such as a characteristic bell-shaped planform during transient response to base-level fall (Figs. 1A and 1B), block-mantled channels, and steep hillslopes (Figs. 1C–1F). Here, we explore how block delivery feedbacks influence canyon form and evolution.

Several processes (e.g., landsliding, debris flows) may influence channel-hillslope coupling. We focused on the delivery of large blocks from hillslopes to channels (Fig. 2), which we hypothesized would control river canyon form. During the early stages of canyon evolution, channel incision causes failure of the cap rock (Ward et al., 2011), which delivers large blocks to the hillslopes and channel. The presence of blocks on the hillslopes inhibits soil erosion, stalling subsequent block release from the cap rock (Glade et al., 2017; Glade and Anderson, 2018). Blocks in the channel, if they are too large to be transported, reduce the river incision rate by armoring the bed and increasing hydraulic roughness (Shobe et al., 2016, 2018). Prolonged inhibition of river incision reduces the rate of hillslope steepening and hence the rate of block delivery to the channel. Thus, as long as hillslopes supply blocks to the channel (Fig. 2), erosion rates both in the channel and on the hillslopes are expected to be highly variable in time. Eventually, the hillslopes retreat far enough from the channel that blocks weather during hillslope transport to a size at which they no longer inhibit river incision (Fig. 2). From then on, the channel lowers at a rate that is unaffected by hillslope-derived blocks.

We tested the conceptual model shown in Figure 2 with a series of numerical experiments coupling models for channel (Shobe et al., 2016) and hillslope (Glade et al., 2017) evolution that incorporate the effects of large blocks of rock (see the GSA Data Repository1). The model domain, designed to represent the layered landscapes in Figures 1 and 2, consists of a horizontal resistant layer of rock overlying softer, more-erodible rock (the domain is 2 km wide × 1 km long, with 5 m resolution). A channel of uniform initial slope, forced with a constant base-level fall rate at its downstream end, incises the weaker, underlying rock. The channel permanently occupies the center of the model domain and has a constant width and discharge. The rest of the model domain operates under hillslope process laws (Glade et al., 2017; Glade and Anderson, 2018). We used this model to investigate erosion dynamics and time scales in a river canyon in which blocks released from the resistant cap rock cause interactions that govern both hillslope and channel evolution. In interpreting model behavior, we focused on the role of the two unique elements of our model: the geometry of horizontal layered rock, and the role of block feedbacks.

The two-dimensional (2-D) model couples a hillslope evolution model (Glade et al., 2017; Glade and Anderson, 2018) with a fluvial incision model (Shobe et al., 2016, 2018), both in the presence of blocks. The hillslope model uses 2-D depth-dependent linear soil diffusion (Johnstone and Hilley, 2015) and the exponential soil production function (Ahnert, 1976; Heimsath et al., 1997). Blocks, which are treated as discrete particles, experience a steady weathering rate and are released and transported down the steepest slope according to a local relief threshold. The fluvial incision model employs a modified shear-stress incision rule that accounts for the inhibition of erosion by blocks that both cover the channel bed and extract momentum from the flow. Block motion in the channel is calculated using a force balance (Larsen and Lamb, 2016), and blocks are abraded in proportion to the shear stress exerted on them. The two models were coupled using the Landlab modeling toolkit (Hobley et al., 2017).

Our model results in characteristic planform and cross-sectional features that agree with first-order field observations (Figs. 1 and 3). At early times, the model reproduces the bell-shaped planform (Figs. 1A, 1B, and 3) observed in the field. The planform shape results from two unique features of our model: (1) hillslope response times in horizontally layered rock, and (2) channel-hillslope block feedbacks. First, the curvature of the planform shape is dictated by both the hillslope and channel adjustment time scales (e.g., Mudd and Furbish, 2007). In layered rocks, the hillslope adjustment time scale is complicated by the geometry of the system, in which the hillcrest is pinned at the elevation of the resistant cap rock; this leads to a hillslope that grows in length and relief through time, and the scarp retreat rate therefore decreases through time. This geometry results in the gentle planform curvature observed in the control case without blocks (Fig. 3C). In the blocky case, in addition to geometric constraints, the blocks modify the hillslope and channel adjustment time scales. This results in a sharp kink in the planform (Fig. 3), which corresponds to the location in the channel that has roughly reached steady state (Figs. DR12 and DR13 in the Data Repository), below which the hillslopes have had sufficient time to lengthen. Our results thus indicate that the planform bell shape of canyons can arise as a result of canyon-escarpment retreat rate declining as hillslopes lengthen and hillslope response times increase, and that a sharply kinked plan-view shape can arise due to block dynamics. This situation is expected to be common in horizontally layered rocks, as well as any canyon system in which hillslope divides can migrate.

Incrementally increasing block size to explore the model parameter space, we found that diagnostic features grow more pronounced as the erosion-inhibiting effects of blocks increase (Fig. 3; Figs. DR14–DR28). Canyon cross sections at early times (Figs. 3A and 3B) show linear-to-concave-up canyon walls mantled with blocks. Channel longitudinal profiles exhibit knickpoints associated with channel steepening in response to base-level fall and block delivery. At later times, when the scarp has retreated far from the channel, blocks mantle only the upper portions of the hillslopes, leaving the lower portions to become convex-up because they are no longer influenced by the blocks (Figs. 3C and 3D). Only at later times do channel profiles become linear, as expected for a steady channel reach with uniform discharge (Fig. 3D).

Tracking the erosion rate through time at three points along the channel and hillslope (Fig. 4) reveals how negative channel-hillslope feedbacks control rates of canyon evolution for long periods of time. Near base level, the channel quickly equilibrates to the imposed base-level fall rate, with steady river erosion occurring until incision triggers block delivery from the hillslopes. At a given location within the canyon, block delivery to the channel, which causes increased bed armoring and hydraulic drag, results in a 100 k.y. period of highly variable, 1000-yr-averaged river erosion rates. This, in turn, causes unsteadiness in the hillslope boundary condition. The 1000-yr-averaged channel erosion rates vary between zero and four times the imposed base-level fall rate, indicating that block dynamics can force erosion rate variations comparable to those caused by several-fold changes in rock uplift rates. Only after 100 k.y. of erosion-rate oscillations driven by negative channel-hillslope feedbacks do the channel and hillslope experience relatively steady erosion. The high-frequency erosion rate fluctuations at later times are due to continued delivery of small blocks that do not perturb the 1000-yr-averaged erosion rates.

Farther up the channel, the response to base-level fall is delayed because block delivery from the hillslopes slows upstream propagation of the base-level signal. The channel erosion rate is minimal while the channel steepens to the point at which it can transport blocks, after which a >200 k.y. period occurs in which both the channel and hillslope erosion rates oscillate about the imposed base-level fall rate. Far upstream of base level, the oversteepened channel erodes more rapidly than the base-level fall rate for >100 k.y. (Fig. 4). This illustrates the potential for blocks to cause counterintuitive prolonged increases in erosion rate in addition to the expected decreases, masking the base-level forcing through internal feedbacks (e.g., Jerolmack and Paola, 2010).

Exploration of model parameter space revealed that greater block delivery increases both the amplitude of erosion-rate perturbations and the time scale required for the canyon to experience near-steady erosion rates (Fig. 3; see the Data Repository). The influence of blocks on river erosion will vary with the scale and hydrology of drainage basins, but our analysis suggests that it may be important in many canyon landscapes. This finding has implications for the applicability of common stream power–type approximations for understanding landscape response to perturbations. When substantial mass is delivered to the channel as blocks, feedbacks among block delivery, river incision, and canyon plan-view evolution significantly alter landscape dynamics in ways not captured by current models. While recent work has shown that stream power theory may be modified to include the effects of hillslope-derived blocks (Shobe et al., 2018), a stronger quantitative framework is needed to extract forcing signals from block-controlled landscapes.

This work presents the first model of canyon evolution capable of matching observations of plan-view form, cross-sectional shape, and the presence of large boulders. While we explored the example of a simple cap-rock canyon, our results apply to any landscape in which blocks are delivered from hillslopes to channels. We observed complex erosion rate dynamics due to negative, two-way feedbacks between incising channels and delivery of blocks from the adjacent hillslopes. We emphasize that model results illustrate a minimum effect of blocks in canyon evolution; many additional elements of reality would amplify the effects described here. For example, rockfall and landsliding would lead to more rapid block delivery. Landscapes where blocks are sourced along the entire hillslope (e.g., in well-fractured granitic rock; Granger et al., 2001) will also experience a greater influx of blocks to the channel than the line source explored in our model.

Our results imply that channel-driven models for canyon evolution may be overly simplistic, even when canyons evolve under a steady external forcing mechanism. Changes to channel incision rates caused by hillslope sediment delivery, in addition to changes in hillslope erosion rate due to unsteady channel incision, set both the shapes of canyons and the time scales of their evolution. Autogenic channel-hillslope feedbacks substantially modify base-level signals for hundreds of thousands of years, and increase the time required for canyon landscapes to equilibrate to an imposed base-level forcing. Bedrock canyon evolution can only be understood as the product of a coupled channel-hillslope system in which large blocks of rock play a critical role.

This work was supported by National Science Foundation grants EAR-1349390, EAR-1331828, EAR-1529284, EAR-1831623, and OAC-1450409. Shobe was supported by a National Defense Science and Engineering Graduate (NDSEG) Fellowship and a University of Colorado Chancellor’s Fellowship. Landlab is an element of the Community Surface Dynamics Modeling System (CSDMS; Boulder, Colorado, USA).We thank Katy Barnhart, Aaron Hurst, Will McDermid, Colin Meloy, and Matt Rossi for helpful discussions. Reviews by Sean Gallen, Sam Johnstone, one anonymous reviewer, and Editor James Schmitt greatly improved the paper.

1GSA Data Repository item 2019229, supplementary methods and results, is available online at, or on request from Data availability: The BlockLab model is archived on GitHub ( and the Community Surface Dynamics Modeling System (CSDMS) model repository ( Model output and figure plotting scripts are permanently available in a Figshare repository (
Gold Open Access: This paper is published under the terms of the CC-BY license.