Influence of Boulders on Channel Width and Slope: Field Data and Theory

Large boulders with a diameter of up to several tens of meters are globally observed in mountainous bedrock channel environments. Recent theories suggest that high concentrations of boulders are associated with changes in channel morphology. However, data are scarce and ambiguous, and process-related studies are limited. Here we present data from the Liwu River, Taiwan, showing that channel width and slope increase with boulder concentration. We apply two mass balance principles of bedrock erosion and sediment transport and develop a theory to explain the steepening and widening trends. Five mechanisms are considered and compared to the field data. The cover effect by immobile boulders is found to have no influence on channel width. Channel width can partially be explained by boulder control on the tools effect and on the partitioning of the flow shear stress. However, none of the mechanisms we explored can adequately explain the scattered width data, potentially indicating a long-timescale adjustment of channel width to boulder input. Steepening can be best described by assuming a reduction of sediment transport efficiency with boulder concentration. We find that boulders represent a significant perturbation to the fluvial landscape. Channels tend to adjust to this perturbation leading to a new morphology that differs from boulder-free channels. The general approach presented here can be further expanded to explore the role of other boulder-related processes.

efficiency with boulder concentration. We find that boulders represent a significant perturbation to the 23 fluvial landscape. Channels tend to adjust to this perturbation leading to a new morphology that differs 24 from boulder-free channels. The general approach presented here can be further expanded to explore the 25 role of other boulder-related processes. Here, we explore the hypothesis that increasing boulder concentration can cause channels to widen and 89 steepen. We collected field data of channel geometry, morphology, and boulder characteristics. The data 90 are based on field and remote sensing observations from the Liwu River in Taiwan, where boulders with 91 diameters of up to 25 meters are ubiquitous on the channel bed along various river reaches, differing in 92 drainage area and geometry. Our goals are to (1) study the relationship between channel slope, width, and 93 boulder concentration and (2) theoretically identify the processes linking these variables. To achieve the 94 second goal, we focus on the influence of boulders on bedrock erosion and sediment transport and describe 95 their effects using conceptual arguments and empirical relations. We establish a suite of mechanisms that 96 relate immobile boulders to bedrock channels' steady-state width and slope. We treat each mechanism 97 independently, discuss the trends it predicts, and compare them against field observations. In the following, 98 we review the literature concerning flow hydraulics induced by boulders, the effects of boulders on 99 Various studies reported links between boulders and channel morphology (e.g., Montgomery and 128  channel steepness (Thaler and Covington, 2016) and slope (Shobe et al., 2020) in bedrock channels. Steady-131 state bedrock channel morphology is assumed to control long-term river bedrock erosion adjustment to the 132 long-term uplift rate (e.g., Whipple and Tucker, 1999). When the uplift rate changes, the erosion rate 133 responds by adjusting the river profile (e.g., Whipple, 2004;Lague et al., 2005) and cross-section geometry 134 (e.g., Turowski et al., 2009a;Yanites, 2018;Turowski, 2020). Assuming immobility of large boulders, 135 Shobe et al. (2020) exploited this notion to argue that boulders hinder erosion by protecting the bedrock 136 channel bed. Their model predicts that a boulder-bed channel would consequently steepen to compensate 137 for the reduced erosion. Immobile boulders were also argued to be consequential for changes in bedrock 138 channel width. Shobe et al. (2020) tested the influence of the proximity of the boulder delivery point (e.g., 139 landslides scars) on the width coefficient, i.e., width normalized by drainage area (e.g., Lague, 2014) and 140 found contrasting results. Accordingly, conclusive data and a general theory of boulder influence on 141 bedrock channel width are still missing. 142 Understanding how boulders influence channel morphology in bedrock rivers requires insights into the 143 process of bedrock erosion and sediment transport. The slope of bedrock channels has been argued to adjust 144 to both bedrock erosion requirement and the mobilization of upstream sediment supply (e.g., Sklar and 145 Dietrich, 2006). However, the degree to which slope adjusts to each of these components remains unclear 146 (Johnson et al., 2009). While channel slope is commonly considered to be the consequence of bedrock 147 erosion and reshaping of the long profile (e.g., Royden and Perron, 2013), recent studies suggest that 148 equilibrium of bedrock channels could be attained by a modification of the slope of sediment overlying the 149 bedrock (Phillips and Jerolmack, 2016;Turowski, 2020Turowski, , 2021. As in alluvial channels, rearrangement of 150 the bed to form a new sediment-bed slope can be achieved via selective deposition and entrainment 151 processes during floods (Mackin, 1948;Schumm and Parker, 1973;Schneider et al., 2015b;Turowski and 152 Hodge, 2017), which relates to sediment transport processes. Furthermore, adjusting sediment-bed slope 153 can be achieved within a timescale of a single flood, significantly faster than the timescale associated with 154 bedrock erosion and the formation of a new bedrock slope (Turowski, 2020). 155 In abrasion-dominated channels, erosion of the bedrock bed and banks are thought to occur during flood 156 events and are driven by impacts of sediment grains, which travel as bedload (e.g., Sklar and Dietrich, 2004;157 Cook et al., 2013;Auel et al., 2017). Channel widening occurs by lateral erosion, which is thought to be a 158 consequence of sediment particles deflected to the sides following encounters with bed roughness elements 159 (e.g., Li et al., 2020). A field study from a bedrock channel gorge in Switzerland showed that wall erosion 160 increases in proximity to roughness elements (Beer et al., 2017). Although recent studies proposed a 161 positive relationship between channel roughness and lateral erosion (Fuller et al., 2016;Turowski, 2018;162 Li et al., 2020; He et al., 2021) the precise nature of this relationship remains to be explored (Turowski, 163 2020). 164 In the light of the above review, adjustment of channel width and slope to perturbations caused by 165 immobile boulders can be expected to be controlled by bedrock erosion and sediment transport processes. 166 Changing channel slope could occur by eroding the bedrock bed and altering sediment cover and depth by 167 sediment deposition and entrainment. In contrast, existing models and observations indicate that widening 168 the channel is only possible by lateral erosion. Due to the estimated long timescales of width adjustment 169 (see below) (Turowski, 2020), a link between steady-state width and boulder concentration can be 170 established if we consider at least one of the following conditions. First, the timescale of bedrock channel 171 width adjustment to boulder input is shorter than the residence time of boulders within a river. Theoretically, 172 the widening of bedrock channels such as the Liwu River is expected to extend to periods of up to thousands 173 of years. Second, boulder supply and boulder degradation balance each other to keep the concentration of 174 boulders steady over the required time scale for width adjustment. These assumptions will be reviewed in 175 the discussion. The Liwu River, Taiwan, exhibits numerous fluvial bedrock reaches hosting huge boulders. This section 184 describes the methods applied for data collection in the Liwu River (Section 2.1) and empirical relations of 185 boulder concentration and channel slope and width based on these data (Section 2.2). 186

Data Collection 187
We documented 20 fluvial reaches along the Liwu River. Field data were collected in field campaigns 188 during the low flow seasons of 2018 and 2019. We selected different fluvial reaches with variable drainage 189 areas and local relief, representing various portions of the drainage basin. Our primary focus was on reaches 190 with a substantial number of boulders, but we also collected data from reaches with lower boulder 191 concentrations. We avoided fluvial reaches with incoming tributaries to ensure a minimal difference in 192 drainage area within a given reach. To avoid lithology differences, we used a geological map of Taiwan to 193 verify that the lithology does not change within the selected reaches. We avoided reaches exhibiting a large 194 spatial variability in channel width. In each channel reach, a drone was used to document the channel at 80 195 -120 m above the channel, constrained by the complexity of the topography and the pilot's location. The 196 channel bed and banks were photographed primarily at vertical and various other angles, with ~80% 197 overlapping area. Due to the steep topography of many bedrock canyon sections, most of the reaches were 198 inaccessible by foot, thus prohibiting emplacement of Ground Control Points (GCPs). We generated point 199 clouds from the photos by using the AGISOFT METASHAPE commercial software. We created 200 orthophotos and DEMs at 5 -25 cm/pixel horizontal spatial resolutions, depending on the site and data 201 quality. To account for the elevation uncertainty associated with the output models, we assume an elevation 202 error of ± 0.5 m for the DEM. 203 The reach area Atot was manually delineated using a digitization process in ArcGIS. First, the upstream 204 and downstream channel reach boundaries were chosen and delineated with straight lines, bounding what 205 we observed as a continuous distribution of boulders (Fig. 2). Second, the channel bank boundaries were 206 identified and tracked by following distinctive bedrock-vegetation contacts. To evaluate the boulder-207 concentration in the channel reach, we manually digitized the map-view area of all of the visible boulders 208 with an average diameter ≥ 2 m (Fig. 2B). A boulder was commonly recognized by observing that it 209 protrudes from either water or a gravel bar. Boulder-concentration was calculated using the relation Γ = 210 Ab/Atot, where Ab is the sum of the areas of all of the boulders and can range between zero and one. To 211 extract boulder diameters from the delineated map-view polygons, we assumed that boulders are circles. 212 Reach-averaged channel width Wb was calculated using two methods: (1) by dividing the reach area by the 213 thalweg length L, the assumed streamwise distance that follows the curvature of the map-view channel 214 banks, and we consider a 5 m uncertainty on the measurement of L. (2) By manually measuring ten bank-215 to-bank lengths, perpendicular to channel banks, along the reach and using the average as a representative. 216 We calculated the Root Mean Square Error (RMSE) value between the two methods to be 3.4 m (5% of the 217 average reach width measurements in all reaches). The standard deviation (STD) of each ten measurements 218 was used as an error on channel width measurements. To calculate reach-scale channel slope Sb in a boulder-219 bed channel, the cross-sections that define the upstream and downstream boundaries of the reach area 220 ( Fig. 2) were extracted from the DEM. The minimum elevation point of each cross-section was subtracted 221 and divided by L. Because a substantial fraction of the bedrock bed is occupied with fine sediments, the 222 slope represents a sediment-bed slope, which might differ from the bedrock-bed slope. Theory and global observations show that to first order and in the absence of other perturbations, 232 channel width increases (e.g., Montgomery and Gran, 2001;Whitbread et al., 2015) and slope decreases 233 (e.g., Wobus et al., 2006) with drainage area. Consequently, to isolate the effects of boulders, the impact of 234 the drainage area needs to be removed. We use a dimensionless width ratio Wb/W, defined as the average 235 width ratio of paired reaches located immediately upstream or downstream from one another, presumably 236 without tributaries joining between them. Wb refers to the measured width of a boulder-bed channel, and W 237 refers to the width of a boulder-free channel. Similarly, we define a slope ratio Sb/S. Each selected pair of 238 reaches shares a similar drainage area and lithology. Calculations of boulder-free width W were performed 239 using two approaches. First, the average of ten measurements exploiting Google-Earth imagery, and 240 second, utilizing a basin-wide scaling relationship between width and drainage area for boulder-free 241 channels (Fig. S2). The first approach can test local width anomalies compared to the standard width 242 derived using the second approach. Data points that show significant difference between the two methods 243 are suspected of experiencing a local effect on width. Channel reaches with a large discrepancy between 244 the two measurements are marked differently in plots. The boulder-free slope was determined by a power-245 law regression between channel slope versus drainage area (

Empirical Relations 249
The collected data include channel reaches with widths ranging between 30 and 120 m, slopes ranging from 250 0.01 to over 0.08, and boulder concentrations that range between ~0 and 0.34 (Table 1 boulder concentration. In both cases, normalization by using the paired boulder-free reach improves the 257 relationship with Γ, as indicated by an increase of the R 2 (Fig. 3). Although the width ratio exhibits scatter 258 for a given Γ, Wb/W is always larger than one for Γ > 0.05. Over a range of 35% variability in boulder-259 concentration, the slope ratio increases from about unity to > 4.

275
In this section, we develop a theoretical framework that yields steady-state analytic solutions that predict 276 the width ratio Wb/W and the slope ratio Sb/S as functions of boulder-concentration, Γ. The geometrical 277 adjustment of a boulder-bed bedrock channel is associated with two aspects of its mass balance. First, 278 Bedrock Rivers evolve by matching their erosion rates to the applied uplift rates. Second, like alluvial rivers 279 (e.g., Mackin, 1948), it has been argued that bedrock rivers evolve to achieve a graded state (Turowski,280 2020), related to the mass balance of river sediments. Aggradation of the bed occurs if the flow is unable 281 to carry the supplied sediments from upstream. Conversely, degradation of the bed arises when the flow's 282 ability to mobilize sediments is larger than the sediment supply. When the channel geometry reflects a 283 condition where the power of the channel to mobilize sediments exactly equals upstream sediment supply, 284 the channel is considered graded. When boulders disturb a bedrock channel, the channel responds by 285 altering its geometry until an erosion-uplift balance and grade are met again. As is shown below, the 286 solutions developed under the erosional balance assumption predict Wb/W as a function of Γ, while 287 predictions derived from the grade assumption yield solutions involve both Wb/W and Sb/S. 288 289

Influence of Boulders on Bedrock Erosion 290
The erosion rate in abrasion-dominated bedrock rivers is thought to be physically driven by the impacts 291 of moving sediment grains during floods (e.g., Dietrich, 1998, 2004;Turowski et al., 2007). 292 When sediment flux increases, more sediments are available to impact the channel bed, causing erosion 293 and contributing to the so-called 'tools effect' (e.g., Cook To account for the effect of immobile boulders in Eq. (1), Cf is defined as the sediment cover due to mobile 303 fine grains only and does not include the cover by large immobile boulders. The fine cover Cf can be 304 calculated from a cross-sectional perspective, by dividing the width which is not covered by sediments 305 ('uncovered width'), with the total width W. To predict steady-state channel width using Eq. (1) we need an 306 assumption about the cover Cf at steady-state. Turowski (2018) suggested that steady-state width can be 307 related to a length scale d [m], which indicates the distance in which a sediment particle is deflected to the 308 side after impacting a roughness element, thereby causing bedrock wall erosion. Bedload deflected towards 309 the sidewalls can cause wall erosion if d is larger than the cover-free channel width. In contrast, no wall 310 erosion occurs when d is smaller than the cover-free width. At some point, the channel width adjusts such 311 that particles almost arrive at the channel wall but do not cause erosion (Turowski, 2018(Turowski, , 2020. In this 312 specific steady-state, d is equal to the uncovered width 313 Substituting (2) into (1) and solving for the width, the model predicts steady-state width to be: 315 The sideward deflection length d is expected to vary in space and time and can expected to depend on 317 channel hydraulics, roughness, and sediment supply (Fuller et  larger than those of fine grains. Here, we assume a cross-section configuration with an immobile boulder 331 in the center and a patch of fine cover that hugs one of the banks (Fig. 4). In contrast to previous works 332 (e.g., Sklar and Dietrich, 2004), we define the riverbed fraction covered by mobile sediments as Cf = Af / 333 (Atot -Ab), where Atot is the reach area, and Af and Ab are the areas covered by fine sediments, and boulders, 334 respectively. The total cover Ctot of the mobile sediments and immobile boulders is then written as 335

Equation (4) can be combined with the equation of vertical erosion rate (1) by replacing (1 -Cf) with (1 -339
Ctot), where both Cf and Γ ranges between zero and one. To illustrate this choice, when Γ is 0.5, half of the 340 channel reach area is covered by boulders, and half is free to accommodate non-stationary, finer sediments. 341 Then, Cf may be adjusted according to the remaining proportion, e.g., Cf = 1 means that the fine sediments 342 cover the remaining bed area, half of the total reach area. In this case, the fine steady-state cover can be 343 described using the definition for the fine sediment cover Cf 344

The Tools Effect 364
In a boulder-bed channel reach, immobile boulders occupy a fraction of the total bed area, thus reducing 365 the bed area exposed to erosion. We assume that impacting sediments acting as erosion tools can 366 concentrate on such reduced exposed bedrock patches. Consequently, for a given cross-sectional geometry, 367 the existence of immobile boulders increases bedload flux per unit exposed (or reduced) width. width portions, including boulders (Wcb), exposed bedrock (WBR), fine cover (Wcf). Here we assume that the particle deflection db equals the exposed bedrock.
According to Eq. (8), for db/d = 1, boulder-bed width increases with boulder-concentration due to the tools 380 effect. The combination of both tools and cover effects into a single model yields 381 In this case, the solution collapses to Eq. (6), for α = 0, or indicates an increase of width with boulder-383 concentration for α > 0. Here, cDB [m] is the length of a typical pile of clustered boulders (Fig. 5), and c (> 1) is a dimensionless 395 parameter, assumed to equal the number of boulders constituting the boulder-island in the cross-section 396 direction. Bedload is considered to be evenly distributed between the in-channels, such that in each of them, 397 the average bedload flux _ ̅̅̅̅̅̅ is given by 398 We assume that steady-state cover adjusts within each in-channel independently so that deflected sediments 400 arrive precisely at the boulder or channel bank but do not cause lateral erosion. In this case, each in-channel 401 width Wic can be approximated using a form of Eq. (3) 402 The total reach width Wb MCE is the sum of all in-channel widths ̅̅̅̅ and the width occupied by boulders 404

Influence on Sediment Transport 434
The concept of grade (Gilbert, 1877;Davis, 1902;Mackin, 1948) stipulates that a channel removed 435 from equilibrium adjusts its system variables to restore the ability to transport the same sediment supplied 436 from upstream. A central paradigm is that channel slope adjusts to achieve grade (e.g., Mackin, 1948 Equation (17) can be replaced by a discharge-based equation for sediment transport (Rickenmann, 2001), 460 which takes the form of (e.g., Turowski, 2021) 461 Here, Q is water discharge, and KBL is a constant describing transport efficiency. The exponent m typically 463 takes values between 1 and 4 ( Barry et al., 2004), while n ranges between 1.5 to 2 (Rickenmann, 2001). 464 The exponent q sets the dependence of bedload transport on channel width and is often assumed to be equal 465 to zero (e.g., Rickenmann, 2001). However, given the unsteady nature of bedload transport and along-466 stream variations in channel width (Cook et al., 2020), the parameter q may differ from zero. Analytically 467 derived end-member approximations have been discussed by Turowski (2021), which give values for q of 468 zero, 0.1, or 2.5. 469 The influence of boulders on sediment transport can be considered via the boulder effects on the various 470 parameters in equations (17)  They showed that fractional transport efficiency K ' BL/KBL, where K ' BL is the reduced transport efficiency 497 coefficient due to roughness, decreases with boulder-concentration. Using digitization of their bedload data 498 (their Fig. 8e) Here, Sb STE /S and Wb STE /W are dependent variables, whereas Γ is independent and q, n, θ, and ν are empirical 506 parameters. Closing equation (21) requires that the width ratio is substituted with either one of the models 507 derived in section 3.1 or supplied with field data. 508 509

The effect of Shear-Stress Partitioning 510
The total shear stress acting on a channel boundary is commonly used as a first-order parameter for 511 prediction of bedload fluxes (e.g., eq. 18; Einstein, 1950 where the geometry is simplified and roughness is considered to be steady. Natural bedrock channels often 514 exhibit bedforms and large grains, which act as obstacles to the flow, altering water velocity gradients and 515 associated shear stresses. Mainly, roughness elements bear a fraction of the total available shear-stress τ, 516 thus decreasing the available shear stress for entrainment of bedload τm. Einstein and Banks (1950) 517 suggested that the total resistance to roughness elements equals the sum of the resistance of each of the 518 individual components. This partitioning approach for transport predictions was further developed for 519 immobile boulders (Yager et al., 2007). We adopt this approach to predict channel width and slope in 520 boulder-bed channels, acknowledging that boulders are roughness elements. Here, Fm = τmAm is the resisting force due to the roughness of the channel bed without boulders, which 527 encompasses both skin friction and drag (Dey, 2014), Fd = τdAd is the resisting force due to drag on boulders, 528 and Atot, Ad and Am are the channel areas upon which the forces are applied, respectively. The skin friction 529 component due to boulders Fs is assumed to be negligible. To facilitate area calculations, we can divide Eq. 530 (22) by the total reach area Atot to obtain 531 In a large flood, the entire bed is submerged, and the mobile area Am upon which drag applies is proportional 533 to the overhead projection area without boulders, i.e., Am/Atot = (1-Γ). However, boulders extend into the 534 flow; thus, drag forces act mostly on their upstream sides and Ad = nDB 2 , with n being the number of 535 boulders in the reach. Using the definitions for boulder-concentration Γ = nDB 2 /WL and for the reach area 536 Atot = WL, we introduce Ad/Atot = Γ. Thus, (23) can be rewritten as: 537 Considering that boulders reduce the total shear stress, we aim to find an expression for the reduced shear 539 stress, τm/τ, which we assume is responsible for fine sediment transport. First, the fractional boulder-drag If only the effect of shear stress partitioning is considered, then the combination of (16) and (17)  Where δ is an exponent relating water velocity to the hydraulic radius Rh and equals ½ for a Darcy-558 Weisbach relation or 2/3 for a Manning-Strickler relation. For δ equals ½, the right-hand side of (29) equals 559 one, and the boulder-bed channel slope Sb equals the boulder-free channel slope, whereas when δ equals 560 2/3, the slope ratio Sb/S depends on the expression on the right-hand side of (29) to the power of 1/7. With 561 such a low exponent, the effect of shear-stress partitioning on the slope ratio is expected to be small and is 562 not likely to reproduce the data. 563 564

571
The mechanisms introduced in Section 4 yield five equations for the width and two for the slope ratio (Table  572 2). These can be tested against the Liwu River data (Fig. 3). Each model contains various parameters, some 573 of which could not be independently constrained. Due to the scatter in the width ratio versus boulder-574 concentration (Fig. 3C), we do not expect a single set of parameters to predict the entire width ratio dataset. 575 Moreover, plotting a single model with specific parameter values requires calibration against field data, 576 which will bias the model towards good performance. Instead, we analyze slope and width ratio sensitivity 577 to parameter changes within the different models by plotting model results for different parameter scenarios 578 while holding other parameters constant. 579 580

The Tools, Cover, 'Tools and Cover,' and Multi-Channel Effect Models 581
The cover (Eq. 6), the tools (Eq. 8), and the combined 'tools and cover' (Eq. 9) effects can be solved 582 explicitly and can therefore be directly compared to the field data. We first note that the cover model is 583 independent of boulder concentration (Fig. 6). Its prediction does not follow the trend observed in the field 584 data and is equivalent to a case where Wb = W. This trivial model yields an RMSE value of 0.48, which 585 forms a benchmark error to which the various width models are compared (Table 2). 586 The tools and the 'tools and cover' effects contain one free parameter, α, which could vary between zero 587 and one, and db/d = 1 is assumed throughout the analysis. In the case of α = 0, both models yield the trivial 588 result of Wb = W. Therefore, we turn to present the results of the two models using α = 1. The tools effect 589 predicts an increase in the width ratio with boulder-concentration (Fig. 6), with a model-data RMSE of 0.40 590 (Table 2). Although a lower RMSE value than the trivial model, the tools effect underpredicts most of the 591 data. The 'tools and cover' model predicts an increase in the width ratio similar to the tools model. Still, it 592 predicts an even smaller exponent and underpredicts the field data, yielding an RMSE value of 0.40. 593 The multi-channel effect, Eq. (14), is implicit for the boulder-free width W, which also appears on the right-594 hand side of the equation. Therefore, to compare the model to data, we assign the measured values of 595 boulder-free width W (using the Google-Earth derived channel width; see Fig. S2), mean boulder diameter 596 DB, and boulder-concentration for each data point. The parameter c can be interpreted as the number of 597 boulders constituting a boulder pile along the cross-section (Fig. 5). Here, to represent the uncertainty in c, 598 we assume that it takes values between 1 and 4. The multi-channel model plots relatively close to the field 599 data (Fig. 7A). Considering the error on Wb/W and the uncertainty in c, the model accounts for 75% of the 600 Liwu width ratio data. The RMSE between the field and best-fit model width ratios is 0.51, which is larger 601 than the prediction for Wb = W. three models: (I) the tools effect (Eq. (8); orange shaded area depicts the model output range when the parameter α is varied between zero and one), (II) the 'tools and cover' effect (Eq. 9; red shaded area, using the same range for α range as in the tools effect), and (III) the cover effect (Eq. 6; black dashed line). The tools and the 'tools and cover' models predict that the width ratio to increase as a response to an increase in boulder-concentration, while the cover effect is constant. All models are plotted using db/d = 1. The notation for the white circles is given in the caption of

Reduction in Sediment Transport Efficiency 628
The reduction of sediment transport efficiency model, equation (21), combines both the width and 629 slope ratios and therefore requires a second equation or independent data to close the system. Furthermore, 630 to solve equation (21), the parameters q, n, θ, and ν need to be constrained. The parameter q was shown to 631 take end-member values of 0, 0.1, 1, and 5/2 (Section 4.2.1). We study the behavior of q on the model since 632 its appearance in Eq. (21) implies a covariant effect of channel width and channel slope. For each q value 633 explored, we iterated and chose random values of the remaining unknown parameters: n, θ, and ν from a 634 specified range of values and selected those that minimized the RMSE value between the model output and 635 Figure 7: The width ratio Wb/W versus boulder-concentration Γ compared between the Liwu river field data and the multi-channel effect model. The width ratio data is plotted versus boulder concentration (blue and white circles). Each gray bar represents the application of the model (Eq. 14) on a single field data point by using the specific boulder-free width of that boulder-bed channel reach. The vertical range of the gray bar represents uncertainty in the parameter c, which is varied between 1 and 4. The notation for the white circles is given in the caption of Fig. 3. the Liwu data. Using the tools model to replace Wb/W and α = 1, when q is low (i.e., equals zero or 0.1), 636 the model captures the increase in Sb/S with Γ (Fig. 8). In contrast, for larger values of q, the model deviates 637 significantly from the data. The model performs best (lowest RMSE) when q is close to or set to zero (Table  638 2), which corresponds to a case where the slope ratio Sb/S is independent of the width ratio Wb/W. exploring a range of β values that fit the width ratio data. Then, β = 1.25 was held constant while Γmax values 661 were varied to study their role in controlling model behavior. 662 Exploring the model, we find that it predicts a non-monotonic trend. At small boulder concentrations, the 663 width ratio is predicted to increase, then it reaches a maximum, after which it predicts a decrease in width 664 ratio with increasing boulder concentration (Fig. 10). For a given Γmax, larger β shifts the width ratio maxima 665 and magnitude towards larger Γ values and larger Wb/W values, respectively (Fig. 9A). A similar behavior 666 is observed when increasing Γmax (Fig.9B). A model with β = 1.0 captures 60% of the width ratio data 667 within one STD error, so does. To test whether the data set can be described by a non-monotonic model, 668 we evaluated Spearman's rank correlation coefficient between the width ratio and boulder-concentration. A 669 calculated value of 0.65 implies that the two variables are positively correlated (for comparison, the rank 670 correlation coefficient between the slope ratio and Γ is 0.76). However, a non-monotonic relationship 671 cannot be ruled out. 672 Considering the effect of the shear stress partitioning on the slope ratio, in section 3.2.2, we showed that 673 the slope ratio depends on boulder concentration to a maximum power of 1/7. Regardless of the choice of 674 the other free parameters, this produces only a weak dependency between the slope ratio and boulder 675 concentration, which makes the shear stress partitioning model inadequate to describe the Liwu slope ratio 676 data (Fig. 10).

Reviewing the Assumptions of Steady-States 712
In our theoretical framework, we have assumed steady-state and tested the resulting equations using 713 field data from the Liwu River. Among the examined models, some have produced higher goodness-of-fit 714 values (e.g., the reduction of transport efficiency effect), while others showed a certain degree of 715 incompatibility compared to the data (Section 4; Table 2), thus requiring an assessment of the applicability 716 of the steady-state assumptions. 717

Steady-State in the Erosional Balance Assumption 718
Under the steady-state assumption in the erosional balance, we assume that (i) neighboring boulder-bed 719 and boulder-free reach incise at the same rate. A substantial difference in incision rates between two 720 adjacent channel reaches would promote a knickpoint between the two reaches. We have not observed any 721 such prominent knickpoint at any of the 20 studied sites. (ii) The channel width is at steady-state with 722 respect to erosion rate and fine cover (Turowski, 2018). This assumption is valid if boulders are present at 723 the same reach-averaged concentration at a particular location for a sufficiently long time. In numerous 724 reaches that we examined, there is direct evidence for a continuous supply of large boulders (Text S3). 725 Hillslopes nearby boulder-bed channels often exhibit scars typical of landslides and rockfalls. However, 726 whether those boulders were delivered to the Liwu tributaries recently or if they were placed a long time 727 ago requires further research. Field evidence from other tectonically active sites demonstrates that boulders 728 may last in rivers for periods of tens of thousands of years (Haviv, 2007;Huber et al., 2020). However, a 729 steady-state width configuration is also dependent on how fast the channel widens in response to boulder 730 input. Direct bedrock erosion measurements from the Liwu river reveal that locally, lateral bank erosion 731 can be as significant at tens of centimeters in a single flood season (Hartshorn et al., 2002;Turowski et al., 732 2008). Hence, channel widening probably occurs much faster in the Liwu River than elsewhere. Ultimately, 733 boulder-bed channel width in the Liwu river may be at a steady-state with respect to uplift, sediment supply, 734 and discharge, but whether width has completely adjusted to boulder input requires further investigations 735 concerning the durability of boulders once they arrive into the river domain. It is also possible that we have 736 not considered different mechanisms responsible for the width anomaly in the Liwu River. 737

Steady-State in the Grade Assumption 738
Under the assumption of a grade steady state, we have assumed that sediment-flux in the boulder-bed 739 and boulder-free reaches are the same. This assumption does not require a continuous supply of boulders 740 into the channel but rather a fast response of the river to change the ability to transport sediment with the 741 same efficiency due to boulder-concentration. Specifically, according to Eq. (15), the 'grade' assumption 742 requires that sediment-bed elevation hs above the bedrock is steady in the long term. Thus, for a channel 743 that has been recently supplied with large boulders, how rapidly can a river restore its sediment transport 744 capacity? We demonstrated that grade conditions could be achieved by adjusting the sediment thickness to 745 form a new channel slope. The rate at which new sediment bed slope forms depends on various hydrological 746 and morphological parameters, such as water discharge, shear stress, and the grain size of the mobile 747 sediment (e.g., Barry et al., 2004). The Liwu river may be a locality in which large variability in water 748 discharge (Lague et al., 2005) and magnitude are expected to promote more sediment transport events in a 749 given flood season (Hartshorn et al., 2002;Dadson et al., 2003). For example, observations from the Liwu 750 River show that the river can remove sediment a few meters in depth following a significant typhoon event 751 (Lague, 2010). Even if boulders disappear quickly once arriving in the fluvial system, the timescale of 752 bedload entrainment and deposition to form a new sediment slope in general and grade conditions, in 753 particular, may correspond to a single flood (Turowski, 2020). Field evidence from various bedrock 754 channels supports recognizing an equilibrium, or 'grade,' in many bedrock river environments (Phillips and 755 Jerolmack, 2016), reinforcing the plausible assumption that the Liwu river is at an approximate sediment 756 transport steady-state, or in grade, at nearly all times. 757 758

Evaluation of the Theoretical Models 759
Under the two steady-state assumptions described in sections 4.1 and 4.2, we formalized five 760 mechanisms presumably underlying the observations of both widening and steepening of boulder-bed 761 channels (Table 2). Below, we examine and discuss the performance of the models to describe the Liwu 762 data channel width and slope predictions using the erosional balance and grade-based mechanisms 763 Five models have been considered for testing the width ratio, Wb/W: the cover, the tools, the 'tools 764 and cover,' the multi-channel effect, and the shear-stress partitioning effect. The first three models are 765 dependent on the ratio of the square root of boulder-bed to boulder-free deflection lengths (Eqs. 6, 8, and 766 9), which we assumed to be one, i.e., db = d. The deflection length likely depends on grain size, hydraulic 767 parameters (Turowski, 2020), the contact angle of the boulder with the mobile particle (Fuller et  deflection is currently unknown, but a positive correlation may account for channel widening for the above-771 discussed models beyond the predictions with db = d. 772 The cover model predicts that Wb/W equals √ ⁄ , meaning that as long as db = d, (I) the width 773 ratio does not independently depend on Γ, and (II) there is no boulder-bed widening with respect to a 774 boulder-free reach. Considering the above, although our theoretical framework of the cover effect does not 775 reproduce the Liwu width ratio data, future advances in our understanding of the relation between deflection 776 length and roughness elements could lead to a modified cover model for channel width that depends on Γ 777 as has been hypothesized for slope (Shobe et al., 2021). 778 The tools (Eq. 8) and the 'tools and cover' (Eq. 9) models predict an increase in the width ratio with 779 boulder-concentration (Fig. 6). The essential difference between the two models is in the exponent, which 780 depends on the parameter α, describing whether bedload particles are routed above boulders (α = 0) or not 781 (α = 1). At the process scale, large boulders protruding into the flow are thought to encourage sediment 782 deposition around them (e.g., Papanicolaou and Kramer, 2006;Tsakiris et al., 2014;Polvi, 2021), which 783 may lead to substantially different protrusion, causing bedload transport to alter significantly (Yager et al., 784 2007). We thus hypothesize that boulder protrusion and hydraulic behavior near boulders have an essential 785 role in controlling α.

786
The multi-channel effect (Eq. 14) predicts an increase in the width ratio with boulder concentration 787 (Fig. 7). For a given boulder-bed channel, the model plots relatively close to the data but commonly 788 overpredicts it, especially for large boulder-concentration values. Given the overall over-prediction of the 789 data and the relatively large RMSE value, we propose that the model with its current assumptions is 790 unsuitable for boulder-bed channels in the Liwu River. We envision three major potential causes for the 791 model-data deviations. 792 The model was derived using three primary assumptions: (I) the channel reach follows a specific 793 geometry, including boulder arrangement (Fig. 5), (II) sediments are redistributed evenly between the in-794 channels, and (III) the overall boulder-bed channel width independently reflects a steady-state configuration 795 for every in-channel. The Liwu boulder-bed channel reaches, however, exhibit a wide range of boulder 796 sizes and inner-reach distributions. Furthermore, at bankfull flows, when the entire bed is submerged, 797 sediments are expected to follow paths set by the flow hydrodynamics-rather than the configuration of 798 boulders-and not be evenly distributed. We believe that some of the scatter of the data concerning model 799 predictions are due to such discrepancies. To better capture the width ratio variability, specific treatments 800 of boulder distributions and sediment paths can be considered in future studies 801 The effect of shear-stress partitioning shows a humped relationship between width ratio and Γ, 802 which can be explained by the non-monotonic log-linear model Eq. (25) used to derive the model. Although 803 the effect captures 60% of the data within the errors, we emphasize two reasons for the model's inadequacy 804 of the Liwu data. First, a simpler, linear model could also account for that fraction captured by the non-805 monotonic model. Second, the parameters β and Γmax need to be different between the different reaches, but 806 independent constraints on their values are missing. As a result, the shear-stress partitioning model cannot 807 predict well the width ratio. While the physical interpretation of β is unclear, we cannot evaluate the extent 808 to which this parameter should vary among the examined reaches. Differences between channel reaches 809 could emerge from contrasts in the boulder size distributions and the streamwise and cross-sectional 810 location distributions. 811 Although the tools and shear-stress partitioning models statistically performed best overall (Table  812 2), the role of multi-channels in shaping boulder-bed channel width cannot be ruled out. Furthermore, 813 additional investigations are awaiting to unravel the role of cover in the relationship between width and 814 boulder presence (e.g., Shobe et al., 2020). However, the overall inability of the models to account for the 815 width anomaly in the Liwu River, and the long timescale of width adjustment implied from the scatter in 816 our data all drive us to suspect that different reaches have adjusted to boulder input to various degrees. 817 Two models have been considered for testing the slope ratio, Sb/S: reduction in transport efficiency and 818 shear-stress partitioning. Both were developed under the assumption of grade. We first note that the shear-819 stress partitioning effect cannot explain the slope ratio increase in the Liwu River (Fig. 10). Thus, this 820 model can be ruled out in explaining the increase in width ratio with boulder-concentration. The prediction 821 of the reduction in transport efficiency model could explain the trend observed in the data (Fig. 8) yet 822 requires calibration of four parameters. Although, according to this model, in the general case, the slope 823 ratio is a function of the width ratio, we find that the best-fit parameters are those that make the slope ratio 824 independent of the width ratio. This outcome points to a steepening effect that relies solely on sediment 825 entrainment and deposition to form a steeper bed and can occur very fast, probably within one or a few 826 floods. This mechanism differs from the one developed by Shobe et al. (2020), which relied on bedrock 827 erosion to accomplish the slope change, and would therefore have a much longer adjustment timescale. The 828 inferred independence of the slope ratio and the width ratio, manifested by the small q (power of the width 829 ratio), may be a consequence of the substantial difference in the adjustment timescales of bedrock width 830 and sediment-bed slope (Turowski, 2020). In other localities with much softer and erodible banks (e.g., 831 Cook et al., 2014), the covariation of slope and width are hypothesized to be more significant. Whereas 832 standard models commonly assume that q is either zero or one, it is also possible that the dependence of 833 sediment flux on channel width is diminished in the long term, thus constraining q to be close to zero 834 (Rickenmann, 2001). Further research on the value q for different timescales of sediment transport is 835 needed. Given the good fit and the general agreement of the model with the data, we attribute most of the 836 steepening of the boulder-bed channel reaches to a necessity to mobilize the upstream sediment supply 837 despite the presence of boulders that inhibit sediment transport efficiency. 838 The reduction in transport efficiency model further predicts a monotonic steepening effect with 839 increasing boulder concentration. However, with increasing boulder concentration, we expect channel slope 840 response to potentially reduce as the channel self-organizes a new bed largely composed of boulders such 841 that boulders are no longer significant roughness elements on the bed. This situation is equivalent to the 842 role of boulder spacing, shown by flume experiments, to strongly influence grade conditions (McKie et al., 843 2021). Each boulder generates a unique zone susceptible to flow reversals and enhanced turbulence 844 (Papanicolaou and Tsakiris, 2017). However, when the spacing is small, the different boulder-influenced 845 zones interact, causing an overall reduction in the total influence zone. Our developed equation does not 846 show this behavior because of the assumption that the transport efficiency reduces monotonically for the 847 entire range of Γ (Eq. 19). 848

Causality Relations between Boulder Concentration and Channel Width 849
The models proposed to explain the observed relation between boulder concentration and the width 850 ratio assume that channel width adjusts (and is, therefore, the dependent variable) to boulder concentration. 851 Notwithstanding, the causality between the two variables can also be presented inversely. Here we pose a 852 hypothesis for a potential dependence of boulder-concentration on channel width. Consider a case in which 853 boulders have an equal probability of arriving at a specific location within the river and assume an initial 854 natural variability in channel width along the river. In wider reaches, the fluid shear-stress deriving both 855 bedload transport, which is responsible for boulder abrasion (Wilson et al., 2013) and boulder 856 transportation, is smaller relative to a narrow channel with otherwise the same parameters. Since bedload 857 transport depends on discharge and erosion depends on bedload, boulders will both abrade and be 858 transported quicker in narrower channels. In such a case, observations would be of a positive scaling of 859 channel width with boulder concentration. However, we specify three reasons why such a case is probably 860 not valid for the Liwu River. First, we do not have direct evidence for initial width variability. Second, an 861 initial spatial variability in channel width between neighborhood reaches is less likely because our paired 862 approach of width normalization should also account for lithological and bank properties. 863 The unknown directionality between boulder-concentration and the channel width is a 'chicken or egg?' 864 problem, where cause and causality between the two variables are potentially bi-directional. Ultimately, we 865 speculate that once boulders arrive into the channel, the effects outlined above interact to form a new steady-866 state width configuration (Turowski, 2018). Nearby failure events can produce new boulders, which in turn 867 contribute to the adjustment of channel width. 868 869 7. Conclusions 870 We studied the controls of large immobile boulders on channel width and slope in bedrock channels. 871 Our data from the Liwu River, Taiwan, show that sediment-bed slope and bedrock width increase in 872 response to higher concentrations of boulders after width and slope are normalized for variations due to 873 drainage area. We invoked rock and sediment mass balance principles to explore possible mechanisms 874 responsible for the observed width and slope increase with boulder concentration. We combined them in a 875 theoretical framework for fluvial bedrock erosion and sediment transport. The theoretical framework yields 876 analytic predictions for the width and slope ratios as a function of boulder concentration. 877 Under the first principle of rock mass balance, we assumed that boulder-bed and boulder-free reaches 878 incise into bedrock at the same rate. We expanded this assumption by considering three effects that boulders 879 impose on the process of bedrock erosion: the cover, the tools, and the multi-channel effect. These models 880 yielded solutions to the width ratio. Under the second principle, we assumed that bedload flux in adjacent 881 boulder-bed and boulder-free reaches is identical under equilibrated grade conditions. Here, two underlying 882 mechanisms were examined for the effect of boulders on bedload transport: a reduction in the efficiency of 883 sediment mobilization and a reduction in the shear stress responsible for sediment mobilization. Both 884 mechanisms yielded solutions for the slope ratios, while the second mechanism provided an independent 885 solution to the width ratio. 886 The width ratio trend was best captured by the tools effect, yet this model underpredicted most of the 887 data. The shear-stress partitioning model can account for a fraction of the width data, but requires 888 knowledge on its numerous parameters, which are not constrained. Under the boulder cover effect, width 889 is insensitive to boulder concentration. Still, an emerging relationship between sediment deflection and the 890 existence of boulders is hypothesized to yield a link between width and boulder-concentration due to a 891 boulder cover effect. The multi-channel effect demonstrated an increase of width but predicted larger width 892 values than mostly seen in the Liwu data. We conclude that the scatter in the data points out to a long 893 timescale adjustment of channel width, and different reaches have adjusted to different extents. 894 The slope ratio was best captured by the effect of a reduction in sediment transport efficiency, with 895 little to no dependence on channel width. In contrast, a shear-stress partitioning effect cannot explain the 896 slope ratio. 897 The theoretical framework presented in this paper is a first attempt to examine and test various physical 898 mechanisms controlling the relationship between bedrock channel morphology and large boulders. We 899 acknowledge two primary key future research questions emerging from our work. First, we have 900 insufficient insights into the dynamics of bedload relating to sediment deflection in vicinity to boulders.

Introduction
The following supporting information describes additional methods, calculations, and extended data analysis that aid this study.

Text S1: Error in Width Calculations
The field data width factor Z = Wb/W includes various uncertainties, including (I) the error on the 3D model-derived orthophoto used to calculate the channel width. Based on one field location in the Liwu river where an orthophoto was compared to numerous measurements using a tape, this error is less than 20 cm and is thus very small in comparison to the actual channel width. (II) Uncertainty related to the channel reach being non-homogeneous in space. In other words, the measurement of a cross-section channel width is dependent on the location along the downstream direction. Here we assume that this is the most significant uncertainty. To calculate it, we utilize an error propagation technique according to the following. Assuming that uncertainties are normally distributed, the uncertainty ΔZ is given by The boulder-bed and boulder-free widths, Wb and W, respectively, were derived using the methods described in the main text (Section 2). The error on the boulder-bed width ∆ was calculated as one standard deviation from the mean of 10 width measurements using a 3D model-derived orthophoto for each boulder-bed channel reach. In contrast, the boulder-free width ΔW was estimated using the maximal error on width derived from Google-earth imagery (Section 2 in the main text) and set a constant of ΔW = 0.2W. With these, each boulder-bed channel width can be assigned with a specific error length.

Text S2: Comparison of Width Calculations Using Basin-scale Relationship and Direct Google-earth-derived Data
Boulder-concentration is generally smaller (~0.1) in smaller upstream reaches and shows large variability of 0 -0.34 in larger drainage areas. Thus, the entire range of Γ is observed in drainage areas > 400 km 2 . The two different approaches for calculating W yield relatively similar width ratios; 18 out of 20 are within 50% error relative to a one-to-one case (Fig. S1). Although the choice of 50% is somewhat arbitrary, the two marked outliers show a large discrepancy between the two applied approaches.

Text S3: Recurrence of Hillslope Failure Events in the Liwu River
The Liwu River forms a steep landscape with elevations drop of ~ 3000 m over a relatively short distance of 40 km (Fig. S1)

Text S4: Testing the Horizontal Errors in Drone-Derived Three-Dimensional Models
To test the errors in boulder-concentration and boulder size analysis, we have performed the following study. A field locality was chosen in Baiyang (24.185005, 121.485815), a touristic trail with a concrete bridge crossing the tributary. Using a standard meter, we have measured by hand the length of seven objects with different lengths. Those objects included the bridge length, its width, and other observable shapes. Additionally, a drone was used to photograph the bridge vicinity using the same method described in Section 2.1. We followed the same procedure (section 2.1) to generate three-dimensional models, from which we generated an orthophoto. Then we measured the length of the same seven objects using the orthophoto and compared them to the lengths estimated by hand at the site. The comparison (Fig. S4) shows that the lengths fall closely to a one-to-one reference line, indicating that our model includes minimal horizontal errors. The RMSE between model and observations is 6.5 cm. Figure S1: A 30 m DEM map of the Liwu River (630 km 2 ) overlain with the surveyed reach locations of boulder-bed channels in light blue circles. The location details are listed in Table 2 in the main text.