# Coupled continuum modeling of size-segregation driven by shear-strain-rate
gradients and flow in dense, bidisperse granular media

Daren Liu222These authors contributed equally to this work., Harkirat
Singh22footnotemark: 2, and David L. Henann888Email address for
correspondence<EMAIL_ADDRESS>
School of Engineering, Brown University, Providence, RI 02912, USA

###### Abstract

Dense granular systems that consist of particles of disparate sizes segregate
based on size during flow, resulting in complex, coupled segregation and flow
patterns. The ability to predict how granular mixtures segregate is important
in the design of industrial processes and the understanding of geophysical
phenomena. The two primary drivers of size-segregation are pressure gradients
and shear-strain-rate gradients. In this work, we isolate size-segregation
driven by shear-strain-rate gradients by studying two dense granular flow
geometries with constant pressure fields: gravity-driven flow down a long
vertical chute with rough parallel walls and annular shear flow with rough
inner and outer walls. We perform discrete element method (DEM) simulations of
dense flow of bidisperse granular systems in both flow geometries, while
varying system parameters, such as the flow rate, flow configuration size,
fraction of large/small grains, and grain-size ratio, and use DEM data to
inform continuum constitutive equations for the relative flux of large and
small particles. When the resulting continuum model for the dynamics of size-
segregation is coupled with the nonlocal granular fluidity model–a nonlocal
continuum model for dense granular flow rheology–we show that both flow fields
and segregation dynamics may be simultaneously captured using this coupled,
continuum system of equations.

##  1 Introduction

Dense granular systems in nature and industry are often non-monodisperse–i.e.,
consisting of particles of disparate sizes. In non-monodisperse granular
systems, the constituent grains segregate based on size during flow, forming
complex patterns (e.g., Shinbrot and Muzzio, 2000; Gray and Thornton, 2005;
Hill and Fan, 2008; Fan and Hill, 2010; Schlick et al., 2015; Gray, 2018;
Umbanhowar et al., 2019). The ability to predict the dynamics of segregation
is important across applications. For example, in geophysics, granular size
segregation can manifest in landslides and debris flows (Johnson et al.,
2012), in which larger grains segregate to the top of the flow, potentially
causing more damage, while in industry, size-segregation can be an undesirable
effect when blending granular constituents of various sizes.

The current understanding is that there are two driving forces for size-
segregation in dense granular flows. The first is pressure-gradients, which
are typically induced by gravity. In pressure-gradient-driven size-
segregation, small particles move more readily through the interstitial spaces
that open and close during flow through a process referred to in the
literature as “kinetic sieving,” leading to a system stratified along the
direction of pressure gradients (Savage and Lun, 1988; Gray and Thornton,
2005; Gray and Chugunov, 2006; Thornton et al., 2012; Fan et al., 2014). While
pressure-gradient-driven segregation has been the focus of significant study,
Hill and coworkers (Hill and Fan, 2008; Fan and Hill, 2010, 2011b; Hill and
Tan, 2014) demonstrated that grains can also segregate in inhomogeneous flows
along directions orthogonal to gravitationally-induced pressure gradients,
driven instead by gradients in the shear-strain-rate. As an example, this
mechanism has been observed in the split-bottom cell experiments of Hill and
Fan (2008). In these experiments, not only do the larger particles segregate
to the top of the cell, but they also segregate perpendicular to the direction
of pressure gradients towards more rapidly shearing regions. Shear-strain-
rate-gradient-driven segregation has received comparatively less attention in
modeling efforts.

Due to the complexity of flow and segregation patterns, developing a general,
predictive, continuum model for coupled size-segregation and flow in dense
granular materials remains an open challenge. Although much progress has been
made over recent decades (e.g., Savage and Lun, 1988; Gray and Thornton, 2005;
Gray and Chugunov, 2006; Fan and Hill, 2011b; Fan et al., 2014; Tunuguntla et
al., 2017; Gray, 2018; Umbanhowar et al., 2019), the development of continuum
models that are capable of simultaneously predicting the evolution of both
segregation and flow fields, based solely on the geometry of the flow
configuration, applied loads, and boundary/initial conditions is still in its
infancy. Instead, most continuum models for size-segregation require some flow
field quantity, such as the velocity or stress fluctuation field, to be
measured first from experiments or DEM simulations and then used as model
input. A crucial reason for the incompleteness of current models is the lack
of a dense granular flow rheology theory that may be coupled to segregation
models. A widely-used class of viscoplastic models for steady, dense granular
flow is based on the $\mu(I)$ rheology (MiDi, 2004; Jop et al., 2005; da Cruz
et al., 2005; Srivastava et al., 2021). One recent work that couples rheology
and segregation in dense granular flows is that of Barker et al. (2021), which
combines a regularized version of the $\mu(I)$ rheology (Barker and Gray,
2017) with a model for gravity-driven segregation. However, it has been well-
documented in the literature (e.g., Kamrin, 2019) that the $\mu(I)$ rheology,
even in its regularized form, can break down in the presence of spatial flow
inhomogeneity, which can be attributed to nonlocal effects not accounted for
in the $\mu(I)$ rheology.

To address this point, significant effort has gone into the development of
size-dependent, nonlocal continuum constitutive theories for dense granular
flow rheology, and coupling a nonlocal rheological model with a segregation
model provides a route to robust, simultaneous prediction of flow and
segregation fields. In this paper, we focus on the nonlocal granular fluidity
(NGF) model of Kamrin and coworkers (Kamrin and Koval, 2012; Henann and
Kamrin, 2013; Kamrin, 2019), which has been successfully applied to predicting
dense flows of monodisperse grains in a wide variety of flow geometries. Then,
the overarching aim of this paper is to formulate a predictive continuum
theory for simultaneous flow and size-segregation in dense granular systems by
integrating the NGF model with a phenomenological size-segregation model. This
is a broad goal, and in this paper, we narrow our focus to several simpler,
quasi-one-dimensional flow configurations. In most real-world flows, both the
pressure-gradient-driven and shear-strain-rate-gradient-driven segregation
mechanisms are present, making it difficult to disentangle them. Therefore,
our plan for this paper is to isolate and examine the shear-strain-rate-
gradient-driven mechanism. Specifically, we study the shear-strain-rate-
gradient-driven segregation mechanism by considering flows of dense,
bidisperse systems of both disks and spheres in flow geometries in which the
pressure field is spatially uniform–namely, vertical chute flow and annular
shear flow. Therefore shear-strain-rate-gradients are the sole drivers of
segregation. In order to inform continuum model development, we perform
discrete element method (DEM) simulations using the open source software
LAMMPS (Plimpton, 1995), which function as “numerical experiments.” The
coupled continuum model that we develop is then validated by comparing its
predictions of the transient evolution of segregation and flow fields against
additional DEM simulation results.

This paper is organized as follows. In Section 2, we discuss the continuum
model that we use to describe flow and size-segregation in bidisperse, dense
granular materials. Specifically, Sections 2.1 and 2.2 introduce the mixture
theory framework used to describe dense, bidisperse granular mixtures, and in
Section 2.3, we briefly revisit the $\mu(I)$ rheology and the NGF model for
monodisperse granular systems and discuss their extension to bidisperse
systems. Then in Section 2.4, we propose a model for shear-strain-rate-
gradient-driven size-segregation. In Sections 3 and 4, we consider granular
diffusion and shear-strain-rate-gradient-driven segregation, respectively, and
independently determine the two dimensionless material parameters that appear
in the size-segregation model for both disks and spheres. Then, in Section 5,
the proposed segregation model is coupled with the NGF model and applied to
both vertical chute flow and annular shear flow to predict the transient
evolution of the segregation dynamics, and the predicted continuum fields are
compared against DEM measurements. In the end, our model demonstrates a level
of fidelity in simultaneously predicting flow and segregation dynamics that
has not been previously achieved. We close with a discussion of the
segregation model and some concluding remarks in Section 6.

##  2 Continuum framework

In this section, we discuss the continuum framework used to describe dense,
bidisperse granular systems and propose constitutive equations for rheology
and size-segregation. Throughout, we utilize a mixture-theory-based approach,
which is common in continuum modeling of dense, bidisperse mixtures (e.g.,
Gray and Thornton, 2005; Gray and Chugunov, 2006; Fan and Hill, 2011b; Gray,
2018; Umbanhowar et al., 2019; Barker et al., 2021; Bancroft and Johnson,
2021; Duan et al., 2021), and we use standard component notation, which
supposes an underlying set of Cartesian basis vectors $\\{{\bf
e}_{i}|i=1,2,3\\}$, and in which the components of vectors, ${\bf v}$, and
tensors, $\bm{\sigma}$, are denoted by $v_{i}$ and $\sigma_{ij}$,
respectively. The Einstein summation convention is employed, and the Kronecker
delta, $\delta_{ij}$, is utilized to denote the components of the identity
tensor.

###  2.1 Bidisperse systems

Figure 1: A representative schematic of a dense, bidisperse granular system
consisting of two-dimensional disks.

We consider granular mixtures consisting of particles with two sizes–large
grains with an average diameter of $d^{\rm l}$ and small grains with an
average diameter of $d^{\rm s}$. We consider both two-dimensional systems of
disks, as illustrated in Fig. 1, and three-dimensional systems of spheres. To
eliminate the effect of density-based segregation (e.g., Tripathi and Khakhar,
2013) and isolate size-based segregation, all particles are made of the same
material with density $\rho_{\rm s}$, which represents the area-density for
disks and the volume-density for spheres. Throughout, we utilize the
notational convention in which we denote large-grain quantities using a
superscript ${\rm l}$ and small-grain quantities using a superscript ${\rm
s}$. The species-specific solid fractions–i.e., the areas occupied by each
species per unit total area for disks and the volumes occupied by each species
per unit total volume for spheres–are $\phi^{\rm l}$ and $\phi^{\rm s}$,
respectively, and the total solid fraction is $\phi=\phi^{\rm l}+\phi^{\rm
s}$. The concentration of each species then follows as $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 mixture grain size is defined as the sizes of both species
weighted by their concentrations, $\bar{d}=c^{\rm l}d^{\rm l}+c^{\rm s}d^{\rm
s}$. We make the common idealization that the total area for dense systems of
disks or total volume for dense systems of spheres does not change (Savage,
1998; Gray and Thornton, 2005; Gray and Chugunov, 2006; Fan and Hill, 2011b),
and therefore $\phi$ is idealized as constant at each point in space and at
each instant in time during the segregation process. We have verified in our
DEM simulations that area (or volume) dilatation at flow initiation occurs
over a much shorter time scale than the process of segregation, so that this
idealization is reasonable. Throughout this study, we use $\phi=0.8$ for
disks, and $\phi=0.6$ for spheres.

Regarding the kinematics of flow, each species has an associated partial
velocity, $v_{i}^{\rm l}$ and $v_{i}^{\rm s}$, and the mixture velocity is
given by $v_{i}=c^{\rm l}v_{i}^{\rm l}+c^{\rm s}v_{i}^{\rm s}$. The mixture
strain rate tensor is then defined using the mixture velocity in the standard
way: $D_{ij}=(1/2)(\partial v_{i}/\partial x_{j}+\partial v_{j}/\partial
x_{i})$, where $D_{kk}=0$ since we have assumed that the mixture area (or
volume) does not change. The equivalent shear strain-rate is defined as
$\dot{\gamma}=(2D_{ij}D_{ij})^{1/2}$.

Then, the relative area (or volume) flux for each grain type $\alpha={\rm
l},{\rm s}$ is defined through the difference between its partial velocity and
the mixture velocity as
$w_{i}^{\alpha}=c^{\alpha}\left(v_{i}^{\alpha}-v_{i}\right),$ so that
$w_{i}^{\rm l}+w_{i}^{\rm s}=0$. Conservation of mass for each species
requires that $D{c}^{\alpha}/Dt+\partial w_{i}^{\alpha}/\partial x_{i}=0$,
where $D(\bullet)/Dt$ is the material time derivative. Due to the fact that
$c^{\rm l}+c^{\rm s}=1$, only one of $c^{\rm l}$ and $c^{\rm s}$ is
independent. Therefore, we will utilize $c^{\rm l}$ as the field variable that
describes the dynamics of size-segregation in the following discussion, and
the evolution of $c^{\rm l}$ is governed by its conservation of mass equation:

$\dfrac{D{c}^{\rm l}}{Dt}+\dfrac{\partial w_{i}^{\rm l}}{\partial x_{i}}=0.$
(2.1)

###  2.2 Stress and the equations of motion

We recognize that the symmetric Cauchy stress tensor $\sigma_{ij}=\sigma_{ji}$
represents the Cauchy stress of the mixture, rather than the partial stress of
either species. Regarding stress-related quantities for the granular mixture,
we define the pressure $P=-(1/3)\sigma_{kk}$, the stress deviator
$\sigma_{ij}^{\prime}=\sigma_{ij}+P\delta_{ij}$, the equivalent shear stress
$\tau=(\sigma_{ij}^{\prime}\sigma_{ij}^{\prime}/2)^{1/2}$, and the stress
ratio $\mu=\tau/P$. The Cauchy stress is then governed by the standard
equations of motion

$\phi\rho_{\rm s}\frac{Dv_{i}}{Dt}=\frac{\partial\sigma_{ij}}{\partial
x_{j}}+b_{i},$ (2.2)

where $\phi$ is the constant total solid fraction, and $b_{i}$ is the non-
inertial body force per unit volume (typically gravitational). In order to
close the system of equations, we require (1) rheological constitutive
equations for the Cauchy stress $\sigma_{ij}$ and (2) a constitutive equation
for the flux $w_{i}^{\rm l}$, each of which are discussed in the following
subsections.

###  2.3 Rheological constitutive equations for bidisperse mixtures

In this section, we discuss the rheology of dense, bidisperse granular
mixtures. Our strategy for formulating rheological constitutive equations for
bidisperse mixtures is to relate mixture-related quantities, such as the
Cauchy stress $\sigma_{ij}$ and the strain-rate tensor $D_{ij}$, instead of
specifying constitutive equations for species-specific partial stresses and
then combining them to obtain the mixture stress.

The starting point of this discussion is the local inertial, or $\mu(I)$,
rheology (MiDi, 2004; Jop et al., 2005; da Cruz et al., 2005), which follows
from dimensional arguments. For a dense, monodisperse system of dry, stiff
grains with mean grain diameter $d$ subjected to homogeneous shearing, the
local inertial rheology asserts that the stress ratio $\mu$ is given through
the equivalent shear strain-rate $\dot{\gamma}$ and the pressure $P$ through
the dimensionless relationship $\mu=\mu_{\rm loc}(I)$, where
$I=\dot{\gamma}\sqrt{d^{2}\rho_{\rm s}/P}$ is the inertial number,
representing the ratio of the microscopic time scale associated with particle
motion $\sqrt{d^{2}\rho_{\rm s}/P}$ to the macroscopic time scale of applied
deformation $1/\dot{\gamma}$. As shown by Rognon et al. (2007) and Tripathi
and Khakhar (2011), the inertial rheology function $\mu_{\rm loc}(I)$ may be
straightforwardly generalized from monodisperse to bidisperse systems by
defining the inertial number for a bidisperse system as
$I=\dot{\gamma}\sqrt{\bar{d}^{2}\rho_{\rm s}/P}$, where the average mixture
grain size for a bidisperse system $\bar{d}$ has been used in place of $d$ for
a monodisperse system. Then, the same local rheology function $\mu_{\rm
loc}(I)$ utilized for the monodisperse system may be used for bidisperse
systems without any changes to the parameters appearing in the fitting
function. This approach neglects potential effects of new dimensionless
quantities that arise in a bidisperse granular system, such as the grain size
ratio $d^{\rm l}/d^{\rm s}$, but has been shown to capture DEM data well
(Rognon et al., 2007; Tripathi and Khakhar, 2011).

Figure 2: (a) Configuration for two-dimensional DEM simulations of bidisperse
simple shear flow. Upper and lower layers of black grains denote rough walls.
Dark gray grains indicate large flowing grains, and light gray grains indicate
small flowing grains. A $10\%$ polydispersity is utilized for each species to
prevent crystallization. (b) The local inertial rheology ($\mu$ versus
$I=\dot{\gamma}\sqrt{\bar{d}^{2}\rho_{\rm s}/P}$) for monodisperse as well as
bidisperse mixtures of disks for grain-size ratios of $d^{\rm l}/d^{\rm
s}=1.5$, 2.0, 2.5, and 3.0 and $c^{\rm l}=0.5$. The solid black line
represents the best fit to the monodisperse DEM data using (2.3) with
$\mu_{\rm s}=0.272$ and $b=1.168$. (c) The local inertial rheology for
monodisperse and bidisperse mixtures of spheres for grain size ratios of
$d^{\rm l}/d^{\rm s}=1.5$ and 2.0 and $c^{\rm l}=0.5$ along with the DEM data
of Tripathi and Khakhar (2011). The solid black curve represents the best fit
to the monodisperse DEM data using (2.4) with $\mu_{\rm s}=0.37$,
$\mu_{2}=0.95$, and $I_{0}=0.58$.

To demonstrate this point, consider DEM simulations of homogeneous, simple
shearing of a dense, bidisperse system of disks, illustrated in Fig. 2(a) for
the case of $d^{\rm l}/d^{\rm s}=1.5$ and a system-wide large-grain
concentration of $c^{\rm l}=0.5$. Details of the simulated granular systems,
including grain interaction properties, for both two-dimensional disks and
three-dimensional spheres are given in Appendix A.1. The large particles are
dark gray, and the small particles are light gray. With the system-wide mean
grain size denoted by $\bar{d}_{0}=c^{\rm l}d^{\rm l}+(1-c^{\rm l})d^{\rm s}$,
the rectangular domain has a length of $L=60\bar{d}_{0}$ in the $x$-direction
and a height of $H=60\bar{d}_{0}$ in the $z$-direction, which is filled with
$\sim 5000$ flowing grains. Shearing is driven through the relative motion of
two parallel, rough walls, which each consist of a thin layer of touching
glued grains, which are denoted as black in Fig. 2(a). The bottom wall is
fixed, and the top wall moves with a velocity $v_{\rm w}$ along the
$x$-direction. Following previous works in the literature (da Cruz et al.,
2005; Koval et al., 2009; Kamrin and Koval, 2012; Zhang and Kamrin, 2017; Liu
and Henann, 2018; Kim and Kamrin, 2020), the $z$-position of the top wall is
not fixed but continuously adjusted using a feedback scheme so that the normal
stress applied by the top wall is maintained at a target value of
$\sigma_{zz}(z=0)=-P_{\rm w}$. Periodic boundary conditions are applied along
the $x$-direction. For homogeneous simple shearing, no segregation will occur
since the flow is homogeneous and no pressure or strain-rate gradients are
present. We utilize the DEM procedures described in detail in Liu and Henann
(2018) in order to extract the relationship between $\mu$ and $I$ for
bidisperse mixtures with grain size ratios of $d^{\rm l}/d^{\rm
s}=1.5,2.0,2.5$, and 3.0 and $c^{\rm l}=0.5$. The simulated relationships are
plotted in Fig. 2(b) using triangular symbols of different colors, along with
the monodisperse data from Liu and Henann (2018) plotted as gray circles. The
relationship between $\mu$ and $I$ for dense systems of disks is observed to
be approximately independent of $d^{\rm l}/d^{\rm s}$. As for the monodisperse
case, the DEM data for bidisperse mixtures of disks may be fit by a linear,
Bingham-like functional form:

$\mu_{\rm loc}(I)=\mu_{\rm s}+bI,$ (2.3)

as shown by the solid line in Fig. 2(b), where $\mu_{\rm s}=0.272$ and
$b=1.168$ are the dimensionless material parameters for monodisperse disks
(Liu and Henann, 2018).

Similarly, we consider DEM simulations of homogeneous, simple shearing of
dense, bidisperse systems of spheres. The simulation domain consists of a
rectangular box of length $L=20\bar{d}_{0}$ in the $x$-direction (i.e., the
shearing direction), width $W=10\bar{d}_{0}$ in the $y$-direction (i.e., the
direction perpendicular to the plane of shearing), and height
$H=40\bar{d}_{0}$ in the $z$-direction. The domain is filled with $\sim$10000
flowing grains, and periodic boundary conditions are applied along both the
$x$\- and $y$-directions. The simulation domain is bounded along the
$z$-direction by two parallel, rough walls, consisted of touching glued
grains, and as for the case of disks, shearing along the $x$-direction and
normal stress along the $z$-direction are applied by the walls. We perform DEM
simulations of steady simple shearing for size ratios of $d^{\rm l}/d^{\rm
s}=1.5$ and $2.0$ for a system-wide large-grain concentration of $c^{\rm
l}=0.5$ as well as for the monodisperse case over a range of top wall
velocities. The $\mu$ versus $I$ relationship extracted from DEM simulations
for these cases along with data from the prior DEM study of Tripathi and
Khakhar (2011) collapse quite well as shown in Fig. 2(c), showing minimal
dependence on $d^{\rm l}/d^{\rm s}$. This relationship for dense systems of
spheres may be fitted using a nonlinear functional form of Jop et al. (2005)
for $\mu_{\rm loc}(I)$:

$\mu_{\rm loc}(I)=\mu_{\rm s}+\frac{\mu_{2}-\mu_{\rm s}}{I_{0}/I+1},$ (2.4)

as shown by the solid curve in Fig. 2(c), where $\\{\mu_{\rm
s}=0.37,\,\mu_{2}=0.95,\,I_{0}=0.58\\}$ are the dimensionless parameters for
frictional spheres. (We note that these parameters are nearly the same as
those determined by Zhang and Kamrin (2017) for monodisperse frictional
spheres.) In this way, one may capture the rheology of bidisperse mixtures of
both disks and spheres in homogeneous simple shearing without introducing
additional fitting functions or adjustable parameters beyond those used for
the monodisperse case.

Despite the successes of the local inertial rheology in capturing steady,
homogeneous shear flow, it has been well established in the literature that a
local rheological modeling approach cannot be applied to a broad set of
inhomogeneous flows, such as annular shear flow (Tang et al., 2018), split-
bottom flow (Fenistein and van Hecke, 2003), and gravity-driven heap flow
(Komatsu et al., 2001). Therefore, to consider inhomogeneous flows of
bidisperse granular systems, it is necessary to generalize a nonlocal
rheological modeling approach to the case of dense, bidisperse mixtures. In
the present work, we focus attention on the nonlocal granular fluidity (NGF)
model of Kamrin and Koval (2012), which has been shown to robustly capture a
variety of inhomogeneous, steady flows of monodisperse granular systems
(Kamrin, 2019). As is standard in the NGF model, we introduce the granular
fluidity $g$, which is a positive, scalar field quantity, and recognize that
$g$ represents the fluidity of the mixture. (See Zhang and Kamrin (2017) and
Kim and Kamrin (2020) for further discussion of the kinematic description of
the granular fluidity field for monodisperse granular systems.) Then, we
utilize the steady-state form of the NGF model, which relates the stress
state, the strain-rate, and the granular fluidity through two constitutive
equations: (1) the flow rule and (2) the nonlocal rheology.

First, invoking the common idealization that the Cauchy stress deviator and
the strain-rate tensor are co-directional (Rycroft et al., 2009), the flow
rule relates the Cauchy stress tensor $\sigma_{ij}$, the strain-rate tensor
$D_{ij}$, and the granular fluidity through

$\sigma_{ij}=-P\delta_{ij}+2\dfrac{P}{g}D_{ij}.$ (2.5)

Taking the magnitude of the deviatoric part of (2.5) and rearranging leads to
the following scalar form of the flow rule:

$\dot{\gamma}=g\mu.$ (2.6)

Second, the granular fluidity of the bidisperse mixture is governed by the
following differential relation:

$g=g_{\rm loc}(\mu,P)+\xi^{2}(\mu)\frac{\partial^{2}g}{\partial x_{i}\partial
x_{i}},$ (2.7)

where $g_{\rm loc}(\mu,P)$ is the local fluidity function and $\xi(\mu)$ is
the stress-dependent cooperativity length. The local fluidity function gives
the granular fluidity during steady, homogeneous shear flow at a given state
of stress and is related to the local inertial rheology function $\mu_{\rm
loc}(I)$. Denote the inverted form of the local inertial rheology function
$\mu_{\rm loc}(I)$ as

$I_{\rm loc}(\mu)=\left\\{\begin{array}[]{cl}\mu^{-1}_{\rm loc}(\mu)&\text{if
$\mu>\mu_{\rm s}$,}\\\\[4.0pt] 0&\text{if $\mu\leq\mu_{\rm
s}$,}\end{array}\right.$ (2.8)

which is a function of the stress ratio $\mu$. Then, consistent with the new
definition of the inertial number involving $\bar{d}$ for a bidisperse
mixture, the local fluidity function is $g_{\rm
loc}(\mu,P)=\sqrt{P/\bar{d}^{2}\rho_{\rm s}}\,I_{\rm loc}(\mu)/\mu$. For the
case of bidisperse disks, using (2.3), the local fluidity function is

$g_{\rm loc}(\mu,P)=\left\\{\begin{array}[]{cl}\sqrt{\dfrac{P}{\rho_{\rm
s}\bar{d}^{2}}}\,\dfrac{(\mu-\mu_{\rm s})}{b\mu}&\text{if $\mu>\mu_{\rm
s}$,}\\\\[8.0pt] 0&\text{if $\mu\leq\mu_{\rm s}$,}\end{array}\right.$ (2.9)

with $\\{\mu_{\rm s}=0.272,b=1.168\\}$, and for the case of bidisperse
spheres, using (2.4), the local fluidity function is

$\displaystyle g_{\rm
loc}(\mu,P)=\left\\{\begin{array}[]{cl}I_{0}\sqrt{\dfrac{P}{\rho_{\rm
s}\bar{d}^{2}}}\,\dfrac{(\mu-\mu_{\rm s})}{\mu(\mu_{2}-\mu)}&\text{if
$\mu>\mu_{\rm s}$,}\\\\[8.0pt] 0&\text{if $\mu\leq\mu_{\rm s}$,}\\\
\end{array}\right.$ (2.12)

with $\\{\mu_{\rm s}=0.37,\,\mu_{2}=0.95,\,I_{0}=0.58\\}$. No additional
adjustable parameters beyond those used to describe the local inertial
rheology for the monodisperse case are introduced in the local fluidity
function.

As discussed in several of our previous works (Henann and Kamrin, 2014; Kamrin
and Henann, 2015; Liu and Henann, 2017), the functional form for the
cooperativity length $\xi(\mu)$ is also connected to the choice of the
$\mu_{\rm loc}(I)$ function. Without going into details here, the functional
forms for the cooperativity length corresponding to (2.3) and (2.4) are

$\xi(\mu)=\frac{A\bar{d}}{\sqrt{|\mu-\mu_{\rm
s}|}}\quad\text{and}\quad\xi(\mu)=A\bar{d}\sqrt{\dfrac{(\mu_{2}-\mu)}{(\mu_{2}-\mu_{\rm
s})|\mu-\mu_{\rm s}|}},$ (2.13)

respectively. In the monodisperse case, the cooperativity length is directly
proportional to the grain size $d$, and in (2.13) for the bidisperse case, it
is taken to be proportional to $\bar{d}$. This is analogous to the process
undertaken above for generalizing the local inertial rheology, in which $d$
for monodisperse grains is replaced by $\bar{d}$ for bidisperse grains.
Moreover, in (2.13), the parameter $A$ is a dimensionless material constant,
referred to as the non-local amplitude, which quantifies the spatial extent of
cooperative effects. Again, following an approach analogous to the
generalization of the local inertial rheology, we utilize values of $A$
previously determined for monodisperse frictional disks and spheres–namely,
$A=0.9$ as determined by Liu and Henann (2018) for monodisperse disks and
$A=0.43$ as determined by Zhang and Kamrin (2017) for monodisperse spheres.
These choices of $A$ for bisdisperse mixtures will be tested in later sections
by comparing flow fields predicted by the NGF model to measured flow fields in
DEM simulations of bidisperse, inhomogeneous flows.

###  2.4 Segregation model

The segregation model consists of the constitutive equation for the large-
grain flux $w_{i}^{\rm l}$. In the present work, we focus on dense flows in
the absence of pressure gradients, and we take the large-grain flux
$w_{i}^{\rm l}$ to be comprised of two contributions: (1) a diffusion flux
$w_{i}^{\rm diff}$ and (2) a shear-strain-rate-gradient-driven segregation
flux $w_{i}^{\rm seg}$, so that

$w_{i}^{\rm l}=w_{i}^{\rm diff}+w_{i}^{\rm seg}.$ (2.14)

First, the diffusion flux acts counter to segregation to mix the species and
is taken to be given in the standard form, in which the diffusion flux is
driven by concentration gradients: $w_{i}^{\rm diff}=-D\left(\partial c^{\rm
l}/\partial x_{i}\right)$, where $D$ is the binary diffusion coefficient
(Utter and Behringer, 2004; Artoni et al., 2021; Bancroft and Johnson, 2021).
Based on dimensional arguments, we expect that

$D=C_{\rm diff}\bar{d}^{2}\dot{\gamma},$ (2.15)

where $C_{\rm diff}$ is a dimensionless material parameter which remains to be
calibrated (Fan et al., 2014; Tripathi and Khakhar, 2013). Therefore, we have
that the diffusion flux is

$w_{i}^{\rm diff}=-C_{\rm diff}\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm
l}}{\partial x_{i}}.$ (2.16)

Second, regarding segregation, a major question is what field quantity drives
the segregation flux in the absence of pressure gradients. Gradients of a
number of kinematic quantities are possible–e.g., strain-rate, velocity
fluctuations, or fluidity. For perspective, we note that recent works (Fan and
Hill, 2011b; Hill and Tan, 2014; Tunuguntla et al., 2016, 2017) have shown
that gradients in kinetic stress, which are defined through the velocity
fluctuations, correlate well with segregation flux. In the present work, we
adopt the simplest approach and hypothesize that the segregation flux is
driven by gradients in the shear strain-rate $\dot{\gamma}$ and take the
segregation flux to be given in the following phenomenological form:

$w_{i}^{\rm seg}=C_{\rm seg}^{\rm S}\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})\dfrac{\partial\dot{\gamma}}{\partial x_{i}}.$ (2.17)

The factor $c^{\rm l}(1-c^{\rm l})$ ensures that segregation ceases when the
bidisperse mixture becomes either all large $(c^{\rm l}=1)$ or all small
$(c^{\rm l}=0)$ grains, and the factor $\bar{d}^{2}$ is present for
dimensional consistency. The quantity $C^{\rm S}_{\rm seg}$ is a dimensionless
material property. While it is possible for $C^{\rm S}_{\rm seg}$ to depend on
the size ratio $d^{\rm l}/d^{\rm s}$, we will demonstrate that this effect is
negligible in the DEM simulations of Section 4 and therefore treat $C^{\rm
S}_{\rm seg}$ as a constant, dimensionless material parameter, which will be
determined by fitting to DEM simulation results for disks and spheres,
respectively.

Combining (2.16), (2.17), and (2.14) with conservation of mass (2.1), we
obtain the following differential relation governing the dynamics of $c^{\rm
l}$:

$\dfrac{D{c}^{\rm l}}{Dt}+\dfrac{\partial}{\partial x_{i}}\left(-C_{\rm
diff}\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm l}}{\partial x_{i}}+C^{\rm
S}_{\rm seg}\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})\dfrac{\partial\dot{\gamma}}{\partial x_{i}}\right)=0,$ (2.18)

where $\\{C_{\rm diff},C^{\rm S}_{\rm seg}\\}$ represent two constant
dimensionless material parameters that remain to be determined.

We close this section by noting that the incompressibility constraint, the
equations of motion (2.2), the nonlocal rheology (2.7), and the segregation
dynamics equation (2.18) represent a closed system of equations for the
velocity field $v_{i}$, the pressure field $P$, the fluidity field $g$, and
the large-grain concentration field $c^{\rm l}$, which may be used to
simultaneously predict flow fields and segregation dynamics in the absence of
pressure gradients.

##  3 Diffusion flux

In this section, we determine values of $C_{\rm diff}$ for dense, bidisperse
systems of frictional disks and spheres. Consider homogeneous simple shear
flow of such a bidisperse mixture, as shown in Fig. 2(a) for disks. Again, no
segregation occurs in this setting, since neither of the segregation driving
forces (pressure gradients or shear-strain-rate gradients) are present
(Tripathi and Khakhar, 2011). During steady, simple shearing, the motion of
individual grains in the direction transverse to flow (the $z$-direction in
Fig. 2(a)) approximates a random walk for both two-dimensional systems of
disks and three-dimensional systems of spheres. Therefore, by measuring the
mean square displacement (MSD) of the system of $N$ particles as a function of
time, we may determine the binary diffusion coefficient $D$ (Natarajan et al.,
1995; Dufty et al., 2002; Utter and Behringer, 2004; Fan et al., 2015; Kharel
and Rognon, 2017; Bancroft and Johnson, 2021) through

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

where $z_{n}(t)$ is the $z$-coordinate of the $n$th grain at time $t$. We
simulate homogeneous, steady simple shear flows of disks for grain-size ratios
of $d^{\rm l}/d^{\rm s}=1.5$, 2.0, 2.5, 3.0, and 4.0 and at various shearing
rates. We also simulate homogeneous, steady simple shear flows of spheres for
the monodisperse case as well as for grain-size ratios of $d^{\rm l}/d^{\rm
s}=1.5$, 2.0, and 2.5 over a range of shearing rates. To avoid wall effects in
the calculation of the MSD (3.1), grains that are initially within
$15\bar{d}_{0}$ of either the top or bottom wall in Fig. 2(a) are excluded
from the system of particles used to calculate the MSD for disks, leaving a
set of $N\approx 2400$ grains. For spheres, particles initially within
$5\bar{d}_{0}$ of the top and bottom walls are excluded, so that a set of
$N\approx 9000$ grains are used to calculate the MSD. Both large and small
grains are included in the calculation of the MSD for the mixture. After a
sufficiently long time, the calculated MSD is linear in time in all cases for
both disks and spheres, allowing one to extract the diffusion coefficient $D$
for each case.

The diffusion coefficient $D$ is plotted against $\dot{\gamma}\bar{d}^{2}$
(with both quantities normalized by $d^{\rm s}\sqrt{P_{\rm w}/\rho_{\rm s}}$)
for disks in Fig. 3(a) and for spheres in Fig. 3(b). 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 size ratios and shearing
rates considered. A best fit of the slopes of the linear relations–the solid
black lines in Figs. 3(a) and (b)–yields:

$C_{\rm diff}=\dfrac{D}{\dot{\gamma}\bar{d}^{2}}=0.20\text{ for disks and
}C_{\rm diff}=\dfrac{D}{\dot{\gamma}\bar{d}^{2}}=0.045\text{ for spheres.}$
(3.2)

We note that our results are consistent with previous results in the
literature. For example, the recent work of Bancroft and Johnson (2021) found
a value of $C_{\rm diff}\approx 0.05$ with a weak dependence on the inertial
number for dense spheres. In order to further assess the fitted value of
$C_{\rm diff}$ in a diffusion-dominated setting, we have performed a
consistency test for disks by considering simple shearing of an initially
fully-segregated cell, which is described in Appendix B. In this case,
diffusion drives remixing of the two species. Using the fitted value of
$C_{\rm diff}$ for disks in (3.2), we are able to quantitatively capture the
diffusive remixing process, which provides confidence in our fitted value of
$C_{\rm diff}$.

Figure 3: The binary diffusion coefficient $D$, calculated using the mean
square displacement (3.1), versus $\dot{\gamma}\bar{d}^{2}$ in homogeneous,
steady simple shear DEM simulations. (a) Simple shearing of bidisperse
mixtures of disks for grain-size ratios of $d^{\rm l}/d^{\rm
s}=1.5,2.0,2.5,3.0$, and 4.0. Both axes are normalized by $d^{\rm
s}\sqrt{P_{\rm w}/\rho_{\rm s}}$. Each symbol represents $D$ calculated from
one DEM simulation of a specified size ratio at one shearing rate. The solid
line represents the best fit of a linear relation with $C_{\rm diff}=0.20$.
(b) Simple shearing of bidisperse mixtures of spheres for the monodisperse
case and for grain-size ratios of $d^{\rm l}/d^{\rm s}=1.5,2.0$, and 2.5. The
solid line represents the best fit of a linear relation with $C_{\rm
diff}=0.045$.

##  4 Shear-strain-rate-gradient-driven segregation flux

Having independently determined the material parameter $C_{\rm diff}$ for both
frictional disks and spheres, we next turn to testing the constitutive
equation for the shear-strain-rate-gradient-driven segregation flux (2.17) and
determining the material parameter $C^{\rm S}_{\rm seg}$ by studying two
representative flow configurations in the absence of pressure gradients: (1)
vertical chute flow and (2) annular shear flow.

###  4.1 Vertical chute flow

Figure 4: (a) Initial well-mixed configuration for two-dimensional DEM
simulation of bidisperse vertical chute flow with $4327$ flowing grains. The
chute width is $W=60\bar{d}_{0}$. As in Fig. 2, black grains on both sides
represent rough walls (only large particles are used as wall grains here). (b)
Segregated configuration after flowing for a total simulation time of
$\tilde{t}=t/\left(d^{\rm s}\sqrt{\rho_{\rm s}/P_{\rm w}}\right)=4.3\times
10^{5}$. (c) Spatiotemporal evolution of the large grain concentration field.
Spatial profiles of (d) the concentration field $c^{\rm l}$ and (e) the
normalized velocity field $(v_{\rm cen}-v_{z})\sqrt{\rho_{\rm s}/P_{\rm w}}$
at three times ($\tilde{t}=4\times 10^{3}$, $4\times 10^{4}$ and $4\times
10^{5}$) as indicated by the horizontal lines in (c).

Consider a dense, bidisperse granular mixture flowing down a long vertical
chute with parallel, rough walls separated by a distance $W$ under the action
of gravity $G$. This flow geometry has been utilized extensively in the
literature to study dense flows of monodisperse, frictional disks (Kamrin and
Koval, 2012; Liu and Henann, 2018) and spheres (Zhang and Kamrin, 2017; Kim
and Kamrin, 2020) as well as flows of bidisperse, frictional spheres (Fan and
Hill, 2011a, b). Beginning with the case of bidisperse disks, the DEM setup is
shown in Fig. 4(a) for $W=60\bar{d}_{0}$, where $\bar{d}_{0}$ is the system-
wide average grain size. In all cases for disks, we take the chute length to
be $L=60\bar{d}_{0}$ and apply periodic boundary conditions along the
$z$-direction. The parallel, rough walls consist of touching glued large
grains, denoted as black in Fig. 4(a). The left vertical wall is fixed, and
the right wall is fixed in the $z$-direction but can move slightly in the
$x$-direction so as to maintain a constant compressive normal stress $P_{\rm
w}$ on the granular material, utilizing the same wall-position control method
described in Liu and Henann (2018). We have verified that the chute length $L$
is sufficiently large, so that it does not affect the resulting flow and
segregation fields and all fields are invariant along the $z$-direction. In
the resulting flow fields, the only non-zero component of the velocity is
$v_{z}$, which only depends on the cross-channel coordinate $x$. A typical
steady velocity field is qualitatively sketched in Fig. 4(a), illustrating
that the shear-strain-rate is greatest at the walls $(x=\pm W/2)$.

In all of our DEM simulations of vertical chute flow of bidisperse disks, we
observe that the stress field quickly becomes independent of time, so that
macroscopic inertia (i.e., the left-hand-side of (2.2)) may be neglected.
Moreover, we observe that the normal stresses are approximately equal, i.e.,
$\sigma_{zz}\approx\sigma_{xx}$. Therefore, due to the force balance along the
$z$-direction, the equivalent shear stress field is
$\tau(x)=|\sigma_{xz}(x)|=|\sigma_{zx}(x)|=\phi\rho_{\rm s}G|x|$, where $x$ is
measured from the centerline of the chute, and due to the force balance along
the $x$-direction, the pressure field is $P(x)=-\sigma_{xx}(x)=P_{\rm w}$. The
stress ratio field then follows as

$\mu(x)=\mu_{\rm w}\left(\dfrac{|x|}{W/2}\right),$ (4.1)

where $\mu_{\rm w}=\phi\rho_{\rm s}GW/2P_{\rm w}$ is the maximum value of
$\mu$, occurring at the walls $(x=\pm W/2)$. We note that while flow is driven
by gravity, the pressure field is constant throughout the chute, and no
pressure gradients are present. Therefore, segregation occurs only due to
shear-strain-rate-gradients, enabling us to consider this effect in isolation.

Apart from the grain interaction properties that are held constant throughout
this work (Appendix A.1), there are four important dimensionless parameters
that fully describe each case of vertical chute flow of dense, bidisperse
granular mixtures: (1) $W/\bar{d}_{0}$, the dimensionless chute width; (2)
$\mu_{\rm w}$, the maximum stress ratio, which occurs at the walls and
controls the total flow rate; (3) $c^{\rm l}_{0}(x)$, the initial large-grain
concentration, which is not necessarily constant but can be a spatially-
varying field; and (4) $d^{\rm l}/d^{\rm s}$, the bidisperse grain-size ratio.
This list of system parameters $\\{W/\bar{d}_{0},$ $\mu_{\rm w},$ $c^{\rm
l}_{0},$ $d^{\rm l}/d^{\rm s}\\}$ specifies the geometry, loads, and initial
conditions of a given case of vertical chute flow. As a representative base
case, we consider the parameter group $\\{W/\bar{d}_{0}=60,\mu_{\rm
w}=0.45,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm s}=1.5\\}$. We then run the
corresponding DEM simulation starting from the well-mixed initial
configuration shown in Fig. 4(a) and observe that after a simulation time of
$\tilde{t}=t/(d^{\rm s}\sqrt{\rho_{\rm s}/P_{\rm w}})=4.3\times 10^{5}$, the
large, dark-gray grains segregate towards the regions near the walls where the
shear-strain-rate is greatest, while the small, light-gray grains gather in
bands just inside these regions, as shown in Fig. 4(b). A well-mixed core
persists along the center of the vertical chute where the shear-strain-rate is
nearly zero. To obtain a more quantitative picture of the segregation process,
we coarse-grain the concentration field $c^{\rm l}$ in both space and time and
plot contours of the spatiotemporal evolution of $c^{\rm l}$ in Fig. 4(c). The
large-grain concentration field evolves quickly in time during flow
initiation. Then, over longer times, the evolution becomes slower. Spatial
profiles of the concentration and velocity fields at three snapshots in
time–specifically, $\tilde{t}=t/(d^{\rm s}\sqrt{\rho_{\rm s}/P_{\rm
w}})=4\times 10^{3}$, $4\times 10^{4}$ and $4\times 10^{5}$ as indicated by
the horizontal lines in Fig. 4(c)–are plotted in Figs. 4(d) and (e). These
three snapshots correspond to early, medium, and late times with respect to
the segregation process. The spatial $c^{\rm l}$ profiles shown in Fig. 4(d)
demonstrate the transition from a well-mixed state to a segregated state with
large-grain-rich and small-grain-rich regions. In Fig. 4(e), the normalized
velocity fields $(v_{\rm cen}-v_{z})\sqrt{\rho_{\rm s}/P_{\rm w}}$, relative
to the velocity at the center of the chute, $v_{\rm cen}=v_{z}(x=0)$, show
that the velocity field rapidly develops into a steady flow field, even while
the segregation process is still ongoing, and the $c^{\rm l}$ field continues
to evolve.

Figure 5: Collapse of $C_{\rm diff}\bar{d}^{2}\dot{\gamma}({\partial c^{\rm
l}}/{\partial x})$ versus $\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})({\partial\dot{\gamma}}/{\partial x})$ for several cases of vertical chute
flow of (a) bidisperse disks and (b) bidisperse spheres. Symbols represent
coarse-grained, quasi-steady DEM field data, and the solid lines are the best
linear fits using (a) $C^{\rm S}_{\rm seg}=0.23$ for disks and (b) $C^{\rm
S}_{\rm seg}=0.08$ for spheres.

At long times, near the end of the simulated time window ($\tilde{t}=t/(d^{\rm
s}\sqrt{\rho_{\rm s}/P_{\rm w}})\gtrsim 3\times 10^{5}$), the concentration
field evolves very slowly, so that $D{c}^{\rm l}/Dt\approx 0$, and the state
of segregation may be regarded as quasi-steady. Therefore, according to
equations (2.1) and (2.14) and the no-flux boundary condition at the walls,
the total flux is approximately zero ($w_{i}^{\rm l}=w_{i}^{\rm
diff}+w_{i}^{\rm seg}\approx 0_{i}$) at each $x$-position, meaning that the
segregation flux is approximately balanced by the diffusion flux in this
quasi-steady flow regime. Then, using the expressions for the two fluxes,
equations (2.16) and (2.17), this observation implies that

$C_{\rm diff}\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm l}}{\partial
x}\approx C^{\rm S}_{\rm seg}\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})\dfrac{\partial\dot{\gamma}}{\partial x}.$ (4.2)

The field quantities appearing in this expression may be obtained by coarse-
graining the DEM data in the quasi-steady flow regime. Therefore, since
$C_{\rm diff}$ has been previously determined, (4.2) may be used to determine
the parameter $C^{\rm S}_{\rm seg}$ as follows. First, we acquire the field
quantities $c^{\rm l}$ (and hence $\bar{d}$) and $v_{z}$ by spatially coarse-
graining $152$ evenly-distributed snapshots in time in the quasi-steady regime
($\tilde{t}>3\times 10^{5}$).111Results are not particularly sensitive to the
quasi-steady regime criterion, and we only provide this value as a guideline.
Then, we arithmetically average these fields in time, yielding fields that
only depend on the spatial coordinate $x$, and take spatial gradients to
obtain ${\partial c^{\rm l}}/{\partial x}$, $\dot{\gamma}={\partial
v_{z}}/{\partial x}$, and ${\partial\dot{\gamma}}/{\partial
x}={\partial^{2}v_{z}}/{\partial x^{2}}$.222Detailed descriptions of our
coarse-graining techniques may be found in Appendix A.2. Next, as suggested by
(4.2), we plot $C_{\rm diff}\bar{d}^{2}\dot{\gamma}({\partial c^{\rm
l}}/{\partial x})$ versus $\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})({\partial\dot{\gamma}}/{\partial x})$ in Fig. 5(a), in which each point
represents a unique $x$-position. A linear relation is observed, supporting
our choice for the form of the constitutive equation for the segregation flux
(2.17). In order to obtain further evidence for this choice, we consider four
additional cases: (1) a lower flow rate case $\\{W/\bar{d}_{0}=60,\mu_{\rm
w}=0.375,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm s}=1.5\\}$; (2) a narrower channel
case $\\{W/\bar{d}_{0}=40,\mu_{\rm w}=0.45,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=1.5\\}$; (3) a more large grains case $\\{W/\bar{d}_{0}=60,\mu_{\rm
w}=0.45,c^{\rm l}_{0}=0.75,d^{\rm l}/d^{\rm s}=1.5\\}$; and (4) a larger size
ratio case $\\{W/\bar{d}_{0}=60,\mu_{\rm w}=0.45,c^{\rm l}_{0}=0.5,d^{\rm
l}/d^{\rm s}=3.0\\}$. Coarse-graining the quasi-steady fields for each case
and including the field data in Fig. 5(a), we observe a strong linear
collapse. Finally, the dimensionless material parameter $C_{\rm seg}^{\rm S}$
may be obtained from the slope of the linear relation in Fig. 5(a) (indicated
by the solid line). We determine the numerical value for disks to be $C^{\rm
S}_{\rm seg}=0.23$ and note that this value indeed appears to be independent
of the grain-size ratio $d^{\rm l}/d^{\rm s}$.

We carry out an analogous set of DEM simulations for dense, bidisperse
mixtures of spheres to determine the value of $C^{\rm S}_{\rm seg}$ for
spheres. The DEM setup for spheres is similar to that shown in Fig 4(a) for
disks with a domain size of length $L=20\bar{d}_{0}$ in the $z$-direction,
width $W$ in the $x$-direction, which is varied in our DEM simulations, and
out-of-plane thickness $H=10\bar{d}_{0}$ in the $y$-direction. Periodic
boundary conditions are applied along both the $y$\- and $z$-directions for
spheres, and a constant compressive normal stress $\sigma_{xx}=-P_{w}$ is
applied using the same feedback scheme as utilized for disks. In DEM
simulations of dense flows of spheres, we observe normal stress differences,
in which the normal stresses $\sigma_{yy}$ and $\sigma_{zz}$ are slightly
different from the prescribed value of $\sigma_{xx}=-P_{w}$, which is a
widely-reported feature of dense flows of spheres in the literature (e.g.,
Srivastava et al., 2021). However, all normal stresses, and hence the pressure
field, are spatially uniform, and segregation occurs only due to shear-strain-
rate gradients. As for disks, the set of system parameters $\\{W/\bar{d}_{0},$
$\mu_{\rm w},$ $c^{\rm l}_{0},$ $d^{\rm l}/d^{\rm s}\\}$ specifies the
geometry, loads, and initial conditions for a given case of vertical chute
flow. Here, the dimensionless parameter $\mu_{\rm w}=\phi\rho_{\rm
s}GW/2P_{\rm w}$ continues to control the total flow rate down the chute. We
consider five different cases: (1) a base case $\\{W/\bar{d}_{0}=60,\mu_{\rm
w}=0.51,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm s}=1.5\\}$; (2) a lower flow rate
case $\\{W/\bar{d}_{0}=60,\mu_{\rm w}=0.46,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=1.5\\}$; (3) a narrower channel and higher flow rate case
$\\{W/\bar{d}_{0}=50,\mu_{\rm w}=0.59,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=1.5\\}$; (4) a more large grains and higher flow rate case
$\\{W/\bar{d}_{0}=60,\mu_{\rm w}=0.58,c^{\rm l}_{0}=0.75,d^{\rm l}/d^{\rm
s}=1.5\\}$; and (5) a larger size ratio and higher flow rate case
$\\{W/\bar{d}_{0}=60,\mu_{\rm w}=0.58,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=2.0\\}$. Each case involves $\sim 20000$ flowing grains. After a
sufficiently long simulation time, a quasi-steady state is attained in each
case, implying the flux balance (4.2). The quasi-steady field quantities
appearing in (4.2) are extracted for 1000 snapshots in time using the coarse-
graining techniques described in Appendix A.2 and then arithmetically averaged
in time. The calculated quasi-steady diffusion flux $C_{\rm
diff}\bar{d}^{2}\dot{\gamma}({\partial c^{\rm l}}/{\partial x})$ at discrete
$x$-positions for each of the five cases is plotted versus the calculated
quantity $\bar{d}^{2}c^{\rm l}(1-c^{\rm l})({\partial\dot{\gamma}}/{\partial
x})$ in Fig. 5(b), and we observe a linear relation. Therefore, the form of
the constitutive equation for the segregation flux (4.2) is also applicable to
dense, bidisperse systems of spheres. The numerical value of the dimensionless
material parameter $C^{\rm S}_{\rm seg}$ for spheres is determined from the
slope of the linear relation to be $C^{\rm S}_{\rm seg}=0.08$.

###  4.2 Annular shear flow

The constitutive equation for the segregation flux (2.17) and the fitted
values of the material parameters should be general across different flow
geometries. To test this for the case of disks, we apply the same process
described in the preceding section for vertical chute flow to a different flow
geometry–annular shear flow. In this flow geometry, flow is driven through
motion of the boundary rather than by gravity, but as in vertical chute flow,
the pressure field is spatially uniform, which eliminates hydrostatic pressure
gradients so that only shear-strain-rate-gradient-driven size-segregation
occurs.

Our DEM simulations of annular shear flow of a dense, bidisperse granular
mixture of disks follow the procedures utilized in prior works in the
literature for monodisperse, frictional disks (Koval et al., 2009; Kamrin and
Koval, 2012, 2014; Liu and Henann, 2018). Consider a dense, bidisperse
granular mixture in a two-dimensional annular shear cell with rough circular
walls of inner radius $R$ and outer radius $R_{\rm o}$, as shown in Fig. 6(a)
for the case of $R=60\bar{d}_{0}$. The inner and outer walls consist of rings
of glued large grains, denoted as black in Fig. 6(a). The circumferential
velocity of the inner wall is prescribed to be $v_{\rm w}$, and its radial
position is fixed. The outer wall does not rotate, and its radius $R_{\rm o}$
fluctuates slightly so as to maintain a constant imposed compressive normal
stress $P_{\rm w}$, utilizing the wall-position control method used throughout
this work (Koval et al., 2009). As in Liu and Henann (2018), we simulate the
full shear cell, as shown in Fig. 6(a), instead of applying periodic boundary
conditions along the circumferential direction to a slice (Koval et al., 2009;
Kamrin and Koval, 2012, 2014). In this flow geometry, all fields are
axisymmetric (i.e., invariant along the $\theta$-direction), and the only non-
zero component of the velocity field is the circumferential component
$v_{\theta}$. Moreover, flow tends to localize near the inner wall with
$v_{\theta}$ rapidly decaying with radial position, as qualitatively
illustrated by the steady velocity field sketched in Fig. 6(a). We choose the
outer radius $R_{\rm o}$ to be sufficiently large so that it does not affect
the resulting flow and segregation fields, and based on our experience, we
take $R_{\rm o}=2R$. Therefore, the role of the outer wall is simply to apply
a far-field pressure, and otherwise, it does not affect the flow and
segregation fields.

Figure 6: (a) Initial well-mixed configuration for two-dimensional DEM
simulation of bidisperse annular shear flow with $40108$ flowing grains. The
inner wall radius is $R=60\bar{d}_{0}$, and the outer wall radius is $R_{\rm
o}=2R$. The inner and outer walls consist of rings of glued large grains,
denoted as black. (b) Segregated configuration after flowing for a total
simulation time of $\tilde{t}=t/\left(R/v_{\rm w}\right)=584$. (c)
Spatiotemporal evolution of the large grain concentration field. Spatial
profiles of (d) the concentration field $c^{\rm l}$ and (e) the normalized
circumferential velocity field ${v}_{\theta}/v_{\rm w}$ at three times
($\tilde{t}=5$, $50$ and $500$) as indicated by the horizontal lines in (c).

Regarding the stress field, in all of our DEM simulations of annular shear
flow of bidisperse disks, we observe that the normal stresses are
approximately equal, i.e., $\sigma_{rr}\approx\sigma_{\theta\theta}$, and
spatially uniform, so that the force balance along the $r$-direction gives
that the pressure field is $P(r)=-\sigma_{rr}(r)=P_{\rm w}$. Since the
pressure field is spatially uniform, all segregation in annular shear flow is
due to shear-strain-rate-gradients. The moment balance gives that the
equivalent shear stress field is
$\tau(r)=|\sigma_{r\theta}(r)|=|\sigma_{\theta r}(r)|=\tau_{\rm w}(R/r)^{2}$,
where $\tau_{\rm w}$ is the inner-wall shear stress. It is important to note
that $\tau_{\rm w}$ is not directly prescribed in our DEM simulations.
Instead, the inner-wall velocity $v_{\rm w}$ is prescribed, and $\tau_{\rm w}$
arises as a result. The stress ratio field is then

$\mu(r)=\mu_{\rm w}(R/r)^{2},$ (4.3)

where $\mu_{\rm w}=\tau_{\rm w}/P_{\rm w}$ is the maximum value of $\mu$,
occurring at the inner wall $(r=R)$.

As for vertical chute flow, there are four important dimensionless parameters
that specify the geometry, loads, and initial conditions for a given case of
annular shear flow of dense, bidisperse granular mixtures: (1)
$R/\bar{d}_{0}$, the dimensionless inner-wall radius; (2) $\tilde{v}_{\rm
w}=(v_{\rm w}/R)\sqrt{\pi\rho_{\rm s}\bar{d_{0}}^{2}/(4P_{\rm w})}$, the
dimensionless inner-wall velocity, which determines $\tau_{\rm w}$ and hence
$\mu_{\rm w}$; (3) $c^{\rm l}_{0}(r)$, the initial large-grain concentration
field; and (4) $d^{\rm l}/d^{\rm s}$, the bidisperse grain-size ratio. We
choose a representative base case of annular shear flow identified by the
parameter set $\\{R/\bar{d}_{0}=60,\tilde{v}_{\rm w}=0.01,c^{\rm
l}_{0}=0.5,d^{\rm l}/d^{\rm s}=1.5\\}$. The well-mixed initial configuration
for the base-case DEM simulation is shown in Fig. 6(a), and the segregated
configuration after driving flow for a total simulation time of
$\tilde{t}=t/\left(R/v_{\rm w}\right)=584$ is shown in Fig. 6(b). The large,
dark-gray grains segregate into a ring near the inner wall, while the small,
light-gray grains form a band just outside this region. Outside of these
bands, where the shear strain-rate is very small, the large and small grains
remain well-mixed. Contours of the spatiotemporal evolution of the coarse-
grained concentration field $c^{\rm l}$ are plotted in Figs. 6(c),
illustrating the time-evolution of $c^{\rm l}$ field. Spatial profiles of the
concentration and velocity fields at three selected snapshots during the
segregation dynamics ($\tilde{t}=t/\left(R/v_{\rm w}\right)=5,50,$ and 500 as
indicated by the dashed lines in Fig. 6(c)) are shown in Figs. 6(d) and (e).
The radial $c^{\rm l}$ profiles in Fig. 6(d) demonstrate the formation of
large-grain-rich and small-grain-rich regions with a persistent well-mixed
far-field. The normalized velocity fields in Fig. 6(e) demonstrate that the
flowing zone is localized near the inner wall with slow creeping flow observed
far from the wall. As in the case of vertical chute flow, the velocity field
quickly develops into a steady flow field, while the large grain concentration
field $c^{\rm l}$ evolves over a longer time-scale before approaching a quasi-
steady state near the end of our simulated time window.

Figure 7: Collapse of $C_{\rm diff}\bar{d}^{2}\dot{\gamma}({\partial c^{\rm
l}}/{\partial r})$ versus $\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})({\partial\dot{\gamma}}/{\partial r})$ for several cases of annular shear
flow of bidisperse disks. Symbols represent coarse-grained, quasi-steady DEM
field data, and the solid line is the best linear fit using $C^{\rm S}_{\rm
seg}=0.23$.

Again, at long times, near the end of the simulated time window
($\tilde{t}=t/\left(R/v_{\rm w}\right)\gtrsim 500$), the concentration field
evolves very slowly, which we identify as the quasi-steady regime. In this
regime, the segregation and diffusion fluxes approximately balance at each
$r$-position, implying that

$C_{\rm diff}\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm l}}{\partial
r}\approx C^{\rm S}_{\rm seg}\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})\dfrac{\partial\dot{\gamma}}{\partial r}.$ (4.4)

As for vertical chute flow, we spatially coarse-grain the DEM data to obtain
the $c^{\rm l}$ and $v_{\theta}$ fields for $144$ evenly-distributed snapshots
in time in the quasi-steady regime ($\tilde{t}\gtrsim 500$), which are
arithmetically averaged in time to obtain the quasi-steady $c^{\rm l}(r)$ and
$v_{\theta}(r)$ fields and then spatially differentiated to obtain the
remaining field quantities in (4.4). Next, we plot $C_{\rm
diff}\bar{d}^{2}\dot{\gamma}({\partial c^{\rm l}}/{\partial r})$ versus
$\bar{d}^{2}c^{\rm l}(1-c^{\rm l})({\partial\dot{\gamma}}/{\partial r})$ in
Fig. 7 with each point representing a unique $r$-position. Finally, this
process is repeated for four additional cases of annular shear flow: (1) a
lower inner-wall velocity $\\{R/\bar{d}_{0}=60,\tilde{v}_{\rm w}=0.001,c^{\rm
l}_{0}=0.5,d^{\rm l}/d^{\rm s}=1.5\\}$; (2) a smaller inner-wall radius
$\\{R/\bar{d}_{0}=40,\tilde{v}_{\rm w}=0.01,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=1.5\\}$; (3) more large grains $\\{R/\bar{d}_{0}=60,\tilde{v}_{\rm
w}=0.01,c^{\rm l}_{0}=0.75,d^{\rm l}/d^{\rm s}=1.5\\}$; and (4) a larger size
ratio $\\{R/\bar{d}_{0}=60,\tilde{v}_{\rm w}=0.01,c^{\rm l}_{0}=0.5,d^{\rm
l}/d^{\rm s}=3.0\\}$, and the coarse-grained, quasi-steady fields are included
in the data plotted in Fig. 7. Collectively, we observe a strong collapse to a
linear relation. Crucially, the slope of the linear relation in Fig. 7
(indicated by the solid line) gives the same value for the dimensionless
material parameter $C_{\rm seg}^{\rm S}$ obtained by fitting to vertical chute
flow data, $C_{\rm seg}^{\rm S}=0.23$. This observation of agreement between
the best fit values of $C_{\rm seg}^{\rm S}$ obtained using two different flow
geometries provides support for our choice of the constitutive equation for
the segregation flux (2.17) and the fitted value of $C^{\rm S}_{\rm seg}$ for
disks. Having established that the parameter $C^{\rm S}_{\rm seg}$ is
independent of the flow geometry, driving conditions, and initial conditions,
we henceforth regard $C^{\rm S}_{\rm seg}$ as a material parameter for a given
dense granular system, analogous to how rheological material parameters such
as $\mu_{\rm s}$ are regarded. Of course, the values of $C^{\rm S}_{\rm seg}$
determined above for disks and spheres likely depend on grain interaction
properties, such as the inter-particle friction coefficient, but elucidating
this dependence is beyond the scope of the present work.

##  5 Validation of the continuum model in the transient regime

In the preceding section, we only used DEM data from the quasi-steady regime
to test the constitutive equation for the shear-strain-rate-gradient-driven
segregation flux (2.17) and to determine the material parameter $C_{\rm
seg}^{\rm S}$. In this section, we compare continuum model predictions of the
transient evolution of segregation and flow fields to the DEM data for both
vertical chute flow and annular shear flow as a validation test of the model.
To obtain continuum model predictions in the transient regime, we couple the
segregation dynamics equation (2.18) with the NGF model, (2.6) and (2.7), and
a use fixed sets of material parameters for disks,

$\\{\mu_{\rm s}=0.272,b=1.168,A=0.9,C_{\rm diff}=0.20,C^{\rm S}_{\rm
seg}=0.23\\},$ (5.1)

and for spheres,

$\\{\mu_{\rm s}=0.37,\mu_{2}=0.95,I_{0}=0.58,A=0.43,C_{\rm diff}=0.045,C^{\rm
S}_{\rm seg}=0.08\\}.$ (5.2)

###  5.1 Vertical chute flow

First, we describe in detail how the continuum model is solved to obtain
predictions for the transient evolution of segregation and flow fields for the
case of vertical chute flow of disks. In vertical chute flow, the stress field
may be straightforwardly deduced from a static force balance, giving that the
pressure field is uniform, $P(x)=P_{\rm w}$, and that the stress ratio field
$\mu(x)$ is given through (4.1), and therefore, the balance of linear momentum
(2.2) is satisfied and does not further enter the solution procedure.
Continuum model predictions are obtained by numerically solving the remaining
governing equations using finite-differences. Summarizing the coupled
boundary/initial-value problem for flow and segregation in the context of
vertical chute flow, the unknown fields are the velocity field $v_{z}(x,t)$
and the accompanying strain-rate field $\dot{\gamma}(x,t)=\partial
v_{z}/\partial x$, the granular fluidity field $g(x,t)$, and the large-grain
concentration field $c^{\rm l}(x,t)$. The governing equations are (1) the flow
rule (2.6)

$\dot{\gamma}=g\mu,$ (5.3)

(2) the nonlocal rheology (2.7)

$g=g_{\rm loc}(\mu,P_{\rm w})+\xi^{2}(\mu)\dfrac{\partial^{2}g}{\partial
x^{2}}$ (5.4)

with $g_{\rm loc}$ and $\xi$ given through (2.9) and (2.13)1, respectively,
and (3) the segregation dynamics equation (2.18)

$\dfrac{\partial c^{\rm l}}{\partial t}+\dfrac{\partial}{\partial
x}\left(-C_{\rm diff}\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm
l}}{\partial x}+C^{\rm S}_{\rm seg}\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})\dfrac{\partial\dot{\gamma}}{\partial x}\right)=0,$ (5.5)

where $\bar{d}=c^{\rm l}d^{\rm l}+(1-c^{\rm l})d^{\rm s}$.

Regarding boundary conditions, we impose Dirichlet fluidity boundary
conditions at the walls (i.e., $g=g_{\rm loc}(\mu_{\rm w},P_{\rm w})$ at
$x=\pm W/2$) as well as no flux boundary conditions at the walls (i.e.,
$w_{x}^{\rm l}=-C_{\rm diff}\bar{d}^{2}\dot{\gamma}(\partial c^{\rm
l}/\partial x)+C^{\rm S}_{\rm seg}\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})(\partial\dot{\gamma}/\partial x)=0$ at $x=\pm W/2$). Due to the time
derivative in (5.5), an initial condition for the concentration field $c^{\rm
l}_{0}(x)=c^{\rm l}(x,t=0)$ is required. In order to account for the
concentration fluctuations inherent in the initial state, we obtain the
coarse-grained $c^{\rm l}$-field from the initial DEM configuration for each
case and utilize this field as the initial condition field $c_{0}^{\rm l}(x)$
in each of the respective continuum simulations.

Then, for a given case identified through a set of input parameters
$\\{W/\bar{d}_{0},$ $\mu_{\rm w},$ $c^{\rm l}_{0}(x),$ $d^{\rm l}/d^{\rm
s}\\}$, we obtain numerical predictions of the continuum model utilizing
finite differences as follows. First, at a given point in time, the
concentration field $c^{\rm l}(x)$ is known, allowing the average grain-size
field to be calculated through $\bar{d}=c^{\rm l}d^{\rm l}+(1-c^{\rm l})d^{\rm
s}$. Using the stress ratio field $\mu(x)$ for vertical chute flow (4.1), the
local fluidity $g_{\rm loc}(\mu,P)$ and the cooperativity length $\xi(\mu)$
(equations (2.9) and (2.13)1) may be calculated at each spatial grid point.
Then, the nonlocal rheology (5.4) may be used to solve for the fluidity field
$g(x)$ at the current step, using central differences in space. The strain-
rate field follows using (5.3), which may be integrated to obtain the velocity
field $v_{z}(x)$. Next, (5.5) is used to determine the concentration field at
the next time step utilizing the forward Euler method and central differences
in space with one modification–the spatial derivatives of $c^{\rm l}$
appearing in the diffusion flux term in (5.5) are treated implicitly in order
to improve numerical stability. This completes one time step, and this process
is repeated to step forward in time and calculate the transient evolution of
the concentration and flow fields, $c^{\rm l}(x,t)$ and $v_{z}(x,t)$. In our
finite-difference calculations, we utilize a fine spatial resolution of
$\Delta x\ll\bar{d}_{0}$, and we have verified that the time-step is
sufficiently small in order to ensure stable, accurate results.

Figure 8: Comparisons of continuum model predictions with corresponding DEM
simulation results for the transient evolution of the segregation dynamics for
three cases of vertical chute flow of disks: (a) Base case
$\\{W/\bar{d}_{0}=60,\mu_{\rm w}=0.45,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=1.5\\}$; (b) Lower flow rate case $\\{W/\bar{d}_{0}=60,\mu_{\rm
w}=0.375,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm s}=1.5\\}$; and (c) Narrower chute
width case $\\{W/\bar{d}_{0}=40,\mu_{\rm w}=0.45,c^{\rm l}_{0}=0.5,d^{\rm
l}/d^{\rm s}=1.5\\}$. Additional cases are shown in Fig. 9. For each case, the
first column shows spatiotemporal contours of the evolution of $c^{\rm l}$
measured in the DEM simulations. The second column shows comparisons of the
DEM simulations (solid black lines) and continuum model predictions (dashed
gray lines) of the $c^{\rm l}$ field at four time snapshots representing
different stages of the segregation process: $\tilde{t}=4\times 10^{3}$,
$2\times 10^{4}$, $1\times 10^{5}$, and $4\times 10^{5}$ in the sequence of
top left, top right, bottom left, bottom right. The third column shows
comparisons of the quasi-steady, normalized velocity profiles at
$\tilde{t}=4\times 10^{5}$ from DEM simulations and continuum model
predictions. Figure 9: Comparisons of continuum model predictions with
corresponding DEM simulation results for the transient evolution of the
segregation dynamics for two cases of vertical chute flow of disks: (a) More
large grains case $\\{W/\bar{d}_{0}=60,\mu_{\rm w}=0.45,c^{\rm
l}_{0}=0.75,d^{\rm l}/d^{\rm s}=1.5\\}$ and (b) Larger size ratio case
$\\{W/\bar{d}_{0}=60,\mu_{\rm w}=0.45,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=3.0\\}$. Additional cases are shown in Fig. 8. Results are organized as
described in the caption of Fig. 8. Figure 10: Comparisons of continuum model
predictions with corresponding DEM simulation results for the transient
evolution of the segregation dynamics for three cases of vertical chute flow
of spheres: (a) Base case $\\{W/\bar{d}_{0}=60,\,\mu_{\rm w}=0.51,\,c^{\rm
l}_{0}=0.5,\,d^{\rm l}/d^{\rm s}=1.5\\}$; (b) More large grains and higher
flow rate case $\\{W/\bar{d}_{0}=60,\mu_{\rm w}=0.58,c^{\rm l}_{0}=0.75,d^{\rm
l}/d^{\rm s}=1.5\\}$; and (c) Larger grain-size ratio and higher flow rate
case $\\{W/\bar{d}_{0}=60,\,\mu_{\rm w}=0.58,\,c^{\rm l}_{0}=0.5,\,d^{\rm
l}/d^{\rm s}=2.0\\}$. Results are organized as described in the caption of
Fig. 8.

We compare predictions of the continuum model against DEM data for all five
cases of vertical chute flow considered in Section 4.1 in order to test the
generality of the model. Figures 8 and 9 summarize the comparisons for these
five cases. The first column of Figs. 8 and 9 shows the spatiotemporal
contours of the evolution of the $c^{\rm l}$ field measured in the DEM
simulations for each case, and the second column shows comparisons of the DEM
simulations (solid black lines) and the continuum predictions (dashed gray
lines) for the $c^{\rm l}$ field at four snapshots in time
($\tilde{t}=t/\left(d^{\rm s}\sqrt{\rho_{\rm s}/P_{\rm w}}\right)=4\times
10^{3}$, $2\times 10^{4}$, $1\times 10^{5}$ and $4\times 10^{5}$) indicated by
the horizontal lines in the first column of Figs. 8 and 9. Based on Figs. 8
and 9, the coupled model generally does a good job capturing the salient
features of the evolution of the $c^{\rm l}$ field across all cases. For
instance, for the narrower chute case shown in Fig. 8(c), the segregation
process nearly completes within the simulated time window with the mixed core
along the center of the chute nearly disappearing, and the continuum model
prediction captures this observation well.

Regarding flow fields, comparisons of the quasi-steady, normalized velocity
fields at $\tilde{t}=4\times 10^{5}$ from the DEM simulations and the
continuum model predictions are shown in the third column of Fig. 8. Since the
velocity field evolves minimally during the segregation process, only the flow
field at long time is shown. We note that the velocity field is very well-
predicted in all cases, including the creeping regions far from the wall. This
favorable comparison provides support of our generalization of the NGF model
to bidisperse granular systems discussed in Section 2.3–in particular, the
choices to use the average grain size $\bar{d}$ in the expression for the
cooperativity length (2.13) and to continue to use the numerical value for the
nonlocal amplitude determined for monodisperse systems, $A=0.9$, without
refitting.

Next, we make comparisons between continuum model predictions and DEM data for
vertical chute flow of bidisperse spheres. The governing equations and
boundary conditions are same as those used above for bidisperse disks. That
is, the fluidity field is governed by the nonlocal rheology (5.4) with
Dirichlet fluidity boundary conditions at the walls ($g=g_{\rm loc}(\mu_{\rm
w},P_{\rm w})$ at $x=\pm W/2$), and the concentration field is governed by the
segregation dynamics equation (5.5) with no flux boundary conditions at the
walls ($w_{x}^{\rm l}=0$ at $x=\pm W/2$). The only difference from the process
described above for bidisperse disks is that the values $P_{\rm w}$ and
$\mu_{\rm w}$ used in the continuum simulations are obtained from the coarse-
grained stress fields in the DEM data for each case, rather than based on the
nominal value of the compressive wall stress $P_{\rm w}$ applied in the DEM
simulation. This is done to account for the normal stress differences that
arise for spheres, which slightly affect the predicted velocity fields and
hence the consequent concentration fields. We consider three different cases:
(1) the base case $\\{W/\bar{d}_{0}=60,\,\mu_{\rm w}=0.51,\,c^{\rm
l}_{0}=0.5,\,d^{\rm l}/d^{\rm s}=1.5\\}$, (2) the more large grains and higher
flow rate case $\\{W/\bar{d}_{0}=60,\,\mu_{\rm w}=0.58,\,c^{\rm
l}_{0}=0.75,\,d^{\rm l}/d^{\rm s}=1.5\\}$, and (3) the larger size ratio and
higher flow rate case $\\{W/\bar{d}_{0}=60,\,\mu_{\rm
w}=0.58,\,c^{l}_{0}=0.5,\,d^{l}/d^{s}=2.0\\}$, which are shown in Figs. 10(a),
(b), and (c), respectively. The leftmost column of Fig. 10 shows the
spatiotemporal contours of the evolution of the $c^{\rm l}$ field from the DEM
data. The middle column shows comparisons between the DEM simulations (solid
black lines) and the continuum model predictions (dashed gray lines) for the
$c^{\rm l}$ fields at four time snapshots in time ($\tilde{t}=t/\left(d^{\rm
s}\sqrt{\rho_{\rm s}/P_{\rm w}}\right)=5\times 10^{2}$, $2\times 10^{3}$,
$1\times 10^{4}$ and $4.5\times 10^{4}$) indicated by the horizontal lines in
the first column of Fig. 10. Lastly, comparisons of the quasi-steady,
normalized velocity fields at $\tilde{t}=4.5\times 10^{4}$ from the DEM
simulations and the continuum model predictions are shown in the third column
of Fig. 10. The continuum model is able to capture the decaying velocity field
quite well in all cases, and therefore, our choice to continue using the value
of the nonlocal amplitude estimated for monodisperse systems, $A=0.43$,
without readjustment works well for spheres. Overall, the coupled continuum
model is capable of quantitatively predicting both the flow fields and the
transient evolution of the segregation dynamics in vertical chute flow of
bidisperse spheres.

Figure 11: Comparisons of the continuum model predictions with corresponding
DEM simulation results for the transient evolution of the segregation dynamics
for three cases of annular shear flow of disks: (a) Base case
$\\{R/\bar{d}_{0}=60,\tilde{v}_{\rm w}=0.01,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=1.5\\}$; (b) Lower inner-wall velocity case
$\\{R/\bar{d}_{0}=60,\tilde{v}_{\rm w}=0.001,c^{\rm l}_{0}=0.5,d^{\rm
l}/d^{\rm s}=1.5\\}$; and (c) Smaller annular shear cell case
$\\{R/\bar{d}_{0}=40,\tilde{v}_{\rm w}=0.01,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=1.5\\}$. Additional cases are shown in Fig. 12. For each case, the first
column shows spatiotemporal contours of the evolution of $c^{\rm l}$ measured
in the DEM simulations. The second column shows comparisons of the DEM
simulations (solid black lines) and continuum model predictions (dashed gray
lines) of the $c^{\rm l}$ field at four time snapshots representing different
stages of the segregation process: $\tilde{t}=5$, $50$, $200$, and $550$ in
the sequence of top left, top right, bottom left, bottom right. The third
column shows comparisons of the quasi-steady, normalized velocity profiles at
$\tilde{t}=550$ from DEM simulations and continuum model predictions. Figure
12: Comparisons of the continuum model predictions with corresponding DEM
simulation results for the transient evolution of the segregation dynamics for
two cases of annular shear flow of disks: (a) More large grains case
$\\{R/\bar{d}_{0}=60,\tilde{v}_{\rm w}=0.01,c^{\rm l}_{0}=0.75,d^{\rm
l}/d^{\rm s}=1.5\\}$ and (b) Larger size ratio case
$\\{R/\bar{d}_{0}=60,\tilde{v}_{\rm w}=0.01,c^{\rm l}_{0}=0.5,d^{\rm l}/d^{\rm
s}=3.0\\}$. Additional cases are shown in Fig. 11. Results are organized as
described in the caption of Fig. 11.

###  5.2 Annular shear flow

We utilize an analogous process to obtain continuum model predictions for the
transient evolution of concentration and flow fields in annular shear flow of
disks. In annular shear, based on the static force and moment balances, the
pressure field is uniform, $P(r)=P_{\rm w}$, and the stress ratio field is
given by (4.3). The governing equations (5.4) and (5.5) are modified to
appropriately account for the divergence and Laplacian operators in
cylindrical coordinates. Also, since $v_{\rm w}$, not $\mu_{\rm w}$, is
specified in our DEM simulations of annular shear, while $\mu_{\rm w}$ is
specified in our continuum simulations, we iteratively adjust the value of
$\mu_{\rm w}$ input into our continuum simulations in order to achieve the
target value of $v_{\rm w}$ in the predicted quasi-steady flow field.
Otherwise, our process for obtaining numerical predictions from the continuum
model is the same. Dirichlet fluidity boundary conditions and no flux boundary
conditions are imposed at the walls, and the initial concentration field is
extracted from the initial DEM configuration for each case.

Then, we compare continuum model predictions against DEM data for the five
cases of annular shear flow discussed in Section 4.2 in order to further
validate the model. Figures 11 and 12 summarize the comparisons for these
fives cases and are organized in the same manner as Figs. 8 and 9. Again, the
coupled, continuum model does a good job capturing the segregation dynamics
and its dependence on the input parameters. Moreover, the quasi-steady
velocity fields are well-predicted by the NGF model in all cases, including
the creeping region far from the inner wall. We reiterate that all continuum
model predictions are obtained using the same set of material parameters for
disks (5.1).

##  6 Discussion and Conclusion

In this paper, we studied coupled size-segregation and flow in dense,
bidisperse granular systems of disks and spheres and developed a
phenomenological continuum model that captures the simultaneous evolution of
both segregation and flow fields. We focused on the shear-strain-rate-
gradient-driven size-segregation mechanism in two configurations in which the
pressure field is uniform–vertical chute flow and annular shear flow–and based
on observations from DEM simulations, we proposed a phenomenological
constitutive equation for the shear-strain-rate-gradient-driven flux. When
combined with a standard model for granular diffusion, the segregation model
involves two dimensionless parameters $\\{C_{\rm diff},C^{\rm S}_{\rm
seg}\\}$, which multiply the two fluxes appearing in the model–the diffusion
and shear-strain-rate-gradient-driven fluxes, respectively. By coupling the
segregation model with the NGF model adapted to bidisperse systems, we may
quantitatively predict both the flow fields and the segregation dynamics for
dense flows of bidisperse disks and spheres for two distinct flow geometries
and under a number of different flow conditions.

Size-segregation in granular materials is a complex and rich problem, so there
remain many avenues for model improvement and unresolved research questions to
be answered. One important question relates to the constitutive equation for
the shear-strain-rate-gradient-driven segregation flux (2.17). Although our
use of a constitutive equation driven by gradients in $\dot{\gamma}$ does a
good job capturing the DEM data, there are other theories in the literature
based on gradients of other field quantities. In particular, Hill and
coworkers (Fan and Hill, 2011b; Hill and Tan, 2014) have proposed that
gradients in the kinetic stress, which is related to velocity fluctuations and
hence the granular temperature, drive segregation. Since Zhang and Kamrin
(2017) have established a connection between velocity fluctuations and the
granular fluidity $g$, it is possible to propose other forms for the
constitutive equation for $w_{i}^{\rm seg}$ based on gradients in $g$. For
example, instead of (2.17), consider the following form for the segregation
flux:

$w_{i}^{\rm seg}=C^{\rm S}_{\rm seg}\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})\dfrac{\partial g}{\partial x_{i}}.$ (6.1)

Then, applying the quasi-steady flux balance condition,

$C_{\rm diff}\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm l}}{\partial
x_{i}}\approx C^{\rm S}_{\rm seg}\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})\dfrac{\partial g}{\partial x_{i}},$ (6.2)

to the quasi-steady DEM data for vertical chute flow and annular shear flow of
disks, we obtain the collapses shown in Figs. 13(a) and (b),
respectively.333The coarse-grained values of $g$ and its gradient are obtained
using the coarse-grained values of $\dot{\gamma}$ and its gradient, calculated
as described in Appendix A.2, along with $\mu$ and its gradient, calculated
using (4.1) and not by coarse-graining. The solid lines represent the best
linear fit using $C^{\rm S}_{\rm seg}=0.08$. The collapses are reasonable but
not as strong as those shown in Figs. 5(a) and 7 for a segregation flux based
on gradients in the shear-strain-rate, leading us to choose to work with the
constitutive equation (2.17) on pragmatic grounds.

Figure 13: (a) Collapse of $C_{\rm diff}\bar{d}^{2}\dot{\gamma}({\partial
c^{\rm l}}/{\partial x})$ versus $\bar{d}^{2}c^{\rm l}(1-c^{\rm l})({\partial
g}/{\partial x})$ for several cases of vertical chute flow and (b) collapse
$C_{\rm diff}\bar{d}^{2}\dot{\gamma}({\partial c^{\rm l}}/{\partial r})$
versus $\bar{d}^{2}c^{\rm l}(1-c^{\rm l})({\partial g}/{\partial r})$ for
several cases of annular shear flow of bidisperse disks. (c) Collapse of
$C_{\rm diff}\bar{d}^{2}\dot{\gamma}({\partial c^{\rm l}}/{\partial x})$
versus $\sqrt{\rho_{\rm s}/P}\bar{d}c^{\rm l}(1-c^{\rm l})({\partial
T}/{\partial x})$ for several cases of vertical chute flow and (d) collapse of
$C_{\rm diff}\bar{d}^{2}\dot{\gamma}({\partial c^{\rm l}}/{\partial x})$
versus $\sqrt{\rho_{\rm s}/P}\bar{d}c^{\rm l}(1-c^{\rm l})({\partial
T}/{\partial x})$ for several cases of annular shear flow of bidisperse disks.
Symbols represent coarse-grained, quasi-steady DEM field data, and the solid
lines represent the best linear fit using $C^{\rm S}_{\rm seg}=0.08$ for (a)
and (b) and $C^{\rm S}_{\rm seg}=0.7$ for (c) and (d).

Additionally, we have tested a possible constitutive equation for the
segregation flux driven by gradients in the granular temperature, defined as
$T=(\delta v)^{2}$, where $\delta v$ is the velocity fluctuation.444The
definitions of the granular temperature $T$ and the velocity fluctuation
$\delta v$ as well as the coarse-graining method used to obtain these
quantities follow Zhang and Kamrin (2017). They are treated as field
quantities just as $\dot{\gamma}$. Following the coarse-graining process
described in (A.2), the instantaneous $T$ field at a bin is the grain-area-
weighted summation of the square of the difference of the grain instantaneous
velocity vector and the coarse-grained instantaneous velocity field at that
grain center, for all grains that are intersected by the bin. Then, consider
the following form for the segregation flux:

$w_{i}^{\rm seg}=C^{\rm S}_{\rm seg}\sqrt{\rho_{\rm s}/P}\bar{d}c^{\rm
l}(1-c^{\rm l})\dfrac{\partial T}{\partial x_{i}},$ (6.3)

where the inertial time $\sqrt{\rho_{\rm s}/P}\bar{d}$ is included in the
prefactor for dimensional reasons. As above, upon applying the quasi-steady
flux balance condition,

$C_{\rm diff}\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm l}}{\partial
x_{i}}\approx C^{\rm S}_{\rm seg}\sqrt{\rho_{\rm s}/P}\bar{d}c^{\rm
l}(1-c^{\rm l})\dfrac{\partial T}{\partial x_{i}},$ (6.4)

to the quasi-steady DEM data, Figs. 13(c) and (d) show the collapses for
granular-temperature-gradient-driven segregation for vertical chute flow and
annular shear flow of disks, respectively. In this case, the solid lines
represent the best linear fit using $C^{\rm S}_{\rm seg}=0.7$. The collapse is
quite good; however, in order to utilize the constitutive equation (6.3) in
practice, an additional constitutive equation that gives the granular
temperature field in terms of other continuum quantities–such as strain-rate,
stress, and fluidity–is needed. Recent work by Kim and Kamrin (2020) offers a
path forward on this point. In summary, we acknowledge that the possibility
for an alternative form for the segregation flux $w^{\rm seg}_{i}$ remains and
that future work involving additional flow geometries will be required to
conclusively judge which constitutive equation is the most predictive.

Another question regarding the constitutive equation for the segregation flux
is the $c^{\rm l}$ dependence of the pre-factor. We simply use a symmetric
dependence $c^{\rm l}(1-c^{\rm l})$, while other researchers (e.g., van der
Vaart et al., 2015; Tunuguntla et al., 2017) invoke expressions that depend on
both $c^{\rm l}$ and $d^{\rm l}/d^{\rm s}$ in a more complex manner.

Finally, in this paper, we have focused on two simple flow geometries that
have two important features: (1) the continuum fields are one-dimensional,
only varying along one spatial direction, and (2) the pressure field is
spatially uniform. In order to apply the proposed continuum model in more
complex flow geometries, such as heap flows or split-bottom flow, two
important steps are necessary. First, this paper solely considered shear-
strain-rate-driven size-segregation. Now that a predictive continuum model for
this mechanism has been established, it remains to return to pressure-
gradient-driven size-segregation in order to incorporate this mechanism by
introducing an additional flux contribution to (2.14) to obtain a more general
model. Second, a robust numerical implementation of the complex system of
coupled equations that is capable of addressing problems in general
geometries, involving multi-dimensional continuum fields, is needed. One
possible approach is to utilize the finite-element method, as in previous work
involving the NGF model (Henann and Kamrin, 2016). These steps will be
addressed in future works.

## Acknowledgements

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

##  Appendix A Discrete-element method simulations and coarse-graining
procedures

###  A.1 Simulated granular systems

We consider two types of simulated granular systems: two-dimensional systems
consisting of a dense collection of circular disks and three-dimensional
systems consisting of a dense collection of spheres. We consider bidisperse
systems and denote the mean diameter of the large particles as $d^{\rm l}$ and
the mean diameter of the small particles as $d^{\rm s}$ for both types. For
both disks and spheres, we take $d^{\rm l}=3$ mm and $d^{\rm s}=2$ mm for the
base case, so that $d_{\rm l}/d_{\rm s}=1.5$. For both large and small
particles, the diameters of individual particles are chosen from a uniform
distribution over the range of $\pm 10\%$ of the respective mean diameters to
prevent crystallization. In the two-dimensional granular system, $\rho_{\rm
s}$ denotes the grain-material area-density, which we take to be
$\rho_{s}=3.26\,{\rm kg/m^{2}}$ for both large and small disks, and in the
three-dimensional granular system, $\rho_{\rm s}$ denotes the grain-material
volume-density, which is taken to be $\rho_{\rm s}=2450\,{\rm kg/m^{3}}$ for
both large and small spheres to eliminate density-based segregation. For two-
dimensional systems, the mass of the large disks is given by $m^{\rm
l}=({\pi}/{4})\rho_{\rm s}(d^{\rm l})^{2}$, and the mass of the small disks is
given as $m^{\rm s}=({\pi}/{4})\rho_{\rm s}(d^{\rm s})^{2}$, so that the
characteristic grain mass is $m=c^{\rm l}_{0}m^{\rm l}+(1-c^{\rm l}_{0})m^{\rm
s}$, where $c^{\rm l}_{0}$ the initial concentration of the large grains.
Similarly, for three-dimensional systems, the mass of the large spheres is
$m^{\rm l}=({\pi}/{6})\rho_{\rm s}(d^{\rm l})^{3}$, and the mass of the small
spheres is $m^{\rm s}=({\pi}/{6})\rho_{\rm s}(d^{\rm s})^{3}$, so that the
characteristic grain mass is $m=c^{\rm l}_{0}m^{\rm l}+(1-c^{\rm l}_{0})m^{\rm
s}$.

For the grain interaction model, the interaction force is given through a
spring/dashpot contact law that accounts for elasticity, damping, and sliding
friction (da Cruz et al., 2005; Koval et al., 2009; Kamrin and Koval, 2014;
Zhang and Kamrin, 2017). The normal contact force $F_{\rm n}$ is given
linearly through the normal component of the contact overlap, denoted by
$\delta_{\rm n}$, with stiffness $k_{\rm n}$ and the relative normal velocity,
denoted by $\dot{\delta}_{\rm n}$, with damping coefficient $g_{\rm n}$ as
$F_{\rm n}=k_{\rm n}\delta_{\rm n}+g_{\rm n}\dot{\delta}_{\rm n}$ whenever
$\delta_{\rm n}\geq 0$ and $F_{\rm n}=0$ whenever $\delta_{\rm n}<0$. The
normal damping coefficient is given by $g_{\rm n}=\sqrt{mk_{\rm n}}(-2\ln
e)/\sqrt{2(\pi^{2}+\ln^{2}e)}$ where $e$ is the coefficient of restitution for
binary collisions and $m$ is the characteristic grain mass discussed above. We
denote the tangential stiffness and damping coefficient as $k_{\rm t}$ and
$g_{\rm t}$, respectively, and take $g_{\rm t}=0$, so that the tangential
force is given as $F_{\rm t}=k_{\rm t}\delta_{\rm t}$ whenever $\delta_{\rm
n}\geq 0$, where $\delta_{t}$ is the tangential component of the contact
displacement. The magnitude of the tangential component of the contact force
is limited by Coulomb friction, which depends on the inter-particle friction
coefficient $\mu_{\rm surf}$. Thus, the parameter set $\\{k_{\rm n},k_{\rm
t},e,\mu_{\rm surf}\\}$ fully describes the interaction properties.
Throughout, the normal stiffness $k_{\rm n}$ is taken to be sufficiently large
so that grains behave as stiff and nearly rigid. For two-dimensional systems,
$k_{\rm n}/P>10^{4}$, and for three-dimensional systems, $k_{\rm
n}/P\bar{d}_{0}>10^{4}$, where $P$ is the characteristic confining pressure
for a given configuration (force per unit length in two dimensions and force
per unit area in three dimensions) and $\bar{d}_{0}=c^{\rm l}_{0}d^{\rm
l}+(1-c^{\rm l}_{0})d^{\rm s}$ is the characteristic grain size. In the stiff
grain regime, the only interaction parameter that significantly affects the
rheology of dense grains is $\mu_{\rm surf}$, which we have kept constant as
$\mu_{\rm surf}=0.4$ throughout. The other parameters, namely the ratio
$k_{\rm n}/k_{\rm t}$ and the restitution coefficient $e$, have negligible
effects on the rheology (Kamrin and Koval, 2014) and thus the segregation
dynamics in the stiff particle regime, so we maintain $k_{\rm t}/k_{\rm
n}=1/2$ and $e=0.1$ throughout. Finally, the open-source software LAMMPS
(Plimpton, 1995) is used to numerically integrate the equations of motion for
each particle, and we restrict the time step for numerical integration to be
$0.01$-$0.1$ of the binary collision time $\tau_{\rm
c}=\sqrt{m\left(\pi^{2}+\ln^{2}e\right)/4k_{\rm n}}$ for stability and
accuracy of the simulation results.

###  A.2 Averaging methods

In this appendix, we describe the spatial and temporal averaging methods
utilized to extract continuum fields from our DEM data. We begin with the
spatial averaging procedure for a given snapshot of DEM data at a time $t$.
All flows considered in this work are one-dimensional, in which the continuum
fields vary only along one direction–i.e., along the $z$-direction in simple
shear flow (Fig. 2), along the $x$-direction in vertical chute flow (Fig. 4),
and along the $r$-direction in annular shear flow (Fig. 6)–and periodic along
all other directions. Here, we describe the procedure for vertical chute flow
of disks in detail, which may be straightforwardly adapted to the other flow
geometries. We utilize a bin-based coarse-graining process, in which we
construct a slender rectangle that spans the simulation domain along the
$z$-direction and is centered at a given $x$-position with a finite width
along the $x$-direction.555For annular shear flow of disks, the bins are thin,
annular rings that are centered at a given $r$-position with a finite
thickness along the $r$-direction. Then, we assign each intersected grain $i$
a weight $A_{i}$, defined as the area of the grain $i$ inside the bin.
Following Tunuguntla et al. (2017) for a bidisperse system, we denote the sets
of large and small grains intersected by the bin as ${\cal F}^{\rm l}$ and
${\cal F}^{\rm s}$, respectively, so that the set of all grains intersected by
the bin is ${\cal F}={\cal F}^{\rm l}\cup{\cal F}^{\rm s}$ with ${\cal F}^{\rm
l}\cap{\cal F}^{\rm s}=\varnothing$. The instantaneous solid area fraction
field for species $\alpha$ is $\phi^{\alpha}(x,t)=(\sum_{i\in{\cal
F}^{\alpha}}A_{i})/A$, where $A$ is the total area of the bin, and the
corresponding concentration field for species $\alpha$ is
$c^{\alpha}(x,t)=\phi^{\alpha}(x,t)/\phi(x,t)$ with $\phi(x,t)=\phi^{\rm
l}(x,t)+\phi^{\rm s}(x,t)$. With the instantaneous velocity of each grain $i$
denoted as ${\bf v}_{i}(t)$, the instantaneous velocity field is ${\bf
v}(x,t)=(\sum_{i\in{\cal F}}A_{i}{\bf v}_{i}(t))/(\sum_{i\in{\cal
F}}A_{i})$.666We note that the definition of the instantaneous velocity field
is consistent with first defining species-specific velocity fields–i.e., ${\bf
v}^{\alpha}(x,t)=(\sum_{i\in{\cal F}^{\alpha}}A_{i}{\bf
v}_{i}(t))/(\sum_{i\in{\cal F}^{\alpha}}A_{i})$–and then calculating the
mixture-level field–i.e., ${\bf v}(x,t)=c^{\rm l}(x,t){\bf v}^{\rm
l}(x,t)+(1-c^{\rm l}(x,t)){\bf v}^{\rm s}(x,t)$. Likewise, with the
instantaneous stress tensor associated with grain $i$ defined as
$\bm{\sigma}_{i}(t)=(\sum_{j\neq i}{\bf r}_{ij}\otimes{\bf f}_{ij})/(\pi
d_{i}^{2}/4)$, 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}(x,t)=(\sum_{i\in{\cal
F}}A_{i}\bm{\sigma}_{i}(t))/A$.

In vertical chute flow of spheres, we utilize a similar spatial coarse-
graining approach. Instead of two-dimensional, rectangular bins, we use three-
dimensional, rectangular-cuboidal bins that span the simulation domain along
the $y$\- and $z$-directions and are centered at a given $x$-position with a
finite width along the $x$-direction. Then, the weight for each grain $i$ is
the volume $V_{i}$ of the grain $i$ inside the bin. The instantaneous solid
volume fraction field for species $\alpha$ is
$\phi^{\alpha}(x,t)=(\sum_{i\in{\cal F}^{\alpha}}V_{i})/V$, where $V$ is the
total volume of the bin; the instantaneous concentration field for species
$\alpha$ is $c^{\alpha}(x,t)=\phi^{\alpha}(x,t)/\phi(x,t)$; the instantaneous
velocity field is ${\bf v}(x,t)=(\sum_{i\in{\cal F}}V_{i}{\bf
v}_{i}(t))/(\sum_{i\in{\cal F}}V_{i})$; 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)$; and the instantaneous stress
field is $\bm{\sigma}(x,t)=(\sum_{i\in{\cal F}}V_{i}\bm{\sigma}_{i}(t))/V$.

Our analysis of the size-segregation process in this paper depends on
obtaining accurate and high-resolution coarse-grained $c^{\rm l}$ fields from
the DEM data. Therefore, the choices of the bin width and the spatial
resolution of the bins are crucial. Throughout, we take a bin width of
$4\bar{d}_{0}$ and a spatial resolution of roughly $0.1\bar{d}_{0}$ for both
disks and spheres. Note that for these choices, adjacent bins overlap. We have
ensured that a bin width of $4\bar{d}_{0}$ is sufficiently small so that the
coarse-grained data is not over-smoothed. Specifically, we have tested that
the coarse-grained velocity and stress fields are insensitive to the choice of
bin width over the range of $0$ to $9.6\bar{d}_{0}$. In the limit that the bin
width goes to zero, the coarse-graining procedure reduces to that utilized
successfully in the literature for the velocity and stress fields for both
disks (e.g., da Cruz et al., 2005; Koval et al., 2009) and spheres (e.g.,
Zhang and Kamrin, 2017; Kim and Kamrin, 2020). However, as shown by Weinhart
et al. (2013), applying a coarse-graining procedure with a small or zero bin
width to the $\phi$ field–and hence the $\phi^{\rm l}$, $\phi^{\rm s}$,
$c^{\rm l}$, and $c^{\rm s}$ fields–will lead to spatial fluctuations due to
particle layering near the walls. We have ensured that a bin width of
$4\bar{d}_{0}$ is sufficiently large so that these layering effects are not
observed in the $c^{\rm l}$ field–i.e., we are within the “plateau range”
(Weinhart et al., 2013) of bin widths that produce bin-width-independent,
coarse-grained continuum fields. We note that when plotting spatiotemporal
contours of the concentration fields, profiles of the concentration fields,
and profiles of the velocity fields (e.g., Fig. 4), we truncate the coarse-
grained DEM data from bins centered within one-half of a bin-width
($2\bar{d}_{0}$) from the walls. Furthermore, to ensure that the DEM data
collapses used to determine the dimensionless material parameter $C_{\rm
seg}^{\rm S}$ (i.e., Figs. 5, 7, and 13) are representative of bulk behavior
and not wall effects, we use a more conservative criterion and do not include
DEM data from bins within $6\bar{d}_{0}$ of the walls.

The collapses of Figs. 5, 7, and 13 are obtained in the long-time regime, in
which the fields evolve slowly in time. Since the flow is quasi-steady, the
instantaneous concentration and velocity fields are simply arithmetically
averaged in time (using 152 instantaneous snapshots for vertical chute flow of
disks, 1000 snapshots for vertical chute flow of spheres, and 144 snapshots
for annular shear flow of disks) to obtain fields that only depend on the
spatial coordinate. Then, the necessary first and second-order spatial
derivatives of the field quantities (e.g., $\partial c^{\rm l}/\partial x$,
$\dot{\gamma}=\partial v_{z}/\partial x$, and $\partial\dot{\gamma}/\partial
x=\partial^{2}v_{z}/\partial x^{2}$ for vertical chute flow) are obtained from
these time-averaged fields. We apply a spatial derivative filter to the time-
averaged DEM fields in order to obtain accurate estimates of the spatial
derivatives. We have tested using both cutoff Gaussian functions and Lucy
functions (Weinhart et al., 2013; Tunuguntla et al., 2016) for the kernel
function of the derivative filter as well as a range of kernel function widths
to ensure that the reported results in this study are independent of these
choices.

Unlike for the steady concentration and velocity fields, which allow for
arithmetic time-averaging, the transient concentration fields for dense flows
of disks (e.g., Fig. 4(d)) are time-averaged in a slightly different manner
using a cutoff Gaussian filter. In particular, this process is performed by
applying a normalized, cutoff Gaussian time filter to the DEM data at each
$x$-position. Denoting the standard deviation of the Gaussian kernel function
as $\sigma_{\rm t}$, so that the cutoff time-width of the Gaussian kernel is
$6\sigma_{\rm t}$, the time-smoothed field quantity at a given $x$-position
and time $t$ is then given by the convolution of the DEM data over a time-
period of $6\sigma_{\rm t}$, centered at time $t$, with the cutoff Gaussian
kernel. We have tested a range of kernel widths $\sigma_{\rm t}$ to ensure
that the coarse-grained concentration fields appearing in this paper are
insensitive to this choice. For the transient concentration fields for dense
flows of spheres (second column of Fig. 10), no additional time-smoothing is
needed, and the concentration fields are simply instantaneous snapshots in
time. This is possible primarily because spatial smoothing is being done over
a greater volume and a larger number of grains, and therefore the
instantaneous concentration fields are relatively more smooth.

##  Appendix B Diffusion flux consistency test

Our process for determining the segregation flux–and hence the material
parameter $C^{\rm S}_{\rm seg}$–is based on the assumption that the
segregation and diffusion fluxes balance in the quasi-steady regime.
Therefore, it is essential that the dimensionless material parameter $C_{\rm
diff}$ appearing in the constitutive equation for the diffusion flux (2.16)
has been accurately determined, so that the coarse-grained diffusion flux is
accurate. In Section 3, we determined $C_{\rm diff}$ for dense flows of
frictional disks to be 0.20 using mean square displacement data from DEM
simulations of simple shear flow of a well-mixed bidisperse granular system.
In this Appendix, we perform an independent consistency check that tests
whether the constitutive equation for the diffusion flux (2.16) using this
fitted value of $C_{\rm diff}$ for disks is capable of predicting the
evolution of the $c^{\rm l}$ field in a diffusion-dominated problem.

Consider homogeneous simple shear flow of an initially-segregated system with
large grains (dark gray) on the bottom and small grains (light gray) on the
top, as shown in Fig. 14(a) for the case of $d^{\rm l}/d^{\rm s}=1.5$. The
rectangular domain has a length of $L=60\bar{d}_{0}$ in the $x$-direction and
a height of $H=120\bar{d}_{0}$ in the $z$-direction. As in Sections 2.3 and 3,
shearing along the $x$-direction and normal stress along the $z$-direction are
applied by the walls. We perform DEM simulations of simple shearing for a
nominal inertial number of $(v_{\rm w}/H)\sqrt{\bar{d}_{0}^{2}\rho_{\rm
s}/P_{\rm w}}=0.1$. We run the DEM simulation starting from the initially-
segregated configuration, and the spatiotemporal evolution of the coarse-
grained $c^{\rm l}$-field is shown in Fig. 14(b) for $d^{\rm l}/d^{\rm
s}=1.5$. We observe that the interface between large and small grains, which
is initially sharp, becomes diffuse with a transition width that grows with
time. We define a transition width at a given point in time as the distance
between the positions at which $c^{\rm l}$ equals 0.1 and 0.9 in snapshots of
the spatial $c^{\rm l}$ profile. This transition width as a function of the
square-root of the dimensionless time $\tilde{t}=t/(H/v_{\rm w})$ is plotted
in Fig. 14(c) as solid curves for grain-size-ratios of $d^{\rm l}/d^{\rm
s}=1.5$ and 3.0, displaying roughly linear behavior–typical of diffusive
behavior–with a slight dependence on $d^{\rm l}/d^{\rm s}$.

Figure 14: (a) Initially-segregated configuration for two-dimensional DEM
simulation of bidisperse simple shear flow with $d^{\rm l}/d^{\rm s}=1.5$ and
$8649$ flowing grains. Upper and lower layers of black grains denote rough
walls. Dark gray grains indicate large flowing grains, and light gray grains
indicate small flowing grains. A 10% polydispersity is utilized for each
species to prevent crystallization. (b) Spatiotemporal evolution of the large
grain concentration field, illustrating the transition width that grows with
time. (c) Normalized transition width versus square root of normalized time
$\tilde{t}=t/(H/v_{\rm w})$.

Next, we apply the continuum model for the evolution of $c^{\rm l}$ (2.18) to
this problem. As for planar shear flow of a well-mixed bidisperse system (Fig.
2), no pressure gradient is present. Therefore, the evolution of $c^{\rm l}$
is governed by (2.18):

$\dfrac{d{c}^{\rm l}}{dt}+\dfrac{\partial}{\partial z}\left(-C_{\rm
diff}\bar{d}^{2}\dot{\gamma}\dfrac{\partial c^{\rm l}}{\partial z}+C^{\rm
S}_{\rm seg}\bar{d}^{2}c^{\rm l}(1-c^{\rm
l})\dfrac{\partial\dot{\gamma}}{\partial z}\right)=0,$ (B.1)

where $\bar{d}=c^{\rm l}d^{\rm l}+(1-c^{\rm l})d^{\rm s}$. In this flow
configuration, the shear-strain-rate is approximately constant. Therefore, the
shear-strain-rate-gradient is approximately zero throughout, and the diffusion
flux is the dominant flux, which acts to remix the flowing grains. This may be
understood in the context of the local inertial rheology. Since the stress
ratio $\mu$ is spatially-constant in homogeneous planar shear, the resulting
inertial number field $I$ is also spatially constant. Since the inertial
number depends on both $\bar{d}$ and $\dot{\gamma}$, spatial variation in
$\bar{d}$ leads to spatial variation in $\dot{\gamma}$. This spatial variation
in $\dot{\gamma}$ is slight, and consequently, the magnitude of the diffusion
flux is at each point in space is much greater than the magnitude of the
segregation flux. While the effect of segregation is small, we still include
the shear-strain-rate-gradient-driven segregation flux with $C_{\rm seg}^{\rm
S}=0.23$ and solve (B.1) when analyzing the problem shown in Fig. 14(a). We
note that due to the small effect of segregation and the dependence of
$\bar{d}$ on $c^{\rm l}$, (B.1) is similar to but not exactly the same as the
linear diffusion equation in one dimension, so the solution is close to but
not exactly an error function.

We obtain predictions for the evolution of the $c^{\rm l}$-field by solving
(B.1) using the fully-segregated initial condition for $c^{\rm l}(z,t=0)$, no-
flux boundary conditions at $z=0$ and $z=H$, a spatially-constant value of
inertial number $I$ consistent with that prescribed in the DEM simulations, a
given grain-size-ratio $d^{\rm l}/d^{\rm s}$, and $C_{\rm diff}=0.20$. We note
that since $I$ is taken to be spatially-constant,
$\dot{\gamma}=(I/\bar{d})\sqrt{P_{\rm w}/\rho_{\rm s}}$ varies slightly in
space due to the $c^{\rm l}$-dependence of $\bar{d}$. We extract the
transition width as a function of time from continuum simulation results for
$d^{\rm l}/d^{\rm s}=1.5$ and 3.0 and include these results in Fig. 14(c) as
dashed lines. The continuum model predictions agree well with the DEM data and
are even capable of capturing the small difference due to the grain-size-
ratio. This result indicates that the expression for the diffusion flux (2.16)
and the fitted material parameter value $C_{\rm diff}=0.20$ are indeed
consistent with DEM data for disks.

## References

  * Artoni et al. (2021) Artoni, R., Larcher, M., Jenkins, J.T., Richard, P., 2021\.  Self-diffusion scalings in dense granular flows.  Soft Matter 17, 2596–2602.
  * Bancroft and Johnson (2021) Bancroft, R.S., Johnson, C.G., 2021\.  Drag, diffusion and segregation in inertial granular flows.  J. Fluid Mech. 924, A3.
  * Barker and Gray (2017) Barker, T., Gray, J., 2017.  Partial regularisation of the incompressible $\mu(I)$-rheology for granular flow.  J. Fluid Mech. 828, 5–32.
  * Barker et al. (2021) Barker, T., Rauter, M., Maguire, E., Johnson, C., Gray, J., 2021.  Coupling rheology and segregation in granular flows.  J. Fluid Mech. 909, A22.
  * da Cruz et al. (2005) da Cruz, F., Emam, S., Prochnow, M., Roux, J., Chevoir, F., 2005.  Rheophysics of dense granular materials: Discrete simulation of plane shear flows.  Phys. Rev. E. 72, 021309\.
  * Duan et al. (2021) Duan, Y., Umbanhowar, P.B., Ottino, J.M., Lueptow, R.M., 2021\.  Modelling segregation of bidisperse granular mixtures varying simultaneously in size and density for free surface flows.  J. Fluid Mech. 918, A20.
  * Dufty et al. (2002) Dufty, J.W., Brey, J.J., Lutsko, J., 2002.  Diffusion in a granular fluid. I. Theory.  Phys. Rev. E 65, 051303\.
  * Fan and 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 and Hill (2011a) Fan, Y., Hill, K.M., 2011a.  Phase transitions in shear-induced segregation of granular materials.  Phys. Rev. Lett. 106, 218301\.
  * Fan and Hill (2011b) Fan, Y., Hill, K.M., 2011b.  Theory for shear-induced segregation of dense granular mixtures.  New Journal of Physics 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.
  * Fan et al. (2015) Fan, Y., Umbanhowar, P.B., Ottino, J.M., Lueptow, R.M., 2015\.  Shear-rate-independent diffusion in granular flows.  Phys. Rev. Lett. 115, 088001\.
  * Fenistein and van Hecke (2003) Fenistein, D., van Hecke, M., 2003\.  Wide shear zones in granular bulk flow.  Nature 425, 256\.
  * Gray (2018) Gray, J.M.N.T., 2018.  Particle segregation in dense granular flows.  Annual Review of Fluid Mechanics 50, 407–433.
  * Gray and 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 and Thornton (2005) Gray, J.M.N.T., Thornton, A.R., 2005\.  A theory for particle size segregation in shallow granular free-surface flows.  Proc. Roy. Soc. London Ser. A 461, 1447–1473.
  * Henann and 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.
  * Henann and Kamrin (2014) Henann, D.L., Kamrin, K., 2014\.  Continuum thermomechanics of the nonlocal granular rheology.  Int. J. Plasticity 60, 145–162.
  * Henann and Kamrin (2016) Henann, D.L., Kamrin, K., 2016\.  A finite element implementation of the nonlocal granular rheology.  Int. J. Numer. Meth. Engng. 108, 273–302.
  * Hill and Fan (2008) Hill, K.M., Fan, Y., 2008.  Isolating segregation mechanisms in a split-bottom cell.  Phys. Rev. Lett. 101, 088001\.
  * Hill and 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.
  * Johnson et al. (2012) Johnson, C., Kokelaar, B., Iverson, R.M., Logan, M., LaHusen, R., Gray, J., 2012\.  Grain-size segregation and levee formation in geophysical mass flows.  Journal of Geophysical Research: Earth Surface 117, F01032.
  * Jop et al. (2005) Jop, P., Forterre, Y., Pouliquen, O., 2005.  Crucial role of side walls for granular surface flows: consequences for the rheology.  J. Fluid Mech. 541, 21–50.
  * Kamrin (2019) Kamrin, K., 2019.  Non-locality in granular flow: Phenomenology and modeling approaches.  Frontiers in Physics 7, 116\.
  * Kamrin and Henann (2015) Kamrin, K., Henann, D.L., 2015\.  Nonlocal modeling of granular flows down inclines.  Soft Matter 11, 179–185.
  * Kamrin and Koval (2012) Kamrin, K., Koval, G., 2012\.  Nonlocal constitutive relation for steady granular flow.  Phys. Rev. Lett. 108, 178301\.
  * Kamrin and Koval (2014) Kamrin, K., Koval, G., 2014\.  Effect of particle surface friction on nonlocal constitutive behavior of flowing granular media.  Comput. Part. Mech. 1, 169–176.
  * Kharel and Rognon (2017) Kharel, P., Rognon, P., 2017\.  Vortices enhance diffusion in dense granular flows.  Phys. Rev. Lett. 119, 178001\.
  * Kim and Kamrin (2020) Kim, S., Kamrin, K., 2020.  Power-law scaling in granular rheology across flow geometries.  Physical Review Letters 125, 088002\.
  * Komatsu et al. (2001) Komatsu, T.S., Inagaki, S., Nakagawa, N., Nasuno, S., 2001\.  Creep motion in a granular pile exhibiting steady surface flow.  Phys. Rev. Lett. 86, 1757\.
  * Koval et al. (2009) Koval, G., Roux, J.N., Corfdir, A., Chevoir, F., 2009\.  Annular shear of cohesionless granular materials: From the inertial to quasistatic regime.  Phys. Rev. E 79, 021306\.
  * Liu and Henann (2017) Liu, D., Henann, D.L., 2017\.  Non-local continuum modelling of steady, dense granular heap flows.  J. Fluid Mech. 831, 212–227.
  * Liu and Henann (2018) Liu, D., Henann, D.L., 2018\.  Size-dependence of the flow threshold in dense granular materials.  Soft Matter 14, 5294–5305.
  * MiDi (2004) MiDi, G.D.R., 2004.  On dense granular flows.  Euro. Phys. Journ. E. 14, 341–365.
  * 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\.
  * Rycroft et al. (2009) Rycroft, C.H., Kamrin, K., Bazant, M.Z., 2009.  Assessing continuum hypotheses in simulation of granular flow.  J. Mech. Phys. Solids 57, 828–839.
  * Savage (1998) Savage, S.B., 1998.  Analyses of slow high-concentration flows of granular materials.  J. Fluid Mech. 377, 1\.
  * Savage and Lun (1988) Savage, S.B., Lun, C.K.K., 1988\.  Particle size segregation in inclined chute flow of dry cohesionless granular solids.  J. Fluid Mech. 189, 311–335.
  * 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 Journal 61, 1524–1534.
  * Shinbrot and Muzzio (2000) Shinbrot, T., Muzzio, F.J., 2000\.  Nonequilibrium patterns in granular mixing and segregation.  Physics Today 53, 25–30.
  * 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.
  * Tang et al. (2018) Tang, Z., Brzinski, T.A., Shearer, M., Daniels, K.E., 2018\.  Nonlocal rheology of dense granular flow in annular shear experiments.  Soft Matter 14, 3040–3048.
  * Thornton et al. (2012) Thornton, A., Weinhart, T., Luding, S., Bokhove, O., 2012\.  Modelling of particle size segregation: calibration using the discrete particle method.  Intl. J. Mod. Phys. C 23, 1240014\.
  * Tripathi and Khakhar (2011) Tripathi, A., Khakhar, D.V., 2011\.  Rheology of binary granular mixtures in the dense flow regime.  Phys. Fluids 23, 113302\.
  * Tripathi and Khakhar (2013) Tripathi, A., Khakhar, D.V., 2013\.  Density difference-driven segregation in a dense granular flow.  J. Fluid Mech. 717, 643–669.
  * 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.
  * Tunuguntla et al. (2017) Tunuguntla, D.R., Weinhart, T., Thornton, A.R., 2017.  Comparing and contrasting size-based particle segregation models.  Comput. Part. Mech. 4, 387–405.
  * Umbanhowar et al. (2019) Umbanhowar, P.B., Lueptow, R.M., Ottino, J.M., 2019.  Modeling segregation in granular flows.  Annual Review of Chemical and Biomolecular Engineering 10, 129–153.
  * Utter and Behringer (2004) Utter, B., Behringer, R.P., 2004\.  Self-diffusion in dense granular shear flows.  Phys. Rev. E 69, 031308\.
  * van der Vaart et al. (2015) van der Vaart, K., Gajjar, P., Epely-Chauvin, G., Andreini, N., Gray, J.M.N.T., Ancey, C., 2015\.  Underlying asymmetry within particle size segregation.  Phys. Rev. Lett. 114, 238001\.
  * Weinhart et al. (2013) Weinhart, T., Hartkamp, R., Thornton, A.R., Luding, S., 2013\.  Coarse-grained local and objective continuum description of three-dimensional granular flows down an inclined surface.  Phys. Fluids 25, 070605\.
  * Zhang and Kamrin (2017) Zhang, Q., Kamrin, K., 2017\.  Microscopic description of the granular fluidity field in nonlocal flow modeling.  Phys. Rev. Lett. 118, 058001\.