River self-organisation inhibits discharge control on waterfall migration

The action of rivers within valleys is fundamentally important in controlling landscape morphology, and how it responds to tectonic or climate change. The response of landscapes to external forcing usually results in sequential changes to river long profiles and the upstream migration of waterfalls. Currently, models of this response assume a relationship between waterfall retreat rate and drainage area at the location of the waterfall. Using an experimental study, we show that this assumption has limited application. Due to a self-regulatory response of channel geometry to higher discharge through increasing channel width, the bed shear stress at the lip of the experimental waterfall remains almost constant, so there was no observed change in the upstream retreat rate despite an order of magnitude increase in discharge. Crucially, however, the strength of the bedrock material exhibits a clear control on the magnitude of the mean retreat rate, highlighting the importance of lithology in setting the rate at which landscapes respond to external forcing. As a result existing numerical models of landscape evolution that simulate the retreat of waterfalls as a function of drainage area with a fixed erodibility constant should be re-evaluated to consider spatial heterogeneity in erodibility and channel self-organisation.

to be studied 2,3 , although the results cannot be directly scaled up to natural rivers in space and 23 time. In particular, processes such as flow acceleration above the knickpoint lip, plunge pool 24 erosion, undercutting of the knickpoint face and channel banks, cantilever failure, erosion and 25 transport of cohesive material by hydraulic shear are present in the experiments (Fig. 3C). The 26 approach has the advantage of being able to vastly reduce the spatial and, therefore, the temporal 27 duration of individual experiments, as well as offering a complete control on boundary conditions. 28 While abrasion by bedload material does not occur in the experiments and would be difficult to 29 reproduce without significantly changing the experimental design, we seize this opportunity to 30 better highlight how knickpoint retreat in an homogeneous cohesive material occurs in response to 31 variations in discharge and lithological strength. It is also important to note that due to the 32 homogenous nature of the material, the impact of some natural phenomenon such as resistant 33 hillslope material is not included within the experiments. 34 Following the 'similarity of process' analogue modelling approach, the cohesive silica paste used 35 as the bedrock material in these experiments is not designed to have material properties that can 36 be directly compared to rocks found in natural environments. However, this does not prevent the 37 exploration of the impact of the relative cohesion of the channel substrate on the rates of fluvial 38 processes in action in the experimental channel (i.e., knickpoint retreat). The knickpoint retreat 39 experiments presented here were performed using three different mixes of two types of silica material: angular silica (54.66 to 65.6% total mass) and spherical silica beads (27.33 to 16.4% total 41 mass), with a constant proportion of water (18% total mass). The impact of the proportion of 42 spherical silica beads on the relative cohesion of the silica paste was tested by placing a 15 x 10 x 43 10 cm slab of material in a tank containing 2.5 cm of water, and monitoring the time taken for the 44 slab to collapse (Fig. S1). A higher proportion of silica beads reduces the relative cohesion of the 45 silica mix, leading to the faster collapse of the slab (Fig. S1C)

50
The 20.5% bead slab after 77 minutes, with 44% of the initial surface collapsed. The proportion of slab 51 collapse was calculated using Digital Elevation Models collected from laser scanner data at 2 minutes 52 intervals. For each pixel of the DEM (pixel size = 2 mm), the pixel was defined to have 'collapsed' when the 53 elevation of the pixel was 1 cm lower than its initial elevation. C. Time taken for a 15 x 10 x 10 cm slab of 54 different silica mixes to collapse after being placed in 2.5 cm of water. From the experiments, the silica mixes 55 containing a greater proportion of silica beads (e.g., 41%) started to collapse earlier, indicating a lower 56 relative cohesion than the silica mixes containing a lower proportion of silica beads.

57
Beyond similarity in hydraulic processes, we note that the relationship between the channel 58 geometry and discharge in these experimental channels is qualitatively similar to natural bedrock 59 rivers (Fig. S2). The relationship between channel width and discharge (W ∝ Q~0 .5 , Fig S2A) is

SI Section 2: Extracting hydraulic information from the experiments 78
To characterise the hydraulics of the experiments, the topography from the laser scanner was 79 coupled with the Floodos hydrodynamic precipiton-based numerical model 6 . All point clouds in 80 each experiment were rasterised with a pixel size of 2 mm, and Floodos was run using the 81 corresponding discharge for each digital elevation model in each experiment. The flow in the 82 experimental channels is laminar in nature with Reynolds numbers typically less than 500 (Table  83 S1), so Floodos was run in the laminar-flow setting. In pure laminar flow, the friction coefficient C is 84 entirely set by the water viscosity μ and should be approximately C~ρg/3μ where ρ is water density 85 and g is gravitational acceleration. This predicts that at 10°C, C ~ 2.5 x 10 6 m -1 .s -1 , and this value 86 has been shown to accurately match the spatial extent of water flow in experimental channels 7 . 87 Performing Floodos simulations is advantageous for this data extraction, as it generates spatial 88 masks of unit discharge and water depth, from which additional hydrological information can be 89 is bed shear stress and S is the water slope (Fig. S3, Table S1). Subsequent calculations of the 92 average value of these parameters for an entire channel mask gives the benefit of a large sample 93 size and the removal of any manual measurement error. Additionally, channel width was extracted 94 from the Floodos output using cross-sections of the wetted area for the entire length of the 95 channel. The length of cross-sections of the wetted area for the reach of the channel unaffected by 96 the input and output conditions (i.e. > 10 cm from each) were extracted from the initial topography 97 (i.e., before a knickpoint was initiated through base level fall). The mean width is shown in Fig. S2, 98 with the error bars indicating one standard deviation.   the Floodos model output. Mean values are provided for all pixels not affected by the inlet or outlet (> 10 cm in distance) and uncertainty is the 108 standard deviation. The shear stress is calculated from a mean of 20 measurements taken one channel width upstream of the knickpoint.

SI Section 4: Relationship between discharge and bed shear stress 112
This section provides the derivation of the relationship between discharge and bed shear stress, 113 allowing the interpretation of the self-organisation of channel geometry and the constant knickpoint 114 retreat rate even under higher discharges. 115 The bed shear stress is calculated using the following equation: 116 where is the bed shear stress, is the fluid density, is the acceleration due to gravity, is the 118 flow depth (m) and is the slope. 119 The discharge Q is calculated according to: 120 where is the flow velocity, is the channel width and is the flow depth. From SI Section 1, we know that ∝~0 .5 and while the relationship between slope and 133 discharge is not consistent between all sets of experiments (Fig. S2) Equations S5.7a and S5.7b show that there is a very weak relationship between shear stress in the 149 channel and discharge, under both laminar and turbulent conditions, consistent with the relatively 150 constant shear stress measured under different discharges in Fig. 4B. Under changing discharge, 151 the channel self-organises through variations in channel width and channel slope. This self-152 organisation prevents higher shear stresses at the knickpoints under higher discharge. As erosion 153 of the silica material is mainly through hydraulic shear, the constancy of shear stress here may 154 explain why knickpoint retreat rate is not higher during the experiments with higher discharges (Fig.  155 3A) as predicted by the stream power incision model. This observation should hold, whether the 156 flow is laminar (as in our experiments, see Table S1) or turbulent. The shear stress is consistent 157 between the experiments using different silica mixes, which could explain the difference in the 158 magnitude of the knickpoint retreat rate. For materials with different cohesions but the same shear 159 stress exerted by the channel, the material with the stronger cohesion should erode slower when 160 erosion is predominantly by shear.

Data extraction and analysis for Roan Plateau 163
The Roan Plateau study area was selected because of the sparse vegetation cover in the area, 164 which allowed the identification of knickpoints and measurement of their width from high-resolution 165 aerial imagery. Knickpoints within 15 catchments were located using Google Earth Pro, and the 166 width of the channel at the knickpoint lip extracted using the measuring tool within Google Earth 167 Pro at each of these locations. The drainage area of the channels at each of these 15 identified