Controls on the hydraulic geometry of alluvial channels: bank stability to gravitational failure, the critical-flow hypothesis, and conservation of mass and energy

The bankfull depths, widths, depth-averaged water velocities, and along-channel slopes of alluvial channels are approximately power-law functions of bankfull discharge across many orders of magnitude. What mechanisms give rise to these patterns is one of the central questions of fluvial geomorphology. Here it is proposed that the bankfull depths of alluvial channels are partially controlled by the maximum heights of gravitationally stable channel banks, which depend on bank 10 material cohesion and hence on clay content. The bankfull depths predicted by a bank-stability model correlate with observed bankfull depths estimated from the bends in the stage-discharge rating curves of 387 U.S. Geological Survey gaging stations in the Mississippi River Basin. It is further proposed that depth-averaged water velocities scale with bankfull depths as a result of a self-regulatory feedback among water flow, relative roughness, and channel-bed morphology that limits depth-averaged water velocities to be within a relatively narrow range associated with Froude numbers that have a weak inverse relationship 15 to bankfull discharge. Given these constraints on channel depths and water velocities, bankfull widths and along-channel slopes consistent with observations follow by conservation of mass and energy of water flow.


Introduction
The bankfull depths, h, widths, w, depth-averaged water velocities, v, and along-channel slopes, S, of alluvial channels exhibit power-law relationships with bankfull discharge, Q: 20 where k ≈ 0.4, b ≈ 0.5, m ≈ 0.1, and z ≈ -0.4 (Leopold and Maddock, 1953). Many studies have proposed a process-based understanding of these patterns (see reviews by Ferguson (1986) and Singh (2003)), but none has achieved widespread acceptance. Schumm (1960) documented that alluvial channels with sand-rich bed and bank material tend to be wide and shallow, while 25 alluvial channels with silt-and-clay-rich bed and bank material tend to be narrow and deep. Schumm's findings have led nearly all subsequent researchers to consider the resistance of bank material to fluvial erosion to be the key factor controlling alluvial channel width (e.g., Parker, 1979;Eaton and Millar, 2004;Dunne and Jerolmack, 2020). Bank retreat, however, is also driven by gravitational failure (ASCE, 1999), a process that limits bank heights to values that depend on bank material cohesion and hence on clay content. The gravitational failure of channel banks may partially control bankfull depths via a self-regulatory 30 mechanism in which channel incision and/or floodplain deposition tend to increase bank height, triggering bank failure when a critical bank height, dependent on bank material cohesion, is exceeded (Andrews, 1982), thus introducing new sediment into the channel bed that, as it is redistributed by fluvial processes, tends to reduce the channel depth back towards a critical value. This paper demonstrates that bankfull depths predicted by a bank-stability model correlate with observed bankfull depths estimated from the bends in the stage-discharge rating curves of 387 U.S. Geological Survey (U.S.G.S.) gaging stations in the 35 Mississippi River Basin (MRB). This analysis supports the hypothesis that the gravitational failure of channel banks partially controls bankfull depths and complements the recent work of Dong et al. (2019) that documented a relationship between bankmaterial texture and the hydraulic geometry of alluvial channels in the Selenga River Delta. Grant (1997) proposed a critical-flow hypothesis in which depth-averaged water velocities are self-regulated via interactions between the water flow and the channel-bed morphology. Grant (1997) argued that, in steep (≳ 0.01 m m -1 ) channels, the 40 Froude number rarely exceeds one for extended periods of time due to interactions between the free surface and the bed that result in an approximate balance between forces that accelerate the flow and forces that extract energy from the flow. Such a balance may also extend to coarse-bedded channels, the water flow in which is prone to wave drag associated with flow around bed sediment grains that protrude above the water surface (Wohl, 2013). Wave drag can be expected to be more common in coarse-bedded channels relative to channels with finer bed sediments both because they have relatively large bed roughness 45 elements that more readily protrude above the water surface and a tendency towards shallower flows as a result of their characteristically large width-to-depth ratios (Schumm, 1960).
Central to the critical-flow hypothesis is the existence of a self-regulatory feedback in which an increase in velocity is met with an increase in drag that tends to reduce the velocity and hence the Froude number. Similarly, a decrease in velocity tends to be met with a decrease in drag that tends to increase the velocity. Here it is hypothesized that such a self-regulatory feedback 50 is not limited to steep channels. In less-steep, sand-bedded channels, faster flow tends to facilitate the development of larger and more well-developed bedforms (which tend to form at Froude numbers ~0.1-1 (e.g., Simons and Richardson, 1966)) that increase relative roughness and hence drag. The mechanisms of self-regulation and the Froude numbers at which steep and less-steep channels may achieve this self-regulation thus differ, but both are likely to have self-regulatory interactions between the flow and the bed that limit the Froude numbers of bankfull discharges. Consistent with this hypothesis, here it is 55 documented that bankfull Froude numbers, and hence the ratio of depth-averaged water velocities to the square root of bankfull depths, tend to be within a relatively narrow range that has a weak inverse relationship to bankfull discharge.
The bankfull widths of alluvial channels are set by the requirement that channels convey geomorphically effective water discharges. Conservation of mass, together with the clay-content control of bankfull depths and the Froude-number control of water velocities, thus constrains bankfull widths. 60 Conservation of energy constrains along-channel slopes. The conversion of gravitational potential energy into the kinetic energy of water leads to a relationship among along-channel slopes, bankfull Froude numbers, bankfull depths, and the size of the largest bed roughness elements, which in gravel-bedded channels tend to be bed sediment grains and in sand-bedded channels tend to be ripples and/or dunes.

Controls on bankfull depths
The maximum stable height, hc, of an alluvial channel bank subject to gravitational shear failure is proportional to bankmaterial cohesion, C (Taylor, 1937;Terzaghi and Peck, 1967;Hunter and Schuster, 1968;Chen et al., 1969;ASCE, 1999): where ρ is the bulk density of the bank material, g is the acceleration due to gravity, and Ns is a stability parameter dependent 70 on the geometry of the potential failure surface (e.g., planar, log-spiral, or circular), the pore pressure of the bank material (which is governed by the water table position if the pore pressure is assumed to be hydrostatic), and the angles of the bank and of internal friction (see Table 1 for a list of variables).
In order to estimate a reference Ns value appropriate for understanding how gravitational stability may influence the scaling of alluvial channel bankfull depths to discharge, a steep bank (i.e., near-vertical at the top of the bank but decreasing to 75 approximately 45˚ near the toe) with an internal friction angle of 35˚ (typical for a loamy or clayey sand), near-saturated conditions, and a log-spiral potential failure surface were assumed. Near-saturated conditions are consistent with the fact that gravitational shear failure has been documented to occur most frequently during the falling limbs of flood discharges when pore pressures tend to be associated with near-saturated conditions (e.g., Casagli et al., 1999;Simon et al., 2000). Chen et al. (1969) derived Ns values for prescribed angles of the bank and of internal friction for unsaturated conditions. For a friction 80 angle of 35˚, Ns values in Table 1 of Chen et al. (1969) decrease with increasing bank angle from Ns = 22 for a 60˚ bank to Ns = 12 for a 75˚ bank and Ns = 7.5 for a vertical bank. Hunter and Schuster (1968) limited their analysis to cases with no internal friction (hence their absolute Ns values are not applicable here) but documented an approximately 3-fold decrease in Ns values from unsaturated conditions (i.e., M = hwγw/hcγ´ ≈ 1, where is hw is the depth to the water table below the top of the bank, γw is the unit weight of water, and γ´ is the submerged unit weight of the bank material) to near-saturated conditions (i.e., M = 0). 85 Combining the results of Chen et al. (1969) and Hunter and Schuster (1968) suggests that Ns values for a saturated bank with an internal angle of friction of 35˚ vary from approximately 7.3 for a 60˚ bank to Ns ≈ 4 for a 75˚ bank and 2.5 for a vertical bank.
Bank material cohesion varies linearly from 0 (for cohesionless sand) to approximately 90 kPa (for pure clay) for moisture contents in the range of 7 to 40% according to a least-squares linear regression of the data from Dafalla (2013) (Fig. S1): 90 ≈ (900 ± 70) c . (3) where the units of C is Pa and pc is percent. The uncertainty in Eq. (3) is the standard error (i.e., one standard deviation) resulting from the regression.
Combining Eqns. (2) and (3) and assuming a bulk density of 1700 kg m -3 and a value of Ns ≈ 6 (corresponding to a nearsaturated bank with an angle of approximately 65˚, i.e., an average angle for a bank that is near-vertical at the top and decreases 95 to an angle of approximately 45˚ near the toe) yields ℎ c ≈ 0.35 c .
(4) Absent site-specific data for bank angles, the largest source of uncertainty in the proportionality coefficient in Eq. (4) as applied to specific locations is likely the bank angle, since relatively modest variations in bank angle (e.g., from 90˚ to 75˚) are associated differences in Ns of approximately a factor of 3, i.e., much larger than other sources of uncertainty such as that 100 between cohesion and clay content as quantified by Eq. (3). Section 4 provides discussion on how uncertainty in bank angle and other factors such as bank vegetation limit the precision of Eq. (4) to specific locations. The primary of objective of this paper, however, is to document an increase, on average, in bankfull channel depth with increasing bank-material clay content: assuming that bankfull depth is approximately equal to the maximum gravitationally stable bank height. 105 To test the tendency for channel depth to increase, on average, with increasing bank-material clay content as predicted by Eq. (5), the bankfull depths for 387 U.S.G.S. gaging stations in the MRB were estimated using the bends in the stage-discharge rating curves (Fig. S2) for each station. Predictions of bankfull depth using Eq. (5) were then compared to the observed bankfull depths derived from the rating curves. The gNATSGO soil database (Soil Survey Staff, 2019) was used to estimate the percent clay content of the floodplain deposits adjacent to each station. This analysis focuses on the MRB because there is no readily 110 available soils database for any other continental-scale river basin that resolves floodplains in comparable detail and is based on a similar richness of field-based soil texture measurements.
The bankfull depth for each U.S.G.S. gaging station was estimated using the intersection of the linear regressions of peak annual gage height (used as a proxy for stage) to peak annual discharge obtained using the five smallest and five largest discharges in each record (shown as gray circles in the example of Fig. S2). This intersection, or bend, in the stage-discharge 115 rating curve can identify the stage and discharge above which overbank flow occurs (Copeland et al., 2000), provided that the slope of the high-flow linear regression is much smaller than the slope of the low-flow linear regression, which is a signature of the abrupt widening of flows as they transition from in-channel to overbank.
The analysis presented here began by including data from all U.S.G.S. gaging stations in the MRB with available peak discharge data (U.S. Geological Survey, 2020). Only those stations for which the slope of the high-flow linear regression is at 120 least five times smaller than the slope of the low-flow linear regression were retained. In addition, only those stations that had at least 20 years of data, have a contributing area larger than 100 km 2 , are not located close to major infrastructure (based on an inspection of each station location in Google Earth imagery), and have a resulting bankfull depth of greater than 2 m were retained. Channels with bankfull depths less 2 m were removed because such channels tend to be associated with low clay contents that are inherently difficult to estimate in the field (see Sect. 3.1 for more detail on the potential bias associated with 125 estimating low clay contents). In order to further filter out stations where the low-flow linear regression is potentially unrepresentative of the hydraulic behavior of in-channel discharges, only those stations for which the extrapolation of the lowflow linear regression is close to zero flow depth at zero discharge were retained. Stations for which the low-flow regression does not extrapolate to a flow stage close to zero at zero discharge may have gage height data that are not an accurate proxy for flow stage and/or have other data quality issues that preclude an accurate estimate of bankfull depth using the stage-130 discharge rating curve. What represents "close" should not be based on an absolute error, e.g., 0.5 m, because such a criterion would require that the low-flow regressions for deep channels be relatively more accurate than those for shallow channels.
Here "close" was defined as being within 50% of the bankfull stage from zero. That is, if the bankfull stage is 5 m, then the extrapolation of the low-flow regression to zero discharge must yield a stage within 2.5 m of zero in order for that station to be retained in the analysis. Similarly, if the bankfull stage is 2 m, then the extrapolation of the low-flow fit to zero discharge 135 must be within 1 m of zero. A total of 387 stations met these criteria.
To estimate the floodplain clay content for each of the 387 U.S.G.S. gaging stations, the depth-averaged percent clay content from soil depths of 5 to 150 cm was computed for every pixel within the MRB using gNATSGO. These depths were chosen to avoid the soil horizon close the surface (typically the O and/or A horizon, which may have clay contents reflective of surficial biological processes that are not representative of the rest of the profile) and because soil properties at depths greater 140 than 150 cm can be inconsistently encoded in U.S. soil databases (Miller and White, 1998) (Archuleta et al., 2017). Bankfull depths predicted by Eq. (5) were then compared to observed bankfull depths using a Pearson correlation coefficient, a root-mean-squared error (RMSE), the percentage of values correctly predicted to within a factor of 2, and a p-value that quantifies the likelihood of the null hypothesis that the predicted and observed bankfull depths may be correlated by chance.
To assess the potential impact of errors in measurements of percent clay content on predictions of the maximum stable bank 150 height and therefore of bankfull depth using Eq. (5), synthetic predictions for bankfull depths, hpred,syn, were generated equal to 0.35 times samples of synthetic percent clay content, pc,syn, drawn from a lognormal distribution designed to mimic the distribution of bankfull depths of U.S.G.S. gaging stations in the MRB, plus a normally distributed random error with a mean of zero and standard deviation of σ: where η is a normally distributed random variable with a mean of zero and a standard deviation of 1. For σ = 0, Eq. (6) produces synthetic data precisely equal to Eq. (5). With finite values of σ, Eq. (6) produces synthetic data with scatter that can be used to assess how errors in percent clay content may affect the relationships between observed and predicted bankfull depths.

Controls on depth-averaged water velocities and bankfull widths 160
The critical-flow hypothesis implies that bankfull Froude numbers, F, are limited to a relatively narrow range of values with an upper limit close to 1 for gravel-bedded channels and a similarly narrow but somewhat lower range of values for sandbedded channels. In Sect. 3.2 it is demonstrated that this variation in F can be quantified using a power-law relationship between F and Q: ∝ .
(7) 165 The power-law exponents reported in this paper were determined via a least-squares linear regression of the logarithms of the data. The definition of Froude number provides a linkage among depth-averaged bankfull water velocities, bankfull Froude numbers, and bankfull depths: Equations (7) and (8) thus constrain the value of m to be 170 The exponent b in the relationship between bankfull width and discharge can be constrained by conservation of mass of water assuming steady, uniform flow, consistent with many previous models for the downstream hydraulic geometry of alluvial channels (e.g., Lindley, 1919;Smith, 1974;Ferguson, 1986;Huang et al., 2004;Julien, 2014): (10) 175

Controls on along-channel slopes
The Darcy-Weisbach equation is based on conservation of energy for steady, uniform flow, i.e., that the gravitational potential energy per unit mass, ghS, 180 produces a depth-averaged water velocity associated with a drag force per unit mass exerted by the channel bed on the water equal to (f/8)v 2 , where f is the friction factor (e.g., Ferguson, 1986). Equation (11) can be rewritten as Here the Variable Power Equation (VPE) of Ferguson (2007) is used to quantify f: where a1 and a2 are coefficients (equal to 6.1 and 2.4 based on the least-squares minimum error for velocity in the calibration performed by Ferguson (2007)), and β is the relative roughness. The VPE equation was chosen because it provides a transition between the Manning-Strickler 1/3 rd power scaling of friction factor to relative roughness for β ≫ 1 (nearly all sand-bedded channels and some gravel-bedded channels) to a quadratic scaling when β is ≲ 1 (channels with very coarse beds) that accords well with available data (Ferguson, 2007). 190 Relative roughness depends on whether or not bedforms are present. Table S2 of Ohata et al. (2017) identifies the range of F and d50 values conducive to ripple and/or dune development in alluvial channels. Of the 3790 data points in the Ohata et al.
(2017) dataset, 1574 (42%) have ripples or dunes, of which 19% are in the field and the remaining 81% are in the laboratory.
By cross-referencing those results with the Dunne and Jerolmack (D&J) (2018) global dataset (i.e. by drawing a curve in F vs. d50 space that separates channels that have ripples and/or dunes from those that do not, and assuming the existence of ripples 195 and dunes in channels of the D&J dataset that have F and d50 values that sit above and to the left (Fig. 3B)  For gravel-(and-coarser)-bedded channels, relative roughness is defined in the calibration of Ferguson (2007) as the ratio of 84 th percentile of bed grain diameter to the hydraulic radius. Since the analysis of this paper relates Eq. (13) to data from the D&J global dataset that uses d50 instead than d84, it is assumed that d50 ≈ d84/2 and, because w > 10h for all points in the D&J global dataset, that the bankfull hydraulic radius is approximately equal to the bankfull depth. The relative roughness for gravel-bedded channels, βg, in the D&J global dataset can, therefore, be approximated as 205 For sand-bedded channels, relative roughness can by estimated as (Bathurst, 1993): where H is the bedform height and α is the ratio of bedform height to length, L.
Combining Eqs. (12)-(15) gives an equation for the along-channel slopes of gravel-bedded channels consistent with 210 conservation of energy: and an analogous equation for sand-bedded channels: The heights and lengths of ripples and dunes can be estimated as (Yalin, 1964): 215 and ≈ 1000 50 .
Equation (17), therefore, can be rewritten in terms of d50/h for a prescribed value of α as  Figure 1 illustrates the tendency for the clay contents of floodplain deposits adjacent to many smaller channels in the MRB to be lower than those of larger channels. For example, the North and South Platte Rivers have typical floodplain clay contents <10%, while the Platte River has typical floodplain clay contents of ≈10-20%, and the Missouri and Lower Mississippi Rivers 225 have typical clay contents of ≈20-30%. There are many relatively small channels, however, that have clay contents > 30% due to clay-rich local bedrock. As such, there isn't a precise, one-to-one correlation between clay content and contributing area or discharge, but rather a general tendency for channels conveying larger discharges to have more clay-rich floodplain deposits.  depths, a RMSE of 1.7 m, and 84% of the data points are within a factor of 2 of the observed bankfull depth. The p value, i.e., the chance that the correlation between hpred and hobs could have occurred by mere chance, is ~10 -17 .

Controls on bankfull depths
Figures 2(b) and 2(c) illustrate the potential impact of errors in percent clay contents on predicted bankfull depths using Eq.
(5) with σ = 0.06 and 0.12 (6 and 12%), respectively. These σ values were chosen because, while gNATSGO does not provide an error estimate, a recent soil property dataset created using machine learning algorithms has an estimated RMSE of 12% 240 (Ramcharan et al., 2018) and a value half that size allows for the effect of error size to be assessed. For a relatively small error (σ = 6%), there is a spread of values around the 1:1 line, with a larger relative spread for smaller clay contents, i.e., a 6% error results in a 100% relative error for a percent clay content of 6% (i.e., hpred,syn ≈ 2 m) but a 50% relative error for a percent clay content of 12% (i.e., hpred,syn ≈ 4 m). As the σ value increases to 12%, the spread of values around the 1:1 line increases as expected but hpred,syn values less than approximately 2 meters also appear to be biased upward (i.e., the geometric mean of 245 hpred,syn values deviate from the 1:1 line).
This upward bias may be associated with the difficulty of measuring very low clay contents in the field. Clay contents estimated in the field can only be constrained to be within the range from 0 to 10% clay. If the actual clay content is close to 0 (e.g., < 1%), the clay content estimate is likely to be overestimated by a much larger fraction than would be the case for a larger clay content (e.g., 5% clay is 400% greater than 1% clay but only 50% lower than 10% clay). The result may be an 250 upward bias in clay content for clay contents less than approximately 10%, which, according to Eq. (5), may be associated with channels less than 2-3 m in bankfull depth. Channels with bankfull depths less than 2 m (Sect. 2.1) were excluded from the analysis due to this potential bias.

Controls on depth-averaged water velocities and bankfull widths
Ripples and dunes tend to form at lower Froude numbers in channels with finer bed sediments (Fig. S3). The curve in Fig. S3 used to identify the range of F and d50 conditions conducive to ripple and/or dune development is reproduced in Fig. 3, where 96% of sand-bedded channels in the D&J global dataset are below the upper limits of F and d50 conducive to ripple and/or 265 dune development identified using the Ohata et al. (2017) dataset. As such, it is assumed for the purposes of this analysis that sand-bedded channels in the D&J global dataset are dominated by bedform roughness while gravel-bedded channels in the D&J global dataset are dominated by grain roughness.  Fig. S3.

Small circles correspond to sand-bedded channels (d50 < 2 mm), large circles correspond to gravel-bedded channels (d50 > 2 mm). The curve in (b) defining channels that likely have ripples and/or dunes is based on the subset of 3791 field studies and experiments compiled by Ohata et al. (2017) and graphed in
Using the D&J global dataset, alluvial channel depths and widths scale with bankfull discharges to the 0.402 ± 0.006 and 0.512 ± 0.007 powers, respectively (R 2 = 0.86 and 0.89) (Fig. 4). If we take the -0.116 ± 0.009 Froude-number-discharge scaling exponent obtained from least-squares regression to the logarithms of the data (Fig. 3(a)) and the 0.402 ± 0.006 depthdischarge scaling exponent as a starting point, Eqs. (9) and (10) predict a width-discharge exponent of 0.51 ± 0.01, i.e., precisely equal to that observed in the D&J global dataset. 280

Controls on along-channel slopes
Along-channel slopes predicted by Eqs. (16) and (20) are consistent with observed values in the D&J global dataset (Pearson correlation coefficient of 0.77) (Fig. 5). For sand-bedded channels, in which ripples and/or dunes are likely to be the dominant roughness elements, the predicted values plotted in Fig. 5 assume τ c /τ 0 ≈ 0 (consistent with the suspended-load-dominated conditions common in sand-bedded channels (Dade and Friend, 2000)) and a representative α value of 0.05 (based on the range 290 0.05-0.1 reported by Guy et al. (1966)). Figure S4 illustrates the sensitivity of the predictions of the along-channel slopes of sand-bedded channels to the presence/absence of bedforms and the assumed value of α. For alternative scenarios in which a) an unrealistically large value α = 0.25 is assumed, and b) bed grains are assumed to be the dominant roughness elements (i.e., no bedforms are present), Eqs. (16) and (20) predict along-channel slopes that are approximately an order of magnitude above and below observed values, respectively. 295

Discussion 300
The model of this paper posits that floodplain deposit clay contents partially control the maximum stable heights of channel banks. This control is unlikely to result in a precise correlation between clay contents and bankfull channel depths for at least two reasons besides data errors. First, the incised depth of a channel that flows through a section of higher clay content to a section of lower clay content may be more strongly controlled by the lower clay content because the downstream reach may act as the local base level of erosion for the upstream reach. Second, alluvial channels can adjust to spatial variations in bank 305 material cohesion by varying the bank angle in addition to the bank height (Knighton, 1974). For a 35˚ angle of internal friction, for example, a bank angle of 60˚ has a stability factor, Ns, that is approximately three times larger than a vertical bank with the same cohesion (e.g., Fig. 3 of Chen et al. (1969)). As such, an alluvial channel that flows through a section with less cohesive bank material compared to neighboring sections may adopt a less steep bank (thus increasing the stability factor, Ns) instead of, or in addition to, becoming wider and shallower in order to minimize variations in channel depth that might otherwise drive 310 large spatial variations in rates of aggradation/incision. Further tests of the model of this paper may require a better understanding of how channel depths and/or bank angles adjust to spatial variations in bank material cohesion through a channel network.
It is also important to consider how potential errors in the data may contribute to the observed scatter in Fig. 2. The estimates of bankfull depths presented here are not exact because gage height (the height of the water above a reference point) is used 315 as a proxy for flow stage. Also, uncertainties in clay content of just 5-10% percent are capable of creating scatter comparable to that in Fig. 2. Despite such large potential errors and the bias they may introduce into predictions of bankfull depths, the analysis presented here rules out the possibility that floodplain deposit clay contents and bankfull depths are related by chance (i.e., p ~ 10 -17 ) and it demonstrates that Eq. (5) predicts bankfull depths to within a factor of 2 of the observed bankfull depths for 84% of the 387 stations included in the analysis. 320 The model of this paper posits that bankfull depth may be self-regulated via a tendency for an increase in bank height caused by channel incision and/or floodplain deposition to trigger bank failure when a critical bank height, dependent on bank material cohesion, is exceeded. An important role for gravitational failure in controlling alluvial channel geometry is consistent with process-based studies of bank retreat, in which gravitational shear failure has been documented to occur frequently during the falling limbs of flood discharges when pore pressures tend to be highest (e.g., Casagli et al., 1999;Simon et al., 2000). An 325 important role for shear failure in bank retreat is also consistent with Li et al.'s (2015) finding that bankfull depth correlates with fluid viscosity because viscosity controls the permeability of bank material and hence pore pressures and therefore susceptibility to gravitational shear failure. Fluvial erosion of the bank toe is still necessary to remove material slumped from the bank into the channel and likely plays an important role in driving bank retreat in channels with cantilevered banks (Pizzuto, 1994). However, if the bank is more than twice the height of the zone of scour, gravitational shear failure nevertheless will be 330 the process by which the majority of material is removed from the bank (Tao et al., 2019). That said, gravitational failure and fluvial shear stresses act in concert in such a way that identifying which process is dominant may be difficult. Fluvial scour transports slumped material away from the bank, likely keeping bank angles higher than they would be without fluvial scour.
Gravitational failure transports material from higher on the bank to the toe (often following fluvial scour, which is maximized at the toe), likely keeping bank angles lower than they would be without gravitational failure. More research is needed, but the 335 relative importance of these two processes may be reflected in the bank angle, i.e., a bank angle persistently much less than vertical may indicate the dominance of gravitational failure while a bank angle persistently at or above 90˚ (i.e., a cantilever or overhang) may indicate the dominance of fluvial scour.
The model of this paper is simplified in at least two specific ways that bear mentioning: it does not account for tension cracks that, if present, can lower the maximum stable height below that predicted by Eq. (2) (Darby and Thorne, 1994), nor 340 does it account for the role of vegetation in bank stability, which can increase bank heights by at least a factor of two over values predicted using bank material texture alone (Huang and Nanson, 1998).
Bank vegetation plays a significant role in controlling bank stability, but it is unlikely that such control is responsible for the scaling relationships that are the focus of this paper, as such scaling relationships exist across climatic regions with very different vegetation characteristics. In addition, the stability of any bank is primarily a function of its weakest portion or layer, 345 i.e., failure is more likely to occur in a zone of low shear strength compared to a zone of high shear strength, all else being equal. Failure of a weaker zone underlying a stronger zone can create a cantilever, but such cantilevers cannot continue to grow indefinitely, hence rates of long-term retreat in stronger and weaker zones of the same bank will tend to be similar. This, together with the fact that vegetation is likely to strengthen just the uppermost approximately 2 m of channel banks (globally, >99% of roots are found in the uppermost 2 m of the soil (Jackson et al., 1996)), suggests that bank material shear strength is 350 unlikely to be controlled primarily by vegetation in channels deeper than 2 m. A dependence of bank retreat rates on vegetation has been documented in many studies (e.g., Ielpi and Lapôtre, 2020). However, such studies may not fully account for the fact that wetter climates with more vegetation also tend to have more clay-rich soils, leaving open the question of whether it is bank vegetation or material texture that is most responsible for bank resistance to erosion.
Four additional assumptions of the model should be emphasized. First, the model does not account for tension cracks that, 355 if present, can lower the maximum stable height below that predicted by Eq. (2) (Darby and Thorne, 1994). Second, the model assumes a maximally incised channel, i.e., a channel that has incised to the point of reaching the threshold of bank stability quantified by Eq. (2). Quaternary climatic changes have driven cycles in which channel aggradation has been followed by a positive feedback of incision and channel narrowing (e.g., Bull, 1991) that have likely made many alluvial channels prone to an incised state. Such considerations, however, do not apply to some types of alluvial channels, e.g., small-scale channels 360 formed in the laboratory and those in which sediment supply is not heavily influenced by climatic changes. Third, the model of this paper involves no explicit constraint on the width-to-depth ratio of alluvial channels despite the fact that such a constraint may play an important role in some cases. An important concept in rill erosion is that unit stream power is maximized for a width-to-depth ratio of ≈ 0.5-3, with the specific value dependent on the cross-sectional functional form (e.g., Moore and Burch, 1986). A similar concept may limit incision in channels with small width-to-depth ratios (Huang and Nanson, 2000), 365 reducing the likelihood of channels with w/h ≲ 1 and hence potentially limiting h to values smaller than that set by Eq. (2) for channels with small discharges. Such a control does not seem likely in the channels of the D&J global dataset (since w/h ≳10), but it may play an important role in some types of channels and should be explicitly considered in future research. Fourth, the analysis assumes that the dominant sources of roughness in the channels of the D&J data are ripples/dunes or gravel clasts.
Long-wavelength topographic features such as bars and meanders are not likely to be dominant roughness/drag-inducing 370 elements given that the presence/absence of the flow separation that tends to dominate drag depends sensitively on the maximum slope of bedforms and other obstacles to the flow, with slopes in excess of 0.2 m/m generally needed to trigger flow separation (e.g., Lefebvre et al., 2014). Vegetation can certainly be a dominant source of roughness on the beds of ephemeral channels, however, and some of the scatter in the analysis of this paper may be a result of vegetation-induced bed roughness.
More data on the relationship between bank cohesion and bank material texture is needed. This study used clay content as 375 a predictive variable for cohesion in part because a transfer function is available between clay content and cohesion based on the work of Dafalla (2013) used in Figure S1. However, silt content also affects cohesion (Huang et al., 2006) and the findings of Schumm (1960) suggest that both the silt and clay content of bank material are relevant to understanding alluvial channel geometry. The specific surface area of bank material may be a more accurate predictive variable for cohesion than either clay content or silt and clay content (weighted equally), as specific surface area includes both the silt and clay contents but weighs 380 the presence of clays more heavily (Huang et al., 2006).
The analysis presented here used a power-law relationship to quantify the relationship between F and Q (Eq. (7)). However, a scaling break appears to exist in the D&J global dataset (Fig. 3(a)), with F values approximately constant for Q values less than ~10 2 m 3 s -1 (i.e., discharges dominated by gravel-bedded channels). This break may be consistent with the critical-flow hypothesis, i.e., steep, predominantly gravel-bedded channels in the D&J global dataset may have F values predominantly in 385 the range of 0.3-1, independent of channel size/discharge, because the general lack of ripples and dunes in gravel-bedded channels requires that the increase in drag near critical flow conditions be caused primarily by wave drag, while in sand-bedded channels the increase in drag near critical-flow conditions is likely to be caused by ripples and/or dunes that form at a range of Froude numbers weakly dependent on bankfull discharge. It would be useful for future research to investigate the relationships among F, Q, bed texture, and the presence/absence of ripples and/or dunes to more fully test this hypothesis and its implications 390 for potential breaks in scaling within Eq. (1).
The model of this paper implicitly includes the sediment-supply control on along-channel slope documented by Li et al. (2005), Pfeiffer et al. (2017) and Blom et al (2017). Drainage basins with higher rates of sediment supply erode faster, resulting in coarser sediments being delivered to channels, all else being equal (Attal et al., 2015). Coarser sediments increase alongchannel slopes because steeper slopes are necessary to achieve critical or near-critical water velocities in alluvial channels with 395 coarser bed sediments (Eqs. (16)&(20)).

Conclusions
This paper proposed that the bankfull depths of alluvial channels may be partially controlled by the maximum heights of gravitationally stable channel banks, which depend on bank material cohesion and hence on clay content. Bankfull depths predicted by a bank-stability model correlate with the observed bankfull depths estimated using the bends in the stage-discharge 400 rating curves for 387 U.S.G.S. gaging stations in the Mississippi River Basin. It was proposed, inspired by the critical-flow hypothesis of Grant (1997), that depth-averaged water velocities scale with bankfull depths as a result of a self-regulatory feedback among water flow, relative roughness, and channel-bed morphology that limits velocities to be within a relatively narrow range associated with Froude numbers that have a weak inverse relationship to bankfull discharge. Given these constraints on channel depths and depth-averaged water velocities, bankfull widths and along-channel slopes consistent with 405 observations follow from conservation of mass and energy of water flow. The model of this paper provides a novel processbased understanding of the hydraulic geometry of alluvial channel networks. Section 4 identified research needs that would enable a more comprehensive test of the hypotheses and a better understanding of the bank-textural controls on the hydraulic geometry of alluvial channels more generally.

Data Availability
The Supplementary Material contains all of the data used in the paper not published elsewhere. 410 Competing interests The author declares that he has no competing interest.