2.1. Simulation methods

The segregation behaviours of size-bidisperse mixtures on an inclined plane are studied using DEM simulations using the open-source software LAMMPS (Thompson et al. Reference Thompson2022). Grains are simulated as spheres of intrinsic density $\rho _{*}$ and diameters $d^s$ and $d^l$ for the small and large grains, respectively. The size ratio $R = d^l / d^s$ is in the range $2 \leq R \leq 7$ in this paper, with particular focus on the cases $R=3$ and $R=6$ . To prevent crystallisation, each species of grain is given a 10 % polydispersity either side of their mean diameter $d^{\nu }$ , with diameters uniformly distributed in the range $(0.9d^{\nu },1.1d^{\nu })$ . The particles are subject to gravitational acceleration $g$ .

Contacts between particles are resolved using a linear spring–dashpot model with static friction (Cundall & Strack Reference Cundall and Strack1979). For the normal force, the coefficient of restitution is $\epsilon =0.8$ and a fixed binary collision time of $t_c=0.005 (d^s/g)^{1/2}$ is used as in Tunuguntla, Weinhart & Thornton (Reference Tunuguntla, Weinhart and Thornton2017). The normal spring stiffness $k_n$ and normal viscoelastic damping constant $\gamma _n$ can then be obtained from $t_c$ and $\epsilon$ (e.g. Thornton et al. Reference Thornton, Weinhart, Luding and Bokhove2012). For the tangential force, the coefficient of friction is $\mu _c=0.5$ , the tangential stiffness is $k_t=2/7k_n$ and the tangential viscoelastic damping coefficient is $\gamma _t=0.5\gamma _n$ . The simulation time step is $t_c/50$ to ensure numerical accuracy (Silbert et al. Reference Silbert, Ertas, Grest, Halsey, Levine and Plimpton2001).

The mixture flows down a rough plane inclined at angle $\theta$ to the horizontal. A Cartesian coordinate system is chosen with $(x,y,z)$ denoting the down-slope, cross-slope and slope-normal directions, respectively. Periodic boundary conditions are enforced in the $x$ and $y$ directions. The results are independent of numerical domain size, verified by varying the size of the periodic box from 12 $d^l$ to 36 $d^l$ in $x$ and 8 $d^l$ to 10 $d^l$ in $y$ . This is 72 $d^s$ to $216d^s$ and 48 $d^s$ to $60d^s$ at size ratio $R=6$ . To minimise slip, the base is made rough by fixing a layer containing a mixture of small and large grains with the same frictional and material properties as those in the flow. This minimises slip at the base whilst also preventing the small particles from filling the base and thus being lost to the flow.

The small-particle fraction $\overline {\phi ^s}$ is defined as the total volume of small particles divided by the total volume of granular material, and is varied between simulations. Since the intrinsic density of both species is the same, $\overline {\phi ^s}$ is equivalent to the mass fraction of small particles. When $\overline {\phi ^s}=1$ the flow consists of entirely small particles, and when $\overline {\phi ^s}=0$ it corresponds to entirely large particles. Small-particle fractions in the range $0.1 \leq \overline {\phi ^s} \leq 0.9$ are considered.

The initial condition consists a layer of small grains lying on top of a layer of large grains, with all grains initially at rest. The simulations are run until a statistically steady state is reached and the resulting steady-state profiles do not depend on the initial arrangement. This independence has been verified by comparing with simulations with the same small-particle fraction $\overline {\phi ^s}$ in which the large particles are initially lying on top of the small particles. The total mass of flowing grains per unit planar area is fixed across all simulations to $60\rho _{*}d^s$ , resulting in a flow height $h \approx 100 d^s$ . A complete list of all DEM simulations on an inclined plane performed in this paper is provided in Appendix A.

Figure 1. Snapshots of 18 bidisperse DEM simulations in steady state on a plane inclined at $25^{\circ }$ . The small (white) particles are rendered partially transparent to improve the visibility of the large (red) particles. Shown are flows at size ratios (a) $R=3$ and (b) $R=6$ . In both rows, the small-particle fraction $\overline {\phi ^s}$ increases from $\overline {\phi ^s}=0.1$ to $\overline {\phi ^s}=0.9$ in each plot from left to right. At $R=3$ the large particles always rise and accumulate at the surface, irrespective of $\overline {\phi ^s}$ . The behaviour at $R=6$ is qualitatively different, with either large or small particles accumulating at the surface, depending on $\overline {\phi ^s}$ .

2.2. Segregation behaviour

Mixtures of size ratio $R=3$ segregate in a conventional way with large particles rising to the surface (figure 1 a). However mixtures with $R=6$ exhibit a new and qualitatively different segregation behaviour (figure 1 b) in which large particles may either rise or sink, or where segregation is strongly suppressed and the grains are almost uniformly mixed throughout the flow depth. The large particles may also form distinct layers which are separated by thin bands of small particles.

At $R=6$ , whether large or small particles accumulate at the surface depends on the small-particle fraction $\overline {\phi ^s}$ . For $\overline {\phi ^s} \gtrsim 0.7$ , the large particles are concentrated near the base of the flow (figure 1 b), while their concentration decreases towards the surface. A surface layer of purely small particles forms, rendering the large particles invisible when viewed from above. For lower $\overline {\phi ^s}\lesssim 0.5$ the segregation behaviour changes as the large particles accumulate at the surface, forming a pure layer devoid of small particles. The transition between these regimes occurs at $\overline {\phi ^s}\approx 0.6$ , where the distribution of particles is almost uniform with depth.

This behaviour differs markedly from the well-described behaviour at $R=3$ where the large particles segregate to the surface at all values of $\overline {\phi ^s}$ , while their concentration decreases with depth (figure 1 a). The rise of large particles (of equal intrinsic density) is typical of segregation in bidisperse mixtures of size ratio three or less (Savage & Lun Reference Savage and Lun1988; Golick & Daniels Reference Golick and Daniels2009; Wiederseiner et al. Reference Wiederseiner, Andreini, Epely-Chauvin, Moser, Monnereau, Gray and Ancey2011; van der Vaart et al. Reference van der Vaart, Gajjar, Epely-Chauvin, Andreini, Gray and Ancey2015). In contrast to the segregation behaviour of grains of only slightly varying size, at size ratio $R=3$ the uppermost part of the flow is composed solely of large particles, with the small particles totally absent. The thickness of this layer decreases as the fraction of small particles $\overline {\phi ^s}$ increases.

The phenomenon of large particles sinking when present in low concentrations and settling close to the base of the flow was termed reverse segregation by Thomas (Reference Thomas2000), who observed this behaviour in deep chute flows of small particles with a small proportion of large grains added ( $\overline {\phi ^s} \geq 0.9$ ). The simulations performed in this paper at $R=6$ show that as the fraction of large particles in the flow is increased (i.e. as $\overline {\phi ^s}$ is decreased), the large particles progressively fill upwards from the base (figure 1 b). At $\overline {\phi ^s}=0.6$ the large particles reach the surface and there is an approximately uniform distribution throughout the depth. For $\overline {\phi ^s}\lt 0.6$ , there are sufficient large particles that they form a pure layer on the surface, riding on top of a mixed region of small and large particles in a configuration that is more similar to segregation at $R=3$ .

Notably, the small-on-top case at $R=6$ ( $\overline {\phi ^s} \gtrsim 0.7$ ) is qualitatively different from simply flipping the corresponding large-on-top case at $R=3$ upside down. In particular, at $R=6$ the large particles do not form a pure basal layer, but remain intermixed with small particles near the base (figure 1 b). Furthermore, in these mixed small- and large-particle regions, large particles are organised into discrete horizontal layers, with each layer being one large particle thick. These layers are separated by thin regions of small particles, and adjacent large particles do not appear to make direct contact. This layering of large particles is discussed in detail in § 3.

The segregation patterns in the images of figure 1 can be represented in continuum concentration profiles, where the concentration $\phi ^{\nu }$ is defined as the local volume fraction of species $\nu$ per unit granular volume, where $\nu =s,l$ for small and large particles, respectively. Continuum concentration fields are constructed using the coarse-graining method of Tunuguntla et al. (Reference Tunuguntla, Thornton and Weinhart2016) with a Gaussian smoothing function as described in Appendix B. The coarse-graining width is set to $c=0.8d^l$ , a value chosen from a plateau on which continuum fields are insensitive to the coarse-graining width (Tunuguntla et al. Reference Tunuguntla, Thornton and Weinhart2016). Choosing this coarse-graining width comparable to the large-particle diameter suppresses oscillations that would occur at smaller coarse-graining widths, due to the particle-scale layering observed in figure 1. This allows the underlying variation in bulk composition to be studied.

Figure 2 shows such continuum volume fraction profiles $\phi ^s(z)$ for the simulations illustrated in figure 1. At $R=3$ (figure 2 a) the small-particle concentration $\phi ^s$ decreases with height over all small-particle fractions, indicating segregation in which the large particles accumulate at the surface. By contrast, at $R=6$ (figure 2 b) there is a distinct switch in behaviour of the profiles when the local concentration $\phi ^s$ reaches a specific ‘neutral’ value $\phi ^s_n \approx 0.6$ (dashed line in figure 2 b), which separates the large-on-top and small-on-top configurations. For $\phi ^s\lt \phi ^s_n$ , the profiles decrease with height resulting in a large-particle layer at the surface. For $\phi ^s\gt \phi ^s_n$ , the concentration profiles increase with height resulting in a small-particle surface layer. Near the neutral concentration $\phi ^s\approx \phi ^s_n$ , the mixture is at an almost uniform concentration with depth. The local concentration $\phi ^s$ therefore plays a crucial role in determining whether large particles sink towards the base, rise to the surface or remain at the same depth.

Figure 2. Profiles of small-particle concentration $\phi ^s$ as a function of slope-normal coordinate $z$ for (a) $R=3$ and (b) $R=6$ . Coloured lines represent different small-particle fractions, $\overline {\phi ^s}$ , increasing from 0.1 (red) to 0.9 (blue) in increments of 0.1. Selected values of $\overline {\phi ^s}$ are indicated on the plots. In (b), the dashed line at $\phi ^s = 0.62$ separates simulations with a surface composed of small particles and those with a large-particle surface.

Figure 3. Snapshots of bidisperse DEM simulations in steady state on a plane inclined at $25^{\circ }$ , with the flow going from left to right within each panel. The small (white) particles are partially transparent to improve the visibility of the large (red) particles. The small-particle fraction is fixed at $\overline {\phi ^s}=0.7$ , and the size ratio increases from $R=2$ to $R=7$ as indicated beneath the images.

The influence of the size ratio $R$ on segregation is illustrated in figure 3. At fixed $\overline {\phi ^s}=0.7$ , increasing $R$ from 2 progressively reduces the thickness of the surface layer of large particles. At $R=5$ this layer is only one large particle thick and the mixture becomes nearly uniform with depth, similar to the behaviour observed for $R=6$ at $\overline {\phi ^s}=0.6$ . At $R=5.4$ , there is no longer a distinct surface layer of large particles, with both small and large particles lying at the surface and a near-homogeneous mixture throughout the depth. As $R$ is increased further (to 6 and 7) the large particles are concentrated towards the flow base, and are entirely absent from the surface.

We hypothesise that the transition from a small-particle surface to a large-particle surface, as the fraction of large particles increases, arises from a competition between two mechanisms: the effective density effect discussed in the introduction and the migration of small particles through voids. At low concentrations of large particles, their higher density relative to the surrounding medium of small particles and void space causes them to sink. This displaces the smaller particles upward, and results in a surface layer composed of entirely small particles. As the concentration of large particles increases, the large particles are no longer isolated but are in close proximity to other large particles, and the average density of the granular mixture surrounding a large particle rises, decreasing the density difference. Once the large particles are densely packed, the mechanisms of kinetic sieving and squeeze expulsion may take place, allowing small particles to fall downwards through the voids between large particles, displacing them upward. At sufficiently high fractions of large particles, this process leads to a layer of large particles at the surface above a mixed region of large and small particles.

2.3. Neutral concentration

Whether the free surface of an inclined-plane flow is made up of small or large particles (that is, whether large particles have a tendency to sink or rise, respectively) depends on both the size ratio $R$ and the small-particle fraction $\overline {\phi ^s}$ . At $R=6$ the transition occurs close to $\overline {\phi ^s}\approx 0.6$ (figure 2 b), while at $R=5.4$ it shifts to $\overline {\phi ^s} \approx 0.7$ (figure 3). Close to the transition, the concentration profiles are approximately uniform with depth. Simulations can therefore be grouped into those with a large-particle surface and those with a small-particle surface, as shown in figure 4. A simulation is defined to have a large-particle surface if its $\phi ^s(z)$ profile moves towards $\phi ^s=0$ as $z \to h$ , as seen for all cases with $\overline {\phi ^s} \leq 0.6$ in figure 2(b), for example. At a given $R$ , the transitional value of $\overline {\phi ^s}$ corresponds to the value of the local concentration $\phi ^s$ at which the profiles become approximately uniform with depth. This is the neutral concentration $\phi ^s_n$ . Therefore the black line sketched in figure 4, which approximately separates flows with a large- and small-particle surface, is the curve $\overline {\phi ^s} = \phi ^s_n(R)$ .

Figure 4. Size ratio $R$ and small-particle fraction $\overline {\phi ^s}$ for which mixtures have a small- or a large-particle surface at $\theta = 25^{\circ }$ . Blue squares indicate a large-particle surface and orange crosses indicate a small-particle surface. The sketched black curve divides the two regimes and represents the neutral concentration $\phi ^s_n(R)$ .

These results connect to the more well-studied limit of a single large intruder. Our DEM simulations on an inclined plane with a single large intruder within a flow of small grains show that the intruder remains at or just below the surface for $R \leq 3.8$ and sinks towards the base for $R \geq 4$ ; this is consistent with previous studies (Thomas & D’Ortona Reference Thomas and D’Ortona2018; Jing et al. Reference Jing, Ottino, Lueptow and Umbanhowar2020). In figure 4, these single-intruder simulations appear with $\overline {\phi ^s}=1$ , since the flow is made up of entirely small particles apart from a single large grain. Therefore the value of $R$ for which $\phi ^s_n \approx 1$ is approximately $R \approx 4$ .

3.1. Inclined-plane simulations

In the inclined-plane simulations, the initial configuration consists of smaller particles above the larger ones. Once flow commences, the small particles begin to filter down towards the base, but no layering is apparent after a short time, as shown in figure 5(a) for $R=6$ and $\overline {\phi ^s}=0.4$ at $t=120 (d^s/g)^{1/2}$ . With continued flow, however, the system transitions into a layered structure as shown in figure 5(b). A video of the transition into a layered state is available in supplementary movie 1 available at https://doi.org/10.1017/jfm.2026.

The overall layered and segregated structure is recovered after $t \approx 600(d^s/g)^{1/2}$ as shown in figure 5(c), in which the small-particle concentration $\phi ^s$ is plotted as a function of $z$ and $t$ . For context, with typical small grain size $d^s=100$ μm, a flow depth $h=1$ cm and gravitational acceleration $g=9.8$ $\mathrm{m\,s^{-2}}$ , this corresponds to a time $t\approx 1.9$ s. To resolve the layers, the simulations are coarse-grained with a width of $c=0.2d^s$ , as indicated by the notation $\phi ^s(t,z;c=0.2d^s)$ (Appendix B). This very small coarse-graining width effectively visualises the vertical distribution of large-particle centres, rather than the volume occupied by the large particles. The width of the red stripes is indicative of the variation in the $z$ position of large-particle centres within a layer, typically $\lesssim 0.4 d^l$ in the lower half of the flow in figure 5(b,c), and increasing towards the surface.

Successive layers form in the direction of the velocity gradient (i.e. along the $z$ axis), with each layer lying in the $x$ – $y$ plane. The lowest layer forms approximately three small particle diameters above the base, and the layered structure extends upward through the flow. The layers are one large particle thick, with large particles separated by thinner regions of small particles typically one to three small-particle diameters thick. The separation between layers, measured as the distance between large-particle centres in adjacent layers, is approximately constant and is equal to $d^s+d^l$ , although this distance increases slightly as the flow dilates (figure 5 c). This layered arrangement extends upward until either the number of large particles becomes insufficient to sustain additional layers or the flow transitions into a nearly monodisperse region of large particles near the surface. In the latter case, the lack of small particles prevents them from effectively separating the large-particle layers, causing the layered structure to break down. Thus, the layering here is present in bidisperse regions but disappears when the flow is approximately monodisperse, indicating a qualitatively distinct mechanism from the crystallisation that can occur in monodisperse flows (Silbert et al. Reference Silbert, Grest, Plimpton and Levine2002). Layers are less distinct in simulations with low small-particle fractions (e.g. $\overline {\phi ^s} = 0.1$ and $0.2$ in figure 1 b), again pointing to a lack of small particles inhibiting layer formation.

Layering may be identified from the small-particle concentration $\phi ^s$ profiles, computed using a fine coarse-graining width $c=0.2d^s$ . The $\phi ^s$ fields are averaged in the flow ( $x$ ) and cross-flow ( $y$ ) directions, so only the variation in the velocity-gradient ( $z$ ) direction, in which the layers form, remains. In these profiles, layers appear as sequence of alternating maxima and minima (figure 5 c). A minimum is classified as a large-particle layer if two criteria are satisfied: (i) the amplitude between the minimum and both its adjacent maxima exceeds 0.3 and (ii) the vertical separation between successive minima is $d^l \leq \Delta z \leq 1.5d^l$ . The small-particle concentration $\phi ^s$ associated with this layer is estimated by coarse-graining the profile at that height using a broader width $c=0.8d^l$ , which smooths over the oscillations. The coarse-graining procedure is described in Appendix B.

Layering does not occur across all size ratios $R$ or local small-particle concentrations $\phi ^s$ . Rather, applying the above definition of layering to all of the simulations in figure 2(b) indicates that, for $R=6$ , layering emerges within the range $0.28 \lesssim \phi ^s \lesssim 0.86$ . Below this range, the proportion of large particles is too high for distinct layers to form. At very high $\phi ^s$ there are too few large particles to distinguish layers.

Some layering is also present at $R=3$ , particularly for the simulations at $\overline {\phi ^s}=0.3$ and 0.4 (figure 1 a). However the layering appears over a narrower range of concentrations, approximately $ 0.32 \lesssim \phi ^s \lesssim 0.76$ at this size ratio. Furthermore, the layers at $R = 3$ are less distinct than at $R = 6$ , suggesting that layering becomes more pronounced and appears over a wider range of $\phi ^s$ as the size ratio increases. Supporting this trend, the set of simulations at $\overline {\phi ^s}=0.7$ in figure 3 display clear layering for $4 \leq R \leq 7$ but not for the smaller size ratios $R=2$ and 3.

The profiles of the down-slope velocity $u$ in figure 6(a) and corresponding small-particle concentration $\phi ^s$ profiles in figure 6(b–d), which are taken from figure 5(c) at the times indicated, reveal a feedback of layering on the flow rheology. After layering occurs, the velocity profiles acquire a stepped structure, with each step corresponding to one layer, while becoming increasingly linear across the entire layered region. By $t=1400(d^s/g)^{1/2}$ the shear rate begins to decrease only in the thin, approximately monodisperse, region of large particles near the flow surface. This near-linear profile may suggest the rheology in layered bidisperse mixtures deviates from the Bagnold scaling (Bagnold Reference Bagnold1954) and, consequently, from constitutive relations based on the $\mu (I)$ rheology (Jop, Forterre & Pouliquen Reference Jop, Forterre and Pouliquen2005; Rognon et al. Reference Rognon, Roux, Naaim and Chevoir2007; Tripathi & Khakhar Reference Tripathi and Khakhar2011).

Figure 6. (a) The down-slope velocity profiles $u(z)$ plotted against the slope-normal coordinate $z$ at times before and after layering occurs, for the bidisperse inclined-plane flow at $R=6$ in figure 5. The times are given in units of $(d^s/g)^{1/2}$ . The corresponding small-particle concentration $\phi ^s$ profiles at times (b) $t=400$ , (c) $t=900$ and (d) $t=1400$ . Note that the vertical scale is $z/d^l$ . All profiles are from the simulation in figure 5 and computed with a coarse-graining width $c=0.2d^s$ .

As layering develops, the flow near the free surface accelerates, an effect that is visible in figure 6(a) and in supplementary movie 1. In contrast, the flow near the base becomes slower. The flow also dilates after layering occurs as figure 6(c,d) shows, with individual layers rising upwards, even though most segregation is already complete. These behaviours are consistent with a feedback between layering and rheology. Figure 6 suggests that the relationship between layering and rheology is rich and potentially complex. A detailed characterisation lies beyond the scope of this paper.

3.2. Lees–Edwards simulations

Figure 7. The time-dependent formation of layers in a Lees–Edwards simulation at $R=6$ . The system is periodic in all spatial directions and so is invariant to a shift in the $z$ axis. (a) The small-particle concentration $\phi ^s$ as a function of velocity-gradient coordinate $z$ and time $t$ . To resolve the layered structure, a coarse-graining width $c=0.2d^s$ is used. The transition into a layered structure is shown in supplementary movie 2. (b) The corresponding velocity profile at $\dot {\gamma }t=132$ .

Figure 8. Layered structure in a Lees–Edwards simulation at $R=6$ . (a) A snapshot of the layering, obtained by averaging 50 frames over a time interval $25/\dot {\gamma }$ . Flow is in the positive $x$ direction. (b) The small-particle concentration as a function of the velocity-gradient coordinate $z$ and cross-flow coordinate $y$ . A coarse-graining width $c=0.2d^s$ is used to resolve the layered structure. (c) The structure in the $x$ – $y$ plane, of the lowest layer in (a,b).

In the inclined-plane simulations, grains are subject to a gravitational body force, which drives segregation. Furthermore there is a rigid rough boundary at $z=0$ . To test whether either gravity or the presence of a boundary is necessary for layering, DEM simulations are performed in LAMMPS (Thompson et al. Reference Thompson2022) with Lees–Edwards boundary conditions (Lees & Edwards Reference Lees and Edwards1972). These boundary conditions impose a simple shear velocity profile $\boldsymbol{u}=(\dot {\gamma }z,0,0)$ with shear rate $\dot {\gamma }$ and are implemented in LAMMPS via the fix deform command which shears the simulation domain and is equivalent to the standard sliding brick representation (Evans & Morriss Reference Evans and Morriss2007). These simulations are performed in a fully periodic, three-dimensional domain with no rigid boundaries. Gravity is absent and so gravity-driven segregation does not occur.

A bidisperse mixture is simulated at a size ratio $R=6$ , a small-particle concentration $\phi ^s=0.5$ and total solids volume fraction $\varPhi =0.63$ , a value representative of the simulations on an inclined plane at $R=6$ . A polydispersity of $10 \,\%$ either side of the mean diameter $d^s$ is retained for the small particles, whilst the large particles are strictly monodisperse with diameter $d^l$ . In these simulations the large particles are given a single diameter to more clearly study the separation between them, though introducing $\pm 10 \,\%$ polydispersity among the large particles does not qualitatively change the behaviour seen. A fixed binary collision time of $t_c=0.00125/ \dot {\gamma }$ is used, ensuring the particles are in the asymptotic hard-particle regime where the rheology is independent of the particle stiffness (Chialvo, Sun & Sundaresan Reference Chialvo, Sun and Sundaresan2013). Other particle interaction parameters are the same as those described in § 2.1. Similarly to the inclined-plane simulations, a Cartesian coordinate system is chosen with $x$ representing the flow direction, $y$ the vorticity direction and $z$ the velocity-gradient direction. The box lengths are $L_x=48d^l$ , $L_y=10d^l$ and $L_z=12d^l$ in the $x$ , $y$ and $z$ directions, respectively. These lengths are varied to ensure the conclusions are independent of domain size. Particles are initially placed at random positions within the domain.

Despite starting in a random, unlayered configuration, once a shear velocity profile has been imposed, the grains transition into a layered state as shown in figure 7, which plots the small-particle concentration as a function of time $t$ and velocity-gradient coordinate $z$ , with a coarse-graining width $c=0.2d^s$ . The layered structure is qualitatively similar to that in inclined-plane simulations (figure 5 c), with the layers emerging from a disordered state and forming along the velocity-gradient direction. Successive layers are again separated by regions of small particles. The evolution into a layered state is shown in supplementary movie 2. Lees–Edwards boundary conditions impose a linear velocity profile, which remains even after layering has occurred; however, the velocity profile becomes stepped around the overall linear variation (figure 7 b). Each step corresponds to a single large-particle layer, consistent with figure 6(a).

By $\dot {\gamma }t = 60$ , the system has evolved into a clearly layered state (figure 7). All Lees–Edwards simulations performed exhibited layering by this time, but the time scale depended somewhat on the streamwise box length $L_x$ . Larger boxes generally layered more slowly, but this dependence plateaued for $L_x \geq 48 d^l$ . For reference, two large particles separated by a distance $d^l$ in the $z$ direction would take a non-dimensional time $\dot {\gamma } t = L_x/d^l$ to shear apart from one another by the streamwise box size $L_x$ and approach each other again due to the periodic boundaries. In the largest domains studied ( $L_x=200d^l$ ), this time is $\dot {\gamma } t = 200$ , well above the observed layering time of $\dot {\gamma } t \approx 60$ . Hence, layering occurs on a time scale shorter than that required for such periodic boundary effects. The layering time showed little sensitivity to the box dimensions in the $y$ and $z$ directions.

Within each layer of large particles, large particles are arranged into straight lines parallel to the flow direction (figure 8). These chains of large particles can be seen in the image in figure 8(a), which is obtained by averaging 50 frames over a time interval of $25/\dot {\gamma }$ once the simulation has layered. Chains are also present in the inclined-plane simulations as shown in figure 5(b) and supplementary movie 1. Each large particle tends to follow directly behind the one ahead, creating extended chains that can span the entire periodic domain, even for the longest streamwise box length of $200d^l$ tested in this study.

Figure 9. The RDF $g(r)$ for large-particle pairs at $R=6$ , as defined in (3.1), with bin width $\Delta r=0.01d^s$ . The blue line represents an unlayered state measured at time $\dot {\gamma }t=22$ and the red line a layered state measured at $\dot {\gamma }t=132$ . The shaded region indicates $r \lt d^l=6d^s$ , where large-particle centres cannot lie due to volume exclusion. The dashed vertical line marks the distance $r = d^l + d^s= 7d^s$ , corresponding to large particles separated by one mean small-particle diameter.

The chains of large particles are further illustrated in figure 8(b) which plots the small-particle concentration $\phi ^s$ as a function of cross-flow coordinate $y$ and velocity-gradient coordinate $z$ , averaged over the box length $48d^l$ in the streamwise direction. The small-particle concentration $\phi ^s$ from the simulation in figure 8(a) is used, with a coarse-graining width of $c=0.2d^s$ . Each red dot therefore represents a large-particle chain. The chains are regularly spaced in the $y$ direction and, combined with the layering in the $z$ direction, they produce the hexagonal pattern visible in figure 8(b). With the box size chosen this arrangement does not form a perfectly regular lattice. The box, of length $12d^l$ in the $z$ axis, results in 11 layers, making the periodicity inconsistent with a regular lattice. Despite this constraint, a lattice-like structure still emerges. A similar, though weaker, hexagonal arrangement is also observed in the $x$ – $y$ structure of a single layer, shown in figure 8(c), which corresponds to the lowest layer (closest to $z=0$ ) in figure 8(a,b).

The self-organisation of large particles into a regular structure (figure 8), with layers forming in the $z$ direction and chains extending in the $x$ direction, indicates that the large particles arrange themselves in a way that minimises collisions between neighbouring large particles. The avoidance of direct contact between large particles is evidenced by the radial distribution function (RDF), $g(r)$ , for large-particle pairs (figure 9). The RDF provides a measure of the probability of finding a large-particle centre at distance $r$ from the centre of a reference large particle, relative to the probability expected for a uniform distribution at the same number density. Peaks and troughs in $g(r)$ indicate preferred separations or avoidance, while, in a homogeneous system, $g(r)$ approaches 1 in the asymptotic limit $r \to \infty$ , indicating particle positions are uncorrelated at large distances. Specifically, the RDF is defined as (Wiacek Reference Wiacek2016)

(3.1) \begin{equation} g(r)=\frac {N(r)}{n^l 4 \pi r^2 \Delta r}, \end{equation}

where $n^l$ is the number density of large particles and $N(r)$ is the number of large particles with centres in a spherical shell of thickness $\Delta r$ , located at radial distance $r$ from the centre of a reference particle. In this study $\Delta r=0.01d^s$ is used. Time-averaging is performed on $g(r)$ using 29 evenly spaced snapshots over a time interval of $14/\dot {\gamma }.$

There is a striking difference between the RDFs before and after layering occurs (figure 9). At $R=6$ , in an unlayered state, $g(r)$ exhibits a peak around $d^l+d^s=7d^s$ indicating that, with these simulation parameters, large particles preferentially position themselves around one mean small-particle diameter apart. There is a second, sharper peak at $r = 6d^s = d^l$ , representing large particles in direct contact. This peak is likely sharper because the large particles here have no polydispersity: when they are in direct contact, their centres are separated by exactly $d^l$ . Values $r\lt d^l$ are not permissible due to volume exclusion, apart from some small overlaps allowed by the soft-sphere DEM. In contrast, once layering has occurred, the peak at $r = d^l$ vanishes, and $g(r = d^l) \approx 0$ , indicating that large particles no longer make direct contact. In fact, $g(r)\lt 0.1$ for $r\lt 6.7d^s$ showing that large particles rarely approach each other at separations much smaller than one small-particle diameter. Instead, they are kept apart by thin regions of small particles, evidenced by the peak in simulation at $r \approx d^l+d^s$ showing that the large particles are positioned approximately one small particle away from their nearest neighbour in this simulation. Qualitatively similar behaviour is observed for mixtures with different $\phi ^s$ once layering develops.

Figure 10. The RDFs $g(r)$ for large-particle pairs from Lees–Edwards simulations at size ratio $R=6$ , small-particle concentration $\phi ^s = 0.5$ and solids volume fraction $\varPhi$ of (a) 0.54, (b) 0.57, (c) 0.6 and (d) 0.63. The dashed vertical line marks the distance $r = d^l + d^s= 7d^s$ , corresponding to large particles separated by one mean small-particle diameter. (e–h) The corresponding vertical $\phi ^s$ profiles, coarse-grained at $c=0.2d^s$ , directly beneath their associated RDFs. The system is periodic in all spatial directions and so is invariant to a shift in the $z$ axis.

The structure of bidisperse flows depends on the solids volume fraction $\varPhi$ , as evidenced by the differences in the RDFs for $0.54 \leq \varPhi \leq 0.63$ (figure 10 a–d). All RDFs are computed at size ratio $R=6$ and small-particle concentration $\phi ^s=0.5$ . At the lowest solids volume fraction of $\varPhi = 0.54$ , there is a pronounced peak at $r=d^l$ , indicating large particles in direct contact. As $\varPhi$ increases, this peak diminishes then disappears as the large particles are separated by thin bands of small particles. For reference, the jamming fraction at $R=6$ and $\phi ^s=0.5$ is approximately 0.7 (Srivastava et al. Reference Srivastava, Roberts, Clemmer, Silbert, Lechman and Grest2021), so all cases considered are well below jamming and are outside of the quasi-static regime.

The changes in the RDFs are reflected in the $\phi ^s(z)$ profiles when coarse-grained with $c=0.2d^s$ (figure 10 e–h). At $\varPhi =0.54$ , the flow is largely disordered and layering does not occur, allowing large particles to make direct contact, although there is some indication of a weak periodic structure. As $\varPhi$ increases to $0.57$ , 0.6 and 0.63, clear layering occurs, with large particles organised into increasingly well-defined bands. These observations indicate that layering not only depends on $R$ and $\phi ^s$ , but also on the solids volume fraction $\varPhi$ , with layering less likely in more dilute flows. Consequently, the phase space governing layering is at least three-dimensional, depending on $R$ , $\phi ^s$ and $\varPhi$ .

The emergence of the layered structure in a gravity-free, fully periodic shear cell suggests that large-particle layering is caused only by shear. In particular, it implies that neither gravity nor a rigid boundary are necessary for layering, although the presence of a rigid base may promote particle ordering (Weinhart et al. Reference Weinhart, Hartkemp, Thornton and Luding2013). Shear-driven ordering has recently been seen in simulations of dense, frictional size-bidisperse suspensions at a size ratio of three (Malbranche, Chakraborty & Morris Reference Malbranche, Chakraborty and Morris2023), where the large particles arranged themselves into distinct layers separated by thin regions of small particles. This behaviour is qualitatively similar to that seen in figure 8 and supplementary movie 2 and suggests that layering is a generic behaviour of size-bidisperse particle mixtures, rather than being due to the particular dry granular friction behaviour considered in this paper. The alignment of large particles into chains is also consistent with the experimental and numerical study of Fraysse et al. (Reference Fraysse, D’Ortona and Thomas2024), who showed that two large intruders in an inclined-plane flow of otherwise monodisperse small particles tend to align with the flow, over a range of flow depths and slope angles. The results in this paper indicate that this alignment is robust and extends to many large particles.

3.3. Experiments

Figure 11. Experimental set-up for the rectangular chute prior to the release of the grains. The chute is inclined at $70^{\circ }$ to the horizontal and is 120 cm long, with transparent, glass sidewalls spaced 18 mm apart. The top and bottom boundaries are roughened with 40 grit aluminium oxide paper and are separated by $53$ mm. At the lower end, grains are held in place by a gate that pivots at the right-hand boundary. When opened, the gate creates an opening at the bottom boundary, allowing particles to exit. The grains fall onto a mass balance that records mass as a function of time. In steady state, the measured mass flux is $0.07$ $\mathrm{kg\,s^{-1}}$ . The full experiment is shown in supplementary movie 3.

Figure 12. Photographs of the experiment in an inclined rectangular chute. (a) The initial disordered state prior to the onset of flow and (b) the system 11 s after flow has commenced, with large particles forming distinct layers. A sketch of the velocity profile is overlaid. Panel (c) reproduces (b), with eight equally spaced 1.9 mm wide rectangles applied from the top; the hue of alternate rectangles has been shifted to blue to highlight the layered structure that is present in (b). The full time evolution of the layering can be seen at full speed and 1/3 speed in supplementary movies 4 and 5, respectively. The top of the images is $27$ mm above the rigid base, placing it approximately halfway between the upper and lower boundary.

To experimentally validate the layer formation present in DEM simulations, an experimental configuration was chosen that subjects a deep region of size-bidisperse grains to sustained simple shear flow. The set-up consists of an enclosed rectangular chute, inclined steeply at $\theta =70^{\circ }$ to the horizontal, with a gate restricting the outflow flux (figure 11). This configuration creates a region of shearing flow at least 16 large particles ( $\approx 100$ small particles) deep. The experiment provides a means to test whether large-particle layering occurs outside the constraints of periodic simulation domains. To our knowledge, the formation of large-particle layers similar to those in figure 5(b) has not previously been observed in real granular flows.

The top and bottom boundaries of the experiment are 53 mm apart, and are roughened to induce shear using 40 grit aluminium oxide paper. The transparent front and back walls are separated by 18 mm and allow direct visualisation of the flow. The chute is 120 cm long and a camera films a portion of the flow, 16 cm from the outflow. A bidisperse mixture of dry spherical glass beads is used, consisting of small white (diameter 200–300 μm) and large red (diameter 1.40–1.65 mm) grains, with a mean size ratio $R \approx 6$ . The intrinsic density of the small and large beads is 1440 and 1510 $\mathrm{kg\,m^{-3}}$ , respectively. The chute is filled with $0.7$ kg of small grains and $0.5$ kg of large grains, which are initially approximately uniformly mixed. The formation of layers is robust and not sensitive to the precise initial conditions. Material exits the chute through a gate at the bottom which pivots around a point at the base of the right-hand boundary, as indicated in figure 11, to create a partial opening for the grains to flow through. The grains land on a balance that records the mass as a function of time. A video of the full experiment is available in supplementary movie 3.

The mass flux exiting the chute quickly stabilises at a steady value of $0.07\,\mathrm{kg\,s^{-1}}$ approximately 1 s after the gate is opened. Within the shearing region, the initially disordered state (figure 12 a) evolves into an ordered, layered structure (figure 12 b). The full transition into the layered state can be seen in supplementary movie 4, and a slow-motion version is available in supplementary movie 5. Layering emerges rapidly once the flow is underway, with visible evidence that large particles have formed layers apparent around 2 s after the gate has opened. The layering persists until the chute is fully emptied, at approximately 17 s.

The layering of grains occurs primarily where significant shear is present in the flow, as evidenced by the velocity profile sketched in figure 12(b). This profile, which is based on the tracked deformation of 13 large particles, illustrates the deformation of an initially straight line element, oriented normal to the slope. The shear rate is reduced near the base due to lateral confinement by the sidewalls (Jop, Forterre & Pouliquen Reference Jop, Forterre and Pouliquen2006; Lloyd et al. Reference Lloyd, Maguire, Mistry, Reynolds, Johnson and Gray2025). Layering is less evident in this region of low shear near the flow base, where the large-particle volume fraction also appears much smaller.

At least eight distinct layers of large particles are evident in the central shearing part of the flow, aligned with the velocity gradient, consistent with the layering observed in DEM simulations. Figure 12(c) reproduces figure 12(b) but the hue of the image has been artificially shifted to blue in four rectangles, each 1.9 mm in height and spaced by 1.9 mm (approximately the diameter of one large particle plus the diameter of between one and two small particles). In the resulting image almost all large grains lie either wholly within a red stripe or wholly within a blue stripe, a pattern which emerges because the large particles are arranged into discrete layers. In particular, it shows that the layers remain contiguous over a considerable distance in the flow direction, in figure 12(c) spanning the entire width of the frame of around 28 large particles. As shown in supplementary movies 4 and 5, the number of large particles visible within each layer can exceed 100, indicating that the layering extends well beyond the limited view captured in the static image.

The experiment confirms that layering is not simply a consequence of periodic boundary conditions or finite simulation domains, but is an experimentally realisable phenomenon that is intrinsic to real size-bidisperse granular flows. Layering was consistently observed across multiple repeated experimental runs giving confidence that the effect is robust. Additionally, the alignment of large particles into long chains and the organisation into discrete layers may both have important consequences for the rheology of such flows; however, this paper focuses on the description and modelling of the segregation behaviour.

4.1. Bidisperse mixture framework

The aim is now to model the segregation behaviour in figure 1 within a continuum framework. Consider a granular material composed of large and small particles of the same intrinsic density $\rho _{*}$ . Each species has diameter $d^{\nu }$ , where $\nu =l,s$ denotes large and small particles, respectively. The solids volume fraction of each species is $\varPhi ^{\nu }$ and the total solids volume fraction of the granular material is $\varPhi =\varPhi ^s+\varPhi ^l$ . In dense monodisperse flows (GDR-MiDi 2004) as well as bidisperse flows up to a size ratio $R=2.5$ (Tripathi & Khakhar Reference Tripathi and Khakhar2011), $\varPhi$ is nearly constant across the flow. Although a bidisperse mixture can pack more densely than monodisperse grains, and these packing effects may become more significant at larger size ratios (Srivastava et al. Reference Srivastava, Roberts, Clemmer, Silbert, Lechman and Grest2021), for simplicity and compatibility with existing segregation theories (e.g. Gray Reference Gray2018; Umbanhower et al. Reference Umbanhower, Lueptow and Ottino2019) $\varPhi$ is here approximated as uniform throughout the flow. Figure 1(b) shows some evidence of volume fraction variation with changes in the free-surface height across different values of $\overline {\phi ^s}$ . The observed variation is modest, however, at only $5\,\%$ – $10\,\%$ , and is therefore neglected. As a result, the bulk velocity $\boldsymbol{u}$ can be approximated as incompressible:

(4.1) \begin{equation} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}=0. \end{equation}

The theory is expressed in terms of the species’ concentration $\phi ^{\nu } \in [0,1]$ , defined, as in § 2, as the volume fraction per unit granular volume,

(4.2) \begin{equation} \phi ^{\nu }=\frac {\varPhi ^{\nu }}{\varPhi }, \end{equation}

such that concentrations sum to unity:

(4.3) \begin{equation} \phi ^s+\phi ^l=1. \end{equation}

Segregation is modelled using the well-established framework of advection–segregation–diffusion equations (see e.g. Gray Reference Gray2018), which governs the spatial and temporal evolution of the concentrations $\phi ^{\nu }$ through

(4.4) \begin{equation} \frac {\partial \phi ^{\nu }}{\partial t}+\underbrace { \boldsymbol{\nabla } \boldsymbol{\cdot }(\phi ^{\nu } \boldsymbol{u})}_{\text{advection}} + \underbrace { \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{F}^{\nu }}_{\text{segregation}}=\underbrace { \boldsymbol{\nabla } \boldsymbol{\cdot }\big (D\boldsymbol{\nabla }\phi ^{\nu }\big )}_{\text{diffusion}}. \end{equation}

In this framework, the flux of each species is assumed to be composed of three distinct components: an advective flux $\phi ^{\nu } \boldsymbol{u}$ that simply sweeps species $\nu$ along with the bulk flow, a local segregation flux $\boldsymbol{F}^{\nu }$ that drives different species to different regions of the flow and a diffusive flux $D \boldsymbol{\nabla }\phi ^{\nu }$ that accounts for the fluctuating motion of the grains (Gray & Chugunov Reference Gray and Chugunov2006), where $D$ is the binary diffusion coefficient.

To close the segregation model, functional forms for the segregation fluxes must be provided, which satisfy a number of constraints. Firstly, in order to satisfy the bulk incompressibility (4.1) the segregation fluxes of large and small particles must sum to zero:

(4.5) \begin{equation} \boldsymbol{F}^s+\boldsymbol{F}^l=\boldsymbol{0}. \end{equation}

This means only a single segregation flux needs to be specified, as the other will be given by (4.5). In this paper, the small-particle segregation flux $\boldsymbol{F}^s$ is specified. A second constraint requires that segregation shuts off when either constituent is in pure phase. Thus $\boldsymbol{F}^s$ must vanish when $\phi ^s=0$ and $\phi ^s=1$ . The simplest flux function satisfying this constraint was postulated by Dolgunin & Ukolov (Reference Dolgunin and Ukolov1995) and has a quadratic dependence on $\phi ^s$ :

(4.6) \begin{equation} \boldsymbol{F}^s=\phi ^s(1-\phi ^s)\boldsymbol{q}, \end{equation}

where vector $\boldsymbol{q}$ has dimensions of velocity and sets the magnitude and the direction of segregation. Experiments, DEM simulations and theoretical analyses of bidisperse mixtures have revealed that $\boldsymbol{q}$ depends on the diameter ratio $R=d^l/d^s$ between the grains (Thornton et al. Reference Thornton, Weinhart, Luding and Bokhove2012), small-particle concentration $\phi ^s$ (van der Vaart et al. Reference van der Vaart, Gajjar, Epely-Chauvin, Andreini, Gray and Ancey2015), shear rate $\dot {\gamma }$ (Bridgwater, Foo & Stephens Reference Bridgwater, Foo and Stephens1985; Savage & Lun Reference Savage and Lun1988; Marks, Rognon & Einav Reference Marks, Rognon and Einav2012; Fan et al. Reference Fan, Schlick, Umbanhowar, Ottino and Lueptow2014), pressure $p$ (Golick & Daniels Reference Golick and Daniels2009; Fry et al. Reference Fry, Umbanhowar, Ottino and Lueptow2019; Trewhela et al. Reference Trewhela, Ancey and Gray2021) as well as gradients of pressure (Gray & Thornton Reference Gray and Thornton2005; Staron & Phillips Reference Staron and Phillips2015), kinetic stress (Hill & Tan Reference Hill and Tan2014) and shear rate (Jing et al. Reference Jing, Ottino, Lueptow and Umbanhowar2021):

(4.7) \begin{equation} \boldsymbol{q}=\boldsymbol{q}(\phi ^s,R,\dot {\gamma },p, \boldsymbol{\nabla }p, \boldsymbol{\nabla }\dot {\gamma }, \ldots ). \end{equation}

A completely general functional form of $F^s$ is still subject to debate, but so far, models of size-bidisperse grains of equal density have only captured the conventional segregation shown in figure 1(a). This paper presents flux functions with a concentration dependence that allows a unified description of both reverse and standard segregation of large particles.

We examine segregation in a free-surface flow over an inclined plane. To facilitate comparison with the DEM simulations, in which periodic boundary conditions are imposed in the down-slope ( $x$ ) and cross-slope ( $y$ ) directions, the flow is assumed to be uniform in the $x$ and $y$ directions, with the only flow gradients occurring in the slope-normal $z$ direction. Under these assumptions the general segregation equation (4.4) simplifies to

(4.8) \begin{equation} \frac {\partial \phi ^s}{\partial t}+\frac {\partial F^s}{\partial z} =\frac {\partial }{\partial z}\!\left (D\frac {\partial \phi ^s}{\partial z}\right )\!, \end{equation}

where $F^s$ is the $z$ component of $\boldsymbol{F}^s$ . This is subject to the boundary condition of no flux of particles through the base at $z=0$ and the free surface at $z=h$ :

(4.9) \begin{equation} F^s = D \frac {\partial \phi ^s}{\partial z} \quad \text{on} \quad z=0, h. \end{equation}

Integrating (4.8) between the base and free surface and applying (4.9) implies

(4.10) \begin{equation} \frac {\partial \overline {\phi ^s}}{\partial t} = 0, \end{equation}

where

(4.11) \begin{equation} \overline {\phi ^s}=\frac {1}{h}\int _0^h \phi ^s \,\mathrm{d} z \end{equation}

is the mean fraction of small particles throughout the flow depth. Seeking steady solutions to the partial differential equation (4.8) simplifies it to an ordinary differential equation (ODE). Integrating (4.8) in $z$ and applying (4.9) then gives

(4.12) \begin{equation} \frac {\mathrm{d} \phi ^s}{\mathrm{d} z}=\frac {F^s}{D}. \end{equation}

This expresses the gradient of the concentration profile as a balance between segregation and diffusion. Integrating (4.12) once more to obtain a unique solution for $\phi ^s(z)$ requires an additional constraint, for which the constant $\overline {\phi ^s}$ is specified. It also requires explicit expressions for $F^s$ and $D$ , which are the focus of § 5.

The ODE (4.12) subject to the integral constraint (4.11) is solved by introducing the potential $\psi$ , defined as

(4.13) \begin{align} \psi (z) = \frac {1}{h}\int _{0}^{z} \phi ^s(z') \: \mathrm{d}z'. \end{align}

This satisfies the ODE

(4.14) \begin{align} \frac {\mathrm{d} \psi }{\mathrm{d}z} = \frac {\phi ^s}{h}, \end{align}

and boundary conditions

(4.15) \begin{align} \psi (0) = 0 \quad \text{and} \quad \psi (h) = \overline {\phi ^s}. \end{align}

Equations (4.12) and (4.14) and boundary conditions (4.15) form a two-point boundary-value problem for $\phi ^s$ , which is solved by shooting.

5.1. Neutral concentration

Recall that at size ratios $R \geq 4$ , large particles tend to sink to the base of the flow when present in low concentrations (reverse segregation). This contrasts with the segregation behaviour at $R \leq 4$ where large particles rise to the surface irrespective of their concentration, leading to the large particles overlaying smaller particles. Current continuum models for segregation have been constructed for $R \leq 4$ ; consequently their small-particle segregation fluxes $F^s$ are always negative to ensure small particles are transported towards the base, and large particles rise to the surface. This renders them incapable of predicting segregation at larger size ratios. To address reverse segregation, a form $F^s$ is proposed that becomes negative for high concentrations of small particles, at size ratios large enough for reverse segregation to occur.

To begin, the small-particle segregation flux $F^s$ , in the limit of a single small or single large intruder, is examined. These limits correspond to the small-particle concentrations $\phi ^s=0^{+}$ and $\phi ^s=1^{-}$ , respectively, where the superscripts ‘ $+$ ’ and ‘ $-$ ’ denote concentrations just above or below 0 and 1.

A single small intruder ( $\phi ^s=0^{+}$ ) will percolate downwards through a shearing matrix of large particles under the action of gravity, resulting in a downward flux of small particles, $F^s\lt 0$ . Conversely, for $R \gtrsim 4$ a single large intruder ( $\phi ^s=1^-$ ) will sink under gravity (Thomas & D’Ortona Reference Thomas and D’Ortona2018). This results in a downward flux of large particles. By the mass conservation condition (4.5), there must be an upward flux of small particles relative to the bulk flow, $F^s\gt 0$ . In particular, there is a change in sign of the segregation flux $F^s$ as the concentration moves from zero to one. Importantly, this suggests that there is some intermediate concentration, denoted the neutral concentration, $\phi ^s_{n}$ , at which $F^s=0$ :

(5.1) \begin{align} F^s(\phi ^s_{n})=0. \end{align}

The above reasoning constrains $F^s$ to have at least one root. However, $F^s$ must also vanish in the monodisperse limits $\phi ^s=0,1$ , where segregation shuts off. The simplest form of $F^s$ that has three roots at $\phi ^s=0$ , $\phi ^s_n$ and 1 has a cubic concentration dependence:

(5.2) \begin{align} F^s=q\phi ^s(1-\phi ^s)(\phi ^s_n-\phi ^s), \end{align}

where $q$ is a negative parameter with dimensions of velocity that sets the rate of segregation, analogously to (4.6). The parameter $q$ may also depend on flow variables and particle properties as per (4.7).

The neutral concentration $\phi ^s_n$ is the concentration at which segregation reverses. It determines whether large particles rise or sink within the flow. For $\phi ^s\lt \phi ^s_n$ the large particles rise to the surface. For $\phi ^s\gt \phi ^s_n$ the large particles sink towards the base. Between these cases there lies a third, intermediate regime where $\phi ^s=\phi ^s_n$ and $F^s=0$ , in which there is no net movement of large particles. This corresponds to $\phi ^s \approx 0.6$ in figure 2(b) where $\phi ^s$ becomes almost uniform with depth.

Figure 13. Schematic plots of (a–c) the small-particle segregation flux, $F^s$ , as a function of small-particle concentration, $\phi ^s$ . Curve (a) corresponds to the symmetric, quadratic model (4.6); curve (b) shows a strictly negative cubic flux based on (5.2) with $\phi ^s_n = 1.1$ ; and curve (c) depicts a bidirectional cubic flux from (5.2) with $\phi ^s_n = 0.6$ . All fluxes plotted here share the same maximum magnitude, $|F^s| = 0.25$ ; the dependence of this magnitude on flow variables is discussed separately in § 5.2. (d–f) The normal segregation velocities of small (5.5) and large (5.6) particles, corresponding to the fluxes in (a–c), respectively. The black circles in (c, f) mark the neutral concentration $\overline {\phi ^s_n}=0.6$ , the concentration at which segregation reverses.

The value of $\phi ^s_n$ will depend upon the size ratio $R$ , since the small-particle concentration at which segregation reverses decreases with $R$ (figure 4, and Thomas & D’Ortona Reference Thomas and D’Ortona2018). In particular, for $R \lesssim 4$ no reversal in segregation has been observed (Thomas & D’Ortona Reference Thomas and D’Ortona2018), so the segregation flux $F^s$ must remain negative on the interval $\phi ^s \in (0,1)$ . From the form (5.2) of $F^s$ , this requires $\phi ^s\lt \phi ^s_n$ throughout this range. Furthermore, given that the direction of segregation does reverse for $R \gtrsim 4$ , the neutral concentration must obey

(5.3) \begin{align} \phi ^s_n \gt 1 \quad &\mbox{for } \quad 1 \leq R \lesssim 4 \end{align}

and

(5.4) \begin{align} 0 \lt \phi ^s_n \leq 1 \quad &\mbox{for } \quad R \gtrsim 4. \end{align}

The introduction of the neutral concentration $\phi ^s_n$ has a significant impact on the segregation flux (figure 13). When $\phi ^s_n\gt 1$ the segregation flux is skewed relative to the symmetric, quadratic model (figure 13 a) with its minimum shifted towards lower small-particle concentrations (figure 13 b). As $\phi ^s_n$ decreases towards unity, $F^s$ becomes increasingly asymmetric. When $\phi ^s_n$ drops below one, a new root appears at $\phi ^s=\phi ^s_n$ (figure 13 c), marking a qualitative change in behaviour. In this regime $F^s$ can take either sign: it is negative for $\phi ^s\lt \phi ^s_n$ , indicating small particles moving downwards, and it is positive for $\phi ^s\gt \phi ^s_n$ , representing small particles being pushed upwards. This is in stark contrast to the case with $\phi ^s_n\gt 1$ , where $F^s$ is always negative to ensure small particles are transported towards the base.

Asymmetric segregation fluxes such as that shown in figure 13(b) have previously been proposed by Gajjar & Gray (Reference Gajjar and Gray2014) in order to model the asymmetry in segregation rates between small and large particles. Specifically, in the regime where large particles rise to the surface, a lone small particle will typically migrate downwards through a matrix of large particles more quickly than a large particle rises through a matrix of small particles (van der Vaart et al. Reference van der Vaart, Gajjar, Epely-Chauvin, Andreini, Gray and Ancey2015). This asymmetry is not captured by the quadratic flux (4.6), which is symmetric about a 50 % concentration of small particles (figure 13 a).

The effect of skewing the flux can be seen in the small- and large-particle segregation velocities, defined by

(5.5) \begin{eqnarray} w^{s}_{\textit{seg}}=w^{s}-w=\frac {F^s}{\phi ^s}=q(1-\phi ^s)(\phi ^s_n-\phi ^s), \end{eqnarray}

(5.6) \begin{eqnarray} w^{l}_{\textit{seg}}=w^{l}-w=-\frac {F^s}{1-\phi ^s}=-q\phi ^s(\phi ^s_n-\phi ^s). \end{eqnarray}

Here, $w^\nu$ is the slope-normal, $z$ , component of the velocity of species $\nu$ and $w=\phi ^sw^s+\phi ^lw^l$ is the bulk $z$ velocity. A positive segregation velocity ( $w^{\nu }_{\textit{seg}} \gt 0$ ) indicates that species $\nu$ is segregating upwards, while a negative value indicates downward segregation. For the symmetric flux function, the segregation velocity for a single small intruder ( $\phi ^s=0^{+}$ ) and the segregation velocity for a single large intruder ( $\phi ^s=1^{-}$ ) are equal in magnitude (figure 13 d). In the skewed case (figure 13 e), however, the segregation velocity of a small intruder is larger in magnitude than that of a large intruder, reflecting the asymmetry in the underlying flux.

Gajjar & Gray (Reference Gajjar and Gray2014) restricted their model to negative $F^s$ to ensure that small particles always segregate downwards, which corresponds to the neutral concentration taking values $\phi ^s_n \geq 1$ . Allowing the neutral concentration $\phi ^s_n$ to take values on the interval $(0,1)$ introduces a bidirectional flux with segregation velocities that change sign at $\phi ^s = \phi ^s_n$ (figure 13 f). In this regime, both small and large intruders exhibit downward segregation, unlike the case where $\phi ^s_n \geq 1$ . Thus, the cubic flux form (5.2) may be considered as a natural extension of the model of Gajjar & Gray (Reference Gajjar and Gray2014) to larger size ratios, that is able to capture both the asymmetry in the segregation rates and the reversal in the segregation direction, depending on the value of the neutral concentration $\phi ^s_n$ . Although a cubic dependence on $\phi ^s$ is the simplest form capable of capturing either the asymmetry in segregation rates or a reversal in segregation direction, the segregation flux $F^s$ may, in general, be a more complicated function of $\phi ^s$ .

The proposed model (5.2) has a mathematical structure similar to combined size and density models that allow for a concentration-dependent reversal in the segregation direction (Gray & Ancey Reference Gray and Ancey2015; Duan et al. Reference Duan, Umbanhower, Ottino and Lueptow2021). In these models, large particles that are denser than the small ones may sink when present in low concentrations but rise when present in high concentrations. This reversal arises because different segregation mechanisms dominate at different concentrations: at low concentrations, buoyancy effects make isolated dense large particles sink, while at high concentrations, the downward migration of small particles through voids drives the dense large particles upward. The analogous mathematical form of the segregation flux may reflect a similar underlying mechanism for the reversal at large size ratios. An isolated large particle, even if it has the same intrinsic density as the small particles, can sink because it is effectively denser, relative to the surrounding mixture of small particles and void space. At higher concentrations of large particles, however, this effective density contrast is reduced, and the segregation is instead controlled by small particles falling downward through voids, which displaces the large particles upward.

5.2. Flow dependence

5.2.1. Diffusion

In steady state, (4.12) implies that the concentration profile is determined by a balance between the local segregation flux and the diffusive flux. Diffusion in granular flows is driven by shear, with regions of higher shear rate experiencing larger fluctuations leading to enhanced diffusion. In monodisperse flows, the grains undergo Fickian (normal) diffusion, and dimensional analysis requires the diffusion coefficient $D$ to be of the form (Bridgwater Reference Bridgwater1980; Savage & Dai Reference Savage and Dai1993; Utter & Behringer Reference Utter and Behringer2004)

(5.7) \begin{equation} D = \dot \gamma d^2 D_0(I), \end{equation}

where $D_0$ is, in general, an arbitrary function of the inertial number

(5.8) \begin{equation} I=\frac {\dot {\gamma }d}{\sqrt {p/\rho _*}}. \end{equation}

For $I \geq 0.1$ the dependence of $D_0$ on $I$ is weak; however, at very low inertial numbers the diffusivity increases more strongly, scaling as $D \propto I^{-1/2}$ (Kharel & Rognon Reference Kharel and Rognon2018).

Cai et al. (Reference Cai, Xiao, Zheng and Zhao2019) extended expressions for self-diffusion in monodisperse flows to bidisperse systems up to a size ratio $R=3$ by replacing the particle diameter $d$ with the volume-weighted mean particle diameter $\bar {d}$ :

(5.9) \begin{equation} \bar {d}=\phi ^s d^s +\phi ^l d^l. \end{equation}

This ensures that when the flow is dominated by small particles ( $\phi ^s=1$ ), $\bar {d}$ reduces to the small-particle diameter $d^s$ , and when the flow is dominated by large particles ( $\phi ^s=0$ ), $\bar {d}$ equals the large-particle diameter $d^l$ . This results in an average diffusion coefficient

(5.10) \begin{equation} D= D_0\dot {\gamma }(\bar {d})^2, \end{equation}

where $D_0\in (0.03,0.05)$ is approximately constant, although there is a weak dependence on the volume fraction (Cai et al. Reference Cai, Xiao, Zheng and Zhao2019). This is consistent with the results of multiple studies using monodisperse grains (Bridgwater Reference Bridgwater1980; Savage & Dai Reference Savage and Dai1993; Campbell Reference Campbell1997; Bancroft & Johnson Reference Bancroft and Johnson2021). In this study, a diffusion coefficient $D_0=0.04$ is used, as in Maguire et al. (Reference Maguire, Barker, Rauter, Johnson and Gray2024).

5.2.2. Segregation flux

Up to this point, only the minimal concentration dependence of the segregation flux required to capture the reverse segregation of large particles has been considered. However the segregation flux may have further dependence on $\phi ^s$ as well as on flow variables and particle properties. Concentration dependence alone has been shown to provide an effective approximation to the segregation profiles at $R=2$ (Wiederseiner et al. Reference Wiederseiner, Andreini, Epely-Chauvin, Moser, Monnereau, Gray and Ancey2011); however, a more complete segregation model requires incorporating additional dependences. These enter the segregation flux through the segregation rate $q$ in (5.2).

Through the use of shear box experiments to examine the single intruder limit of segregation, Trewhela et al. (Reference Trewhela, Ancey and Gray2021) determined that $q$ is linear in the shear rate $\dot {\gamma }$ and inversely proportional to the pressure $p$ . This relationship can be primarily attributed to the drag on segregating particles (Bancroft & Johnson Reference Bancroft and Johnson2021; Jing et al. Reference Jing, Ottino, Lueptow and Umbanhowar2022). In free-surface flows over an inclined plane, which is the configuration considered in this paper, the normal segregation force due to the hydrostatic pressure gradient dominates over that due to the shear rate gradient (Jing et al. Reference Jing, Ottino, Lueptow and Umbanhowar2021). These findings suggest that the segregation flux can be expressed as

(5.11) \begin{equation} F^s=\frac {\dot {\gamma } (\bar {d})^2}{p} \frac {{\rm d}p}{{\rm d}z} f(\phi ^s,R), \end{equation}

where the hydrostatic pressure for a flow of height $h$ on a plane inclined at angle $\theta$ is given by

(5.12) \begin{equation} p=\rho _{*}\varPhi g\cos \theta (h-z), \end{equation}

and where $\bar {d}$ is chosen as a reference length scale (Trewhela et al. Reference Trewhela, Ancey and Gray2021). The function $f(\phi ^s,R)$ encapsulates the concentration dependence of the segregation flux, as was discussed in § 5.1. It is expected that $f$ is positive for $R \lesssim 4$ to ensure that small particles migrate downwards. To model reverse segregation, $f$ must be at least a cubic in $\phi ^s$ , though higher-order forms are possible. In this paper a cubic form is primarily considered, as it provides the simplest model with the necessary properties:

(5.13) \begin{equation} f_3(\phi ^s,R)=A\phi ^s(1-\phi ^s)(\phi ^s_n-\phi ^s). \end{equation}

This function vanishes at $\phi ^s=0,1$ as required, and is characterised by two dimensionless parameters: a constant overall amplitude $A$ which scales the strength of segregation and the neutral concentration $\phi ^s_n$ which sets the point of reversal. Both parameters depend on the size ratio $R$ and are determined by fitting to the data from DEM simulations. For comparison, a quintic form of $f$ is also considered:

(5.14) \begin{equation} f_5(\phi ^s,R)=B_1\phi ^s(1-\phi ^s)(B_2-\phi ^s)(1+B_3\phi ^s+B_4{\phi ^s}^2), \end{equation}

which allows quantitative concentration profiles to be captured more accurately. However, this comes at the expense of using four dimensionless parameters compared with two for the cubic model. These parameters are denoted $B_1$ , $B_2$ , $B_3$ and $B_4$ and depend on the size ratio $R$ . The amplitude of the flux function is controlled by $B_1$ , while $B_2$ sets the concentration at which segregation reverses. The remaining parameters, $B_3$ and $B_4$ , adjust the curvature of the flux function, enabling a closer match to DEM profiles.

5.3. Steady state

Substituting the forms of the segregation flux and diffusive flux into the equation for the steady-state concentration profiles (4.12) results in

(5.15) \begin{equation} \frac {\mathrm{d} \phi ^s}{\mathrm{d}z}=\frac {-f(\phi ^s,R)}{D_0(h-z)}. \end{equation}

Since both the segregation flux $F^s$ and the diffusive flux are linear in the shear rate, the shear rate dependence is removed in steady state, and thus the steady-state profiles are not dependent on the rheology. This equation is integrated between the base $z=0$ and the free surface $z=h$ subject to the constraint (4.11) on the depth-averaged concentration of small particles.

By inspection, one solution of (5.15) is $\phi ^s = \phi ^s_n$ , corresponding to well-mixed particles with a uniform vertical distribution. The model therefore predicts the complete suppression of segregation in such mixtures, corresponding to the vertical dashed line in figure 2(b), which is approximately equal to the results of the DEM simulation at $\overline {\phi ^s}=0.6$ through the majority of the flow depth. This divides other solutions of (5.15) into two classes: those for which $\phi ^s \lt \phi ^s_n$ for all $z$ and those for which $\phi ^s \gt \phi ^s_n$ for all $z$ . For both classes, there is a singularity in (5.15) at $z=h$ , which arises because the pressure vanishes at the free surface, causing the segregation flux $F^s$ to diverge as $z$ approaches $h$ . This divergence in $F^s$ forces each species into pure phase at the free surface ( $\phi ^s=0$ or $\phi ^s=1$ at $z=h$ ), a prediction in good agreement with the results of the DEM simulations in figure 2.