Role of gravity or confining pressure and contact stiffness in granular rheology

The steady shear rheology of granular materials is investigated in slow quasi-static states and inertial flows. The effect of the gravity field and contact stiffness, which are conventionally trivialized is the focus of this study. Series of Discrete Element Method simulations are performed on a weakly frictional granular assembly in a split-bottom geometry considering various gravity fields and contact stiffnesses. While traditionally the inertial number, i.e., the ratio of stress to strain-rate timescales describes the flow rheology, we find that a second dimensionless number, the ratio of softness and stress timescales, must also be included to characterize the bulk flow behavior. For slow, quasi-static flows, the density increases while the macroscopic friction decreases with respective increase in particle softness and gravity. This trend is added to the $\mu(I)$ rheology and can be traced back to the anisotropy in the contact network, displaying a linear correlation between macroscopic friction and deviatoric fabric in the steady state. Interestingly, the linear relation holds when the external rotation rate is increased for a given gravity field and contact stiffness.


Introduction
Gravity is a critical factor in many natural (granular) phenomena like avalanches, landslides, sand-piles and even in some industrial applications [1,2]. Avalanches and debris flows play an important role in the transportation of mass existing at the surface of earth. Gravity-driven flows have also been observed on other planetary bodies of our Solar System and are of particular importance in understanding the geology of other planets and asteroids as well as for the human exploration of the Moon and Mars in the coming decades [3]. Currently, surface features found on Mars [4], Venus [5], and the Moon [6] are hypothesized to be the results of avalanches of granular material.
One of the important aspects of granular shear flows is the dependence of stress on external driving. Various experimental and numerical studies have shown that for slow-dense, quasi-static flows, the ratio of shear to compressive stress (effective friction coefficient) is independent of the imposed driving rate [7,8,9]. However, very little is known regarding the same in the presence of very weak gravity fields or low confining stress. Shear tests performed on parabolic flights have shown an increase in the friction coefficient at low confinement [10,11,12,13,14,15]. Brucks et al [16] also obtained similar trend using centrifuge experiments at gravity levels larger than Earth's gravity. Despite these studies, the effect of external compression (gravity) on granular flows is still poorly understood, which leads to issues like exploration vehicles getting stuck in the Martian soil.
Soft materials like hydrogel and elastomer, which can support large deformation are of increasing importance in engineering and biological applications such as tissue scaffolding, biosepration and micro-and-nano-printing [17]. While inertial number has been relatively successful in understanding the dynamics of rigid particles [18], elasticity becomes relevant for soft particles [8,19]. The deformability of the soft particles has been shown to affect the force network close to the jamming transition [20]. Recent study by Vaart et al [21] has shown different rheological behavior for hard and soft particle suspensions. Despite the increasing importance, the models for soft deformable particles have been largely ignored.
We claim here that these two factors, i.e., gravity and softness are two aspects of the same phenomenon. We aim to test this claim by answering the following questions: (1) How does gravity and softness of particles affect the bulk flow behavior? (2) Is there a unique law that can describe the flow behavior on Earth, Mars and the Moon for both soft and rigid particles?
In this paper, we address the above questions with a focus on dense, frictional, quasistatic granular flow. Using Discrete Element Method (DEM), we simulate cohesionless frictional granular material in a split-bottom ring shear cell. An important aspect of this setup is that the shear rate is given solely by external rotation rate and the geometry. At the same time, in this geometry the local strain rate does not depend strongly on the external compression [22], unlike the inclined plane and rotation drum where gravity has a strong effect on the local strain rate [16,23,24]. To study the effect of gravity and particle softness, we independently vary both parameters by two orders of magnitude. A change in particle softness provides an adjustment on the microscopic scale, while gravity is a macroscopic (field) modification. We find that they have similar effect at the mesoscopic (local) scale. The bulk behavior can be described well using a dimensionless parameter, defined as the ratio between the time scales due to gravitational compression and contact stiffness. Further, by increasing the external rotation rate, we study the dependence of macroscopic friction and contact network anisotropy (deviatoric fabric) on the inertial number. The dependence of macroscopic friction and deviatoric fabric on pressure is added to µ(I) rheology. Additionally, we find some non-local effects in our results due to the presence of gradients in both stress and strain rate, which are quantified by following an approach similar to Koval et al [25].
The outline of this paper is as follows: In section 2 we explain our numerical setup and methodology. We present our results for quasi-static state in section 3. In section 4, we provide results on the rheology and combine it with the results from quasi-static state to present new rheological laws. We then close the paper with a discussion and conclusion of our results in section 5 along with a possible outlook for future research.

Discrete Element Method (DEM)
We present our numerical simulation scheme and setup in sections 2.1 and 2.2 respectively. Section 2.3 briefly presents our averaging methodology and definitions of the tensorial quantities. We summarize various time scales associated with the system in section 2.4.

Model
Our computational techniques are based on the soft-sphere DEM simulations as developed by Cundall and Strack [26], Walton [27] and Luding [28,29,30]. The normal force between the particles in contact is modeled by a Hookean spring as f n = −k n δ n − η n v n , where k n , δ n , η n , and v n are the normal stiffness, particle overlap, normal viscous damping coefficient, and relative velocity in the normal direction, respectively. Similarly, the tangential force is modeled as f t = −k t δ t − η t v n , where k t = 2k n /7, δ t , η t = η n /4, and v t are the tangential stiffness, relative displacement in tangential direction, tangential viscous damping coefficient, and relative velocity in tangential direction, respectively. We also introduce Coulomb's friction between the particles, where the tangential force f t is switched to the sliding force f s = −µ p |f n |, µ p being the particle friction coefficient when f t exceeds the critical value, i.e. |f t | > µ p |f n | (with µ p = 0.01) [30].
To study the effect of particle softness on macroscopic behaviors, we explore a range of normal contact stiffness k n , from 10 Nm −1 ≤ k n ≤ 10 4 Nm −1 . While changing k n , the time step for numerical integration δt is adjusted to be 1 50 times the contact duration to ensure accurate dynamic integration [30]. When k n is changed, k t and η n are changed as well, to ensure that the coefficient of restitution remains unchanged. Figure 1 is a sketch of our numerical setup [31,32,33,34]. In this figure, the inner, split, and outer radii are given by R i , R s , and R o , respectively, where the concentric cylinders rotate relative to each other at a rate Ω around the symmetry axis (the dotdashed line). The ring shaped split at the bottom separates the moving and static parts of the system, where a part of the bottom and the outer cylinder rotate at the same rate. The system is filled with N ≈ 3.7 × 10 4 polydispersed spherical particles with density ρ = 2000 kgm −3 = 2 gcm −3 up to height H. The average size of particles is a = 1.1 mm, and the width of the homogeneous size-distribution (with a min /a max = 1/2) is 1 − A = 1 − a 2 / a 2 = 0.18922. The cylindrical walls and the bottom are roughened due to some (about 3% of the total number) attached/glued Figure 1. (Color online) A sketch of our numerical setup consisting of a fixed inner part (light blue shade) and a rotating outer part (white). The white part of the base and the outer cylinder rotate with the same angular velocity, Ω, around the symmetry axis (dot-dashed line). The inner, split, and outer radii are given by R i = 8 a , R s = 40 a , and R o = 52 a , respectively, where each radius is measured from the symmetry axis. The gravity, g, points downwards as shown by an arrow.

Split-bottom ring shear cell
particles [33,35,36]. The initial state of the system is prepared frictionless. Friction is switched on at the onset of shear.
When there is relative motion at the split, a shear band initiates at the bottom from the split position R s and propagates upwards and inwards, remaining far away from the cylindrical walls in most cases [32,34]. The qualitative behavior is governed by the ratio H/R s and three regimes can be observed as reported in [32,34]. We keep H/R s < 0.5, such that the shear band reaches the free surface and stays away from inner wall.
To understand the effect of gravity on the bulk behavior, gravity is varied in the range 0.5 ms −2 ≤ g ≤ 50 ms −2 . The details regarding rotation rate of the system are reported in table 1. The total simulation time is chosen such that the cylinder completes half a rotation in that time.

Local averaging
One of the goals of current research in the granular community is to derive macroscopic continuum theory based on the given micro-mechanical properties [37,38,39]. We assume translational invariance in the azimuthal θ−direction. The averaging is thus performed over toroidal volumes over many snapshots in time, leading to fields Q(r, z) as function of the radial and vertical positions. Here, the averaging is performed with spacings ∆r and ∆z around two particles diameter in radial and vertical directions (averaging procedure for two and three dimensions is discussed in detail in Refs. [35,38] and hence not discussed further here).

Macroscopic (tensorial) quantities
Here, we present general definitions of the averaged microscopic and macroscopic quantities -including strain rate, stress and fabric (structure) tensors.
The strain rate is calculated by averaging velocity gradients of particles over cells and is given bẏ where Greek letters represent components radial distance r, azimuthal angle θ, and height z, while ▽ represents the gradient in cylindrical coordinate [35]. The stress tensor is given by with particles p, contacts c, mass m p , velocity v p , force f c and branch vector l c , while Greek letters represent components r, θ, and z [35]. The first term is the sum of kinetic energy fluctuations, and the second term involves the dyadic product of the contact-force with the contact-branch vector. The quantity which describes the local configuration of a granular assembly is the fabric tensor [40,39] and is given by where V p represents the particle volume which lies inside the averaging volume V , n c is the normal unit branch-vector pointing from the center of particle p to contact c.

Isotropic and deviatoric parts
Any given tensor Q can be uniquely decomposed into isotropic and deviatoric parts as with Q V = 1 3 Tr Q and the traceless deviator Q D . The latter contains information about the eigensystem of Q, that is identical to the eigensystem of Q D itself.
Let us assume Q 1 , Q 2 , and Q 3 are the eigenvalues of Q D sorted such that Q 1 ≥ Q 2 ≥ Q 3 . Based on the normalization, we use following definition to quantify the anisotropy of any tensor Q D using a single scalar quantity: the deviatorsǫ dev , σ dev , and F dev refer to strain rateǫ αβ , stress σ αβ and fabric F αβ respectively. Other definitions of Q dev are reported in [41]. γ =ǫ dev quantifies the local strain rate magnitude, which is very close to the form defined in cylindrical coordinates [42] as tested in [35]. The pressure p is the isotropic hydrostatic stress, while τ = σ dev quantifies objectively the shear stress. The deviatoric stress ratio, µ = τ /p, quantifies the "stress anisotropy"or the macroscopic friction coefficient. The volumetric fabric 3F v represents the contact number density, while the deviatoric fabric F dev quantifies the anisotropy of the contact network (as studied in detail in [43]). Table 1. Table showing the values of g (units of ms −2 ) and particle stiffness k n (units of Nm −1 ) used in our simulations, and various time scales associated with the system, as discussed in the main text (in units of s). The values of Tγ and T p are the average values reported at z = 2 d , H/2, H − 2 d in the center of the shear band.

Time scales
We characterize the dynamics of the system looking at different time scales. At first, we define two microscopic time scales related to the contact duration and the viscous damping coefficient between two particles in contact, respectively, as where m is the mass of a particle with mean diameter. T k and T η can be related to . Next, two time scales associated with external forces, i.e.
the gravity and external rotation rate, can be introduced as respectively, where T g is the time taken by a particle with zero initial velocity to fall a distance d /2. g Ω 2π The time scales, (6) and (7), are functions of material constants and applied external forces, hence, are constants throughout the system. In this sense, the time scales, T c , T η , T g , and T Ω , are global. On the other hand, two local time scales related to the local shear rateγ and pressure p, as introduced in [18,44] are: As shown in the following sections, the spatial distributions of pressure and shear rate are inhomogeneous due to gravity and shear band localization. Therefore, in contrast to the global time scales, T˙γ and T p are local field variables that depend on spatial position. The time scales can be combined to formulate dimensionless numbers that give indications of dominance of one of the time scales. For example, the inertial number, as introduced by [44,45,46], provides an estimate of the local rapidity of the flow. For I ≪ 1, the flow is quasi-static, where particles interact via enduring contacts and inertial effects are negligible. For I ∼ 1, the flow is in the dense inertial regime, and for I ≫ 1, the flow is in the rapid, collisional gas like state. A variant of inertial number I k =γ d / p kin /ρ can be defined by using only the kinetic pressure instead of the total pressure. Other expressions such as Savage or Coulomb numberγ 2 d 2 ρ/p is simply the square of the inertial number I [47]. A dimensionless parameter global compressibility can be introduced as the ratio between T g and T k : providing a global measure of compressibility of the bulk material. A high κ g signifies that the bulk material is compressible, which comes from either very high confinement by strong gravity or from low contact stiffness at particle level. On the other hand, when κ g is small, the average overlap is very small compared to the particle diameter, which means that the bulk material is closer to being the rigid limit. A similar dimensionless parameter can be introduced as: which estimates the compressibility of the material at the local scale. Table 1 shows typical values of timescales in our simulations, and table 2 reports various dimensionless numbers. We observe that for flows with a rotation rate Ω/2π = 0.01 s −1 and the gravity g ≥ 1 ms −2 , the inertial number I is well below 1 for all values of k n . For lower values of gravity g = 0.5ms −2 , the rotation rate is chosen to be Ω/2π = 0.005 s −1 , such that I stays in the same range. From this table we infer that the systems for wide a range of g and k n can be safely assumed to be in the rate-independent quasi-static state. However, we observe that with increase in rate of rotation Ω/2π, I begins to increase and the flow behavior enters into the rate dependent inertial regime.

Quasistatic rheology
In this section, we present our results on the analysis of macroscopic rheology in the quasi-static state. At first, in section 3.1 we study the steady state and critical state and the amount of time system requires to reach the critical state. We explore the effect of gravity and stiffness on the macroscopic friction coefficient in section 3.2. We also show the results of local volume fraction in section 3.3, and connect the rheology to the microscopic structure tensor in section 3.4. We will extend our analysis to dense inertial flows in section 4.

Steady state and critical state
Shearing leads to dilation and a build-up of shear stress and anisotropy in fabric. Additionally, anisotropy in the fabric of the granular medium needs a finite amount of strain to build up. Since, we are interested in the steady state flow properties, we need to ask: what is the minimum time or equivalently strain required to reach the steady state flow regime?
We define ϕ = 180 2π Ω∆t (where ∆t is the simulation time) as a shear parameter, which is the same as the amount, in degrees, through which the system rotates. In general, the time required for stabilization can depend on the considered variable [25]. In the following, we consider the relaxation of various quantities to the steady state (both locally and globally averaged). As we expect this relaxation to be slowest for small Ω, we perform these tests for Ω/2π = 0.01 s −1 that is the smallest Ω we explore in our simulations.
At first, we analyze the global quantities (averaged over the whole system) like the averaged kinetic energy and average number of contacts. We find that they relax very fast (ϕ ∼ 5) and are not reliable since they only represent the fastest relaxation in the system.
Next, we analyze the relaxation of local quantities. In order to estimate the relaxation time for local quantities, we analyze local velocity profiles. We find that such a relaxation requires (ϕ ∼ 30 or more) which is in agreement with Ries et al [22]. Consequently we consider that the condition to approach the steady state is ϕ ≥ 30. Experimental results from Wortel et al [48] also show that the system needs to be rotated by an angle of approximately 30 degrees in order to reach a reasonable steady state. When the steady state is reached, the local averaging is performed over almost 600 snapshots distributed over ϕ ≥ 30.
The consistency of the local averaged quantities also depends on the accumulated local shear strain during the procurement of data. We concentrate our interest in the region where the system can be considered in a critical state, which occurs at large enough shear strain γ =γ∆t ≥ 1. In other words, the system can be assumed to be in the critical state. The critical state is a unique state reached after long shear, where material deforms with applied strain without any change in normal stress, shear stress and volume fraction, and the system forgets its preparation history [49]. Figure 2 shows the local shear stress τ (r, h) plotted against the local pressure p(r, h). We observe that for a given pressure, τ is higher for largerγ, however forγ >γ c (witḣ γ c ≈ 0.08 s −1 ), τ becomes almost independent of the local strain rate. This means that τ /p is almost constant for all data points with strain rateγ >γ c . In other words, if the dimensionless shear length l γ = t avγ [35] exceeds unity (where t av is the time over which averaging of the data is performed), layers can be assumed to be sheared by almost a particle diameter. For l γ ≥ 1, the shear deformation can be considered to be fully established, and the system is assumed to be in the critical state [35].
In the same setup, Ries et al [22] and Szabó et al [50] also found that after long enough shear, the material inside the shear band reaches the critical state and can be characterized by estimating local accumulated strain γ ≥ 1. Our previous works [35,51] also showed that for rotation rate Ω/2π = 0.01 s −1 ,γ c ≈ 0.1 s −1 is the shear rate above which the shear band is well established. Since we are interested in the flow behavior of the material, as default in the rest of the paper we focus only on the data well inside the shear band with local strain rate, unless specified otherwise.

Friction coefficient
We will now turn our attention towards the effect of gravity and softness on the macroscopic friction coefficient τ /p in the quasi-static state. In previous studies, the friction coefficient has been assumed to be independent of both the particle stiffness and gravity. However, particles used in these were extremely rigid. Few studies were performed systematically investigating the dependence of the flow behavior on gravity [10,11,12,13,14,15]. Figures 3(a) and 3(b) show shear stress-pressure curves for different values of normal stiffness k n and external gravity g, respectively. In these plots, for the sake of clarity, both shear stress and pressure are plotted only at the center position R c of the shear band (R c being the the position at which the local shear rate is maximum). For a better comparison, both shear stress and pressure are normalized with the maximum pressure p max reached in the simulation with particular k n and g (so that both abscissa and ordinates are of the same order). We observe that both softness of the particles (interpreted as opposite of contact stiffness) and external gravity drastically affect the shear stress. Moreover, they act in the same direction as for a given pressure τ decreases with increase in either particle softness or the external gravity. We also see that with increasing softness or gravity, the relation between τ and p becomes non-linear.

Linear approximation
To understand the dependence of the macroscopic friction coefficient in a quasi-static state on the softness and gravity, we estimate it as the slope where µ global 0 is the global friction coefficient which depends neither on the shear rate nor on the pressure. Figure 4(a) displays the global friction coefficient µ global 0 plotted against gravity g for different values of the normal stiffness k n , as given in the legend. We observe that µ global 0 decreases with increasing gravity, while it increases with increasing k n . Figure  4(b) shows the global friction coefficient with different values of the normal stiffness k n and gravity g, where all results of µ global 0 are collapsed if we plot them against κ = κ 2 g (κ g is given by (10)).
In figure 4(b), the solid line is given by where µ global r = 0.17 ± 0.01 is the global friction coefficient in the rigid particle limit and the exponent and characteristic global compressibility are given by α ≃ 0.5 and κ 0 ≃ 2.01, respectively.
Previous microgravity parabolic flight and centrifuge experiments [10,12,13,15,16,14] showed a similar decreasing trend of the macroscopic friction coefficient on gravity. Few authors [12,15] attributed this dependence to the fact that at low gravity, the body forces become weak and the electrostatic cohesive forces begin to dominate. Klein et al [12] also argued that a load-dependent interparticle friction coefficient might be responsible for this behavior. However, no cohesive force or load-dependent friction was implemented in any of the DEM simulation data presented here. Hence, we claim that  there should be an additional mechanism responsible for this interesting behavior. In order to gain a better understanding in the following, we have a closer look by studying the system locally.

3.2.2.
Local compressibility Since our system is heterogeneous in both stress and strain rate fields, a local description of the system is highly desirable. In the shear stresspressure (τ − p) curves for different softness and gravity (figure 3), the dependence of shear stress on pressure slightly "bends" with increasing softness and gravity, i.e., the friction coefficient has a dependence on the pressure and the shear stress becomes a nonlinear function of pressure as: where µ local 0 (p, k n ) is a local friction coefficient which depends on pressure and contact stiffness. Figure 5 shows the local friction coefficient with different values of the normal stiffness and gravity, where all results of µ local 0 (p, k n ) collapse if we introduce the local compressibility, defined as square of the ratio between two time scales, T k and T p . p * can also be interpreted as non-dimensional average overlap (scaled with mean particle diameter).
In this figure, we scanned through a wide range of p * by systematically varying g and k n . We observe that for p * < 5 × 10 −4 , µ local 0 (p * ) is almost constant, while for higher values, µ local 0 (p * ) decreases with p * up to p * ≈ 0.1. This dependence can be written in the form, where µ local r = 0.172 ± 0.01 is the value of macroscopic friction in the rigid limit, which is in fair agreement with contact dynamics shear simulations for the same particle friction coefficient [52]. The exponent is found to be β 1 ≈ 0.5 ± 0.04 and the characteristic local compressibility is p * σ = 10.1 ± 0.2. As one extreme of p * , for p * = 0.1 the average overlap is almost 10% relative to the mean particle diameter, i.e., the soft particle limit. The upper bound of µ local 0 (p * ) is the low compression case, i.e., the rigid limit, where the average overlap is much smaller (1/10000) compared to the particle diameter and µ local 0 (p * ) is almost double as large as for p * ≈ 0.1.

Local volume fraction
In figure 6, the local volume fraction ν(r, h) is plotted against the local compressibility, p * . Because of slow quasi-static flows, no strong dilation is observed, i.e., no strong dependence of ν on local shear rate is observed. We observe that the packing is rather loose for lower p * and tends to a critical value ν c = 0.642 ± 0.002, as observed in [53]. The data can be fitted well by the functional form with a * = 0.48 ± 0.02 (a * can be further expressed in terms of volumetric fabric as reported in [40,54]). Interestingly, no significant difference in volume fraction ν is observed for p * < 10 −3 , while for p * > 10 −3 within the fluctuations, ν increases almost linearly with p * (the curvature is due to the logarithmic p * axis). The relation between ν and p * is well established in the case of static packings [40,54,55]. Here we show that the same relation holds for a slow granular flow, involving considerable finite but small strain rates.

Local structure
Shearing of a granular assembly always leads to the buildup of a contact anisotropy in the system [56,57,58]. To study contact anisotropy in our system, we analyze the deviatoric fabric as defined in (3). We use the second invariant or F dev (5) of the deviatoric fabric tensor to quantify anisotropy of the contact network in the system. Figure 7 displays the local deviatoric fabric, F dev (r, h), plotted against the local compressibility p * , where F dev (r, h) for different values of the particle stiffness and gravity is found to collapse on a unique curve (solid line). This dependence can be written in a similar fashion to (17),

Anisotropy
where F r dev is the anisotropy of contact network in the rigid limit, the exponent is found to be β 2 ≈ 0.5 ± 0.03, and p * F ≈ 26.3 ± 0.6. The decrease in F dev with increasing p * can be explained in terms of the increasing volume fraction ν(r, h) with increase in p * . In the case of a denser packing, particles have less free space to re-arrange (and thus align along preferential direction) and build up anisotropy in response to the local shear, compared to a relatively loose packing (at low p * ).
In figures 5 and 7, we observe that the the local shear resistance and the local anisotropy of the contact network in the quasi-static state show a similar trend. In figure 8, we plot µ local 0 (p * ) against F dev (p * ) for different values of κ, where a linear correlation can be inferred as, where µ iso = 0.01 ± 0.01(≈ 0) is the friction coefficient in the (extrapolated) limit of the isotropic contact network (F dev = 0) and b = 1.38 ± 0.02 is a constant of proportionality. This clearly shows that the shear resistance seems to accompany the anisotropy in the contact network in the critical state. It is worthwhile to mention that no information about such a relation in the transient regime can be inferred from our analysis, while it is strongly supported by our data in the critical state. It also links the increase in the friction coefficient with decreasing gravity (as observed in figure 4(a)) to the change in the contact network configuration of the material at lower g.

Shape factor
In this section, we compare the shape factor ( Q 2 Q 1 ) for stress and fabric tensors, where Q 2 , and Q 1 are the eigenvalues of the deviatoric tensors as defined in section 2.3.2. The shape factor quantifies the out of shear plane neutral eigenvalue, compared to the maximum eigenvalue. Note that for strain rate tensorǫ 2 /ǫ 1 ≡ 0 due to geometry and symmetry, i.e., we have plain-strain shear.
In figure 9(a), we plot the shape factor for the stress tensor. Different symbols represent different values of compressibility κ as given in the legend of figure 6, while black circles show the data in the center of shear bands for these simulations. We observe that the shape factor is highly fluctuating for the two extremes of p * , while in the range 10 −4 < p * < 10 −2 , it is significantly below zero. This implies that the stress in the shear plane is higher as compared to the axial stress along the neutral direction orthogonal to the shear plane. The sign means that this axial stress is reduced perpendicular −1 and enhanced +1/2 within the shear plane [39]. The dashed line is the data from [39], where authors studied the flow behavior on an inclined plane. We observe that the sign for both the cases is negative, while the magnitude is different, which might be due to the difference in interparticle friction. In figure 9(b), the shape factor of the fabric tensor fluctuates (strongly) around zero. It is important to mention that the fluctuations in the data are from a single simulation.
These two observations suggest that the fabric and stress tensors behave differently even though they are proportional in magnitude (norm) as shown in figure 8. The fabric tensor is in a planar state like the strain rate tensor, whereas the induced stress state is triaxial, as expected for a solid-like material [39]. F 2 /F 1 tends to positive values for larger p * , further establishing the difference between structure and stress tensors. However, in order to have a clear picture for the fabric tensor, the strong and weak subnetworks should be studied separately, since only the strong subnetwork carries almost all of the fabric anisotropy [43,59].
As discussed in section 3.1 the cutoff shear rateγ c can depend on the simulation time or the averaging time. In this section, we focused on the data only inside the shear band, which are in the critical state and have forgotten their initial configuration due to large strain. However, the velocity gradients in the setup are smooth, which implies that part of the system outside the shear band is also flowing, but only very slowly. If the simulation runs longer, the cutoff can be lowered, eventually if the simulation would run infinitely long, no cutoff on the local strain rate is needed. If we lower the cutoff on local strain rate (see next section), by settingγ c (Ω) ≡ Ω 2π , we observe much wider variation of data in our results, especially for the local friction coefficient, the deviatoric fabric and the volume fraction. However, the shape factors are not strongly influenced by a change inγ c , by nearly an order of magnitude.

Complete rheology
The previous section showed that in the quasi-static state the friction coefficient and deviatoric fabric are strongly correlated in the critical state. Motivated by this, we check if this correlation also works in the rate-dependent inertial regime. To test the correlation, high inertial number data are generated by varying the external rotation rate Ω for a fixed gravity and contact stiffness. In the following, we will explore the evolution of the local macroscopic friction coefficient µ and deviatoric fabric with I and combine it with the dependence of both µ and deviatoric fabric on p * , to propose the complete rheological law. It is important to mention that in this section the cutoff on local strain rate is lowered toγ c (Ω) ≡ Ω 2π , so as to capture the maximum data and present a unique rheology law outside and inside the shear band.

Friction law
The local friction coefficient µ is plotted against inertial number I in figure 10. Different symbols show data from different rates of rotation as given in the inset; the black solid circles represent the data in the center of the shear band. The solid black line shows the friction law, as proposed in [45]: with µ local 0 (p * ), as given in (17). We observe that data from simulation using a single gravity and contact stiffness does not give a wide variation in µ and µ local 0 = 0.14, , which is similar to the frictional law proposed in [18,44]. Two different trends emerge, i.e., the shear band center data can be very well fitted by (21) and for I ≥ 0.01 data collapse on a unique curve. On the other hand, for lower values of I, deviations from this relation are observed, depending on the external rotation rate. The friction coefficient in slow flows (steady state) becomes smaller than µ local 0 , i.e. in our system the granular material is able to flow below µ local 0 . The deviation of our data from the main law (21) is consistent with observations in [25,60], where this deviation was explained based on the heterogeneity in the stress  field (due to strain rate). In our system, we have gradients in stress arising due to gradients in both strain rate and pressure (due to gravity).
In order to quantify the deviation from (21), we fit the data with: where α is a constant and I * is the characteristic inertial number when µ ∼ = µ local 0 . This relation is similar as proposed in [25] for two-dimensional ring shear cell data. As the above relation was derived for a two-dimensional setup with constant pressure, we fit it to our data at three different heights (i.e. pressure levels), close to top, at mid-height and close to bottom. In figure 10, different colored dashed lines represent this fit at the mid-height of the system. We observe that the prediction is in close agreement with the data, even though our system has different dimension and boundary conditions. Appendix A shows the data and corresponding fits for different heights (pressures). We find that both α and I * do not show a clear dependence on pressure.

Fabric anisotropy
In order to look for the connection between anisotropic fabric and macroscopic friction coefficient in the inertial regime, here we explore the dependence of F dev on I. In figure  11, we plot the local F dev as obtained by simulations with different rates of rotation against I. We observe that like µ, F dev varies strongly against I and its dependence on I can be represented as:  Figure 11. (Color online) The local fabric anisotropy F dev plotted against the inertial number I for results from simulations with different rates of rotation. Different symbols represent different rates of rotation as given in the legend. Black circles represent the data in the center of the shear band. Green, red, and black lines are fit to (23) for pressure levels p = 100, 200, and 400 Nm −2 respectively, with fit parameters given in Table 3. The arrow shows increasing pressure.
with F 0 dev being the fabric anisotropy in the quasi-static state, F dev is the limiting fabric anisotropy, and I F 0 is the typical inertial number, which is an order of magnitude different from I σ 0 . Green, red and black lines show the fit to above relation at pressure levels 100, 200 and 400 Nm −2 , respectively, with points inside the shear band highlighted (black circles). Fit parameters to these results are presented in table 3. It is noticeable that unlike µ, I alone is not able to describe F dev , with the effect of pressure being prominent in case of slow flows i.e., low I. In contrast, for fast flows, the anisotropy seems to become independent of pressure.
The increase in the contact anisotropy with inertial number is in accordance with the recent studies [39,58]. It is important to mention that for even higher rates of rotation of the system, i.e. inertial number I > 0.1, F dev shows a different behavior as predicted by (23) and a decreasing trend is observed (as reported in [61]), which is beyond the scope of this work. This might be due to the fact that for I > 0.2 the packing becomes very loose (ν ≤ 0.55). Also for such high rates of rotation, the centrifugal force on grains due to rotation becomes comparable to the gravitational force. As a result, the top surface is not flat anymore, instead the surface develops a dip in the middle, as observed previously [61,62,63]. In this range, also the kinetic and contact contributions of the local macroscopic friction µ become comparable.
Starting from both variations of macroscopic friction and fabric anisotropy as a function of inertial number I, it is tempting to ask the question if the correlation in (20) holds for the inertial regime as well. The result is displayed in figure 12. Solid line shows (20), which fits well the shear band center data being highlighted by black circles. Here again, we find data outside the shear band showing a different trend. Red dashed line separates the data into quasi-static and dense inertial regimes. It is interesting to note   that (20) (as shown by black solid line) very well captures the trend at the onset of dense inertial regime. However, for faster flows, where rate dependence becomes prominent, a different trend is observed which can also be fitted by a slightly different linear relation (dot-dashed line). This shows that the macroscopic friction and fabric anisotropy are correlated even in the dense inertial regime.

Discussion and Conclusion
To summarize, this work is an exploration of the flow behavior for 3D granular shear flow using DEM simulations. We particularly focused on the effect of external compression (gravity) and the particle softness on the flow behavior, considering both stress and structure.
Quasi-static flows Our study shows that the shear strength (macroscopic friction coefficient µ) of the material decreases with increase in either gravity or particle softness for the quasi-static flows. We find that the data for different gravity and particle softness can be expressed as a unique power law, when analyzed in terms of a control parameter called local compressibility p * (as defined in (16)). p * can be interpreted as the non-dimensional local average overlap (scaled with mean particle diameter) and can be used to quantify the softness of the bulk material relative to the local compression level. Low values of p * signify rigid particles, while a high p * implies soft, deformable particles. Both softness and gravity are also found to affect the local microstructure, i.e., the anisotropy of the contact network, which is quantified by the deviatoric fabric (F dev ). We show that the deviatoric fabric can also be expressed as a power law of p * with an exponent similar to that of the shear strength. This points out that the local anisotropy of the contact network (deviatoric fabric) and the shear strength of the material are highly correlated in the slow, quasi-static flows and the shear strength follows the anisotropy of the contact network. These results are highly significant for planetary studies regarding the shear strength of the granular material in extraterrestrial bodies such as Moon or Mars. Significant amount of experimental work using parabolic flight have shown the increase in shear strength of the material with decreasing gravity [10,11,12,13,14,15] without proper explanation. We propose that the anisotropy, i.e., the rearrangement in the contact network, is the key parameter for this dependence if no new forces are involved. With decreasing gravity, the packing becomes loose (due to decrease in body force acting on the particles), which in turn provides more free space to the particles to rearrange (and thus align) in response to the local shear. The fact that the particle softness and gravity have been shown to have similar effects on the local flow behavior, makes this work equally relevant for soft particles, that find their applications in many engineering and biological systems [17].
Since it is extremely difficult and expensive to perform in situ experiments on the Moon (or the parabolic flight), the 'compensation' effect we find with the ratio of gravity g and particle stiffness k n allows to mimic variable gravity by tweaking/tuning the particle stiffness.
Inertial flows Further, to study the rheology of the system for gravity g = 10 ms −2 , the rotation rate of the system is increased. For faster flows the system enters into a rate dependent inertial regime, consistent with previous studies [18,44,46]. We find that both the macroscopic friction µ and contact anisotropy show a similar increasing behavior with inertial number I. This shows that macroscopic friction and contact anisotropy are correlated also in the dense inertial regime. The increase of µ with I as observed in the inertial regime is accompanied by the evolution of the microstructure (contact anisotropy) with increasing inertial number. This picture is consistent with the recent study of Azéma et al [58].
Open issues We find that the frictional laws obtained from homogeneous shear flows [44] can be applied locally in the inertial regime, while they fail to predict the behavior of the material in the slower, quasi-static regime. The local rheology laws can be applied to our data in the center of the shear bands, where the strain-rate and stress gradients are zero, hence the material can be assumed to be homogeneously sheared. However, away from the center of the shear bands, in the quasi-static regime, we observe a nearly identical range of µ values corresponding to a completely different range of I. We find that in our system the material is able to flow even below µ local 0 , but only very slowly. We have shown that some observations can be explained by using an approach similar to Koval et al [25]. These deviations might be well captured using the non-local models by Kamrin and coworkers [60,64,65]; this work is in progress. Another related issue which remains untouched is the effect of particle softness and external compression (gravity here) on the non-locality. A study of effect of gravity on primary and secondary velocity fields, as done recently in [66,67] also deserve a further study.
Conclusion The macroscopic friction (shear strength) of the material is found to be affected not only by the local shear rate, but also by external compression (gravity) and softness of the particles. While traditionally the inertial number, the ratio of stress to strain-rate time-scales, is dominating the flow rheology, we find that a second dimensionless number, the ratio of softness and stress time scales, must be involved to characterize the bulk flow behavior. For very slow shear rate the former can be ignored, while the latter two affect the shear strength by decreasing it with increase in either gravity (and thus local pressure) or particle softness. For faster flows, the macroscopic friction is found to increase in general with increasing shear rate. However, the tails of shear bands feature an anomalously small macroscopic friction -as observed previously [25,60,65]. For the dependence of macroscopic friction on above three quantities, the change in local microstructure (contact anisotropy) is found to be a key parameter, that was not often considered yet.
Looking towards the future, we are now in a position to address various important issues, such as unexpectedly high shear strength of the material at low (confining) stress or reduced gravity and a direct relation between the contact anisotropy and the shear strength of the material. These issues are vital for a better explanation of the macroscopic behavior of the granular systems from a microscopic observation. The current study dealt with a dense system with small interparticle friction (µ p = 0.01), where the effect of softness on the macroscopic behavior is more direct than for large µ p . However, an issue which remains unanswered and will be an extension of this study is whether the same effect can also be observed for relatively loose system (with higher interparticle friction). Further, the question whether the correlation between contact anisotropy and shear strength is just a consequence of relatively low interparticle friction or if it will also hold for a more realistic material (with higher interparticle friction) remains to be answered.  "Jamming and Rheology" project of the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is financially supported by the "Nederlandse Organisatie voor Wetenschappelijk Onderzoek" (NWO), is acknowledged. (c) Figure B1. (Color online) µ plotted against F dev for different local pressures in the system (a) p = 100, (b) p = 200, and (c) p = 400 Nm −2 . The solid line represent the corresponding fit to (20), while the dashed line is the best fit to the data.