Does the Threshold of Sediment Motion Constrain the Width of an Incising Laboratory River?

A physically rational model for river width is critical to predict macroscopic landscape evolution driven by fluvial sediment transport. Growing evidence suggests that rivers widen until the stress exerted by the fluid on the bed surface is close to the critical entrainment stress of the bank material. In this study, we test the limits of this model as a closure assumption in dynamically evolving river systems. We consider a simple laboratory channel with a fixed water discharge, monodisperse bed material, no sediment supply, and an initial relief that was sufficiently large to guarantee a finite transport capacity. Over time, the transport capacity approaches zero through changes in channel morphology. Concurrent measurements of width and sediment load illustrate how theory developed for equilibrium rivers fails to predict transient channel processes. Similarity to first order trends in natural rivers suggest that many rivers could be far from equilibrium.

Over the last four decades, numerous studies have advanced the notion that river channels adjust their width W so the stress exerted by the fluid on the bed surface at the bankfull flow condition τ is close to the critical stress needed to mobilize material that makes up the banks and bed, τ c (Phillips et al., 2022, references therein). This basic reasoning is encapsulated by the "1 + ɛ" model, that is: where ɛ ≪ 1. Despite its broad appeal, this model lacked rigorous theoretical justification until recently due to the inherent difficulty of modeling hydraulics and sediment transport over arbitrarily sloping beds (e.g., Parker et al., 2003;Seminara et al., 2002). Popović et al. (2021) present a comprehensive analysis of this problem that ultimately supports the 1 + ɛ hypothesis. Although the detailed derivation and testing depends on the assumption that flow is laminar, the underlying physical reasoning is also thought to be valid for turbulent flows. Critically, their model describes a stable configuration where width is constant, there are no overbank fluxes, and the sediment load Q s is spatially uniform (∂Q s /∂x = 0, where x is the longitudinal position).
What remains unclear is how channels achieve this stable state following a perturbation. This is problematic because rivers are continuously responding to climatic, tectonic, and anthropogenic perturbations over a wide range of scales. Indeed, finite sediment loads are essentially an outcome of perturbation; the sediment load in any given reach integrates longitudinal gradients in transport rate over the contributing drainage area, precluding globally uniform transport conditions. Tight empirical scaling relationships suggests that rivers can quickly adjust their width in response to perturbations (Phillips et al., 2022), and it may be reasonable to neglect small longitudinal gradients in sediment load for the purpose of predicting channel width in some cases. In other words, rivers appear to be able to adjust their width much faster than they can regrade their longitudinal profile following a perturbation. However, rapidly aggrading (Kim & Jerolmack, 2008) and incising (Croissant et al., 2017) channels exhibit dynamic changes in width that are not explained by Equation 1 but directly influence locations and rates of geomorphic change. This observation prompts an important question: what are the limits of the 1 + ɛ model as a closure assumption in disequilibrium rivers? Specifically, is the modeling framework presented by Phillips et al. (2022) appropriate for predicting the evolution of longitudinal river profiles?
Presently, we investigate this question using a simple laboratory experiment. Our objective was to create an incising alluvial channel (∂Q s /∂x > 0) and document the transient process of adjustment toward the stable state (∂Q s /∂x = 0). To this end, we carved a straight channel in uniform cohesionless substrate and imposed a constant water discharge. The initial relief was sufficiently large to guarantee a finite transport capacity but no sediment was supplied at the inlet. Because water discharge and grain size are fixed, the system achieves a stable state through morphological changes that reduce the sediment transport capacity. Concurrent measurements of width (from overhead photos) and sediment discharge (from repeat scans of topography) provide a direct test of the questions outlined above.

Setup
The experiment was conducted in a 1.2 m wide by 2.9 m long rectangular stream table elevated over a standing basin of water ( Figure 1c). A wire mesh enclosure was placed at the upstream end and filled with gravel (D 50 ≈ 2 cm) to dissipate fluid energy. This enclosure spanned the entire width of the stream table, allowing the channel to self-select its inlet width. The stream table was then filled with 0.3 m 3 of sand with a measured density of ρ s = 2.66 g/cm 3 , and an approximately lognormal size distribution with a median diameter of D 50 = 0.4 mm, and base-2 logarithmic standard deviation of σ ϕ = 0.72. The sand was leveled to a constant thickness of 12 cm except at the downstream end where the bed surface sloped downward toward the edge of the stream table. Then, a straight, triangular channel with banks at the angle of repose was carved down the center of the experimental domain.
Water was pumped from the standing basin into a head tank, supplying a constant discharge of 0.16 L/s. The reported discharge is an average of three measurements obtained by filling a plastic container at the inlet for a known duration, weighing, and assuming a fluid density of 1 g/cm 3 . Turbulence conditions were transitional with estimated Reynolds numbers ranging from 800 to 1,400. After filling the gravel enclosure, water flowed through the channel to the outlet and spilled over the edge of the metal platform into the standing basin. In effect, this setup fixes the base level position. Time t is measured from the instant that the water reached the downstream end of the experimental domain.
Lateral channel migration was unrestricted except by the sides of the stream table. The channel briefly touches one side but never touches both sides simultaneously; we argue that width was never restricted by the dimensions of the stream table. Water discharge measured at the inlet and outlet differed by less than 5%, indicating that groundwater flow was negligible. Overhead photos and topography scans were collected throughout the experiment to quantify changes in channel width and sediment discharge. The overhead images were also used to ensure the water surface position in the head tank remained constant. The experiment was ended when there were no observable changes in overhead photos over a period of 24 hr. This took approximately 340 hr (14 days). Only one experiment was performed. Below, we provide an overview of the data acquisition and analysis methods used in this study. Detailed processing workflows and a video are available with the published data set (Ashley & Strom, 2022).

Measurements of Sediment Discharge
Sediment discharge is calculated from repeat scans of topography (Figure 1a) obtained using a SICK Ranger E50 3-D laser scanning system mounted to a moving cart (see also : Hoyal & Sheets, 2009;Hamilton et al., 2015). The scan interval was increased from 30 min to 8 hr throughout the experiment. Each scan comprises 9 overlapping passes with the laser system, producing between 3 and 9 measurements of bed elevation for every 1 × 1 mm region of the bed. Measurements were merged using a custom algorithm to obtain final elevation models at 1 cm resolution.
The average sediment discharge ( , 0, 1) past the longitudinal position x between two topography scans obtained at t 0 and t 1 is estimated by integrating the Exner equation (Paola & Voller, 2005) in space and time. Careful accounting leads to the following expression (Ashley & Strom, 2022): Here, V(x, t) is the volume of sediment stored upstream of x at time t. To plot temporal trends, we interpret the average sediment discharge ( , 0, 1) as a proxy for the instantaneous ensemble average at = √ 0 1 , the midpoint of t 0 and t 1 in log(t). In total, we report 676 measurements of sediment discharge corresponding to 26 different times and 26 different longitudinal positions.
We note that topography scans were not corrected for refraction through the water surface. This is acceptable because the potential magnitude of uncertainty in Q s (x, t) introduced through systematic biases in measured bed elevation is small.

Measurements of Channel Width
Channel width was estimated from overhead photos ( Figure 1b) collected at an interval that increased from 1 to 16 min throughout the experiment. No images were acquired from approximately t = 17.2 to t = 23.2 hr due to a camera malfunction. Images were also not acquired for approximately 6 min during each topography scan.
An automated channel identification algorithm was used to extract 2,680 width measurements (at 1 mm longitudinal resolution) each from 3,498 images. The main element of this algorithm is a logistic regression model (fit using three manually digitized channel masks) that predicts the probability that each pixel is part of the channel as a function of lightness-corrected red, green, and blue bands. Probabilities are predicted for each pixel, smoothed in space, and thresholded at p = 0.5. Finally, width is calculated by summing the total number of channel pixels at each longitudinal position and multiplying by the pixel width. This approach correctly differentiates between wet and dry pixels with 97% accuracy in the training data set. 68% of the predicted widths lie within a factor of 1.2 of the measured width, and all widths lie within a factor of 1.9 of the measured width. The algorithm was spot-checked to ensure results are reasonable outside the training data set. Finally, average widths are computed at a spatiotemporal resolution that matches the 676 reported sediment discharge measurements to enable a direct comparison of these quantities.

Phases of Adjustment
The data described above document transient river adjustment at an unprecedented level of detail. In this section, we integrate data and observations to identify key processes associated with morphologically mediated relaxation toward the threshold state. Our goal is to establish a phenomenological context for the comparison of data and theory presented in Section 5.
Measurements of sediment discharge (Figure 2, insets) and qualitative changes in channel behavior observed in timelapse images (Ashley & Strom, 2022) indicate that morphological evolution is fastest at the beginning of the experiment and then slows over time. We identify three distinct intervals, referred to throughout the remainder of this paper as "initial," "exponential," and "autogenic" phases of adjustment. Most (75%) of the morphological change occurs during the exponential phase which begins at t ≈ 1.25 hr and ends at t ≈ 90 hr. This phase is characterized by a highly regular pattern of spatial and temporal variability in the average sediment discharge, while the initial and autogenic phases are characterized by deviations from this pattern ( Figure 2). Approximately 15% of the total geomorphic change occurs during the initial phase and 10% occurs during the autogenic phase.
During the exponential phase, the sediment discharge Q s (x, t) at each cross-section (i.e., each x) decreases exponentially over time (Figure 2a, Exponential Phase), that is, is well described by: Here, ′ ( ) is a parameter that is analogous to an initial sediment discharge at x (see below) and T c is a characteristic exponential timescale. A single value of T c = 11 hr provides a good fit for every x, and the longitudinal profile of the sediment load maintains a constant shape, varying in uniform proportion to ′ ( ) . The following expression constrains this shape ( Figure 2b): Here, L is the length of the experimental domain. Thus, the sediment discharge at any x and t can be predicted from four parameters: ′ (0) , ′ ( ) , T c , and L. Note that a finite best-fit value of ′ (0) reflects a boundary effect characterized by enhanced bank erosion where flow accelerates at the inlet. This effect is visible as an incipient erosional wedge on the left side of Figure 1. Because the boundary effect is relatively small, assuming ′ (0) = 0 provides a good fit except near the upstream boundary (Figure 2b).
Prior to t = 1.25 hr, measured sediment discharges are significantly higher than predicted from the best-fit exponential function (Figure 2a, Initial Phase). Consequently, ′ ( ) is not the true sediment discharge at t = 0; rather, it is a parameter that scales the sediment discharge during the exponential phase.
We suggest that the spatially uniform exponential timescale is indicative of a quasi-stable condition where local state variables are dictated entirely by the global configuration of the system. In other words, the parameters ′ ( ) , and T c are governed by global boundary conditions like Q w , D, L, and the initial relief. Recognizing that the initial channel was straight and had a triangular shape with banks close to the angle of repose, we believe the enhanced transport during the initial phase reflects rapid morphological changes that lead to this quasi-stable configuration. This conclusion is supported by rapid changes in width in the first 20 min of the experiment (Figure 3) and increases in sinuosity observed in timelapse images. Although spatially uniform exponential relaxation is not predicted by existing physical theory, we note that the total average longitudinal flux Q s /1.2 m (i.e., averaging over the width of the experimental domain including the channel and floodplain) follows a solution for one-dimensional landscape evolution driven by linear diffusion. Our simple fluvial landscape exhibits rich similarity to diffusion-dominated landscapes like hillslopes, obliquely supporting a hypothesis proposed by Reitz et al. (2014).
Significant departures from exponential relaxation are observed after t = 90 hr. During this time, autogenic pulses of incision are separated by periods of relative stasis. Broadly, we believe this behavior is an outcome of nonlinearity in transport rate as a function of shear stress close to the threshold of sediment motion. Departures may also reflect a failure of the assumption that short temporal averages are a good proxy for ensemble averages. The overall effect is that intermittent transport processes appear to persist for much longer than expected given the exponential trend observed during the previous phase. . Note that q s0 is the sediment flux corresponding to ɛ = 0.2 and is computed after Wong and Parker (2006) for our experiment. This is appropriate because flow is turbulent and inferred stresses cannot suspend bed material. (c) Comparison of field (Dunne & Jerolmack, 2018) and laboratory (Métivier et al., 2017) data with Equation 12.

First-Order Theory of Alluvial River Width
In this section, we propose a model for channel width that encapsulates recent theoretical advances proposed by (Popović et al., 2021). This model is identical to the 1 + ɛ model except at very low sediment discharges. Although the theoretical development of this model is based on an assumption of longitudinal equilibrium (∂Q s /∂x = 0), it is common to assume it is appropriate for predicting the evolution of disequilibrium river profiles because longitudinal gradients in transport rate are typically very small (e.g., Phillips et al., 2022;Wickert & Schildgen, 2019). Our primary objective is to test whether this assumption is reasonable for our disequilibrium experimental channel.
The concept of channel "stability" implies topographic invariance through time. Presently, we suggest that this conceptual definition includes uniformly aggrading and eroding landscapes as long as the rate of change of bed elevation η(x, y, t) is spatially uniform. Formally, a stable landscape satisfies: where x and y are spatial coordinates, t is time, and Ψ is the rate of bed elevation change. Neglecting uplift and subsidence, this condition is met when the divergence of the sediment flux is spatially uniform (Paola & Voller, 2005): Here, q s,x is the longitudinal component of the flux, q s,y is the lateral component of the flux, p is the bed porosity and Ψ is a spatially uniform rate of bed elevation change ∂η/∂t. To predict channel width from physical theory, researchers seek solutions to Equation 6 that incorporate simplified models for hydraulics and sediment transport.
Recognizing that all rivers exhibit stochastic fluctuations in reach-averaged state variables like width, it is often implicitly assumed that such solutions describe a deterministic expectation associated with a probabilistic ensemble of possible channel geometries conditional on appropriately defined boundary conditions.
The "threshold" channel is a trivial solution to this problem where q s,x = q s,y = Ψ = 0. It follows that channel geometry balances fluid, friction, and gravitational forces acting on particles resting on the bed surface. Neglecting lateral diffusion of fluid momentum and assuming (a) the water discharge discharge Q w is fixed, (b) the bed material is uniform and cohesionless with a representative diameter D and (c) all of the flow resistance comes from grains, granular forces are balanced across the entire channel if the cross-sectional shape h(y) follows a cosine function with a constant aspect ratio of where W is the width, H is the mean flow depth, and μ ≈ 0.7 is the critical friction angle of the bed material, that is (Devauchelle et al., 2011;Glover & Florey, 1951;Savenije, 2003;Seizilles et al., 2013): To accommodate a specified water discharge, width and slope S must satisfy and Here, the subscript 0 denotes values of W and S predicted for the threshold channel, * = ∕ √ 5 is a dimensionless water discharge, τ * c = τ c /ρgRD is a dimensionless critical stress for sediment motion, = √ ∕ is a flow resistance factor, U is the mean flow velocity, R is the submerged specific gravity of the sediment, g is gravitational acceleration, and α = (8/π) 1/4 . While the precise value of α depends on the flow model, defensible choices evaluate to α = O(1). Throughout this paper, we assume fixed values of C f = 0.1, R = 1.66, g = 9.81 m/s 2 , and τ * c = 0.04 such that W 0 and S 0 depend only on Q w and D.
The threshold solution was recently extended to include laminar channels with finite sediment loads by incorporating lateral diffusion of fluid momentum and sediment and assuming that the sediment load is constant in the longitudinal direction (∂Q s /∂x = 0) (Popović et al., 2021). This implies that Ψ = 0 and q s,y = 0 on the channel margins (i.e., there is no erosion or aggradation and no overbank fluxes) and allows Equation 6 to be expressed as an ordinary differential equation that can be solved numerically. In this case, diffusive fluxes of sediment and fluid momentum exert first-order control on the dimensionless excess stress ɛ = τ/τ c − 1 which, in turn, scales the width averaged sediment flux q s = Q s /W per an appropriate sediment transport formula (e.g., Wong & Parker, 2006). Mathematical, numerical, and experimental tests reveal two distinct scaling regimes that are delineated by Q s . At low Q s , width is effectively constant and changes in sediment supply are accommodated through changes in ɛ. At high Q s , ɛ, and q s saturate at fixed values ɛ 0 and q s0 , and changes in sediment supply are accommodated through changes in width. Popović et al. (2021) find that ɛ 0 ≈ 0.2, supporting the hypothesis that rivers cannot maintain large excess stresses (i.e., ɛ ≪ 1 as proposed by Parker, 1978). We propose the following expression to approximates this behavior over the full range of sediment loads: Here Q s * = Q s /[q s0 W 0 ] and k is a parameter that scales the smoothness of the transition between regimes at Q s * = 1. The fit of this expression is shown in Figure 3b against laminar channel data reported by Abramian et al. (2020). We find that k = 5 provides a reasonable fit, though the value of k may vary in different settings. This expression is identical to the 1 + ɛ model except at very low transport rates (Q s * < 1). We emphasize that the laminar flow assumption simplifies their mathematical derivation; however, the underlying physical reasoning is generally thought to be valid for turbulent flows.
The relationship between width and slope can be predicted by substituting a sediment transport formula into Equation 11. Then, the normal flow assumption and a flow resistance equation are used to relate the stress to the discharge and slope. In effect, the result is identical to the 1 + ɛ model except when S ≈ S 0 . Neglecting nonlinearity close to S 0 leads to (Dunne & Jerolmack, 2020): In other words, excess width is expected to increase in proportion to excess slope. This expression does not depend on the sediment transport formula used as long as the dimensionless sediment discharge ∕ √ 3 is expressed as a monotonic function of ɛ (e.g., Wong & Parker, 2006).
In summary, Popović et al. (2021) show that the 1 + ɛ model constrains the relationship between W and Q s when ∂Q s /∂x = 0. This relationship is captured by Equation 11 and leads to Equation 12 under normal flow conditions. A key question is whether this relationship remains valid when ∂Q s /∂x ≠ 0. This question is crucial because longitudinal gradients in transport rate are explicitly tied to macroscopic landscape evolution. Below, we present a direct empirical test of Equation 11 using concurrent measurements of Q and W in our incising laboratory channel. Then, we reinterpret empirical trends in a global database of alluvial rivers (Dunne & Jerolmack, 2018) in light of our experimental results and Equation 12.

First-Order Trend
Our simplified approximation of the model presented by Popović et al. (2021) (Equation 11), predicts that widths should increase in proportion to Q s when Q s > q s0 W 0 . For the specified water discharge and assumed values of fixed parameters, Equation 4 predicts W 0 = 9.1 cm. Throughout the experiment, measurements of Q s vary by a factor of O(10 5 ), leading to predicted widths that are as large as O(10 2 ) × W 0 . However, measured widths remain within a factor of O(10 1 ) of W 0 (Figure 3a). Because widths do not increase with Q s * as expected when Q s * > 1 (Figure 3b), we conclude that the 1 + ɛ model and the model of Popović et al. (2021) are insufficient to predict channel width in disequilibrium river systems with high sediment loads.

Fluctuations
High resolution measurements of width (Figure 3a) reveal deviations from the first-order trend described above. Width fluctuates rapidly, and though the magnitude of fluctuations increases in the downstream direction, every cross-section periodically returns to W 0 . The upper limit of fluctuations generally increases in the downstream direction and the frequency of fluctuations decreases over time.
In plan view, fluctuations are reminiscent of periodic oscillations between channelized transport and sheet flow described in an aggrading fan delta by Kim and Jerolmack (2008). These authors document coupled changes in channel width and slope that are linked to autogenic pulses of proximal and distal deposition and argue that periods of narrowing are initiated by a lateral instability mechanism that is analogous to the channelization instability on hillslopes (Izumi & Parker, 2000;Loewenherz-Lawrence, 1994;Schorghofer et al., 2004). Notably, they identify a clear link between vertical channel motion and width; periods of narrowing are coupled with periods of incision and sediment delivery to the delta toe.

Comparison With Field Data and Other Laboratory Rivers
Our results mirror trends observed in bankfull alluvial rivers (Figure 3c). Dunne and Jerolmack (2018) show that rivers remain close to W 0 across a wide range of conditions and settings, even when they are significantly steeper than the S 0 . Additionally, the range of W/W 0 values observed in our experiment is comparable to the range of values observed in bankfull alluvial rivers on Earth's surface. This result is significant because hydraulic considerations require Thus, if W ≈ W 0 and S ≫ S 0 , then τ ≫ τ c . This relationship is highlighted by isocontours of constant τ/τ c in Figure 3c, and may be derived for a specified Q w and C f assuming normal flow conditions (i.e., τ/τ c = HS/H 0 S 0 , where H 0 = W 0 2μ/π 2 per Equation 7). For reference, the same isocontours are plotted in Figure 3b assuming q s = f(τ) per Wong and Parker (2006). Note that these isocontours are plotted for illustrative purposes; the essential outcome is not sensitive to the choice of bedload transport formula used here.
Several mechanisms have been proposed to explain this trend including enhanced bank strength (e.g., Dunne & Jerolmack, 2018, references therein) and stress partitioning (Francalanci et al., 2020). While these are sufficient to explain elevated stresses with respect to τ c , they are not predictive because they introduce additional model parameters that preclude mathematical closure (e.g., the critical stress of the bank material). Our experiment reproduces this first-order trends observed in bankfull alluvial rivers using uniform sand, evincing a fundamental organizing principle that (a) does not depend on secondary factors like fine cohesive sediment or vegetation to increase the fluid entrainment stress of the bed material (b) is valid even when ∂Q s /∂x ≠ 0.

Discussion and Conclusions
This paper investigates dynamic fluctuations in channel width and sediment discharge throughout the process of incisional relaxation toward the threshold state. After an initial period of rapid adjustment, the sediment discharge decays over a spatially uniform exponential timescale. At long timescales, this pattern breaks down and sediment transport occurs primarily through intermittent, autogenic events.
Throughout the experiment, width remains close to the width of the threshold channel (within a factor of 10) across a wide range of sediment loads. As a result, changes in sediment supply are accommodated primarily through changes in the unit sediment discharge rather than changes in width. This result illustrates an important limitation of existing physical theory derived assuming ∂Q s /∂x = 0 (e.g., Phillips et al., 2022;Popović et al., 2021).
Recognizing that incision and narrowing appear to be related (Croissant et al., 2017;Kim & Jerolmack, 2008), we suggest that vertical channel motion associated with longitudinal gradients in sediment load may play an important role in maintaining large excess stresses in disequilibrium river systems.
Our experimental river departs from the 1 + ɛ model at high sediment loads in a manner that is consistent with observations of bankfull alluvial rivers. This result is surprising because recent studies have attributed large excess stresses to enhanced bank strength. We propose an alternative explanation: many rivers are far from equilibrium. While bank strength may be needed to maintain large stresses in equilibrium channels (i.e., ∂Q s /∂x = 0), our results suggest that vertical channel motion allows rivers that are far from equilibrium to maintain large excess stresses while they approach the near-threshold state. This may explain the absence of very wide, shallow rivers predicted by the 1 + ɛ model when S ≫ S 0 .

Data Availability Statement
A GitHub repository containing all data and detailed processing workflows can be accessed at https://github. com/tashley/Relaxing-Channel-Data-and-Workflows. Static versions of the processing workflows are viewable without a GitHub account. This repository is also preserved with a persistent identifier through Zenodo (Ashley & Strom, 2022).