Sediment flux‐driven channel geometry adjustment of bedrock and mixed gravel–bedrock rivers

Sediment supply (Qs) is often overlooked in modelling studies of landscape evolution, despite sediment playing a key role in the physical processes that drive erosion and sedimentation in river channels. Here, we show the direct impact of the supply of coarse‐grained, hard sediment on the geometry of bedrock channels from the Rangitikei River, New Zealand. Channels receiving a coarse bedload sediment supply are systematically (up to an order of magnitude) wider than channels with no bedload sediment input for a given discharge. We also present physical model experiments of a bedrock river channel with a fixed water discharge (1.5 l min−1) under different Qs (between 0 and 20 g l−1) that allow the quantification of the role of sediment in setting the width and slope of channels and the distribution of shear stress within channels. The addition of bedload sediment increases the width, slope and width‐to‐depth ratio of the channels, and increasing sediment loads promote emerging complexity in channel morphology and shear stress distributions. Channels with low Qs are characterized by simple in‐channel morphologies with a uniform distribution of shear stress within the channel while channels with high Qs are characterized by dynamic channels with multiple active threads and a non‐uniform distribution of shear stress. We compare bedrock channel geometries from the Rangitikei and the experiments to alluvial channels and demonstrate that the behaviour is similar, with a transition from single‐thread and uniform channels to multiple threads occurring when bedload sediment is present. In the experimental bedrock channels, this threshold Qs is when the input sediment supply exceeds the transport capacity of the channel. Caution is required when using the channel geometry to reconstruct past environmental conditions or to invert for tectonic uplift rates, because multiple configurations of channel geometry can exist for a given discharge, solely due to input Qs. © 2020 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd


Introduction
The understanding of the processes that set the morphology of bedrock river channels, and the feedbacks and interactions that drive their evolution, are crucial to develop a detailed understanding of rapidly evolving landscapes (Whipple, 2001;Turowski et al., 2006;DiBiase et al., 2015). The elevation of the river channel acts as the local base level for adjoining hillslopes. Vertical and lateral bedrock channel dynamics directly impact the hillslope sediment flux, through processes such as bank erosion and the undermining of hillslopes, and associated geohazards such as landsliding (Golly et al., 2017). In response to variations in climate or tectonic uplift rate, it has been shown that bedrock channels evolve dynamically through adjustments to both width and slope. This was observed across contrasting settings, including Nepal (Lavé and Avouac, 2001), California (Duvall et al., 2004), Italy (Whittaker et al., 2007), Taiwan  and Scotland (Whitbread et al., 2015), as well as in experiments (Turowski et al., 2006) and predicted by mechanistic numerical modelling approaches (e.g. Stark, 2006;Turowski et al., 2007Turowski et al., , 2008Lague, 2010;Nelson and Seminara, 2011;Turowski, 2018;Li et al., 2020). However, landscape evolution models often assume that channel width scales simply with water discharge (Q) or drainage area (e.g. Whipple and Tucker, 1999;Shobe et al., 2017) or with a constant ratio between the channel width and channel depth (e.g. Finnegan et al., 2005), failing to capture dynamic patterns of channel geometry identified in the settings above. This restricts the applicability of these models to account for changes in channel geometry due to climate or tectonic forcing that are important for setting river sediment transport capacity (Q sc ) , bank erosion (e.g. Turowski et al., 2008), channel sinuosity (Stark, 2006;Turowski, 2018) and the undermining or stabilization, in the case of channel narrowing, of adjoining hillslopes (e.g. Golly et al., 2017).
The key factors in setting the pattern and rate of erosion and the resulting morphology of bedrock river channels include the strength of the bedrock material (e.g. Montgomery and Gran, 2001;Yanites et al., 2017;Baynes et al., 2018a), sediment availability or 'tools' for erosion by abrasion (e.g. Sklar and Dietrich, 2004;Turowski et al., 2007;Lamb et al., 2008;Cook et al., 2013Cook et al., , 2014Brocard et al., 2016), the proportion of bedrock exposed on the bed or the 'cover effect' (e.g. Turowski et al., 2007;Chatanantavet and Parker, 2008;Johnson et al., 2009;Lague, 2010), the degree of rock fracturing that promotes erosion by plucking or block toppling (e.g. Lamb and Dietrich, 2009;Dubinski and Wohl, 2013) and the magnitude and frequency of flood events (e.g. Hartshorn et al., 2002;Baynes et al., 2015). Of these factors, recent work has focused on the processes driving vertical incision and the resulting channel geometries by sediment impacts, including the relative importance of suspended and bedload sediment (Scheingross et al., 2014) and the role of bed roughness and sediment flux (Finnegan et al., 2007;Whipple, 2007, 2010). It has been observed, first by Gilbert (1877) and many times subsequently in the field (e.g. Sklar and Dietrich, 2001;Turowski et al., 2007;Johnson et al., 2009;Whitbread et al., 2015) and in laboratory experiments (e.g. Baynes et al., 2018b), that high sediment supply relative to river transport capacity can protect the bed from vertical incision , with potentially increased lateral erosion rates as a result of sediment impacts being more focused on the channel banks (Turowski et al., 2008;Beer et al., 2017). Both the availability and nature (i.e. grain size) of bedload sediment that acts as 'tools' are important for controlling rates of erosion through abrasion in controlled experiments (e.g. Sklar and Dietrich, 2001), but the related impact on channel morphology remains relatively understudied. Finnegan et al. (2017) identified the impact of grain size and sediment supply on steady-state bedrock channel slopes in a small (30km 2 ) catchment in the Santa Cruz Mountains (California, USA). Tributary channels are three times steeper for a given drainage area where they transport coarse bedload but the main channel widens, shallows and shifts to more complete alluvial cover downstream of where it starts to transport coarse bedload material . In fully alluvial channels, increasing sediment supply (Q s ) is crucial in setting the channel geometry and planform morphology, including an increase in width above that predicted by threshold theory (Lacey, 1930;Glover and Florey, 1951;Parker et al., 2007;Phillips and Jerolmack, 2016) and the transition to a braided river pattern at high Q s compared to Q sc (see review in Métivier et al., 2017). However, there remain outstanding questions on the link between Q s , vertical and lateral erosion rate and the channel geometry (width and slope) in bedrock and mixed bedrock-gravel bed rivers across multiple scales Lague, 2014;Finnegan et al., 2017;Turowski, 2018) and how they respond to changes in Q s (Turowski and Hodge, 2017;Yanites, 2018).
This study presents a combined field and experimental study to explore outstanding questions related to (i) the relationship between bedrock and mixed gravel-bedrock channel morphology (width and slope) and bedload sediment availability and (ii) the morphodynamic processes that drive bedrock and mixed gravel-bedrock river channel evolution under different sediment supply regimes. To answer these questions, we first explore trends in channel geometry of river channels in the Rangitikei catchment (~3000km 2 ) of the North Island of New Zealand, a natural laboratory where the impact of bedload sediment availability on channel geometry can be isolated and the nature of the bedload measured. We then use a physical model coupled with a numerical hydrodynamic model to explore the processes and interactions between sediment dynamics and channel incision in a controlled laboratory setting, allowing the processes leading to small changes in channel geometry to be identified and quantified. In particular, we focus on the impact of bed cover on lateral erosion and the changing distribution of shear stress within channels under different sediment supply regimes.

Geomorphological setting
The channel geometry dataset presented here is from the Rangitikei River and surrounding catchments located in the central southern part of the North Island of New Zealand ( Figure 1). The rivers are, at present, incised bedrock and mixed gravel-bedrock channels, with a geomorphological history characterized by climatically controlled cycles of aggradation and incision. The modern-day trunk channel of the Rangitikei is deeply incised into a strath terrace, topped by fluvial deposits, that was abandoned following the Last Glacial Maximum (LGM) (Milne, 1973;Pillans et al., 1993;Litchfield and Rieser, 2005;Bonnet et al., 2019). In the middle reach of the Rangitikei River, the dating of recent terraces formed during post-LGM incision indicates a recent (post-10 kyear) incision rate of 1.5-2.0mmyear À1 (Bonnet et al., 2019).
Importantly, for the purposes of this study, much of the trunk valley network is incised into weak Plio-Pleistocene marine mudstones to sandstones of the Whanganui Basin (Proust et al., 2005), but the lithology of the headwater regions has a strong spatial variation ( Figure 1). The eastern headwaters drain the Ruahine mountain range, supplying metamorphosed Triassic to early Cretaceous greywacke material as coarse bedload sediment to the river network while, in contrast, western headwaters have a very limited supply of such material due to the lithological composition of these areas ( Figure 1). Therefore, due to the uniform weak lithology in the trunk channels, which are fed by tributaries both with and without a supply of coarse bedload sediment, the Rangitikei catchment is a natural laboratory where the role of bedload sediment supply in setting channel geometry can be assessed directly.

Methods
Channel geometry and drainage area Measurements of the active channel width were collected from across the catchment using DigitalGlobe imagery in Google Earth (resolution 0.5m). Measurements were collected from the trunk channel and tributaries both with and without a supply of coarse-grained bedload material from headwaters containing greywacke, to cover the full variability in channel geometry across the Rangitikei catchment ( Figure 1). All measurements of width were taken from locations where the local lithology is either sandstone or mudstone ( Figure 1) and were measured orthogonally to the river flow direction at distances of approximately 1km. Measurements were made, as far as possible, from straight channel reaches rather than the apex 3715 SEDIMENT FLUX-DRIVEN CHANNEL GEOMETRY ADJUSTMENT of meander bends to minimize the natural variability in the dataset. Active width was defined in the images as the distance orthogonal to the flow direction between either bedrock banks (i.e. cliffs) or the limit of vegetated bars on the edge of the channels, to avoid uncertainty introduced by using the extent of the water coverage on the imagery. Where channels were multithreaded, the width measurement did not include the width of vegetated islands between channel threads. To test the level of natural variability in the channel width, 50 measurements were taken along a 500m reach at Mangaweka, upstream of a large meander bend ( Figure 2F; 39.80932°S, 175.80811°E). The mean channel width from these measurements was 62.8m with a standard deviation of 3.2m, giving an uncertainty of approximately 5% in the width measurements as a result of the natural variability of the channel. Uncertainty related to human measurement error and the resolution of the Google Earth imagery was estimated by measuring the channel width in the same location within this reach 20 times, with a mean measurement of 63.8m and a standard deviation of 1.2 m (~1.8%).
From each of the width measurements, the upstream drainage area was extracted using a 90m-resolution digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM) using the r.water.outlet tool in QGIS version 2.18. The lithology of the upstream drainage area was extracted for each of the width measurements using the geological map (Lee et al., 2011;Heron, 2014) and the percentage of the catchment draining basement greywacke rock used as a proxy for the supply of hard coarse-grained material.
Rainfall and river discharge Drainage area at each location was converted to mean daily discharge, to remove any regional orographic rainfall control in the relationship between channel width and drainage area. Based on the difference between average annual rainfall from four stations at high elevations in the Ruahine range and nine stations in the main catchment (see Table S1 in the online supporting information for annual rainfall data from locations across the catchment; provided by Horizons Regional Council), we identified an orographic rainfall pattern with a ratio of 1:1.54 between the western edge of the study area and the eastern limit of the measurements on the western side of the Ruahine range. To convert the drainage area to mean daily discharge (Q m ) for the western catchments, we assumed the orographic rainfall increased at a constant rate from west to east and used the following equation: where R is the value of the orographic rainfall value between 1 and 1.54 as a function of the relative distance from the measurement location between the western edge of the study and the location of peak rainfall, A is the catchment drainage area and K is the scaling constant between the drainage area and the mean daily discharge calculated using the discharge record from the gauging station at Mangaweka (39.80932°S, 175.80811°E): For the three catchments located to the east of the Ruahine, in the absence of any rainfall data, we assumed a decrease in the orographic rainfall pattern of 1.54:1 from the peak rainfall location used above to the eastern limit of the study area and calculated the mean discharge for each location according to the relative position along this orographic gradient.

Sediment grain size distribution
The catchment-wide channel geometry and catchment lithology data were supplemented by field measurements of sediment grain size at locations across the catchment where sediment was stored within the channel and easily accessible for measurement (Figure 1). A random walk method of sampling the a-axis of 125 grains along a 20m transect was performed at each location (Wohl et al., 1996). Multiple transects at the same location were collected where clear spatial differences in grain size were apparent on active bars within the channel compared to bars on the channel banks (see Section S2 in the online supporting information for details of the number of transects per site, and photos of the sites).

Field data
Channel width The presence of hard, coarse-grained, bedload sediment is clear within the channels that drain part of the Ruahine mountain range (Figure 2). These channels are wider than channels that do not receive a supply of coarse-grained bedload material from the headwaters (Figures 2A-C). The width of the Kawhatau River ( Figure 2D) is an order of magnitude greater than the Makohine River (Figure 2A), despite having a drainage area only three times greater and flowing through similar weak lithology.
Analysis of the catchment-wide channel width dataset ( Figure 3A; Table 1) demonstrates that the channel width increases with discharge, consistent with observations from Figure 3. (A) Channel width against drainage area and discharge, for sites located within mudstone or sandstone, coloured by proportion of catchment-draining greywacke. There is a clear distinction in the channel width between catchments that drain greywacke (grey/black) and those that do not drain greywacke (open symbols). (B) Channel width normalized by the square root of discharge. For a given discharge, the channels that contain coarse sediment have a higher normalized channel width than those that do not contain sediment. (C) Relationship between normalized channel width and the proportion of the upstream drainage area that drains greywacke. Data plotted is the mean normalized channel width (error bars indicate ±2 standard deviations) within bins of 5% range.
other natural settings (e.g. Lague, 2014). There is also a clear relationship between channel width and whether the catchment drains areas of greywacke lithology, used as a proxy here for coarse-grained bedload sediment. Channels that do not drain any greywacke (in white in Figure 3A) are systematically narrower than channels that do have a supply of coarse bedload sediment from the headwaters (shades of grey/black in Figure 3A) across two orders of magnitude variation in mean daily discharge. Normalization of the channel width by Q 0.5 (the common scaling for bedrock channels reported in the literature; see Lague, 2014) further demonstrates a dependency between the channel width and the proportion of the catchment supplying coarse sediment ( Figure 3B). Further analysis of the normalized channel width according to the proportion of the upstream catchment that drains greywacke ( Figure 3C) also highlights the role of bedload sediment supply on the channel morphology, with the normalized channel width broadly increasing in a linear relationship (r 2 = 0.48, p-value = 0.009). Channels without a bedload supply have very low normalized channel widths (~0.01), while channels with a supply of coarse bedload sediment have normalized channel widths ranging from 0.02 to 0.09. The relationship between the coarse bedload sediment supply and the channel width is most prominent in small tributaries close to the Ruahine range front where a larger upstream proportion of the catchment is greywacke (mean daily discharge between 0.5 and 10 m 3 s À1 ) compared to the trunk channel where the discharge is largest (>50 m 3 s À1 ) and the normalized channel widths are generally lower ( Figures 3B and C).

Grain size
We have investigated how the grain size varies with transport distance rather than the catchment drainage area due to the spatially variable nature of the bedload sediment supply into the headwater tributaries. For example, the drainage area of the Rangitikei (stars in Figure 4A) increases threefold between the upstream and downstream locations (Table S1), but there is no coincident increase in supplied sediment volume as the majority of the additional tributaries do not drain the source of the greywacke sediment; the Ruahine mountain range ( Figure 1). The transport distance is therefore a more representative parameter in this setting for exploring sediment grain size variability and any possible impact of the grain size on the channel geometry within the study catchments. Transport distance was measured as the longitudinal distance between the edge of the greywacke areas to the grain size measurement site ( Figure 1). Where multiple tributaries drain the greywacke (e.g. in the Rangitikei, downstream of the confluence with the Kawhatau; Figure 1), we used the shortest transport distance to the greywacke contact but the transport distance from the greywacke contact along the trunk channel is also provided in Table S2.  Other than a single outlying transect from the Pohangina River, there is no clear trend in downstream fining of the median (D 50 ) or coarsest (D 90 ) sediment grain size fractions over a transport distance of up to 70km from the Ruahine mountain range. D 90 is consistently between~100 and 300 mm ( Figure 4). The outlying transect of the Pohangina (D 90 = 510mm) was measured from the outside of a meander bend near a bridge, and D 90 from another transect measured from a gravel bar in the middle of the active channel at the same location (40.05237°S, 175.93847°E) was 153mm. While it wasn't the only survey site with a nearby bridge, we suggest the local conditions of the outlying transect on the outer meander bend immediately downstream of the bridge pile have led to a non-representative grain size from this location and we do not consider this transect further. Given the relatively consistent grain size across the catchments located close (<25km) to the Ruahine range front, which typically have low mean daily discharge (Table S2), we suggest that variations in the sediment grain size cannot solely be the cause of the large variability in normalized channel width (between 0.02 and 0.09; Figure 4B) for these rivers.
In the trunk channel of the Rangitikei, there is no clear downstream fining trend in the coarse fraction of the sediment grain size distribution ( Figure 4A). The normalized channel widths in the Rangitikei are typically lower than the catchments close to the mountain range front, despite a similar proportion of the catchment draining greywacke (normalized channel width 0.01-0.02; Figure 4B), possibly the result of the slightly smaller grain size in the trunk channel. However, the trend in the sediment grain size data is not strong enough to suggest that the variability in the normalized channel width with discharge is solely the result of sediment grain size variability within the channels ( Figure 4B).

Physical Modelling of the Impact of Sediment Flux on Channel Geometry
The results from the Rangitikei River clearly demonstrate that bedload sediment flux has a significant impact on the channel geometry at the catchment scale but is limited by the lack of data showing the dynamics of channel geometry variability at the reach scale. Here, we complement the data from the Rangitikei with a suite of controlled laboratory experiments performed at the reach scale, allowing additional data to be collected related to the impact of sediment flux on the hydraulics within the channel (e.g. shear stress) and also monitoring of the evolution of the channel slope to take place.

Experimental approach
Flume setup Analogue physical model experiments were performed in the 80×30cm Bedrock River Experimental Incision Tank at Université de Rennes 1 ( Figure 5; see Turowski et al., 2005Turowski et al., , 2006Turowski, 2007;Baynes et al., 2018a, b for more information on experimental apparatus). While no formal scaling of the experiments with a particular natural location was sought, the similarity of the processes represented in the analogue model allows the exploration of the impact of different boundary conditions on channel morphology and the evolution of the system over spatial and temporal scales that would otherwise be impossible using strict scaling approaches (Hooke, 1968;Paola et al., 2009;Baynes et al., 2018c). The experiments were designed to isolate the role of sediment flux (Q s ) on bedrock channel geometry in a tectonically quiescent setting and identify the processes that explain the observations of channel geometry in natural rivers. Seven experiments (Table 2) were performed under a constant water discharge using a mix of homogenous cohesive silica material (grain size = 45μm) with uniform strength (Baynes et al., 2018a). During an initial channel adjustment phase, the channel outlet was dropped by 2cm to generate a knickpoint that formed an incised channel with cohesive 'bedrock' banks rather than unchannelized flow over a flat surface. Following this initial base level fall to channelize the flow, the elevation of the channel outlet was fixed throughout the rest of the experiment. Channel geometry has been shown to be a function of the tectonic uplift rate (e.g. Lavé and Avouac, 2001;Duvall et al., 2004;Whittaker et al., 2007), but is not considered in the experiments presented here as the channel outlet elevation was fixed following the initial channelization phase. The experiments therefore represent a different type of tectonic setting to many existing studies, where the role of tectonics in controlling channel geometry is a main focus (e.g. Turowski et al., 2006). To represent bedload sediment material, coarse sand (grain size = 250μm) was added at the channel inlet using a Gericke Infinite Screw feeder system. The sand was collected beyond the channel outlet to avoid any deposition in the channel as a result of the channel outlet configuration ( Figure 5). The channel width and slope were free to adjust in response to Q s through lateral erosion of the banks, or deposition of material or erosion of the channel bed when exposed to clear water flow. Channel geometry evolution during the experiments was monitored using a Leica ScanStation 2 terrestrial laser scanner, programmed to collect high-resolution (sub-millimetre) point clouds of the topography every 2min. The ScanStation 2 uses a green laser that penetrates water. The point clouds were gridded to generate DEMs at 2mm resolution, and the resulting DEMs were used as input for the Floodos numerical hydrodynamic model . The Floodos model generates a water depth map over the topography for a given discharge, which was used to automatically extract channel geometry information as well as calculate hydraulic information such as bed shear stress or Froude number (see Baynes et al., 2018b for further description of the Floodos model parameter calibration for predicting hydraulic information in the experimental channels).
Relevance of laboratory experiments to the Rangitikei catchment Direct comparisons of the laboratory results to particular locations within the Rangitikei study area are limited. Instead, the experiment is designed to replicate the key behaviour of the target system, allowing the spatial and temporal scales to be reduced, and allowing greater flexibility in a range of parameters to be tested. There is a key difference between the dominant erosion process in the Rangitikei (abrasion) and the experiments (hydraulic shear by clear water flow), negating the possible comparison of the 'tools effect' between the experimental results and the natural setting. However, we suggest that the other key ingredients of the experiments replicate the processes in the Rangitikei well: a weak bedrock material that can easily be eroded by hydraulic shear (aided in the Rangitikei by wetting-drying cycles that weather the bedrock) and the coarse cover of sediment that can protect the channel bed from vertical incision, the 'cover effect' (Sklar and Dietrich, 2004;Turowski et al., 2007;Chatanantavet and Parker, 2008;Johnson et al., 2009;Lague, 2010). Despite the lack of the 'tools effect' in the experiments, we suggest that the important process for studying the effect of sediment supply of bedrock channel geometry is the cover effect, with channel widening occurring when the bed is protected from vertical erosion and flow is directed towards the channel walls (e.g. Turowski, 2018). We therefore acknowledge the limitations of the experimental approach for the lack of the 'tools effect' but suggest that the similarity of the other processes, particularly the presence of the 'cover effect', allows insights to be drawn on the role of sediment supply in setting the width and slope of bedrock channels.

Experimental results
Evolution of channel geometry Channel width evolution through time during the experiments is shown in Figure 6, with the initial phase at the start of the experiment where the channelization of the flow took place highlighted in red. All subsequent reporting of results does not include this initial phase when the channel geometry is affected by the incisional nature of the channel as well as the input Q s . It should be noted that the end of this initial phase does not necessarily indicate the time when the channel geometry has stabilized, but rather when the initial phase of channelization and generation of a bedrock-constrained channel has been completed (i.e. the knickpoint has reached the channel inlet and the channel is fully channelized). At this time, the channel width is still adjusting to the balance of the input Q s , including the reduction in supply of material eroded from the bed (i.e. by the retreating knickpoint). The fluctuations in the channel width measurements can be estimated through the error bars in Figure 6, which indicate two standard deviations from the mean value of the width of the cross-sections from each of the DEMs. For channels with low Q s ( Figures 6A-C), the natural variability in the channel is low. The channels become more dynamic with increasing Q s (e.g. Figures 6D-G), due to the presence of a migrating channel thalweg and the generation of small terraces. Figure 7 provides summary statistics for the channel geometry for all of the DEMs during each experiment, excluding the channel initialization phase (see Figure 6). For example, the violin plot for Q s = 0gl À1 contains all of the channel geometry measurements after the initial incision phase (t= 92 min) until the end of the experiment (t = 240min). The data was extracted from the whole channel reach, with the exception of the upper and lower 20cm of the flume to negate any possible impact of the channel inlet and channel outlet. The width (W), slope (S) and width-to-depth ratio (W/D) of the experimental channels varies with the input Q s (Figure 7), with channels typically wider and steeper under higher Q s . Q s appears to have the largest impact on the channel slope, and the spread with the plots in Figure 7 indicates the degree of variability in W, S and W/D under the different sediment regimes. The median channel width and slope values at low Q s are similar, but the wider spread of the channel widths at Q s = 2.5 and 5gl À1 compared to Q s = 0gl À1 indicates the impact on channel morphology. However, some of the extreme values that make the tails of the plots may be indicators of the state of the channel shortly after the initial channelization phase when the channel has not necessarily fully stabilized ( Figure 6).
Channel form at representative conditions In addition to the impact on mean W, S and W/D, Q s also has an important role in setting the morphology and characteristic dynamics of the channels. To directly compare the results of typical channel geometry under the different Q s regimes, we extracted data from individual DEMs that are representative of, or as close as possible to, channels that have fully evolved to the input Q s . We define a channel as being 'representative' at the time when the channel geometry has had plenty of time to evolve to the input conditions and is not undergoing widespread changes (Figure 6). The timing of the DEMs that we selected for each Q s regime is shown by the blue lines in Figure 6. Typically, these 'representative' conditions are at the end of the experiment, giving the full 240min for the channel width and slope to stabilize to the input Q s . The mean channel width data for each DEM presented in Figure 6 indicate the temporal evolution of the channel through time and the standard deviation error bars indicate the variability in the channel width for each DEM. Typically, this variability is reduced as the channel stabilizes through time (e.g. Figures 6A and B), but there is some variability in the width throughout the length of the experiments for the higher sediment fluxes. Under these conditions, we selected the DEM closest to the end of the experiment that appeared visually stable and did not contain any DEM artefacts generated during the interpolation process.
Given this continued small degree of variability, we refer to these as 'representative' channels rather than 'steady state' to reflect the potential for not having reached perfectly stable conditions, and also to avoid confusion with the common usage of the term 'steady state' to refer to the balance between vertical erosion and uplift rate. Additionally, we chose to consider all the cross-sections for each DEM to demonstrate the natural variability in the channel width, rather than the temporal evolution of a single cross-section which may not fully represent the channel morphology at any given time. However, we suggest that using the channel geometries from these extracted 'representative' times allows a better comparison between Q s regimes than using data extracted from the entire time period, including immediately after the initialization phase ( Figure 6). Due to these limitations discussed above, we exercise caution when directly comparing the experimental results to the Rangitikei natural example, but we suggest that the experiments are still informative, showing the relative impact of Q s on the channel geometry for a fixed discharge within the same system (from Q s = 0to 20gl À1 ). The experiments with a low Q s (0 and 2.5gl À1 ) were typically characterized by a single active channel, with a stable channel bed and partial bed cover by sediment for 2.5 and 5 gl À1 (Figures 8A-H; Figure S8). The channels at higher Q s deposited material on the whole channel bed, increasing the elevation of the channel bed and adjusting the channel slope (Figure 7). It is important to note that during the Q s = 20gl À1 experiment, the banks of cohesive 'bedrock' were buried by sediment deposited on the channel bed ( Figures 8M and N). This led to dynamic channels that migrated across the width of the flume, sometimes reaching the edges (e.g. Figure 8M) and typically behaving more like an alluvial system rather than a bedrock or mixed gravel-bedrock channel. This experiment therefore represents a comparison dataset between the dynamics of an 'alluvial' channel and bedrock or mixed gravelbedrock channels in the experimental setup.
Bed shear stress In addition to the water depth mask (Figure 8), the Floodos hydrodynamic model output  was used to calculate the bed shear stress for each pixel within the channels, τ =ρgDS, where τ is bed shear stress, ρ is density of water (1000kgm À2 ), D is water depth and S is hydraulic slope. The distribution of shear stress values within the channel at 'representative' conditions for each of the Q s regimes is shown in Figure 9. Typically, the mean and median τ within the channel across the different Q s regimes are relatively similar, with all values above the critical shear stress (τ c ) (see below). For experiments where Q s > 0gl À1 , the range of τ values was greater than for the Q s = 0gl À1 experiment, with a higher peak shear stress (τ max ) but also more values closer to 0Pa (Figure 9). From Figure 9, there is no apparent link between increasing Q s and increasing τ for individual pixels within the channel, but due to variability in channel width (Figures 7 and 8) this does not necessarily mean that there is no increase in total sediment transport potential with higher Q s . The critical shear stress (τ c ) was calculated using the following equation: where τ Ã c is the Shield's criterion (typically 0.03-0.06; Buffington and Montgomery, 1997), ρ s is the density of the sediment (2750kgm À3 ), ρ w is the density of water (1000kg m À3 ), g is the acceleration due to gravity (9.81ms À2 ) and D 50 is the median grain size of the sediment (250μm). Using these values, we calculate τ c to be~0.15-0.25Pa, approximately three to four times smaller than the mean and median τ (Figure 9), which is a similar ratio to many natural gravel-bedded rivers such as those located in western North America . Following Meyer-Peter and Muller (1948) and Lajeunesse et al. (2010), we compute the 'excess bedload shear stress' (τ bl ) for each pixel during the representative width conditions to explore the potential for sediment transport within the channel: We calculated τ bl for each pixel within the studied reach of the experiments, and extracted the cumulative τ bl across 20 channel cross-sections at equal intervals (18mm) from the studied reach (see Figure 8 for excluded areas at channel inlet and outlet). A higher cumulative τ bl within the channel indicates a greater potential for sediment transport, as there are more pixels where τ > τ c . The mean total τ bl (across 20 cross-sections) increases with Q s when Q s > 5gl À1 (Figure 10) under a range of τ c values (0.1-0.3), to take into account a range of values for the Shield's criterion discussed above. Where Q s ≤ 5gl À1 , the total τ bl is approximately the same for all of the experiments.
Sediment flux, bed shear stress and river lateral mobility Importantly, despite an increase in τ bl with increasing Q s , the distribution of τ on the channel bed does not follow the same pattern for each Q s regime (Figures 9 and 11). For each experiment, the two columns of Figure 11A show examples of two DEMs, with the 20 cross-sections plotted next to each other (same cross-sections used for the τ bl calculation in Figure 10). The normalized shear stress at a single point (τ/τ max ) is coloured at each location along the cross-section. Probability density functions (pdf) of τ/τ max across the cross-sections are shown in Figure 11B for the first two DEMs as well as three additional DEMs. Where τ is close to being uniformly distributed across the width of the channel, τ/τ max in a cross-section would all be close to 1 (darker colours in Figure 11A and the pdf show that most values lie towards high values of τ/τ max ). Where τ is distributed unevenly, there would be sections of the cross-section where τ/τ max is close to 1, and other areas where τ/τ max is close to zero (lighter colours in Figure 11A and a pdf dominated by values that lie close to low values of τ/τ max ).
Where Q s = 0gl À1 , τ is distributed uniformly across the whole width of the channel, and the planform morphology is stable through time (e.g. the top row of Figure 11 and the small range of values in Figure 9). When Q s > 0gl À1 , the distribution of τ/τ max becomes increasingly complex with a network of 'active' channel threads that contain higher τ than other parts of the cross-section, shown by the pattern of light and dark patterns in the bottom panels of Figure 11. Additionally, the pdfs show that with an input Q s , the shape of the curve shifts such that more values lie towards a low value of τ/τ max . At Q s = 0gl À1 , there is a unimodal peak with a relatively narrow standard deviation centred at τ/τ max~0 .75, while at Q s = 2.5, 6.66 and 7.5gl À1 a broader peak is centred at τ/τ max0 .5. At Q s = 10, 20 and, to some extent, 5gl À1 , there is a general decay in the frequency of high τ/τ max , indicating that the majority of the channel does not have τ close to τ max . At Q s = 20gl À1 , where the channel is alluvial in nature, the distribution of τ is dynamic and evolves quickly through time (see the difference in the spatial map of τ/τ max in the bottom panels of Figure 11 showing DEMs that are separated by 2min). Therefore, under high Q s , the channel morphology evolves to accommodate the extra material but the sediment transport does not occur at a uniform rate across the whole width of the channel. A similar behaviour can be seen, to a lesser extent, for Q s = 7.5 and 10gl À1 , suggesting that mixed bedrock-gravel channels, with cohesive 'bedrock' banks, are also dynamic systems that evolve their morphology according to the input Q s .

Discussion
Comparison to theoretical predictions of channel geometry Constant width-to-depth ratio Finnegan et al. (2005) proposed, on the basis of Manning's formula, that river width is a function of the lithology into which the channel is developed, with a constant ratio between the width and the depth for a given substrate. The experimental results presented here (Figure 7) demonstrate that W/D can be highly variable, as a function of Q s . Therefore, we suggest that while the use of constant W/D for the modelling of bedrock channel width gives a first-order estimate of the channel geometry through the consideration of the role of discharge, this simple approach fails to capture the impact of the complexities of erosion-sedimentation processes. In particular, the experiments highlight the role of sediment cover on the bed, favouring erosion of the channel banks that, in turn, drives increased rates of lateral erosion and channel widening (e.g. Turowski et al., 2009;Turowski, 2018Turowski, , 2020Li et al., 2020). This behaviour is consistent with observations from Boulder Creek, California, USA, where the main channel was wider and shallower downstream of the introduction of coarse bedload material .
Threshold theory According to Lacey's law (Lacey, 1930;Glover and Florey, 1951;Métivier et al., 2017), when the conditions within a channel are just above the threshold for entrainment of the sediment, the channel geometry (width and slope) is set by the balance between gravity and fluid friction, and the width scales with the square root of the discharge: where μ is the friction angle, θ t is the threshold Shield's param- eter, C f is the turbulent friction coefficient, K 1 2 is a transcendental integral with value of~1.86 and Q Ã ¼ Q= ffiffiffiffiffiffiffiffiffiffiffi gD 5 50 q (Glover and Florey, 1951;Métivier et al., 2017). The width of channels normalized by the sediment grain size in the Rangitikei is typically larger than the width predicted by the Lacey relationship for channels at the threshold of entrainment for the corresponding discharge and sediment grain size ( Figure 12). Note that in this analysis we use three different scaling factors between bankfull discharge (Q bf ) and mean daily discharge (Q m ). In Section S3 of the online supporting information we estimate, based on analysis of gauge data and a survey of river stage, that Q bf =4.5Q m ( Figure 12B) for the Rangitikei but we provide the other plots to demonstrate the impact of the Q bf scaling factor on Q * . For the experimental channel, the grain size used (D 50 ), when Q s = 0, 2.5 and 5g l À1 , where the bed was not or only partially covered, was 45μm (for the silica 'bedrock' material). For all other experiments, the grain size used was 250μm (for the sand). The choice of D 50 is important for the calculation of Q * and W/D 50 and is why the channels with low Q s have larger values of Q * and W/D 50 in Figure 12 than channels where Q s > 5g l À1 . Despite this, it is still possible to compare the relative  Figure 11. Visualization of the distribution of normalized shear stress within the channel. For each sediment flux, 20 cross-sections were extracted from each of five DEMs which were collected at 2-min intervals. Normalized shear stress is calculated as the value of shear stress relative to the peak shear stress along each cross-section. The DEMs represent the five closest to 'representative' conditions. (A) The 20 cross-sections (normalized by width) for two of the DEMs (Time 1 and Time 2) for each of the sediment fluxes. Only two of the DEMs are plotted to save space, with the data from all five DEMs plotted in the pdfs in (B). Where the majority of the points are coloured black, with little variability along a cross-section (e.g. when Q s = 0gl À1 ), the distribution of shear stress is relatively uniform. Where there is a large variability in the colours across a cross-section (e.g. Q s = 20gl À1 ), the distribution is not uniform within the channel and multiple narrower channels are active across the channel bed, which can be mostly inactive. In the case of Q s = 20gl À1 , the location of these active channels is dynamic, with the pattern of distributed shear stress changing between the DEMs difference between observed W/D 50 for the experimental channels with the predictions of Lacey's law ( Figure 12). Using an extensive compilation of previous experiments of alluvial channels, Métivier et al. (2017) found a relationship of increasing channel width Q s before a transition in channel morphology from single-thread to braided river conditions occurs. In this compilation, single-thread channels typically have a lower Q s , and the width is closer to the threshold width predicted by Lacey's law (Figure 12). A similar relationship is found in the experimental data in this study and data from Baynes' et al. (2018a) experiments without Q s but at variable discharge. Experiments in Baynes et al. (2018a) are close to the predicted width at threshold for entrainment conditions, while the experiments in this study with an input Q s plot further from threshold theory ( Figure 12). Differences in the channel geometry from a clear single thread to a more variable channel planform with increasing Q s are also observed in the Rangitikei (Figure 3), in other field locations that contain strath terraces where narrower channels have incised within wider braid plains following reductions in sediment supply (e.g. Finnegan and Balco, 2015) and in our experiments (Figures 8 and 11), although a threshold Q s where this behaviour begins is hard to identify beyond being in the region of 5-7.5gl À1 . Importantly, despite the experiments in this study being performed with cohesive bed material, they follow the same pattern of a transition to different planform morphologies and distribution of shear stress at very high Q s .
The sensitivity of Q * to the Q bf scaling factor controls how aligned the Rangitikei channels are to the compilation dataset of Métivier et al. (2017) and threshold theory, but does not impact the pattern within the Rangitikei dataset, where channels with a higher bedload sediment supply typically plot further from the threshold theory and can experience a change in the channel behaviour (e.g. Figure 2). We therefore suggest that a Q s -dependent channel form is not just a feature of alluvial rivers typically considered in the compilation of Métivier et al. (2017), but also rivers constrained by banks made of relatively weak bedrock material (e.g. the mudstone of the Rangitikei or the experimental substrate).
Adaptation of channel geometry and bed shear stress in response to sediment flux Figure 13 demonstrates the range of possible channel geometry configurations (i.e. W, S) at a constant discharge under different Q s . For comparison, we have included a reference dataset from similar experiments, where Q s = 0 but under different discharges (from Baynes et al., 2018a). For a constant discharge, the impact of Q s on the width and slope of channels is greater than the impact of changing the discharge when Q s = 0gl À1 , with the variability in the Q s dataset larger than the reference  Stebbings, 1963;Schumm et al., 1987;Ikeda et al., 1988;Sapozhnikov and Foufoula-Georgiou, 1996;Warburton, 1996;Métivier and Meunier, 2003;Tal and Paola, 2007;Peakall et al., 2007;Brauderick et al., 2009;Dijk et al., 2012;Ashmore, 2013;Reitz et al., 2014). Black solid line is the threshold theory calculated when θ t = 0.05, C f = 0.1, as used by Métivier et al. (2017). Grey dashed lines indicate the uncertainty in the calculation of the threshold width due to uncertainty in the values used for the critical shear stress and friction coefficient. Upper grey dashed line calculated using θ t = 0.03, C f = 0.3; lower grey dashed line calculated using θ t = 0.3, C f = 0.02. Panels A-C indicate a sensitivity analysis based on different values of the scaling factor to calculate the bankfull discharge (Q bf ) from the mean daily discharge (Q m ). [Colour figure can be viewed at wileyonlinelibrary.com] 3727 SEDIMENT FLUX-DRIVEN CHANNEL GEOMETRY ADJUSTMENT dataset that covers an order of magnitude of discharges ( Figure 13). Here we explore the processes that could explain the geometric adaptation to these different Q s scenarios, linking back to the Rangitikei.
When Q s is low, the characteristic bedrock channel conditions are detachment-limited and bedrock erosion processes dominate the evolution of the channel geometry (Lague, 2010). For an incising channel under detachment-limited conditions, the channel geometry is set by the incision rate and independent of Q s (Lague, 2010). In our experimental channels where incision rate I = 0 and Q s = 0, the channel geometry is dependent on τ c and there is a uniform distribution of τ bl within the channel (Figures 10 and 11). A small increase of Q s (2.5 and 5gl À1 ) does not seem to drive an increase in τ bl (Figure 10), suggesting that the channel transport capacity (Q sc ) was higher than or equal to 5gl À1 .
For Q s > 5gl À1 , the channels widen and steepen to increase their transport capacity. By definition, if a channel changes its geometry to accommodate an increase in sediment supply, it is no longer under detachment-limited conditions (e.g. Lague, 2010). The channel can then evolve towards strict transport-limited conditions (Q s = Q sc ) or hybrid conditions (Q s < Q sc ) (Lague, 2010) for which the geometry is set by a combination of the upstream Q s and local incision rate. Because I = 0 in our experimental setup and the discharge is constant, the fraction of the bed covered by sediment is an indirect indication of these conditions: if the channel bed is fully covered by sediment, and in particular with a thick immobile sediment cover, it is under transport-limited conditions. If the bed is only partially covered with immobile sediment, this means that the local Q s is equal to or larger than the local Q sc , but the complete channel section could be in hybrid or transport-limited conditions (Q s ≤ Q sc ). This is specific to the experimental configuration for which discharge is constant. In natural rivers, the analysis of the long-term conditions of mixed bedrock-alluvial channels based on the present-day sediment cover at low flow is more difficult owing to stochastic fluctuations in discharge and sediment supply (Lague, 2010). At 'representative width' conditions, the experimental channel was fully covered with sediment when Q s ≥ 6.66gl À1 , with partial bed coverage for Q s = 2.5 and 5gl À1 (Figure S8), suggesting that it is very likely that the experimental channel transitions to transport-limited conditions when Q s > 5gl À1 .
Under constant discharge conditions, wider channels would decrease τ as the flow is spread over a larger wetted area and thus has a reduced water depth. However, the coincident increase in the channel slope (Figures 7 and 13), as well as the morphodynamic transition towards an actively migrating braided system of narrow threads with high local τ (Figure 11), contributes to the overall increase of channel transport capacity. This behaviour of the experimental bedrock-gravel channels appears consistent with the behaviour and characteristic morphology of alluvial rivers. The presence of an input Q s leads to a larger offset from the theoretical predictions of Lacey's law ( Figure 12) and the transition from single-thread to multiple-threaded channels after reaching a threshold sediment supply and shift to transport limited conditions (Q s > 5g l À1 ).

Geometry of the rivers in the Rangitikei region
The scaling between width and discharge for most of the channels in the study area (Table 1) follows the common relationship found in bedrock rivers across multiple settings (W = kQ 0.5 , where k is a constant; Lague, 2014). However, three of the channels (Hautapu, Pohangina and Makaroro) are notable exceptions due to a negative exponent. The Pohangina ( Figure S6) and Makaroro both drain the Ruahine mountain range, are heavily laden with coarse sediment, and in-channel bars and multiple active threads are evident. We suggest that the negative exponent in these channels may be related to the high sediment supply or local variabilities in lithological strength, especially in the Lower Pohangina, where the bedrock includes a thick succession of coarse terrigeneous clastic deposits (Rees et al., 2017). Further detailed field study of these channels would highlight these controls. The Hautapu contains no coarse sediment and close to the confluence with the Rangitikei trunk channel, the valley is more incised and narrower ( Figure 2C) than upstream ( Figure 2B). We therefore suggest that the general narrowing of the channel with discharge in the Hautapu is related to the response of the tributary to vertical incision in the Rangitikei, with the incision wave yet to fully migrate upstream. Despite these three exceptional channels, it is clear that the availability of bedload sediment has an important role in the bedrock channel width in the Rangitikei by affecting the magnitude of k in the scaling between width and discharge (Table 1; Figures 2 and 3).
There exist two components to the 'tools' effect: (i) sediment availability and (ii) the nature of that sediment (i.e. grain size or lithology; Sklar and Dietrich, 2001). In the Rangitikei, we suggest that sediment availability is the dominant component of the 'tools' effect as the sediment grain size is relatively consistent throughout the channels (Figure 4). We propose that in the channels with a high availability of coarse-grained bedload, the 'cover' effect protects the bed from erosion and drives an increased volume of impacts on the channel banks, driving lateral erosion and widening the channel.
The transition from detachment-limited to transport-limited conditions and different behaviour of the channel system due to bedload sediment flux is evident in the example from the Rangitikei. Tributaries that do not receive a supply Figure 13. Mean channel width against channel slope closest to 'representative width' for (i) the experiments in this study (squares), where sediment flux is variable but all points with a constant discharge of 1.5lmin À1 and (ii) a reference dataset where Q s = 0gl À1 under different discharges between 0.5 and 3lmin À1 (white circles). Labels indicate the value of Q for each of these points (data from Baynes et al., 2018a). The channel width for the 20gl À1 experiment was potentially limited by the width of the experimental box due to the unconstrained banks and dynamic migrating channel, and the true width may be larger than plotted (indicated by *). The channel at 2.5gl À1 remained constrained within banks of silica material and thus was not limited by the experimental box. Mean values of width and slope are from the five DEMs leading up to the representative conditions, with the error bars indicating ±2 standard deviations from the mean. 3728 E. R. C. BAYNES ET AL.
of coarse-grained bedload material represent the detachment-limited channel configuration for a given discharge (Figures 2 and 3), while any departure from this geometry due to coarse-grained bedload sediment supply represents a transport-limited or hybrid channel configuration. The channels that are not supplied with coarse-grained bedload material from greywacke areas plot close to the threshold theory ( Figure 12), similar to the Q s = 0gl À1 experiments, while the tributaries that drain a high proportion of greywacke plot further from the threshold theory, consistent with channels that have a high bedload sediment flux ( Figure 12). These channels also actively migrate laterally (Lague et al., 2013;Bonnet et al., 2019) as a meandering bedrock channel, and sometimes contain multiple active threads (Figure 2), whereas the channels with low or no bedload sediment supply are more uniform in nature, consistent with the distribution of shear stress within the experimental channels where Q s = 0gl À1 (Figure 11). We therefore suggest that the findings from the experiments, where channels become less uniform in terms of the value and distribution of shear stress after a transition to transport-limited or hybrid conditions, are applicable to natural environments such as the Rangitikei region.

Conclusion
Sediment is a critical component of bedrock channel systems, with the supply of hard, coarse-grained material as bedload having a direct impact on the width of channels in the Rangitikei River, New Zealand. The presence of a bedload sediment supply acts to protect the channel bed from vertical erosion, increasing impact on the channel banks and leading to increased channel widths for a given river discharge. Under experimental conditions, varying sediment supply leads to a range of possible configurations of channel width and slope for a given discharge. Increased sediment supply leads to a complex and evolving pattern of shear stress within the channel as the pattern of flow becomes more dynamic, while maintaining a constant mean shear stress just above the critical shear stress for the motion of sediment, replicating a similar effect identified in laboratory and natural alluvial channels. We therefore suggest that sediment flux should be an important parameter to consider when using preserved bedrock channel geometries to reconstruct past conditions (e.g. palaeodischarge), or estimating uplift rates from channel geometry, as a range of possible configurations may be possible for a given discharge. The feedback between sediment supply and channel geometry should also be included in numerical models of landscape evolution. Further work is required to identify the response of bedrock channel geometries to variable sediment supply through time. Table S1: Rainfall data provided by Horizons Regional Council, used in calculation of orographic rainfall pattern. In the station name column, (R) indicates a station located in the Ruahine and (MC) indicates a station located in the main catchment. Table S2: Sediment grain size information, including site location, site characteristics and grain size data. Figure S1: Sediment grain size measurement site on the Kawhatau River. Figure S2: Sediment grain size measurement site on the Kiwitea River. Figure S3: Sediment grain size measurement site on the Mangawharariki River. Figure S4: Sediment grain size measurement site on the Mangoira River. Figure S5: Selected sediment grain size measurement sites on the Oroua River. Figure S6: Selected sediment grain size measurement sites on the Pohangina River. Figure S7: Selected sediment grain size measurement sites on the Rangitikei River. Figure S8: Images from the mounted camera at 'representative width' conditions for each of the experimental conditions, with the exception of Q s = 10gl À1 where the camera failed, and no images were collected. 3731 SEDIMENT FLUX-DRIVEN CHANNEL GEOMETRY ADJUSTMENT