The importance of threshold in alluvial river channel geometry and dynamics

Civil and Environmental Engineering, Utah State University Earth and Planetary Sciences, Washington University in Saint Louis School of Geography and the Environment, University of Oxford, Oxford, UK Earth, Environmental, and Planetary Sciences, Rice University Civil and Environmental Engineering, University of Florence Civil, Environmental and Architectural Engineering, University of Padua Earth and Environment, Franklin and Marshall College Universite de Paris, Institut de Physique du Globe de Paris, CNRS, F-75005, Paris,


Introduction
In the limit of no sediment supply, experiments with laminar and turbulent flow demonstrate that channels evolve to a 151 threshold condition with a cross section in the shape of a cosine 46, 48, 114 -where fluid and gravitational stresses everywhere 152 on the bed are balanced by friction. This reduces the solution for the stable threshold channel to a hydraulic problem: with 153 expressions for fluid-mass conservation and flow resistance, the shape and slope of the channel can be predicted if one imposes 154 values for: discharge, sediment properties (D 50 , R) and flow resistance (C f ) (Box 2). In the absence of bed forms, flow 155 resistance arises primarily from grain-scale roughness and hence C f may be estimated from D 50 53, 114 . Alluvial rivers are not 156 canals of course; they regularly transport their bed sediment, and therefore experience fluid stresses in excess of threshold. Yet 157 many alluvial rivers maintain stable banks (on average), which would appear to require fluid stresses at or below threshold.

158
Gary Parker referred to this problem as the 'stable channel paradox' 50 , and presciently stated " [such] paradoxes are often 159 resolved in terms of singular perturbation analysis". This suggests that sediment transport may be treated, conceptually and 160 mathematically, as a perturbation to the threshold state; and that the corresponding average stress condition is τ = (1 + ε) τ c , 161 where ε 1. We refer generically to the model class based on a perturbation approach as the "1 + ε model". Indeed, trend 162 lines in hydraulic geometry scaling of alluvial rivers follow predictions of the threshold channel theory (Fig. 3a), but with an 163 offset that indicates a formative fluid stress that is above threshold 44 . 164 Parker's 50 original 1 + ε model built directly on the hydraulic stable canal theory, and assumed ideal conditions including: 165 a straight channel, constant imposed discharge and C f , and uniform grain size along the bed and banks. It was formulated for 166 gravel rivers, in which sediment moves purely by bed load. Parker proposed that lateral diffusion of momentum, from the 167 channel center toward the margins due to turbulence, is the perturbation that solves the stable channel paradox. The solution 168 describes a channel with stable banks (τ ≤ τ c ), and active sediment transport in the channel center (τ > τ c ). This model predicts 169 a width-averaged formative shear stress τ ≈ 1.2τ c , i.e., ε ≈ 0.2. It is important to note that the value for ε depends on specific 170 model choices, such as the turbulent closure scheme and flow resistance relation. All reasonable choices, however, would 171 produce a near-threshold channel.
In Parker's formulation, channel geometry is fundamentally determined by hydraulics. An alternative model was recently proposed by researchers from the Institut de Physique du Globe de Paris (IPGP), in which channel shape, under laminar flow conditions, is adjusted to achieve a balance between lateral diffusion of bed-load flux toward the channel margins, and inward sediment motion due to gravity 51 . In this formulation, raising the imposed sediment discharge drives increases in channel 176 aspect ratio (W /H) and slope, away from the threshold state associated with no sediment flux (Fig. 3e). Experiments show, 177 however, that a channel will not tolerate a large increase in sediment discharge; instead, a single channel destabilizes into 178 multiple near-threshold threads akin to a braided river (Fig. 3d) 53,[115][116][117][118] . In this manner, the threshold state is like the critical 179 angle of a sandpile 119 : alluvial rivers can adjust their slope and channel geometry when driven by an imposed sediment load, 180 but they always remain close (i.e., 1 + ε) to the threshold state.

181
The two 1 + ε models above propose different perturbations to the threshold state to explain alluvial rivers. However, both 182 models are founded on the observation that an alluvial river adjusts its cross section to a state in which bed-load transport is 183 accommodated in the active channel center where fluid stress slightly exceeds threshold, and diminishes laterally to zero on the 184 banks. We believe that these models complement, rather than contradict, each other; Parker's approach neglects lateral bed-load 185 flux arising from a concentration gradient in sediment transport, while the IPGP model 51 neglects lateral fluid-momentum flux 186 arising from a gradient in stress. A more complete solution, which includes both effects, will likely have greater power for 187 predicting alluvial channel shape 51 . Nevertheless, the near-threshold constraint that τ = (1 + ε) τ c is sufficient to close the set 188 of governing equations for a first-order model of channel geometry. As we shall see, this model has surprising explanatory 189 power when applied judiciously.

190
Modifications and generalizations of near-threshold models 191 One person's boundary condition is another person's model. The near-threshold models above typically impose the following 192 variables as fixed boundary conditions: grain size, sediment discharge, threshold-fluid stress and flow resistance (among 193 others). In natural rivers, however, all of these parameters -where measured -can and do adjust to achieve a channel geometry 194 consistent with the near-threshold state. Here we briefly summarize relevant studies that explicitly examine these adjustments,

196
The fluid entrainment threshold is typically described by the dimensionless Shields stress, representing the ratio of fluid 197 forces over the submerged weight of a particle: τ * c = τ c /((ρ s − ρ)gD 50 ). For loose and non-cohesive grains, τ * c is primarily a 198 function of near-bed turbulence and its mean value may be estimated from the Shields curve 120, 121 . Variation in flow resistance 199 can result in apparent changes in τ * c , if an appropriate form drag correction is not applied when estimating the boundary channel" model: it posits that alluvial rivers adjust their geometry to the threshold fluid entrainment stress of the most resistant 211 material lining the channel, i.e., τ = (1 + ε)τ c max . In practice, gravel-bed rivers are adjusted to τ c of the gravel bed, while 212 sand-bed rivers are adjusted to τ c of the muddy banks (when present). This empirically validated model explains how sand-bed 213 rivers maintain stable banks, even though boundary shear stresses are far in excess of τ c for bed material (Fig. 4).

214
The importance of flow resistance, in terms of channel hydraulics and sediment transport, has long been recognized 45, 47 .

215
The boundary stress available to transport sediment is only a fraction of the total fluid stress; the rest, termed form drag, is 216 dissipated by turbulence arising from channel roughness at all scales -from grain, to bed form, to bank curvature 143, 144 .

229
A broader and more pervasive class of models, based on "extremal hypotheses", has been proffered as the primary alternative 230 to near-threshold models for explaining hydraulic geometry scaling. There is some physical basis for proposing an extremal 231 hypothesis as a closure scheme: often in problems that can be cast in terms of conservation of energy, there is a unique system 232 configuration that minimizes energy or maximizes entropy 145, 146 . In classical physics problems, this configuration may be 233 formally derived from a well-posed mechanical or thermodynamic constraint 147 . For rivers, the entrainment threshold is one 234 such mechanical constraint; yet, models that invoke extremal hypotheses do not formally apply this constraint. Researchers have predictions, and observed hydraulic geometry, appear to be due to error arising from these two issues 66 .

246
The importance of parameterization 247 We first consider gravel-bed rivers; based on the threshold-limited channel model, we assume that bank composition may be 248 neglected to first order 54 . The bankfull Shields stress (τ * b f ) values in the global database cluster around τ c predicted using 249 the Shields curve; the scatter around this trend, however, is more than an order of magnitude (Fig. 4a). These data would 250 appear, at first blush, to suggest that some gravel-bed rivers sustain bankfull shear stresses of almost 10τ c -conditions for 251 which bed material could be suspended -while others fall below the entrainment threshold at bankfull. Hydraulic geometry 252 scaling is correctly predicted by the 1 + ε model, but with similarly large scatter around the trend (Fig. 3ab). There is mounting 253 evidence 6, 66 that these discrepancies arise primarily from mis-estimates of τ c . Determining the threshold entrainment stress where the appropriate τ c could be measured or estimated, however, observed τ b f and hydraulic geometry scaling of sand-bed 271 rivers are in remarkably good agreement with predictions of the 1 + ε model 52 (Fig. 4ab). the sensitive dependence of turbulent momentum dissipation on complex boundaries is of fundamental importance for river hydraulics -but is also clearly beyond the scope of a first-order model for hydraulic geometry.

283
One question that arises in the application of the 1 + ε model is whether channel slope is an input parameter or a model 284 output. Both the Parker 50 and IPGP 51 models for gravel-bed rivers derive stationary solutions for channel slope, width and 285 depth as functions of water and sediment discharge. However, solutions for width and depth can be rearranged to depend only 286 on hydraulic factors -and not sediment discharge -if channel slope is imposed as an input parameter (Box 2). Hydraulic 287 geometry data show that width and depth are well predicted by hydraulic conveyance, while the large scatter in slope (Fig. 3c) 288 suggests an additional unmeasured factor -presumably sediment discharge -is required. Another possible factor is time, 289 which of course is neglected in stationary solutions. Sediment transport models, that couple channel geometry to long-profile 290 evolution via sediment mass conservation, predict that the timescale of slope adjustment may be on the order of millenia -291 much larger than the decadal timescales of width and depth adjustment 5, 77, 111, 157, 158 . This separation of scales suggests that 292 slopes of many natural rivers are not stationary; i.e., they may still be adjusting to modern water and sediment loads. This 293 change may be slow enough, however, to be considered quasi-steady in terms of hydraulic geometry; i.e., width and depth 294 may adjust in lockstep with changes in slope. Practically, this means that on engineering timescales slope should be treated as 295 an input parameter to the 1 + ε model 52 ; it is certainly easier to measure than sediment discharge. On geologic timescales, 296 however, alluvial rivers set their own slope through regrading of valleys and meandering.

297
The importance of averaging 298 Given a constant imposed water discharge above the entrainment threshold, a channel will develop a (statistically) stationary 299 geometry that just contains this flow 51, 114, 159 . Natural rivers, however, experience a wide range of discharges; most are well 300 below bankfull, while occasional floods can be well above 6 . This raises a fundamental question: is bankfull discharge merely a 301 useful reference point for hydraulic geometry comparisons, or is it a channel-forming flow condition with physical significance?

302
The seminal work of Wolman and Miller 30 provided an elegant conceptual framework for answering this question. They 303 reasoned that channels are adjusted to the flow of 'maximum geomorphic work': the stress whose product of frequency of 304 occurrence, and intensity of sediment transport, moves the most sediment in the long-time limit. Large floods have high 305 transport intensity but low frequency, while frequent low flows that do not exceed threshold do no work in moving material; it 306 is intermediate stresses, with low transport intensity and moderate frequency, that do the most work in shaping the channel.

307
Empirically, it has been demonstrated that Q b f also corresponds to the stress of maximum geomorphic work; in other words, 308 the bankfull flow indeed appears to generally be the channel-forming discharge 6, 160-162 . How is this achieved? To answer 309 this question one must understand how water discharge is converted into boundary shear stress (Box 1). Discharge may be discharge; they follow a thin-tailed distribution often well-described by an exponential function whose average value coincides with the bankfull discharge 6, 52, 165 . This occurs because the boundary stress that results from an imposed water discharge is 319 determined by channel shape and flow resistance; i.e., Q is imposed by watershed hydrology, but τ is an intrinsic property of the 320 channel. For flow within the channel (Q < Q b f ) we expect that flow depth, and hence τ, increases nearly linearly with Q. Once

321
Q exceeds Q b f , however, flow spreads across the floodplain and τ increases much more slowly with Q (Box 1). This results in 322 a rapid decline in the frequency of high stresses as flows exceed bankfull. The threshold constraint on channel organization is 323 central here: river banks destabilize, and widen the channel, with increases in boundary stress above threshold -producing a 324 negative feedback that keeps the channel in a near-threshold state 159 . The transformation of widely varying discharges into 325 a common thin-tailed distribution of excess shear stresses has been termed the 'critical filter' 6 . It is a logical consequence 326 of the organization of alluvial rivers to a near-threshold state, and justifies the use of a single bankfull discharge value in the 327 application of the 1 + ε model for hydraulic geometry.

328
The above should not be interpreted to mean that rivers do not respond to flows larger or smaller than bankfull, or experience

341
We now consider averaging over spatial variability in channel morphology. Dunes, bars and meander bends create systematic 342 variations in channel width, depth, slope and grain size -variations absent within a first-order hydraulic geometry model.

343
The length scales of these features should inform the spatial scales required for averaging 52, 168 . Despite these sources of 344 variability, a first-order model of channel geometry can still provide useful information. For example, measured channel widths 345 of a meandering river exhibited a wide statistical distribution, but the modal value was well predicted by the 1 + ε model 52 . complexities, the question of channel adjustment must first be addressed empirically. In recent decades, the same USGS gage 369 data discussed above has begun to be utilized to examine changes in cross-sectional geometry and hydraulic conveyance ( fig.   370 5a-c) 107, 182, 183 . Although not intended for this purpose, long-term gage records -extending decades and, in some cases, over 371 a century -are the best available data for relating channel size to hydroclimate 184, 185 .

372
A general observation is that, to first order, the hydraulic geometry of alluvial rivers is more-or-less adjusted to modern 373 hydroclimate regimes. This can be inferred from the tight scaling relations between bankfull channel width, depth and 374 discharge (Fig. 1e), and the general success of the stationary 1 + ε model -which involves averaging timescales on the order 375 of a decade or longer 6 . This result implies that statistically significant changes in hydroclimate -i.e., the frequency and 376 magnitude of discharge events -may be expected to result in detectable changes in hydraulic geometry (Fig. 5 a-d). Indeed, 377 multi-decadal trends in river channel form are widespread 108 . An analysis of almost a thousand, minimally-disturbed USGS 378 gauges across the USA found that 2/3 of the stations displayed significant temporal trends in riverbed elevation (Fig. 5), with 379 disproportionately higher rates of change in drier regions 107 . Results imply that at least some of these trends may be attributable 380 to anthropogenic-driven climate change, although no formal attribution analysis has yet been performed.

381
The sensitivity of alluvial river geometry to climate change is only just beginning to be explored. Hydroclimatic variability -along channel banks, on bars, and within the upstream catchment -could contribute significantly both in changing the 393 entrainment threshold, and also by potentially introducing lags and hysteresis in channel response 193,194 .

394
The considerations above offer some reassurance regarding the utility of the first-order 1 + ε model for interpreting observed 395 changes in mean hydraulic geometry as a consequence of climate change. They also, however, underscore a major shortcoming: 396 the assumptions of stationary water and weather regimes are broken 195 and many channels may be adjusting to a changing 397 climate, climatic oscillations, significant shifts in land use, and/or water management (Fig. 5). The close agreement between 398 channel width and discharge (Fig. 1e) indicates that the 1 + ε model can be used to predict channel size following adjustment, 399 however the rates and modes of adjustment cannot be predicted from a stationary model. To move forward with an empirical 400 approach, the next logical step is to consider the information contained in the higher-order moments of channel geometry 401 data. In other words, one possibility may be to consider river morphology and hydroclimate as probability density functions 402 that evolve together over time. What would such an approach look like? The cross-sectional river width, for example, can stresses that line the valley bottoms (Fig. 6). Although the strong perturbations to hydrology and sediment supply associated 422 with mill production have been removed, the rivers have not returned to their original form. Rather, they are relatively deep and 423 narrow, meandering channels 25, 196 that appear to be pattern stable (Fig. 6a-c). We posit that both channel states, pre-European -the exhumed Pleistocene cobbles (Fig. 6d). Application of the threshold-limited model provides close predictions of the 429 modern channel width (Fig. 6d). This case study reveals how the history of a landscape, embedded in sedimentary deposits, 430 becomes a substrate that can exert a first-order control on channel geometry through the entrainment threshold. This idea 431 has practical consequences: currently, major restoration projects are underway in some of these Northeastern rivers, based 432 on the premise that channel geometry may be returned to pre-disturbance conditions by altering the substrate. More broadly, 433 dam removal projects are rapidly growing in number around the world, with the similar goal of returning rivers to a natural 434 state 198-201 . The ultimate success of these projects will also be a measure of the success of the threshold-limited channel model.

435
Understanding and predicting threshold 436 If the 1 + ε model is at least a sturdy vessel for encapsulating our current understanding of first-order channel patterns, it is 437 anchored to a shifting bottom: the entrainment threshold. For all intents and purposes, τ c of in-situ river sediments -whether 438 gravel beds, or cohesive river banks -cannot be predicted from existing models to better than a factor of ten 52, 125, 152 . Given

446
The final challenge is that many of the aforementioned factors influencing threshold are spatially and temporally variable.

447
As a result of these challenges, we advocate that researchers and practitioners collect site-specific, in-situ measurements 448 of τ c for the most resistant material. There are currently so few measurements of cohesive bank-toe material that no general 449 trends can be reported 52 ; but our hope is that recent methodological improvements 214 will allow for broader data collectiong.
. Velocity shows considerably more scatter than either width or depth data, but is also not independently measured. The fitted trend lines are remarkable in that they provide first order approximations for river hydraulic geometry across 10 orders of magnitude in discharge.

30/36
Figure 2. Schematic illustration of the proposed orders of channel behavior. a| Schematic channel cross sections representing the orders of channel behavior from left to right: a threshold channel, near-threshold channel, channel with morphology, and meandering or braided morphology. The red dashed box represents the 1st order near-threshold channel that is the focus of this review. Light red regions represent regions of the channel at or above the threshold of motion. b| Planform or map view of the channel cross sections. c| Photographs and satellite images of (from left to right) a grass lined canal in Sweden, a cobble river in Alaska, the Eel River with alternating bars in California, and meandering and braided rivers in Indiana (USA) and New Zealand, respectively. The photo of the canal is courtesy of B. Neilson, and the alternate bar, meandering and braided rivers are from Google Earth.

31/36
Figure 3. The width of natural and laboratory alluvial rivers follow near threshold predictions. a| Natural (1,581 gray points 110 ) and laboratory alluvial rivers 53 dimensionless width (W /D 50 ) and discharge scaling compared with threshold theory (red line represents a threshold channel with τ * c = 0.05 and C f = 0.1). Both data sets sit slightly offset from the threshold theory. The shaded area denotes uncertainty within the possible parameter estimates for a threshold channel. b| Dimensionless depth (H/D 50 ) against dimensionless discharge. Fine-grained rivers are significantly shallower than threshold theory predicts. c| River slope against dimensionless discharge. The threshold channel is less steep than coarse-grained rivers and significantly lower gradient than fine-grained rivers. d| Schematic of the evolution of a transient experimental channel illustrating the transition from threshold, to increasing sediment flux to the point of channel instability. The early experiments of Stebbings 115 illustrate the end member conditions of single thread alluvial channels and the importance of sediment flux. e| Recent experimental efforts under laminar flow conditions directly measure the influence of increasing sediment flux on channel geometry 51 . From left to right, under no sediment flux (q s = 0) and constant discharge the channel cross section nearly exactly matches the cosine prediction from threshold theory, increasing sediment flux (middle and right) drives a stark increase in channel aspect ratio (W/H) and a steepening of the channel banks. . Shields diagram and illustrations of the near-threshold and threshold limiting models. a| Variation of bankfull dimensionless shear stress (Shields stress) with dimensionless grain size. The compiled rivers create two clouds of data between coarse and fine-grained rivers (dotted vertical line is 2.5 mm). The coarse grained rivers cluster near the threshold of motion as defined by the Shields curve, while fine-grained rivers cluster about the average threshold for clay sand mixtures (τ c−cs , dashed diagonal line). b| Illustration of the threshold-limited model for the Mullica River (yellow circle), a fine-grained river with cohesive banks. The black line represents the surveyed cross section. The cohesive banks (mud) are near the threshold at the channel center (τ b f /τ c = 1.13), while the fine-grained bed material (sand) is well above its respective threshold of motion (τ b f /τ c = 17). c| Illustration of the near-threshold model for a surveyed cross section (black line) the Salmon River, a cobble lined river in rural Idaho. Surveyed points below the bankfull flow are shaded according to the bankfull transport capacity (average τ b f /τ c = 1.17). Shear stresses are computed via the depth slope product using the hydraulic radius due to the narrow channel aspect ratio for the Mullica River and using the local depth for the Salmon River. Figure 5. River channel size grows and shrinks in response to hydroclimatic cycles and persistent changes in hydroclimate. a| Daily mean discharge (Q) records from 1955-2020 for Little Cedar River near Ionia, Iowa, USA (Gage No. 05458000). At this gage the annual peak flow (gray circles) has increased over time (dashed trend line). Periodicity within the discharge record is correlated with the Arctic Oscillation 109 . b| Mean bed elevation measurements over time showing a gradual degradation of the bed. The rolling average (red line) highlights periods of persistent scour or fill relative to the start of the record (dashed black line). c| Changes in flood stage channel capacity (∆CC, m 3 /s) over time. Increases in discharge were accommodated by channel bed degradation resulting in increased channel capacity since the 1950s. The rolling average (black line) represents periods of increased/decreased channel capacity relative to the average stage-discharge rating curve. d| Observed change in discharge frequency between the initial (1955-65, light gray) and final (2010-20, dark gray) ten years of the record. The probability density functions (PDF) are represented by the kernal density estimates of the natural log-transformed discharge data and show an overall shift in the discharge distribution. e| Potential statistical changes within the discharge record (non-stationarity) as a result of changes in landuse or hydroclimate include a shift in the mean and/or the frequency of extreme values. These changes result in differing forcing scenarios for channel adjustment. f| Schematic showing increases in channel capacity (conveyance) through degradation, increased widening and/or declining roughness (i.e., τ b f ), and decreases through aggradation, narrowing and/or increasing roughness. Creek showing buried precolonial wetland sediment characterized by relic tree stumps and wetland vegetation. The increased mobile sediment following the breach of the dam resulted in rapid incision down to the coarser periglacial sediment below. Modern inset gravel bars are a result of current sediment mobility. c| Longitudinal profiles of the modern river bed (blue) and two profiles (black solid and red dashed lines) showing the elevation up and downstream of the milldam. d| Modeled bankfull width for Mountain Creek (red circle) and coarse-grained rivers (D 50 > 5 mm, blue points) from the data compilaiton 110 . The 1 + ε model provides an accurate prediction of the modern channel width based on the periglacial sediment diameter (D 50 = 68 mm) indicating that the current channel is well described by the near-threshold model. Modeled predictions follow from Box 2, with a coarse-grained river average C f =7 and τ c = τ b f /1.2. The misalignment at low and larger widths for the compilation is a consequence of the use of a single value for C f .