River-bed armouring as a granular segregation phenomenon

River bed-load transport is a kind of dense granular flow, and such flows are known to segregate grains. While gravel-river beds typically have an “armoured” layer of coarse grains on the surface, which acts to protect finer particles underneath from erosion, the contribution of granular physics to river-bed armouring has not yet been investigated. Here we examine these connections in a laboratory river with bimodal sediment size, by tracking the motion of particles from the surface to deep inside the bed, and find that armour develops by two distinct mechanisms. Bed-load transport in the near-surface layer drives rapid, shear rate-dependent advective segregation. Creeping grains beneath the bed-load layer give rise to slow but persistent diffusion-dominated segregation. We verify these findings with a continuum phenomenological model and discrete element method simulations. Our experiments suggest that some river-bed armouring may be due to granular segregation from below—rather than fluid-driven sorting from above—while also providing new insights on the mechanics of segregation that are relevant to a wide range of granular flows.


INTRODUCTION
River-bed grain size controls the exchange of solutes, nutrients and fine particulates across the sediment-fluid interface [1,2], and determines the flood magnitude required to initiate motion [3][4][5].Grain size, however, also evolves over a series of floods as particles are sorted longitudinally and vertically during transport [3,[6][7][8].A ubiquitous pattern observed in gravel-bed rivers is armoring, in which the median grain size of the surface is significantly larger than that of the subsurface (Fig. 1 A).Laboratory experiments designed to simulate gravel rivers -i.e., bed-load transport of heterogeneous grain sizes -reproduce the phenomenon, but disagree on its origins.Three competing mechanisms have been proposed: (i) kinetic sieving, in which smaller particles migrate downward through the void spaces between larger particles during motion [9]; (ii) 'equal mobility', whereby the proportion of large and small surficial grains adjusts to achieve a spatially constant entrainment stress [6]; and (iii) sediment supply imbalance, in which the transport capacity of the flow locally exceeds the upstream supply and results in surface coarsening [7].All of them assume that gravel in transport only mixes with the substrate over a small "active layer" that is one to several grain diameters deep.
Recently, sediment transport experiments have revealed that granular motion extends much deeper into the subsurface [10,11].In particular, grains transition continuously from rapid bed-load motion at the surface to slow creeping motion far below the surface [11,12].
Both kinds of motion also occur in dry granular systems, where bed-load corresponds to a dense granular flow, and creep is characteristic of quasi-static deformation of disordered granular packs [11,12].The former is known to produce robust vertical size segregation by kinetic sieving [13][14][15].Phenomenological continuum models based on this premise [15][16][17][18] produce vertical segregation that is consistent with experimental observations [16,19] and discrete element method (DEM) simulations [20,21].Segregation by creep is unexplored; while reports of slow coarsening do exist [22], its connection to creep has not been demonstrated.
The contribution of granular physics to river-bed armoring has only begun to be examined.Frey and Church [9] showed with laboratory experiments that bed load drives segregation by kinetic sieving that is qualitatively similar to dense granular flows.Here we investigate granular segregation and quantify its contribution to armoring using an ideal- computed as τ = ηU f /h f [11] where U f and h f are the top-plate speed and flow depth, respectively.
The red curve shows the long-term-averaged streamwise particle velocity u x (z), where I and II correspond to the bed load and creep zones, respectively.The directions x and z are indicated.
ized laboratory river experiment.Our setup [Fig.1B] is designed to: eliminate the disruptive influence of flume boundaries by using an annulus; image particle motion from the sedimentfluid interface to deep in the subsurface, away from the wall; isolate granular contributions by simplifying particles to bimodal spheres and eliminating fluid turbulence with a viscous fluid; and explore a range of transport conditions from near threshold to vigorous bed load.
These experiments demonstrate how river-bed armor can develop due to bottom-up motion of subsurface grains, while revealing new insight on granular segregation mechanisms in a system where rapid and slow granular flows co-exist.Results are compared with predictions from a modified phenomenological segregation model, and with DEM simulations of a dry-granular bed under shear.

EXPERIMENTS
Experiments were conducted in a closed-top annular flume (Fig. 1B); details of the apparatus have been described previously [11,12].The channel walls are smooth to allow slip between grains and the boundary, in order to approximate an infinitely deep and wide channel.The flume is filled with a bidisperse granular bed of acrylic spherical grains with small and large diameters, d s = 1.5 mm and d l = 3.0 mm, respectively, and density ρ p = 1.19 g mL −1 ; the ratio of total small to large grain volume in the channel is V small V large = 2.The grains are submerged in a fluid of viscosity η = 72.2mPa s and density ρ = 1.05 g mL −1 .A fluid gap is sheared from above by rotating the lid of the flume to apply a constant fluid-boundary shear stress, τ (Fig. 1C); it is reported here as dimensionless Shields number for the small grains, defined as τ * s = τ (ρp−ρ)gds , where g is gravity.The associated Shields stress for large grains as τ * l = ds d l τ * s .For reference, Shields numbers for each experiment are compared to the critical Shields number, τ * c , that is classically used to identify the onset of sediment transport.Our previously determined critical Shields number for a monodisperse bed of small grains, τ * cs 0.1 [11], is used here as the reference critical stress, recognizing that the actual value may differ in this bidisperse system [24].We determined empirically for the present experiments that the range of Shields numbers τ * cs ≤ τ * s ≤ 5 τ * cs corresponds to bedload transport: a thin surface layer of moving grains in frequent contact with, and supported by, an underlying granular bed that is creeping [11].We report data from experiments conducted at five Shields numbers, τ * s = [2.7,3.8, 4.1, 4.4, 4.7]τ * cs .All flows were laminar (Reynolds number ≤ 4) and grain collisions were viscously damped (Stokes number < 1) (see SI, section 1).
The bed at the start of each experiment was composed of sedimented particles forming an approximately flat granular bed (see Methods).At the beginning of an experiment (t = 0 s) fluid shear was initiated at the specified Shields stress, and applied for a duration of 24 hr or longer.We image a cross-section of particles from the bed surface (z s ) to the bottom of the channel through time using refractive-index matched scanning [25] (Figs.1B; S2).Vertical profiles of streamwise particle velocity (u x (z)) for experiments at all Shields stresses were determined from averaging pixel strips in the streamwise direction over all time, using image cross-correlation (SI, section 4).Velocity profiles confirm the existence of two distinct regions of particle motion (Fig, 1C).Zone I corresponds to bed load, where velocity decays rapidly with depth below the surface; below this is zone II associated with creep, characterized by a much slower decay [11].All runs show a qualitatively similar evolution of the bed through time: a coarse surface "armor" layer develops as large grains are delivered from below; first more rapidly by bed load, and then more slowly by creep (Fig. 1, SI Movies 1 & 2).This is explored in more detail below.In order to probe the size-and depth-dependent behavior of grain motion, and its contribution to vertical segregation, we construct trajectories of all imaged grains using the particle tracking method [11] (see Methods) for a representative experiment at Shields stress τ * s = 4.1 τ * cs .Profiles of average vertical velocity for large (u z,l ) and small (u z,s ) grains, computed from these trajectories, show a striking pattern: they are anti-correlated in the bed-load regime, with net upward (positive) velocity for large grains and net downward (negative) velocity for small grains.Although there are deviations in the near-surface (within 1d s of z s ) due to intermittent saltation, below this region the velocity of large grains decreases with depth and reaches approximately zero at the transition from bed load to creep (z c ) (Fig. 2A).The decay rate of u z,l is roughly exponential, and coincides with the decay of the bulk streamwise velocity u x (Fig. 2A-inset).This suggests that the observed vertical advection of larger grains is linked to horizontal granular shear in the bed-load zone.
Grains in the creep zone have a small but detectable vertical velocity.To determine the dominant modes of particle motion in bed load and creep, we inspect the scaling of the mean-square displacements (MSD) versus time.For the same experiment at τ * s = 4.1 τ * cs we compute the vertical MSD as a function of depth for the large grains as MSD(∆t) ≡ 2 over a duration of 20 minutes; the brackets indicate ensemble averaging over grains and the reference time t, and z is the particle's vertical position (Fig. 2; see Methods).A distinction can be made between grains above and below the depth associated with the transition from bed load to creep.Grains in the bed-load zone exhibit MSD growth at short times that approaches ballistic motion, and is consistent with the advection described earlier (Fig. 2A).The strength of the advection behavior diminishes at larger timescales where it perhaps transitions to super-diffusive behavior.In contrast, grains in the creep zone appear to exhibit caged dynamics in which MSD grows slowly or not at all at short timescales.Motion transitions toward diffusive and sometimes super-diffusive dynamics at longer times.The crossover timescale indicates the average lifetime of cages, and it increases with depth into the creeping zone.This behavior is similar to what has been observed in slow granular flows [26][27][28], and indicates that particle movement in creep is related to creation and destruction of the granular contact network [28,29].
To visualize the resulting development of surface armor, we examine the spatio-temporal concentration map of large grains; φ l represents the streamwise-averaged areal fraction at a given depth and time (see Methods).The development of surface armor is seen as a highconcentration surficial layer that thickens through time (Figs.2C; S4-S8).We quantify the thickness of the armor layer (Fig. 2 C) as z sa −z i , where z sa and z i are the position of the top and bottom surfaces that define the armor layer, respectively (SI, section 5) for all five Shields stress experiments.The data suggest the existence of two stages in the creation of armor, anticipated by the granular dynamics described above (Fig. 3 A).First is rapid segregation (duration of 10 2 -10 3 s), as large grains are delivered up from the shallow subsurface.The rate of segregation shows a strong dependence on the driving Shields number, consistent with shear-rate dependent segregation of bed load.Once the bed-load zone is depleted of large grains, there follows a slow but persistent segregation that continues for the duration of the run (∼ 24 hr).We interpret the slow stage of segregation as creep driven.Interestingly, the rate of segregation in this stage is insensitive to the driving Shields number, suggesting that creep segregation does not depend strongly on the driving shear rate.
Armor development in our experiments results from a vertical flux (z direction) of coarse grains toward the bed surface.We quantify this segregation flux, J, as the time derivative of the number of large grains in the armored layer: where A is the cross-sectional area of the armor interface in the x−y plane.The variation of segregation flux density (J/A) with time (Fig. 3B) clearly shows the existence of two stages of armor formation.We introduce a dimensionless time t/t adv , where the characteristic advection time t adv = h bl u z,l ∼ h bl aU sf ; u z,l , U sf and h bl are the average large grain velocity, the average surficial grain velocity and the thickness of the bed-load layer, respectively, and a = u z,l /U sf ∼ 10 −3 is measured for the experiment at τ * s = 4.1 τ * cs , but the value for a collapses all data later in Fig. 3C.We also define a dimensionless segregation flux J/J(0) where J(0) is the initial value for J at the start of each experiment.Utilizing the dimensionless time and flux variables produces a reasonable collapse of the data (Fig. 3C).
For all experiments J/J(0) decays to a value of 1/e at a characteristic dimensionless time of O(1).

ADVECTION-DIFFUSION SEGREGATION MODEL
Sediment transport produces armoring that appears similar to reported granular segregation experiments [13,16], implying that the presence of a viscous fluid has little influence beyond determining the shear rate of surficial grains.In particular, some previous experiments in dry granular flows suggested that segregation rate depends on the granular shear rate [30,31], consistent with our findings for bed load (although another study found otherwise [19]).In addition, a recent study found that particle diffusion was shear-rate dependent for rapid granular flows but independent of shear rate for creep [32], similar to our experiments.
Because the exact mechanism of segregation is still a subject of debate [18,21,33], there is no universally agreed upon continuum theory.Nonetheless, one-dimensional (1D) continuum models generally describe the vertical evolution of concentrations of binary mixtures through time with a phenomenological advection-diffusion equation [16,20].Here we develop and apply a modified version of one such model, the Gray-Thornton model [15,34].
The model requires specification of: vertical advection and diffusion coefficients, usually assumed to be constant [15]; the vertical granular velocity profile; and the initial concentration profile.It then solves for the temporal segregation of large and small grains subject to mass conservation constraints.Two new ingredients must be included to account for the granular dynamics observed in our experiments: (i) for the bed-load regime, both advection and diffusion depend on shear rate; and (ii) for the creep regime there is no advection, and diffusion is independent of shear rate [32].Our modified advection-diffusion segregation model, written in terms of the evolution of the large-grain concentration φ l , becomes: Equation 2 is written in terms of dimensionless elevation ẑ = z/H and time t = tU sf /L, where H and L are the height of the granular pack and the length of the centerline of the annular flume.The flux function F (φ l ) determines the dependence of the segregation flux (S r F (φ l )) on φ l .Although there are ongoing debates on the mathematical form of the flux function [34,35], we implement the simplest choice: a quadratic function that is symmetric about φ l = 0.5, which assumes that small and large grains behave identically but in opposite directions.The original Gray-Thornton model assumed a non-dimensional advective segregation velocity S r that is independent of shear rate.We introduce a depth-dependent parameter, S rn , in order to redistribute the non-dimensional advective segregation velocity, S r , according to the depth-dependent grain velocity, u x (ẑ), normalized by the vertical average of grain velocities, u x (ẑ) (Eq.3).
The form of the normalized velocity is determined by a fit to the bed-load velocity profile such that: where β is the exponential decay constant of the bed-load velocity profile.Accordingly, for our analysis we define the parameter S r = L H ux(ẑ) q, where q is the maximum bulk advective segregation velocity, i.e., that associated with the start of the experiment (t = 0; see SI section 6; Fig. S3).Similarly, we introduce a dimensionless and vertically-varying diffusivity D rn that has the same exponential decay as the velocity profile characterized by β (Eq.5).The parameter D r = DL H 2 ux(ẑ) is a non-dimensional diffusive-remixing constant, where D is the dimensional diffusivity.
To apply the new model Eq. 2 to our experiments requires specification of several parameters, determined from each experimental run (see Methods and SI).The input velocity profile u x (z) is determined by fitting two exponential functions to the time-averaged velocity profiles of the bed-load and creep zones, respectively (see Methods; Fig. S6).The input value for q is computed as the upward migration velocity of the center of mass for the large particles at the start of each experiment (see Methods).Note that the advective segregation term in Eq. 3 decays with decreasing velocity (and depth) in the bed-load zone, and is set to zero in the creep zone (z < z c ).The diffusivity also decays with velocity (and depth) in the bed-load zone, and is constant for creep (Fig. S6 E).We take the dimensionless diffusion constant D r as a fitting parameter.In particular, the ratio S r /D r is estimated by fitting the position of the armor interface through time for each experiment.We find that a constant ratio of S rn /D rn ∼ 318 for the bed-load zone and S rn /D rn = 0 for the creep zone is sufficient to describe the development of armor for all Shields numbers.We use the profile of φ l at the start of our experiments (t = 0) as the initial concentration profile for the continuum model (Fig. S6F).
A visual comparison of armor development for the example condition τ * s = 4.1τ * cs shows that the advection-diffusion segregation model (Eq.2) captures the experimental behavior well (Fig. 2D).A more quantitative comparison of the thickness of the armored layer through time (Fig. 3 A) demonstrates good agreement between the model and data for all Shields numbers.Importantly, the model correctly captures the initial fast and subsequent slow stages of segregation.The large ratio S rn /D rn ∼ 318 for z > z c confirms the idea that the rapid stage of armor development is driven by shear-rate dependent advection associated with bed load.The fact that the ratio S rn /D rn remains constant for all experiments suggests that the model results are robust.The bulk kinetics can be related to particle-scale advection and diffusion by noting that S rn /D rn = Sr Dr β = P e H d β, where the particle-scale ratio of advection to diffusion is given by the Peclet number P e = u z d l /D z .For the experiment with τ * s = 4.1τ * cs we determined from measurements that u z = 1.51 mm/s and D z = 3.38 mm 2 /s, which leads to P e = 1.3 and S rn /D rn = 140; the latter is the same order of magnitude as the ratio used in the continuum simulations.The creeping zone is characterized by a constant value for D r , and a lack of advection (S r = 0), for z < z c .This supports the notion that the slow stage of armoring results from diffusion by creeping grains that is independent of local shear rate.

DISCRETE ELEMENT MODELING
The analysis presented thus far shows how explicit accounting for the kinematics of granular motion in bed load and creep can produce a reasonable continuum description of armor development.In order to demonstrate that the observed armoring in experiments is entirely a consequence of granular physics, we now turn to DEM simulations in which the velocity profile and segregation dynamics arise spontaneously from grain-grain interactions.Simulations are performed with LIGGGHTS, an open-source granular modeling package based on LAMMPS (http://lammps.sandia.gov).Details of model implementation are available in Methods and Supplementary Information.In accord with the low Stokes number of our laboratory experiments, the restitution coefficient is chosen to be very small (e n = 0.01) such that collisions are highly damped (Table S2).Otherwise, there is no treatment of the viscous fluid in DEM simulations.
The model domain is constructed to have a geometry, grain size and size-volume ratio that are the same as the experimental setup (Fig. 4A).The system is driven by a layer of large grains deposited at the surface and moving at constant velocity u top in the x direction.
Simulations are run for a duration that is equivalent to ∼ 60 s due to their computational expense; however, we show below that this is sufficient time to observe the fundamental dynamics.Simulations show behavior that is qualitatively comparable to the fluid-driven experiments of armoring, confirming the existence of two stages of segregation (Fig. 4).First is fast segregation within the rapid granular flow regime (first few grain diameters from the surface).Then, once grains are depleted from this "bed load" zone (Fig. 4C), armoring transitions to a slow stage driven by creep from deeper layers.
For a more direct comparison, we examine the growth of armor thickness through time (see Methods) for the previous simulation and an additional run with u top = 0.05 (m/s) which corresponds to τ * s = 2.7τ * cs .For both runs the agreement of DEM simulations with the experiments is reasonably good (Fig. 4D).This agreement is especially encouraging given that: simulations neglect fluid flow entirely; the initial concentration of large grains in the experiments was difficult to control and not uniform; and there was no tuning or calibration done for the DEM runs, beyond adjusting the velocity of surface grains to match experiments.

DISCUSSION
Even though our flows were laminar, experiments and theory have shown that laminar bed load is similar to its turbulent counterpart in many respects [24,[36][37][38][39].Our results show armoring dynamics that are qualitatively similar to previous experiments [40,41] conducted under conditions more representative of gravel rivers -i.e., poly-disperse and natural-shaped particles (average grain diameter d ∼ 1 cm) in turbulent flows with driving stress τ * ≈ 2τ * c .Those studies [40,41] found a Shields-stress dependent armoring rate with a relatively rapid initial stage (a few hours) followed by slower stage.While data on particle motions were not reported, we can perform a scale analysis of the expected bed-load armoring timescale due to granular segregation, t adv ∼ h bl aU sf , by assuming: h bl = (3−5)d [42]; U sf ∼ 1 cm/s [43,44]; and our experimentally-determined value a ∼ 10 −3 .This analysis yields t adv ∼ [1 − 2] hrs, within the observed range of experiments [40,41], and may be a reasonable bed-load armoring timescale for natural gravel rivers.Translation to the field, however, may need to account for the presence of bed and bar forms that can influence armor formation [45].[7,40,41] attributed armor development to a lack of sediment supply to the channel, which they hypothesize resulted in winnowing of fines and concentration of coarse grains -in other words, sediment-supply imbalance.Our experiments, however, showed no significant size-selective transport at the surface and, more importantly, there were no supply limitations because the flume is annular.We can thus rule out sediment-supply imbalance for our experiments.Our results support the kinetic sieving model, on which the phenomenological Gray-Thornton equation is based.An important new finding, however, is that segregation does not occur only in the "active layer".If the bedload zone corresponds to the active layer, then the associated sorting is important but occurs rapidly.Creep delivers grains from far below the bed-load zone to the surface, contributing to persistent armor development that was not previously recognized.The agreement of DEM simulations and experiments confirm the contention of Frey and Church [9] that riverbed armoring is a granular segregation phenomenon, and suggest minimal influence of the fluid beyond determining the surface grain velocity.We point out that sediment-supply imbalance may still be important for armoring in some rivers; in particular, under sediment limitations such as downstream of dams where river beds experience net erosion that may preferentially remove finer grains [7].The granular segregation dynamics revealed here, however, would operate in all environments regardless of sediment supply and may therefore be more prevalent.Future experiments with more 'natural' flow and particle conditions, that control for sediment supply while also examining precise granular motion from the surface to the bottom of the granular pack, would be helpful for assessing the relative importance of these different mechanisms.A potential field confirmation of the armoring mechanism proposed here would be the observation of a zone underneath the armor layer that is almost entirely depleted of large grains (Figs.2C and D).Size-selective surficial transport would not be expected to influence the concentration of coarse grains beneath the armor layer.

Authors of previous experiments
Our work sheds new light on the mechanics of granular segregation.Experiments clearly show that vertical advection of large grains is shear-rate dependent.Explicit accounting of this dependence, and also of shear-rate dependent diffusion, is needed in order to explain observed segregation rates for the rapid granular flow regime (i.e., bed load).Moreover, data and models demonstrate that creep contributes to segregation, and that its mechanism is distinct from rapid granular flows.Large grains in the creep zone show no preference for upward-or downward-directed motion.Their long-time motion may be modeled as vertically-isotropic and constant diffusion.Short-time dynamics show that creeping grains are caged, and indicate that their motion is likely induced by long-range transmission of forces through the granular contact network [46,47].This may be why creep motion is independent of shear rate, at least for the range of Shields stresses examined here.It is intriguing that isotropic diffusion in creep can give rise to a net upward flux of large grains.
Based on our results, we hypothesize that this flux arises because large grains that cross the boundary into bed load are then advected to the surface.If correct, this implies that a purely creeping granular pack (no bed load) should not produce armoring.
The experimental and modeling results presented here are a first step in assessing the contribution of fast and slow particle motion to vertical segregation.Our sediment mixture was bi-disperse in order to establish connections between granular shear segregation and river-bed armoring, but many systems of interest (including rivers) have a polydisperse grain size distribution that may exhibit different behavior.Such a distribution would challenge the application of continuum models, but is amenable to experimental and DEM simulation approaches.River-bed armoring in our experiments and models was found to be driven by bottom-up granular segregation, rather than top-down surficial sorting driven by the fluid.Our findings show how information from the surface, in terms of fluid-driven shear, is transmitted deep into the subsurface through grain-grain interactions that are typically neglected in sediment transport models.Granular motion in the subsurface transmits information back to the surface through the delivery of coarse grains, linking surface dynamics to subsurface structure.By examining the river bed as a discrete medium, we were able to link the macroscopic pattern of armor development to the physics of sheared granular systems.
Our results add to a growing body of evidence that sediment transport systems belong to a broader class of granular flows [9,11,12,48], and show how examining geophysical flows through the lens of granular physics can reveal novel insights for both fields.

MATERIALS AND METHODS
Details of the experimental setup and shearing protocol are described in SI section 1.The bed preparation protocol was inspired by Golick and Daniels [30].Grains were initially deposited in an inverse-segregated state, with large grains at the bottom, and then subject to a driving stress equivalent to τ * s = 20τ * cs for ∼ 1 minute to fully suspend and mix the large and small populations.Fluid shear was halted and the suspension left for ∼ 30 minutes to allow sedimentation, relaxation and compaction of the granular bed to reach completion (Fig. S1).
Implementation of the continuum model is described in SI section 6 (Figs.S3-S8; Table S1).Details of the DEM model and parameters appear in SI section 7 (see Table S2).The local concentration of large grains is defined as, φ l (z) = A l x A l +As x where A l and A s are large and small grains areas, respectively, and • x indicates pixel-wise streamwise integration.
unidirectional assumption is justified based on the small ratio of radial viscous stress to the azimuthal viscous stress for our experimental conditions: where h f 3d s , R is the flume radius and c 0.06 is an estimated coefficient [36] that is only weakly dependent on the flow aspect ratio.

Detection of the bed surface
In order to detect the surface position, first the concentration profile C(z) for a given configuration of particles is determined from a processed binary image, valued at zero outside of particles and one inside of particles.For each elevation z, the concentration is determined as the pixel-wise average in the x direction.This concentration is the 1D analogue of packing fraction, the fraction of space occupied by the particles.The surface is defined as the position z s at which the concentration crosses fifty percent of its saturated value [11,49].We use a fixed threshold of 0.35 to define the surface position, as the saturated value does not vary significantly from experiment to experiment.We define z s after averaging the concentration for a ∆t = 100 s at the beginning of each experiment.The time duration is sufficiently long for the flux convergence time as observed in our earlier study [11].Slow granular compaction [50,51] and slow dilation due to segregation [30] approximately counterbalance such that the surface position remains almost constant as the armoring experiments progress.The bed surface position is used for calculating the Shields stress at each experiment.

Imaging technique and particle detection/tracking
Using a Nikon DSLR 5100 digital camera, we record the real-time positions of single particles by acquiring the fluorescence intensity from a laser dye (concentration ≈ 1µM ) dispersed in the fluid that is suitable for long data acquisition without significant photobleaching.To sample fast dynamics near the surface, where the relevant timescale is the settling time of particles over their own diameter d/v sed = 0.68s, we acquire images continuously at 24 Hz for 10 − 20 mins.To sample slow dynamics in the bed, we acquire single images at a rate of one every 15s for 24hr or longer.Figure S2A shows a sample raw image at the start of an experiment.To detect the positions of the particles with subpixel accuracy, we find particle positions to pixel accuracy by peak-finding above a threshold.The details of the background correction process are described in the supplementary materials of our previous publication [11].The process is designed specifically to handle both long-wavelength background intensity variations and intermediate wavelength fluctuations (stripes) due to slight mismatches in the index of refraction of the particles relative to the fluid.After removing the background, we determine sub-pixel positions using the radial symmetry method [52].
We use this method to identify radially symmetric features at different probing length scales.
Then, features are identified as peaks in the three-dimensional image (x, y, R), over a range of radii that is sampled linearly in log-radius.Peaks are obtained by two passes, first a pixel-scale search, then a local quadratic fit to obtain subpixel positions.A snapshot of detected particles with this method is superimposed on a gray-scale raw image and is shown in Figure S2B.The same detected particles are also shown in binary format in Figure S2C.
A fixed diameter threshold of d = 1.38mm is used for separating large and small particles in all experiments as shown in Figure S2D.Finally, a snapshot of identified large and small particles using this threshold is shown in Figure S2E.
We stitch positions at different frames into tracks using Lagrangian particle tracking [53].

Velocity profiles
For each experiment, a 10 minute video capture with frame rate 24 fps at the start of the experiments is converted and processed into consecutive binary images following the procedure described in imaging section above.The consecutive binary images, I(t) and I(t + ∆t) are then used as the input of pixel-wise cross-correlation analysis along the x direction at each pixel elevation z.The position of the central peak in the cross-correlation between I(t) and I(t + ∆t) corresponds to the average streamwise distance traveled by gains at elevation z during ∆t without regard to small and large particle species, i.e. for all particles.The result is averaged over the full duration of the video capture.This technique yields a time-averaged streamwise velocity profile u x (z) for all particles.Note, large particles are weighted more heavily than small particles.The results are in agreement with velocity profiles determined from the particle-tracking method for all experiments.For the case of Figure 2A, the velocity profile of large grains is computed using the particle tracking method described in the imaging technique section above.

Determination of armor thickness
The top surface of the armor layer, z sa , is characterized as the position where the streamwise averaged concentration of large grains equals φ l = 0.9.The interface (bottom) of the armor layer with the rest of the granular bed, z i , is calculated as the location where the gradient of φ l reaches a minimum below the surface.The surface and interface positions time-series are smoothed using a running average of temporal window size 8.33 s for images obtained from video capture conversions and temporal window size 833 s for image captures.These are shown with dashed lines in Fig. 2C.The thickness of the armored layer is defined as z sa − z i .

Implementation of continuum model
The variables used to compute the advection and diffusion parameters for each experiment are reported in Table [S1  The DEM model consists of a shear cell with sizes 0.027 × 0.025 m in the y × z directions, and has a length 0.2m in x direction where periodic boundary conditions are applied.The lateral sides in the x − z plane and the lower boundary in the x − y plane are smooth and frictional walls, with the same mechanical and frictional properties as the grains (Table S2).The initial concentration of large grains is uniform in the simulation domain.The grains are free to move in other directions (e.g. to dilate) in order to resemble a free-surface and shear-driven system.The grains are modeled as compressible spheres of diameter d s,l that interact when in contact via the Hertz-Mindlin model [54][55][56]: where the first term is total normal force, F n , and the second term is total tangential force, F t .In Equation 7, k n and k t are normal and tangential stiffness respectively, δ is the overlap between grains, γ n is the normal damping, v is the relative grain velocity, n ij is the normal vector at grain contact, t ij is the tangential vector at grain contact, γ t is the tangential damping.The full model implementation is available on the LAMMPS/LIGGGHTS webpage and several references [57][58][59].In accord with the low Stokes number of our laboratory experiments, the restitution coefficient is chosen to very small (e n = 0.01) such that collisions are highly damped (Table S2).Otherwise, there is no treatment of the viscous fluid in DEM simulations.The DEM model system is frictional, meaning that the coefficient of friction, µ, is the upper limit of the tangential force through the Coulomb criterion F t = µF n .
The tangential force between two grains grows according to non-linear Hertz-Mindlin contact law until F t /F n = µ and is then held at F t = µF n until the grains lose contact.The values of density, grain diameter, Poisson's ratio and acceleration due to gravity are chosen to match the experimental conditions.The values for coefficient of restitution and friction coefficient are chosen to mimic the effects from interactions with the fluid.The Young's modulus of the particles used here is chosen to be low (Table S2), MPa rather than GPa, in order to increase the calculation time step and decrease computational cost; however, since the system is not under significant confining pressure, a softer grain-grain interaction will not have considerable effect on the results, and the simulation remains in the hard-sphere limit.The particles are sufficiently hard that we find no additional rescaling of time is necessary.The damping coefficients γ n and γ t are determined within the implementation of LIGGGHTS from the chosen value for the restitution coefficient, e n .Movie S2

FIG. 1 .
FIG. 1. Phenomenology and setup.(A) Bed sediment of the River Wharfe, U.K., that shows a pronounced surface armor.Photo courtesy D. Powell [23].(B) Sketch of the experiment, showing position of the camera and laser plane used for imaging inside the granular bed.(C-E) Snapshots during armor development for τ * s = 3.8τ * cs .Also shown is the fluid boundary stress, which is

FIG. 2 .
FIG. 2. Experimental particle and segregation dynamics.(A) Vertical velocity profile for small and large grains for the interval ∆t = [0 : 20] minutes at the beginning of shearing at Shields number τ * s = 4.1τ * cs .The elevations of the bed surface (z s ) and transition to creep (z c ) are indicated.Inset shows the horizontal streamwise velocity u x profile for all grains and vertical velocity for large grains u z,l in logarithmic scale for the bedload zone.The streamwise particle velocity measurements have error-bar (one standard deviation/mean) value ∼ 0.3%, whereas vertical particle velocity measurements have error-bar (one standard deviation/mean) value ∼ 1%.(B) Vertical meansquare displacement (MSD) for large grains as a function of time ∆t.MSD shown for various depths of the granular bed, defined with a colorbar.The top and bottom red-dashed lines indicate the limiting behaviors of advection and diffusion, respectively.Boundary between bed load and creep is indicated.Note that near-surface grains are advective at short times; creeping grains show no change in MSD at short times indicating caged dynamics, but transition to diffusive behavior at longer times.(C) 1D (x-averaged) concentration map of large grains over time for shear stress τ * s = 4.1τ * cs .The red dashed lines show the positions of the armor surface (z s a) and the bottom of the armored layer (z i ), which are used to calculate the thickness of the armored layer in Figure 3 A. (D) Same as (C) for the advection-diffusion model, using velocity profiles and initial conditions that correspond to the shear stress τ * s = 4.1τ * cs in (A).

FIG. 3 .
FIG. 3. Armor thickness and segregation flux through time.(A) Temporal evolution of the thickness of the armored layer at different Shields number.Legend indicates the Shields number associated with each curve, and applies to (B) and (C) also.The brighter continuous lines are predictions from the advection-diffusion model Eq. 2. Note the first rapid stage of armoring which is dependent on Shields number and is associated with bed-load transport, and the second slower stage that exhibits a nearly constant rate for all Shields numbers and is the result of creep.(B) The variation of segregation flux density (J/A) with time.(C) Normalized flux against dimensionless time; data are reasonably collapsed.
].The maximum bulk segregation velocity, q, for each experiment is measured from relative displacement of the vertical z component of the center of mass position of large particles Z CM,l relative to small particles Z CM,s .The data for the relativeZ CM displacementfor all five stresses is presented in Figure S3.The initial concentration profile of large particles φ l (0) is determined from the first 10 s of each experimental run; a simplified version of it is used as another input to the PDE model (Figs.S4 to S8).The time-averaged streamwise velocity profiles for each experiment are reported in Figs.S4 to S8 and are used to estimate the value of β.We use a numerical implementation of the method of lines solution to solve the PDE equations and use N ts = 10000 time steps for the full experimental time (∼ 10 5 s).Comparisons of the concentration maps for the PDE model and experiments, for the four additional driving stresses, are presented in the supplementary materials (Figs.S4 to S8).

FIG. 4 .
FIG. 4. DEM simulation of a dry, sheared granular bed with u top = 0.08 m/s equivalent to the fluid-driven sheared granular bed at τ * s = 4.1τ * cs .The top layer of large grains that drives particles underneath is shown in Fig S9. (A) Model domain and initial conditions.Granular pack shown after (B) 5 s and (C) 50 s of shearing.Note rapid segregation, and depletion of the near-surface zone of large grains, as a consequence of bed load.(D) Evolution of armor thickness for the DEM model and experiments at two equivalent driving stresses indicated in the legend.
FIG. S1. (A)A schematic for the preparation protocol used for preparing the initial bidisperse mixture of particles.Small and large particles are initially deposited in an inversely segregated manner (panel B).The initial deposition has been performed by gently pouring first large particles and then small particles uniformly but manually from a very close distance to the bottom of the chamber.The system is then subjected to a rotation of Ω = 22 r.p.m. that is a driving stress equivalent to τ * s = 20τ * cs for ∼ 1 minute in order to mix the large and small populations (panel C).This shear stress was sufficient to fully suspend all particles in the channel.Fluid shear was then stopped completely and the suspension was left for ∼ 30 minutes to allow time for sedimentation, and relaxation and compaction of the granular bed, to reach completion.The final random packed layer at the end of preparation protocol has a thickness ∼ 15.5d s for all experiments.This state is used for the five Shields number armoring experiments reported in the main manuscript.(C) A hypothetical fully segregated state that one could obtain if continuing shearing the initial sample at the large preparation shear stress of τ * s = 20τ * cs for about 3 minutes.
FIG. S9.A snapshot from the initial conditions of the numerical DEM simulation that shows the layer of large grains deposited at the surface and moving at constant velocity u top .

TABLE S1 .
Parameters of the PDE simulations for five Shields numbers studied in the main manuscript.

TABLE S2
The top side in the x − y plane is open.The cell is filled with N =38812 grains that are initially inserted randomly in the cell with a desired volume fraction of 0.45.It is then equilibrated under gravitational forces for 10 million time steps equivalent to ∆t = 20 s.