# Anti-plane segregation and diffusion in dense, bidisperse granular shear
flow

Harkirat Singh and David L. Henann222Email address for correspondence:
<EMAIL_ADDRESS>
School of Engineering, Brown University, Providence, RI 02912, USA

###### Abstract

Many dense granular systems are non-monodisperse, consisting of particles of
different sizes, and will segregate based on size during flow. This phenomenon
is an important aspect of many industrial and geophysical processes,
necessitating predictive continuum models. This paper systematically studies a
key aspect of the three-dimensional nature of segregation and diffusion in
flowing, dense, bidisperse granular mixtures—namely, segregation and diffusion
acting along the direction perpendicular to the plane of shearing, which we
refer to as the anti-plane modes of segregation and diffusion. To this end, we
consider discrete-element method (DEM) simulations of flows of dense,
bidisperse mixtures of frictional spheres in an idealized configuration that
isolates anti-plane segregation and diffusion. We find that previously-
developed constitutive equations, calibrated to DEM simulation results from
flows in which both the segregation and diffusion processes occur within the
plane of shearing, do not capture aspects of the anti-plane segregation
dynamics. Accordingly, we utilize DEM simulation results to inform and
calibrate constitutive equations for the segregation and diffusion fluxes in
their anti-plane modes. Predictions of the resulting continuum model for the
anti-plane segregation dynamics are tested against additional DEM simulation
results across different cases, while parameters such as the shear strain rate
and mixture composition are varied, and we find that the calibrated model
predictions match well with the DEM simulation results. Finally, we suggest a
strategy for generalizing the constitutive forms for the segregation and
diffusion fluxes to obtain three-dimensional constitutive equations that
account for both the in-plane and anti-plane modes of the segregation and
diffusion processes.

##  1 Introduction

Dense granular mixtures consisting of particles of different sizes tend to
segregate, or demix, based on size during shear flow. This phenomena is an
important aspect of geophysical flows as well as industrial processes
involving blending of granular mixtures. Continuum models capable of
predicting the evolution of the spatially inhomogeneous composition of a
granular mixture are valuable in understanding or designing processes
involving size segregation. Consequently, substantial effort has gone into
developing continuum constitutive equations for segregation and diffusion
fluxes in flowing, dense granular mixtures based on both simulated and
experimental granular systems (e.g., Gray & Thornton, 2005; Gray & Chugunov,
2006; Fan & Hill, 2011; Marks et al., 2012; Hill & Tan, 2014; Fan et al.,
2014; Gajjar & Gray, 2014; Tunuguntla et al., 2014; Jones et al., 2018; Fry et
al., 2018; Trewhela et al., 2021; Liu et al., 2023; Singh et al., 2023),
primarily focusing on bidisperse mixtures made up of grains of two different
sizes.

The vast majority of works on continuum modeling of size segregation in dense,
bidisperse granular flows—including all of those cited in the previous
paragraph—have dealt with two-dimensional flow and segregation scenarios. More
precisely, denoting the plane of shearing as the $x$-$y$-plane, the components
of the velocity field are of the form $v_{x}(x,y)$, $v_{y}(x,y)$, and
$v_{z}=0$, and the vector describing the relative flux of large and small
grains—which accounts for both the segregation and diffusion processes—lies
entirely within the $x$-$y$-plane. Throughout this paper, we refer to the
segregation and diffusion processes in this scenario as their “in-plane”
modes. Common examples from the literature include flow down an inclined plane
(e.g., Gray & Thornton, 2005; Gray & Chugunov, 2006; Marks et al., 2012;
Gajjar & Gray, 2014; Tunuguntla et al., 2014; Singh et al., 2023), two-
dimensional bounded heap flow (e.g., Fan et al., 2014; Jones et al., 2018;
Schlick et al., 2015), and flow in a rotating drum (e.g., Hill & Tan, 2014;
Liu et al., 2019; Barker et al., 2021; Yang et al., 2021). However,
segregation and diffusion in many bidisperse shear flows—e.g., split-bottom
flow (Hill & Fan, 2008; Fan & Hill, 2010)—are not in their in-plane modes. It
has not been established whether flux constitutive equations proposed and
calibrated based on measurements from discrete simulations or experiments
involving the in-plane modes of segregation and diffusion remain predictive in
more general, three-dimensional scenarios, which serves as the motivation for
this work.

The purpose of this paper is to systematically study the three-dimensional
nature of the continuum constitutive equations for the segregation and
diffusion fluxes in flowing, dense, bidisperse granular mixtures, using
discrete-element method (DEM) simulations of frictional spheres. In
particular, we consider an idealized flow configuration, in which a semi-
infinite layer of dense granular material is sheared within the $x$-$y$-plane
but gravity and composition gradients act along the $z$-direction, so that the
segregation and diffusion fluxes act along the direction perpendicular to the
plane of shear. Throughout this paper, we refer to the segregation and
diffusion processes in this scenario as their “anti-plane” modes. We then test
whether a continuum model for the in-plane modes of segregation and diffusion
in a dense, bidisperse mixture is capable of capturing the evolution of the
spatially inhomogeneous composition of the granular mixture in this anti-plane
scenario. We focus on the continuum model proposed in our prior work (Singh et
al., 2023), which can predict the dynamics of segregation across different
cases of flow down an inclined plane as well as planar shear flow with
gravity. Observations of size segregation in both of these flow configurations
are examples of the in-plane mode of segregation. We find that the in-plane
segregation model fails to capture several features of the segregation
dynamics in the anti-plane mode—namely, the rate of diffusion and segregation
at sufficiently low shear strain rates (quantified through the inertial
number) and the scaling of the segregation flux with pressure. To remedy these
points, we first isolate the anti-plane mode of diffusion by removing
segregation due to gravity to motivate a modified constitutive equation for
the diffusion flux. Second, we return to the anti-plane mode of segregation
driven by a gravitational pressure gradient perpendicular to the plane of
shearing to inform a modified constitutive equation for the segregation flux.
When combined, the modified continuum model is capable of capturing the anti-
plane segregation dynamics across different combinations of shear strain rate,
mixture composition, and layer thickness.

This paper is organized as follows. We first discuss the continuum framework
used to describe dense, bidisperse granular systems in Section 2. Then, we
introduce the configuration for anti-plane segregation in simple shear flow in
Section 3 and illustrate the deficiencies of the planar model of Singh et al.
(2023). Next, based on DEM simulation data, we propose continuum constitutive
equations for the anti-plane modes of the diffusion and segregation fluxes in
Sections 4.1 and 4.2, respectively. The proposed continuum model is solved to
predict the transient evolution of the concentration field in anti-plane
segregation in simple shear flow, and continuum model predictions are tested
against DEM simulation results in Section 5. Revisiting anti-plane diffusion,
in Section 6.1, we characterize the diffusion flux through the mean square
displacement (MSD) of the granular mixture and contrast the result with the
diffusion flux characterized through the dynamics of mixing of an initially
segregated granular mixture from Section 4.1. Lastly, we describe a strategy
for generalizing the constitutive equations for the segregation and diffusion
fluxes to obtain three-dimensional constitutive forms that are capable of
switching between the respective in-plane and anti-plane modes in Section 6.2.
We conclude the paper with some closing remarks in Section 7.

##  2 Continuum framework

In this section, we discuss the continuum framework used to describe
segregation and diffusion in dense, bidisperse granular mixtures. Regarding
notation, we use component notation, in which the components of vectors,
$\bm{v}$, and tensors, $\bm{\sigma}$, relative to a set of Cartesian basis
vectors $\\{\bm{e}_{i}|i=1,2,3\\}$ are denoted by $v_{i}$ and $\sigma_{ij}$,
respectively, and the Einstein summation convention is employed.

We consider bidisperse granular mixtures consisting of grains of two different
sizes—large grains with diameter $d^{\rm l}$ and small grains with diameter
$d^{\rm s}$. All grains are made of same material with mass density $\rho_{\rm
s}$, thereby eliminating density-based segregation and isolating size-based
segregation. To describe the dynamics of segregation in dense, bidisperse
granular systems at the continuum level, we utilize a mixture-theory approach,
which is common in the literature (e.g., Gray, 2018; Umbanhowar et al., 2019).
The solid volume fractions of each species are denoted as $\phi^{\rm l}$ and
$\phi^{\rm s}$ for the large and small grains, respectively. For the dense
shear flows considered in the present work, volume changes are minimal, and
any volume change at flow initiation occurs over a much shorter time scale
than the process of segregation. Accordingly, we adopt the idealization that
the total solid volume fraction $\phi=\phi^{\rm l}+\phi^{\rm s}$ is constant
and, based on observations from DEM simulations, approximately equal to
$\phi=0.6$ for dense spheres. The concentrations of each species are $c^{\rm
l}=\phi^{\rm l}/\phi$ and $c^{\rm s}=\phi^{\rm s}/\phi$, so that $c^{\rm
l}+c^{\rm s}=1$. The average grain size of the mixture is then defined as
$\bar{d}=c^{\rm l}d^{\rm l}+c^{\rm s}d^{\rm s}$.

The mixture velocity $v_{i}$ is defined through the partial velocities of each
species, $v_{i}^{\rm l}$ and $v_{i}^{\rm s}$, as $v_{i}=c^{\rm l}v^{\rm
l}_{i}+c^{\rm s}v^{\rm s}_{i}$. Then, the strain-rate tensor for the mixture
is $D_{ij}=(1/2)\left(\partial v_{i}/\partial x_{j}+\partial v_{j}/\partial
x_{i}\right)$, and the equivalent shear strain rate is
$\dot{\gamma}=(2D_{ij}D_{ij})^{1/2}$. Due to the constant-volume idealization,
the mixture velocity is divergence-free, $\partial v_{i}/\partial x_{i}=0$,
and the mixture strain-rate tensor is deviatoric, $D_{kk}=0$. The relevant
stress-related fields for the mixture are the symmetric Cauchy stress tensor
$\sigma_{ij}=\sigma_{ji}$ and the pressure $P=-(1/3)\sigma_{kk}$. The inertial
number for a bidisperse system is defined in terms of the average grain size
as $I=\dot{\gamma}\bar{d}\sqrt{\rho_{\rm s}/P}$ (Rognon et al., 2007; Tripathi
& Khakhar, 2011).

The relative volume flux of species $\nu={\rm l},{\rm s}$ is defined as
$w^{\nu}_{i}=c^{\nu}(v^{\nu}_{i}-v_{i})$, so that $w^{\rm l}+w^{\rm s}=0$. We
utilize $c^{\rm l}$ as the field variable that describes the segregation
dynamics, and under the constant-volume idealization, the mass conservation
equation for the large grains governs the evolution of $c^{\rm l}$ as follows:

$\frac{Dc^{\rm l}}{Dt}+\frac{\partial w^{\rm l}_{i}}{\partial x_{i}}=0,$ (2.1)

where $D(\bullet)/Dt$ is the material time derivative. To close the
segregation model, a constitutive equation for the relative volume flux of
large particles $w^{\rm l}_{i}$ is needed. We take the flux to be comprised of
two contributions:

$w^{\rm l}_{i}=w^{\rm diff}_{i}+w^{\rm P}_{i},$ (2.2)

where $w^{\rm diff}_{i}$ is the diffusion flux that acts to remix the species
and $w^{\rm P}_{i}$ is the pressure-gradient-driven segregation flux that acts
to demix the species. In this paper, we focus on simple shear flows, in which
the equivalent shear strain rate $\dot{\gamma}$ is spatially uniform, so the
shear-strain-rate-gradient-driven segregation of Liu et al. (2023) does not
appear in (2.2).

##  3 Anti-plane segregation in simple shear flow

###  3.1 Discrete-element method simulations

In this section, we examine the anti-plane mode of segregation by considering
simple shear flow of a layer of a dense, bidisperse granular mixture within
the $x$-$y$-plane. The layer has a finite thickness $H$ along the
$z$-direction, which is perpendicular to the plane of shearing, and gravity
$G$ acts along the $z$-direction. We consider simple shear flows of dense,
bidisperse spheres with an inter-particle friction coefficient of $\mu_{\rm
surf}=0.4$ using discrete-element method (DEM) simulations. Further details of
the simulated granular system are given in Appendix A. The DEM setup for the
case of a well-mixed system of bidisperse spheres and $H=30\bar{d}_{0}$ is
shown in Fig. 1(a), where $\bar{d}_{0}$ is the system-wide average grain size.
Shearing in the $x$-$y$ plane is achieved through Lees-Edwards boundary
conditions (Lees & Edwards, 1972). The dimensions of the DEM simulation domain
are taken to be $L=20\bar{d}_{0}$ along the $x$-direction and
$W=10\bar{d}_{0}$ along the $y$-direction. The bottom boundary is a flat and
frictionless wall that does not impart any shear traction within its plane,
and the top surface is traction-free. This setup achieves a homogeneous state
of strain rate, in which the only non-zero components of the strain-rate
tensor are $D_{xy}=D_{yx}=\dot{\gamma}/2$, and eliminates segregation due to
strain-rate gradients (e.g., Liu et al., 2023), isolating pressure-gradient-
driven segregation in its anti-plane mode. As discussed in Section 1, we use
the term “anti-plane” to designate that pressure-gradient-driven size
segregation acts along the direction (the $z$-direction) perpendicular to the
plane of shearing (the $x$-$y$ plane), which contrasts with the “in-plane”
mode examined in our prior work (Singh et al., 2023), where segregation acts
along a direction that is within the plane of shearing.

Figure 1: Representative base case for anti-plane segregation in simple shear
flow. (a) Initial well-mixed configuration for three-dimensional DEM
simulation of bidisperese spheres. (b) Segregated state after a simulation
time of $\tilde{t}=\dot{\gamma}t=1200$. (c) Spatiotemporal evolution of the
$c^{\rm l}$ field. (d) Spatial profiles of the $c^{\rm l}$ field at three
different time instants ($\tilde{t}=60$, 300, and 1200) as indicated by the
three vertical lines in (c).

The important dimensionless parameters that characterize anti-plane
segregation in simple shear flow are (1) the dimensionless layer thickness
$H/\bar{d}_{0}$, (2) the dimensionless strain rate
$\dot{\gamma}\sqrt{\bar{d}_{0}/G}$, (3) the initial large-grain concentration
field $c^{\rm l}_{0}(z)$, and (4) the bidisperse grain-size ratio $d^{\rm
l}/d^{\rm s}$. DEM simulation results for a representative base case
corresponding to the parameter set
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$ are summarized in Fig 1. The
segregated state after a simulation time of $\tilde{t}=\dot{\gamma}t=1200$ is
shown in Fig. 1(b), showing that the large grains (dark gray) indeed segregate
to the top free surface, leaving small grains (light gray) on the bottom of
the layer. To further quantify the process of segregation, we coarse-grain the
large-grain concentration field $c^{\rm l}$ and plot spatiotemporal contours
of the evolution of $c^{\rm l}$ in Fig. 1(c). (Coarse-graining procedures are
described in Appendix A.) The growing width of the dark region indicates the
accumulation of large grains at the top of the layer. Finally, in Fig. 1(d),
we plot spatial profiles of the $c^{\rm l}$ field along the $z$-direction at
three snapshots in time, associated with short, moderate, and long times
($\tilde{t}=60,300,1200$) with respect to the process of segregation,
indicated by the vertical lines in Fig. 1(c).

###  3.2 Continuum modeling based on in-plane constitutive equations

Next, we test the segregation model described in Singh et al. (2023) for the
in-plane modes of segregation and diffusion in a bidisperse mixture in the
context of anti-plane segregation and diffusion. First, we briefly recap the
constitutive equation for the relative volume flux $w^{\rm l}_{i}$, which
comprises contributions due to diffusion and segregation. The diffusion flux
$w^{\rm diff}_{i}$ acts to remix the two species and is taken to be driven by
concentration gradients as $w^{\rm diff}_{i}=-D(\partial c^{\rm l}/\partial
x_{i})$, where $D$ is the binary diffusion coefficient which scales with the
mean grain size $\bar{d}$ and the shear strain rate $\dot{\gamma}$ as
$D\sim\bar{d}^{2}\dot{\gamma}$ (Utter & Behringer, 2004; Fry et al., 2019;
Bancroft & Johnson, 2021; Liu et al., 2023). Accordingly, the diffusion flux
is

$w^{\rm diff}_{i}=-C_{\rm diff}^{\rm in}\bar{d}^{2}\dot{\gamma}\dfrac{\partial
c^{\rm l}}{\partial x_{i}},$ (3.1)

where $C_{\rm diff}^{\rm in}$ is a dimensionless, constant material parameter
that characterizes the in-plane diffusion flux. Regarding the segregation
flux, in the present work, we consider homogeneous shear flows, in which the
shear strain rate $\dot{\gamma}$ is spatially uniform. Therefore, shear-
strain-rate gradients (e.g., Liu et al., 2023) do not drive the segregation
process, and we only consider segregation driven by pressure gradients (e.g.,
Singh et al., 2023) and adopt the following constitutive equation for the
segregation flux $w^{\rm P}_{i}$:

$w^{\rm P}_{i}=-C^{\rm P,in}_{\rm seg}\dfrac{\bar{d}^{2}\dot{\gamma}}{P}c^{\rm
l}(1-c^{\rm l})(1-\alpha+\alpha c^{\rm l})\dfrac{\partial P}{\partial x_{i}},$
(3.2)

where $C^{\rm P,in}_{\rm seg}$ and $\alpha$ are dimensionless, constant
parameters that characterize the in-plane segregation flux. The scaling of the
pre-factor in (3.2) with $\bar{d}^{2}\dot{\gamma}/P$ follows prior works in
the literature (Trewhela et al., 2021; Barker et al., 2021; Singh et al.,
2023) and ensures dimensional consistency. Moreover, this scaling captures the
requirement that the segregation flux is zero when there is no flow
$(\dot{\gamma}=0)$ as well as the attenuation of segregation with increasing
pressure (Golick & Daniels, 2009; Fry et al., 2018). The dependence of the
pre-factor in (3.2) on $c^{\rm l}$ through the function $f(c^{\rm l})=c^{\rm
l}(1-c^{\rm l})(1-\alpha+\alpha c^{\rm l})$ follows from Gajjar & Gray (2014)
and captures the asymmetric dependence of the flux on $c^{\rm l}$, quantified
by the dimensionless parameter $\alpha\in[0,1]$. For $\alpha=0$, the
dependence of $f(c^{\rm l})$ is symmetric about its maximum at $c^{\rm
l}=0.5$, but for $\alpha\in(0,1]$, the maximum of $f(c^{\rm l})$ skews from
$c^{\rm l}=0.5$ towards $c^{\rm l}=1$, while maintaining that segregation
ceases when the bidisperse mixture becomes either all large ($c^{\rm l}=1$) or
all small ($c^{\rm l}=0$) grains. While it is well established in the
literature that the pressure-gradient-driven segregation flux depends on the
grain-size ratio $d^{\rm l}/d^{\rm s}$ (e.g., Schlick et al., 2015; Jones et
al., 2018; Trewhela et al., 2021), we focus on a single grain-size ratio of
$d^{\rm l}/d^{\rm s}=1.5$ throughout this work, and therefore, this dependence
does not appear in (3.2). The relative volume flux $w^{\rm l}_{i}$ is then
taken to be the sum of (3.1) and (3.2): $w^{\rm l}_{i}=w^{\rm diff}_{i}+w^{\rm
P}_{i}$, where the material parameters associated with the in-plane
segregation model for a given $d^{\rm l}/d^{\rm s}$ are $\\{C_{\rm diff}^{\rm
in},C^{\rm P,in}_{\rm seg},\alpha\\}$.

Utilizing the flux constitutive equations (3.1) and (3.2) in the balance of
mass equation (2.1), we obtain a governing equation for the $c^{\rm l}$ field
in the absence of shear-strain-rate gradients. In the context of anti-plane
segregation in simple shear flow (Fig. 1(a)), in which the concentration and
pressure fields only vary along the $z$-direction, the large-grain
concentration field $c^{\rm l}(z,t)$ is governed by the following partial
differential equation (PDE):

$\frac{\partial c^{\rm l}}{\partial t}+\frac{\partial}{\partial
z}\left(-C_{\rm diff}^{\rm in}\bar{d}^{2}\dot{\gamma}\frac{\partial c^{\rm
l}}{\partial z}-{C}^{\rm P,in}_{\rm
seg}\frac{\bar{d}^{2}\dot{\gamma}}{P}c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha
c^{\rm l})\frac{\partial P}{\partial z}\right)=0.$ (3.3)

In this configuration, the strain rate $\dot{\gamma}$ is prescribed and is
spatially and temporally uniform. Moreover, the stress field does not vary in
time, and the out-of-plane normal stress component $\sigma_{zz}(z)$ may be
determined from a static force balance along the $z$-direction to be
$\sigma_{zz}(z)=-\phi\rho_{\rm s}Gz$. The in-plane normal stresses
$\sigma_{xx}(z)$ and $\sigma_{yy}(z)$ are equal and, in coarse-grained DEM
data, are observed to also be linearly proportional to $z$; however, due to
the normal stress differences that arise in dense flows of spheres (Depken et
al., 2007; Srivastava et al., 2021), the magnitude of the slope of this linear
dependence is observed to be slightly higher than $\phi\rho_{\rm s}G$. As a
result, the pressure field $P(z)$ varies linearly in $z$ but with a constant
slope $\partial P/\partial z$ that is slightly higher than the nominal value
of $\phi\rho_{\rm s}G$. To control for this effect, when obtaining continuum
solutions, we utilize the value of the slope of the pressure field $\partial
P/\partial z$ obtained from the coarse-grained stress fields in the DEM data,
rather than the nominal value of $\phi\rho_{\rm s}G$. Furthermore, to avoid a
singularity in the pressure-gradient-driven segregation flux at the free
surface where $z=0$, we add a small constant to the pressure field
corresponding to the weight of a layer of $(1/4)\bar{d}_{0}$ thickness, i.e.,
$(1/4)(\partial P/\partial z)\bar{d}_{0}$. We have verified that continuum
model predictions are insensitive to the exact choice of this constant, so
long as it is sufficiently small. Therefore, with $P(z)=(\partial P/\partial
z)z+(1/4)(\partial P/\partial z)\bar{d}_{0}$, we have that $(1/P)\partial
P/\partial z=1/(z+(1/4)\bar{d}_{0})$ in the last term of (3.3).

When accompanied by boundary and initial conditions, (3.3) may then be used to
obtain predictions for the transient evolution of the large-grain
concentration field $c^{\rm l}(z,t)$ during anti-plane segregation in simple
shear flow. Regarding boundary conditions, we impose no-flux boundary
conditions at the top free surface as well as at the bottom frictionless wall,
i.e., $w^{\rm l}_{z}=-{C}_{\rm diff}^{\rm in}\bar{d}^{2}\dot{\gamma}({\partial
c^{\rm l}}/{\partial z})-{C}^{\rm P,in}_{\rm
seg}({\bar{d}^{2}\dot{\gamma}}/{P})c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha
c^{\rm l})({\partial P}/{\partial z})=0$ at $z=0$ and $H$. For the initial
condition for the large-grain concentration field $c^{\rm l}_{0}(z)=c^{\rm
l}(z,t=0)$, we obtain the coarse-grained field from the initial DEM
configuration at $t=0$ and use this field as the initial condition $c^{\rm
l}_{0}(z)$ in the corresponding continuum simulation. This is done to control
for spatial concentration fluctuations in the initial DEM state. Then, for a
given case of anti-plane segregation in simple shear flow, described by the
parameter set $\\{H/\bar{d}_{0},\dot{\gamma}\sqrt{\bar{d}_{0}/G},c^{\rm
l}_{0},d^{\rm l}/d^{\rm s}\\}$, we obtain numerical solutions using finite
differences. Central differences are used to discretize the spatial
derivatives in (3.3). The Euler method is used for temporal discretization, in
which the spatial derivative of $c^{\rm l}$ appearing in the diffusion flux
term in (3.3) is treated implicitly for improved numerical stability while the
pre-factors of both flux terms in (3.3) are treated explicitly. We use
sufficiently small spatial and temporal resolutions in the finite-difference
scheme to ensure stable and accurate results.

We compare predictions of the in-plane-based continuum model (3.3) against DEM
data for the representative base case of anti-plane segregation in simple
shear flow from Fig. 1,
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$, using the in-plane parameters from
Singh et al. (2023), i.e., $\\{{C}_{\rm diff}^{\rm in}=0.045,{C}^{\rm
P,in}_{\rm seg}=0.34,\alpha=0.4\\}$. Comparisons of the continuum model
predictions (dashed gray lines) with the DEM simulations (solid black lines)
are summarized in Fig. 2 at four different time instants during the
segregation process: $\tilde{t}=60$, 300, 600, and 1200. The in-plane model
falls short in capturing the transient evolution of the segregation dynamics
in the anti-plane mode of segregation. First, predictions of the segregation
dynamics lag behind the DEM simulation results in time, and second, a sharp
transition zone near the top is predicted by the continuum model at early
times $(\tilde{t}\lesssim 300)$, which is contrary to DEM simulation
observations that show a comparatively wider transition zone. The main
objective of this paper to adapt the constitutive equations for diffusion and
pressure-gradient-driven segregation to capture segregation dynamics in the
anti-plane mode. To do so, we examine the anti-plane mode of diffusion in
Section 4.1 and the anti-plane mode of pressure-gradient-driven segregation in
Section 4.2.

Figure 2: Comparisons of predictions of the in-plane continuum model of Singh
et al. (2023) with corresponding DEM simulation results for the transient
evolution of the segregation dynamics for the base case of anti-plane
segregation in simple shear flow
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$. On the left are spatiotemporal
contours of the evolution of $c^{\rm l}$ measured in the DEM simulation, and
on the right are comparisons of DEM simulation results (solid black lines) and
continuum model predictions (dashed gray lines) for the $c^{\rm l}$ field at
four different time instants during the segregation process: $\tilde{t}=60$,
300, 600, and 1200 in the sequence of top left, top right, bottom left, bottom
right and indicated by the vertical lines on the contour plot.

##  4 Constitutive equations for anti-plane diffusion and segregation

###  4.1 Anti-plane diffusion flux

In this section, we characterize the anti-plane mode of diffusion—i.e., when
the diffusion process occurs along the direction perpendicular to the plane of
shearing. To isolate the anti-plane diffusion flux, steady simple shear flow
at a uniform strain rate $\dot{\gamma}$ is prescribed within the $x$-$y$-plane
at a constant pressure $P$, while the system is initially segregated along the
$z$-direction (the anti-plane direction). This is in contrast to the situation
in Section 3, where the pressure field was not uniform but hydrostatic. A
snapshot of the initial DEM configuration is shown in Fig. 3(a), where large
grains (dark gray) are on the top and small grains (light gray) are on the
bottom, and the grain-size ratio is $d^{\rm l}/d^{\rm s}=1.5$. We employ Lees-
Edwards boundary conditions within the $x$-$y$-plane to prescribe a uniform
strain rate $\dot{\gamma}$, and the dimensions of the simulation domain along
the $x$\- and $y$-directions are $L=20\bar{d}_{0}$ and $W=10\bar{d}_{0}$,
respectively. The size of the simulation domain along the $z$-direction is
$H=60\bar{d}_{0}$, and periodic boundary conditions are not employed along
this direction. Instead, the bottom boundary consists of a flat, frictionless
wall, and a constant compressive normal stress is applied to the top through a
mass of large grains (not pictured in Fig. 3(a)) subjected to a gravitational
body force along the $z$-direction. The grains in the region shown in Fig.
3(a) are not subjected to any body forces. The result of these boundary
conditions along the anti-plane direction is that the normal stress component
$\sigma_{zz}$ is spatially and temporally constant and equal to $-P_{\rm
top}$. Moreover, the normal stresses $\sigma_{xx}$ and $\sigma_{yy}$ are also
uniform and compressive but of slightly higher magnitude due to normal stress
differences that arise in dense flows of spheres. The spatially uniform
pressure field is then $P=-(1/3)(\sigma_{xx}+\sigma_{yy}+\sigma_{zz})$ and is
slightly greater than $P_{\rm top}$.

Since the strain rate and pressure are uniform in this flow configuration, no
size segregation takes place, and only diffusive mixing occurs. A DEM snapshot
of an intermediate state after a simulation time of
$\tilde{t}=\dot{\gamma}t=590$ is shown in Fig. 3(b), illustrating that the
initially sharp transition between large and small grains in Fig. 3(a) becomes
diffuse due to mixing. To further quantify the diffusion process, we coarse-
grain the $c^{\rm l}$ field from the DEM simulations, and the spatiotemporal
contours of the evolution of the $c^{\rm l}$ field are plotted in Fig. 3(c)
for the case of a nominal inertial number of
$\dot{\gamma}\bar{d}_{0}\sqrt{\rho_{\rm s}/P}=0.02$, where
$\bar{d}_{0}=(d^{\rm l}+d^{\rm s})/2$ and $P$ is the pressure obtained from
the coarse-grained DEM results to control for the normal stress differences.
We observe that the width of the transition region between large and small
grains, which is initially sharp, grows in time as the diffusion process
proceeds and the grains mix.

Figure 3: Representation of the anti-plane mode of diffusive mixing in simple
shear flow. (a) Initially segregated configuration for three-dimensional DEM
simulation of bidisperse spheres with $\sim 20~{}000$ grains. (b) Intermediate
remixed state after shearing for a simulation time of
$\tilde{t}=\dot{\gamma}t=590$. (c) Spatiotemporal evolution of the large-grain
concentration field for a nominal inertial number of
$\dot{\gamma}\bar{d}_{0}\sqrt{\rho_{\rm s}/P}=0.02$. (d) Estimated dependence
of the anti-plane diffusion parameter $\hat{C}_{\rm diff}^{\rm anti}$ on the
nominal inertial number $\dot{\gamma}\bar{d}_{0}\sqrt{\rho_{\rm s}/P}$.
Symbols are estimates based on fitting the continuum model (4.2) to DEM data
for different values of the nominal inertial number; the solid black curve
represents the best fit of (4.3) with $\\{C_{0}=0.104,C_{\infty}=0.057,I_{\rm
tr}=0.00744,m=3.23\\}$.

In the continuum model for anti-plane diffusion, we continue to take the
diffusion flux $w^{\rm diff}_{i}$ to scale with the average grain size
$\bar{d}$, the strain rate $\dot{\gamma}$, and the concentration gradient
$\partial c^{\rm l}/\partial x_{i}$ as in (3.1), so that

$w^{\rm diff}_{i}=-\hat{C}_{\rm diff}^{\rm
anti}\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm l}}{\partial x_{i}}.$ (4.1)

The key difference from (3.1) is that we allow the dimensionless parameter
associated with diffusion $\hat{C}_{\rm diff}^{\rm anti}$ to depend on
additional dimensionless quantities—in particular, the inertial number
$I=\dot{\gamma}\bar{d}\sqrt{\rho_{\rm s}/P}$. Plugging (4.1) in the mass
conservation equation (2.1) in the context of anti-plane diffusion in which
the concentration field only varies along the $z$-direction, we obtain the
following governing equation for the diffusion process:

$\frac{\partial c^{\rm l}}{\partial t}-\frac{\partial}{\partial
z}\left(\hat{C}_{\rm diff}^{\rm anti}\bar{d}^{2}\dot{\gamma}\frac{\partial
c^{\rm l}}{\partial z}\right)=0.$ (4.2)

To estimate the parameter $\hat{C}_{\rm diff}^{\rm anti}$ and its functional
dependence on $I$, we perform DEM simulations of anti-plane diffusive mixing
for 11 different values of the nominal inertial number
$\dot{\gamma}\bar{d}_{0}\sqrt{\rho_{\rm s}/P}$ spanning from approximately
0.001 to 0.1. We then calculate predictions of (4.2) for the concentration
field $c^{\rm l}(z,t)$ for different values of $\hat{C}_{\rm diff}^{\rm anti}$
over the range $\hat{C}_{\rm diff}^{\rm anti}\in(0.04,0.12)$. By minimizing
the difference between the continuum model predictions and the coarse-grained
concentration field from each DEM simulation, the best-fit value of
$\hat{C}_{\rm diff}^{\rm anti}$ is determined for each value of the nominal
inertial number $\dot{\gamma}\bar{d}_{0}\sqrt{\rho_{\rm s}/P}$. In the fitting
process, we only consider sufficiently short dimensionless times $\tilde{t}$,
corresponding to a transition width that remains substantially smaller than
the layer thickness $H$, so that the diffusion process is unaffected by the
boundaries. We note that this fitting process neglects the spatial variation
in the inertial number field $I=\dot{\gamma}\bar{d}\sqrt{\rho_{\rm s}/P}$ that
arises due to the spatial variation of $c^{\rm l}$ and hence $\bar{d}$.
However, this simplification has a minimal effect on the estimates of
$\hat{C}_{\rm diff}^{\rm anti}$ since, for a given nominal inertial number,
the inertial number field varies by a factor of $d^{\rm l}/d^{\rm s}=1.5$
between the small grains at the bottom of the layer and the large grains on
the top of the layer, while the inertial-number range-of-interest spans two
orders of magnitude. Moreover, this simplification is only invoked during the
fitting process, and, in all subsequent calculations, $\hat{C}_{\rm diff}^{\rm
anti}$ is taken to be a function of the spatially varying inertial number
field and not a global, nominal value of the inertial number.

The best-fit values of $\hat{C}_{\rm diff}^{\rm anti}$ are plotted as a
function of the nominal inertial number
$\dot{\gamma}\bar{d}_{0}\sqrt{\rho_{\rm s}/P}$ as light gray circles in Fig.
3(d). The estimated $\hat{C}_{\rm diff}^{\rm anti}$ plateaus to a lower value
at high values of the nominal inertial number and to a higher value at low
values of the nominal inertial number, and there is an intermediate range in
which $\hat{C}_{\rm diff}^{\rm anti}$ transitions from the higher plateau to
the lower plateau. We use the following phenomenological functional form to
fit the estimated values of the diffusion parameter over the range
$0.001\lesssim I\lesssim 0.1$:

$\hat{C}^{\rm anti}_{\rm diff}(I)=C_{0}+\dfrac{C_{\infty}-C_{0}}{(I_{\rm
tr}/I)^{m}+1},$ (4.3)

where $C_{0}$ is the plateau value at low inertial numbers, $C_{\infty}$ is
the plateau value at high inertial numbers, $I_{\rm tr}$ determines the
position of the transition, and the exponent $m$ determines the width of the
transition with greater values of $m$ corresponding to sharper transitions.
Fitting (4.3) to the DEM-based estimates in Fig. 3(d) gives the parameter set
associated with the anti-plane mode of diffusion for frictional spheres to be
$\\{C_{0}=0.104,C_{\infty}=0.057,I_{\rm tr}=0.00744,m=3.23\\}$. The fitted
functional form is shown by the solid black line in Fig. 3(d). In the semi-
logarithmic plot of Fig. 3(d), (4.3) corresponds to hyperbolic-tangent-type
behavior.

The dependence of the anti-plane diffusion parameter $\hat{C}^{\rm anti}_{\rm
diff}$ on the inertial number $I$ is in contrast to the in-plane diffusion
parameter, which, in our experience (Liu et al., 2023; Singh et al., 2023),
may be idealized as independent of the inertial number, based on both mean
square displacement DEM data and continuum-model predictions of the evolution
of the $c^{\rm l}$ field. We note that the plateau value at high inertial
numbers of $C_{\infty}=0.057$ is similar to the in-plane diffusion parameter
$C_{\rm diff}^{\rm in}=0.045$ determined in Liu et al. (2023) and used in
Section 3.2, so that diffusion is nearly isotropic for $I\gtrsim 0.01$.
However, we find that out-of-plane diffusion becomes comparatively stronger
for low values of the inertial number with $\hat{C}^{\rm anti}_{\rm diff}$
increasing to a plateau value of $C_{0}=0.104$, indicating anisotropy in
diffusion at low values of the inertial number. The dependence of
$\hat{C}^{\rm anti}_{\rm diff}$ on $I$ is also in contrast to the mean square
displacement DEM data of Bancroft & Johnson (2021), who found the diffusion
coefficient along the anti-plane direction to be nearly constant with
$\hat{C}^{\rm anti}_{\rm diff}\approx 0.04$ over a range of inertial numbers.
We revisit this point in Section 6.1.

###  4.2 Anti-plane segregation flux

Next, we return to anti-plane segregation in simple shear flow as described in
Section 3 to characterize the anti-plane mode of segregation. We take the
pressure-gradient-driven segregation flux to be given in the following form:

$w^{\rm P}_{i}=-\hat{C}^{\rm P,anti}_{\rm
seg}\frac{\bar{d}^{2}\dot{\gamma}}{P}c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha
c^{\rm l})\dfrac{\partial P}{\partial x_{i}},$ (4.4)

where $\hat{C}^{\rm P,anti}_{\rm seg}$ is a dimensionless quantity, which can
depend on additional dimensionless parameters, and the pre-factor
$({\bar{d}^{2}\dot{\gamma}}/{P})c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha c^{\rm
l})$ has the same form as that used in the in-plane model summarized in
Section 3.2, wherein $\alpha$ is a constant parameter.

We hypothesize that the quantity $\hat{C}^{\rm P,anti}_{\rm seg}$ depends on
the inertial number $I$ (as does $\hat{C}^{\rm anti}_{\rm diff}$) and an
additional dimensionless parameter $J={P}/({\bar{d}|\partial P/\partial{\bf
x}|})$, which quantifies the ratio of the pressure to the change in pressure
over the mean grain size at a point. Inclusion of $J$-dependence in the
constitutive equation for the pressure-gradient-driven segregation flux (4.4)
allows for potential nonlinear dependence of the flux on the magnitude of the
pressure gradient. The procedure to estimate the dimensionless parameter
$\hat{C}^{\rm P,anti}_{\rm seg}$ and its functional dependence is as follows.
First, we run a given DEM simulation of anti-plane segregation in simple shear
flow to steady state, so that the bidisperse mixture is fully segregated, and
the total flux at each $z$-position is zero, i.e., $w^{\rm l}_{z}=w^{\rm
diff}_{z}+w^{\rm P}_{z}=0$. Therefore, the flux balance

$-\hat{C}_{\rm diff}^{\rm anti}(I)\dfrac{\partial c^{\rm l}}{\partial
z}=\hat{C}^{\rm P,anti}_{\rm seg}(I,J)\frac{1}{P}c^{\rm l}(1-c^{\rm
l})(1-\alpha+\alpha c^{\rm l})\dfrac{\partial P}{\partial z}$ (4.5)

is valid at every $z$-position. We find that the steady-state concentration
field, and hence the flux balance (4.5), is independent of the dimensionless
strain rate $\dot{\gamma}\sqrt{\bar{d}_{0}/G}$ when the other parameters
$\\{H/\bar{d}_{0},c^{\rm l}_{0},d^{\rm l}/d^{\rm s}\\}$ are fixed. Therefore,
(4.5) must be independent of the inertial number, and as a result, the
function $\hat{C}^{\rm P,anti}_{\rm seg}(I,J)$ should depend on the inertial
number in the same manner as $\hat{C}_{\rm diff}^{\rm anti}(I)$. Accordingly,
we take $\hat{C}^{\rm P,anti}_{\rm seg}(I,J)$ to be given by

$\hat{C}^{\rm P,anti}_{\rm
seg}(I,J)=\tilde{C}(J)\left[1-\left(\frac{C_{\infty}-C_{0}}{C_{\infty}+C_{0}}\right)\dfrac{(I_{\rm
tr}/I)^{m}-1}{(I_{\rm tr}/I)^{m}+1}\right],$ (4.6)

where $\tilde{C}(J)$ is a function of only $J$. The $I$-dependent function in
brackets in (4.6) is the same as (4.3) but normalized by
$(C_{\infty}+C_{0})/2$, and the parameters $\\{C_{0},C_{\infty},I_{\rm
tr},m\\}$ are the same as those estimated in Section 4.1.

Then, plugging (4.6) into (4.5) and rearranging, we have that, at steady
state, $\tilde{C}(J)$ is given by

$\tilde{C}(J)=\frac{-\left(C_{0}+C_{\infty}\right)P\,\partial c^{\rm
l}/\partial z}{2c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha c^{\rm l})\partial
P/\partial z}.$ (4.7)

The field quantities appearing on the right-hand-side of this expression may
be obtained by coarse-graining steady-state DEM data, which may be used to
determine the function $\tilde{C}(J)$. For a given case of anti-plane
segregation in simple shear flow, we consider $z$-positions within the
intermediate transition region of the steady-state concentration field where
$c^{\rm l}\in(0.2,0.8)$ and calculate the steady-state concentration and
pressure fields and their gradients at a given $z$-position from DEM data as
described in Appendix A. From the coarse-grained, steady-state DEM data,
$\tilde{C}$, using (4.7), and $J={P}/({\bar{d}|\partial P/\partial z|})$ are
calculated at each $z$-position and plotted against one another in Fig. 4,
where each symbol represents a unique $z$-position. We apply this procedure to
steady-state DEM data from 13 different cases of anti-plane segregation in
simple shear flow by changing the size of the system $H/\bar{d}_{0}$, the
strain rate $\dot{\gamma}\sqrt{\bar{d}_{0}/G}$, and the initial conditions
$c^{\rm l}_{0}$, and the data from all cases are included in Fig. 4. (The
grain-size ratio is maintained at $d^{\rm l}/d^{\rm s}=1.5$ throughout.) We
consider values of $\alpha$ within the range $\alpha\in(0,0.8)$ and find that
the best collapse is obtained for $\alpha=0.4$, which is the same as the value
determined for the in-plane mode of segregation in Singh et al. (2023). The
resulting steady-state data collapse for $\alpha=0.4$ across all 13 cases is
shown in Fig. 4. The data can be fitted using a power-law function, i.e.,
$\tilde{C}(J)=CJ^{\rm n}$, so that the anti-plane segregation parameter
becomes

$\hat{C}_{\rm seg}^{\rm P,anti}(I,J)=CJ^{\rm
n}\left[1-\left(\frac{C_{\infty}-C_{0}}{C_{\infty}+C_{0}}\right)\dfrac{(I_{\rm
tr}/I)^{m}-1}{(I_{\rm tr}/I)^{m}+1}\right],$ (4.8)

where the associated parameters are estimated to be $\\{C=0.22,n=0.58\\}$. The
fitted power-law function is shown by the solid black curve in Fig. 4. (An
alternative fitting function $\tilde{C}(J)=C_{1}/(1+C_{2}/J)$ with constants
$\\{C_{1},C_{2}\\}$ based on the works of Trewhela et al. (2021) and Barker et
al. (2021) is discussed in Appendix B and shown by the dash-dotted black curve
in Fig. 4.) The role of $J$-dependence in the segregation flux (4.4) and its
effect on the predicted concentration fields is discussed further in Section
5.

Figure 4: Steady-state collapse of the parameter $\tilde{C}$ (4.7) versus the
non-dimensional parameter $J={P}/({\bar{d}|\partial P/\partial z|})$ over $13$
different DEM simulations of anti-plane segregation in simple shear flow.
Symbols represent coarse-grained, steady-state DEM field data. The solid black
curve is the fit function $\tilde{C}=CJ^{n}$ with $\\{C=0.22,n=0.58\\}$, and
the dash-dotted black curve is the fit function $\tilde{C}=C_{1}/(1+C_{2}/J)$
with $\\{C_{1}=3.2,C_{2}=28.62\\}$.

##  5 Validation tests in the transient regime

In Section 4.1, we studied anti-plane diffusion in simple shear flow by
considering an initially segregated mixture and introduced four parameters
$\\{C_{0},C_{\infty},I_{\rm tr},m\\}$ to characterize the rate-dependence of
the anti-plane mode of diffusion. Then, in Section 4.2, we used steady-state
DEM data from simple shear flow in the presence of an anti-plane gravitational
pressure gradient to inform and assess a constitutive equation for the
pressure-gradient-driven segregation flux in the anti-plane mode. We
introduced three parameters $\\{C,n,\alpha\\}$ to characterize the anti-plane
mode of the segregation flux. In this section, we compare continuum model
predictions of the transient evolution of the large-grain concentration field
against transient DEM data for several cases of anti-plane segregation in
simple shear flow. Throughout, we continue to use the following fixed set of
parameters for frictional spheres:

$\\{C_{0}=0.104,C_{\infty}=0.057,I_{\rm
tr}=0.00744,m=3.23,C=0.22,n=0.58,\alpha=0.4\\}.$ (5.1)

First, we briefly recap how the continuum model is solved to obtain the
transient evolution of the $c^{\rm l}$ field in anti-plane segregation in
simple shear flow. The large-grain concentration field $c^{\rm l}$ is governed
by the following PDE:

$\frac{\partial c^{\rm l}}{\partial t}+\frac{\partial}{\partial
z}\left(-\hat{C}_{\rm diff}^{\rm anti}\bar{d}^{2}\dot{\gamma}\frac{\partial
c^{\rm l}}{\partial z}-\hat{C}^{\rm P,anti}_{\rm
seg}\frac{\bar{d}^{2}\dot{\gamma}}{P}c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha
c^{\rm l})\frac{\partial P}{\partial z}\right)=0,$ (5.2)

where $\hat{C}_{\rm diff}^{\rm anti}(I)$ is given by (4.3) and $\hat{C}^{\rm
P,anti}_{\rm seg}(I,J)$ is the dimensionless function depending on $I$ and $J$
given in (4.8), and $\bar{d}=c^{\rm l}d^{\rm l}+(1-c^{\rm l})d^{\rm s}$. As
discussed in Section 3, the pressure field is hydrostatic (varying linearly
with $z$), so that $\partial P/\partial z$ is constant, and to avoid a
singularity at the top free surface, a small constant $(1/4)(\partial
P/\partial z)\bar{d}_{0}$ is added to the pressure field. Therefore, we have
that $(1/P)\partial P/\partial z=1/(z+(1/4)\bar{d}_{0})$ in the last term of
(5.2), and moreover, the dimensionless quantity $J$ is
$J=(z+(1/4)\bar{d}_{0})/\bar{d}$. The pressure appearing in the inertial
number $I=\dot{\gamma}\bar{d}\sqrt{\rho_{\rm s}/P}$ is obtained from the
coarse-grained DEM data to control for the normal stress differences, and the
strain rate $\dot{\gamma}$ is prescribed and spatially uniform. Regarding the
boundary conditions, we impose no-flux boundary conditions at the top free
surface as well as on the bottom frictionless surface, i.e., $w^{\rm
l}_{z}=-\hat{C}_{\rm diff}^{\rm anti}\bar{d}^{2}\dot{\gamma}({\partial c^{\rm
l}}/{\partial z})-\hat{C}^{\rm P,anti}_{\rm
seg}({\bar{d}^{2}\dot{\gamma}}/{P})c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha
c^{\rm l})({\partial P}/{\partial z})=0$ at $z=0$ and $H$. For the initial
condition $c^{\rm l}_{0}(z)=c^{\rm l}(z,t=0)$, we use the coarse-grained field
from the initial DEM configuration for each case and use it as the initial
condition in the corresponding continuum simulation.

Figure 5: Comparisons of continuum model predictions with corresponding DEM
simulation results for the transient evolution of the segregation dynamics for
three cases of anti-plane segregation in simple shear flow of bidisperse
spheres: (a) Base case
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$; (b) Lower strain-rate case I
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.03,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$; (c) Lower strain-rate case II
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.015,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$. For each case, spatiotemporal
contours of the evolution of $c^{\rm l}$ measured in the DEM simulations are
shown on the left. Comparisons of the DEM simulations (solid black lines) and
continuum model predictions (dashed gray lines) of the $c^{\rm l}$ field are
shown on the right at four time snapshots representing different stages of the
segregation process, indicated by the vertical lines on the contour plots for
each case.

Next, we compare continuum model predictions of the transient evolution of the
segregation process against DEM data for different cases of anti-plane
segregation in simple shear flow. In the first set of comparisons, we vary the
applied strain rate and consider (1) the base case $\\{H/\bar{d}_{0}=30$,
$\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06$, $c^{\rm l}_{0}=0.50$, $d^{\rm
l}/d^{\rm s}=1.5\\}$, (2) lower strain-rate case I $\\{H/\bar{d}_{0}=30$,
$\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.0$3, $c^{\rm l}_{0}=0.50$, $d^{\rm
l}/d^{\rm s}=1.5\\}$, and (3) lower strain-rate case II $\\{H/\bar{d}_{0}=30$,
$\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.015$, $c^{\rm l}_{0}=0.50$, $d^{\rm
l}/d^{\rm s}=1.5\\}$, as shown in Figs. 5(a), (b) and (c), respectively.
Spatiotemporal contours of the evolution of the $c^{\rm l}$ field measured in
each DEM simulation are shown on the left, and comparisons between continuum
model predictions and DEM simulation data are plotted on the right for four
different snapshots in time corresponding to different stages of the
segregation process (short time, moderate time, long time, and steady-state),
as indicated by the vertical lines on the contour plots for each case. The
dimensionless times for each snapshot are given in the bottom right corner of
each plot in Fig. 5(a) and correspond to total shear strains of
$\tilde{t}=\dot{\gamma}t=60$, 320, 640, and 1300 in all three cases. The
steady-state concentration fields at $\tilde{t}=1300$ across all three cases
are nearly identical, which is why it was enforced that the steady-state flux
balance (4.5) should be independent of the strain rate. While the overall
segregation process proceeds faster with increasing strain rate, this
dependence is not simply linear. This may be seen by comparing the
concentration fields across all three cases at a fixed non-steady-state shear
strain. It is clear that the lowest strain-rate case (Fig. 5(c)) evolves
fastest in shear strain to its steady state. This is indicative of the
increased rate of diffusion (and hence segregation) at low inertial number
(Fig. 3(d)). Since this effect has been accounted for in the constitutive
equations for the anti-plane modes of diffusion (through (4.3)) and
segregation (through (4.8)), the continuum model predicts this characteristic
quite well, as can be seen in the transient comparisons.

Figure 6: Comparisons of continuum model predictions with corresponding DEM
simulation results for the transient evolution of the segregation dynamics for
three additional cases of anti-plane segregation in simple shear flow of
bidisperse spheres: (a) Higher large-grain concentration case
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.60,d^{\rm l}/d^{\rm s}=1.5\\}$; (b) Lower large-grain concentration
case I $\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.40,d^{\rm l}/d^{\rm s}=1.5\\}$; (c) Lower large-grain concentration
case II $\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.20,d^{\rm l}/d^{\rm s}=1.5\\}$. Results are organized as described in
the caption of Fig. 5. Figure 7: Comparisons of continuum model predictions
with corresponding DEM simulation results for the transient evolution of the
segregation dynamics for two additional cases of anti-plane segregation in
simple shear flow of bidisperse spheres: (a) Deeper layer case I
$\\{H/\bar{d}_{0}=40,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$; (b) Deeper layer case II
$\\{H/\bar{d}_{0}=60,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$. Results are organized as described in
the caption of Fig. 5.

Next, we vary the initial concentration of large grains and compare continuum
model predictions against DEM simulation data. We consider the following three
cases: (1) higher large-grain concentration case $\\{H/\bar{d}_{0}=30$,
$\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06$, $c^{\rm l}_{0}=0.60$, $d^{\rm
l}/d^{\rm s}=1.5\\}$, (2) lower large-grain concentration case I
$\\{H/\bar{d}_{0}=30$, $\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06$, $c^{\rm
l}_{0}=0.40$, $d^{\rm l}/d^{\rm s}=1.5\\}$, and (3) lower large-grain
concentration case II $\\{H/\bar{d}_{0}=30$,
$\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06$, $c^{\rm l}_{0}=0.20$, $d^{\rm
l}/d^{\rm s}=1.5\\}$, as shown in Figs. 6(a), (b), and (c), respectively.
Spatiotemporal contours of the evolution of the $c^{\rm l}$ field measured in
each DEM simulation are shown on the left, and comparisons of the continuum
model predictions against DEM simulation data are plotted on the right for
four time instants, as indicated by the vertical lines on the contour plots
for each case. Starting with the higher large-grain concentration case, the
steady-state transition zone between the large-grain region and the small-
grain region is pushed farther beneath the free surface, as can be seen on the
bottom right of Fig. 6(a). As the initial concentration of large grains is
reduced, the steady-state transition zone moves closer to the free surface, as
seen in the bottom right plots of Figs. 6(b) and (c). Finally, we vary the
depth of the layer and consider two additional cases: (1) deeper layer case I
$\\{H/\bar{d}_{0}=40$, $\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06$, $c^{\rm
l}_{0}=0.50$, $d^{\rm l}/d^{\rm s}=1.5\\}$ and (2) deeper layer case II
$\\{H/\bar{d}_{0}=60$, $\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06$, $c^{\rm
l}_{0}=0.50$, $d^{\rm l}/d^{\rm s}=1.5\\}$, as shown in Figs. 7(a) and (b).
Spatiotemporal contours of the evolution of the $c^{\rm l}$ field measured in
the DEM simulations are again shown on the left, and profiles of the continuum
model predictions compared against DEM simulation data at four time instants
are shown on the right. For the deeper layer cases, the steady-state
transition zone is farther beneath the free surface compared to the base case
of Fig. 5(a), and it takes more time for the segregation process to reach its
steady state. In fact, for deeper layer case II, the steady state is not
reached within the simulation time window considered in Fig. 7(b). In sum,
increasing the initial large-grain concentration or the thickness of the layer
both result in the steady-state transition zone being farther beneath the free
surface, and the continuum model can capture both the transient evolution of
the $c^{\rm l}$ field and its steady-state quite well across all of these
cases.

A collective view of Figs. 6 and 7 reveals that the width of the steady-state
transition region increases slightly with its depth beneath the free surface.
This is because the pressure-gradient-driven segregation flux attenuates with
pressure while the diffusion flux does not, so that the influence of the
diffusion flux is comparatively stronger than that of the segregation flux as
the depth beneath the free surface increases, leading to a wider transition
width. The $J$-dependence of the segregation flux in (4.8) is crucial to
capture this effect quantitatively and actually serves to make the segregation
flux less dependent on the pressure. Without $J$-dependence, the segregation
flux (3.2) scales as $P^{-1}$ for a fixed magnitude of the pressure gradient,
while with the $J$-dependent prefactor (4.8), the segregation flux scales as
$P^{-(1-n)}$, where we have determined that $n=0.58$ for frictional spheres
based on the steady-state DEM data of Fig. 4, so that the segregation flux
scales as $P^{-0.42}$. Neglecting $J$-dependence would lead to the segregation
flux attenuating too strongly as the pressure increases with depth beneath the
free surface, resulting in predicted transition widths that are too sharp near
the free surface and too wide deep beneath the free surface. This effect is
evident in the transient evolution of the $c^{\rm l}$ field predicted by the
$J$-independent, “in-plane” version of the continuum model in Fig. 2 for the
base case, in which the predicted transition width is too sharp at early times
when the transition region is near the free surface and too wide at late times
when the transition region has moved deeper beneath the free surface. In
conclusion, the proposed continuum model, which incorporates $I$\- and
$J$-dependence in the constitutive equations for the anti-plane fluxes, is
capable of quantitatively capturing the segregation dynamics in anti-plane
segregation in simple shear flow, including the rate at which the
concentration field approaches steady state and the width of the transition
region, across various cases subject to changes in strain rate, initial
conditions, and layer size.

##  6 Discussion

###  6.1 Diffusion flux based on mean square displacement data

In Section 4.1, we characterized the diffusion flux by starting with a fully
segregated granular mixture and examining the dynamics of mixing under
conditions in which diffusion is dominant. However, it is common in the
literature (e.g., Bancroft & Johnson, 2021) to characterize the diffusion flux
by calculating the mean square displacement (MSD) of the granular mixture. In
this section, we follow this process to characterize the diffusion flux and
assess whether the resulting flux constitutive equations can capture the
dynamics of the $c^{\rm l}$ field in the anti-plane mode of diffusion and
segregation.

To this end, we consider both the in-plane and anti-plane modes of diffusion
in DEM simulations of simple shear. In our prior work (Liu et al., 2023), the
in-plane MSD was calculated from DEM data for steady, simple shear flows of
the same dense granular system (Appendix A) to determine the binary diffusion
coefficient $D$ for both bidisperse mixtures of spheres and for the
monodisperse case. We follow an analogous procedure here to calculate the
diffusion coefficient in the anti-plane mode. DEM simulations of simple shear
flow are performed as described in Section 4.1 for a well-mixed bidisperse
granular system with $d^{\rm l}/d^{\rm s}=1.5$ as well as for the monodisperse
case. The MSD of a system of $N$ particles along the anti-plane direction
(i.e., the $z$-direction) is calculated as a function of time (e.g., Natarajan
et al., 1995; Campbell, 1997; Utter & Behringer, 2004; Fry et al., 2019;
Bancroft & Johnson, 2021) according to

$\text{MSD}(t)=\dfrac{1}{N}\sum_{n=1}^{N}(z_{n}(t)-z_{n}(0))^{2}=2Dt,$ (6.1)

where $z_{n}(t)$ is the $z$-coordinate of the $n$th grain at time $t$. The
binary diffusion coefficient $D$ may then be inferred from the slope of the
MSD in time. To avoid boundary effects in the calculation of the MSD, grains
that are initially within $10\bar{d}_{0}$ of either the top or bottom boundary
are excluded from the system of $N$ particles used to calculate the MSD,
leaving a set of $N\approx 3000$ grains. The normalized binary diffusion
coefficient $D$ is plotted versus $\dot{\gamma}\bar{d}^{2}$ (with both
quantities normalized by $d^{\rm s}\sqrt{P_{\rm top}/\rho_{\rm s}}$) for
different strain rates in Fig. 8(a) for both modes of diffusion and both
bidisperse and monodisperse granular systems. The DEM data for the in-plane
mode of diffusion in Fig. 8(a) is the data reported in Fig. 3(b) of Liu et al.
(2023). The DEM data for the binary diffusion coefficient collapses to a
nearly linear relation with $D\sim\dot{\gamma}\bar{d}^{2}$ across the range of
strain rates considered and for both monodisperse and bidisperse systems in
in-plane as well as anti-plane modes of diffusion. The solid black line in
Fig. 8(a) is a linear fit that estimates the diffusion parameter as $C_{\rm
diff}=0.048$, which is very close to the estimated value of $C_{\rm diff}^{\rm
in}=0.045$ determined in Liu et al. (2023).

Figure 8: (a) The binary diffusion coefficient $D$, calculated using the mean
square displacement, plotted versus $\dot{\gamma}\bar{d}^{2}$ in homogeneous,
steady simple shear DEM simulations. Each symbol represents $D$ calculated
from a distinct DEM simulation. Both axes are normalized by $d^{\rm
s}\sqrt{P_{\rm top}/\rho_{\rm s}}$. The solid line is the best fit of a linear
relation with a slope of $C_{\rm diff}=0.048$. Comparisons between continuum
model predictions (dashed gray lines) and DEM simulation results (solid black
lines) of the $c^{\rm l}$ field at four time snapshots using a constitutive
equation for the diffusion flux based on mean square displacement data for (b)
the base case of anti-plane segregation in simple shear flow
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$ and (c) the lower strain-rate case II
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.015,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$.

The diffusion parameter $C_{\rm diff}$ estimated via the MSD approach in Fig.
8(a) is approximately rate-independent, i.e, it is only weakly dependent upon
the strain rate (or inertial number), which contradicts the previously
observed $I$-dependence at the continuum scale discussed in Section 4.1 and
illustrated in Fig. 3(d). Therefore, we test this constant value of $C_{\rm
diff}=0.048$ in the segregation dynamics equation (5.2). In doing so,
$I$-dependence is neglected in both the diffusion and segregation parameters,
so that the segregation dynamics equation becomes

$\frac{\partial c^{\rm l}}{\partial t}+\frac{\partial}{\partial
z}\left(-{C}_{\rm diff}\bar{d}^{2}\dot{\gamma}\frac{\partial c^{\rm
l}}{\partial z}-CJ^{n}\frac{\bar{d}^{2}\dot{\gamma}}{P}c^{\rm l}(1-c^{\rm
l})(1-\alpha+\alpha c^{\rm l})\frac{\partial P}{\partial z}\right)=0,$ (6.2)

with the reduced parameter set $\\{C_{\rm
diff}=0.048,C=0.22,n=0.58,\alpha=0.4\\}$. This approach maintains the steady-
state collapse of Fig. 4. We then numerically calculate the evolution of the
large-grain concentration field using (6.2) accompanied by the same boundary
and initial conditions as in the preceding sections. First, we compare
predictions of the continuum model with DEM simulation data for the base case
of anti-plane segregation in simple shear flow (Fig. 5(a)) as shown in Fig.
8(b). The continuum model predictions using the constant diffusion parameter
capture the DEM simulation results reasonably well for this case. However, for
the lower strain-rate case II (Fig. 5(c)), as shown in Fig. 8(c), continuum
model predictions using the constant diffusion parameter lag behind the DEM
data in time. This indicates that the diffusion parameter is indeed rate-
dependent, and at low strain rates, its value is higher than what is inferred
from MSD data. This is consistent with the estimates of the diffusion
parameter shown in Fig. 3(d). We leave this discrepancy between the MSD-based
estimate and the continuum-model-fitting-based estimate of the diffusion
parameter as an open problem for future work. However, we have established in
this study that the anti-plane diffusion parameter is rate-dependent in
contrast to the in-plane diffusion parameter, which is nearly rate-
independent. This implies that the diffusion process in general is mode-
dependent and to account for the mode-dependence, we propose a strategy for
generalizing the diffusion parameter in the next section.

###  6.2 Generalized three-dimensional constitutive relations

In Section 4, we proposed constitutive equations for the anti-plane modes of
the diffusion flux and the pressure-gradient-driven segregation flux. Both the
anti-plane diffusion and segregation fluxes take different forms from the in-
plane modes studied in our prior works (Liu et al., 2023; Singh et al., 2023)
and described in Section 3.2. In this section, we propose a strategy for
generalizing the constitutive forms for both fluxes, which takes into account
both the in-plane and the anti-plane modes of diffusion and pressure-gradient-
driven segregation.

In order to propose a generalized constitutive form for the diffusion flux, we
denote the unit vector along the direction of the concentration gradient as
$n^{\rm c}_{i}=({\partial c^{\rm l}/\partial x_{i}})/{|\partial c^{\rm
l}/\partial{\bf x}|}$ and the flow direction tensor as $N_{ij}={D_{ij}}/{|{\bf
D}|}$. Next, we introduce a combined scalar invariant $K^{\rm c}=n^{\rm
c}_{i}N_{ij}N_{jk}n^{\rm c}_{k}$, which is used to quantify whether diffusion
is in its in-plane or anti-plane mode. Straightforward calculations described
in Appendix C show that for the in-plane mode of diffusion, the scalar
invariant is $K^{\rm c}=1/2$, and for the anti-plane mode of diffusion, the
scalar invariant is $K^{\rm c}=0$. This feature indicates that $K^{\rm c}$ is
a suitable invariant for capturing mode-dependence, and we define a
generalized diffusion parameter through a linear combination of $\hat{C}_{\rm
diff}^{\rm anti}$ and $C_{\rm diff}^{\rm in}$ based on the invariant $K^{\rm
c}$ as follows:

$\hat{C}_{\rm diff}(I,K^{\rm c})=\hat{C}^{\rm anti}_{\rm diff}(I)+2K^{\rm
c}\left(C^{\rm in}_{\rm diff}-\hat{C}^{\rm anti}_{\rm diff}(I)\right),$ (6.3)

where $C^{\rm in}_{\rm diff}$ is the constant diffusion parameter for the in-
plane mode of diffusion determined in Liu et al. (2023) for frictional
spheres. In (6.3), $\hat{C}_{\rm diff}=C^{\rm in}_{\rm diff}$ for the in-plane
mode when $K^{\rm c}=1/2$, and $\hat{C}_{\rm diff}=\hat{C}^{\rm anti}_{\rm
diff}$ for the anti-plane mode when $K^{\rm c}=0$. The generalized diffusion
flux may then be written as $w^{\rm diff}_{i}=-\hat{C}_{\rm diff}(I,K^{\rm
c})\bar{d}^{2}\dot{\gamma}({\partial c^{\rm l}}/{\partial x_{i}})$.

Similarly, in order to propose a generalized constitutive equation for the
pressure-gradient-driven segregation flux, we denote the unit vector along the
direction of the pressure gradient as $n^{\rm P}_{i}=({\partial P/\partial
x_{i}})/{|\partial P/\partial{\bf x}|}$, and introduce another combined scalar
invariant $K^{\rm P}=n^{\rm P}_{i}N_{ij}N_{jk}n^{\rm P}_{k}$. The generalized
pressure-gradient segregation parameter is then defined through a linear
combination of $\hat{C}^{\rm P,anti}_{\rm seg}$ and $C^{\rm P,in}_{\rm seg}$
as

$\hat{C}^{\rm P}_{\rm seg}(I,J,K^{\rm P})=\hat{C}_{\rm seg}^{\rm
P,anti}(I,J)+2K^{\rm P}\left(C_{\rm seg}^{\rm P,in}-\hat{C}_{\rm seg}^{\rm
P,anti}(I,J)\right),$ (6.4)

where $C_{\rm seg}^{\rm P,in}$ is the constant in-plane segregation parameter,
so that the generalized pressure-gradient-driven segregation flux is

$w^{\rm P}_{i}=-\hat{C}^{\rm P}_{\rm seg}(I,J,K^{\rm
P})\frac{\bar{d}^{2}\dot{\gamma}}{P}c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha
c^{\rm l})\frac{\partial P}{\partial{x}_{i}}.$ (6.5)

Analogous to the discussion of the preceding paragraph, for the in-plane mode
of segregation, $K^{\rm P}=1/2$ and $\hat{C}^{\rm P}_{\rm seg}=C_{\rm
seg}^{\rm P,in}$, and for the anti-plane mode of segregation, $K^{\rm P}=0$
and $\hat{C}^{\rm P}_{\rm seg}=\hat{C}_{\rm seg}^{\rm P,anti}$.

Putting the generalized flux constitutive equations together, the generalized
form of the segregation dynamics equation can be written as

$\frac{Dc^{\rm l}}{Dt}+\frac{\partial}{\partial x_{i}}\left(-\hat{C}_{\rm
diff}(I,K^{\rm c})\bar{d}^{2}\dot{\gamma}\frac{\partial c^{\rm l}}{\partial
x_{i}}-\hat{C}^{\rm P}_{\rm seg}(I,J,K^{\rm
P})\frac{\bar{d}^{2}\dot{\gamma}}{P}c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha
c^{\rm l})\frac{\partial P}{\partial x_{i}}\right)=0,$ (6.6)

where $\hat{C}_{\rm diff}$ and $\hat{C}^{\rm P}_{\rm seg}$ are given by (6.3)
and (6.4), respectively. The diffusion parameter $\hat{C}_{\rm diff}$ is
described by the parameter set $\\{C_{\rm diff}^{\rm
in},C_{0},C_{\infty},I_{\rm tr},m\\}$, and the segregation parameter
$\hat{C}^{\rm P}_{\rm seg}$ is described by the parameter set $\\{C^{\rm
P,in}_{\rm seg},C,n\\}$. In the limits of pure in-plane and pure anti-plane
modes, the generalized segregation dynamics equation (6.6) will reduce to the
appropriate forms that have been tested in their respective modes, and future
work will test the generalized model (6.6) in more complex flow configurations
that involve mixed-mode diffusion and segregation, such as split-bottom flow
(Hill & Fan, 2008; Fan & Hill, 2010).

##  7 Concluding remarks

In this paper, we have studied diffusion and size segregation in their anti-
plane modes in dense, bidisperse granular mixtures of spheres. In DEM
simulations, we observed that the diffusion and pressure-gradient-driven
segregation fluxes are mode-dependent in nature, and the previously-developed,
in-plane size-segregation model of Singh et al. (2023) is unable to capture
the segregation dynamics in the anti-plane mode of segregation. Guided by DEM
simulations, we proposed phenomenological constitutive equations for the anti-
plane diffusion and segregation fluxes. We found that the anti-plane diffusion
parameter $\hat{C}_{\rm diff}^{\rm anti}$ is rate dependent and depends on the
inertial number $I$, and we proposed a phenomenological fitting function (4.3)
to capture the $I$-dependence. Moreover, we found that the pressure-gradient-
driven segregation flux depends on the magnitude of the pressure gradient in a
nonlinear fashion, and we proposed a constitutive equation for the anti-plane
segregation parameter $\hat{C}^{\rm P,anti}_{\rm seg}$ (4.8) that captures
this dependence using an additional dimensionless parameter
$J={P}/{(\bar{d}|\partial P/\partial\textbf{x}|)}$. The newly-developed
phenomenological constitutive equations predict the transient evolution of the
segregation dynamics quite well across different cases of anti-plane
segregation in simple shear flow as the strain rate, initial conditions, and
layer size are varied. Finally, we have proposed a strategy for generalizing
the constitutive equations for the diffusion and pressure-gradient-driven
segregation fluxes that synthesizes the anti-plane flux constitutive equations
proposed in this work with the in-plane flux constitutive equations of Singh
et al. (2023), which is suitable for more complex, three-dimensional flows.

Although the proposed continuum model performs well in predicting the
segregation dynamics in anti-plane segregation in simple shear flow, several
important questions remain. First, the present work has focused on a single
grain-size ratio of $d^{\rm l}/d^{\rm s}=1.5$. It has been established that
the in-plane pressure-gradient-driven segregation flux should depend upon the
grain-size ratio (e.g., Tunuguntla et al., 2014; Schlick et al., 2015; Jones
et al., 2018; Trewhela et al., 2021), and it remains to probe this dependence
in the anti-plane mode. Second, future work is required to investigate the
discrepancy between the MSD-based estimate and the continuum-model-fitting
based estimate of the diffusion parameter in the anti-plane mode of diffusion,
as discussed in Section 6.1. Finally, the generalized constitutive equations
that account for both in-plane and anti-plane modes of diffusion and
segregation still need to be tested in more complex flow configurations—for
example, annular shear flow with gravity, split bottom flow (Hill & Fan, 2008;
Fan & Hill, 2010), or flow in a blade mixer (Yang et al., 2022). To achieve
this, it is necessary to integrate the proposed generalized segregation model
with rheological constitutive equations (e.g., the nonlocal granular fluidity
model (Kamrin & Koval, 2012; Henann & Kamrin, 2013)) to develop a coupled
continuum model that is able to predict the flow and segregation dynamics
simultaneously and to develop a robust numerical formulation to solve the
coupled system of continuum equations in complex geometries. These points will
be pursued in future work.

## Acknowledgements

This work was supported by funds from NSF-CBET-1552556.

##  Appendix A Simulated granular system and averaging methods

We consider granular systems consisting of dense, bidisperse mixtures of
spheres. The mean large-grain diameter is $d^{\rm l}=3\,{\rm mm}$, and the
mean small-grain diameter is $d^{\rm s}=2\,{\rm mm}$, so that the grain-size
ratio $d^{\rm l}/d^{\rm s}=1.5$ is held constant throughout. We impose a
polydispersity of $\pm 10\%$ to the mean diameter of each species to prevent
crystallization. The grain-material density is $\rho_{\rm s}=2450\,{\rm
kg/m}^{3}$. The grain-grain interaction force is given through a contact law
that accounts for Hookean elasticity, damping, and sliding friction (Silbert
et al., 2001). The relevant interaction properties are the normal contact
stiffness $k_{\rm n}$, the tangential contact stiffness $k_{\rm t}$, the
coefficient of restitution $e$, and the inter-particle friction coefficient
$\mu_{\rm surf}$. We restrict our attention to the nearly-rigid particle
regime and set the normal contact stiffness to be sufficiently high compared
to the confining pressure throughout, i.e., $k_{\rm n}/P\bar{d}_{0}>10^{4}$.
The other parameters are set as $k_{\rm t}/k_{\rm n}=1/2$, $e=0.1$, and
$\mu_{\rm surf}=0.4$ throughout. Numerical integration of the equations of
motion for each grain is performed using the open-source software LAMMPS
(Plimpton, 1995), and the time step of numerical integration is specified to
sufficiently small to ensure accuracy and stability of simulation results.

To extract coarse-grained concentration and pressure fields from DEM
simulations, we use a bin-based approach, as described in our prior work (Liu
et al., 2023). Briefly, this approach utilizes rectangular-cubiodal bins of
width $\Delta$ that are centered at a $z$-position and span the simulation
domain along the $x$\- and $y$-directions. For a snapshot at time $t$, each
grain $i$ intersected by the bin is given a weight $V_{i}$, which is equal to
the volume of grain $i$ inside the bin. The sets of large and small grains
intersected by the bin are denoted as ${\cal F}^{\rm l}$ and ${\cal F}^{\rm
s}$, respectively (Tunuguntla et al., 2016). The instantaneous solid volume
fractions for large and small grains are $\phi^{\rm l}(z,t)=(\sum_{i\in{\cal
F}^{\rm l}}V_{i})/V$ and $\phi^{\rm s}(z,t)=(\sum_{i\in{\cal F}^{\rm
s}}V_{i})/V$, respectively, where $V$ is the volume of the bin. The
concentration field for the large grains is calculated as $c^{\rm
l}(z,t)=\phi^{\rm l}(z,t)/(\phi^{\rm l}(z,t)+\phi^{\rm s}(z,t))$. The
instantaneous stress tensor associated with grain $i$ is
$\bm{\sigma}_{i}(t)=(\sum_{j\neq i}{\bf r}_{ij}\otimes{\bf f}_{ij})/(\pi
d_{i}^{3}/6)$, where ${\bf r}_{ij}$ is the position vector from the center of
grain $i$ to the center of grain $j$, ${\bf f}_{ij}$ is the contact force
applied on grain $i$ by grain $j$, and $d_{i}$ is the diameter of grain $i$.
The instantaneous stress field is $\bm{\sigma}(z,t)=(\sum_{i\in{\cal
F}}V_{i}\bm{\sigma}_{i}(t))/V$, where ${\cal F}={\cal F}^{\rm l}\cup{\cal
F}^{\rm s}$ is the set of all grains intersected by the bin, and the
instantaneous pressure field is $P(z,t)=-(1/3){\rm tr}(\bm{\sigma}(z,t))$.
Throughout, we take a bin width of $\Delta=4\bar{d}_{0}$ and a spatial
resolution of about $0.2\bar{d}_{0}$. For the steady-state collapse of Fig. 4,
we consider DEM data from the spatial region where $c^{\rm l}\in(0.2,0.8)$.
When calculating the quantities appearing in the collapse of Fig. 4, the
instantaneous concentration and pressure fields are spatially smoothed and
differentiated using Lucy’s quartic kernel function (Tunuguntla et al., 2016)
with a kernel width of $4\bar{d}_{0}$. To obtain the relevant steady-state
field quantities in (4.7), we consider a time window within the steady-state
regime and generate $N=1000$ snapshots of the instantaneous, smoothed fields
and arithmetically average these fields over all snapshots to obtain fields
that only depend on the spatial coordinate.

##  Appendix B Alternative constitutive equation for the segregation flux

The works of Trewhela et al. (2021) and Barker et al. (2021) suggest a
constitutive form for the pressure-gradient-driven segregation flux that
depends on $J$ in an alternative manner, adapted to the present work as
follows:

$w^{\rm
P}_{i}=\dfrac{C_{1}}{1+C_{2}/J}\left[1-\left(\frac{C_{\infty}-C_{0}}{C_{\infty}+C_{0}}\right)\dfrac{(I_{\rm
tr}/I)^{m}-1}{(I_{\rm
tr}/I)^{m}+1}\right]\dfrac{\bar{d}^{2}\dot{\gamma}}{P}c^{\rm l}(1-c^{\rm
l})(1-\alpha+\alpha c^{\rm l})\dfrac{\partial P}{\partial x_{i}},$ (B.1)

where the flux depends on $J$ through the factor
$\tilde{C}(J)=C_{1}/(1+C_{2}/J)$ rather than $\tilde{C}(J)=CJ^{n}$, and
$\\{C_{1},C_{2}\\}$ are constant fitting parameters. The prefactor
$(C_{1}/(1+C_{2}/J))(\bar{d}^{2}\dot{\gamma}/P)$ in (B.1) may be rewritten as
$C_{1}\bar{d}^{2}\dot{\gamma}/(P+C_{2}\bar{d}|\partial P/\partial{\bf x}|)$,
and the term in the denominator involving $C_{2}\bar{d}|\partial
P/\partial{\bf x}|$ was included by Trewhela et al. (2021) and used by Barker
et al. (2021) to prevent a singularity when the pressure equals zero, such as
at a free surface. (The flux constitutive equation using the factor $CJ^{n}$
does exhibit a singularity when the pressure equals zero, and in Section 5, we
circumvented this issue by adding a small top pressure equal to the weight of
a layer of $(1/4)\bar{d}_{0}$ thickness, i.e., $(1/4)(\partial P/\partial
z)\bar{d}_{0}$, to the pressure field.) Here, we consider the alternative
segregation flux constitutive equation (B.1) and assess its ability to capture
segregation dynamics in the anti-plane mode of segregation.

With $\tilde{C}(J)=C_{1}/(1+C_{2}/J)$, the parameters $\\{C_{1},C_{2}\\}$ may
be determined by fitting to the steady-state DEM data from 13 different cases
of anti-plane segregation in simple shear flow, and the best fit is shown by
the dash-dotted black line in Fig. 4. The estimated parameters are
$\\{C_{1}=3.2,C_{2}=28.62\\}$. We note that the value of $C_{2}=28.62$ is
quite large compared to the value of 0.2712 used by Trewhela et al. (2021) for
the purpose of preventing a singularity, so in the anti-plane context, this
term plays a key role in capturing the dynamics of segregation and is not
simply included to regularize the free-surface singularity. Then, we plug the
alternative constitutive form for the pressure-gradient-driven segregation
flux in the segregation dynamics equation (5.2) and compare numerical
predictions of the dynamics of the $c^{\rm l}$ field against the DEM
simulation results. We consider four different cases: (1) the base case, (2)
lower strain-rate case II, (3) lower large-grain concentration case II, and
(4) deeper layer case II, as shown in Figs. 9(a), (b), (c), and (d),
respectively. The predictions of the alternative constitutive form are shown
using dash-dotted gray lines and are as capable of capturing the segregation
dynamics as predictions of the prior constitutive form, which are shown using
dashed gray lines.

Figure 9: Comparisons of continuum model predictions with corresponding DEM
simulation results for the transient evolution of the segregation dynamics for
four different cases of anti-plane segregation in simple shear flow of
bidisperse spheres: (a) Base case
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/g}=0.06,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$; (b) Lower strain-rate case II
$\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/g}=0.015,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$; (c) Lower large-grain concentration
case II $\\{H/\bar{d}_{0}=30,\dot{\gamma}\sqrt{\bar{d}_{0}/g}=0.06,c^{\rm
l}_{0}=0.20,d^{\rm l}/d^{\rm s}=1.5\\}$; and (d) Deeper layer case II
$\\{H/\bar{d}_{0}=60,\dot{\gamma}\sqrt{\bar{d}_{0}/G}=0.06,c^{\rm
l}_{0}=0.50,d^{\rm l}/d^{\rm s}=1.5\\}$. DEM simulation results are shown
using solid black lines. Theory-I, shown using dashed gray lines, corresponds
to predictions of the continuum model using the prior constitutive form for
the anti-plane segregation flux (4.8). Theory-II, shown using dash-dotted gray
lines, corresponds to predictions of the continuum model using the alternative
constitutive form (B.1).

##  Appendix C Generalized diffusion and segregation parameters

In this appendix, we give additional context to the combined invariants used
in Section 6.2 to obtain the generalized diffusion parameter (6.3) and the
generalized segregation parameter (6.4). Here, we focus the discussion on the
constitutive equation for the generalized diffusion parameter, but an
analogous process may be applied to the generalized segregation parameter.
Between our prior work (Liu et al., 2023) and Section 4.1 of the present work,
we have observed that the diffusion process depends on the mode of diffusion
relative to the plane of shearing (in-plane or anti-plane) and the inertial
number $I$. To account for this dependence in a constitutive equation, we
consider the flow direction tensor ${\bf N}={\bf D}/|{\bf D}|$, which is
deviatoric ${\rm tr}({\bf N})=0$ and of unit magnitude $|{\bf N}|=1$, and the
unit vector along the concentration gradient ${\bf n}^{\rm c}=(\partial c^{\rm
l}/\partial{\bf x})/|\partial c^{\rm l}/\partial{\bf x}|$. Assuming that the
diffusion flux and the concentration gradient act in opposite directions along
the same line of action, a generalized constitutive equation for the diffusion
flux may then be written in terms of ${\bf N}$ and ${\bf n}^{\rm c}$ as
follows:

$w^{\rm diff}_{i}=-\check{C}_{\rm diff}\left(I,{\bf N},{\bf n}^{\rm
c}\right)\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm l}}{\partial x_{i}},$
(C.1)

where $\check{C}_{\rm diff}\left(I,{\bf N},{\bf n}^{\rm c}\right)$ is a scalar
function that accounts for the mode-dependent nature of diffusion. The
function $\check{C}_{\rm diff}\left(I,{\bf N},{\bf n}^{\rm c}\right)$ must be
an isotropic function of its arguments, so it may be written in terms of the
combined invariants associated with the unit flow direction tensor ${\bf N}$
and the unit vector ${\bf n}^{\rm c}$ as follows:

$\begin{split}&K_{1}={\rm tr}(\textbf{N}),\,\,\,K_{2}={\rm
tr}(\textbf{N}^{2}),\,\,\,K_{3}={\rm tr}(\textbf{N}^{3}),\\\
&K_{4}=\textbf{n}^{\rm c}\cdot\textbf{n}^{\rm c},\,\,\,K_{5}=\textbf{n}^{\rm
c}\cdot\textbf{N}\textbf{n}^{\rm c},\,\,\,K_{6}=\textbf{n}^{\rm
c}\cdot\textbf{N}^{2}\textbf{n}^{\rm c},\end{split}$ (C.2)

where the invariants $K_{1}=0,\,K_{2}=1$ and $K_{4}=1$ are constants due to
the properties of ${\bf N}$ and ${\bf n}^{\rm c}$, and the remaining invariant
set reduces to $\\{K_{3},K_{5},K_{6}\\}$. (The invariants are denoted as
$\\{K_{1},K_{2},K_{3},K_{4},K_{5},K_{6}\\}$ to avoid confusion with the
inertial number $I$ and the dimensionless quantity $J$ introduced in Section
4.2.) An isotropic representation of $\check{C}_{\rm diff}\left(I,{\bf N},{\bf
n}^{\rm c}\right)$ may then be written as $\hat{C}_{\rm
diff}\left(I,K_{3},K_{5},K_{6}\right)$, and the flux constitutive equation
(C.1) reduces to

$w^{\rm diff}_{i}=-\hat{C}_{\rm
diff}\left(I,K_{3},K_{5},K_{6}\right)\bar{d}^{2}\dot{\gamma}\dfrac{\partial
c^{\rm l}}{\partial x_{i}}.$ (C.3)

To further reduce (C.3), we consider a situation, in which the concentration
field only varies along the $z$-direction, so that $[{\bf n}^{\rm
c}]=\pm[0\,\,\,\,0\,\,\,\,1]^{\top}$, and the velocity field corresponds to a
steady, viscometric flow of the following form: $v_{x}(y,z)$, $v_{y}=0$, and
$v_{z}=0$. The velocity gradient tensor ${\bf L}={\rm grad}\,{\bf v}$ and the
strain-rate tensor ${\bf D}=(1/2)({\bf L}+{\bf L}^{\top})$ are then

$[\textbf{L}]=\begin{bmatrix}0&\partial v_{x}/\partial y&\partial
v_{x}/\partial z\\\ 0&0&0\\\
0&0&0\end{bmatrix}\quad\text{and}\quad[\textbf{D}]=\frac{1}{2}\begin{bmatrix}0&\partial
v_{x}/\partial y&\partial v_{x}/\partial z\\\ \partial v_{x}/\partial y&0&0\\\
\partial v_{x}/\partial z&0&0\end{bmatrix}.$ (C.4)

This steady, viscometric flow field involves shear within two planes: the
$x$-$y$-plane and the $x$-$z$-plane. The concentration gradient lies within
the $x$-$z$-plane, so we denote $\partial v_{x}/\partial z=\dot{\gamma}_{\rm
in}$. Similarly, the concentration gradient is perpendicular to the
$x$-$y$-plane, so we denote $\partial v_{x}/\partial y=\dot{\gamma}_{\rm
anti}$. The equivalent shear strain rate is $\dot{\gamma}=\sqrt{2}|{\bf
D}|=\sqrt{(\partial v_{x}/\partial y)^{2}+(\partial v_{x}/\partial
z)^{2}}=\sqrt{\dot{\gamma}_{\rm anti}^{2}+\dot{\gamma}_{\rm in}^{2}}$. The
unit flow direction tensor ${\bf N}={\bf D}/|{\bf D}|$ may then be written as

$[{\bf N}]=\frac{\sqrt{2}[{\bf
D}]}{\dot{\gamma}}=\frac{1}{\sqrt{2}}\begin{bmatrix}0&\dot{\gamma}_{\rm
anti}/\dot{\gamma}&\dot{\gamma}_{\rm in}/\dot{\gamma}\\\ \dot{\gamma}_{\rm
anti}/\dot{\gamma}&0&0\\\ \dot{\gamma}_{\rm
in}/\dot{\gamma}&0&0\end{bmatrix},$ (C.5)

and the tensor ${\bf N}^{2}$ becomes

$[{\bf N}^{2}]=\frac{1}{2}\begin{bmatrix}\frac{\dot{\gamma}_{\rm
anti}^{2}+\dot{\gamma}_{\rm in}^{2}}{\dot{\gamma}^{2}}&0&0\\\
0&\frac{\dot{\gamma}_{\rm anti}^{2}}{\dot{\gamma}^{2}}&\frac{\dot{\gamma}_{\rm
anti}\dot{\gamma}_{\rm in}}{\dot{\gamma}^{2}}\\\ 0&\frac{\dot{\gamma}_{\rm
anti}\dot{\gamma}_{\rm in}}{\dot{\gamma}^{2}}&\frac{\dot{\gamma}_{\rm
in}^{2}}{\dot{\gamma}^{2}}\end{bmatrix},$ (C.6)

so the reduced set of invariants is $K_{3}={\rm tr}(\textbf{N}^{3})=0$,
$K_{5}={\bf n}^{\rm c}\cdot{\bf N}{\bf n}^{\rm c}=0$, and $K_{6}={\bf n}^{\rm
c}\cdot{\bf N}^{2}{\bf n}^{\rm c}=\dot{\gamma}_{\rm
in}^{2}/2\dot{\gamma}^{2}$. Therefore, the invariants $K_{3}$ and $K_{5}$ are
always zero in this steady, viscometric flow scenario and cannot differentiate
between the in-plane and anti-plane modes of diffusion, so we neglect the
dependence of the scalar function $\hat{C}_{\rm diff}$ on $K_{3}$ and $K_{5}$.
The remaining invariant $K_{6}$ is bounded within the range $K_{6}\in[0,1/2]$.
In the in-plane mode of diffusion, $\dot{\gamma}_{\rm anti}=0$,
$\dot{\gamma}=|\dot{\gamma}_{\rm in}|$, and $K_{6}=1/2$, and in the anti-plane
mode of diffusion, $\dot{\gamma}_{\rm in}=0$, $\dot{\gamma}=|\dot{\gamma}_{\rm
anti}|$, and $K_{6}=0$, so $K_{6}$ is a suitable invariant for capturing the
mode-dependence of the diffusion flux. Renaming this invariant $K_{6}={\bf
n}^{\rm c}\cdot{\bf N}^{2}{\bf n}^{\rm c}=K^{\rm c}$, the flux equation (C.3)
may be written as

$w^{\rm diff}_{i}=-\hat{C}_{\rm diff}(I,K^{\rm
c})\bar{d}^{2}\dot{\gamma}\frac{\partial c^{\rm l}}{\partial x_{i}},$ (C.7)

and when a simple linear dependence of $\hat{C}_{\rm diff}$ on $K^{\rm c}$ is
adopted, we recover (6.3) from Section 6.2.

Following an analogous set of arguments, we denote the unit vector along the
pressure gradient as ${\bf n}^{\rm P}=(\partial P/\partial{\bf x})/|\partial
P/\partial{\bf x}|$ and introduce the combined scalar invariant $K^{\rm
P}={\bf n}^{\rm P}\cdot{\bf N}^{2}{\bf n}^{\rm P}$. We then write the
generalized form for the pressure-gradient-driven segregation flux as follows:

$w^{\rm P}_{i}=-\hat{C}^{\rm P}_{\rm seg}(I,J,K^{\rm
P})\frac{\bar{d}^{2}\dot{\gamma}}{P}c^{\rm l}(1-c^{\rm l})(1-\alpha+\alpha
c^{\rm l})\frac{\partial P}{\partial x_{i}},$ (C.8)

where $J=P/(\bar{d}|\partial P/\partial{\bf x}|)$. In the in-plane mode of
segregation, the scalar invariant is $K^{\rm P}=1/2$, and $K^{\rm P}=0$ in the
anti-plane mode of segregation. Adopting a linear dependence of $\hat{C}^{\rm
P}_{\rm seg}$ on $K^{\rm P}$ gives (6.4) from Section 6.2.

## References

  * Bancroft & Johnson (2021) Bancroft, R.S.J. & Johnson, C.G. 2021 Drag, diffusion and segregation in inertial granular flows. J. Fluid Mech. 924, A3.
  * Barker et al. (2021) Barker, T., Rauter, M., Maguire, E.S.F., Johnson, C.G. & Gray, J.M.N.T. 2021 Coupling rheology and segregation in granular flows. J. Fluid Mech. 909, A22.
  * Campbell (1997) Campbell, C.S. 1997 Self-diffusion in granular shear flows. J. Fluid Mech. 348, 85–101.
  * Depken et al. (2007) Depken, M., Lechman, J. B., van Hecke, M., van Saarloos, W. & Grest, G. S. 2007 Stresses in smooth flows of dense granular media. Europhys. Lett. 78, 58001.
  * Fan & Hill (2010) Fan, Y. & Hill, K.M. 2010 Shear-driven segregation of dense granular mixtures in a split-bottom cell. Phys. Rev. E 81, 041303.
  * Fan & Hill (2011) Fan, Y. & Hill, K.M. 2011 Theory for shear-induced segregation of dense granular mixtures. New J. Phys. 13, 095009\.
  * Fan et al. (2014) Fan, Y., Schlick, C.P., Umbanhowar, P.B., Ottino, J.M. & Lueptow, R.M. 2014 Modelling size segregation of granular materials: the roles of segregation, advection and diffusion. J. Fluid Mech. 741, 252–279.
  * Fry et al. (2018) Fry, A.M., Umbanhowar, P.B., Ottino, J.M. & Lueptow, R.M. 2018 Effect of pressure on segregation in granular shear flows. Phys. Rev. E 97, 062906.
  * Fry et al. (2019) Fry, A.M., Umbanhowar, P.B., Ottino, J.M. & Lueptow, R.M. 2019 Diffusion, mixing, and segregation in confined granular flows. AIChE J. 65, 875–881.
  * Gajjar & Gray (2014) Gajjar, P. & Gray, J.M.N.T. 2014 Asymmetric flux models for particle-size segregation in granular avalanches. J. Fluid Mech. 757, 297–329.
  * Golick & Daniels (2009) Golick, L.A. & Daniels, K.E. 2009 Mixing and segregation rates in sheared granular materials. Phys. Rev. E 80, 042301.
  * Gray (2018) Gray, J.M.N.T. 2018 Particle segregation in dense granular flows. Annu. Rev. Fluid Mech. 50, 407–433.
  * Gray & Chugunov (2006) Gray, J.M.N.T. & Chugunov, V.A. 2006 Particle-size segregation and diffusive remixing in shallow granular avalanches. J. Fluid Mech. 569, 365–398.
  * Gray & Thornton (2005) Gray, J.M.N.T. & Thornton, A.R. 2005 A theory for particle size segregation in shallow granular free-surface flows. Proc. R. Soc. A 461, 1447–1473.
  * Henann & Kamrin (2013) Henann, D.L. & Kamrin, K. 2013 A predictive, size-dependent continuum model for dense granular flows. P. Natl. Acad. Sci. USA 110, 6730–6735.
  * Hill & Fan (2008) Hill, K.M. & Fan, Y. 2008 Isolating segregation mechanisms in a split-bottom cell. Phys. Rev. Lett. 101, 088001\.
  * Hill & Tan (2014) Hill, K.M. & Tan, D.S. 2014 Segregation in dense sheared flows: gravity, temperature gradients, and stress partitioning. J. Fluid Mech. 756, 54–88.
  * Jones et al. (2018) Jones, R.P., Isner, A.B., Xiao, H., Ottino, J.M., Umbanhowar, P.B. & Lueptow, R.M. 2018 Asymmetric concentration dependence of segregation fluxes in granular flows. Phys. Rev. Fluids 3, 094304.
  * Kamrin & Koval (2012) Kamrin, K. & Koval, G. 2012 Nonlocal constitutive relation for steady granular flow. Phys. Rev. Lett. 108, 178301\.
  * Lees & Edwards (1972) Lees, A.W. & Edwards, S.F. 1972 The computer study of transport processes under extreme conditions. J. Phys. C: Solid State Phys. 5, 1921–1929.
  * Liu et al. (2023) Liu, D., Singh, H. & Henann, D.L. 2023 Coupled continuum modelling of size segregation driven by shear-strain-rate gradients and flow in dense, bidisperse granular media. J. Fluid Mech. 976, A16.
  * Liu et al. (2019) Liu, Y., Gonzalez, M. & Wassgren, C. 2019 Modeling granular material segregation using a combined finite element method and advection–diffusion–segregation equation model. Powder Technol. 346, 38–48.
  * Marks et al. (2012) Marks, B., Rognon, P. & Einav, I. 2012 Grainsize dynamics of polydisperse granular segregation down inclined planes. J. Fluid Mech. 690, 499–511.
  * Natarajan et al. (1995) Natarajan, V.V.R., Hunt, M.L. & Taylor, E.D. 1995 Local measurements of velocity fluctuations and diffusion coefficients for a granular material flow. J. Fluid Mech. 304, 1–25.
  * Plimpton (1995) Plimpton, S. 1995 Fast parallel algorithms for short-range molecular dynamics. J. Comp. Phys. 117, 1–19.
  * Rognon et al. (2007) Rognon, P.G., Roux, J.-N., Naaïm, M. & Chevoir, F. 2007 Dense flows of bidisperse assemblies of disks down an inclined plane. Phys. Fluids 19, 058101.
  * Schlick et al. (2015) Schlick, C.P., Fan, Y., Isner, A.B., Umbanhowar, P.B., Ottino, J.M. & Lueptow, R.M. 2015 Modeling segregation of bidisperse granular materials using physical control parameters in the quasi-2D bounded heap. AIChE J. 61, 1524–1534.
  * Silbert et al. (2001) Silbert, L.E., Ertaş, D., Grest, G.S., Halsey, T.C., Levine, D. & Plimpton, S.J. 2001 Granular flow down an inclined plane: Bagnold scaling and rheology. Phys. Rev. E 64, 051302.
  * Singh et al. (2023) Singh, H., Liu, D. & Henann, D.L. 2023 Continuum modeling of size-segregation and flow in dense, bidisperse granular media: Accounting for segregation driven by both pressure gradients and shear-strain-rate gradients, arXiv: 2307.08942.
  * Srivastava et al. (2021) Srivastava, I., Silbert, L.E., Grest, G.S. & Lechman, J.B. 2021 Viscometric flow of dense granular materials under controlled pressure and shear stress. J. Fluid Mech. 907, A18.
  * Trewhela et al. (2021) Trewhela, T., Ancey, C. & Gray, J.M.N.T. 2021 An experimental scaling law for particle-size segregation in dense granular flows. J. Fluid Mech. 916, A55.
  * Tripathi & Khakhar (2011) Tripathi, A. & Khakhar, D.V. 2011 Rheology of binary granular mixtures in the dense flow regime. Phys. Fluids 23, 113302\.
  * Tunuguntla et al. (2014) Tunuguntla, D.R., Bokhove, O. & Thornton, A.R. 2014 A mixture theory for size and density segregation in shallow granular free-surface flows. J. Fluid Mech. 749, 99–112.
  * Tunuguntla et al. (2016) Tunuguntla, D.R., Thornton, A.R. & Weinhart, T. 2016 From discrete elements to continuum fields: Extension to bidisperse systems. Comput. Part. Mech. 3, 349–365.
  * Umbanhowar et al. (2019) Umbanhowar, P.B., Lueptow, R.M. & Ottino, J.M. 2019 Modeling segregation in granular flows. Annu. Rev. Chem. Biomol. Eng. 10, 129–153.
  * Utter & Behringer (2004) Utter, B. & Behringer, R.P. 2004 Self-diffusion in dense granular shear flows. Phys. Rev. E 69, 031308.
  * Yang et al. (2021) Yang, L.Y.M., Zheng, Q.J., Bai, L. & Yu, A.B. 2021 Continuum modelling of granular segregation by coupling flow rheology and transport equation. Powder Technol. 378, 371–387.
  * Yang et al. (2022) Yang, L., Zheng, Q. & Yu, A. 2022 A continuum model for the segregation of bidisperse particles in a blade mixer. AIChE J. 68, e17734.