Supplement of Bedrock river erosion through dipping layered rocks: quantifying erodibility through kinematic wave speed

Abstract. Landscape morphology reflects drivers such as tectonics
and climate but is also modulated by underlying rock properties. While
geomorphologists may attempt to quantify the influence of rock strength
through direct comparisons of landscape morphology and rock strength
metrics, recent work has shown that the contact migration resulting from the presence of mixed lithologies may hinder such an approach. Indeed, this work counterintuitively suggests that channel slopes within weaker units can sometimes be higher than channel slopes within stronger units. Here, we expand upon previous work with 1-D stream power numerical models in which we have created a system for quantifying contact migration over time. Although previous studies have developed theories for bedrock rivers incising through layered stratigraphy, we can now scrutinize these theories with contact migration rates measured in our models. Our results show that previously developed theory is generally robust and that contact migration rates reflect the pattern of kinematic wave speed across the profile. Furthermore, we have developed and tested a new approach for estimating kinematic wave speeds. This approach utilizes channel steepness, a known base-level fall rate, and contact dips. Importantly, we demonstrate how this new approach can be combined with previous work to estimate erodibility values. We demonstrate this approach by accurately estimating the erodibility values used in our numerical models. After this demonstration, we use our approach to estimate erodibility values for a stream near Hanksville, UT. Because we show in our numerical models that one can estimate the erodibility of the unit with lower steepness, the erodibilities we estimate for this stream in Utah are likely representative of mudstone and/or siltstone. The methods we have developed can be applied to streams with temporally constant base-level fall, opening new avenues of research within the field of geomorphology.



S1: Extracting and analyzing stream profile data
We used TopoToolbox v2 (Schwanghart and Kuhn, 2010;Schwanghart and Scherler, 2014) to extract the stream profile data shown in Figs. 1, 14, and 15. We extracted these data from 10-m digital elevation models (DEMs) provided by the United States Geological Survey. For each reach assessed (i.e., the five streams with steepness values in Fig. 1a, which includes Tank Wash), we first used the getnal function in TopoToolbox to extract the elevations, flow distances, and drainage areas for a flow path along the filled DEM (filled using ArcGIS). We then selected a minimum drainage area for the flow path, removing all nodes upstream of the minimum value. Minimum drainage areas were manually selected in slope-area plots as values that (1) exceeded 0.1 km 2 and (2) occur at a transition towards decreasing channel slopes with increasing drainage area (Montgomery and Foufoula-Georgiou, 1993). To account for the noise common within DEMs (Wobus et al., 2006), we smoothed river profile elevations over 55 nodes (where node spacing is either 10 m or ~14.1 m). We chose this smoothing interval because it minimized scatter on slope-area plots while still preserving the shape of the profile.
We used two approaches for measuring channel steepness (ksn) in this study: (1) one approach for the data shown in Fig. 1 and (2) one approach for both our estimation of kinematic wave speed with Eq. 13 and our preliminary analysis of Tank Wash (Figs. 14-15). The channel steepness data shown in Figs. 1b and 1c are calculated as (Wobus et al., 2006): where S is the channel slope calculated with the smoothed channel elevations (where the slope at each node is calculated in the downstream direction), A is drainage area [L 2 ], and m/n is taken as 0.5. This same approach was used for the channel steepness values shown in Fig. 1a, although Fig. 1a displays the average steepness values within 250-m intervals along the channels. We used this "at-a-point" channel steepness for Fig. 1 because that figure is only used for demonstrative purposes (i.e., those data are not used in further analyses).
Besides Fig. 1, we used ksn values for both (1) calculating kinematic wave speed with Eq. 13 (Figs. 9-10) and (2) our preliminary analysis of Tank Wash . In these applications, we calculated steepness values as (Perron and Royden, 2013): where Δz is the change in elevation (L) from upstream to downstream across the selected reach (defined by contact locations), Δχ is the change in transformed river distance (L) (Eq. 4) from upstream to downstream across the reach, A0 is reference drainage area (1 km 2 ), and m/n is taken as 0.5. Although we only assess an m/n value of 0.5 in the main text of this study, the variability of this parameter (and the channel concavity it is often thought to influence) is a significant issue in the field of geomorphology Supplement (Mudd et al., 2018). We present examples of simulations using different m/n values in Figs. S6-S10 and S12-S16 below.

S2: Approach for setting erodibility values
Because we examine different n values here, we set the range of erodibility by considering slope patch migration rates (Royden and Perron, 2013): where χsp is the χ value of a slope patch and U is the rock-uplift rate (or base-level fall rate) at the time t when the slope patch was generated at the basin's outlet (constant over time here; U ≠ f(t)). For the sake of concision, we express slope patch migration rates as Cχ. We calculate the strong layer's erodibility so that the Cχ for the strong layer is a certain fraction (33, 50, 67, 75, or 90 %) of the Cχ for the weak layer.
Using this approach, we set the weak layer's erodibility (KW) to a reference value (Table 1) and calculate the strong layer's erodibility (KS) as: where nS and nW are the n values of the strong and weak layers and CχS and CχW are the Cχ values (Eq. S3) for the strong and weak layers, respectively. We calculate different KS values by varying (CχS / CχW) from 0.33 to 0.9 (i.e., the percentages listed above). Note that Eq. S4 was calculated by expressing (CχS / CχW) with Eq. S3 and then rearranging for KS. Although these slope patch migration rates are calculated with U, these rivers do not always erode at U. The rock-uplift rate is still a physically representative rate, however (i.e., erosion rates will vary around U). Most importantly, this approach to scaling erodibility is consistent across different n values. Note that while we include nS and nW in Eq. S4, we only assess simulations where the strong and weak layers have the same slope exponent n values. If nS = nW, then KS simplifies to KW × (CχS / CχW) n .
We use three reference weak erodibilities (KW) for each n value (Table 1). We chose the low, medium, and high reference KW values so that they allow slope patch migration rates equal to those for n = 1, U = 0.15 mm yr -1 , and K values of 5×10 -7 yr -1 , 10 -6 yr -1 , and 2×10 -6 yr -1 , respectively. This range of KW for n = 1 was based on values reported in the literature (e.g., (Armstrong et al., 2021)).
In the simulations using three rock types (scenario 2; Table 1), the erodibilities for the medium unit are scaled to be between those for the weak and strong units. For each combination of weak and strong erodibilities (KW and KS), three erodibilities are evaluated for the medium rock type (KM): (1) the KM that provides the average of the Cχ (Eq. S3) in the strong and weak units (for rock-uplift rate U); (2) the KM that provides a Cχ value that is three times closer to that of the strong unit (for rock-uplift rate U); and (3) the KM that provides a Cχ value that is three times closer to that of the weak unit (for rock-uplift rate U). For example, if the KS value provides a Cχ that is 50 % of the Cχ for the weak unit, KM would be varied to provide slope patch migration rates that are 62.5, 75, and 87.5 % of those in the weak unit for three separate model simulations.

S3: Recording contact migration rates
To track the movement of specific contacts in our numerical models, there must be a system for determining if each contact was observed in the last time step evaluated or if it is a new contact. Our approach for this involves assigning each contact a designation (e.g., contact #1, #2, etc.). When our automated system proceeds to the next time step, it then determines if the first contact (contact with the lowest elevation) is the same as the first contact from the last time step. This comparison considers which units are above and below the contact; if the first contact in the previous time step had the weak unit beneath it, but the first contact now has the strong contact beneath it, then this first contact is a new contact (i.e., the first contact from the previous time step has migrated further upstream, and the newly identified contact has only just entered the region where we record contact migration rates). We also use the recorded χ values at each contact. Consider, for example, if the first contact encountered sits over the same rock type as the first contact encountered in the previous time step but the χ value of the newly encountered contact is lower. In that scenario, two new contacts have just entered the region where contact migration is recorded. Note that this approach requires that contacts cannot migrate too fast relative to the time interval over which contacts are tracked (every 25 kyr, larger than model timestep dt = 25 yr). For example, a contact cannot migrate so fast that it crosses the entire profile over the 25 kyr interval.
The aspect of the contact tracking system we stress here is that recording contact migration rates near the basin outlet is inadvisable because of base-level perturbations (Perne et al., 2017). We describe how we avoid the problems caused by base level perturbations further below, but first we provide some background on the base-level perturbations themselves. Previous work (Forte et al., 2016;Perne et al., 2017;Darling et al., 2020) shows that contact migration dynamics can cause much of the profile to have erosion rates that are not equal to the rock-uplift rate (U). Near the outlet, however, the constant forcing of the base level fall rate at this hard boundary means that an equilibrated reach eroding at U is constantly being created at the basin outlet. This reach with E = U is at odds with the rest of the profile, where E is instead a function of contact migration dynamics. As the units are steadily raised by rock uplift, the equilibrated reach near the outlet is perturbed when a different rock type is exposed beneath it. For example, a reach in the strong unit may be connected to the outlet (causing it to have E = U), only for the weak unit to appear beneath it. Then, an equilibrated reach in the weak unit with E = U begins to form, and the previously equilibrated reach in the strong unit is disconnected from the basin outlet. The erosion at the base of the previously equilibrated reach in the strong unit will then change in order to conform with the dynamics of the contact's migration. This change in erosion rate can, in certain situations, create a consuming knickpoint (Royden and Perron, 2013) that migrates upstream. These knickpoints are the reason why recording contact positions near the basin outlet is problematic for our comparison of measured contact migration rates and estimated kinematic wave speeds; one would record pronounced variations in contact migration rates near the outlet as reaches initially have E = U and are then perturbed by knickpoints once they are disconnected from the basin outlet. To test the equations for kinematic wave speed we examine here, we must examine these equations where they are applicable (i.e., not where knickpoints from base level effects dominate). The expression of these knickpoints will, however, sufficiently diminish at some distance upstream from the basin outlet. To characterize the upstream extent of these knickpoints, Perne et al. (2017) derived the damping length scale λ. This damping length scale defines the χ distance from the outlet at which the effects of such disturbances have largely disappeared: where HS is the layer thickness of the stronger unit. Note that λ is calculated as the χ-distance from the outlet at which a knickpoint migrates from the bottom to the top of a reach underlain by the strong layer (Perne et al., 2017). We avoid such base-level perturbations by only recording and analyzing contact migration rates upstream of the location where χ = λ. Depending on the parameters used in a simulation, however, λ can sometimes be larger than the highest χ value. Otherwise, if λ is too close to the highest χ value, this may leave only a few stream nodes for analysis. Tracking contact migration over such a small portion of the stream profile can cause issues for our approach (e.g., a reach extending across one rock type may not fit within the area). We therefore only record contact migration if a minimum proportion of the profile can be used (25 % of the profile for scenarios 1-3 and 10 % of the profile for scenario 4; these choices were based on model observations). We do not utilize results for simulations that do not meet these criteria, but this loss is alleviated by the large number of simulations we assess.
Quantifying contact migration in some of our numerical models requires additional considerations, however. Such models include (1) those with three rock types instead of two and (2) those where contacts dip downstream. When we evaluate numerical models with three rock types (weak, medium, and strong), the issue is that damping length scale λ is calculated using the parameters for only two rock types. When we use three rock types, we therefore calculate λ between the weak and medium layers and between the medium and strong layers. We use the larger of these two λ values to set the starting position for contact tracking. Although this approach for simulations with three rock types is simple, there are far more considerations for numerical models with contacts dipping in the downstream direction. For example, the conceptual model in Fig. 2d has contacts dipping in the downstream direction. Based on previous work (Darling et al., 2020), we anticipate that locations where channel slope is less than the contact slope will have contacts that migrate in the downstream direction. We do not focus on contact migration in the downstream direction in this study; when contacts migrate in the downstream direction, they migrate in the opposite direction of base level signals. As a result, these contacts (1) do not have the same capacity for driving feedbacks between erosion and contact migration and (2) do not offer the same opportunities for extracting information regarding rock strength in the manner we propose in this study. Downstream-migrating contacts may still offer other kinds of valuable information, but we restrict this study's focus to upstream-migrating contacts.
Because of the change in contact migration direction when contacts dip downstream (ϕ > 0), the position where χ equals the damping length scale λ (Eq. S5) is not always appropriate place to set the starting position for tracking contact migration. To define an appropriate starting position, we run the simulation for 50 Myr (so that a dynamic equilibrium is achieved) and find the lowest position where the absolute value of average channel slope exceeds the contact slope. Because there are variations in slope within each rock type, we use the average steepness of the profile to find where the average slope exceeds contact dip. We used the average steepness because it will not change significantly over time, once a dynamic equilibrium is achieved. Instead of using the point where average channel slope equals the contact slope as a starting position, however, we advance the starting position further upstream. We made this decision because the exact area where a new unit is exposed varies slightly, which can be problematic for our contact tracking approach. To avoid such complications near the point where new units are exposed, utilize the expected χ-distance a reach extending across both the strong unit and a weak unit would span if the dip was zero (ϕ = 0°). The distance of a reach extending across one rock type Δχ can be solved using the theory developed by Perne et al. (2017): where ΔχS is the χ-distance required for a reach spanning the strong unit when ϕ = 0°. Erosion rate ES can be solved with Eq. 10. To solve for ΔχW with Eq. S6, the variables HS, ES and KS would simply have to be changed to HW, EW, and KS. The length scale we use to adjust the starting position is then ΔχS + ΔχW. This distance was chosen because it varies with simulation input such that it is never too large or small relative to the profile's elevation range (e.g., one fixed distance would not be applicable across the full range of simulations).
Through a process of trial and error, we found that simulations using lower reference weak erodibilities (KW) require a larger distance between the point where new units are exposed and the starting position for contact tracking. With lower KW, there is more variability in the locations where new units are exposed. When the starting position for contact tracking is too close to areas where new units are exposed, the variations in erosion rate can cause Eqs. 12 and 13 (Fig. 10) to become less accurate. To avoid such issues, we define the starting position for contact tracking when ϕ > 0° as follows: where ζ is the χ value defining the starting position for contact tracking, χexposure is the χ value closest to the point where average channel slope is equal to the contact slope (where new units will likely be exposed), and F is a factor set by the user through trial and error. In simulations with the low, medium, and high KW, we use F values of one, three, and four, respectively. These F values allowed for accurate estimations of kinematic wave speed using Eqs. 12 and 13 (Fig. 10).
When the contact dip is low enough (or the erodibility values are high enough), however, there will be no inflection point (χexposure = 0 m) and all contacts will migrate upstream. To address such situations, we evaluate the starting positions described by both λ (Eq. S5) and ζ (Eq. S7) and use whichever starting position is farther upstream.
The results for some of our simulations were not used because the damping length scale λ (Eq. S5; (Perne et al., 2017)) occurred at a χ value that was too high (e.g., beyond the maximum χ). The damping length scale is used to define the distance from base level at which the influence of constant base level forcing disappears (i.e., consuming knickpoints migrating through the strong layer). Upstream of this location, the dynamics of contact migration are fully expressed. One might wonder if real streams can fail to fully express the influence of contact migration in a similar manner; can the length scales required for these dynamics to develop be too large to occur in a real landscape? Such complications would hinder our ability to extract rock strength information from landscape morphology. One then must wonder, however, where a hard boundary in base level forcing would occur in a real landscape. Our numerical models have a basin outlet that acts as a hard boundary, and the constant base level fall at this boundary is the reason why we use the damping length scale. Perhaps such a boundary could effectively occur where a bedrock-dominated tributary flows into a larger, transport-limited river. Such a relationship between base level forcing and a tributary's response remains to be demonstrated, however.

S4: Erodibilities for the simulations in Figs. S6-S10 and S12-S16
The simulations in Figs. S6-S10 and S12-S16 use m/n values of 0.3, 0.5, or 0.7. Although we use a reference weak erodibility (KW) to calculate the medium and strong erodibilities (KM and KS) for our main simulations (scenarios 1-4; Table 1), such an approach is less convenient in simulations where m/n changes (i.e., the dimensions of K depend on m). In the simulations in Figs. S6-S10 and S12-S16, we therefore chose to calculate the KW value as one that would provide a fluvial relief of 1500 m (if only the weak layer was present and E = U): where zmax ref is the reference maximum elevation (L) (zmax ref = 1500 m) and χmax is the maximum χ value (L). The erodibility for the strong layer (KS) is then calculated from KW in the same manner described in Sect. S2 (Eq. S4). Because we present a limited number of examples varying m/n, however, we use specific (CχS / CχW) values (Eq. S4): (CχS / CχW) = 0.33 for n = 0.67 and (CχS / CχW) = 0.5 for n = 1.5. We made these choices because those (CχS / CχW) values were used in Figs. 3, 7, and 8 (for the corresponding n values). While the six simulations in Figs. S6-S10 use ϕ = 0°, the six in Figs. S12-S16 use ϕ = -2.5°. All simulations in Figs. S6-S10 and S12-S16 use a rock-uplift rate U of 0.15 mm yr -1 and a layer thickness H of 100 m. Figure S1. Variations in maximum elevation over time t (zmax(t)) normalized by the final maximum elevation (zmax(tmax)) for the simulations in (a) Fig. 3a and (b) Fig. 3b. Note how the maximum elevation can gradually increase or decrease before settling upon a consistent range. (c) Normalized maximum elevations over time for all simulations assessed in scenario 1 (two layers with contact dip ϕ = 0°; Table  1). Although there are variations in maximum elevation, the range of elevations is constant over time (after the 50 Myr periods of initialization). These data show that the 50 Myr initialization period allowed all simulations assessed in scenario 1 to achieve a dynamic equilibrium such that the range of elevations was constant with time. Larger variations in maximum elevation occur for larger contrasts in erodibility. Figure S2. Variations in maximum elevation over time t (zmax(t)) normalized by the final maximum elevation (zmax(tmax)) for the simulations in (a) Fig. 6a and (b) Fig. 6b. Note how the maximum elevation can gradually increase or decrease before settling upon a consistent range. (c) Normalized maximum elevations over time for all simulations assessed in scenario 2 (three layers with contact dip ϕ = 0°; Table  1). Although there are variations in maximum elevation, the range of elevations is constant over time (after the 50 Myr periods of initialization). These data show that the 50 Myr initialization period allowed all simulations assessed in scenario 2 to achieve a dynamic equilibrium such that the range of elevations was constant with time. Larger variations in maximum elevation occur for larger contrasts in erodibility. Figure S3. Variations in maximum elevation over time t (zmax(t)) normalized by the final maximum elevation (zmax(tmax)) for the simulations in (a) Fig. 7a and (b) Fig. 8a. Note how the maximum elevation can gradually increase or decrease before settling upon a consistent range. (c) Normalized maximum elevations over time for all simulations assessed in scenario 3 (two layers with contact dip ϕ < 0°; Table  1). Although there are variations in maximum elevation, the range of elevations is constant over time (after the 100 Myr periods of initialization). These data show that the 100 Myr initialization period allowed all simulations assessed in scenario 3 to achieve a dynamic equilibrium such that the range of elevations was constant with time. Larger variations in maximum elevation occur for larger contrasts in erodibility. Figure S4. Variations in maximum elevation over time t (zmax(t)) normalized by the final maximum elevation (zmax(tmax)) for the simulations in (a) Fig. 7b and (b) Fig. 8b. Note how the maximum elevation can gradually increase or decrease before settling upon a consistent range. (c) Normalized maximum elevations over time for all simulations assessed in scenario 4 (two layers with contact dip ϕ > 0°; Table  1). Although there are variations in maximum elevation, the range of elevations is constant over time (after the 50 Myr periods of initialization). These data show that the 50 Myr initialization period allowed all simulations assessed in scenario 4 to achieve a dynamic equilibrium such that the range of elevations was constant with time. Larger variations in maximum elevation occur for larger contrasts in erodibility. Figure S5. Longitudinal profile for a simulation using n = 1. Here, the strong layer (dark gray) has erodibility KS = 5×10 -7 yr -1 , while the weak layer (light gray) has KW = 10 -6 yr -1 . Here, the rock-uplift rate U is 0.15 mm yr -1 , and the layer thickness H is 100 m. At a distance upstream of the outlet sufficient to obscure the influence of base-level effects, the stream's slopes are either zero or infinite (zero in the flat reaches, infinite at the steps).       Fig. 6b. The slope exponent (n) and weak, medium, and strong erodibilities (KW, KM, and KS) are shown in the upper left of each subplot. Note that contact dip (ϕ) here is 0°. Each subplot has dotted, dash-dot, and dashed lines for the kinematic wave speeds estimated for the weak, medium, and and strong layers, respectively, if the erosion rate (E) was equal to rock-uplift rate (U). The erosion rates in each layer do not conform to this assumption, however, and instead vary so that kinematic wave speed is maintained at a moderate value between these different lines. Figure S12. Longitudinal profiles (a-c) and χ-plots (d-f) for three different simulations. All simulations have n = 0.67, contact dips ϕ of -2.5°, rock-uplift rates U of 0.15 mm yr -1 , and layer thicknesses H of 100 m. The simulations have different m/n values, however (0.3, 0.5, or 0.7). The erodibility K values have different dimensions, impeding their direct comparison. Each simulation has the same K * value (Eq. 7c), however, allowing similar variations in erosion rate. Note that the distribution of peaks in erosion rate varies with m/n. With a lower m/n value of 0.3, the peaks in erosion rate have a smaller range and rate of change along the profile. Conversely, with a higher m/n value of 0.7 the peaks in erosion rate have a larger range and rate of change along the profile. Figure S13. Longitudinal profiles (a-c) and χ-plots (d-f) for three different simulations. All simulations have n = 1.5, contact dips ϕ of -2.5°, rock-uplift rates U of 0.15 mm yr -1 , and layer thicknesses H of 100 m. The simulations have different m/n values, however (0.3, 0.5, or 0.7). The erodibility K values have different dimensions, impeding their direct comparison. Each simulation has the same K * value (Eq. 7c), however, allowing similar variations in erosion rate. Figure S14. Contact migration rates (dxcontact/dt) vs drainage area (A) for the (a) the simulation shown in Figs. S12a and S12d, (b) the simulation shown in Figs. S12b and S12e, and (c) the simulation shown in Figs. S12c and S12f. Note that there is no lower error bar on the point for the lowest drainage area bin in subplot (c). The measured dxcontact/dt for that bin is about 7.84×10 -5 m yr -1 while the standard deviation (represented by error bars) is about 9.77×10 -5 m yr -1 , causing the lower error bar to extend to -1.92×10 -5 m yr -1 . Figure S15. Contact migration rates (dxcontact/dt) vs drainage area (A) for the (a) the simulation shown in Figs. S13a and S13d, (b) the simulation shown in Figs. S13b and S13e, and (c) the simulation shown in Figs. S13c and S13f. Figure S16. Contact migration rates measured in the simulations in Figs. S12-S13 (dxcontact/dt) vs. kinematic wave speeds (CH) estimated using Eqs. 12-13. Subplot (a) shows CH estimated with Eq. 12 while subplot (b) shows CH estimated with Eq. 13. Figure S17. Contact migration rates (dxcontact/dt) vs drainage area (A) for the four simulations shown in Figs. 7 and 8. In contrast with Fig. 9, this presentation of the data uses 20 drainage area bins instead of 10. The performances of Eqs. 12 and 13 are similar but not exactly the same (e.g., the R 2 for Eq. 12 in subplot (b) is 0.41 instead of 0.37). All simulations here have layer thickness H = 100 m. The simulations have slope exponent (n) and weak and strong erodibility values (KW and KS) of (a-b) n = 0.67, KW = 6.83×10 -6 m 0.33 yr -1 , and KS = 3.25×10 -6 m 0.33 yr -1 and (c-d) n = 1.5, KW = 4.44×10 -8 m -0.5 yr -1 , and KS = 1.57×10 -8 m -0.5 yr -1 . Figure S18. Contact migration rates measured in our models (dxcontact/dt) vs. kinematic wave speeds (CH) estimated using Eq. 13. Note that this is a version of Fig. 10 that uses all ksn measured over the 10 Myr duration for each simulation (instead of only the final model timestep). Subplot (a) shows results for scenario 3 (contacts dipping upstream, ϕ < 0°) and subplot (c) shows results for scenario 4 (contacts dipping downstream, ϕ > 0°). Dashed lines show the minimum and maximum values for most values in each subplot, with labels denoting the corresponding relationships between observed contact migration rate and estimated kinematic wave speed. Figure S19. Residuals for the three-dimensional regression of EW / U shown in Fig. 11 (n = 1.5, contacts dipping upstream). Figure S20. Variations in the average erosion rate in the weak layer (EW) normalized by rock-uplift rate (U) with both the logarithm of the absolute contact dip in χ-space (log(|ϕχ|)) and the enforced K * (Eq. 7c) for simulations with n = 0.67 and contacts dipping upstream (ϕ < 0º). Note that points are colored by ϕ and have shadows directly beneath them. A gray plane is situated at EW / U values of one to highlight that points at high log(|ϕχ|) values (> -2) generally have EW / U values of about one. A regression is fit to all data (R 2 = 0.64): EW / U = (1.4×10 -3 ln(|ϕχ|) 5 ) + (-2.9×10 -3 ln(|ϕχ|) 4 K * ) + (1.3×10 -2 ln(|ϕχ|) 4 ) + (9.7×10 -4 ln(|ϕχ|) 3 K * 2 ) + (-3.0×10 -2 ln(|ϕχ|) 3 K * ) + (4.8×10 -2 ln(|ϕχ|) 3 ) + (-6.1×10 -4 ln(|ϕχ|) 2 K * 3 ) + (1.3×10 -2 ln(|ϕχ|) 2 K * 2 ) + (-8.2×10 -2 ln(|ϕχ|) 2 K * ) + (1.2×10 -1 ln(|ϕχ|) 2 ) + (-5.9×10 -3 ln(|ϕχ|) K * 3 ) + (9.1×10 -2 ln(|ϕχ|) K * 2 ) + (-3.4×10 -1 ln(|ϕχ|) K * ) + (4.3×10 -1 ln(|ϕχ|)) + (-2.4×10 -3 K * 3 ) + (3.5×10 -2 K * 2 ) + (-1.4×10 -1 K * ) + 1.1×10 0 . Note that points are colored by ϕ and have shadows directly beneath them. The red dashed line represents the erosion rates expected if the contact dip was 0° (Eq. 10). Figure S21. Residuals for the three-dimensional regression of EW / U shown in Fig. S20 (n = 0.67, contacts dipping upstream). Figure S22. Variations in the average erosion rate in the weak layer (EW) normalized by rock-uplift rate (U) with both the logarithm of the absolute contact dip in χ-space (log(|ϕχ|)) and the enforced K * (Eq. 7c) for simulations with n = 1.5 and contacts dipping downstream (ϕ > 0º). A regression is fit to all data (R 2 = 0.32): EW / U = (-6.0×10 -3 ln(|ϕχ|) 4 ) + (-5.6×10 -1 ln(|ϕχ|) 3 K * ) + (1.9×10 -1 ln(|ϕχ|) 3 ) + (3.7×10 0 ln(|ϕχ|) 2 K *2 ) + (-9.9×10 0 ln(|ϕχ|) 2 K * ) + (3.6×10 0 ln(|ϕχ|) 2 ) + (3.0×10 1 ln(|ϕχ|) K *2 ) + (5.9×10 1 K *2 ) + (-5.3×10 1 ln(|ϕχ|) K * ) + (1.8×10 1 ln(|ϕχ|)) + (-8.6×10 1 K * ) + 2.8×10 1 . Note that points are colored by ϕ and have shadows directly beneath them. The red dashed line represents the erosion rates expected if the contact dip was 0° (Eq. 10). Figure S23. Residuals for the three-dimensional regression of EW / U shown in Fig. S22 (n = 1.5, contacts dipping downstream). Figure S24. Variations in the average erosion rate in the weak layer (EW) normalized by rock-uplift rate (U) with both the logarithm of the absolute contact dip in χ-space (log(|ϕχ|)) and the enforced K * (Eq. 7c) for simulations with n = 0.67 and contacts dipping downstream (ϕ > 0º). A regression is fit to all data (R 2 = 0.92): EW / U = (3.3×10 -2 ln(|ϕχ|) 3 ) + (-2.6×10 -2 ln(|ϕχ|) 2 K * ) + (4.5×10 -1 ln(|ϕχ|) 2 ) + (-2.5×10 -1 ln(|ϕχ|) K * ) + (2.0×10 0 ln(|ϕχ|)) + (-2.7×10 -1 K * ) + 3.6×10 0 . Note that points are colored by ϕ and have shadows directly beneath them. The red dashed line represents the erosion rates expected if the contact dip was 0° (Eq. 10).  Figure S26. Variations in the average erosion rate in the strong layer (ES) normalized by rock-uplift rate (U) with both the logarithm of the absolute contact dip in χ-space (ln(|ϕχ|)) and the enforced K * (Eq. 7c) for simulations with n = 1.5 and contacts dipping upstream (ϕ < 0º). Note that symbol size represents the reference weak erodibility (KW), with smaller points corresponding with higher KW values. Also note that the ES / U and ln(|ϕχ|) values here are the mean values taken within logarithmically spaced drainage area bins (e.g., Fig. 9). Points are colored by ϕ and have shadows directly beneath them. The red dashed line represents the erosion rates expected if the contact dip was 0° (Eq. 10). A gray plane is situated at ES / U values of one. Figure S27. Variations in the average erosion rate in the strong layer (ES) normalized by rock-uplift rate (U) with both the logarithm of the absolute contact dip in χ-space (log(|ϕχ|)) and the enforced K * (Eq. 7c) for simulations with n = 0.67 and contacts dipping upstream (ϕ < 0º). Note that points are colored by ϕ and have shadows directly beneath them. The red dashed line represents the erosion rates expected if the contact dip was 0° (Eq. 10). The positions of maximum and minimum ES / U values change with the reference weak erodibility (KW). Figure S27:

Explanation of
The complex patterns in Fig. S27 occur because stretch zones (Royden and Perron, 2013) form when n < 1 and contacts dip upstream. A stretch zone is a gap between slope patches; an example of a stretch zone that many would recognize is the broad, convex-upwards knickzone that sometimes occurs along transient bedrock rivers. Gaps between slope patches occur when n < 1 because the slope patches for higher rates of base level fall ("high-E" slope patches) migrate upstream at a slower rate than the slope patches for lower rates of base level fall ("low-E" slope patches). Higher erosion rates within a weak unit can undercut the strong unit, generating high-E slope patches and a stretch zone. Examples of these dynamics are in Fig. 7c. Near the basin outlet, the strong unit has a high steepness and the weak unit has a low steepness. Moving upstream, a steep reach begins to form at the top of each reach within the weak unit, and reaches within the strong unit have a curved, convex-upwards morphology. These convex reaches are the stretch zones. Moving further upstream, eventually each reach within the weak unit has a high steepness, while each reach within the strong has a low steepness (approaching the behavior for a contact dip of 0°). These longitudinal changes in stream morphology are the physical representation of the trends in Fig. S27; the maximum erosion rates in Fig. S27 (ES / U) occur right before the reaches within the strong unit begin to be dominated by convex stretch zones. The erosion rates decrease with ln(|ϕχ|) past this point because the influence of contact migration becomes less significant at higher contact dips.
The reason the positions of maximum and minimum erosion rates change with the reference weak erodibility (KW) in Fig. S27 is because erodibility influences the turning point at which stretch zones begin to dominate each reach within the strong unit (also the point at which the weak unit begins to have a pronounced section with high steepness). Figure S28. Variations in the average erosion rate in the strong layer (ES) normalized by rock-uplift rate (U) with both the logarithm of the absolute contact dip in χ-space (ln(|ϕχ|)) and the enforced K * (Eq. 7c) for simulations with n = 1.5 and contacts dipping downstream (ϕ > 0º). Note that symbol size represents the reference weak erodibility (KW), with smaller points corresponding with higher KW values. Also note that the ES / U and ln(|ϕχ|) values here are the mean values taken within logarithmically spaced drainage area bins (e.g., Fig. 9). Points are colored by ϕ and have shadows directly beneath them. The red dashed line represents the erosion rates expected if the contact dip was 0° (Eq. 10). A gray plane is situated at ES / U values of one. Figure S29. Variations in the average erosion rate in the strong layer (ES) normalized by rock-uplift rate (U) with both the logarithm of the absolute contact dip in χ-space (ln(|ϕχ|)) and the enforced K * (Eq. 7c) for simulations with n = 0.67 and contacts dipping downstream (ϕ > 0º). Note that symbol size represents the reference weak erodibility (KW), with smaller points corresponding with higher KW values. Also note that the ES / U and ln(|ϕχ|) values here are the mean values taken within logarithmically spaced drainage area bins (e.g., Fig. 9). Points are colored by ϕ and have shadows directly beneath them. The red dashed line represents the erosion rates expected if the contact dip was 0° (Eq. 12). A regression is fit to all data (R 2 = 0.91): ES / U = (-3.7×10 -4 ln(|ϕχ|) 3 ) + (3.4×10 -4 ln(|ϕχ|) 2 K * ) + (7.1×10 -4 ln(|ϕχ|) 2 ) + (-1.4×10 -3 ln(|ϕχ|) K * 2 ) + (1.9×10 -2 ln(|ϕχ|) K * ) + (3.7×10 -3 ln(|ϕχ|)) + (1.8×10 -3 K * 2 ) + (-4.0×10 -2 K * ) + 9.7×10 -1 . Figure S30. Residuals for the three-dimensional regression of ES / U shown in Fig. S29 (n = 0.67, contacts dipping downstream). Figure S31. Comparison of best-fit K values in our numerical models to the weak erodibility (KW) used in each simulation. Note that this is a version of Fig. 12 where we calculate Eq. 13 estimates of kinematic wave speed using all ksn values recorded over the entire 10 Myr duration for each simulation (rather than only the final timestep). (a-b) Χ 2 Misfit Function values for kinematic wave speeds (CH) estimated using Eq. 12, the enforced contact dip (ϕ), and a wide range of K values (200 points spaced logarithmically from 10 -9 to 10 -4 m 1-2n θ yr -1 , where θ = 0.5) relative to the Eq. 13 estimates of CH. (c-d) Comparison between the best-fit K and the KW enforced in the simulations. Subplots (a) and (c) show results for n = 0.67, while subplots (b) and (d) show results for n = 1.5. Figure S32. Comparison of best-fit K values in our numerical models to the strong erodibility (KS) used in each simulation. Note that this is a version of Fig. 13 where we calculate Eq. 13 estimates of kinematic wave speed using all ksn values recorded over the entire 10 Myr duration for each simulation (rather than only the final timestep). (a-b) Χ 2 Misfit Function values for kinematic wave speeds (CH) estimated using Eq. 12, the enforced contact dip (ϕ), and a wide range of K values (200 points spaced logarithmically from 10 -9 to 10 -4 m 1-2n θ yr -1 , where θ = 0.5) relative to the Eq. 13 estimates of CH. (c-d) Comparison between the best-fit K and the KS enforced in the simulations. Subplots (a) and (c) show results for n = 0.67, while subplots (b) and (d) show results for n = 1.5.