A general(ized) local rheology for wet granular materials

We study the rheology of dry and wet granular materials in the steady quasistatic regime using the discrete element method in a split-bottom ring shear cell with focus on the macroscopic friction. The aim of our study is to understand the local rheology of bulk flow at various positions in the shear band, where the system is in critical state. We develop a general(ized) rheology, in which the macroscopic friction is factorized into a product of four functions, on top of the classical μ ( I ) rheology, each of which depends on exactly one dimensionless control parameter, quantifying the relative importance of different micro-mechanical machanisms. These four control parameters relate the time scales of shear rate t γ ˙ , particle stiffness tk, gravity tg and cohesion tc, respectively, with the governing time scale of confining pressure tp. While t γ ˙ is large and thus of little importance for most of the slow flow data studied, it increases the friction in critical state, where the shear rate is high and decreases friction by relaxation (creep) where the shear rate is low. tg and tk are comparable to tp in the bulk, but become more or less dominant relative to tp at the extremes of low pressure at the free surface and high pressure deep inside the bulk, respectively. The effect of wet cohesion on the flow rheology is quantified by tc decreasing with increasing cohesion. Furthermore, the proposed rheological model predicts well the shear thinning behavior both in the bulk and near the free surface; shear thinning rate becomes less near the free surface with increasing cohesion.


Introduction
The ability to predict a material's flow behavior, its rheology (like the viscosity for fluids) gives manufacturers an important product quantity. Knowledge on material's rheological characteristics is important in predicting the pourability, density and ease with which it may be handled, processed or used. The interrelation between rheology and other product dimensions often makes the measurement of viscosity the most sensitive or convenient way of detecting changes in flow properties. A frequent reason for the measurement of rheological properties can be found in the area of quality control, where raw materials must be consistent from batch to batch. For this purpose, flow behavior is an indirect measure of product consistency and quality.
Most studies on cohesive materials in granular physics focus on dry granular materials or powders and their flow [15,39]. However, wet granular materials are ubiquitous in geology and many real-world applications where interstitial liquid is present between the grains. Many studies have applied the ( ) m I -rheology to flows of dry materials at varying inertial numbers I [40,41,43,45,49]. Studies of wet granular rheology include flow of dense non-Brownian suspensions [3,13,14,21]. Here, we study partially wetted system of granular materials, in particular the pendular regime, which is also covered in many studies [35,38,51]. While ideally, unsaturated granular media under shear show redistribution of liquid content among the contacts [28,36], we assume a simplistic approach of homogeneous liquid content for liquid bridges of all contacts. One of the important aspects of partially wetted granular shear flows is the dependence of shear stress on the cohesive forces for wet materials. Various experimental and numerical studies show that addition of liquid bridge forces leads to higher yield strength. The yield stress at critical state can be fitted as a linear function of the pressure with the friction coefficient of dry flow m o as the slope and a finite offset c, defined as the steady state cohesion in the limit of zero confining pressure [35]. This finite offset c is constant in the high pressure limit. However, very little is known regarding the rheology for granular materials in the low pressure limit. Depending on the surrounding conditions, granular flows phenomenon are affected by appropriate time scales namely, t p : time required for particles to rearrange under certain pressure,ġ t : time scale related to strain rateġ , t k : related to the contact time between particles, t g : elapsed time for a single particle to fall through half its diameter under the influence of gravity and t c : time scale for the capillary forces driving the flow are primarily hindered by inertia based on particle density. While various time scales, as related to the ongoing mechanisms in the sheared bulk of the material, can interfere, they also can get decoupled, in the extremes of the local/ global condition, if one time scale gets way smaller in magnitude than the other. A detailed description of this time scales are given in section 3. While t k , t g and t c are global, other time scalesġ t and t p depends on local field variables strain rateġ and pressure p respectively. We restrict our studies to the quasi-static regime (˙ g t t p ) as the effect of cohesion decreases with increasing inertial number due to the fast decrease in coordination number [1]. Moreover, the quasistatic regime observed for non-cohesive particles also persist for cohesive particles, while the inertial regime of noncohesive particles bifurcates into two regimes: rate-independent cohesive regime at low shear rates and inertial regime at higher shear rates [11]. In the present work, we shed light on the rheology of non-cohesive dry as well as cohesive wet granular materials at the small pressure limit, by studying free surface flow. While the inertial number I [19], i.e. the ratio of confining pressure to strain-rate time scales, is used to describe the change in flow rheology from quasi-static to inertial conditions, we look at additional dimensionless numbers that influence the flow behavior. (i) The local compressibility * p , which is the squared ratio of the softness and stress time scales (ii) the inverse relative pressure gradient * p g , which is the squared ratio of gravitational and stress time scales and (iii) the Bond number Bo [48] quantifying local cohesion as the squared ratio of stress to wetting time scales are these dimensionless numbers. We show a constitutive relation based on these dimensionless numbers in sections 4-6 of this paper. Additional relevant parameters are not discussed in this study, namely granular temperature or fluidity. All these dimensionless numbers can be related to different time scales or force scales relevant to the granular flow.
Granular materials display non-Newtonian flow behavior for shear stresses above the so called yield stress while they remain mostly elastic like solids below this yield stress. More specifically, granular materials flow like a shear thinning fluid under sufficient stress. When dealing with wet granular materials, it is therefore of fundamental interest to understand the effect of cohesion on the bulk flow and yield behavior. Recently, the majority of investigations of non-Newtonian flow behavior focused on concentrated colloidal suspensions. Shear thickening is often observed in those flows due to the formation of flow-induced density fluctuations (hydroclusters) resulting from hydrodynamic lubrication forces between particles [46]. Similar local clusters (aggregates) can also be found in strongly cohesive wet granular materials, especially near to the free surface, where attractive forces dominate their repulsive counterparts [39]. However, the strong correlations observed between particles of close proximity in suspensions seem to be irrelevant in wet granular systems, where the range of force interactions is much more limited. On the other hand, Lin et al [22] show that contact forces dominate over hydrodynamic forces in suspensions that show continuous shear thickening. Fall et al [7] propose that discontinuous shear thickening of cornstarch suspensions is a consequence of dilatancy: the system under flow attempts to dilate but instead undergoes a jamming transition because it is confined-a phenomenon that was recently also explained by a moving jamming point [20]. Another possible cause for shear thickening is the large stress required to maintain flow due to particle-particle friction above a critical stress as in [6,29]. This is more likely to happen in charge stabilized colloidal suspensions. Here we only intended to speculate the flow behavior of cohesive granular materials in relevance to micro scale analogy for shear thickening in suspensions and section 7 of this paper is devoted to understand more on the behavior of wet granular materials with increasing cohesion.  [42,50], an open-source implementation of the discrete particle method, to simulate a shear cell with annular geometry and a split bottom plate, as shown in figure 1. Some of the earlier studies in similar rotating set-ups include [37,47,52]. The geometry of the system consists of an outer cylinder (outer radius R o =110 mm) rotating around a fixed inner cylinder (inner radius R i =14.7 mm) with a rotation frequency of Ω=0.01 revolutions per second. The granular material is confined by gravity between the two concentric cylinders, the bottom plate, and a free top surface. The bottom plate is split at radius R s =85 mm. Due to the split at the bottom, a narrow shear band is formed. It moves inwards and widens towards the flow surface. This set-up thus features a wide shear band away from the bottom and the side walls which is thus free from boundary effects. The filling height (H=40 mm) is chosen such that the shear band does not reach the inner wall at the free surface.

Model
In earlier studies [33,39,40], a quarter of this system ( • 0  f  • 90 ) was simulated using periodic boundary conditions. In order to save computation time, here we simulate only a smaller section of the system with appropriate periodic boundary conditions in the angular coordinate, unless specified otherwise. We have observed no noticeable effect on the macroscopic behavior in comparisons between simulations done with a smaller ( • 30 ) and a larger ( • 90 ) opening angle. Note that for very strong attractive forces, agglomeration of particles occur. Then, a higher length scale of the geometry is needed and thus the above statement is not true anymore.

Contact model and parameters
The liquid bridge contact model is based on a combination of an elastic-dissipative linear contact model for the normal repulsive force and a nonlinear irreversible liquid bridge model for the non-contact adhesive force as described in [35]. The adhesive force is determined by three parameters; surface tension σ, contact angle θ which determine the maximum adhesive force and the liquid bridge volume V b which determines the maximum interaction distance between the particles at the point of bridge rupture. The contact model parameters and particle properties are as given in ). To study the effect of inertia and contact stiffness on the non-cohesive materials rheology, we compare our data for non-cohesive case with data from simulations of [40] for different gravity as given below:  We also compare the effect of different rotation rates on the rheology for the following rotation rates: The liquid capillary force is estimated as stated in [51]. It is observed in our earlier studies [35] that the shear stress τ for high pressure can be described by a linear function of confining pressure, p, as t m = + p c o . It was shown that the steady state cohesion c is a linear function of the surface tension of the liquid σ while its dependence on the volume of liquid bridges is defined by a cube root function. The friction coefficient m o is constant and matches the friction coefficient of dry flows excluding the small pressure limit. In order to see the effect of varying cohesive strength on the macroscopic rheology of wet materials, we vary the intensity of capillary force by varying the surface tension of the liquid σ, with a constant volume of liquid bridges The first case, s = 0.0 N m −1 , represents the case of dry materials without cohesion, whereas s = 0.50 N m −1 corresponds to the surface tension of a mercury-air interface. For s > 0.50 N m −1 , smooth, axisymmetric shear band formation is not observed and the materials agglomerate to form clusters as shown in figure 2, for our particle size and density. Hence, σ is limited to maximum of 0.50 N m −1 .

Averaging methodology
To extract the macroscopic properties, we use the spatial coarse-graining approach detailed in [24][25][26]. The averaging is performed over a grid of 47-by-47 toroidal volumes, over many snapshots of time assuming rotational invariance in the tangential f-direction. The averaging procedure for a three-dimensional system is explained in [24,26]. This spatial coarse-graining method was used earlier in [26,33,39,40,52]. We do the temporal averaging of non-cohesive simulations over a larger time window from 30 to 440 s with 2764 snapshots to ensure the rheological models with enhanced quality data. All the other simulations are run for 200 s and temporal averaging is done when the flow is in steady state, between 80 and 200 s with 747 snapshots, thereby disregarding the transient behavior at the onset of the shear. In the critical state, the shear band is identified by the region having strain rates higher than 80% of the maximum strain rate at the corresponding height. Most of the analysis explained in the later sections are done from this critical state data at the center of the shear band.

Macroscopic quantities
The general definitions of macroscopic quantities including stress and strain rate tensors are included in [40]. Here, we define the derived macroscopic quantities such as the friction coefficient and the apparent viscosity which are the major subjects of our study.
The local macroscopic friction coefficient is defined as the ratio of shear to normal stress and is defined as m t = p.
The magnitude of strain rate tensor in cylindrical polar coordinates is simplified, assuming u r =0 and u z =0: The apparent shear viscosity is given by the ratio of the shear stress and strain rate as: whereġ is the strain rate.

Critical state
We obtain the macroscopic quantities by temporal averaging as explained in section 2.3. Next we analyze the data, neglecting data near walls ( < » r r 0.045 figure 3. Further, the consistency of the local averaged quantities also depends on whether the local data has achieved the critical state. The critical state is defined by the local shear accumulated over time under a constant pressure and constant shear rate condition. This state is reached after large enough shear, when the materials deform with applied strain without any change in the local quantities, independent of the initial condition. We focus our attention in the region where the system can be considered to be in the critical state and thus has a well defined macroscopic friction. To determine the region in which the flow is in critical state,˙( ) g z max is defined to be the maximum strain rate for a given pressure, or a given height z. The critical state is achieved at a constant pressure and strain rate condition over regions with strain rate larger than the strain rate˙( ) g z 0.1 max as shown in figure 3 corresponding to the region of shear band. While [40] showed that for rotation rate 0.01 rps, the shear band is well established above shear rateġ > 0.01 s −1 , of our analysis shown in the latter sections are in the shear band center is obtained by˙˙( ) g g > z 0.8 max at different heights in the system. This is defined as the region where the local shear stress τ becomes independent of the local strain rateġ and t p becomes constant. We also extend our studies to the shear-rate dependence in critical state which is effective for critical state data for wider regions of shear band (section 4.4). This shear rate dependence is analyzed in the regions of strain rate (ġ ) larger than the˙( ) g z 0.1 max at a given height z. These data include the region from the center to the tail of the shear band, with typical cut-off factors s c =0.8 or 0.1, respectively, as shown in figure 3, and explained in section 4.4.

Time scales
Dimensional analysis is often used to define the characteristic time scales for different physical phenomena that the system involves. Even in a homogeneously deforming granular system, the deformation of individual grains is not homogeneous. Due to geometrical and local parametric constraints at grain scale, grains are not able to displace as affine continuum mechanics dictates they should. The flow or displacement of granular materials on the grain scale depends on the timescales for the local phenomena and interactions. Each time scale can be obtained by scaling the associated parameter with a combination of particle diameter d p and material density ρ. While some of the time scales are globally invariant, others are varying locally. The dynamics of the granular flow can be characterized based on different time scales depending on local and global variables. First, we define the time scale related to contact duration of particles which depends on the contact stiffness k as given by [40]: In the special case of a linear contact model, this is invariant and thus represents a global time scale too. Two other time scales are globally invariant, the cohesional time scale t c , i.e. the time required for a single particle to traverse a length scale of d 2 p under the action of an attractive capillary force and the gravitational time scale t g , i.e. the elapsed time for a single particle to fall through half its diameter d p under the influence of the gravitational force. The time scale t c could vary locally depending on the local capillary force f c . However, the capillary force is weakly affected by the liquid bridge volume while it strongly depends on the surface tension of the liquid σ. This leads to the cohesion time scale as a global parameter given by: with surface tension σ and capillary force ps The corresponding time scale due to gravity which is of significance under small confining stress close to the free surface is defined as: The global time scales for granular flow are complemented by locally varying time scales. Granular materials subjected to strain undergo constant rearrangement and thus the contact network re-arranges (by extension and compression and by rotation) with a shear rate time scale related to the local strain rate field: Finally, the time for rearrangement of the particles under a certain pressure constraint is driven by the local pressure p. This microscopic local time scale based on pressure is: As the shear cell has an unconfined top surface, where the pressure vanishes, this time scale varies locally from very low (at the base) to very high (at the surface). Likewise, the strain rate is high in the shear band and low outside, so that also this time scale varies between low and high, respectively. Dimensionless numbers in fluid and granular mechanics are a set of dimensionless quantities that have a dominant role in describing the flow behavior. These dimensionless numbers are often defined as the ratio of different time scales or forces, thus signifying the relative dominance of one phenomenon over another. In general, we expect five time scales (t g , t p , t c ,ġ t and t k ) to influence the rheology of our system. Note that among the five time scales discussed here, there are ten possible dimensionless ratios of different time scales. We propose four of them that are sufficient to define the rheology that describes our results. Interestingly, all these four dimensionless ratios are based on the common time scale t p . Thus, the time scale related to confining pressure is important in every aspect of the granular flow. All the relevant dimensionless numbers in our system are discussed in brief in the following two sections of this paper for the sake of completeness, even though not all are of equal significance.

Rheology of dry granular materials 4.1. Effect of softness in the bulk of the materials
We study here the effect of softness on macroscopic friction coefficient for different gravity in the system. Thus the pressure proportional to gravity is scaled in dimensionless form * p [40] given by: This can be interpreted as the square of the ratio of time scales, = * p t t k p 2 2 , related to contact duration and pressure respectively. Figure 4 shows the macroscopic friction coefficient as a function of the dimensionless pressure * p and the dashed line is given by: * p o denotes the limiting dimensionless pressure around the correction due to softness of the particles, where the correction is not applicable anymore, since  f 0 [27]. We have used this fit, as our data range is too limited to derive the functional form of the fit. This is shown by the solid line in figure 4 with the plotted data from our present simulation (&) and with data for different gravity in the system [40] which we use to describe other corrections for dry non-cohesive materials. Despite the deviation of data for different gravity from the trend for small * p , the agreement with our data is reasonable. The dashed line represents the softness correction as proposed by [40]. The effect of softness is dominant in regions of large pressure where the pressure time scale t p dominates over the stiffness time scale t k and thus the data in plot are corresponding to higher than a critical pressure ( > * p 4 g , explained in section 4.3). Here, the compressible forces dominate over the rolling and sliding forces on the particles, the flow being driven by squeeze. Thus, the macroscopic friction coefficient decreases with softness.

Effect of inertial number
For granular flows, the rheology is commonly described by the dimensionless inertial number [30]:   [27]. Note that these fitting constants change with the range of I that are included in the fitting. Given that we do not have data

Effect of gravity close to the free surface
In this section, we investigate the effect of the another dimensionless number * p g on local friction coefficient, given by: This can be interpreted as the square of the ratio of time scales, = * p t t g g p 2 2 , related to gravity and pressure respectively. The effect of inertial number and softness correction are eliminated by scaling μ by the correction factors m I and f p respectively and studying the effect of * p g on the scaled friction coefficient. Figure 6 shows μ scaled by m f I p as a function of dimensionless pressure * p g for different gravity g (different * p ) and different rotation rates Ω (different I), including our data for g=9.81 ms −2 and W = 0.01 rps which is also in agreement with other data set. The data for different slower rotation rates and different gravitational accelerations g agree well with our new data set, while the higher rotation rates deviate. Note that the higher rotation rates are in a different regime where kinetic theory works and hence agreement with the generalized rheology is not expected strictly. All the data for different gravity and slower rotation rates collapse and these can be fitted by the solid line given by the correction ( ) where, ¢ » a 0.71 is the relative drop in friction coefficient at go is the dimensionless pressure at which the friction coefficient drops below m 0.74 o and ( ) * f p g g is the correction corresponding to the dimensionless pressure * p g . Due to lack of confining stress close to the free surface ( < * p 4 g ), the macroscopic friction coefficient exponentially decreases with decrease in * p g . Here, the gravity time scale t g dominates over the pressure time scale t p . Thus, while the effect of gravity close to the free surface is dominant for is the critical pressure above which the effect of softness p * is significant as explained in section 4.1.

Shear rate dependence in critical state flow
After having quantified the dependence of the macroscopic friction on inertial number and softness, another correction was proposed in [40], taking into account a reduced, relaxed friction correction in very slow quasistatic flow.  is a representation of this correction f q (I) where: 4.85 1.08 10 5 for very small inertial numbers (  * I I ) and a =  0.48 0.07 1 . This correction is in inspiration with [24] where I * scales linearly with the external shear rate and thus is proportional to the local strain-rate and the granular temperature. Although the data represented in ) which corresponds to the data given by red • which are invariant to the effect of small inertial number which allows us to assume ( ) » f I 1.0 q . Hence, in the following sections, we do not take into consideration the correction f q (I), though we mention it.

Bond number
The Bond number (Bo) is a measure of the strength of the adhesive force relative to the compressive force. A low value of Bo (typically much less than 1) indicates that the system is relatively unaffected by the attractive forces; high Bo indicates that the attractive force dominates in the system. Thus Bo is a critical microscopic parameter that controls the macroscopic local rheology of the system. While the conventional way of defining the Bond number as the ratio of the time scales t c and t g [48] is appropriate for single particles, or close to the free surface, we define the local Bond number relative to the confining force: where, p mean is the mean pressure in the system. This is an experimentally measurable quantity and is related to quantifying the system as a whole. The global Bond number corresponding to surface tension of liquid defined in equation (3) is given by:

Effect of local bond number
The properties of the particles and the interstitial fluid strongly affect the macroscopic behavior of granular materials. The local macroscopic friction is studied as a function of local Bond number Bo for different wet cohesion intensity. Figure 8 shows the macroscopic friction coefficient as a function of the local Bond number Bo for different wet cohesion. It is evident that the friction coefficient increases with local Bond number with a constant value m o in the low Bond number limit. For frictionless wet cohesive materials, the rheology can be defined by a linear fitting function given by: is the macroscopic friction coefficient in the high pressure limit [35] and » a 1.47. This is shown by the solid line in figure 8. However, it is observed that the data deviate from the solid fitting line in the high Bond number or low pressure limit. This deviation is explained by the small pressure correction ( ) * f p g g as explained in section 4.3 and discussed in details in the next section. Figure 6 shows the dependence of the local friction coefficient on the local scaled pressure * p g for dry noncohesive materials and this effect is small in the high pressure limit. With an attempt to separate the effect of Bond number on the rheology of cohesive materials, we plot the local friction coefficient μ scaled by the Bond number correction f c and other corrections m I and f p , as a function of scaled pressure * p g as shown in figure 9. The solid line is given by equation (23), where the non-cohesive function fits for the wet data as well.

Rheological model
We studied the rheology of dry and wet granular materials in terms of different dimensionless numbers. The trends are combined and shown to collectively contribute to the rheology as multiplicative functions given by: The proposed general(ized) multiplicative rheology function for the macroscopic friction coefficient is dependent on four dimensionless numbers * * p p I , , g and Bo. Table A1 in the appendix gives the summary and details of our proposed rheological model. This rheological model is based on constant liquid bridge volume at all contacts and we do not take into account liquid redistribution among contacts [28,36]. This is a simplified approach to establish the generalized rheology and we are working further on liquid redistribution and will analyze its effect on the rheology. However, the cohesion time scale is only weakly affected by the liquid bridge volume and mainly depends on the surface tension of the liquid. Preliminary results using a liquid redistribution model show that in this state, 40% of the contacts in the shear band center become dry, resulting in a higher probability of dry contacts with microcontact local Bond number = Bo 0. This results in a lower local Bond number in the shear band center. Our present rheological model is shown to be valid for a wide range of Bond number and thus use of a liquid redistribution model is expected to shift data further, towards the lower Bond numbers but is expected to follow the same trends.
For a full constitutive law, one also needs to take into account the solid volume fraction also. For dry granular shear flow [27,40], the constitutive relations for the volume fraction given by corrections (to first order) based on dimensionless numbers * p and I as follows: where, f » 0.65 c is the critical or the steady state density under shear, in the limit of vanishing pressure and inertial number. = I 0.85 c is the inertial number corresponding to strain rate when the dilation turns to fluidization. = * p 0.33 c is the typical pressure for which softness leads to huge densities. Though the volume fraction in an inhomogeneous system is a field (fluctuating around a mean value), its local values are captured by the above equation in terms of the local dimensionless numbers. The above relation shows that the volume fraction decreases (and the friction increases) when the quasi-static regime is exceeded. However, the generalized rheology is expected to be valid everywhere in the inhomogenous system where the system has been sheared long enough to reach the critical state, irrespective of their different volume fraction. The volume fraction increases with increase in confining stress as shown in [27,40]. In ongoing research [34], we show that inter-particle cohesion has a considerable impact on the compaction of the soft materials. Cohesion causes additional stresses, due to capillary forces between particles, leading to an increase in volume fraction due to higher compaction. This effect is not visible in a system of infinitely stiff particles. On the other hand, we observe a general decrease in volume fraction due to increased cohesion, which we attribute to structural changes in the bulk material.

Local apparent viscosity
For unsaturated granular materials, being heterogeneous systems, it is not relevant to define their viscosity. Nevertheless, we introduce the local apparent viscosity η of granular materials which is barely the ratio of the shear stress to the strain rate as an alternative to μ. To see the combined effect of pressure and strain rate on the local apparent viscosity, we analyze them as functions of the inertial number. For a given pressure, the inertial number is proportional to the shear rate. Thus, the analysis of local apparent viscosity as a function of the inertial number for small pressure ranges can be interpreted as the analysis of apparent viscosity versus strain rate. We Since we here focus on the data in the center of the shear band, the dependence on shear rate in the critical state flow which includes data outside the shear band center can be neglected ( ( )»  * f I I 1 q ) and thus the rheological model for the local friction coefficient given by equation (22) is simplified by: The dimensionless variable h* can be related to three time scales namely, contact duration t k , strain rate related time scaleġ t and pressure related time scale t p asḣ m = g * t t t k p 2 . Alternatively, the flow rules of granular materials can be approximated as that of a power-law fluid as given by: where, m = a -* K p I is the flow consistency and α is the flow behavior index. The flow rules of granular materials are pretty straightforward at high pressures with a » 0. However, deviations are observed from the power-law behavior at small pressures. More details on the flow rules at large and small pressure are explained in sections 7.1.2 and 7.1.3 respectively. Figure 10 shows the local apparent viscosity h* as a function of the inertial number I for different global Bond numbers. The data shown correspond to all the data close to the shear band center for different heights. The inertial number is lowest at an intermediate height, and increases towards surface and base. With increasing inertial number, the apparent shear viscosity decreases, indicating that granular materials flow like non-Newtonian fluids, specifically shear-thinning fluids. It is also evident from the figure that the flow behavior is different at large and small confining pressure. We use this simplified constant strain rate to predict the apparent viscosity near the surface of the shear cell where the pressure is very small. The apparent shear viscosity for wet cohesive materials confined to small pressure is more intricate and is predicted by from equations (24) and (25) Figure 10 shows the prediction of apparent viscosity at small pressure as given by the green solid lines. Noncohesive materials upto weakly cohesive materials ( < Bo 0.60 g ), at low pressure, are less viscous than those at high pressure, as shown in the figure. For global Bond number = Bo 0.60 g , materials for a given inertial number have the same apparent viscosity independent of pressure. For even higher cohesion ( > Bo 0.60 g ), the flow behavior changes qualitatively. Though, the apparent viscosity decreases with the inertial number (h µ d -* I ), even for cohesive materials, the qualitative decay power δ decreases towards zero (d  0). For a given inertial number, the material near the surface has higher apparent viscosity than in the bulk and at the base. Materials confined by small pressure thus display reduced shear thinning with increase in cohesion. This is represented by the direction of black arrows marked with Bo g in figure 10. Thus, granular materials have different shearthinning properties depending on local confining pressure and Bond number.

Analytical prediction of apparent viscosity
We extract the position and the width of the shear band R c and W respectively from the fit function in equation (28). Both position and width of the shear band depend on the height in the system and the position moves inwards with increasing height (decreasing pressure). Predictions of the position of the shear band center as a function of height is given in [44]. Since the analytical prediction discussed here is not significantly affected by this varying position of the shear band, we use the mean shear band positionR c for our prediction. The shear band moves inward with increase in global Bond number [39]. Thus the mean shear band positionR c decreases with increasing Bo g (not shown here).
The width of the shear band is predicted as function of height as given by [32]: where b = 0.6 for non-cohesive materials and b < < 0.5 0.7 for cohesive materials are fitted well by our data. Assuming the pressure varying hydrostatically and the bulk density as r r » 0.6 b , we translate equation (34) to W as a function of p. Substituting equations (31) and (34) in (13) and rearranging, we get the inertial number I max in the shear band center as a function of the local pressure p. Further, by substituting p, we get h* max in the shear band center and thus obtain a quantitatively accurate prediction of h* max (I max ), plotted as blue solid lines and cyan dashed lines in figure 10.
The results show that the analytical solution is in good agreement with our numerical results. Focusing on the slope of the small pressure line, we observe that it changes with increasing cohesion in the same way as shown by numerical data. It is observed from the analytical solution that this change in slope is governed by m. Thus, the shear-thinning rate for materials under small pressure depends on local friction coefficient, which depends on the corrections ( ) * f p g g and f c (Bo).

Eliminating the effect of cohesion and gravity
Under larger confining pressure (as stated in section 7.1.2), with increase in cohesion, the apparent viscosity of the granular fluid increases, however, the flow behavior remains qualitatively the same even for very high cohesion. For materials confined to large pressure, where p is slowly varying, the apparent viscosity is inversely proportional to the strain rate and approximately also to the inertial number. At smaller pressure, the materials are more free only under the effect of gravity, with less dominant forces due to particle contacts. Therefore, cohesion is relatively more dominant for higher local Bond numbers, resulting in the qualitative change in shear thinning rate (α). Thus the flow of materials is affected by both dimensionless numbers Bo and * p g at the same time. Then, the granular fluid appears to no longer behave like a power-law fluid. Several of these rheological correction factors make the flow behavior even more nonlinear under small pressure. In order to see the rheology of the granular fluid under small pressure, which is devoid of the effect of these dimensionless numbers, we rescale the local dimensionless apparent viscosity h* by ( ) f Bo c and ( ) * f p g g and analyze it as a function of inertial number. Figure 11(a) shows the dimensionless apparent viscosity h* scaled by f c (Bo) as a function of inertial number for different cohesion. All the data for different cohesion collapse to a single plot for the triad of different pressure scales. Further, we rescale and plot it as a function of inertial number for different cohesion as shown in figure 11(b). The fitted solid line corresponding to the data at large pressure is given by equation (26) with a = 0 and » K 0.01. Furthermore, the fitted dashed line corresponding to the data at small pressure is given by equation (26) with a = -1 and »´-K 5.6 10 6 . This is explained theoretically by substituting p * in equation (13) and using equation (33) with constant friction coefficient m 0 yielding: Thus, for slowly varying strain rate at small pressure, h* is proportional to -I 2 and is represented by equation (26) with a = -1. This eventually explains the earlier observations in [25].
Thus, the flow behavior for granular materials in a simple hypothetical case with high confining stress constant friction coefficient can be approximated by that of a power-law fluid flow behavior. However, for more realistic systems, e.g., unit operations at low stress, several other factors influence the flow rheology, e.g., near to the free surface. Thus, under small pressure, granular materials behave more interestingly and complex than a power-law fluid.

Discussions and conclusions
The rheology of dry as well as wet granular materials (in the pendular regime) has been studied by simulations using the discrete element method in steady state shear. Our results show that the conventional ( ) m I rheology must be modified to take into account other factors such as cohesion, contact softness, corrections at small pressures where gravity dominates, and a generalized inertial number dependence for very slow quasi-static flow (creep) in the tails of the shear bands. The trends are combined and shown to collectively contribute to the rheology as multiplicative functions, i.e. ignoring one contribution can lead to inconsistent results. This new generalized rheological model applies to a wide range of parameters from dry non-cohesive to strongly cohesive materials, and contains also both the small and the large pressure limits. Note that additional contributions from viscous forces should be included in case of rapid flow. Our ongoing work shows that the generalized rheology is independent of system configuration, pressure or volume control, in the critical state and is applicable for both homogeneous simple shear and inhomogeneous systems like the split-bottom shear cell. Given this is justified, the shear thinning behavior for granular materials is valid for every locally reached critical state, irrespective of the system configuration for moderate to low pressure and in the dense regime. Furthermore, we study the apparent viscosity as a function of inertial number for granular fluids of varying cohesive strength. Most strikingly, the cohesive strength not only increases the magnitude of the apparent viscosity, but also decreases the shear thinning rate, but only for material under small confining pressure e.g. close to the free surface. This variable shear thinning behavior of granular materials, close to a free surface, is attributed to the higher local Bond number i.e. it is a low pressure effect. Thus, the flow rheology (friction and apparent viscosity) is predicted by the proposed rheology model for dry and wet granular materials under both low and high confining stress. Further, we develop an analytical solution for the apparent viscosity using the proposed rheology (with some simplifications) and show that the results are in good agreement with our numerical analysis. Materials have less shear thinning with an increase in cohesion as quantified by the high Bond numbers under small confining pressure.
Finally, it is shown that the effect of each of the dimensionless numbers can be eliminated by rescaling, and thus the scaled apparent viscosity of a simple system with a (small) constant friction coefficient is predicted as that of a power-law fluid with Bagnold type scaling withI .
As an outlook, we aim to implement the generalized rheological model in a continuum description of the splitbottom shear cell geometry. A successful implementation is only the first step for validation and paves the way to use this rheological model in industrial applications for material flow descriptions. We aim to also include higher order effects of the Bond number in the generalized rheology. We included the small pressure (free surface) correction in the rheology, as an effect of gravity. It is to be noted that even in a micro-gravity system, both pressure and gravity change identically and thus the corresponding correction term remains the same as in a system with high gravity. Thus this correction corresponds to an effect active at interfaces or at the free-surface. Next step is to perform the micro-structural analysis [39] also for our system and in particular close to the free surface in order to understand the change of the shear thinning rate. Another open question concerns the creep correction and its relation to the microstructure and granular temperature. Last, the present rheology has to be merged to kinetic theory in the rapid, collisional flow regime [45], which presents another open challenge. Table A1. List of rheological correction functions for application in a continuum model, see equation (22).

Dimensionless numbers
Corrections Coefficients from fits Coefficients in [27] Inertial