Steady state rheology of homogeneous and inhomogeneous cohesive granular materials

This paper aims to understand the effect of different particle/contact properties like friction, softness and cohesion on the compression/dilation of sheared granular materials. We focus on the local volume fraction in steady state of various non-cohesive, dry cohesive and moderate to strong wet cohesive, frictionless-to-frictional soft granular materials. The results from (1) an inhomogeneous, slowly sheared split-bottom ring shear cell and (2) a homogeneous, stress-controlled simple shear box with periodic boundaries are compared. The steady state volume fractions agree between the two geometries for a wide range of particle properties. While increasing inter-particle friction systematically leads to decreasing volume fractions, the inter-particle cohesion causes two opposing effects. With increasing strength of cohesion, we report an enhancement of the effect of contact friction i.e. even smaller volume fraction. However, for soft granular materials, strong cohesion causes an increase in volume fraction due to significant attractive forces causing larger deformations, not visible for stiff particles. This behaviour is condensed into a particle friction—Bond number phase diagram, which can be used to predict non-monotonic relative sample dilation/compression.

Material's bulk responses like density/shear resistance are influenced strongly by different particle properties such as frictional forces, as well as dry or wet cohesion. How the microscopic particle properties influence the granular rheological flow behaviour is still a mystery and thus attracted more attention in the last few decades [15-17, 20, 30, 36, 42, 51, 52, 63]. Although the influence of individual particle properties is better understood now, there is still very little known about the combined effect of particle friction and cohesion on the rheological behaviour of granular flows [48].
For dry cohesive granular media, one needs to account for the dominant van der Waals force between particles when they are not in contact. This attractive force can be modeled by a reversible linear long range interaction [54], if the interaction energy and dissipation are matched with the true non-linear force.
In contrast, in wet granular media, particles attract each other, as liquid bridges cause capillary forces [66,74,78]. The capillary bridge force becomes active at first contact, but then is active up to a certain cut-off distance where the bridge ruptures. This gives asymmetric loading/unloading behaviour. Recent studies based on Discrete Element Method (DEM) confirmed that the specific choice of capillary bridge models (CBMs) has no marked influence on the hydrodynamics of granular flow [11,24] for small volumes of interstitial liquids. In the present paper, our focus is thus to investigate the effect of two qualitatively different cohesion mechanisms, dry vs. wet, on the bulk density of cohesive granular materials. We quantify the cohesion associated with different contact models by the dimensionless Bond number, Bo, which is defined as the ratio of time scales related to confining stress and adhesive force. In this way, we generalize the existing rheological models with cohesion dependencies.
For dry, non-cohesive but frictional granular materials, the bulk density under steady shear decreases [17]. Cohesive grains are more sensitive to stress intensity as well as direction, and exhibit much larger variations in their bulk densities [55]. Dry granular media with median particle size below 30 μ m flow less easily and show a bulk cohesion due to strong van der Waals interactions between particle pairs [60]. Unsaturated granular media with interstitial liquid in the form of liquid bridges between particle pairs can also display bulk cohesion which can be strongly influenced by the flow or redistribution of the liquid [47,53]. Fournier et al. [10] observed that wet systems are significantly less dense than dry granular materials even for rather large particles. The packing density is only weakly dependent on the amount of wetting liquid, because the forces exerted by the liquid bridges are very weakly dependent on the bridge volume [78]. In general, one expects that the steady state local volume fraction decreases with increasing either particle friction or cohesion [16,69]. This decrease of volume fraction is related to both structural changes and increasing bulk friction of the materials [64]. However, Roy et al. [51] reported an opposing increase of local steady state volume fraction proportional with cohesive strength. Thus whether the particle cohesion enhances or suppresses the frictional effects that causes dilatancy or compression remains still a debate. Therefore, the second focus of this paper is to understand the interplay between particle friction and cohesion on the steady state volume fraction of the granular systems under shear.
Dilatancy, or shear dilatancy in classical plasticity and soil mechanics, has a specific definition as the ratio of incremental volumetric strain to the shear strain. While the "relative dilatancy" discussed in this work is a dynamic effect that is distinct from the classical dilatancy since in the steady simple shear flow the volumetric strain is zero. Choosing different dynamic steady states based on different inter-particle cohesion, it represents the bulk density change between two different steady states.
The work is structured as follows: in Sect. 2.1, we provide information on the two simulation geometries; in Sect. 2.2, the two cohesive contact models are introduced; in Sect. 2.3 the important dimensionless number and their related time scales are elaborated; and in Sect. 2.4, the input parameters are given. Section 3 is devoted to the discussion of major findings of rheological modelling with a focus on the combined influences of several particle parameters, while conclusions and outlook are presented in Sect. 4.

Simulation methods
The Discrete Particle Method (DPM) or Discrete Element Method (DEM) is a family of numerical methods for simulating the motion of large numbers of particles. In DPM, the material is modeled as a finite number of discrete particles with given micro-mechanical properties. The interactions between particles are treated as dynamic processes with states of equilibrium developing when the internal forces balance. As previously stated, granular materials are considered as a collection of discrete particles interacting through contact forces. Since the realistic modeling of the shape and deformations of the particles is extremely complicated, the grains are assumed to be non-deformable spheres which are allowed to overlap [32]. The general DPM approach involves three stages: (1) detecting the contacts between elements; (2) calculating the interaction forces among grains; (3) computing the acceleration of each particle by numerically integrating Newton's equations of motion after combining all interaction forces. This three-stage process is repeated until the entire simulation is complete. Based on this fundamental algorithm, a large variety of codes exist and often differ only in terms of the contact model and some techniques used in the interaction force calculations or the contact detection. After the discrete simulations are finished, there are two popular ways to extract relevant continuum quantities for flow description such as stress, bulk density from discrete particle data. The traditional one is ensemble averaging of "microscopic" simulations of homogeneous small samples, a set of independent RVEs [19,27,40]. A recently developed alternative is to simulate an inhomogeneous geometry where dynamic, flowing zones and static, high-density zones coexist. By using adequate local averaging over a fixed volume (inside which all particles can be assumed to behave similarly), continuum descriptions in a certain parameter range can be obtained from a single set-up [30,35,64,65].

Geometries
In this study, we use MercuryDPM [68,76,77], an opensource implementation of the Discrete Particle Method (DPM) to simulate granular flow in two geometries: a homogeneous stress controlled simple shear box (RVEs) and an inhomogeneous split bottom shear cell, where the local stress is given by the weight of the particles above.

Stress controlled simple shear box (SS)
The first simple geometry is a cuboid shear volume with Lees-Edwards periodic boundaries in x and y directions [25] and normal periodic boundaries in z direction as shown in Fig. 1. The initial length of each box side is set to L, but the box side length in y-direction, L y , can vary in time. The polydispersed granular sample contains 4096 soft particles. The system is sheared along the x direction with a constant shear rate ̇= 2V x ∕L y 1 by moving all the particles at each timestep to achieve homogeneous shearing. Meanwhile, the normal stress along y direction, yy is kept constant by allowing L y to dilate/compact, so that it smoothly reaches its steady state [59,61]. Using this setup, one can keep both shear rate and normal stress constant while measuring the shear stress responses in both transient and steady state. Although this setup is not achievable in reality, it still represents typical lab experiments, sand or/and powders sheared in different shear cells. This setup allows the user to explore two variables (shear rate and stress) independently with low computational cost. Therefore, in the current study, we systematically vary the shear rate ̇= V x ∕ y and confining stress yy to understand the shear flow in both quasi-static and dynamic regimes, as well as the influence of particle Softness.

Split-bottom ring shear cell (SB)
In our present study, we also simulate a shear cell with annular geometry and a split bottom, as described in detail in [50,64]. The geometry of the system consists of an outer cylinder (radius R o = 110 mm ) rotating around a fixed inner cylinder (radius R i = 14.7 mm ). We vary the rotation frequency from = 0.06 to 4.71 rad/s (0.01 to 0.75 rotation/s). The granular material is confined by gravity between the two concentric cylinders, the bottom plate, and a free top surface. The gravity varies from 1 to 50 ms −2 . 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 free surface. The filling height ( H ≈ 40 mm) is chosen such that the shear band does not reach the inner wall. Figure 2 shows the 3D schematic presentation and side view of the split bottom shear cell geometry with colors blue to red indicating low to high kinetic energy of the particles. It is visible that a wide shear band is formed away from the walls, which is thus free from boundary effects. In earlier studies [64,65], a quarter of this system ( 0 • ≤ ≤ 90 • ) was simulated using periodic boundary conditions. All the data corresponding to different gravities and different rotation rates belong to this system. In order to save computation time, we simulate only a smaller section of the system ( 0 • ≤ ≤ 30 • ) with appropriate periodic boundary conditions in the angular coordinate. No noticeable differences in the macroscopic flow behavior were observed between simulations performed using a smaller ( 30 • ) and a larger ( 90 • ) section.

Force models
In our non-cohesive simulations, we use the linear viscoelastic frictional contact model between particle contacts in both geometries mentioned above: Simple Shear Box (SS) and Split Bottom Shear Cell (SB) [31].
The dry cohesive particles are simulated in the geometry SS using a linear reversible adhesive visco-elastic contact model as shown in Fig. 3a , where we combine Simulation of a 3D system of polydisperse particles simple shear using Lees-Edwards periodic boundary conditions, with L y controlled to keep the normal stress yy constant a linear attractive non-contact force with a linear viscoelastic frictional contact force [31] as the sum of Eqs. (2) and (3). This reversible attractive model is chosen as a (linearized) representation of the van-der-Waals attractive non-contact force.
And in the geometry SB, we simulate wet cohesive particles using an irreversible adhesive non-contact force with a linear visco-elastic frictional contact force as shown in Fig. 3b, which only differs in the attractive non-contact part compared to the dry cohesive simulations as shown in Eq. (4).

Classical linear contact force
When two particles i, j are in contact, the overlap can be computed as, with positions i and j , radii a i and a j , for the two primary particles, respectively and the unit vector Thus the normal force f n between two particles is simply computed when > 0 (Fig. 3) as with a normal stiffness k n , a normal viscosity n and the relative velocity in normal direction v n [31].
On top of this simple elastic contact force law, we separately add two types of adhesive forces: reversible adhesive force and jump-in liquid bridge capillary force. The details of these two adhesive forces are elaborated in the following.

Reversible adhesive force
For the dry granular material, we assume a (linear) van-der-Waals type long-range adhesive force, f a , as shown in Fig. 3a. The adhesive force law between two primary particles i and j can be written as,  with the range of interaction a ∶= f max a ∕k a c , where k a c is the adhesive strength of the material and f max a is the (constant) adhesive force magnitude, active for the overlap > 0 in addition to the contact force.
The adhesive force f a is active when particle overlap is greater than − a , when it starts increase/decrease linearly with slope k a c independent of the two particles approaching and separating. In the current study, this contact model is applied in the case of homogeneous stress controlled simple shear box simulations.

Jump-in liquid bridge capillary force
The capillary attractive force between two particles in a wet granular bulk is modeled with a discontinuous, irreversible attractive law as shown in Fig. 3b. The jump-in capillary force can be simply written as: is the normalized separation distance, S = − being the separation distance, r the reduced radius and V b the volume of liquid bridges. The maximum capillary force at contact ( S = 0 ) is given by f max a = 2 r s cos , with the surface tension of the liquid, s , and the contact angle, [54].
There is no attractive force before the particles come first into contact; the adhesive force becomes active and suddenly drops to a negative value, −f max a , when = 0 (the liquid bridge is formed). Note that this behavior is defined here only during first approach of the particles. We assume the model to be irreversible: the forces will not follow the same path, i.e. the attractive force is active until the liquid bridge ruptures at = −S c . This attractive force is following Willet's capillary bridge model [78], as explained in [54]. Similar to the reversible model, we combine the attractive capillary adhesive model with the linear visco-elastic contact model defined in Sect. 2.2.1 and use this combined model to simulate wet granular materials under shear in the split bottom shear cell.

Friction
Introducing the additional cohesive force in the normal direction between the two particles also influences the frictional force in the tangential direction. The tangential force is calculated in a similar fashion like the normal force through a "spring-dashpot" model ( f t = k t t ) and coupled to the normal force through Coulomb's law, f t ≤ f C s ∶= s f rep , where for the sliding case one has dynamic friction with In this study, we use p = s = d . For a non-adhesive contact, f rep = f n > 0 , and the tangential force is active and calculated based on the history in tangential direction. For an adhesive contact, Coulomb's law has to be modified in so far that f rep is replaced by f n − f a . In the current model, the reference for a contact is no longer the zero force level, but it is the adhesive, attractive force level [31,32].

Time scales and dimensionless numbers
Dimensional analysis is often used to define the characteristic time scales for different physical phenomena that involved in the system. 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/particle density p . 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 n as given by [65]: In the special case of a linear contact model, this is invariant and thus represents a global time scale. Two other time scales are globally invariant, the cohesive time scale t c , i.e. the time required for a single particle to traverse a length scale of d p ∕2 under the action of an attractive 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 s . This leads to the cohesion time scale as a global parameter given by: with maximum adhesive force f max a between two particles. 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 shear rate field: the shear rate is high in the shear band and low outside, so that also this time scale varies between low and high, respectively.
The time for rearrangement of the particles under a certain pressure constraint is driven by the total 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).
Finally, one has to also look at the granular temperature time scale due to the velocity fluctuations ( 2 ) among all the particles which is of significance of determining how "hot" the system is [13,34,72,75]: 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 four time scales ( t p , t c , ṫ and t k ) to influence the rheology of our system. Note that among the four time scales discussed here, there are six possible dimensionless ratios of different time scales. We propose three of them that are sufficient to define the rheology that describes our results. We do not expect a strong influence of the timescale related to the granular temperature on such dense granular flow and thus the t T is not included in the dimensionless study. The first dimensionless number is the Inertial number, 2 which is the ratio between t p and ṫ, The second Softness characterizes how "soft" the system is and is the square of the ratio between contact collision timescale t k and pressure timescale t p , The last one studied in this paper which is important to characterize the cohesiveness in the whole system is Bond number. It is also the square of the ratio between pressure timescale t p and cohesion timescale t c , Interestingly, all three 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 [52,65].
Even though not used, for completeness, we also propose the dimensionless granular temperature, T * = (t p ∕t T ) 2 , which is much smaller than 1 in all data analyzed further on.

Fixed particle parameters
The fixed input parameters for the two different simulation set-ups are summarized in Table 1. In the liquid bridge capillary contact model, there are two extra input parameters compared to the dry adhesive model: contact angle and liquid bridge volume. These two parameters are kept constant for all the simulations. Apart from that, we use restitution coefficient, e = 0.8 to compute the normal viscosity in the dry adhesive model but directly input the normal viscosity in the liquid bridge capillary contact model, which gives the same damping in the two contact models.

Control parameters
Apart from the fixed particle parameters, there are also several control parameters which we vary to explore their influences. For instance, varying shear rate, pressure, gravity and cohesion, we explore our rheology towards different inertia, softness and cohesion regimes. The details of parameter ranges are summarized in Table 2.
Alternative Inertial numbers could be defined using the bulk density, so that I � = I √ =̇d p ∕ � p∕ p , or the strain rate ( ̇= 1∕2̇ ), will lead to small differences in fitted characteristic numbers. For sake of brevity, we use the definition in Eq. (10).
To study the effect of inertia and contact stiffness on the dry-non-cohesive materials rheology in the geometry SS, we vary the shear rate, ̇ , between 0.0001 and 1 s −1 , as well as the pressure, p, between 125 and 4000 Pa. Thus, we obtain Inertial numbers, I, between 0.001 and 1; Softness, p * , between 0.001 and 0.1.
For the dry-cohesive simulations in the same set-up, we keep the shear rate, ̇= 0.005 s −1 , pressure, p = 500 Pa, and vary only the maximum adhesion force, f max a , between 0 and 0.01 N. This leads to the range of Bond number, Bo, between 0 and 5.
For the dry non-cohesive simulations using the geometry SB, we vary the rotation rate, , between 0.063 to 4.712 rad/s (0.01 and 0.75 rotation/s), of the outer wall to vary the system inertia and therefore match the case of the simple shear box. And to change the pressure in this free surface system, we achieve low or high pressure by varying the gravity, g, between 1 and 50 m s −2 .
In the case of wet granular material in SB, we vary the intensity of the maximum capillary force f max Note that for all the simulations performed using geometry SB, the inter-particle friction is kept constant at p = 0.01 , while for the simulations performed in geometry SS, we vary p from 0 to 1. Therefore, for the comparison between the two geometries, we choose only the data with p = 0.01 , and in each simulation performed, the Coulomb friction static and dynamic are always kept the same and referred to by the inter-particle friction, p = s = d .

Rheology
Our first goal is to check the consistency of bulk densities from two different geometries: SS and SB when using the same non-cohesive contact model. And the second goal is to investigate how different cohesive contact models influence  In addition, we also check the validity of existing frictional rheology and consummate it towards cohesive materials. Finally, in order to achieve our second goal, we explain the interplay between particle friction and cohesion in Sect. 3.3.

Non-cohesive granular materials
For dry granular materials, which have rather large particle size and negligible attractive forces ( Bo ≈ 0 ). The rheology (i.e. the equations of state for volume fraction and macroscopic friction) depend mainly on the Inertial number, I, and the Softness, p * , which are the competition between the time scales t p ∕ṫ and (t k ∕t p ) 2 , respectively. The dependence of the macroscopic friction coefficient = ∕p on p * and I has been studied earlier [52,65] as summarized in Appendix 1. In order to complete the rheology for soft, compressible particles, a relation for the volume fraction, , as function of pressure and shear rate is still missing for cohesive frictional materials. In previous work [35], the following dependency was reported: with the critical or steady state volume fraction under shear, i.e., the limit for vanishing pressure and Inertial number, 0 | p =0.01 = 0.64 . Luding et al. [35] found that the typical shear rate for which shear dilation would turn to extreme fluidization is Î � = 0.85 (for more details on number, see Eq. (10) and Table 3), and the typical pressure level for which Softness leads to too huge densities is p * = 0.33 . Note that both correction functions are first order Taylor expansions with respect to full, higher order functions, i.e. they are valid only for sufficiently small arguments. In slow, quasi-static flows in the split bottom shear cell simulations, weak dilation is observed, i.e., no strong dependence of on the local shear rate [35]. On the other hand, large Inertial numbers fully fluidize the system so that the rheology should be that of a granular fluid, for which kinetic theory applies. Small pressure, i.e. small overlaps have little effect, while too large pressure would lead to enormous overlaps, for which the contact model and the particle simulation with pair forces become questionable. Therefore, we focus on the linear expansion as shown in Eq. (13), valid in the quasi-static to (13) (I, p * ) = 0 g I (I)g p (p * ) ≈ 0 (1 − I∕I )(1 + p * ∕p * ) moderate inertial regime ( I < 0.25 ≈ I ∕10 ) and low to moderate overlaps between particles ( p * < 0.03 ≈ p * ∕10 ) with only few data outside this regime [35,59].
Note that in Eq. (13), we do not include inter-particle friction p as functional variable although the volume fraction depends on p . We consider inter-particle friction as a material micro-parameter like polydispersity or restitution coefficient, which is different from the dimensionless state variables in our functions, e.g. shear rate and pressure. However, we did explore this dependency in detail in "Appendix 2".
In Fig. 4, we plot the volume fraction in steady state as a function of Inertial number for both the stress controlled simple shear box (SS) as well as the split bottom shear cell (SB) with the same friction coefficient p = 0.01 . For the case of SS, when we we keep Softness p * constant and vary only shear rate (black squares and dashed lines), the volume fractions decrease with Inertial number, as the increase of shear rate leads to higher collisional/dynamic pressure. Correspondingly in the case of SB, when we vary the rotation velocity (red triangles), the volume fraction follows the same trend and the data collapse well with the data from SS but with slightly more scatter due to the local, small volumes used for averaging. The red triangles are not following any single of the dashed lines since both I and p * are small but neither is constant. When we fix shear rate and vary only pressure (blue circles and solid lines), the larger pressures lead to considerable compaction. In order to focus on the dependency of on Softness p * , in Fig. 5, we plot the same data as in Fig. 4 against Softness. For the case of SS, when the shear rate ̇ is kept constant while Softness p * is varying (blue circles), we observe an increase of volume fraction with Softness. For the corresponding case in SB, when we keep the rotation velocity constant and vary gravity (brown crosses), the volume fractions follow the same trend as in the case of SS and increase with Softness in a linear relation as specified in [65]. Different gravity data from the split bottom shear cell collapse also well with the simple shear box data but with slightly more scatter. The increase in Softness p * can be due to the pressure increase in the system or the particles becoming softer [61,72]; in either case the particles will overlap more thus leading to an increase in the volume fraction.
In the same Figs. 4 and 5, we also plot the prediction of Eq. (13) fitted using data from SS and with the fitting parameters summarized in Table 3. The coefficients for p = 0, 0.01 and 0.5 are from the fitting whereas the others are from estimation of the curve fitting, which are further tested and confirmed in this paper. Since the two systems have similar rheological behaviour, one would expect the rheological model developed from SB to capture also the behaviour of SS and vice versa. The prediction of fitted function looks promising and it captures well both the dependencies on inertia and Softness as shown in dashed and solid lines, respectively. All the fitting deviations are within 5%, even for data with relatively high inertia ( I > 0.25 ). However, when the volume fraction gets too low (below 0.5), the predictions of our rheological model deviate from the simulation data because the system goes from the dense state towards a loose granular fluid, which our model does not take into account. Alternative rheological models for this transition region are reviewed in [72], while the appropriate model for the liquid/gas state is standard kinetic theory as studied in [1,5,12,18]. This indicates the limits of using the rheological model Eq. (13); it predicts well moderate to dense granular flows with finite contact stresses, but not dynamic, less dense granular gases. Note that here in Table 3, we have also included the I = 2.11 from [34,35,52], 3 which is very close to our I = 2.246 obtained from the SS setup.

Cohesive granular materials
For cohesive granular materials, attractive forces enhance the local compressive stresses acting on the particles, thus leading to an increase in the volume fraction. On the other hand, rough and frictional particles will display a stronger dilatancy, i.e, a reduced volume fraction under shear at steady state, compared to their initial (over-consolidated) volume fraction. We systematically vary the inter-particle friction ( p from 0 to 1) and then for each inter-particle friction value we vary cohesion (Bo from 0 to 5) to study the volume fraction variations in steady state. Then we introduce a generalized rheological model involving cohesion, as based on evidence from the simulations.
The Bond number quantifies the competition between cohesion and pressure. Low Bond numbers refer to weak cohesion, or t p ≪ t c , while high Bond numbers indicate the  Table 3 The fitting coefficients of the non-cohesive rheological model for in Eq. (13) for various p Note that in this study, we use shear rate ̇ and particle density p in the Inertial number, whereas strainrate ̇ and bulk density b = p are used in other literature [34,35,52] for p = 0.01 . All the coefficients are fitted using data with > 0.5 , I < 0.25 and p * < 0.03 . Note the * symbol indicating the different definition as discussed in the main text and Eq. (10) 14 Page 10 of 20 opposite. The variation of Bond number in the current study covers a wide range of cohesion intensities by comparing the maximum pressure contributed by cohesion, f max a ∕d 2 p , to the actual pressure.

Effects of cohesion for different particle friction
In Fig. 6, we plot the steady state volume fraction against the Bond number for samples with different inter-particle friction using the SS geometry. As expected, for very low inter-particle friction ( p ⩽ 0.01 ), we observe a continuous increase in the volume fraction with Bo. However, for higher inter-particle friction ( p ⩾ 0.05 ), increasing dilatancy is observed with increasing Bo up to values around 2, an effect which we will discuss later. When the Bond number becomes larger, the volume fraction increases with increasing Bo. When we fix cohesion (Bo ≈ const.) and vary the particle friction, the volume fraction always decreases with increasing particle friction; this indicates that particle friction leads to shear dilation, not to compression. Thus, we assume that the change in volume fraction is due to both compression of soft particles in presence of cohesion and the dilatancy due to structural changes in presence of friction. We also added here the predictions of our proposed rheological models as in Eq. (16) and Eq. (1B) with solid and dashed lines, respectively. The predictions of Eq. (16) are in very good agreements with the actual data, but the lines of including Eq. (1B) over-predict the volume fractions at the steady state and deviate more towards higher Bond numbers. The details of how the rheological model is modified/generalized will be discussed in the following.
In cohesive flows, attractive forces enhance the compressive pressure acting on the particles. This can be quantified as follows: the net pressure can be split into two components, p = p rep − p coh , 4 denoting the respective contributions of repulsive and cohesive contact forces. The ratio between the total cohesive contribution and the total pressure is given by the local Bond number, Bo = p coh ∕p , and thus p rep = (1 + Bo)p , which represents the dimensionless overlap. As the geometrical compression (deformation at each contact) is related to the repulsive stress, it is the repulsive pressure p rep that has to be considered in the Softness factor g p to correct this cohesion induced bulk density change. However, this (1 + Bo) correction is not considered inside the Inertial number I, because the dominating timescale here is the pressure time scale t p that correlates to the total pressure p in the system not the repulsive part p rep . Thus, using Eq. (13), the modified Softness correction for cohesive systems is given as: For non-cohesive systems, Bo = 0 , g p is consistent with the non-cohesive rheological function. This modified pressure is similar in spirit to the modified Inertial number as presented in [4], which takes into account the cohesive contribution to stress. A similar modification of pressure in the Inertial number would require an additional parameter but would only reduce I, and thus is not considered here.  . 7 The scaled volume fraction, i.e. the cohesion term, g c = ∕( 0 g I g p ) as a function of the Bond number Bo for different inter-particle friction coefficients as shown in the legend. The lines are the predictions of the fitted Eq. (15) and the data are the same as in Fig. 6 In Fig. 7, we plot the same data as in Fig. 6 to obtain the functional form of the cohesive correction g c : the normalized volume fraction, g c = ∕[ 0 g I (I)g p (p * , Bo)] as function of the Bond number, Bo, in order to isolate the effect of cohesion on the dilatancy of soft particles. The normalized local volume fraction helps isolate the effect of cohesion only. From low to moderate Bo, the correction g c decreases with increase of Bo. An increase of the repulsive force between two contacting particles also leads to an increase of the sliding limit tangential forces, resulting in an enhancement of the role of friction. For Bo ⩾ 3 , g c increases from its minimum at Bo ≈ 3 , since the attractive force is so high that sample dilation is compensated by compression. Note that compression is prevailing for soft particles but is negligible when p * ≈ 0 , in the limits of low confining stress or infinite stiffness, when the local volume fraction is expected to monotonically decrease with Bo. Some previous studies have confirmed the negligible effect of confining stress p * when using stiff particles [4,21].
Here we assume that the frictional contributions are the same for both cohesive and non-cohesive materials. Therefore, we use the coefficients of the non-cohesive material: 0 , I and p * from Table 3. The local volume fraction is scaled by g I as in Eq. (13) and g p as modified in Eq. (14). In such a way, we hope to remove the effects of particle inertia as well as particle Softness.
Thus, cohesion contributes to the initial decrease and the subsequent increase in the volume fraction correction of g c for sheared material. In order to quantify this dependence, the correction function g c of the Bond number, Bo, is given by a fourth degree polynomial function: where, p 1 , p 2 , p 3 and p 4 depend on p . The actual fitted values for g c with different inter-particle friction, p , are summarized in Table 5. With increasing p , the polynomial constants p 1 , p 3 decrease, but p 2 , p 4 increase. The lines in Fig. 7 represent Eq. (15), where a promising collapse between this prediction and the simulation data is observed. For low particle friction, e.g. p = 0.01 , the volume fraction decreases gradually with Bo, while for higher particle friction, e.g. p = 0.5 , the volume fraction decreases more sharply with Bo. We note that also frictionless particles ( p = 0 ) under shear, the correction g c decreases with increasing Bo, which suggests that cohesion causes the structural changes that relate to shear dilation. In other words, on top of the effect in g p , cohesion could make frictionless particles stick together and form clusters, which results in an increase of the bulk shear resistance (data not shown here). The sample has to dilate to reduce its shear resistance to compensate this effect. Thus the role of contact friction in enhancing dilation is not a solitary effect, also cohesion is affecting the system behaviour: one should not neglect the role of either. Frictional particles need at least four mechanical contacts to form a stable structure while frictionless particles need six contacts. The higher the friction, the higher the probability that a stable packing can be formed at lower volume fractions with less average contacts/coordination number.
Comparing to the previous works [52,65], a more complete rheology including the cohesion/Bond number influence on the volume fraction is introduced as: In case of rigid particles, p * → 0 , Eq. (16) reduces to stiff = (I, 0, Bo) = 0 g I g c . The details of each contribution are summarized in Table 4.

Generalized rheology
Now we go further and include data from SB with a different contact model to compare with the data from SS. As discussed earlier in Sect. 2.2, we use two different contact models for the attractive forces in the two systems, SS and SB. Thus, the question arises whether the two models are influencing both systems in the same way. Since we discuss the granular rheology for soft particles, our second focus would be the effects of cohesion (Bond number) on the volume fraction, which are the increase due to the (15) g c (Bo) = 1 + p 1 Bo + p 2 Bo 2 + p 3 Bo 3 + p 4 Bo 4 (16) (I, p * , Bo) = 0 g I (I)g p ((1 + Bo)p * )g c (Bo)  compression and decrease due to the structural change enhanced by friction. First, we illustrate our data in a slightly different way in Fig. 8a, plotting the scaled volume fraction ∕[ 0 g I (I)g c (Bo)] as a function of the repulsive pressure (1 + Bo)p * . Unlike the previous section, here we use the function g c from Eq. (15) for the scaling. In such a way, dilation due to an increasing Bond number and inertial effects due to increasing shear rate are removed. For the case of SS, the scaled volume fraction increases with Softness, and the same trend is observed for the case of SB, regardless of non-cohesive or cohesive data.
The scaled local volume fractions ∕[ 0 g I (I)g c (Bo)] increase with the repulsive pressure (1 + Bo)p * , and all the data are found to collapse on a single increasing trend that is g p . Therefore we plot the predictions of our rheological functions as red and black solid lines in Fig. 8a. We see the two lines differ slightly from each other and we associate this to the small difference in the particle size distributions used in the two geometries, the simple shear box has a sample with polydispersity w = 3 , while the split bottom shear cell has w = 2 . Nevertheless, this difference is minimal and the global trend is well captured using the proposed rheological model. Note here, for the cohesive SB data, the Inertial number I ranges from 0.0004 to 0.0013 and the Softness p * ranges from 0.002 to 0.012, which differs from the cohesive SS data with constant I and p * . If one expects the correction function g c to be strongly affected by different I and p * , then these influences would show up in the data of SB where a wider range of Inertial numbers and Softness is covered. However, the influence/difference is not observed, therefore, the correction function g c if at all, is weakly affected by the dimensionless numbers I and p * , apart from the dominant dependence on the Bond number Bo.
When we look closer at the cohesive SS data as shown in Fig. 8b, we observe that the data for p > 0.01 deviate from the prediction for p = 0.01 (black line), which is captured using the proper p * ( p ) from Table 3 ( p = 1 is shown as dashed line). As we control all the samples having the same p * for the data shown here, the cause of the increase in volume fraction can only be from the increase of cohesion (Bo), which is the increase of the attractive forces among the soft particles leading to compression. The repulsive component of pressure is solely contributed by p * for non-cohesive materials. For cohesive material, the net repulsive contribution is due to the cohesive (tensile) as well as the pressure (compressive) forces, and this effect is found to be independent of inter-particle friction. For Bond number up to 2, all the data are collapsing on a single trend with less than 2% deviations relative to their absolute values, while the high Bond number data are deviating more from this trend, but still within 5% deviations up to Bond number 5. Here, the model predictions from this work are given as lines for two values of friction: black solid line for p = 0.01 and black dashed line for p = 1 . Both lines predict the corresponding dataset very well. The red line is slightly off due to a different polydispersity in [35].  Fig. 8 a The scaled volume fraction, i.e. the repulsive stress term, g p = ∕( 0 g I g c ) as a function of the cohesive softness p * (1 + Bo) and b zoom in to the simple shear data with different inter-particle friction coefficients and Bond numbers. Different symbols represent different particle friction, where small black diamonds represent local data from the split bottom shear cell (SB) for different Bond numbers using a wet cohesive capillary bridge force model (Fig. 3b) and p = 0.01 , while all other are from the homogeneous stress controlled simple shear box (SS). The red circles are the simple shear data using the normal visco-elastic contact model with no cohesive forces involved but varying shear rate and confining stress. The rest of the SS data are performed with different inter-particle friction coefficients and Bond numbers using the aforementioned linear reversible adhesive contact model (Fig. 3a). The black solid and dashed lines are the predictions of our rheological model Eq. (14) fitted using SS data with p = 0.01 and 1, respectively, while the red line is based on SB data with the same inter-particle friction as reported in [35] Page 13 of 20 14

The combined effect of inter-particle friction and cohesion
As mentioned earlier, cohesion among particles can lead to either compression or dilation of steady granular flow relative to its non-cohesive limit. Which mechanism dominates depends on the interplay between inter-particle friction and cohesion. In order to investigate how cohesion works with friction, we choose the SS set-up and systematically vary both inter-particle friction and cohesion to check their combined influences. For each inter-particle friction, we use the steady state volume fraction ( p , I, p * , Bo = 0) as reference and subtract it from the steady state volume fraction of finite Bond numbers, defining the cohesive effect on the steady state volume fraction: where positive indicates compression and negative indicates dilation of the sample relative to the cohesionless case.
In Fig. 9a, we plot the dependence of on both interparticle friction and Bond number. Using all points with colors indicating relative compression and dilation, we further interpolate , and obtain the phase diagram as shown in Fig. 9b. For weak friction and low to moderate cohesion, Bo < 0.1 , volume fraction changes are inconspicuous ( ≈ 0 ). For very high cohesion, compression is observed ( > 0 ), while relative dilation is observed ( < 0 ) at intermediate Bond numbers and strong enough inter-particle friction. However, when the role of friction is almost negligible ( p < 0.05 ), we observe a purely compressive effect with  Fig. 9 a The change of steady state volume fraction as shown in the color bar, plotted against Bond number Bo and inter-particle friction p ; b continuous phase diagram created from data in a with red and blue indicating relative compression and dilation, respectively. The Inertial number I = 0.02 and Softness p * = 0.01 are kept constant for all the simulations shown here (a) (b) Fig. 10 The differences of the steady state volume fractions plotted against a Bond number, Bo, for datasets with fixed inter-particle friction, p , and b inter-particle friction, p , for datasets with fixed Bond number, Bo, as shown in the legend. Lines are the predictions of our proposed model in Eq. (16) increasing Bo. This further confirms that sufficient interparticle friction leads to relative dilation of cohesive flow, where cohesion plays an important role in enhancing the dilation effect, but turning to compaction if large enough.
Although the phase diagram offers us a nice overview of the relative compression and dilation behaviour, we further study the data more quantitatively by fixing one parameter and looking at the influence of the second parameter. Therefore, in Fig. 10a and b, we plot the volume fraction for constant inter-particle friction and Bond number, respectively. When we fix the inter-particle friction and increase Bo, we observe negligible volume fraction change ( ≈ 0 ) up to Bo = 0.1 . If we increase the Bond number further, for low inter-particle friction ( p = 0.01 ), we see a monotonically increasing trend, while for moderate to high inter-particle friction ( p > 0.05 ), we get first a decrease then an increase, where the point of inflection depends on the inter-particle friction: higher inter-particle friction moves this point towards higher Bond number. This could be explained by the influence of friction: high inter-particle friction with moderate cohesion leads to strong relative dilation, but one needs even stronger cohesion to compensate the dilation effect and turn the bulk back into compaction. The predictions of the rheological model are plotted as lines. Instead of using the coefficients from p dependency fitting in "Appendix 2", here the model parameters (Tables 3 and 5) from the actual data fitting are used. For all the inter-particle frictions presented here, the predictions of the model are quite good up to Bond numbers around 1, then the model over-predicts the differences of the volume fractions, and the over-prediction increases with the inter-particle friction. Nevertheless, our model still predicts the volume fraction at steady state quite well, the highest deviation between the model prediction and the actual data point is around 4% of the actual volume fraction value measured from the simulation. Because we focus on a much smaller scale of the relative here (-0.06 to 0.08), the differences between the line and the actual data points only appear large.
When we fix the Bond number and increase the interparticle friction as shown in Fig. 10b, we observe the turning point more clearly. If the cohesion between particles is very small ( Bo = 0.1 ), the volume fraction is not changing much with inter-particle friction (relative to the Bo = 0 case). However, for stronger cohesion between particles ( Bo = 0.1, 0.5 and 1), we see some relative compression (positive values) for low inter-particle friction, which turns into strong relative dilation (negative values) with increasing inter-particle friction. Similarly, the predictions of the rheological model are also plotted here. As we pointed out earlier, for large Bond number, the prediction quality gets lower. Our model tends to over-predict the steady state volume fraction at higher Bond numbers, which results in larger . Note that here we focus on an even smaller variation of volume fraction ( −0.04 < < 0.01 ) and thus the deviations are acceptable.
Although we have shown our data in the above mentioned two possible ways, it is still interesting to look further into this inconspicuous volume fraction change (green) region in Fig. 9b and find out where the transition between relative compression and dilation happens. In Fig. 11, we narrow down this transition by applying a two-color map. Using the data shown in Fig. 11a, a sharp transition line between relative compression and dilation is reconstructed and highlighted in Fig. 11b. In the parameter range ( 0 ⩽ p ⩽ 1 and 0 ⩽ Bo ⩽ 5 ), three regimes can be identified: i) pure relative compression for p ≲ 0.027 ; ii) relative compression followed by relative dilation for increasing Bo (a) (b) Fig. 11 a The same data as in Fig. 9a plotted with a different colorbar and b the phase diagram with a sharp red/blue-relative compression/dilation transition. The vertical lines are the boundaries between zones 1-2 and 2-3, at p = 0.027 and 0.4, respectively and then relative compression again for 0.027 ≲ p ≲ 0.40 ; iii) pure relative dilation for p ≳ 0.40 (this value depends on the maximal Bo ≈ 5 chosen). In regime i), the friction effect is almost negligible such that cohesion dominates the flow behaviour leading only to relative compression. In regime ii), both friction and cohesion affect the flow behaviour: When cohesion is low, it has an almost negligible relative compression effect. Increasing cohesion, the relative dilation due to friction is enhanced by cohesion. But when cohesion is strong ( Bo > 1 ), the relative compression effect dominates again. In regime iii), the flow shows almost only relative dilation, mostly contributed by friction, up to Bo ⩽ 5 . If cohesion is extremely strong, ( Bo > 5 ), we expect that relative compression contributed by cohesion dominates the system again. However, this is beyond the scope of this study since our simulations are not homogeneous with such high cohesion.

Conclusion and outlook
We have extended an existing rheological model that predicts the relation between volume fraction, Inertial number and Softness [52], including friction and cohesion dependencies. We have calibrated this extended model without/ with cohesion using two different simulation geometries: a homogeneous stress controlled simple shear box and an inhomogeneous split bottom shear cell. Furthermore, two geometries are featured with two different cohesive contact models (dry, reversible vs. wet, irreversible). We systematically varied the inter-particle friction and cohesion in the stress controlled simple shear box with shear rate and normal stress fixed. We report an interesting interplay between the inter-particle friction and cohesion, as represented in phase diagrams. This allows the prediction of compression-dilation behaviour relative to non-cohesive reference systems in steady states.
Besides extending the rheological model towards cohesive-frictional granular media, we had two main goals: (1) to understand the relation between microscopic properties such as inter-particle friction or/and cohesion and macroscopic, bulk properties such as volume fraction; (2) to check the validity of our rheological model in both systems i.e., to confirm if the homogeneous representative elementary volumes (REVs) represent the center of a shear band in an inhomogeneous system (and maybe even the data away from the center). For completeness, the rheological model for the macroscopic friction is given in "Appendix 1".
Furthermore, we introduce a reversible van der Waals type cohesive force in the simple shear box and an irreversible liquid bridge type cohesive force in the split bottom shear cell. Independent of the type of cohesive model, if the two systems are in same steady state, e.g., same Inertial number, Softness, friction and Bond number, the steady state volume fractions from the two geometries agree well in the range of parameters studied. Our extended rheological model is thus applicable in different systems. The fact that different contact models can be unified by the Bond number indicates that the macroscopic rheological behaviour of the steady state volume fraction depends on cohesion intensity but not on the microscopic origin of cohesion between particles, whereas macroscopic friction dose, see Appendix 1.
Interestingly, when investigating the effects of cohesive models, we discovered that cohesion can either contribute to a decrease or an increase of the steady state volume fraction of sheared materials, relative to the non-cohesive reference case, enhancing or counter-acting the dilation effect of the inter-particle friction. Using the extended rheological model, we can successfully distinguish the two micro-mechanical mechanisms of cohesion: compression from the increase in the normal contact forces and relative dilation from enhancing both frictional forces as well as structural stability.
A phase diagram reveals how the combinations of these two particle parameters lead to sample compression or dilation in steady state shear, relative to a non-cohesive case. In addition, a sharp interface between compression and dilation on our phase diagram allows to categorize the explored parameter space into three regimes: i) pure relative compression for p ≲ 0.027 ; ii) non-monotonic behaviour with Bo: relative compression followed by relative dilation for higher Bo and then relative compression again, for 0.027 ≲ p ≲ 0.40 ; and iii) pure relative dilation for p ≳ 0.40.
The present paper is an extension of former works [52,65] on rheological modeling, but with a deeper insight on the influence of friction and cohesion. It could be enriched by exploring more closely how the micro-structure [64] is influenced by the combination of inter-particle friction and cohesion or vice-versa. Furthermore, extending the rheological model towards the intermediate to low volume fraction regime, where most dense rheological models fail but kinetic theory works well, is still a great challenge [72]. This will involve the granular temperature results in at least one more relevant dimensionless number in the rheology. Moreover, comparing the stress controlled system used here to a volume controlled system as in [72] is still ongoing research and will be addressed in the future.

Compliance with ethical standards
Conflict of Interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Macroscopic friction coefficient
For a full constitutive law of steady state rheology, one also needs to take into account the shear resistance as quantified by the macroscopic friction coefficient, = ∕ yy . This has been already developed in a previous work, in which the classical − I rheology on hard spheres was generalized for soft, cohesive granular flows [52], involving both dry [64] and wet [52,65] granular materials. The trends caused by different mechanisms are combined and shown to collectively contribute to the rheology [52] as multiplicative functions given by: where (I) = 0 f I (I) (see Table 6 or [52]) and all functions depend on one dimensionless number, while their coefficients depend on the particle friction coefficient, see Table 7.
The rheological function Eq. (1A) is slightly simpler than the general form in previous work [52], since we have taken out the contributions of f g and f q due to absence of free boundaries and too small Inertial numbers in the simple shear box. Remarkably, all the functions were (to leading order) linear in ratios of time scales and in particular, the macroscopic friction increases linearly to first order with the Bond number [52]. However, based on the dry, cohesive data (from SS) shown in the main text for , we propose here a new modification for the Bond number contribution as f Bo = 1 + 1 Bo 1 , where we claim that the increase of (1A) ∶= ( p , I, p * , Bo) = (I)f p (p * )f Bo (Bo) macroscopic friction is ruled by a power law. The details of each contribution are given in Tables 6 and 7 and then explained next. Note that here we limit our fitting range as following: > 0.5 , 0.005 < I < 0.5 and p * < 0.1 . The only different range used here compared with fitting is the limits of the inertial number I, due to the granular temperature and its consequent creep effect, as described in detail in [52].

Non-cohesive slightly frictional material
First we compare the data for non-cohesive slightly frictional material using the above mentioned two shear cell setups and show the results in Fig. 12. We see an inverse trend of macroscopic friction compared to the trend of volume fraction in Fig. 4. The macroscopic friction increases  [52] f c = 1 + aBo with Inertial number but decreases with Softness. The lines are the prediction of Eq. (1A) with Bond number equal to zero. The prediction is accurate when the Inertial number I is less than 0.5, but deviating for larger I. This can be explained by the considerable dilation of our constant stress shear box setup. When the system goes to very high Inertial numbers, the stress contribution from the kinetic part increases substantially and the granular bulk dilates in order to keep the stress constant; a much reduced bulk density is the consequence. Furthermore, we also compare the results between the two different setups: SS and SB, the center of shear band data of SB agrees very well with SS data, which confirms that the shear stress (macroscopic friction) responses are also the same when using the same material (same contact model) independent of the boundary conditions.

Cohesive materials
In Fig. 13, we plot the macroscopic friction, , against Bond number, Bo, with the predictions of our proposed rheological model. For the sake of clearness, we only plot data of three inter-particle friction coefficients instead of all the data we have. With increase of cohesion, the macroscopic friction stays almost constant for Bo < 1 , but increases when the cohesion gets stronger. High attractive forces between particles lead to a more stable packing, thus higher resistance to shear. We have also included the predictions from both our modified rheological model and the model from previous work [52], shown as solid and dashed lines, respectively. The latter can only capture the behaviour at low Bond number ( Bo < 1 ), while our modified model agrees very well with the SS data up to Bo < 5 . A possible reason that the previous linear model agrees well with the split bottom shear cell data is that they were fitting a data cloud, see Fig. 13c, which could introduce higher deviations. The data from simple shear box are cleaner and allow us to look at the influence of cohesion more closely, Fig. 13c nevertheless shows a significant qualitative difference between wet and dry cohesive models. All the model parameters are summarized in Table 7, the difference in the coefficients (related to cohesion) from the two setups are due to different contact models.
If we look at all the friction data (not shown), we observe that the macroscopic friction increases first with interparticle friction coefficient, but then decreases slightly for p > 0.3 . This non-monotonic trend might be caused by the friction induced, cohesion enhanced micro-structural anisotropy. The shear stress also shows similar trends as the macroscopic friction, whereas the pressure is changing continuously.
The granular temperature time scale, t T is comparable to the shear time scale ṫ but slightly smaller, while they  [52]. b The same data as in a but relative to 0 and in logarithmic scale. c The comparison between SS and SB data for p = 0.01 only are both much larger than the pressure time scale t p , which indicates that neither granular temperature nor shear are the dominating mechanisms for the friction trend, for the range of parameters studied here.

Appendix 2: The influence of inter-particle friction on the coefficients in the rheological model of volume fraction
For completeness, one has to also look at the dependency of the inter-particle friction p on the fitting coefficients obtained in our proposed rheological model as presented in Tables 3 and 5. The function used to describe this dependency is: where Q represents the coefficients obtained from the rheological model, Q 0 and Q ∞ represent the values of the coefficient at zero and infinite p , respectively. While 0 p refers to the typical inter-particle friction at which the change of coefficient starts saturating with the increase of p . For all the model coefficients, our proposed dependency laws fit well most of the data points with very small values of Root-Mean-Square (RMS) residuals, see the summarized details in Table 8. For the sake of brevity, we only plot here two coefficients out of seven due to the similar behaviour among all the coefficients. The results are shown in Fig. 14 with the effects of p on 0 and p 1 in (a) and (b), respectively. In the case of 0 , it decreases with characteristic 0 p around 0.28, then the decrease becomes weaker and saturates towards large p . The characteristic 0 p for the changes of I is smaller but larger for p * . For the cohesive coefficient p 1 , the decrease becomes weaker at p around 0.29, and this saturation turning points are very similar for all the other (1B) Q = Q ∞ + (Q 0 − Q ∞ )e (− p ∕ 0 p ) Table 7 The fitting coefficients of the rheological model Eq. (1A) for macroscopic friction. The numbers from [52] have additional correction terms for small I, but ignored for our fitting here. inf and I are not independent from each other and thus can vary strongly dependent on the available range of data All the coefficients are fitted using SS data with > 0.5 , 0.005 < I < 0.5 and p * < 0.1 p 0 inf (a) (b) Fig. 14 The effect of inter-particle friction, p , on a the fitting coefficient of the dry cohesive model, 0 , in Eq. (13); b the fitting coefficient from the cohesive model, p 1 , in Eq. (15). The lines are the predictions of Eq. (1B) and details are explained in the main text and summarized in Table 8  cohesive coefficients p 2 , p 3 and p 4 (see Table 8). The values of Q ∞ and Q 0 show a very wide variation for I , and rather small variation for 0 and p * . The cohesive coefficients might be approximated all by Q 0 = 0 and Q ∞ (p i ) ≃ 0.16(−1) i ∕i! , but our data quality is not good enough to support such a guessed series.