# Rayleigh-Taylor instability in size-bidisperse isodense granular flow down an incline

Umberto d'Ortona, Richard Lueptow, Nathalie Thomas

► To cite this version:

Umberto d'Ortona, Richard Lueptow, Nathalie Thomas. Rayleigh-Taylor instability in size-bidisperse isodense granular flow down an incline. 2025. hal-05233035

HAL Id: hal-05233035

https://hal.science/hal-05233035v1

Preprint submitted on 1 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Rayleigh-Taylor instability in size-bidisperse isodense granular flow down an incline

Umberto D'Ortona*

*Aix-Marseille Université, CNRS, Centrale Med., M2P2, Marseille, France*

Richard M. Lueptow

*Department of Mechanical Engineering, Northwestern University, Evanston, Illinois 60208, USA and*
*Department of Chemical and Biological Engineering,*
*Northwestern University, Evanston, Illinois 60208, USA*

Nathalie Thomas

*Aix-Marseille Université, CNRS, IUSTI, Marseille, France*

(Dated: September 1, 2025)

Size-bidisperse granular material flowing down a rough incline, or chute, may develop a Rayleigh-Taylor instability even when the two particles species have the same density. Unlike fluid instability, which occurs under quiescent conditions, the granular material must be flowing for the instability to occur. Initially, size segregation results in a monodisperse layer of large particles above a monodisperse layer of small particles with an interfacial layer of mixed particles that has a higher volume fraction and, hence, is denser than the monodisperse layers above and below it. These two layers destabilize via a Rayleigh-Taylor instability and form dense descending plumes of mixed particles and ascending plumes of less dense pure small particles. As a result of these plumes, the upper layer of large particles breaks and accumulates above the descending plumes; the small particle plumes reach the free surface. The instability evolves into persistent longitudinal rolls corresponding to streamwise bands, or stripes, of small and large particles at the free surface. The appearance of the instability depends on the particle and mixture properties as well as the flow conditions. In all cases, the propensity for the appearance of the instability and subsequent band formation can be traced back to the total particle volume fraction, or packing density, in the mixed particle layer, which depends on the particle size ratio, fraction of each particle species, and thickness of the flow.

## I. INTRODUCTION

The Rayleigh-Taylor instability occurs at the horizon-
tal interface between a layer of dense fluid above a layer
of less dense fluid in a gravitational field [1–4], resulting
in initial deformation of the interface between the layers
leading to downward plumes of heavy fluid and upward
plumes of light fluid. The instability can be observed over
a wide range of scales from that in a laboratory to the
explosion of a supernova, where acceleration plays the
role of gravity [5, 6]. This instability has been studied
experimentally [7], numerically [8, 9], and theoretically
over a wide range of conditions using miscible or immis-
cible fluids [10–12], quiescent or flowing systems [13], fi-
nite or infinite systems [14–16], and bilayer or multilayer
systems [13, 17].

The existence of a granular Rayleigh-Taylor instabil-
ity in flows of dry non-cohesive granular media has only
recently been described [18, 19]. When a layer of dense
particles is placed on top of a layer of less dense particles
of the same size, a Rayleigh-Taylor instability (RTI) oc-
curs provided that the particles are mobilized via a flow,
for example down on a rough inclined plane. The insta-
bility is evident as a spanwise wavy deformation of the
interface between the layers in the plane perpendicular
to the direction of flow, and it becomes visible at the

surface of the flow as a series of alternating streamwise
bands of dense and light particles.

More surprisingly, if the two particle species differ in
size in addition to density, with the dense particles being
larger, a self-induced RTI can occur [18]. In this case,
the two particles species are initially a mixed homoge-
neous system with no layers, and therefore a priori sta-
ble. When set into motion, say flowing down an incline,
large dense particles will form a layer at the surface due
to granular size segregation, which will then destabilize
in a second step via the RTI. Whether starting from an
initially homogeneous mixture of large dense and small
light particles or from an initially size-segregated bilayer
configuration, the alternating bands are maintained by
convective rolls similar to Rayleigh-Bénard cells, where
the driving force for convection is granular segregation
rather than temperature.

Even more surprisingly, a granular RTI can occur even
in an *isodense* size-bidisperse system, where the two
species of particles have the same density and differ only
in size. This counter-intuitive phenomenon comes about
because particles of two different sizes pack more densely
than either pure species. Upon size segregation during
flow, a layer of mixed particles forms at the interface be-
tween the segregated upper layer of large particles and
the lower layer of small particles. UDO July 24: When
we start with a bilayer system, which is mainly the case
here, it is not segregation that forms the mixed layer, but
mixing, and segregation limits it. So we propose: As a

<footnote>* umeonto.dortona@cnrs.fr</footnote>
result of a balance between size segregation and mixing during flow, a layer of mixed particles forms at the interface between the segregated upper layer of large particles and the lower layer of small particles. A granular RTI results from the dense mixed particle interface layer above the less dense monodisperse layer of small particles. In fact, granular axial size segregation and band formation in a partially-filled cylindrical rotating drum (or tumbler), a phenomenon discovered over 80 years ago but never having a satisfactory mechanism to explain it, can be understood in terms of the granular RTI [20].

The objective of this article is to better understand the granular RTI by considering flow down an incline, which is a simpler geometry than the rotating drum. Unlike the classical Rayleigh-Taylor instability with a static layer of dense fluid above a static layer of less dense fluid, the dense layer in the granular case is directly connected to the dynamics of the flow along the incline. The density and thickness of the unstable dense mixed-particle layer at the interface between flowing size-segregated layers of the two species of particles depends on the combined effects of mixing due to collisional diffusion induced by the granular flow and segregation induced by differences in the sizes of the two particle species. This dependence of the dense layer on the dynamics of the flow is unique for a RTI.

This article is organized as follows. Section II presents the Discrete Element Method (DEM) [21, 22], and Section III explains the granular Rayleigh-Taylor instability phenomenon for isodense size-bidisperse particles. Section IV provides an in-depth examination of the origin of the dense layer of mixed-size isodense particles at the origin of the granular RTI, and Section V explores how parameters such as flow thickness, particle volume fraction, particle size ratio, and friction coefficient affect the RTI instability and subsequent spanwise segregation into bands of large and small particles. The article concludes with a discussion of the findings and future perspectives.

## II. NUMERICAL METHOD

The open-source DEM code LIGGGHTS [23] is used to explore size-bidisperse granular flow on a rough incline using periodic boundary conditions in the streamwise $x$ and spanwise $y$ directions of the computational domain. The small particles (blue) have a diameter $d = 2$ mm, while the size of the large particles, $d_l$ (red), depends on the size ratio, $1.25 \le d_l/d \le 3.25$. The two particle species have identical properties that correspond to cellulose acetate: density $\rho = 1308$ kg m$^{-3}$ and restitution coefficient $e = 0.87$ [24–26]. Particles have uniform size distributions of 0.95-1.05 times $d$ or $d_l$ to prevent ordering in the flow. These properties are identical to our previous study on axial segregation in a rotating drum [20] to facilitate comparison. Particle interactions are calculated using a linear contact model with Young’s modulus $E = 10^7$ Pa and Poisson’s ratio $\nu = 0.3$ for the normal

forces, and tangential forces include a Coulomb friction with a friction coefficient varied as $0.2 \le \mu \le 1.5$ and a tangential linear spring that deforms at first contact between particles. As the simulations are lengthy, we often use $d_l/d = 3$ with $\mu = 0.3$, because spanwise segregation appears rapidly under these conditions. The integration time step is $5 \times 10^{-6}$ s to provide stable results.

The rough inclined plane is made up of randomly placed particles of diameter $d$ fixed to a smooth planar surface in a single layer with a coverage of about 57%, which is well below hexagonally close packed particle coverage of about 91% to ensure a high level of roughness. The numerical domain is typically $L = 200d$ long and $W = 200d$ wide in the streamwise and spanwise directions, respectively. The width may be increased to as much as $W = 400d$ to measure the wave length more accurately or reduced to $W = 30d$ to prevent any spanwise segregation. When numerous or long-term simulations are required, the length is reduced to $L = 60d$ to save computing time. The flowing particles are initially placed on the inclined plane in either a bilayer configuration with the layer of large particles above the layer of small particles, or as a single layer of homogeneously mixed particles.

To measure the growth, and, potentially, the decay, of the instability, a sinusoidal perturbation can be introduced at the interface between large and small particles. In this scenario, the width of the simulation box is reduced to approximately one spanwise wavelength to permit the appearance of only one pair of streamwise bands, one of mostly large particles and one of mostly small particles. A $t = 0$, gravity is set to $g = 9.81$ m s$^{-2}$ at an angle of 23°, corresponding to the incline of the plane. The flow depth, $H$, which is measured from the top of the $d = 2$ mm particles forming the rough surface of the inclined plane, is varied from 12d to 60d. The velocity is zero at the no-slip surface of the inclined plane and increases nearly linearly with distance above the surface for $H = 12d$ but is more Bagnold-like for larger $H$.

To quantify the evolution of the instability, we define spanwise and vertical segregation indices, $I_{\text{seg}}^s$ and $I_{\text{seg}}^v$. The vertical segregation index,

$$I_{\text{seg}}^{v} = 2 \frac{\text{CM}_l - \text{CM}_s}{H}, \quad (1)$$

computes the difference in the vertical positions of the centers of mass ($\text{CM}_l$) for the large and small particles. It varies from 0 for a homogeneously mixed system to 1 for a perfectly vertically segregated system of two distinct, non-overlapping layers. This range of values is independent of the concentration of large and small particles. The spanwise segregation index is defined as

$$I_{\text{seg}}^{s} = \frac{1}{4F_{l}F_{s}WHL} \int_{0}^{H} \int_{0}^{W} \int_{0}^{L} |f_{l}(x, y, z) - F_{l}| \quad (2) \\ + |f_{s}(x, y, z) - F_{s}| \, dx \, dy \, dz,$$

where $F_l$ and $F_s$ are the global volume concentrations of large and small particles ($F_s + F_l = 1$), and $f_l$ and $f_s$ are
the local volume concentrations ($f_s + f_l = 1$). $I_{seg}^s$ ranges from 0 for no spanwise segregation to 1 for perfectly segregated streamwise bands of each particle species, and it is equivalent to the axial segregation index, $I_{seg}$, used previously for axial segregation in a rotating drum [20]. In the case of an idealized spanwise sinusoidal perturbation of the large-small particle interface of amplitude $B$ in the $z$-direction, $B$ is directly linked to $I_{seg}^s$ as

$$ \frac{B}{H} = \pi F_l F_s I_{seg}^s. \qquad (3) $$

As a consequence, for an exponential growth of the amplitude $B = B_0 \exp \gamma t$, the growth rate $\gamma$ can be directly measured by plotting $I_{seg}^s$ versus $t$ on a semi-log scale.

### III. THE PHENOMENON

#### A. Formation of segregated bands

Figure 1 highlights the phenomenon of interest: the establishment of spanwise segregated alternating streamwise bands of large ($d_l = 6$ mm, red) and small ($d = 2$ mm, blue) isodense particles ($d_l/d_s = 3$), whether starting from layers of the two particle species or from a homogeneous mixture of the two particle species. The $H = 24d = 48$ mm thick flow is in the $x$-direction, from left to right, down a rough surface inclined at $\theta = 23^\circ$. The simulation domain is $L = 200d$ in the periodic direction of flow, and $W = 200d$ in the periodic spanwise direction. The global volume concentrations are $F_l = 0.2$ and $F_s = 0.8$.

Figure 1(a-e) shows successive instants of the system starting with a segregated bi-layer of pure large red particles above pure small blue particles. As time progresses, streamwise bands begin to emerge between $t = 20$ s and $30$ s, and become more pronounced by $t = 40$ s. At steady state (around $t \approx 80$ s), a spanwise pattern of alternating streamwise bands composed of nearly pure large particles and bands composed predominantly of small particles but mixed with a few large particles is observed at the free surface.

It is also possible to start with a homogeneously mixed system of isodense size-bidisperse particles [Fig. 1(f)] and achieve a similar set of streamwise bands of large and small particles [Fig. 1(g, h)]. Starting from a homogeneous configuration, granular size segregation first results in a layer of large particles over a layer of small particles on the inclined surface that will later develop the alternating band pattern. This occurs more slowly than for the initially two-layer case, because the two layers must first form (compare (d) and (g) at $t = 40$ s). But the initially homogeneous case eventually reaches a stationary regime of spanwise-alternating bands of large and small particles [Fig. 1(h)], comparable to the results for a two-layer initial condition [Fig. 1(e)].

The band formation in Fig. 1 is similar to that in the flowing layer of isodense size-bidisperse particles in the

<figure><img src="image_2.png" /><figcaption>FIG. 1. Typical examples of spanwise segregation in a $W = L = 200d$ periodic domain forming streamwise bands of large and small particles in a $H = 24d = 48$ mm thick flow on an inclined plane ($\theta = 23^\circ$) with particle size ratio $d_l/d = 3$ and global species concentrations $F_l = 0.2$ and $F_s = 0.8$. (a-e) The system is initially two separate layers with large red particles above. (f-h) The system is initially homogeneously mixed particles. Flow is from left to right.</figcaption></figure>

flowing layer of a partially filled rotating drum, which eventually leads to axial (spanwise) bands of large and small particles due to the RTI [20]. Although the granular material in a rotating drum typically starts in a homogeneous state, size segregation is highly efficient, so that within a few drum rotations a flowing layer of pure large particles above a layer of pure small particles results, with a mixed layer in between. Since the characteristic timescales for layer formation and band formation in a rotating drum are so different, it is reasonable to consider the initial two-layer state on an inclined plane as analogous to the initial condition for axial segregation in a rotating drum. However, axial segregation in a rotating drum emerges much more slowly than that in Fig. 1(a-e), requiring hundreds of drum rotations or more, probably because particles are constantly exiting and reentering the flowing layer. The advantage of the inclined plane for studying the granular RTI phenomenon is that it is not contaminated by this problem, thereby isolating the Rayleigh-Taylor instability.
## B. The granular RTI mechanism

<figure><img src="image_3.png" /><figcaption>FIG. 2. Large particle species concentration, $f_l$ (red-blue colorbar), and particle volume fraction, $\phi$ (yellow-black colorbar), in a cross-section normal to the flow at (a, d) $t = 5$ s, (b, e) $t = 20$ s and, (c, f) $t = 80$ s. Same conditions as Fig. 1(a-e).</figcaption></figure>

To more clearly demonstrate the granular RTI mechanism, consider Fig. 2, which shows a vertical ($y, z$)-cross-section, perpendicular to the flow, of an initially bi-layer granular flow for the conditions in Fig. 1(a-e). In this figure and similar figures later in this article, the local values for concentration $f_l$ and overall particle volume fraction, or packing density, $\phi$, are based on averages over streamwise length $L$ and over 1 s. At $t = 5$ s in Fig. 2(a), there is a thin layer of large particles (red) above a thicker base layer of small particles (blue) with a thin layer of mixed particles (white) between them. The mixed layer of particles has a larger volume fraction than the small particles below it, as shown in Fig. 2(d), where darker color corresponds to high volume fraction. The larger density in the mixed layer over the lower density of the bottom pure layer initiates the Rayleigh-Taylor instability, and perturbations of the mixed layer appear in Fig. 2(b, e) by $t = 20$ s as variations in the thickness of the mixed and upper layers. The density difference between the mixed layer composed of particles of different sizes and the lower layer of small particles results in the initiation of sinking dense plumes of mixed particles and rising light plumes of pure small particles. The rising plumes eventually reach the free surface, forming alternating bands of large and small particles at the surface of the flow [Fig 2(c, f)]. The large particle bands are shallow with small particles below them, while the small particle bands include some large particles and extend from the free surface to the base of the flow. We note the

appearance of a layering pattern of large particles, evident as faint horizontal strips below the surface bands of large particles. This is a separate ordering phenomenon that occurs for flows of size-bidisperse particles at large size ratios [28, 29].

Compared to a classical Rayleigh-Taylor instability, no inversion of the two phases is observed, nor is there the characteristic mushroom-shaped pattern. Here, the dense phase plumes that descend remain "frozen" at the surface level. A similar phenomenon is observed in the case of axial segregation in a rotating cylinder. This is related to the fact that the dense phase is not a fluid in the strict sense of the term, but rather results from a balance between mixing and segregation of the two other phases.

Before going into detail about the excess density at the interface between the large particle layer and the small particle layer that is intrinsically linked to size-disperse granular media, it is useful to describe the conditions for the onset of the RTI in two- and three-layer fluid systems, inspired by the work of Batchelor and Nitzche on multi-layer fluid systems [15]. The textbook version of the RTI

<figure><img src="image_4.png" /><figcaption>FIG. 3. Mechanisms governing the Rayleigh-Taylor instability for two layers (a, g) and for a variety of situations with three layers having a dense middle layer (b-f), where gravity is downward. Horizontal dashed lines above the top light layer and below the bottom light layer suggest the presence of interfaces, although these layers may be unbounded. In case (g), the light upper layer is quite thin or non existent.</figcaption></figure>

is typically explained in terms of a two-layer system of a dense fluid layer over a light fluid layer in a gravitational field [Fig. 3(a)]. However, the concept readily extends to a three-layer system [Fig. 3(b)] where the central layer may be denser or less dense than the other two layers, which in our case, have the same density. If the central layer is more dense, the lower interface is unstable. For a perturbation which induces a thickness variation of the dense layer [Fig. 3(c)], consider the pressures at points 1 and 2, which are at the same height. The pressure at point 1 is higher than at point 2 if the layer of dense
fluid is thicker above point 1 [Fig. 3(c)]. The system is therefore unstable and the perturbation will amplify, assuming no surface tension and negligible diffusion between the phases, which both act to stabilize the system. When this occurs, dense fluid will additionally drain toward point 1 (large white arrows) to further thicken the dense layer at this point. Whether the situation is stable or unstable depends on the potential for drainage in the denser layer which allows thickness variations during the perturbations. Viscous fluids or thin layers may have a limited drainage, and when the lower interface deforms, the upper interface deforms in the same way [Fig. 3(c)]. If the dense layer does not locally thicken, as in Fig. 3(d), the system remains stable because the pressures at points 1 and 2 are equal.

In a flowing bilayer granular system with large particles over small particles, the intermediate denser layer dynamically arises between the two layers by mixing and segregation [Fig. 3(e)]. The layer of mixed particles has a higher volume fraction than monodisperse particles, thereby forming a denser layer. Contrarily to the 3 layers fluid case, the denser layer, and its thickness, are due to a balance between collisional mixing and segregation [Fig. 3(e)]. More segregation, respectively more mixing, induces a thinning, respectively a thickening of the mixed layer. It is likely that drainage plays only a minor role in the process. As segregation and mixing both depends on the vertical location in the flow, the stability of the system depends on the behavior of the dense mixed layer undergoing a vertical perturbation. The descent of the dense layer (above point 1) can be associated to a thickening leading to the unstable case similar to Fig. 3(c); or to a thinning or a constant thickness leading to the stable case similar to Fig. 3(d). Both segregation and mixing intensities depend also on the flow dynamics leading to various possible behaviors of this bidisperse system. Finally, if the mixed layer's thickness decreases where it sinks, the pressure at point 1 can become lower than that at point 2, and the perturbation will be damped.

One particular scenario must be considered. If the upper layer of large particles is very thin, the number of large particles available to form the mixed dense layer may be limited. The packing of the dense layer by mixing is less efficient above as large particles are missing than above point 1 where the upper layer is thicker. When a perturbation of the lower surface of the dense layer occurs, as shown in Fig. 3(g), a pressure difference between points 1 and 2 always results in an instability.

A final comment is in order. For the sake of simplicity, we have presented the mixing of the two species—leading to the formation of a dense layer—and the resulting Rayleigh-Taylor instability as two successive mechanisms. We will see that in most cases, this assumption is justified; however, in certain situations, the instability may begin while diffusion is still ongoing. These cases will be discussed as they are encountered.

### C. Recirculation cells and segregated bands

A striking feature of the isodense granular RTI is that each pair of bands (one band of large particles and one band of small particles) is associated with four recirculation cells at steady-state, as shown in Fig. 4(a), where the velocity field is obtained by averaging the velocities of all particles in each cell of a grid perpendicular to the flow over 10 s of flow at steady state. A pair of large streamwise-oriented recirculation rolls occurs in the region of small particles (blue), and a pair of weaker small rolls occurs in the region of large particles (red). The pairs of rolls in Fig. 4(a) rotate in opposite directions in each region with a net upward flow in between the rolls at the center of each region and a net downward flow between the rolls in the tilted interface region of mixed particles (white) at the interface between the large (red) and small (blue) particle regions, which is associated with an excess density of mixed particles. This recirculation flow is very small, only a few millimeters per second, compared to a mean velocity of 0.6 m/s in the main flow down the incline. The flow in the tilted mixed particle region drives mixed material under the descending plumes, in a similar way to the drainage motion in a dense fluid layer, shown in Fig. 3(c).

The recirculation cells for each species are clearer in Fig. 4(b) where the velocity vectors are plotted for each species after multiplying by $f_i$ to account for the concentration of the species. A schematic of the recirculation rolls is shown in Fig. 4(c). It appears that the recirculation is driven by the density excess at the tilted interface between the rolls. Simultaneous with this tilted buoyant driven flow, size segregation is occurring in the tilted mixed region, evident as a very subtle deviation in the vectors for the two species in Fig. 4(b). The red vectors (large particles) point slightly upward compared to the blue ones (small particles) at the same location at the mixed interface. This indicates that size segregation occurs everywhere in the flow, whenever both species are present. However, the small differences between vectors indicates that the segregation remains small compared to the velocity of the recirculation cells.

The flow induced by the excess density in the mixed zone combined with size segregation that tends to restore less dense pure zones is what drives the rolls and reinforces the segregated bands. For small particles in the mixed particle zone, excess density adds to size segregation to induce a strong recirculation, while for large particles, the two effects oppose each other and the recirculation is weaker. At the free surface, bands of large particles are nearly pure. In contrast, bands of small particles also include some large particles, depending on parameters such as the particle size ratio and flow thickness.

Recirculation cells also occur for axial segregation with isodense size-bidisperse particles in rotating drums [20]. In this case, the recirculation rolls and axial bands evolve more slowly because particles enter the static bed after
<figure><img src="image_5.png" /><figcaption>FIG. 4. Typical example of recirculation cells with $F_l = 0.15$, $d_l/d = 2.5$, $\mu = 0.7$, $H = 32d$, $L = 100d$ and $W = 60d$. The upper thick black vectors in (a) and (b) correspond to 5 mm/s. (a) Transverse velocity field of all particles superimposed with the local concentration of large particles, $f_l$ (see color bar). (b) Velocity field of large (red vectors) and small (blue vectors) particles rescaled by their local concentration, $f_i$. (c) Schematic of the recirculation cells for an isodense granular RTI.</figcaption></figure>

each pass through the flowing layer. Furthermore, the recirculation rolls must be indirectly inferred from the integration of particle trajectories after each pass through the flowing layer in the drum. In fact, it is this issue that motivated the current study of the simpler case of flow down an incline. However, there are significant differences between the recirculation rolls here and those forming axial bands in a rotating drum. Although the small particle rolls are similar in both cases with an upward flow between rolls at the middle of a small particle band, the large particle rolls are in opposite directions in the two cases. In the granular drum case, it seems that segregation at the tilted mixed interface drives the small particle roll downward at the interface, which is the same as that for flow down an incline, and the large particle

roll upward, which is opposite that for flow down an incline. The difference may lie in the relative strengths of the tilted buoyant driven flow of the mixed layer and the segregation-driven upward individual movement of large particles, both occurring at the tilted mixed particle interface. The buoyant driven flow seems to dominate for the incline, whereas the segregation dominates in the rotating drum. However, the differences are subtle and deserve further investigation.

Unlike the two adjacent bands corresponding to four rolls in isodense size-bidisperse systems on inclines (Fig. 4) and for rotating drums [20], a pair of bands corresponds to only two rolls for two particle species having different sizes and densities on inclined planes [18]. Alternating descending plumes of mainly dense large particles and ascending plumes of mainly light small particles form pairs of rolls, similar to Rayleigh-Bénard cells. Contrarily to the iso-dense case, descending plumes of dense large particles between a pair of rolls extend from the upper dense, large particle-rich region downward to the bottom of the region of small less dense particles, and the adjacent upward plume extends from the small particle-rich region upward to the denser region near the surface. Slight mixing occurs during the downward motion that entrains surrounding small particles, and size segregation mainly occurs during the horizontal surface motion expanding above the ascending plumes. Consequently, the surface pattern appears as pure large particle bands always alternating with a mixture of mostly small particles and few large ones.

## IV. CHARACTERISTICS OF THE MIXED LAYER

To focus on the dense mixed layer itself, we utilize a narrow system with spanwise periodic width $W = 25d$ or $W = 30d$ for $d_l/d > 3$. This allows us to examine the profile of the particle volume fraction, $\phi$, corresponding to the monodisperse layers and the mixed layer between them, but prevents the development of the RTI, or at least delays it sufficiently, because the system is too narrow for the perturbation to grow.

### A. Mixed layer excess density

The excess density of the mixed layer can be quantified in terms of the total particle volume fraction, which depends on several factors. Consider first the particle size ratio, $d_l/d$. Figure 5 presents the volume fraction profiles from the surface of the incline at $z = 0$ to the free surface of the flow at $z/d \approx 32$ for $d_l/d = 1.5, 2, 2.5$, and 3. The solid curves indicate the steady-state total particle volume fraction profiles, $\phi$, averaged over $x$ and $y$ for a period of 10 s, while the dotted and dashed curves show the small and large particle volume fractions, $\phi_s$ and $\phi_l$. The two layers are evident with a very low
<figure><img src="image_6.png" /><figcaption>FIG. 5. Particle volume fraction, $\phi$ (solid curves), small particle volume fraction, $\phi_s$ (dotted curves), and large particle volume fraction, $\phi_l$ (dashed curves), versus height for various particle size ratio $d_l/d$, with $H = 32d$, $F_l = 0.3$, $\mu = 0.7$, and $30 s < t < 40 s$.</figcaption></figure>

large particle volume fraction for $0 \le z/d \lesssim 15$, rising to $\phi_l \approx 0.6$ for $27 \lesssim z/d \lesssim 32$, where the upper layer is nearly pure large particles. The small particle volume fraction starts at $\phi_s \approx 0.56$ in the lower portion of the flow for $0 \le z/d \lesssim 15$ and decreases to zero in the large particle layer for $z/d \ge 27$. The oscillations in $\phi$ and $\phi_s$ at small $z/d$ is related to particle ordering near the wall.

The dense mixed layer of both small and large particles occurs in the region where the $\phi_s$ and $\phi_l$ curves are both non-zero ($10d \lesssim z \lesssim 27d$) with a maximum at $z/d \approx 25$, evident as a bump in the total volume fraction, $\phi$. As the size ratio increases, the smaller particles fit more effectively into the voids between the larger ones and the excess density in the mixed layer becomes more pronounced. The oscillations in $\phi$ at $d_l/d = 3$ appear due to the preferential layering of large particles during flow. Also note that the more compact packing at large $d_l/d$ (magenta curve) reduces the overall height of the free surface, located around $z/d \approx 32$. We further note that the maximum in $\phi$ occurs where locally $f_l \approx 0.7$, which corresponds to the value of $f_l$ for the maximum static packing of size-bidisperse particles for $d/d_l = 2$ and 3, although the maximum packing in a static bed after tapping is slightly higher [30].

Figure 6 presents a magnified view of the excess density bumps in the volume fraction profiles associated with the mixed particle layer. The increase in $\phi$ in the mixed particle layer above that for the primarily small particle layer for $z/d \lesssim 20$ is relatively small. To obtain a consistent measure of the maximum volume fraction, $\phi_{max}$, we fit a vertically shifted Gaussian function $\phi = a + b \exp(-(z-m)^2/\sigma^2)$ to the bumps of the profile (dashed curves in Fig. 6). The maximum of the bump is given by $\phi_{max} = a + b$. A fourth-degree polynomial was also tested and yielded similar results. Figure 7 presents the maximum volume packing as a function of the particle size ratio $d_l/d$ for two friction coefficients,

<figure><img src="image_7.png" /><figcaption>FIG. 6. Particle volume fraction $\phi$ of the mixed layer versus height (measured in small particles diameter d) for various particle size ratio $d_l/d$, with $H = 32d$, $F_l = 0.3$ and $\mu = 0.7$. Solid curves depict $\phi$, and dashed curves correspond to Gaussian fits.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Maximal particle volume fraction $\phi_{max}$ given by Gaussian fits versus particle size ratio $d_l/d$ for $H = 32d$, $F_l = 0.3$ and two friction coefficients, $\mu = 0.3$ (blue squares) and $\mu = 0.7$ (red circles).</figcaption></figure>

$\mu = 0.3$ and $0.7$. The maximum particle volume fraction increases nearly linearly with $d_l/d$, with a slight decrease in slope at large $d_l/d$. A lower friction coefficient results in a higher volume fraction, as would be expected.

## B. Mixed layer thickness

The width of the Gaussian fits in Fig. 6 suggest that particle size ratio has little effect on the mixed layer thickness, even though the value of $\phi_{max}$ increases with $d_l/d$. On the other hand, the concentrations of the two particle species have a much larger effect on the mixed layer thickness, as shown in Fig. 8 for $0.1 \le F_l \le 0.5$. For $F_l = 0.1$, the mixed particle layer is localized near the free surface. Increasing $F_l$ broadens the bump in $\phi$ associated with the mixed particle layer, as more large
<figure><img src="image_9.png" /><figcaption>FIG. 8. Particle volume fraction $\phi$ versus height for various volume fraction of large particles $F_l$, with $H = 32d$, $d_l/d = 2.5$ and $\mu = 0.7$. Dashed curves correspond to a Gaussian fit of $\phi$. The inset illustrates how the thickness of the mixed particle layer is based on the standard deviation of $\phi$ for the case $F_l = 0.2$.</figcaption></figure>

particles are available for the mixed particle layer, and shifts the bump toward lower $z$, as the layer of large particles at the surface thickens.

As $F_l$ increases, $\phi_{max}$ for the Gaussian fits remains nearly unchanged (not accounting for individual peaks associated with preferential layer of the large particles), but the mixed layer thickens and shifts to a lower $z/d$. Note that for the lowest values of $F_l$, the Gaussian fits extend above the free surface. Since the standard deviation of a such a fit would overestimate the thickness of the mixed layer for cases like this, the thickness of the bump is directly computed from the standard deviation of the hatched portion of the volume fraction profiles lying above a baseline defined by a Gaussian fit (dotted line), as illustrated in the inset of Fig. 8. Figure 9 shows how the thickness of the density excess layer, as quantified by the standard deviation, increases with $F_l$, before reaching a plateau. The thickness of the mixed layer also increases with increasing flow depth $H$. Both of these results can be interpreted in terms of the availability of large particles to form the mixed layer. Increasing the number of large particles in the system by increasing the total volume of particles (by increasing $H$) or by increasing the concentration of large particles, $F_l$, both result in a thicker mixed particle layer.

### C. Bilayer versus homogeneous initial conditions

We now return to the two initial conditions mentioned earlier: an initially segregated two-layer system of large particles above small particles and homogeneously mixed small and large particles. Figure 10 highlights the temporal evolution of the initially two-layer (red) system and the initially homogeneously mixed (blue) system in terms of the vertical segregation index ($I_{seg}^v$) and the span-

<figure><img src="image_10.png" /><figcaption>FIG. 9. Thickness of the excess density layer, measured as the standard deviation, $\sigma$, of the Gaussian fit to $\phi$ at the excess density bump in $\phi$, versus large particle concentration $F_l$ for different values of flow depth $H$ for $d_l/d = 2.5$ and $\mu = 0.7$.</figcaption></figure>

wise segregation index ($I_{seg}^s$). The vertical segregation index for the bilayer case (red dashed curve) starts from $I_{seg}^v = 1$ for perfect vertical segregation but decreases between $t = 0$ s and $10$ s due to the diffusion between the layers of large and small particles and the formation of an intermediate mixed layer (see vertical cross-section at $t = 10$ s). The decrease in $I_{seg}^v$ slows around $t = 10$ s but then decreases further around $t \approx 40$ s, which is associated with a rise in the spanwise segregation index (red solid curve), indicating the onset of spanwise segregation. Bands of small particles reach the surface, while bands of large particles become thicker than the initial upper layer (see vertical cross-section at $t = 40$ s). Spanwise segregation continues to evolve, intensifying the band pattern, and a quasi-steady state is reached around $t \approx 70$ s. We call this quasi-steady because over extended times, bands may reorganize, either through the merging of two bands or the disappearance of an initially narrow band.

It is also worth noting that the final state does not resemble that for a classical Rayleigh-Taylor instability for dense and light fluid layers, where after the formation of a wavy interface, protrusions into the opposing phase result in downward heavy fluid plumes and upward light fluid plumes that transport fluid until a stable state is reached in which the two layers are inverted from the initial condition. Instead, here a system of longitudinal bands emerges where bands of large particles protrude into the lower layer but whose descent is frozen, which is related to the recirculation cells driven by the excess density of the mixed zone reinforced by size segregation at the tilted interface between the large and small particle bands, as described previously (Fig. 4).

The initially homogeneous case starts quite differently but exhibits similarities to the bilayer case after some time. The vertical segregation index (dashed blue curve) starts at $I_{seg}^v = 0$, corresponding to a homogeneously mixed system, and increases to $I_{seg}^v \approx 0.7$ due to granular segregation, which results in a layer of nearly pure
<figure><img src="image_11.png" /><figcaption>FIG. 10. Time evolution of the vertical segregation index ($I_{seg}^v$, dashed curves) and spanwise segregation index ($I_{seg}^s$, solid curves) for initially bilayer (red) and homogeneous (blue) systems. Simulation parameters are $H = 24d$, $F_l = 0.2$, $d_l/d = 3$, and $\theta = 23^\circ$. Data points on the curves correspond to successive images of vertical cross-sections, perpendicular to the flow, for the cases of initially (b) bilayer and (c) homogeneous systems.</figcaption></figure>

large particles at the free surface above a layer composed mainly of small particles with an intermediate mixed layer. $I_{seg}^v$ decreases as the RTI begins, associated with collisional diffusion to form a mixed layer and then a rise of the spanwise segregation index $I_{seg}^s$ (blue solid curve). By $t = 40$ s, the pattern becomes similar to the initially bilayer system at shorter times. As vertical segregation decreases, spanwise segregation develops, but at a slower rate compared to the initially bilayer case. Nevertheless, the timescale for segregation and band formation in the initially homogeneous case (approximately 60 s) is not much longer than that for the initially bilayer case (approximately 40 s). The steady states reached in the two cases are very similar (see vertical cross-sections at $t = 80$ s), although the segregation indices differ slightly because they depend on both the number of bands and on their vertical and horizontal extents.

Since the timescales for the initially bilayer case are shorter than for the initially homogeneously mixed case,

we consider only the initially bilayer case in the remainder of the paper. Also note that for the homogeneous initial condition, simulations with periodic boundary conditions are equivalent to experiments with an infinitely long incline. During flow, the global concentration of each species remains constant while the streamwise flux of large particles increases due to size segregation as large particles rise to regions of faster flow. Alternately, for an experiment with equal concentrations of each species introduced at the upstream end of an incline, the streamwise flux of each species remains constant during segregation, but the large particles move faster in the upper portion of the flow, and the faster-moving large particle layer becomes thinner compared to the slower small particle layer. In the case of an initially bi-layer system, simulations and experiments on a long enough incline are essentially equivalent.

To put these results in context, it is helpful to convert these timescales into the distance traveled by particles near the free surface of the flow. In 60 s, for $H = 24d$, $F_l = 0.2$ and $\theta = 23^\circ$, the distance covered is $\sim 20000d$, which is $\sim 40$ m for the particle size used in the simulation. As a result, it is difficult to observe spanwise segregation on an inclined plane for size-bidisperse particles at the laboratory scale, as even with $d = 200$ µm particles, an inclined plane of approximately 4 m in length would be required. However, spanwise, or axial, segregation can be observed at laboratory scale in a rotating drum. For example, axial segregation is achieved in approximately 300 s in a 0.18 m diameter cylinder rotating at 15 rpm that is half-filled with equal volumes of 2 mm and 6 mm particles [20]. Assuming that particles pass through the 0.18 m long flowing layer approximately twice per drum rotation, particles at the flowing layer surface travel approximately 27 m, or 13500d, which is comparable to the inclined plane case investigated here.

## V. PARAMETRIC ANALYSIS FOR BAND FORMATION

While the influence of parameters including the concentration of large particles and the particle size ratio were explored in the context of the appearance of a dense mixed layer necessary for the granular Rayleigh-Taylor instability in previous sections, here we consider the development of the instability itself, subsequent band formation, and their dependence on various parameters of the flow.

### A. Flow thickness and band formation

To consider band formation for various flowing layer thicknesses, the simulation domain width is increased to $W = 400d$ for more accurate measurements of the segregated band wavelength, and its length is reduced to $L = 60d$ to maintain a reasonable computational time,
noting that the segregation pattern is invariant in the flow direction. Figure 11 shows space-time diagrams of the local large particle concentration averaged in the *x* and *z* directions, $f_l(y, t)$, for various flow thicknesses, $H$, for an initial bilayer condition at global large particle concentration $F_l = 0.2$. Gray regions correspond to small particle bands, while red regions correspond to large particle bands, which have relatively low large particle concentration when averaged through the depth of the flowing layer. In all cases, bands become visible within the first 20 s. Typically, during the initial band formation, there are more bands and their spanwise wavelengths are shorter compared to the quasi-steady regime after about 50 s, when some band merging has occurred, which, of course, increases the wavelength.

At steady-state, the bands are well-differentiated, so the wavelength is calculated by counting the number of bands. However, for the initial band formation, counting bands is more subjective. Hence, we use the objective measure of the spanwise segregation index reaching the value of $I_{\text{seg}}^s = 0.1$. As shown later (see Fig. 14), this value for $I_{\text{seg}}^s$ corresponds to the onset of band formation very early in the spanwise segregation process. The spanwise concentration profile at this time is extracted, and a Fourier transform is used to determine the dominant wavelength. When the same procedure is applied after bands merge and steady-state is achieved, the results are consistent with simply counting bands to determine the wavelength.

<figure><img src="image_12.png" /><figcaption>FIG. 11. Space-time diagrams of $f_l(y, t)$ showing the spanwise segregation of a size bidisperse inclined plane flows of thicknesses $12d \le H \le 50d$, with $d_l/d = 3$, $F_l = 0.2$, $\mu = 0.3$, $W = 400d$, and $L = 60d$. Plots are 500 s long for $H < 32$ and 250 s otherwise.</figcaption></figure>

The dependence of the band wavelength on the flow thickness $H$ during both the initial band formation and at quasi-steady-state are shown in Fig. 12. Error bars reflect the uncertainty in basing the wavelength on counting an

integer number of bands in the finite width periodic domain. For instance, if 8 pairs of large and small particle bands appear in the $W = 400d$ wide domain ($\lambda = 50d$), the error bar represents the range of 7.5 to 8.5 pairs of bands ($\lambda$ of 53d and 47d, respectively) that could occur if the domain width was not artificially constrained.

For initial band formation, $\lambda$ increases nearly linearly with $H$ with a relation of $\lambda \sim 2H$, similar to the granular RTI with density difference [18]. This is consistent with approximately circular rolls of small particles that have a diameter of about $H$, probably due to the constraints of the bottom and free surfaces.

Quasi-steady results are quite different for small flow thicknesses ($H \lesssim 16d$). The quasi-steady wavelength at $H = 12d$ is quite large and it quickly decreases with increasing $H$. To understand this result, consider the species concentration perpendicular to the flow direction shown in Fig. 13 for a 100d-wide portion of the domain containing one or more large particle bands. For $H = 12d$, the large particles form a thin layer at the free surface and do not penetrate far into the flow depth, unlike for larger $H$ where large particles penetrate deeper into the flow. As a consequence, spanwise segregation converges to a pattern of bands with a large wavelength for $H = 12d$. For $H \ge 16d$, the depth of the large particle bands are less constrained resulting in an approximately linear correlation of $\lambda \sim 3H$, consistent with band merging that occurs for the quasi-steady regime for particles with a density difference [18]. Of course, as noted earlier, there is the possibility of further band merging after much longer times. Here, the simulations are stopped after twice the time from the last band merge.

Finally, it is important to note that, in both size-bidisperse and density-varying granular flows on inclined planes, the wavelength is comparable to the flow depth $H$, resulting in a band spacing wavelength $\lambda$ of approximately $2H$ to $3H$. This result is different from that in rotating drums, where the wavelength is on the order of the drum diameter [34, 35]. For a 50% filled drum, the flowing zone thickness is typically a quarter of the radius, leading to a wavelength approximately eight times the flowing zone thickness. This difference is likely due to the more complex geometry and flow in a rotating drum.

We also consider the temporal evolution of the band formation in terms of the spanwise segregation index in Fig. 14. The evolution of the segregation is quite different for shallow and deep flowing layers. For $H \le 24$, $I_{\text{seg}}^s$ increases monotonically to its steady-state value, which is relatively large, demonstrating strong spanwise segregation. Although the initial growth is also monotonic for $H \ge 32d$, once $I_{\text{seg}}^s$ reaches a maximum value, it decreases before converging to a quasi-steady value that decreases with increasing $H$. This time evolution of $I_{\text{seg}}^s$ is also evident in the space-time plots of Fig. 11, where for large $H$ the light colored bands corresponding to small particles are more pure for $20 \lesssim t \lesssim 40$, indicating stronger spanwise segregation, than at larger $t$.

Figure 13 is helpful to explain the lower quasi-steady
<figure><img src="image_13.png" /></figure>

FIG. 12. Band wavelengths $\lambda$ at initial destabilization (red squares) and at steady state (blue circles) as a function of the flow thickness $H$ for $d_l/d = 3$, $F_l = 0.2$, $\mu = 0.3$, $W = 400d$, and $L = 60d$. Quasi-steady states are obtained for times longer than those presented in Fig. 11, as discussed in the text. Error bars reflect the uncertainty in the wavelength based on counting an integer number of bands in a finite-width periodic domain.

spanwise segregation at large $H$. As evident in the figure, the small particle regions are less pure as $H$ increases. At large $H$, large particles from the mixed zone of the downward plumes that form after the original RTI are entrained and contaminate the small particle regions, leading to a decrease in $I_{\text{seg}}^s$ at steady state. This effect is amplified with increasing $H$ because diffusion increases due to stronger shear in a deeper flowing layer. For small $H$, the shear rate is smaller, so diffusion is reduced. As a result, the large particles segregate efficiently rather than being entrained into the bands of small particles, thereby maintaining a high level of $I_{\text{seg}}^s$.

## B. Mixture composition

The effect of the mixture composition on quasi-steady bands for an initially bi-layer flow is shown qualitatively in Fig. 15 for $F_l = 0.1$ to 0.4. Spanwise segregation is clearly visible for $F_l \le 0.2$, but only one band of small particles has emerged at the surface for $F_l = 0.3$, and no bands appear for $F_l = 0.4$ or $F_l = 0.5$ (not shown). The lack of band formation is in contrast to the formation of bands in a rotating drum for $F_l = 0.5$ [20]. However, this difference can be explained by the flow geometry for segregation in a rotating drum. After the initial radial segregation in a drum filled with a $F_l = 0.5$ mixture, the large particles are in the upper 30% of the flowing layer thickness and the small particles are in the lower 70%, assuming a linear velocity profile. This corresponds to $F_l = 0.3$ for the bilayer initial condition in the incline system studied here, where bands form, as shown in Fig. 15.

To better understand the conditions for the onset of instability and spanwise band formation, Fig. 16 presents

<figure><img src="image_14.png" /></figure>

FIG. 13. Concentration of large particles $f_l$ in a plane perpendicular to the flow for $H = 12d$, $24d$ and $40d$ with $d_l/d = 3$, $F_l = 0.2$, $\mu = 0.3$, $W = 400d$, and $L = 60d$. Only 100d-wide portions of the domain are shown. Snapshots are taken at $t = 200$ s, when the systems are quasi-static.

<figure><img src="image_15.png" /></figure>

FIG. 14. Time evolution of the spanwise segregation index $I_{\text{seg}}^s$ for different flow thicknesses, $H$, with $F_l = 0.2$, $d_l/d = 3$, $\mu = 0.3$, $W = 400d$, and $L = 60d$.

the temporal evolution of the spanwise segregation index, $I_{\text{seg}}^s$, for a wide range of $F_l$. In these simulations, the width has been reduced to $W = 60d$, and a sinusoidal perturbation with wavelength $\lambda = W = 60d$ and amplitude of $B_0 = 0.1H$ is introduced at the interface between large and small particles. This simulation protocol allows us to monitor the growth of the instability
<figure><img src="image_16.png" /><figcaption>FIG. 15. Impact of the large particle concentration ($F_l$ = 0.1 to 0.4) on band formation in a $W = L = 200d$ periodic domain at $t = 160$ s for $H = 24d$, $d_l/d = 3$, and $\mu = 0.3$. The black arrow indicates the flow direction.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 16. Time evolution of the spanwise segregation index $I_{seg}^s$ for a range of large particle concentrations $F_l$ with $H = 32d$, $d_l/d = 2.5$, $\mu = 0.7$, $W = 60d$, and $L = 100d$. An initial sinusoidal perturbation of amplitude $B_0 = 0.1H$ is used, which results in an initial value of $I_{seg}^s$ that depends on $F_l$ based on Eq. 3.</figcaption></figure>

or, alternatively, the decay of the initial perturbation. In all cases, the initial value for $I_{seg}^s$ is based on Eq. 3, which is near $B_0 = 0.1$ for large $F_l$ but increases substantially for $F \le 0.2$.

From Fig. 16, it is evident that spanwise segregation due to the instability grows faster and to a larger degree for small values of $F_l$. On the other hand, for $F_l > 0.3$, the spanwise segregation decreases, indicating that the initial perturbation is damped, resulting in no spanwise bands forming. For the case of $F_l = 0.1$, the initial strong growth followed by decay is related to downward plumes that form after the original RTI carrying large particles downward that eventually contaminate small particle regions, similar to the situation for large $H$ in Fig. 14.

To understand the mechanisms involved, again con-

<figure><img src="image_18.png" /><figcaption>FIG. 17. Local large particle concentration, $f_l$ (red-blue colorbar), and particle volume fraction, $\phi$ (yellow-black colorbar), for global large particle concentrations $F_l = 0.15$ (a-b) and $F_l = 0.5$ (c-d) for $H = 32d$, $L = 100d$, and $W = 60d$ with $\mu = 0.7$ and $d_l/d = 2.5$. The species are initially positioned in two layers with a sinusoidal perturbation of the interface having an amplitude $B_0 = 0.1H$, and the fields are measured after 10 s.</figcaption></figure>

sider $W = 60d$ domain width with an initial sinusoidal perturbation of the interface having an amplitude of $B_0 = 0.1H$. Figure 17 shows the local concentration of large particles and the volume fraction in a vertical cross-section perpendicular to the flow for two different concentrations of large particles after 10 seconds of flow, allowing sufficient time to establish the unstable mixed layer.

In the case of a low concentration of large particles [Fig. 17(a, b)], $F_l = 0.15$, the thickness of the excess density layer (black) is quite small, about $1d_l$, when it is near the free surface in the region $5d \lesssim x \lesssim 25d$, whereas it thickens in the region where the interface sinks, $30d \lesssim x \lesssim 60d$. Consequently, large particles in
the thicker part of the surface layer enter the downward plume, amplifying the instability, as shown schematically in Fig. 3(g). Conversely, in the case of a larger concentration of large particles, $F_l = 0.5$, the thickness of the excess density region is more uniform across the entire width of the flow [Fig. 17(c, d)]. The spanwise variation density is weak and insufficient to amplify the perturbation, which is eventually damped out, corresponding to the schematic in Fig. 3(f).

The results in Fig. 17 can also be understood in terms of Fig. 8, considering the profiles of $\phi$, not the Gaussian fits. When the concentration of large particles is small, for example, $F_l = 0.2$ (blue curve), a perturbation in the mixed layer will create spanwise regions with higher and lower concentrations of large particles. For lower $F_l$ cases, the excess density decreases sharply near the free surface, corresponding to the thinner mixed layer. Spanwise variations in $F_l$ create spanwise mixed layer thickness and unstable density variations like that in Fig. 17(a, b), where the thin dense layer corresponds to $F_l = 0.1$ and the thick dense layer corresponds to $F_l = 0.3$. This is equivalent to the schematic in Fig. 3(g), which is unstable. Applying the same reasoning to the case of $F_l = 0.4$ (violet curve in Fig. 8), the two neighboring curves, black and red, have similar excess densities and thickness, making it less likely that a spanwise variation of pressure occurs, so the likelihood of an instability is reduced, similar to Fig. 17(c, d) and shown schematically in Fig. 3(f).

<figure><img src="image_19.png" /><figcaption>FIG. 18. Phase diagram of the stable and unstable regions for the RTI dependence on $H$ and $F_l$. The magnitude of the growth rate is indicated by the color bar. With $\mu = 0.3$, $d_l/d = 3$, $W = 200d$, and $L = 100d$.</figcaption></figure>

Figure 18 shows a phase diagram of the stable and unstable regions as a function of the flow thickness $H$ and the concentration of large particles $F_l$ for $d_l/d = 3$, $W = 200d$, and $L = 100d$ (to assure adequate length to assure isotropy for the roughness of the incline). As the depth of the flow increases, the instability can develop for larger $F_l$. The color code indicates the initial growth rate of the RTI. The growth rate $\gamma$ is readily obtained from a semi-log plot of $I_{seg}^s$ versus $t$ for large values of $\gamma$, but for the low growth rate values (yellow), diffusion and other

factors make estimating the growth rate more difficult, so the values are only indicative. Further note that for large values of $H$, where $I_{seg}^s$ grows quickly at first and then decreases to its steady-state value (Fig. 14), $\gamma$ reflects only the initial growth associated with the RTI. With these caveats, it is evident that the growth rate increases with increasing $H$ and to a lesser extent with decreasing $F_l$, as discussed in the context of Figs. 14 and 16.

### C. Particle size ratio

As shown in Fig. 6 for $H = 32d$, as the particle size ratio, $d_l/d$, increases, the particle volume fraction increases, leading to a higher density in the mixed particle layer. In the fluid RTI, the growth rate increases with the density ratio, so we would expect that the particle size ratio would similarly affect the granular RTI instability and subsequent formation of spanwise bands. Figure 19 demonstrates this result qualitatively for the granular RTI as spanwise segregation and subsequent bands depend on the particle size ratio, with no bands forming for $d_l/d = 2$ and narrower large particle bands for the largest size ratio $d_l/d = 3$. The differences likely result from a complex interplay between the excess density in the mixed particle layer that increases with the particle size ratio, and size segregation, which tends to plateau for $d_l/d > 2$ [33].

The time evolution of the spanwise segregation index in Fig. 20 shows how the growth rate increases with the particle size ratio $d_l/d$ with a maximum in the range $2.5 \le d_l/d \le 3$ for $H = 24d$. Figure 20 confirms that for $d_l/d \le 2$ no instability occurs and bands do not form. Note that even though the density excess due to mixing increases nearly linearly with the particle size ratio (see Fig. 7), there is a clear threshold for the instability and

<figure><img src="image_20.png" /><figcaption>FIG. 19. Impact of the size ratio $d_l/d$ on band formation in a $W = L = 200d$ periodic domain at $t = 160$ s for $F_l = 0.2$, $H = 24d$, and $\mu = 0.3$. The black arrow indicates the flow direction.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIG. 20. Time evolution of the spanwise segregation index $I_{seg}^s$ for particle size ratios ranging from $d_l/d = 1.75$ to 3.25 with $F_l = 0.2$, $H = 24d$, $\mu = 0.3$, and $W = L = 200d$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>FIG. 21. Growth rate $\gamma$ as a function of particle size ratio $d_l/d$ for $F_l = 0.2$, $H = 32d$, $\mu = 0.3$, $W = 60d$, and $L = 100d$ with an initial perturbation of 0.1H to favor the instability.</figcaption></figure>

band formation.

To more clearly show the growth rate, $\gamma$, we consider $W = 60d$ with an initial sinusoidal perturbation. As shown in Fig. 21, $\gamma$ increases linearly up to a size ratio $d_l/d \simeq 3$, consistent with the nearly linear increase in $\phi_{max}$ in Fig. 7, but then it saturates for $d_l/d > 3$. For $H = 32d$ used here, the instability appears at a lower value of $d_l/d$ than for $H = 24d$ in Fig. 20 showing how multiple factors play a role in the granular RTI.

Figure 22 shows a phase diagram of the stable and unstable regions for the RTI as a function of the particle size ratio $d_l/d$ and the concentration of large particles $F_l$. The growth rate of the instability, indicated by color, is largest for $0.2 \le F_l \le 0.3$ and $d_l/d \ge 3$. It decreases at higher $F_l$ until it disappears for $F_l \ge 0.4$. The instability also disappears at low values of $d_l/d$ as would be expected because of the reduced tendency for size segregation. Based on Fig. 18, which is for somewhat different simulation conditions, we would expect the region that is unstable to expand slightly for larger $H$ and contract

<figure><img src="image_23.png" /><figcaption>FIG. 22. Phase diagram of the stable and unstable regions for the RTI dependence on $d_l/d$ and $F_l$ with $\mu = 0.3$, $H = 32d$, $W = 60d$, $L = 100d$, and an initial perturbation of 0.1H to favor the instability. The magnitude of the growth rate $\gamma$ is indicated by the color bar.</figcaption></figure>

slightly for smaller $H$. We also note that the results here do not match perfectly with Fig. 16 where the maximum value of $\gamma$ occurs at $F_l = 0.10$ for $d_l/d = 2.5$, but this a result of different parameter values.

D. Friction coefficient

<figure><img src="image_24.png" /><figcaption>FIG. 23. Time evolution of the spanwise segregation index for different friction coefficients $\mu = 0.2$ to 1.5 for $F_l = 0.2$, $H = 24d$, $d_l/d = 3$, $W = 200d$, and $L = 60d$.</figcaption></figure>

Finally, we consider the impact of the friction coefficient because of the dependence of the RTI on the density of the mixed layer, which depends on collisional diffusion, which, in turn, is affected by friction. Figure 23 shows the time evolution of $I_{seg}^s$ for various sliding friction coefficients, $\mu$, with $F_l = 0.2$. Increasing $\mu$ from 0.2 to 1.5 reduces the growth rate from $\gamma = 0.16$ to 0.075, indicating the significant role of $\mu$ on the instability. As shown in Fig. 7 in terms of the maximum volume frac-
<figure><img src="image_25.png" /><figcaption>FIG. 24. Time evolution of the spanwise segregation index for different friction coefficients $\mu = 0.3$ to $1.5$ at $F_l = 0.3$, with $H = 24d$, $d_l/d = 3$, $W = 200d$, and $L = 60d$.</figcaption></figure>

tion, $\phi_{max}$, increasing $\mu$ reduces the excess density of the mixed layer. Nevertheless, the instability develops regardless of $\mu$ for $0.2 \le \mu \le 1.5$, although the degree of spanwise segregation at the steady state is reduced for $\mu = 1.5$.

Increasing the large particle concentration to $F_l = 0.3$, which corresponds to the system being closer to the stable-unstable transition (see Fig. 18), suppresses the instability for the highest friction considered, $\mu = 1.5$. As shown in Fig. 24, the initial exponential growth is clear for small $\mu$, but the instability develops more slowly or not at all as $\mu$ increases. This comes about from the combined effects of $\mu$ on reducing $\phi_{max}$ and $F_l$ being closer to the stability limit (see Fig. 18), thereby reducing the likelihood of the mixed layer becoming unstable.

## VI. CONCLUSIONS

We demonstrate in this article that a granular Rayleigh-Taylor instability (RTI) in the flow of dry non-cohesive granular material of size-bidisperse particles down an incline leads to spanwise surface bands of large and small particles, even when both particle types have the same density. The instability is a consequence of the coexistence of size segregation and collisional mixing due to the shear, which induces the formation of a three-layer system. The intermediate layer is a mixture of particles of both sizes, making it more densely packed than the upper monodisperse layer of large particles and the lower monodisperse layer of small particles. This results in a granular RTI between this denser interfacial mixed layer and the less-dense monodisperse lower layer. The granular RTI may develop or not depending on the relative thicknesses of the three layers. When the upper layer of segregated large particles is thin enough, RTI develops. On the other hand, high concentrations of large particles leads to a stable three-layer system. The development of the instability is initially exponential with a growth rate

that increases due to increasing particle volume fraction $\phi$ for $d_l/d \le 3$.

Compared to the fluid RTI where the heavy and light layers of fluid either invert to a stable situation with the light fluid above the heavy fluid or mix due to the upward and downward plumes of the light and heavy fluids, the granular RTI gives rise to recirculation rolls with descending plumes of mixed particles topped with large particles and ascending plumes of small particles such that the regions of dense and light materials in the granular RTI reorganize and segregate into streamwise bands of small and large particles along the incline, even when starting from a fully mixed condition. A band of large particles visible at the free surface corresponds to a region of large particles accumulated at the top of a descending plume of a mixture of large and small particles denser than the pure regions. By contrast, a band composed mainly of small particles visible at the free surface has roughly the same composition through the entire depth of the flow. The wavelength of bands resulting from the granular RTI scale with the thickness of the flow down the incline as $\lambda \approx 2H$, similar to that for the RTI in fluids, as we previously demonstrated [19] based on applying the RTI dispersion relation in fluids [36] to a granular system.

The RTI and subsequent formation of streamwise-oriented bands of large and small particle bands in isodensity inclined plane flow described in this article have several similarities to that for axial band formation in a rotating drum [20]. In both cases, alternating streamwise-oriented bands of large and small particles appear at the surface of the flow, while the bottom layer (or the core in the case of a rotating drum) is mainly composed of small particles. For flow down an incline, there are two counter-rotating rolls associated with each band with an upward plume of mainly small particles at the center of a small particle band and a downward plume of large particles at the center of the large particle band. These rolls are driven first by the downward motion of the mixed region at the tilted mixed interface, and second by size segregation in the mixed region. In this mixed particle tilted interface region, large particles move slightly upward at the interface, driving small rolls of large particles so that descending large particle plumes do not extend to the bottom wall. At the same time, small particles at the mixed particle interface segregate downward out of the mixed layer, driving large rolls of small particles that extend from the bottom to the free surface of the inclined plane, resulting in the formation of bands of small particles. Thus, for each pair of bands, four rolls are observed.

Although the RTI mechanism initiates the process leading to recirculation rolls and band formation in the three cases of isodense size-bidisperse flow down an incline, isodense size-bisperse flow in a tumbler, and bidisperse flow down an incline of particles having different densities, the number and nature of the rolls for these cases are different. For the isodense flows down an incline,
cline studied here and in a rotating drum [20], there are four rolls for each pair of bands. In the case of RTI with particles having different densities, there are two rolls for each pair of bands and plumes extending through the entire flow thickness, looking similar to RB convection [18]. The differences seem to lie in the relative impacts of segregation, collisional diffusion, and flow, as discussed in detail in Section III C.

Perhaps one of the most interesting results here is that even a homogeneous mixture of particles with two sizes and the same density spontaneously tends towards the RTI as a consequence of size segregation initially forming a layer of large particles above a layer of small particles. The subsequent balance between mixing and size segregation creates a gravitationally unstable mixed layer between these two layers, which is a prerequisite for the

RTI to develop. We call this a self-induced RTI [18].

It is evident from the results presented here and in our previous papers [18, 20] that the granular nature of the flowing material and the existence of size and density segregation results in a rich range of new phenomena that can occur related to the granular RTI. These phenomena reflect mechanisms related to mixing, segregation, and packing that are not usually associated with the Rayleigh-Taylor instability, and, hence, suggest many opportunities for further research.

## ACKNOWLEDGMENTS

Centre de Calcul Intensif d'Aix-Marseille University is acknowledged for granting access to its high performance computing resources.

[1] Lord Rayleigh (J. W. Strutt), Proc. R. Soc. London **14**, 170-177 (1883).

[2] G. I. Taylor, Proc. R. Soc. London Ser. A **201**(1065), 192-196 (1950).

[3] S. Chandrasekhar, *Hydrodynamic and Hydromagnetic Stability*, (Oxford University Press, London, 1961).

[4] F. Charru, *Hydrodynamic Instabilities*, (Cambridge University Press, Cambridge, 2011).

[5] D. H. Sharp, Physica D **12** 3-18 (1984).

[6] A. W. Cook and D. Youngs, Scholarpedia **4**(2), 6092 (2009).

[7] A. Banerjee, J. Fluids Eng. **142**, 120801 (2020).

[8] A. Banerjee and M. J. Andrews, Int. J. Heat Mass Transfer **52**, 3906-3917 (2009).

[9] K. Guo, R. Chen, Y. Li, W. Tian, G. Su, and S. Qui, Prog. Nucl. Energy **109** 130-144 (2018).

[10] A. W. Cook and P. E. Dimotakis, J. Fluid Mech. **443**, 69-99 (2001).

[11] S. I. Abarzhi, A. Gorobets, and K. R. Sreenivasan, Phys. Fluids **17**, 081705 (2005).

[12] H. Liang, X. Hu, X. Huang, and J. Xu, Phys. Fluids **31** 112104 (2019).

[13] P. Suchandra and D. Ranjan, J. Fluid Mech. **974**, A35 (2023).

[14] J. R. Lister and R. C. Kerr, J. Fluid Mech. **202** 577-594 (1989).

[15] G. K. Batchelor and J. M. Nitsche, J. Fluid Mech. **227** 37 (1991).

[16] S. A. Piriz, A. R. Piriz, N. A. Tahir, and F. Cobos, Phys. Fluids **37** 022130 (2025).

[17] J. W. Jacobs and S. B. Dalziel, J. Fluid Mech. **542**, 251-279 (2005).

[18] U. D'Ortona and N. Thomas, Phys. Rev. Lett. **124**, 178001 (2020).

[19] U. D'Ortona, D. Martinand, and N. Thomas, Rayleigh-

Taylor instability in two-layer granular flows, Phys. Rev. E, accepted for publication (2025).

[20] U. D'Ortona, R. M. Lueptow, and N. Thomas Phys. Rev. Research **6**, L032038 (2024).

[21] P. A. Cundall and O. D. L. Strack, Geotechnique **29**, 47 (1979).

[22] C. Kloss, C. Goniva, A. Hager, S. Amberger, and S. Pirker, Prog. Comput. Fluid Dyn. **12**, 140 (2012).

[23] LIGGGHTS documentation website: https://www.cfdem.com/media/DEM/docu/Manual.html, [Accessed: 25-August-2024].

[24] T. G. Drake and R. L. Shreve, J. Rheol. **30**, 981 (1986).

[25] S. F. Foerster, M. Y. Louge, H. Chang, and K. Allia, Phys. Fluids **6**, 1108 (1994).

[26] J. Schäfer, S. Dippel, and D. E. Wolf, J. Phys. I (France) **6**, 5 (1996).

[27] J. W. Vallance and S. B. Savage, Particle segregation in granular flows down chutes, in IUTAM Symposium on Segregation in Granular Flows: Proceedings of the IUTAM Symposium, pp. 31-51. Springer Netherlands, 2000.

[28] N. Thomas, Phys. Rev. E **62** 961 (2000).

[29] N. Thomas and U. D'Ortona, Phys. Rev. E **97**, 022903 (2018).

[30] S. Pillitteri, G. Lumay, E. Opsomer, and N. Vandewalle, Sci. Rep. **9** 7281 (2019).

[31] GDR MiDi, Euro. Phys. J. E **14**, 341 (2004).

[32] B. Utter, and R. P. Behringer, Phys. Rev. E **69**, 031308 (2004).

[33] R. P. Jones and A. B. Isner and H. Xiao and J. M. Ottino and P. B. Umbanhowar and R. M. Lueptow, Phys. Rev. Fluids **3**, 094304 (2018).

[34] C. R. J. Charles and Z. S. Khan and S. W. Morris, Granular Matter **8**, 1 (2006).

[35] S. J. Fiedor and P. Umbanhowar and J. M. Ottino, Phys. Rev. E **73**, 041303 (2006).

[36] H. Ramberg, Phys. Planet. Interiors **1**, 63 (1968).