Steady state , continuity , and the erosion of layered rocks

Considerations of landscape steady state have substantially informed our understanding of the relationships between landscapes, tectonics, climate, and lithology. Topographic steady state, where topography is fixed in time, is a particularly important tool in the interpretation of landscape features, such as bedrock channel profiles, within a context of uplift patterns and rock strength. However, topographic steady state cannot strictly be attained in a landscape with layered rocks with non-vertical contacts. Using a combination of analytical solutions, stream erosion simulations, and full landscape evolution simulations, we 5 show that an assumption of channel continuity, where channel retreat rates in the direction parallel to a contact are equal above and below the contact, provides a more general description of steady state landscapes in layered rocks. Topographic steady state is a special case of the steady state derived from continuity. Contrary to prior work, continuity predicts that channels will be steeper in weaker rocks in the case of subhorizontal rock layers when the stream power erosion exponent n < 1. For subhorizontal layered rocks with different erodibilities, continuity also predicts larger slope contrasts than would be predicted by 10 topographic steady state. Continuity steady state is a type of flux steady state, where uplift is balanced on average by erosion. If uplift rate is constant, continuity steady state is perturbed near base level. These perturbations decay over a length scale that is an inverse function of the contrast in erodibility, such that erodibility contrasts of more than approximately a factor of three lead to rapid decay. Though examples explored here utilze the stream power erosion law, continuity steady state provides a general mathematical tool that can be used to explore the development of landscapes in layered rocks using any erosion model. 15


Introduction
The formation of landscapes is driven by tectonics and climate, and often profoundly influenced by lithology, the substrate on which tectonic and climate forces act to sculpt Earth's surface.Much of our interpretation of landscapes, and their relationship to climatic and tectonic forces, employs concepts of landscape equilibrium, or steady state.Though there are a variety of types of landscape steady state (Willett and Brandon, 2002), topographic steady state, in which topography is constant over time, is perhaps most often used in the interpretation of landscapes.In particular, topographic steady state is used ubiquitously within studies of bedrock channel morphology.
Bedrock channels are of particular geomorphic interest because they span most of the topographic relief of mountainous terrains (Whipple and Tucker, 1999;Whipple, 2004), providing the pathways through which eroded material is routed to lowlands and a primary means by which the landscape is dissected and eroded.Therefore, bedrock channels exert important controls on the relief of mountain ranges and set the pace at which mountainous landscapes respond to changes in climate or tectonic forcing.Research on bedrock channels has driven new understanding concerning the coupling between mountain building, climate, and erosion (Molnar and England, 1990;Anderson, 1994;Whipple et al., 1999;Willett, 1999).
The elevation profiles of bedrock channels provide a primary means for analyzing landscapes for evidence of transience, contrasts in rates of tectonic uplift, or the influence of climate (Stock and Montgomery, 1999;Snyder et al., 2000;Lavé and Avouac, 2001;Kirby and Whipple, 2001;Lague, 2003;Duvall et al., 2004;Wobus et al., 2006;Crosby and Whipple, 2006;Bishop and Goldrick, 2010;DiBiase et al., 2010;Whittaker and Boulton, 2012;Schildgen et al., 2012;Allen et al., 2013;Prince and Spotila, 2013).Within this analysis, erosion rates are typically assumed to scale as power law relations with drainage area and slope, as given by the stream power erosion law (Howard and Kerby, 1983;Whipple and Tucker, 1999), where E is erosion rate, K is erodibility, A is upstream drainage area, S is channel slope, and m and n are constant exponents.
While the stream power erosion law has known limitations (Lague, 2014), it remains the most frequently used tool for channel profile analysis.Under steady climatic and tectonic forcing, channels are typically assumed to adjust toward topographic steady state (Hack, 1960;Howard, 1965;Willett and Brandon, 2002;Yanites and Tucker, 2010;Willett et al., 2014), where uplift and erosion are balanced and topography is constant with time.This assumption, in conjunction with the stream power erosion law, enables interpretation and comparison of stream profiles to identify spatial contrasts in uplift rates or transient responses to changes in tectonic or climatic forcing.
Topographic steady state has also been used to explain channel responce to substrate resistance, generally leading to a conclusion that channels are steeper within more erosion resistant bedrock and less steep within more erodible rocks (Hack, 1957;Moglen and Bras, 1995;Pazzaglia et al., 1998;Duvall et al., 2004).However, this result depends on an implicit assumption of vertical contacts between strata as in Fig. 1A.Strictly speaking, topographic equilibrium does not exist when channels incise layered rocks with different erodibilities and non-vertical contacts (Howard, 1988;Forte et al., 2016).In the case of nonvertical contacts, the contact positions shift horizontally as the channel incises, resulting in topographic changes as shown in Fig. 1B,C.Studies of bedrock channel morphology have primarily focused on regions with active uplift, where rock layers are often deformed and tilted from horizontal.However, a substantial percentage of Earth's surface contains subhorizontal strata.
Many of these settings also contain bedrock channels, with examples including the Colorado Plateau, the Ozark Plateaus, and the Cumberland and Allegheny Plateaus.In such settings, intuition developed from assumptions of topographic equilibrium does not necessarily apply.Here we explore the limitations of topographic equilibrium in landscapes formed in layered rocks and develop a generalized mathematical framework to predict the steady state form of a landscape developed within layered rocks that have an arbitrary dip.

Landscape continuity and steady state
To understand erosion and morphology near a lithologic contact, we begin with an ansatz that the land surface will strive to maintain topographic continuity at the contact.This is equivalent to an assumption that erosion rates above and below the contact will be equal in the direction parallel to the contact plane.This assumption is also identical to topographic equilibrium in the case of vertical contacts.The ansatz of continuity is supported by physical reasoning.If the upper layer erodes slower than the lower layer, this produces a steep, or possibly overhanging, land surface at the contact (Fig. 2A).This steepening will lead to faster erosion in the upper layer and drive the system towards continuity (Fig. 2B).Similarly, if the upper layer erodes faster, this produces a flattened zone near the contact (Fig. 2B) that will reduce erosion rates in the upper layer and also push the system toward continuity.Therefore, the same types of negative feedback mechanisms that drive landscapes to topographic steady state (Willett and Brandon, 2002) will also drive landforms near a contact into a state that maintains continuity.We refer to this type of equilibrium as continuity steady state.There are cases in natural systems where continuity is not maintained at all times.For example, caprock waterfalls are similar to the case in Fig. 2A.However, even in this case the discontinuity cannot grow indefinitely.If the waterfall reaches a steady size then the system has once again obtained a state where continuity is maintained in a neighborhood near the contact.Numerical landscape evolution models do not typically allow cases such as Fig. 2A-B.Therefore, numerical models are likely to maintain continuity even more rigidly than natural landscapes.
Using the constraint of continuity, one can write a very general relationship between surface erosion rates and slopes at a contact between two rock types, where E i and S i are erosion rates and slopes, respectively, and the index refers to rock types 1 and 2. S c is the slope of the rock contact.This relationship results from an assumption of equal retreat rate at the contact within both rock layers in a direction parallel to the rock contact plane, illustrated in Fig. 1C and Fig. A1.A similar relationship is used by Imaizumi et al. (2015) to examine the parallel retreat of rock slopes.If we consider the more specific case of stream power erosion through a pair of weak and strong rocks, this leads to where K w is the erodibility of the weaker rock, K s is the erodibility of the stronger rock, S w = tan θ w and S s = tan θ s are the slopes of the channel bed in each rock type, and the contact slope is S c = tan φ (derivation in Appendix A).
When the contact slope is much larger than the channel slopes (S c S w , S s ), the right hand side of Eq. ( 3) is approximately one, and vertical erosion rates in both lithologies are roughly equal.Rock uplift can thus be balanced by erosion in both segments, and the standard relationship between channel steepness in the two lithologies, normally derived from topographic equilibrium, is recovered, with If the contact slope is in this steep limit, but not vertical, the contact position and topography will gradually shift horizontally with erosion and vertically with uplift, while still obeying this relation derived from topographic equilibrium.
For the case where channel slopes are much greater than the slope of the contact (S w , S s S c ), as would be common in subhorizontal rock layers, Eq. (3) simplifies to (5) In this case, continuity results in roughly the same rate of horizontal retreat in both rocks at the contact, Fig. 1B.This contrasts with the standard assumption of equal rates of vertical erosion, and leads to unexpected behavior.Specifically, if n < 1 then higher slopes are predicted in weaker rocks, which is in strong contrast to intuition developed from the perspective of topographic equilibrium.This results because the rate of horizontal retreat within a given rock layer (dx/dt ∝ K i S n−1 i ) is a decreasing function of slope if n < 1.Since horizontal retreat rate is an increasing function of erodibility, continuity requires that increases in erodibility are offset by increases in slope.For subhorizontal contacts with n > 1, higher slopes are once again predicted in stronger rocks.
The slope ratio (S w /S s ) is depicted for the vertical and horizontal limits in Fig. 3A as a function of n for an erodibility contrast of K w = 2K s .In general, contrasts in the slopes within the two strata in the subhorizontal case as in Eq. ( 5) are larger than would be predicted using the standard formulation for vertical contacts in Eq. ( 4).In subhorizontal rocks, channel slopes may become sufficiently high or low to be driven to values outside the range of validity of the stream power erosion law, particularly for cases of n ≈ 1.Perhaps the most common value of n used within landscape evolution models is n = 1, therefore it is also notable that the continuity relation for subhorizontal strata contains a singularity at n = 1 (Fig. 3).The slope ratio (S w /S s ) diverges for n → 1 − and approaches zero for n → 1 + .This suggests strong dependence of channel behavior on n when n is close to 1.The singularity results because for n = 1 the horizontal retreat rate is independent of slope and solely a function of erodibility.Therefore the channel cannot maintain continuity by adjusting steepness.

Continuity steady state and stream profiles
The channel continuity relations above apply to channels within the neighborhood of a contact.Though there are clear longterm constraints on the relative retreat rates of any two contacts, these are not sufficient to determine an entire profile.Here, we examine whether the continuity relation applies along an entire profile using simulations of channel and landscape evolution in horizontally layered rock.

Methods for one-dimensional simulations and analysis
We solve the stream power erosion law using a first order explicit upwind finite difference method.This method is conditionally stable, and the timestep was adjusted to produce a stable Courant-Friedrich-Lax number of CFL = 0.9.The explicit upwind scheme has commonly been used for prior studies, though it is also known to produce smoothing of channel profiles near knickpoints (Campforts and Govers, 2015).Experimentation with resolution suggests that the conclusions presented here are not dependent on the numerics.The simulations employed 2000 spatial nodes.For simplicity, basin area was held fixed over where k a = 6.69 m 0.33 and h = 1.67.These parameter values are representative of natural drainage networks (Hack, 1957;Whipple and Tucker, 1999).Simulations were run with n = 2/3, n = 1, and n = 3/2.The value of m in the stream power erosion law was adjusted according to the choice of n to assure that the concavity m/n = 0.5, which is typical of natural channels (Snyder et al., 2000).Both high uplift (2.5 mm y −1 ) and low uplift (0.25 mm y −1 ) cases were run.Simulation parameters were adjusted to provide a similar number of rock contacts in each case.For the high uplift cases, rock layers were 50 m thick, whereas for the low uplift cases rock layers were 10 m thick.Longitudinal distances were also adjusted with the high uplift cases simulating 50 km long profiles and the low uplift cases simulating 200 km long profiles.Specific parameter values are provided in Table 1.
Simulation results are most easily visualized in χ space (Perron and Royden, 2013;Royden and Perron, 2013), where the horizontal coordinate x is replaced with a transformed coordinate χ: One advantage of this transformation is that the effect of basin area is removed such that equilibrium channels that evolve according to the stream power erosion law appear as straight lines in this transform space.The relation predicted by Eq. ( 5) is invariant under the transformation to χ space, and therefore the relation also holds if slope is replaced with slope in χ-elevation space.Throughout this work, we use a value of A 0 = 1 m 2 in the χ transforms.

Comparison of continuity steady state and simulated profiles
For simulations where n = 1, channel profiles far from base level approach a steady configuration, in which channel slope in χ space is a unique function of rock erodibility, and the profiles exhibit straight line segments in each rock type (Figs. 4,5).
For the horizontally layered case, channel profiles evolve towards a state in which they are maintaining the same shape in χ space while retreating horizontally into the bedrock.For small changes in basin area, this is equivalent to a channel maintaining constant horizontal retreat rates.For non-horizontal rocks, profile shapes will gradually change in χ space, as the slope of the contact plane in χ space changes with basin area.Animations of the simulations depicted in Figs. 4 and 5 are provided in the online supplementary material.
For n = 1 there is no one-to-one relation between erodibility and slope, and the profiles do not exhibit straight-line segments in each rock type.The n = 1 case produces this result because the horizontal retreat rates are independent of slope and purely a function of erodibility and basin area.Consequently, adjustments of slope cannot produce equal horizontal retreat rates along the channel.Instead, segments within weaker rocks will retreat more quickly than those within stronger rocks.This produces "stretch zones" as a channel crosses from weak to strong rocks and "consuming knickpoints" as a channel crosses from strong to weak rocks (Royden and Perron, 2013;Forte et al., 2016).The channels in the simulations ultimately reach a steady stepped  ,5C) in which weak rock layers retreat until they intercept the contact with strong layers; however, this state will contain near-vertical channels that are sufficiently steep to negate the applicability of the stream power erosion law.Within simulations, the nature of such profiles may be strongly dependent upon the numerical algorithm employed (Campforts and Govers, 2015).
The continuity relation (3) predicts a slope ratio rather than absolute values of slope in each rock type.The predicted slope ratio matches the slopes in the simulation at sufficient distances from base level.Notably, the counterintuitive prediction that profiles would be steeper in weaker rocks for n < 1 is confirmed by the simulations (Figs.4A,5A).However, absolute slopes, and therefore entire profiles, can be predicted by realizing that continuity steady state is actually a type of flux steady state (Willett and Brandon, 2002), where the rate of uplift of rock into the domain is equal to the rate of removal of material by erosion.First, it must be noted that the weak and strong rocks experience different rates of vertical incision in the equilibrium state (Forte et al., 2016).The rates of vertical incision within two alternating layers at equilibrium can be calculated from an assumption of flux steady state, or, equilvalently, an assumption that the time-averaged incision rate through both rock types is equal to the uplift rate.This leads to a relation for erosion rate in a given layer, where E 1 is the erosion rate of one rock layer, H i is the thickness of the ith layer measured in the vertical direction, and U is the uplift rate (see derivation in Appendix B).Entire theoretical profiles can be constructed using this relationship, in combination with the stream power erosion law and the continuity relation (Eq.5), which provides the slope ratio.These profiles closely match the simulations in cases where n = 1 at a sufficient distance from base level (Figs.4A-B,5A-B).Therefore, continuity state is a type of flux steady state.In addition to describing behavior near contacts, it also describes portions of the profile that are distant from contacts.For sub-horizontal rocks this often produces a landscape that is quite different from that which would be predicted by topographic steady state (Fig. 3).

Dynamics of base level perturbations
Continuity steady state is perturbed near base level, where a constant rate of base level fall is imposed.This results because vertical incision occurs at different rates in each rock in continuity steady state, whereas uplift is constant.Despite this disequilibrium near base level, theoretical profiles produced using Eq.(3) and Eq. ( 8) closely match the shapes of the profiles for the cases where n is not one.Therefore, these perturbations decay rapidly away from base level in the simulated cases.However, a question remains as to what controls this decay length scale, and how typical the cases are that we have simulated.
In a horizontally layered rock sequence, a segment of stream profile with erosion rate equal to uplift is continuously developing at base level.The slope of this base level segment in χ-space is given by The difference between this slope, and the continuity steady state slope produces a knickpoint that propagates upstream with a celerity given by As the knickpoint crosses into the other rock type, continuity demands that C does not change, because C is identical to horizontal retreat rate and continuity requires this to be equal across a horizontal contact.Since celerity is a monotonic increasing function of erodibility, knickpoints formed at base level in the stronger rock are slower than those formed in the weak rock.
Therefore, the weak rock knickpoints catch up to the strong rock knickpoints, and the profile damps toward equilibrium as the two interact.Consequently, we can estimate the damping length scale as the χ distance at which the knickpoints generated in weak rock at base level catch up to the knickpoints generated in strong rock at base level.
The strong rock knickpoint begins with a head start equal to the χ distance spanned by the strong rock segment, which we call χ s and is given by The strong rock knickpoint will travel an additional distance χ s,+ before the weak rock knickpoint catches up, and these distances are related by where C s and C w are the knickpoint celerities in the strong and weak rocks, respectively.The damping length scale, λ = χ s,0 + χ s,+ , is the distance from base level over which the weak rock knickpoint catches the strong one and can be solved for by combining Eqs. 10, 11, and 12, leading to To generalize the damping behavior of the base level perturbations it is useful to analyze a dimensionless version of λ, which is normalized by χ s,0 , It can be seen that the damping length scale is primarily a function of the relative erodibilities of the two rock types.When the contrast is large, damping occurs rapidly, whereas when the contrast is small the damping length scale is large.However, in this latter case there is also very little contrast in steepness, since the erodibilities are similar.We show λ * as a function of the erodibility ratio for several choices of n in Figure 6.Here it can be seen that if the erodibility ratio is greater than about two or three damping occurs within two cycles through the rock layers away from base level.If the erodibility ratio is greater than about ten, then damping occurs within a single cycle.
To illustrate this damping behavior, we run two simulations with somewhat longer damping length scales.Both simulations have profile lengths of 500 km, uplift rates of 2.5 mm y −1 , repeating rock layers with a 50 m thickness, and weak rock layers Profiles are shown for these simulations in Figure 7. Fast knickpoints catch the slow knickpoints at roughly the calculated length scale (Fig. 7C,D).This is visualized more clearly in animations in the supplementary material that depict the damping length scale.Beyond this damping length scale, some minor perturbations remain, and one can see fast and slow knickpoints migrating through the upper parts of the profile as the system evolves.However, beyond λ the theoretical profiles derived from continuity and flux steady state are good approximations to profile shape.

Full landscape simulations
To determine whether continuity steady state is obtained within whole landscape models, or whether addition of hillslope processes might eliminate it, FastScape V5 (Braun and Willett, 2013) was used to simulate stream power erosion coupled to an entire landscape model.All simulated cases employ a constant rock uplift rate and horizontal rock layers with alternating high and low erodibility.
The stream power erosion law used in FastScape has the form where Φ is discharge, calculated as the product of the drainage area and the precipitation rate P .Each of the three presented model runs uses two different erodibility coefficients, K f w for the weak rock and K f s for the strong rock, in place of K f .For each one of them, a grid of 3000 × 3000 pixels representing 100 × 100 km is simulated.The initial condition used is a slightly randomly perturbed flat surface at base level.The boundary condition is open on all sides.15000 m of uplift is simulated in 60000 timesteps.The weaker rock is exposed for the first 10800 m of the uplift, allowing an initial drainage network to establish.Afterwards, a layered rock structure starts to be exposed, with alternating layers of 200 m of the stronger rock and 300 m of the weaker rock.The different bed thicknesses enables testing of whether any of the previous theoretical results require the layered rocks to have equal thickness.The main difference between the model runs is in the slope exponent n, with cases using n = 2/3, n = 1, and n = 3/2.A listing of numerical parameters is provided in Table 2.The necessary timestep was calculated from the uplift rate and the ratio of total uplift to the number of timesteps.
Floating point digital elevation models (DEMs) were produced for the final time step for each FastScape simulation.Using the Landlab landscape evolution model (Tucker et al., 2013) to calculate flow routing, channel profiles were extracted from the FastScape DEMs for each case of n.Landlab was extended to enable calculation of χ values for each channel.χ-plots were then generated for 50 channels in each simulation and are shown in Fig. 8.The continuity equilibrium state described above is also reflected within the full landscape evolution model, and plots of elevation versus χ for channels within each model demonstrate similar relationships as displayed in Figure 4A,C,E.

Discussion and conclusions
The implications of topographic equilibrium for channel form break down within layered rocks as the layers approach horizontal (Howard, 1988;Forte et al., 2016).However, under constant forcing, the stream power erosion law drives channels to a different type of steady state, which we refer to as continuity steady state.Continuity steady state is also a type of flux steady state, where time averaged incision in the channel at any given horizontal position is equal to uplift.In this state, channels maintain a constant normalized steepness (gradient in χ space) within each layer for most values of n, the exponent in the stream power erosion law.In the horizontal case, a singularity occurs at n = 1, in which case continuity cannot be maintained through the adjustment of slope.
If constant uplift occurs, the channel does not equilibrate at base level, because continuity steady state results in different vertical incision rates in each rock type.However, the perturbations introduced by this disequilibrium at base level rapidly decay over a length scale that is primarily a function of the ratio of rock erodibilities, with larger erodibility contrasts resulting in shorter decay lengths.Practically speaking, for rocks that have erodibilities sufficiently different to have a strong effect on the profile, base level perturbations decay after a couple rock contacts are passed.
When contacts between rocks dip at slopes much greater than the channel slope, then the vertical contact limit from Eq.
(4) applies and standard conceptions derived from topographic steady state hold.The considerations introduced here become important as rock dips approach values comparable to or less than channel steepness.This subhorizontal limit, given by Eq. ( 5), is most likely to apply for rocks that are very near horizontal and/or channels that are very steep.Therefore, these considerations are most applicable in cratonic settings, in headwater channels, or when considering processes of scarp retreat in subhorizontal rocks (Howard, 1995;Ward et al., 2011).In the subhorizontal limit, slope contrasts are larger than would be predicted by topographic steady state.In the case of n < 1, slope patterns in continuity steady state are also qualitatively different than those predicted by topographic steady state, with steeper channel segments in weaker rocks.For n ∼ 1 slope contrasts become extreme, which is particularly important since n = 1 is the most common value used in landscape evolution models.In this case, large slope contrasts at contacts may accentuate numerical dispersion.It also must be realized that n = 1 is quite a special case in subhorizontal rocks, and the rest of the parameter range for n results in substantially different dynamics and steady state.
The framework presented here provides a possible hint concerning the generation of caprock waterfalls, a feature that has long fascinated geologists (Gilbert, 1895).Caprock waterfalls, such as Niagara Falls, have a resistant caprock layer that is underlain by a weaker rock.The waterfall has the caprock at its lip, followed by a vertical, or often overhanging, face within the weak rock.This is a case of a very steep channel within a highly erodible rock, which would not be predicted from topographic equilibrium and stream power erosion.Such a state is predicted by the continuity relation developed here for subhorizontal layers with n < 1.Values of n might be expected to fall in this range for erosion processes active in the weak rock layer, such as plucking (Whipple et al., 2000).Furthermore, caprock waterfalls typically form in relatively horizontal strata, and are common within steep headwater channels.The stream power law arguably does not apply to waterfalls (Lamb and Dietrich, 2009;Haviv et al., 2010;Lague, 2014), and a variety of erosion mechanisms that are independent of stream power can act in such an oversteepened reach.However, starting from an initial condition of low relief, topographic equilibrium would not predict a channel to evolve toward the caprock waterfall state.In contrast, the framework presented here naturally produces features resembling caprock waterfalls from considerations of landscape equilibrium.While further work would be needed to test this hypothesis, it remains plausible that caprock waterfalls are the result of channels steepened within weaker rocks to maintain continuity.Once the channel becomes sufficiently steep, stream power erosion may no longer provide a good approximation to erosion rates, though the concept of continuity could also be applied to other, more mechanistic, erosion models.
Though the focus of this work is on bedrock channel profiles in layered rocks, the concepts of continuity and flux steady state can be applied in general to any mathematical model for erosion.Much like topographic steady state, both continuity and flux steady state result from negative feedbacks within the uplift-erosion system that drive it toward steady state as uplift and erosion become balanced.Though topographic steady state has been a powerful theoretical tool to understand landscapes, the generalized concept of continuity may prove more useful in interpreting steep landscapes in subhorizontal rocks.

Appendix A: Derivation of the continuity relation
Here we detail how the constraint of channel continuity can be used to derive the relationship given in Eq. ( 3).Consider a planar contact between lithologies with different erodibilities.We label the downstream and upstream erodibility with K 1 and K 2 .Downstream and upstream slopes are S 1 and S 2 , the slope of the contact is S c , and their respective slope angles are θ 1 , θ 2 , and φ, see Fig. A1.
In this section we use the subscript i to denote either 1 or 2, as the relationships are valid for the channel within both rock types.Erosion at a rate E i in the vertical direction, as is calculated by the stream power erosion law, can be transformed to an erosion rate B i that is perpendicular to the channel bed using the slope of the channel bed, θ i , with B i = E i cos θ i , see Fig. A1.
The contact and the channel intersect at angle θ i + φ, thus the rate of exposure of the contact plane is For the case where θ i + φ > π/2 the diagram changes, but these same relationships can be recovered using sin(π − θ i − φ) = sin(θ i + φ).Continuity of the channel bed requires that the contact exposure rates R 1 and R 2 are equal, which gives Using a trigonometric identity for angle sums leads to Simplifying the fractions and multiplying both sides of the equation with cos φ we get Solving for the ratio of erosion rates in the two rock types and converting to slopes rather than angles, the relation becomes If erosion rates are given by the stream power law, then it follows that which is identical to Eq. ( 3) with the general subscripts 1 and 2 replaced with s and w for strong and weak.

Appendix B: Derivation of the erosion relation
Using the stream power erosion law, erosion rates in two channel segments above and below a contact are where A is the recharge area.Taking the ratio of both equations at an arbitrary basin area, we get We define H 1 and H 2 to be the thicknesses of the rock layers measured in the vertical direction.If flux steady state is assumed, then the average erosion rate equals the uplift rate U .Therefore, the time needed to uplift a distance equal to the sum of the two layers' thicknesses equals the sum of the times needed to erode through the two layers: Combining Eq. (B2) and Eq.(B3) gives an expression for the erosion rate in a given rock: While flux steady state seems like a reasonable assumption, simulations also confirm that the erosion rates predicted by Eq.
(B4) are approached within a few contacts above base level.Similarly, simulations that alternate uplift rate over time to match the erosion rate of the rock type currently at base level, as given by Eq. ( B4), obtain straight line slopes in χ-elevation space all the way to base level.This confirms that the disequilibrium seen in the profiles in Fig. 4B-D   Simulation Earth Surf.Dynam.Discuss., doi:10.5194/esurf-2016-41,2016   Manuscript under review for journal Earth Surf.Dynam.Published: 18 August 2016 c Author(s) 2016.CC-BY 3.0 License.time and was computed as a function of longitudinal distance, with Earth Surf.Dynam.Discuss., doi:10.5194/esurf-2016-41,2016   Manuscript under review for journal Earth Surf.Dynam.Published: 18 August 2016 c Author(s) 2016.CC-BY 3.0 License.shape(Figs.4C
is produced by the difference between the constant uplift rate and the equilibrium incision rates experienced in each layer.Earth Surf.Dynam.Discuss., doi:10.5194/esurf-2016-41,2016 Manuscript under review for journal Earth Surf.Dynam.Published: 18 August 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 2 .
Figure 2. Erosional continuity.(A) If the upper layer at a contact erodes slower, this produces a discontinuity at the contact and the resulting steepening in the upper layer will drive the system toward erosional continuity.(B) If the upper layer erodes faster, this produces a low slope zone near the contact, which will reduce erosion rates in the upper layer (C) In general, the system will approach a state where continuity is maintained.

Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .
Figure3.Channel slope response at a subhorizontal contact from an assumption of continuity.The ratio of slope within the weaker rock (Sw) and the slope within the stronger rock (Ss) near a horizontal contact (solid line) with differing values of the power n in the stream power erosion law.Erodibility in the weaker rocks (Kw) is twice that of the stronger rocks (Ks).The dashed line displays the traditional relationship, which applies for cases where the contact slope is much larger than the channel slope.

Figure A1 .
Figure A1.Geometric relationships used to derive the equation for continuity of the channel at a contact between two rock types.Note that the slope of the contact plane (φ) is positive when the contact dips in a direction opposite that of channel bed slope.

Table 2 .
Parameters used in the FastScape model runs.