Abstract

Landscapes developed over heterogeneous stratigraphy exhibit a spectrum of landforms from dramatic cliffbands to hogbacks, depending on the dip and spacing of the layers. In deeply incised landscapes, a single cliffband may consist of multiple resistant layers, whereas similar stratigraphy elsewhere is separated by strike valleys into individual cuesta benches or hogbacks. This paper presents a geometric analysis, informed by a numerical landscape model, to explain the conditions for development of a strike valley floored by erodible rocks. The results define a threshold incision rate below which strike valleys are more likely to form; this threshold incision rate is proportional to the stratigraphic spacing of cliff-forming layers and a trigonometric function of dip angle. The analysis also yields a time scale for the adjustment of structural landforms to changes in regional incision rate, which is a function of dip angle and the coupling between cliff retreat rate and escarpment height. In example landscapes of the Colorado Plateau, this time scale is likely much longer than that of documented variations of incision rates due to late Quaternary climate and land-use changes. The transitional state of escarpments in layered rock may therefore contain information about regional downcutting rates over time scales different from those recorded by the fluvial network. The utility of such features will require better understanding of the coupling between incision of a foot slope and the retreat rate of the cliff above in different kinds of rocks.

INTRODUCTION

Motivation

In landscapes with uniform rocks and a steady rate of rock uplift or base-level fall, topography evolves to a quasi-steady form that equalizes erosion rates spatially across the landscape (Gilbert, 1877; Tucker and Bras, 1998). Even where rocks differ spatially in erodibility, a quasi-steady topography can form under conditions of uniform base-level fall, where landforms on more resistant rocks evolve to have steeper slopes to attain the same erosion rates as those on weaker rocks (Gilbert, 1877; Hack, 1960). This relationship implies that, given appropriate information about material and erosion processes, uplift rates and gradients thereof may be deduced from topography (Wobus et al., 2006).

In contrast, recent work has highlighted how the heterogeneity in rock strength can lead to the indefinite persistence of transient topography and to varying local erosion rates in landscapes experiencing steady base-level fall. For example, where fluvial networks cross dipping stratigraphy, lithologic knickpoints that form on harder units will migrate upstream, which may be updip or downdip. This process results in a local, internal base-level control that behaves differently from the external base level (Cook et al., 2009; Berlin and Anderson, 2009; Forte et al., 2016; Perne et al., 2017). Interpretation of erosion rates, or tectonics, in these landscapes based on landform is therefore complicated by multiple sources of variability: changes to dip, stratigraphic thickness, erodibility, and internal sources of transience in the downcutting rates (Darling and Whipple, 2015).

Nonfluvial landforms that commonly develop on variably erodible stratigraphy include cliffbands, cuestas, homoclinal ridges, hogbacks, and flatirons. These “structurally controlled” landforms are typically classified based on somewhat arbitrary ranges of dip angles of the resistant rock capping the feature (Twidale and Campbell, 2005). There is not, however, any unifying theory that relates stratigraphic and lithologic properties under various rates of base-level fall to specific structurally controlled landforms. It is my purpose here to explain a key morphological distinction along this spectrum of landforms, which is the presence or absence of a strike valley floored by more erodible rocks. The result lends insight into the steady and transient states of landscapes developed on heterogeneous stratigraphy.

Colorado Plateau

Structurally controlled landscapes occur on sedimentary rocks worldwide. Prominent examples include the Valley and Ridge region of eastern North America; the Atlas Mountains of northern Africa; the Jura Mountains of Europe; the Arabian Desert and Negev Desert; the Flinders and Macdonnell Ranges of Australia; and the Qinling and Dalou Mountains of China. Some of the most iconic landscapes of this type are found in the Colorado Plateau region of the United States.

In all of these areas, the landforms are clearly controlled by the stratigraphic spacing of resistant units. For example, resistant layers spaced closely together commonly form a compound cliffband, wherein the foot slope below an upper cliff drains directly over the lower cliff (Fig. 1). Where resistant layers are spaced more widely, broad benches may form on the weaker rock between them, separating the cliffs and routing drainage parallel to strike. Dip is another primary control: Steeply dipping layers form homoclinal ridges or hogbacks separated by strike valleys, whereas more horizontal layers of the same rock may form compound cliffs. The Laramide monoclines of the Colorado Plateau provide classic examples of this stratigraphic and dip control on the topography of escarpments (Fig. 1).

Because the field-measurable terms of interlayer spacing and structural dip contribute to the variability of landforms, the relationship between landform and base-level fall is not as straightforward as in uniform substrates (Forte et al., 2016). In the example of the Colorado Plateau, the cliff and cuesta landforms developed during a time of significant documented spatial and temporal variability in downcutting rates (Darling et al., 2012; Jochems and Pederson, 2015). The most prominent transience in the fluvial networks is related to deep, late Cenozoic entrenchment of the trunk streams due to drainage integration through the Grand Canyon. This transience resulted in several sets of major (∼200-m-relief) knick zones along the river system, some of which may reflect a combination of base-level and lithological effects (Cook et al., 2009; Darling et al., 2012; Bursztyn et al., 2015) or differential tectonic uplift rates (e.g., Crow et al., 2014). Superimposed on these major, long-term incision events, there are 102 to 104 yr cut-and-fill cycles that are typically of 1 to 10 m in magnitude (Harvey et al., 2011; Pederson et al., 2013b; Sheehan and Ward, 2018). These episodes of incision of different magnitudes and durations imply that variability in the surrounding landforms may be due to their current states of response to the variability in downcutting, with their specific morphology conditioned by the bedrock template.

Reference Case: Coal Cliffs Cuesta, Utah, USA

As noted already, steeply dipping cap rocks that form hogbacks are commonly separated by strike valleys, but it is less common for strike valleys to form among subhorizontal rocks. An example of this less common case is the Coal Cliffs cuesta (Fig. 2), located on the western flank of the San Rafael Swell in central Utah. The San Rafael Swell is a Laramide-aged upwarp, ∼60 km across, with a steeply dipping (20°–80°) eastern flank and a shallowly dipping (2°–10°) western flank. The stratigraphy of the western flank is composed of Triassic through Cretaceous terrestrial and marine sedimentary rocks (Cashion, 1973). The sandstones tend to act as cliff-forming units, and the weak shale and mudstone interbeds form foot slopes and strike valleys.

The Coal Cliffs escarpment is developed in the Cretaceous Mancos Formation. The 10–20-m-thick Ferron Sandstone Member of the Mancos Formation acts as the cap rock for the Coal Cliffs (Fig. 2B). Below, the ∼100-m-thick Tununk Shale Member forms a bench or strike valley, ranging from a few hundreds of meters to just over 1 km wide atop the Dakota Formation sandstone, which itself forms the cap of a lower cliffband. Above the Ferron member, there is the Blue Gate Shale Member of the Mancos Formation, which is several hundred meters thick and forms a kilometers-wide strike valley between the Coal Cliffs and the base of the Mesaverde Group sandstones at the edge of the Wasatch Plateau to the west. These layers dip consistently 3°–4° to the west.

Climate-driven incision cycles have resulted in short-term erosion rates of 1–3 mm/yr in the vicinity of the Coal Cliffs (Sheehan and Ward, 2018). The escarpment is crosscut by multiple tributaries to the San Rafael River and Dirty Devil River, which join the main-stem Green River and Colorado River, respectively, ∼100 km downstream. Near its junction with the San Rafael River, the Green River has been incising at 0.45 mm/yr over the last 100 k.y. (Pederson et al., 2013a). Upstream, in Desolation Canyon, the Green River has incised more slowly, at ∼0.04 mm/yr, over the last 1.8 m.y. (Darling et al., 2012). The faster rates downstream possibly represent a transient acceleration in downcutting rates within the Green-Colorado network (Cook et al., 2009; Roberts et al., 2012; Pederson et al., 2013a), but it is not clear how far into the San Rafael Swell such a signal has propagated. Considering these factors, rates of longer-term fluvial incision in the Coal Cliffs are likely on the order of 0.1–0.5 mm/yr.

The Coal Cliffs represent a clear end-member expression of low dip angles, thin resistant and thick weak layers, and well-developed strike valleys that drain to escarpment-crossing transverse streams, some of which have a catchment area above the Ferron sandstone and some of which do not.

Pour Points and Benches

The upper segments of any escarpment that exposes two resistant units must ultimately drain across the lower resistant unit. On a classic hogback, or a cuesta like the Coal Cliffs, the foot slope of the upper cliffband drains to a bench underlain by the more erodible rock. Channels on the bench drain parallel to strike, until reaching a larger transverse stream or break in the lower cap rock.

Where a bench or strike valley is present, the width over which it is floored by more erodible rocks is correlated with the distance from the upper cliffline to the nearest pour point through the lower cap rock (Fig. 2C). Unsurprisingly, the farther the pour point lies outboard of the upper cliffline, the wider is the area in which the more erodible rocks are exposed in the strike valley. Where the pour point distance becomes less than the length of the steep foot slope of the upper cliffband, the strike valley disappears (Fig. 2A), and the cliffs transform to a compound form. Drainage along strike becomes limited, and small channels of the foot slope of the upper cliff drain directly over the edge of the next cliff below. In well-developed compound cliffs, the pour points through the lower cliff are positioned within the bounds of the upper cliffline, as in the classical “rule of V’s.” Similarly, where transverse streams are superimposed on a compound cliff form, they cross the lower cap rock at pour points that are inboard of the upper rim of the escarpment. In turn, the pour point positions are related to the relative entrenchment of the local drainage network.

I posit that this geometry of pour point location dictates the form of the escarpment by controlling the route by which base-level changes are communicated to the upper portion of the escarpment. Next, I derive the conditions under which the lower cap rock pour point on an escarpment-crossing stream would be aligned along strike with the forefront of the upper cliffband. The goal is to relate the stratigraphic spacing of resistant layers, the structural dip, and the regional downcutting rate to the presence or absence of a strike valley, under conditions of steady base-level lowering and steady topographic relief. I then discuss some salient aspects of the transient evolution of these systems as base-level lowering rates change through time.

ANALYTICAL MODEL

Model Geometry

Following the Coal Cliffs example, consider an escarpment exposing two resistant units separated by a weaker unit of thickness Hi, and dipping upstream relative to larger transverse streams that cross the escarpment approximately perpendicular to strike (Fig. 3). This analysis seeks the conditions of stratigraphy, dip angle α, and base-level lowering rate (ζ) under which, upon reaching steady topographic relief, the transverse stream’s pour point through the lower cap rock is aligned along strike with the top of the upper cliff (Fig. 3D). Geometrically, this occurs when the steady-state upper cliff height Hs is equal to Hv, defined as the apparent thickness of the weak rock interlayer in the vertical direction (i.e., Hv = Hi/cos α). Hs < Hv is the bench-forming case (Figs. 3A and 3C), and Hs > Hv favors drainage of the upper cliff over the lower cliff and a compound cliff (Figs. 3B and 3E).

For simplicity, and based on the Coal Cliffs example, the transverse streams are assumed to have significant upstream area (>>Abench). They therefore maintain low slopes (here taken to be zero slope), even where crossing the resistant layers, and respond quasi-instantaneously to base-level change. This assumption can be relaxed if necessary with a geometric correction to account for the gradient of the transverse stream.

The rate of change of the upper cliff’s height H can be written: 
graphic
where H is the instantaneous height of the upper cliff, ζ is the base-level lowering rate, assumed to be applied synchronously along the entire transverse stream, R is the horizontal retreat rate of the upper escarpment, and α is the dip angle. Assuming steady relief has been attained, dH/dt = 0, and 
graphic
Unless the retreat rate R is a function of escarpment height, there can be no steady-relief condition. This analysis is predicated on the assumption that such a condition is attainable and therefore that such a relationship exists. There is not any universal theory for the form of R(…), but there is some empirical and model support for taller escarpments experiencing larger retreat rates (cf. Howard and Selby, 2009; Haviv et al., 2010), particularly where the retreat is driven by undermining of the resistant layer by downcutting of the foot slope (e.g., Koons, 1955; Howard, 1995). An increase in retreat rate with height is furthermore consistent with the scaling of hillslope relief with uplift/incision rate in linear and nonlinear-diffusive settings (Roering et al., 2001). The simplest such form of R(H) is linear with proportionality coefficient c1 [1/T]: 
graphic

The c1 term here encapsulates the mechanics of the backwasting process, scaling the creation of relief by the channel network of the foot slope to the size and time distribution of rockfall failure events.

For an escarpment as illustrated in Figure 3, H cannot physically exceed Hv, but we may still consider that there is a theoretical steady height for the upper cliff. Inserting Equation 3 into Equation 2, and solving for height, which per the assumptions of Equation 2 must be the steady-state height Hs, we obtain: 
graphic
Where Hs = Hv, the pour point through the lower cliff is aligned along strike with the salient edge of the upper cliff (Fig. 3D). The system is at the point of transition from bench to compound cliff, as streams flowing along the bench must cross the lower cap rock before reaching the more competent transverse stream. 
graphic
Defining the incision rate that results in this configuration as ζc, and rearranging Equation 5, we obtain: 
graphic

Therefore, for a given layer spacing Hi and dip angle α, whether a weak-floored bench or strike valley develops depends on a threshold incision rate that increases linearly with interlayer spacing and is proportional to tan α / cos α (Fig. 4A).

Dimensionless Height

Nondimensionalization can reveal the key relationships between terms that govern the overall behavior of a system. Equation 4 can be divided by Hv as a reference length, yielding a dimensionless number H*: 
graphic
H* < 1 means that the equilibrium cliff height is less than the vertical apparent thickness of the interbed, and therefore a cuesta bench may form. H* > 1 represents conditions likely to form a compound cliff. The ratio H* is linearly dependent on the collection of terms: 
graphic
which is the dimensionless ratio between the base-level lowering rate and the cliff retreat rate at a cliff height equal to the interlayer spacing. The interlayer spacing here takes the role of a scale height for the behavior of the escarpment against the regional downcutting rate.

Analytical Results

Equations 6 and 7 posit a nonlinear relationship for resistant layer spacing, dip angle, and incision rate that dictates whether the landscape will evolve toward a compound cliffband or cuesta/hogback form as it approaches equilibrium relief (Fig. 4). At low dip angles, the apparent vertical thickness of weak units between the cap rocks is only slightly exaggerated, so the product of this vertical thickness and the retreat constant c1 is outstripped by even fairly slow incision rates. Thus, for a given value of c1, low dip angles require low erosion rates for benches to form, and even widely spaced resistant layers will tend to form a compound cliff. Conversely, vertical dips result in an infinite apparent thickness of the weak interbed, and so the cap-rock layers must be separated by a strike valley under any incision rate (in that case, parallel hogbacks).

NUMERICAL MODEL

Model Description

The simplicity of the previous solution relies on the assumption that R = c1H. To determine the functional form of this relationship for a set of common erosion rules, I applied a two-dimensional (2-D) finite-difference landscape evolution model derived from that of Ward et al. (2011). The model simulates stochastic failure of a hard-capped cliff within a detachment-limited fluvial landscape carved into more erodible rock. Hillslopes evolve according to a linear-diffusion scheme (Tucker and Bras, 1998), and bedrock incision is based on unit stream power (Whipple and Tucker, 1999). Stochastic rockfall from the cliffband is simulated by instantaneously reduction of the height of a selected cliff-edge cap-rock grid cell to that of its lowest downslope neighbor (following Howard, 1995). The equivalent volume of the cell is added to a model layer that tracks rockfall debris thickness and is distributed across the landscape below this cell by relaxing the surface of the debris layer to the angle of repose. Debris is continuously redistributed to the angle of repose as the surrounding landscape evolves, or until it erodes away. Further details of the model implementation are included in the GSA Data Repository material.1

A single cap rock was configured among uniformly weaker rocks within a long model domain oriented parallel to the dip direction. Model stratigraphy dipped 3.4° in the +y direction, based on the example of the Coal Cliffs (Fig. 2). The resistant layer was 20 m thick and 200× more resistant to erosion than the surrounding rock. The cliff was initially set to fail by the rockfall process when the upper edge of the cliff reached a slope exceeding 60° (in the model, corresponding to a particular elevation difference between cells as a proxy for the height of the vertical portion of the cliff; the effects of this threshold are discussed below). The angle of repose of rockfall debris was 25°. To isolate the intrinsic relationship between height and backwasting rate, rockfall debris was assigned a high erodibility, such that its effect on cliff retreat was minor and identical for all experiments. Harder debris would create an armoring effect on the foot slope, which is a potentially important factor in the height–retreat rate relationship (Ward et al., 2011; Glade et al., 2017). The strike-parallel model boundaries were open, with periodic boundaries on the transverse edges. I imposed a range of rock uplift rates (10−5 to 10−4 m/yr) and allowed the cap rock to be exhumed, develop an escarpment, and retreat until retreat rates became quasi-steady.

Model Results: Retreat Rate Proportional to Height

Once the cap rock was exhumed, the escarpment height increased over the duration of each model run until exhumation of the downdip edge of the cap rock resulted in rapid erosion of the entire landform via dissection of the backscarp. During the increasing-relief phase, the height of the escarpment edge increased most rapidly at first and then slowed until approaching a steady height (Fig. 5). More rapid uplift rates resulted in taller escarpments. Exhumation of the backscarp occurred more rapidly with higher uplift rates, requiring longer model runs and longer domains to capture the growth phase of increasing relief and approach to steady relief that are relevant for these experiments.

The instantaneous rate of retreat during the growth phase was found by extracting the average y position and height of the clifftop through each model run (Fig. 5B). Retreat rates were noisy, due to the stochastic, cellular nature of the rockfall process, but systematically increased through model time as the escarpment grew taller. Dividing the instantaneous retreat rate by the instantaneous escarpment height revealed a constant proportionality, as was asserted heuristically in Equation 3. This proportionality constant took on a value (c1 ≈ 3.5 × 10−6 yr−1) that did not vary between runs with different uplift rates. Sensitivity testing showed that this constant is not strongly affected by dip angle, cap-rock thickness, or the fluvial erodibility of either rock type. Instead, it is an approximately linear function of the specified gradient of initiation for the rockfall process in the model (see Data Repository item). In other words, it dictates how steep the upper extent of the drainage network must become to cause rockfall, effectively coupling the stream network and the retreat rate via the cliff height.

Transient Approach to Steady Height

The approach to a steady escarpment height following either emergence of the cliff or a change in base-level lowering rate, given Equations 1 and 3, has an exponential secular-equilibrium form: 
graphic
where H0 is the initial height of the escarpment. This can be seen in the model results between the times of cliff emergence and backscarp dissection (Fig. 5B). The e-folding time scale defined by τ = (c1 tan α)−1 dictates the time needed for the cliff height and retreat rate to respond to a new uplift/incision rate. In the model configuration used here, this time scale is 4.8 m.y., so nearly full adjustment requires on the order of 10 m.y. In rocks with steeper dips, this response time is reduced as 1/tan α.

Following Equation 4, at steady state, τ = Hs/ζ. This relationship allows an estimate in the field of τ, and therefore also c1, in the rare circumstances where steady-state relief can be assumed. For example, if the ∼100 m height of the Ferron escarpment is assumed to be in equilibrium with a long-term average downcutting rate of 0.1 mm/yr at the Coal Cliffs site, then c1 = 1.4 × 10−5 yr−1, and τ = ∼1 m.y. Even if the escarpment there is equilibrated to a much faster long-term mean downcutting rate of 1 mm/yr, then c1 = 1.4 × 10−4 yr−1, and τ = ∼100 k.y. Although this escarpment has not been shown to be in relief steady state, these example calculations suggest that the landform would be expected to equilibrate over a time scale longer than that of late Pleistocene–Holocene environmental changes in downcutting rate.

MODEL RUNS WITH MULTIPLE CAP ROCKS AND TRANSVERSE STREAMS

Model Setup

The same numerical model was used to examine the transient response of simulated landscapes with more than one cap rock that are near the threshold erosion rate expected from the analytical treatment above. A second, 10-m-thick cap rock was added, 30 m below the 20-m-thick cap rock of the single-layer runs. The dip-parallel boundaries were changed to open, fixed-elevation boundaries and are equivalent to the transverse streams.

Here, I report the results of simulations given two uplift histories that illustrate the key behavior. In both cases, there was an equilibration period in which uplift rate was steady at 0.1 mm/yr for 1.8 million model years, until a near-steady topography was reached over a homogeneous substrate. This was followed by 0.2 m.y. of zero uplift to allow the initial relief to decay somewhat before the cliff was exhumed, preventing the existing topography from draining over the upper cliff as it emerged.

Following the 2 m.y. equilibration period, the first scenario maintained a constant 0.008 mm/yr until 10 m.y. and then increased to 0.03 mm/yr for the remainder of the model run (Fig. 6A). These values lie near to but on either side of the H* = 1 line for this stratigraphic geometry. The second uplift scenario consisted of a series of brief pulses of 0.1 mm/yr uplift spaced irregularly (Fig. 6E). The pulses resulted in long-term uplift rates that averaged 0.03 mm/yr for the next 4.4 m.y. of the model run and 0.008 mm/yr for the final 3.8 m.y. This time series more closely idealizes the episodic, transient forcing documented on the Colorado Plateau, and it allows examination of the system response over different time scales.

Model Results: Transient Behavior with Two Cap Rocks

Simple Step Change in Uplift Rates

In this scenario, during the slow uplift (H* < 1) phase, a bench of constant width formed between the two layers. When the uplift rate increased (H* > 1), there was a period of transient adjustment that lasted ∼5 million model years, as first the upper cliff increased in height until the lower cap rock was exposed beneath it, and then the bench of weak rock disappeared. As the height of the upper cliff above the lower became limited by the vertical projection of the interlayer thickness, its backwearing rate also became limited. This allowed the lower cliff to “catch up,” combining the two cliffs into a compound cliff (Figs. 6B–6D).

Episodic Uplift

Throughout the model run, the cliffs began to separate to form a bench during the intervals between uplift pulses and then converged again when uplift resumed. During the portion of the model run with the faster average uplift rate (H* > 1), the cliffband generally behaved as a compound cliff. In this phase, drainage from the upper foot slope over the lower cliff enhanced lower cliff retreat, creating alcoves and bringing the two cliffs closer during uplift pulses.

During the phase with the slower average uplift rate (H* < 1), the pulses of uplift were too infrequent to maintain convergence of the two cliffs, and a persistent cuesta bench formed. Erosion transients reached the upper cliff via small scarps that propagated from the transverse drainages along strike on the erodible rock of the bench. (Equivalent features are observed in the field; see Sheehan and Ward, 2018.) The upper cliff at the midpoint between the two transverse streams was the last area to respond to each uplift pulse, and salients formed in the cliffband there (Figs. 6F–6H).

DISCUSSION

Expected Steady Form and Landform Adjustment Time Scales

The analytical and numerical model results presented here describe how a spectrum of landforms can result from a common base-level forcing due to the nonlinear interplay among dip, stratigraphic bed thickness, and incision rate (cf. the “multiple modes of adjustment” described by Phillips, 2003). The analysis presented in Figure 4 shows what escarpment form should result on approach to a state of steady relief: Closely spaced resistant layers, shallow dip angles, and rapid downcutting rates each favor a compound cliff morphology. If background downcutting rates remain slower than a threshold defined by Equation 6, cuestas or hogbacks separated by strike valleys will form so long as there are streams cutting the escarpment that connect the strike valley to base level.

Equilibration time scales to reach these states following a change in incision rate are geologically significant and may be long enough that the spatially variable dip of the beds, further changes to downcutting rate, or dissection of the escarpment via the backscarp would be encountered before steady state is reached (Howard and Selby, 2009). In the numerical model runs, the e-folding time scale for adjustment was 4.8 m.y.; this was the outcome of an arbitrary parameterization. For highly erodible escarpments in the Colorado Plateau, such as the Coal Cliffs, conservative estimates of the response time scale are ∼100 k.y., but it is more likely on the order of 1 m.y. Different landforms may therefore be in different states of adjustment to the same changes in incision rate. For example, because of the dip dependence of the response time scale, portions of the steep limb of the San Rafael Swell dipping at 25° would be expected to attain steady relief ∼7× faster than the shallowly dipping limb at 4° dip (Eq. 9).

The response of the escarpment form locally affects fluvial incision rates. Because separation between cuesta cliffbands expands and contracts during transients in the incision rate, the stream networks experience a fluctuating configuration of upstream area (see Movies SM1 and SM2 in the Data Repository). These locally driven fluctuations of incision rate would be further modulated by the dip control on local incision rates described by Berlin and Anderson (2009) and Forte et al. (2016).

Nonetheless, the transient numerical model results indicate that the dominant form of the escarpment reflects the downcutting rate averaged over the equilibration time scale, including any superimposed shorter periods of more rapid or less rapid incision. In the context of the Colorado Plateau, this means that escarpment landforms with longer equilibration time scales likely reflect the long-term regional rock uplift rates, or long-term downcutting trends due to river integration over the late Cenozoic, while faster-responding landforms may change substantially in response to hydrologically driven cycles of downcutting related to glacial-interglacial climate variability.

Retreat Rate as a Function of Height: Role of Lithology and Debris Armoring

These relationships between downcutting rate and landform are primarily due to large-scale geometry, and they are surprisingly independent of erosion process and the detailed lithology of the rocks involved. As long as cliff retreat rate is a monotonic function of escarpment height, the systematics should be grossly similar to those described here. In the model framework presented here, the coupling is assumed to be linearly proportional to the constant c1, which is validated by the numerical model (given its own set of erosion rules and corresponding assumptions). A nonlinear, but monotonic relationship (e.g., RHp, where p is a real number) would result in a different steady height for a given erosion rate, and therefore a different threshold incision rate for form transition, but the same general framework relating dip, spacing, and incision rate would persist. If retreat rate is a superlinear function of height (p > 1), the landscape adjustment time scale would be shorter than in the linear case; in the case of a sublinear relationship (p < 1), it would be longer. This emphasizes the importance of understanding the coupling between backwearing rate and escarpment height in a given setting.

For this analysis, I explicitly assumed that the retreat of the cap-rock edge at the top of a cliff is controlled by erosion of the weaker rocks beneath, and not the properties of the cap rock itself. The presence of lithologically different cap rocks with different backwasting mechanics might invalidate these assumptions or alter the form of the relationship between height and backwearing rate. For example, on the Colorado Plateau, the backwearing process of cliffs in massive Jurassic eolianites without exposed basal shales may be dominated by rock mass strength and fracture patterns, and not by basal undermining (e.g., Kimber et al., 2002; Howard and Selby, 2009).

An additional lithologic effect on the height–retreat rate coupling is due to the presence of debris shed from the hard cap rock onto the foot slope. Depending on the longevity, size distribution, and thickness of this talus, it may create a significant negative feedback on downcutting of the foot slope, and therefore on the retreat rate, due to channel-armoring (Bryan, 1940; Ward et al., 2011) and modification of hillslope sediment fluxes (Glade et al., 2017, 2019). As the talus may vary significantly depending on the character of the cap rock from which it is shed, and its thickness and residence time may be independently modulated by climate and vegetation (Bull, 1991; Gutiérrez et al., 1998), it represents a significant factor in the relationship between height and retreat rate.

Role of Transverse Streams

The role of transverse streams in both the analytical treatment and numerical model is to provide a path of communication of base-level change into the region between two exposed cap-rock layers. If the bench or strike valley cannot respond quickly to a base-level change, then the response of the upper cliffband is delayed as well, favoring loss of the bench to the retreating lower escarpment. It is important to the formation of a strike valley that the transverse stream is erosionally efficient relative to the other streams that drain the upper foot slopes. However, it is not strictly required that the transverse streams have significant area upstream of the upper cliff (as in the Coal Cliffs), or even that they cross the upper escarpment. If a large portion of the upper foot slope is routed through a single notch in the lower cliff, a similar relationship between the position of this pour point and the escarpment form ensues (see Movie SM3 in the Data Repository for an example).

The distance between transverse streams or other pour points along an escarpment affects the longitudinal variability in the escarpment form. Base-level changes reach the upper cliff quickly near the transverse stream, but it takes time for these changes to propagate along the bench or strike valley (Sheehan and Ward, 2018). If the transverse streams are widely spaced, the response of the upper cliff near the midpoint between them will be significantly delayed; this introduces another time scale of adjustment into the problem that is separate from (and additive to) the retreat response time scale (Eq. 9). The along-strike response to incision therefore exhibits a complex behavior, the exploration of which is beyond the scope of this paper.

CONCLUSIONS

Here, I have quantitatively explored how the field-measurable terms of interlayer spacing and structural dip contribute to the variability of structural landforms under different downcutting rates. For a given dip angle and local stream downcutting rate, wide layer spacings tend to result in strike valleys or cuesta benches once steady topographic relief is reached. For a given layer spacing and incision rate, more steeply dipping rocks are more likely to form strike valleys, resulting in the classic hogback landscape. Faster incision rates promote the loss of the strike valleys and the development of compound cliffs, where multiple cap rocks are exposed in a single escarpment face. Under conditions of a linear relationship between escarpment height and backwearing rate, the threshold erosion rate below which a strike valley will emerge increases linearly with interlayer spacing and is proportional to the tangent over the cosine of the layer dip α. Therefore, the threshold incision rate decreases rapidly with increasing dip over low dip angles and becomes less sensitive to dip at high dip angles.

The approach to steady-state relief has an e-folding time scale that scales with 1/tan α, and this is possibly several million years given reasonable values for the height–retreat rate coupling coefficient c1. This response time scale can be found as the ratio of escarpment height to downcutting rate in a setting known to be in relief steady state. However, given the range of dips and rock types found on the Colorado Plateau, it is expected that various cliff and cuesta landforms are in different states of response to documented episodes of downcutting.

The most difficult term to constrain for a particular cliffband is the coupling between escarpment height and backwearing rate. This coupling may vary widely depending on the mechanics of erosion and may not everywhere (or anywhere) be linear in form. This strongly motivates a physical approach, via understanding undermining mechanics and debris dynamics, as well as further chronometric work to determine rates of cliff retreat in response to documented base-level changes.

ACKNOWLEDGMENTS

This work was supported by U.S. Army Research Office grant 67195-EV-YIP. Chris Sheehan, Missy Eppes, Brian Yanites, and Bob Anderson are gratefully acknowledged for insightful discussions. Critical comments by two anonymous reviewers greatly improved the presentation. Model code is available via the Community Surface Modeling Dynamics System (https://csdms.colorado.edu).

1GSA Data Repository Item 2019305, Numerical model description and Movies SM1–SM3, is available at http://www.geosociety.org/datarepository/2019 or on request from editing@geosociety.org.
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.