Saeed Khamseh
email: s.khamseh-1@utwente.nl

Jean-Noël Roux
email: jean-noel.roux@ifsttar.fr

François Chevoir
email: francois.chevoir@ifsttar.fr

Flow of wet granular materials: a numerical study

Keywords: numbers: 83.80.Fg, 47.57.Qk Keywords: granular flow, capillary cohesion, Mohr-Coulomb criterion Multiscale Mechanics (MSM), CTW, MESA+,

We simulate dense assemblies of frictional spherical grains in steady shear flow under controlled normal stress P in the presence of a small amount of an interstitial liquid, which gives rise to capillary menisci, assumed isolated (pendular regime), and to attractive forces, which are hysteretic: menisci form at contact, but do not break until grains are separated by a finite rupture distance. The system behavior depends on two dimensionless control parameters: inertial number I and reduced pressure P * = aP/(πΓ), comparing confining forces ∼ a 2 P to meniscus tensile strength F0 = πΓa, for grains of diameter a joined by menisci with surface tension Γ. We pay special attention to the quasi-static limit of slow flow and observe systematic, enduring strain localization in some of the cohesion-dominated (P * ∼ 0.1) systems. Homogeneous steady flows are characterized by the dependence of internal friction coefficient µ * and solid fraction Φ on I and P * . We also record normal stress differences, fairly small but not negligible, and increasing for decreasing P * . The system rheology is moderately sensitive to saturation within the pendular regime, but would be different in the absence of capillary hysteresis. Capillary forces have a significant effect on the macroscopic behavior of the system, up to P * values of several units, especially for longer force ranges associated with larger menisci. The concept of effective pressure may be used to predict an order of magnitude for the strong increase of µ * as P * decreases but such a crude approach is unable to account for the complex structural changes induced by capillary cohesion, with a significant decrease of Φ and different agglomeration states and anisotropic fabric. Likewise, the Mohr-Coulomb criterion for pressure-dependent critical states is, at best, an approximation valid within a restricted range of pressures, with P * ≥ 1. At small enough P * , large clusters of interacting grains form in slow flows, in which liquid bonds survive shear strains of several units. This affects the anisotropies associated to different interactions, and the shape of function µ * (I), which departs more slowly from its quasistatic limit than in cohesionless systems (possibly explaining the shear banding tendency).

I. INTRODUCTION

Over the last decade, constitutive modelling of dense granular flows has been proposed [1,2] in terms of internal friction laws directly applying to normal stresscontrolled steady shear flows, for which the internal state of the material is characterized by a single dimensionless number, the inertial parameter I [3]. Number I might be regarded as a reduced, dimensionless form of shear rate γ = ∂v 1 /∂x 2 , related to the stress σ 22 normal to flow direction x as I = γ m aσ22 , m denoting the particle mass and a its diameter. The constitutive law relating the effective internal friction coefficient, µ * , defined as a stress ratio, µ * = σ 12 /σ 22 , to inertial number I should be supplemented with a similar relation of solid fraction Φ to I [1,[4][5][6]. I characterizes dynamical effects, and the quasistatic limit is that of vanishing I. In this limit of I → 0, the material is in the so-called critical state of soil mechanics [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF], i.e., quasistatic plastic shear flow at constant solid fraction Φ c , under constant stresses and effective internal friction µ * c . In various experimental and numerical studies, the constitutive law, suitably generalized, was shown to apply to different grain shapes and flow geometries [START_REF] Jop | [END_REF][9][10]. On regarding inertial number I as the sole state parameter in a granular material in shear flow, it is implicitly assumed that small contact deflections due to the finite elastic stiffness of the grains are irrelevant -this is the rigid limit. A major advantage of the "µ * (I) and Φ(I)" approach is its ability to deal with both the quasistatic limit and the rigid limit without any divergence or singularity.

In the presence of attractive forces between neighboring grains, contacts are endowed with a finite tensile strength F 0 , whence a new dimensionless parameter, P * , a reduced pressure comparing the applied confining stress P (say, the controlled normal stress value σ 22 in shear flow) to force scale F 0 , as P * = a 2 σ 22 F 0 (similarly a "cohesion number" η = 1/P * was defined in Ref. [11]). Under small P * , cohesion stabilizes loose structures [12,13], which collapse upon increasing P * [14,15]. In steady shear flow, generalizing rheological laws to the cohesive case involves expressing internal friction coefficient and density as functions of both numbers I and P * or η [16].

In wet granular materials cohesion arises from capillary forces due to small liquid bridges joining particles touching or in close vicinity to each other [17,18]. The effect of such forces has been investigated in quasistatic deformation [19][20][21], and some of its consequences in terms of microstructure were discussed [22]. In the pendular regime of saturation [17,18] those bridges are small enough and do not merge, so that capillary forces are pairwise additive. Those attractive forces act as a source of cohesion, and are also characterized by a small range and some dependence on intergranular distance, as a liquid meniscus might join grains that are not in contact.

A traditional approach of partially saturated granular materials in geomechanics [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF], which has been investigated in recent DEM studies [START_REF] O'sullivan | Particulate discrete element modeling, a geomechanics perspective[END_REF][START_REF] Scholtès | [END_REF], is to resort to the concept of effective stresses, or stresses such that, if applied to the dry material, would produce the same deformation and plastic flow of the granular skeleton as the ones observed in the wet material. Proposed definitions of such an effective stress tensor in the unsaturated case generalize the Terzaghi principle [START_REF] Terzaghi | Soil mechanics in engineering practice[END_REF] applying to saturated media, and involve a correction of the average pressure related to saturation and capillary pressure [START_REF] Lu | [END_REF].

On the macroscopic scale, the effect of adhesive forces are sometimes described in the quasistatic limit of slow flow by the phenomenological Mohr-Coulomb law [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF][START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF][START_REF]Granular Media: Between Fluid and Solid[END_REF],

σ 12 = c + µ * 1 σ 22 , (1)
characterized by macroscopic cohesion c and internal friction coefficient µ * 1 . The present paper investigates the constitutive laws applying to wet model granular materials, in the pendular regime, and discusses the influence of capillary effects on macroscopic behavior and microstructure. Similarly to refs. [11,16], the rheology and micromechanical aspects are studied for varying P * and I values (with special emphasis on the quasistatic limit of I → 0). As in dry granular systems and in previous studies on 2D cohesive materials the material rheology is described in terms of apparent friction coefficient (stress ratio) µ * and solid fraction Φ as functions of I and P * , and the applicability of a Mohr-Coulomb relation is tested. Rheological and microstructural features, such as normal stress differences and formation of large clusters bonded by liquid bridges, are also investigated.

In the following, we first introduce (Sec. II) the microscopic ingredients of the model material, and report then, in Sec. III, on the conditions in which homogeneous steady states are observed in shear flows, enabling material constitutive laws to be deduced. Such laws are measured, depending on the relevant dimensionless parameters and on some features of the microscopic model, in Sec. IV. Next, in Sec. V, we investigate the role of capillary forces and distant interactions in the material rheology, and revisit the traditional concepts of effective stress and Mohr-Coulomb cohesion. Additional studies of microstructural and micromechanical aspects follow: force distributions (Sec. VI), agglomeration effects (Sec. VII), structural anisotropy (Sec. VIII). The results are discussed and put in perspective in the final, conclusive section IX.

II. MODEL MATERIAL AND SIMULATION SETUP

We consider a granular assembly composed of N equalsized spherical beads of diameter a, made of a material with Young modulus E and Poisson ratio ν. The contacts are frictional, satisfying Coulomb's law with friction coefficient µ. The granular flow is set by imposing a uniform shear rate γ to a rectangular parallelipipedic cell with edge lengths (L α ) 1≤α<3 . In order to avoid wall effects periodic boundary conditions are used in all three directions. The periodicity, in the direction of the flow gradient (direction 2), is applied with the Lees-Edwards procedure [START_REF] Allen | Computer simulations of liq-uids[END_REF] and in the two other directions the boundary condition is simple periodic. The system size L 2 is allowed to fluctuate in order to keep normal stress σ 22 constant, equal to a prescribed value P , while L 1 and L 2 are fixed [5].

A. Force model

Elastic and frictional forces are jointly implemented in contacts as in Ref. [START_REF] Agnolin | [END_REF], in which a simplified Hertz-Mindlin-Deresiewicz force model is used for the elastoplastic contact behaviour. This model combines the normal Hertz force F N , depending on contact deflection h as

F N = Ẽ√ a 3 h 3/2 , (2)
in which we introduced notation Ẽ = E/(1 -ν 2 ), with a tangential elastic force F T , to be evaluated incrementally in each time step of the simulation. The simplification of tangential elasticity adopted is that of Ref. [START_REF] Agnolin | [END_REF], involving a constant ratio (2-2ν)/(2-ν) of tangential (K T ) to normal (K N ) stiffnesses in contacts, both depending on F N , as, from (2), one has

K N = dF N dh ∝ F 1/3
. Caution should be exercised to avoid spurious creation of elastic energy with such laws, and K T should be suitably rescaled in cases of decreasing normal force and deflection. For the details of the elastic model, for the enforcement of the Coulomb condition ||F T || ≤ µF N , and for the objective implementation of the force law, with due account of all possible motions of a pair of contacting grains, the reader is referred to [START_REF] Agnolin | [END_REF].

An estimate of the typical contact deflection under confining stress P defines a dimensionless parameter, stiffness number κ [START_REF] Radjaï | Discrete-element modeling of granular materials[END_REF], such that h/a ∝ κ -1 . For a Hertzian contact, one may use [START_REF] Agnolin | [END_REF]

κ = ( Ẽ P ) 2/3 . ( 3
)
Two values of κ, 8400 and 39000, used in this study, respectively correspond to glass beads with E = 70 GPa and ν = 0.3 under pressures 100 kPa and 10 kPa, as in Ref. [5]. Finally, the force model of [START_REF] Agnolin | [END_REF] which is used here may also comprise a viscous damping term opposing normal relative motion of contacting grains, chosen to correspond to a restitution coefficient close to zero in normal collisions. We do not comment this feature, as it was shown [1,5] to have very little influence in the slow shear flows of the present study.

The presence of a small amount of an interstitial wetting liquid introduces additional capillary forces, transmitted between contacting or neighboring grains by a liquid bridge, or meniscus, as sketched in Fig. 1. We consider a perfectly wetting liquid, with contact angle θ equal to zero. As in Ref. [17], we assume that the menisci only form for touching particles, but break for gaps larger than a certain rupture distance D 0 , as observed in [START_REF] Kohonen | [END_REF]. D 0 relates to meniscus volume V as D 0 V 1/3 [34][35][36][START_REF] Maeda | Proc. Natl. Acad. Sci. USA[END_REF].

For the attractive force between particles separated by distance h ≤ D 0 we adopt the Maugis approximation [START_REF] Maugis | [END_REF], which is appropriate for small enough meniscus volume, for its simplicity. The maximum attractive force (tensile strength) is reached for contacting particles, and equal, according to this model, to F 0 = πaΓ (Γ is the liquid surface tension) independently of the meniscus volume. The capillary force varies with the distance h between particle surfaces as

F Cap =            -F 0 h ≤ 0 -F 0 [1 - 1 1 + 2V πah 2 ] 0 < h ≤ D 0 0 h > D 0 (4)
h < 0 corresponds to an elastic deflection of the particles in contact. This formula is a simpler, analytical form of the toroidal approximation with the "gorge method" [34] for the capillary force in a meniscus, which describes the meniscus as limited by circular arcs in a plane containing the two sphere centers. An alternative form was given by Willett et al. [35], while Soulié et al. [20,39,40] proposed a parametrized numerical solution. Fig.

B. Saturation range of pendular state

The morphology of partially saturated granular materials depends on the liquid content [18,41]. The present study, like a number of previous ones [19,21,40], is restricted to the pendular state of low saturations, in which the wetting liquid is confined in bonds or menisci joining contacting grains. Liquid saturation S is defined as the ratio of liquid volume Ω l to interstitial volume Ω v . It is related to meniscus volume V , solid fraction Φ = 1 -Ω v /Ω and wet coordination number z (the average number of liquid bonds on one grain) as:

S = Ω l Ω v = 3z π Φ 1 -Φ V a 3 . (5)
In our study, we fix the value of meniscus volume V . Such a choice does not conserve the total liquid volume, which is proportional to the varying coordination number z of liquid bonds. Its consequences have to be assessed, and we shall check that the results are not significantly affected, within the range of investigated material states. The pendular state to which our model applies is only relevant in some limited saturation range. On the one hand, a minimum liquid volume is necessary for menisci to form at contacts, as the liquid will first cover the grain surface asperities. This minimum saturation S min for bridges to form might be roughly estimated upon introducing a roughness scale δ, assuming a layer of thickness δ covers the surface of the grains, as

S min = 6Φδ (1 -Φ)a . ( 6
)
For Φ = 0.5 and δ ∼ 10 -4 a the minimum value for saturation is slightly below 10 -3 , comparable to experimental observations [17]. Using (5), and typical values for z (5 or 6) and Φ (0.5 or 0.6), this sets a lower bound to meniscus volume, of order 10 -4 a 3 . On the other hand, the upper saturation limit for the pendular state corresponds to the merging of the menisci pertaining to the same grain, which, considering a triangle of spherical grains in mutual contact, happens as soon as filling angle ϕ (see Fig. 1) reaches π/6. The analytical formula for V [34], within the toroidal approximation, as a function of ϕ (setting h = 0, and θ = 0), then yields V a 3 8.10 -3 , corresponding, using (5), to a maximum saturation between 0.05 and 0.1, similar to experimental observations [17,18].

C. Choice of parameters

Tab. I gives the values of parameters employed in our simulations. While stiffness number κ and friction coefficient µ are kept fixed, reduced pressure varies from the dry case P * = ∞ down to the lowest value 0.1, for which cohesive effects are strong, while the investigated range of I values allows us to approach the quasistatic limit with some accuracy, as well as assess the effects of inertia in faster flows (although rapid, strongly agitated flow are not studied here). The meniscus volume is chosen as V = 10 -3 a 3 in most simulations. A few tests are carried out with different values (as given within brackets) of the number of particles, the stiffness number and the preset meniscus volume. Taking Γ = 7.3 × 10 -2 J.m -2 κ 8400 (occasionally 39000) µ 0.3 N 4000 (8000) I from 10 -4 to 0.562 by factors of √ 10 P * 0.1 ; 0.436 ; 1 ; 2 ; 5 ; 10 ; ∞ V /a 3 10 -3 (10 -2 ; 5 × 10 -3 ; 2 × 10 -4 ; 10 -6 ) TABLE I: List of parameter values: N particles of diameter a, interacting with friction coefficient µ, forming menisci of volume V at contacts, are subjected to normal stress-controlled shear flow for which inertial number I, reduced pressure P * (evaluated with normal stress σ 22 ) and stiffness parameter κ take values as prescribed. Attractive forces fall to zero at distance

D0 = V 1/3 .
for water, and a = 0.1 mm, P * = 1 corresponds to σ 22 = (πΓ/a) 2.3 kPa -the pressure, under gravity, below a granular layer with a thickness of a few tens of centimeters. In our study, stiffness number κ is chosen so as to approach the limit of rigid grains (see, e.g., Ref. [42] for the effect of softer contacts).

III. HOMOGENEITY AND STATIONARITY

A. Steady states, macroscopic measurements

Starting from a dense initial configuration, with solid fraction close to the random close packing value (Φ RCP 0.64), we impose a constant shear rate γ and wait until a steady state is reached before measuring constitutive relations for stresses and solid fraction, which are identified as averages over time series. Stresses are measured using the standard formula, for all coordinate index pairs α, β

σ αβ = 1 Ω   i mv α i v β i + i<j F α ij r β ij   , (7)
involving a kinetic contribution with a sum over grains i, of velocities v i and a sum over all pairs with center-to center vector r ij , interacting with force F ij , Ω denoting the sample volume.

The evolution of solid fraction Φ with strain γ is shown in Fig. 3: Φ decreases until it approaches its steady state value for γ ≥ 5 in this case. Fig. 4 shows the evolutions of σ 22 and σ 12 with γ. (Note that σ 12 is negative with our sign convention). We thus check that normal stress σ 22 is well controlled since its fluctuations about its prescribed value P are very small. Shear stress σ 12 exhibits a fast increase and an overshoot at small strain and then decreases, approaching its steady state value, after a few strain units, over a strain interval (a few units) similar to the one corresponding to the transient evolution of Φ. Stresses and solid fraction fluctuate in the steady state, and a careful evaluation of measurement errors on their time averages is required (especially for shear stress, for which fluctuation levels reaching about 20% of the mean value are apparent in the example of Fig. 4). We use the blocking technique of Ref. [43] to estimate error bars on averages over finite time series.

B. Shear localization

Velocity profile

Instantaneous velocity profiles v s (≡ v 1 (x 2 ) s ) are computed on averaging particle velocities along the mean flow direction within slices of thickness 0.01L 2 . We observe a strong shear strain localization for the smallest value of the reduced pressure, P * = 0.1, for both slow and fast flows. As represented in Fig. 5a, the velocity gradient, initially uniform, gradually concentrates within a shear band of thickness H 3a which may move vertically but persists for all values of strain γ > 250. Localization tendencies in slow flow of dry granular materials are sometimes reported [1,5,44], although, for uniform strain rates, usually not observed as an enduring, systematic phenomenon. In the present study persistent localization profiles are also detected in nearly quasistatic flows, with a shear band thickness between 5a and 10a. However, for the intermediate values of the inertial number (10 -2 I 10 -1 ) this effect diminishes and strongly localized profiles are less frequent. Measurements of constitutive laws are then limited to intervals of time for which the velocity profile is devoid of shear banding (a criterion is specified below). For all P * ≥ 0.436, localization is not frequently observed and remains temporary, even in the quasistatic limit. The velocity profiles for P * = 0.436 and I = 10 -3 , as represented in Fig. 5b, are nearly linear and on average the flow is homogeneous.

Local solid fraction

Similarly we record solid fraction profiles, on averaging the solid contents of slices orthogonal to the velocity gradient (splitting the volume of one grain between different slices if necessary). Fig. 6 shows the velocity (v s ) and solid fraction (Φ s ) profiles for two different values of shear strain, γ = 1 and γ = 352, which belong to the simulation of Fig. 5a. We see that in the homogeneous flow the distribution of mass in the system is almost uniform, but when the localization occurs Φ s strongly decreases within the shear bands, to a value below 0.2. It slightly increases outside the shear band, especially in its vicinity. A slighter decrease of density within thicker shear bands in quasistatic flow is observed. For instance, when I = 10 -3 , Φ s decreases from 0.47 to about 0.4 inside the shear band of thickness H ≈ 7a.

Deviation from linear profile

The deviation from the linear profile is characterized by parameter ∆:

∆(t) = 12 L 3 2 γ2 L 2 /2 -L 2 /2 (v 1 (x 2 ) -γx 2 ) 2 dx 2 . ( 8
)
The normalization by L 2 as if two solid blocks were sliding on each other. As defined in Eq. ( 8), parameter ∆ is not affected by a vertical shift in the velocity profile, due to the periodic boundary condition in direction x 2 . If the strain rate is homogeneous within a shear band of thickness H and vanishes outside, one observes

X 2 v 1 Φ s (b) -10 -8 -6 -4 -
∆ = (1 - H L 2 ) 2 (9)

Occurence of shear banding

Large values ∆ > 0.8 for P * = 0.1 in the faster flows (I ≥ 0.17) indicate strong localization in this range. At I = 0.1 ∆ drops down to small values, typically below 0.1, but in the quasistatic limit it increases again: at I = 10 -3 it mainly fluctuates between 0.4 and 0.8. For P * ≥ 0.436 the shear rate is much more homogeneous. ∆ almost vanishes in faster flows, increases somewhat in the quasistatic limit, but rarely exceeds 0.2, even for the smallest inertial numbers.

Simulations carried out with a larger stiffness number (κ = 39000), for the two smallest values of P * and for all values of I in Tab. I, do not record any significant influence of κ on the homogeneity of the flow. The influence of the sample size is studied by simulating some samples with height L 2 twice as large as in the standard sample, containing 8000 grains, with P * = 0.436 and different values of I. The size dependence in formula 9 implies then larger values of ∆, should the shear strain tend to localize, temporarily or permanently, within a region of fixed thickness. In our tall, 8000 grain systems, as the quasistatic limit is approached, ∆ reaches peak values above 0.4 but continuously evolves and no persistent localization pattern is detected.

Consequently, our results reveal a strong localization tendency at P * = 0.1 for both small (below 0.03) and large (above 0.3) values of the inertial number. We performed some measurements at P * = 0.1 for intermediate values of I, over strain intervals for which values of inhomogeneity parameter ∆ averaged below 0.1, as in the first part of the graph of Fig. 7. The influence of cohesive forces on shear band widths was also reported (in a different context of 'forced localization') in Ref. [45].

A systematic fluid depletion in shear bands was reported in [46] -this requires a model for liquid migration between menisci, which we did not introduce in the present study. Leaving detailed investigations of flow localization phenomena for future works, we limited here our results to the issue of whether shear banding occurs for given P * and I values. The remainder of the paper deals with homogeneous flows, for the values of P * for which no evidence of enduring localization effects is observed.

IV. MACROSCOPIC BEHAVIOR AND CONSTITUTIVE RELATIONS

We now deduce macroscopic constitutive relations from the simulations of homogeneous flows.

A. Shear stress and solid fraction

Friction coefficient µ * and solid fraction Φ, depending on I for various P * values, are shown in Fig. 8 for the parameter choice adopted in most simulations. We fit the

* (D0 = 0.1a = V 1/3 ).
following power law functions to those data, denoting as µ * 0 and Φ 0 the quasistatic limits (critical state) values of the macroscopic friction coefficient and of the solid fraction:

µ * = µ * 0 + cI α Φ -1 = Φ -1 0 + eI ν (10)
Tabs. II and III give the values of the fitting parameters introduced in Eqs. 10. While the increase of µ * and the decrease of Φ as functions of I are familiar trends, similar to observations made with dry grains [1,[4][5][6], some other features are remarkable. The quasistatic limit is quite nearly approached for I ≤ 0.01, and is strongly influenced by capillary forces. Internal friction coefficient µ * , compared to the dry, cohesionless value (0.332±0.004), already shows a notable increase at P * = 10, reaching values as high as 0.6 for P * = 1 (i.e., as cohesive and confining forces are of the same order), and nearly 0.9 for P * = 0.436, about 2.3 times the cohesionless value.

Our partial results for P * = 0. Such a strong influence of cohesive (capillary) forces contrasts with the results of Refs. [11,16], in which similar deviations between cohesionless and cohesive systems are not observed until P * decreases to much lower values, of order 0.01. Such 2D results were however obtained with a different attractive force law, of vanishing range beyond contact.

B. Normal stress differences

The first and the second normal stress differences are defined as

N 1 = σ 11 -σ 22 N 2 = σ 22 -σ 33 (11)
Note that those definitions coincide with the one used in complex fluid or suspension rheology [START_REF] Guazzelli | A Physical Introduction to Suspension Dynamics[END_REF], but that we use the opposite sign convention for normal stresses. Signs of N 1 and N 2 should thus be reversed for comparisons to this literature. Most often, considering dense flows of dry granular materials, those differences, deemed small, are ignored or neglected [START_REF]Granular Media: Between Fluid and Solid[END_REF]. We find it worthwhile to record their values nevertheless, since, as shown in Fig. 9, where N 1 and N 2 are plotted versus I for different values of P * , they are strongly influenced by capillary forces. The first normal stress difference is very small in the quasistatic limit and for large values of the reduced pressure. It increases with I and for decreasing values of P * , going through a transition from small negative values to positive values, between I = 0.01 and I = 0.1, for P * ≥ 5. N 1 variations with I are nearly parallel for different P * values. The second normal stress difference N 2 also increases for faster flows and for decreasing reduced pressure P * . In the quasistatic limit, it is considerably larger then N 1 .

C. Sensitivity to capillary force model and saturation

Capillary force model

We tested the effect of the capillary force model by replacing the Maugis approximation, Eq. 4, with the more accurate parametrized capillary force law proposed by Soulié et al. [20,[START_REF] Gras | Youssoufi[END_REF], for V = 10 -3 a 3 . Although the difference in the force models is appreciable on a plot of F Cap versus h (with the Soulié force about 10% smaller at contact, see Fig. 2), the macroscopic results are very close: the difference in stress ratio µ * and solid fraction Φ increases with I but does not exceed 2%.

Meniscus volume and force range

Changing the meniscus volume amounts to changing the distance at which the attractive force vanishes, rupture distance D 0 = V 1/3 , as well as the gap dependence of the capillary force F Cap (h) (Fig. 2). Fig. 10 shows internal friction coefficient µ * (I) to be sensitive to meniscus volume for the lowest P * values. For a meniscus volume of 10 -6 a 3 , as compared to the standard value 10 -3 a 3 , µ * decreases by about 20%. Actually, for such a small meniscus volume, the decay of the attractive force (Fig. 2) is so fast that, as we checked, results are hardly changed on setting the meniscus rupture distance to zero. The effect on the solid fraction Φ remains small. To explore the rheological properties throughout the pendular regime, we varied the meniscus volume, and recorded the solid fraction and the friction coefficient in the quasistatic limit for the smallest studied P * value, as indicated in Tab. IV, thus fully covering the corresponding saturation range (see Sec. II B). Saturation S, by relation (5), is related to the wet coordination number, z, whose values are also provided in the table. While the change in solid fraction does not exceed 0.01, the variation of the macroscopic friction coefficient is about 20% in the pendular regime (up to 50% upon extending the numerical study to unrealistically small menisci, V = 10 -6 a 3 ). We therefore predict a moderate variation of rheological properties within the simulated pendular regime of the partially saturated granular assembly. Returning to the basic assumptions of our model, one of its drawbacks is that it ignores liquid volume conservation. Within the granular sample, the total liquid volume is proportional to coordination number z. As z varies with I, we should in principle correct the meniscus volume to maintain a constant product zV for different shear rates. However, the V dependence of macroscopic properties is so slow (µ * varies by 20% as V is multiplied by 20) that the resulting correction on V (as z changes, typically, from 6 to 4 at most) should be hardly notable. We explicitly checked it on running two series of simulations, for P * = 0.436 and for P * = 10, starting from the standard value V = 10 -3 a 3 for I = 0.001, and correcting V so that the total liquid content was kept constant at larger inertial numbers. This required slight increases of V , only exceeding 1% for I or order 0.1. Most quantities were indistinguishable between the unchanged and the corrected V values. In particular, all flow characteristics in the quasistatic limit were quantitatively unaltered. Only µ * varied by about 3% at most for I ≥ 0.1.

Hydraulic hysteresis

Another feature of the meniscus model the role of which should be explored is the hysteresis of the attractive force, which appears at contact, and vanishes at distance D 0 . As the number of interacting grains increases with the force range (see the values of z in Tab. IV), one may expect a strongly enhanced influence of distant interactions (as reported, e.g., in Ref. [49]) if menisci are assumed to form as soon as non-contacting grains approach below distance D 0 . Fig. 11 compares internal friction and solid fraction, for different values of P * and I, in the standard, hysteretic model, and without the capillary force hysteresis, assuming a force range D 0 between approaching grains that were not previously in contact. Without hystreresis, Φ notably increases, especially for small values of I. The internal friction µ * , for I 0.1, is close to the standard case, but larger values are obtained as I decreases. Even for the smallest values of I investigated (I = 0.001), the material properties still depend on shear rate and no proper critical state appears to be approached in our simulations. The decrease of µ * as a function of I in interval 0.001 ≤ I ≤ 0.01 should trigger shear-banding instabilities, as discussed in [50,51]. A slightly decreasing trend of µ * versus I was also apparent in Fig. 10, for very small D 0 . For the standard value V = 10 -3 a 3 adopted in this study (as one for which laboratory observations should be possible), the friction coefficient does increase with I, albeit slower and slower as P * decreases (see coefficient c in Tab. II). The stabilizing effect of this growing variation is weaker as cohesion gets stronger, consistently with the systematic shear banding behaviour at P * = 0.1, and might be jeopardized on tampering with the capillary force model.

V. RHEOLOGICAL EFFECT OF CAPILLARY FORCES

We now seek to explain the strong influence of capillary forces on the macroscopic material rheology. The roles of different interactions, in the force network and in the stresses are investigated. We first collect information on coordination numbers and neighbor distances (Sec.V A). Simple relations to average forces are recalled in Sec. V B. We split the stresses into several contributions, in order to appreciate the importance of different types of forces. This decomposition (Sec. V C) suggests an attempt to relate the rheology of wet grains to that of dry ones, in terms of some "effective pressure" approach, in the quasistatic limit, which we present in Sec. V D.

A. Coordination numbers and near neighbor distances

Fig. 12 shows the I dependence, for different P * values, of coordination numbers z c , for pairs of grains in contact, and z d , for pairs of grains attracting each other without contact at a distance lower than D 0 . The average number of contacts per grain, z c , decreases for larger inertial numbers, as previously observed in cohesionless systems [1,5] and in cohesive ones [16], slower for smaller P * , as in [16] too. z c also increases as P * decreases at constant I, as previously observed as well [16]. Note that this latter trend is opposite to that of the solid fraction (Fig. 8): as the importance of adhesion, relative to confinement stresses, increases, looser systems are obtained, yet better coordinated. Grains tend to stick to one another, and may form loose aggregates, as in static or quasistatically compressed assemblies, for which little correlation is also observed [13,15] between density and coordination number. On the other hand, the variations of the coordination number of distant interactions, z d , with both parameters I and P * , are in the opposite direction to those of z c . As I increases, so does z d : contacting pairs tend to separate, but some remain bonded by liquid bridges. And for stronger cohesion (smaller P * ), z d is correlated with the system density. The faster approach to quasistatic limit at smaller values of P * is apparent in both figures. The fraction of rattlers (beads carrying no force [5]) in non-cohesive systems is about 5%. In the cohesive case, due to the attractive forces, nearly all of the particles are bonded to others and the number of rattlers tends to zero, as observed in 2D simulation of cohesive powders [13,15]. z d tends to compensate the changes of z c , so that the total coordination number z = z c + z d , throughout the investigated range of I and P * values, exhibits rather small variations (see Fig. 13). Within the investigated parameter range, the maximum change in z, between 6.8 and 4.8, corresponds to a correction of internal friction µ * , should we change the meniscus volume to maintain the total liquid volume constant, below 5% (see discussion in Sec. IV C 2). The contact coordination number does not change much with the force range or the meniscus volume. Setting D 0 = 0 (instead of the standard value 0.1 used in the present study) or decreasing the volume of the meniscus from its standard value V = 10 -3 a 3 down to 10 -6 a 3 , merely leads to a small decrease of z c , from 5 to 4.7 in the quasistatic limit, when P * = 0.436. However, it has a strong influence on z d . Compared to the standard case, for P * = 0.436 and small values of I, it decreases from 0.9 down to a value below 0.3 when we set D 0 = 0.01a (still with V = 10 -3 a 3 ), or down to about zero when we set V = 10 -6 a 3 (with

D 0 = V 1/3 ).
It is interesting to compare the number of distant, interacting pairs to the total number of neighbor pairs at distance below D 0 . The coordination number, z(h), of neighbor grains at distance below h (such that z(0) = z c ) grows with h as depicted in Fig. 14, corresponding to I = 10 -3 (quasistatic limit). z(h), like the contact coordination number, is a decreasing function of P * for small h/a (below about 2.5 × 10 -3 , see the inset in Fig. 14). It increases with P * , like the density, beyond that distance. In denser systems grains have more neighbors on average, but this is only true if neighbors at some distance are included in the count, and does not apply to contacts (a situation reminiscent of some observations in static packings of cohesionless grains [START_REF] Agnolin | [END_REF]). Up to meniscus rupture distance D 0 , equal to 0.1a in the present case, each grain has on average z(D 0 )-z c non-contacting neighbors, among which z d are joined by a liquid bridge. Values of ratio z d /(z(D 0 ) -z c ) for different P * and I are given in Tab. V. The proportion of the neighbors within range D 0 that are bonded by a liquid bridge varies from 0.61 to 0.71 for I = 10 -3 and between 0.68 and 0.79 for I = 10 -1 -slightly larger than the proportion ∼ 50% reported by Kohonen et al. [START_REF] Kohonen | [END_REF] in static grain packs. If the meniscus forms as soon as grains approach at distance D 0 , rather than at contact, the number of contacts hardly changes (z c increases by about 5% for P * = 0.436 in the quasistatic limit), but the increase in the number of menisci is larger than expected from the data of Tab. V, from a simple count of pairs within range D 0 : z d is multiplied by 1.7 at small P * and I.

I = 10 -3 I = 10 -1 P * z d z d /(z(D0) -zc) z d z d /(z(

B. Pressure and average normal forces

From Eq. 7, neglecting the deflection of contacts in comparison to grain diameter a, and ignoring the kinetic term, one may relate [START_REF] Agnolin | [END_REF] the average pressure, P = trσ/3, to the average normal force F N for all interactions, and to the average, F N h d , over pairs in distant interactions, of the product of force by distance h ≤ D 0 :

P = Φz πa 2 F N + Φz d πa 3 F N h d (12)
Due to normal stress differences, ratio P σ 22 is only slightly different from 1 (about 0.95) at small I. We checked that formula ( 12) is very accurate for all P * values, and found its second term to be negligible, contributing less than 2% of the pressure.

A Bond number is defined in Ref. [45] as Bo = F 0 F N to compare cohesion and confinement forces. From (12), it is roughly proportional to 1/P * .

C. Contributions to stresses

The contribution of the kinetic term in Eq. 7 to stresses is quite small. Even for the fastest flow in our simulation (I = 0.562), this contribution does not exceed 2% of the shear stress or 5% of the normal stress components, and for I = 0.178 it is nearly zero for all stress components. Therefore, in this section we only discuss the contributions of forces to the stress components, for the different values of the control parameters, P * and I. These contributions may be split in different ways, on distinguishing different forces.

Contact forces and distant capillary attraction

First one may consider the total stress as a sum of the contributions of the contacts and of the distant interacting pairs, as

σ αβ = σ c αβ + σ d αβ . (13)
Our results show that the contribution of contact forces dominate in the shear stress. It is larger than 90%, regardless of the values of P * and I. The contribution of distant interactions to σ 12 , as represented in Fig. 15 a, although not negligible, hardly reaches 10% of the total shear stress for the smallest values of P * . The contribution of distant interactions to σ 22 is displayed in Fig. 15b. Capillary forces being attractive, σ d 22 is a tensile stress. For P * = 0.436, in the quasistatic limit, this contribution increases up to 20% in magnitude. Consequently, the positive contribution of contact forces to σ 22 reaches about 1.2σ 22 for P * = 0.436.

The relative importance of the contributions of contacts and distant capillary forces to σ 11 and σ 33 is similar: in the quasistatic limit and for P * = 0.436, one has σ d 11 /σ 11 -0.16 and σ d 33 /σ 33 -0.25. Accordingly, the contribution of normal contact forces is the dominant one in normal stress differences N 1 , N 2 (with a notable contribution of tangential forces to N 2 , typically 20% at low P * ).

Elastic-frictional forces and capillary forces

An alternative decomposition of the stress tensor is:

σ αβ = σ cap αβ + σ Ne αβ + σ T αβ , (14)
in which σ cap αβ is the contribution of capillary forces (either in the contacts or for distant interacting pairs), σ Ne αβ is the contribution of normal elastic forces and σ T αβ is the contribution of tangential forces.

The normal elastic forces contribute more than 90% of the shear stress, whatever P * and I.

The contribution of tangential forces to the normal (diagonal) elements of the stress tensor is negligible, but that of capillary forces is very important: for P * = 0.436, negative terms σ cap αα (1 ≤ α ≤ 3) are very large in magnitude: one observes σ cap αα < -2σ αα for small P * , as shown in Fig. 16. This large negative contribution is compensated by that of the repulsive normal elastic forces, σ Ne αα > 3σ αα . Such a large negative contribution of capillary forces to pressure implies that the particles are strongly pushed against one another, which increases the sliding threshold for tangential contact forces. Fig. 17 shows the contribution of tangential forces to the total shear stress. As P * is decreased to P * = 0.436, the ratio σ T 12 /σ 12 increases to 0.18. Capillary forces contribute to the shear stress with the opposite sign (σ cap 12 is positive while σ 12 is negative). Fig. 17 shows that ratio σ cap 12 /σ 12 is always negative and decreases down to -0.12 for P * = 0.436. Similarly to the case of normal stresses, the largest contribution is that of elastic normal forces, the (negative) contribution of which to σ 12 compensates the (positive) term σ cap 12 .

D. Discussion

One important clue to understand the enhanced shear strength of the cohesive material, as compared to the cohesionless, dry granular assembly, is the large tensile contribution of capillary force to normal stress:

σ cap 22 = -βσ 22 , (15)
with a coefficient β ranging, in the quasistatic limit, from about 0.15 (P * = 10) to 2.1 (P * = 0.436). Upon including the result for P * = 0.1 and I = 0.01 (the intermediate value of inertial number, assumed close to the quasistatic limit, for which measurements are possible in homogeneously sheared systems), β reaches about 7.2. This coefficient, and its variations with P * , can be approximately predicted from the values of solid fraction and coordination numbers. Contacts (z c , on average, per grain) carry capillary force -F 0 , and distant forces (z d per grain) average to a fraction of -F 0 . Relation ( 12) can be used to evaluate the capillary contribution to pressure P, as -ΦzF 0 πa 2 ≤ P cap ≤ -Φz c F 0 πa 2 . (This relation between P cap and contact tensile strength F 0 is sometimes referred to as the Rumpf formula, especially in the context of a prediction of rupture conditions [13,19,52]). Dividing by σ 22 , one obtains:

- Φz πP * ≤ P cap σ 22 ≤ - Φz c πP * . ( 16
)
Ignoring the small difference between P and σ 22 , ( 16) provides an estimate of coefficient β defined in (15). Thus the value of β for reduced pressure P * = 0.436 is predicted between 1.9 and 2.3 (and for P * = 0.1, it should reach about 8). Thus, quite unsurprisingly, the (negative) relative contribution of capillary forces to normal stress if of order (1/P * ) ∝ F 0 /P , with a coefficient that may be deduced from Φ and coordination numbers, according to (16).

It is tempting to invoke a classical concept in geomechanics, that of effective pressure, to describe the effect of capillary forces on the shear resistance of the material: the attractive forces create larger repulsive elastic reactions in the contact, corresponding to an effective pressure equal to (1 + β)P. Furthermore, the local Coulomb condition in the contacts is to be written with those enhanced normal repulsive forces. Capillary forces also contribute to shear stress, but, as apparent in Fig. 17, in comparison to their influence on normal stresses, this is a small effect, and one may ignore it in a first approach. One assumes then that the shear behavior of the material is identical to that of a dry material under such effective normal stress σ eff 22 . This approach leads to the following prediction for the P * -dependent quasistatic friction coefficient µ * 0 :

µ * 0 = (1 + β)µ ∞ 0 , (17)
in which µ ∞ 0 denotes the quasistatic internal friction coefficient for dry grains, P * = ∞. Remarkably, if we further assume, as suggested by (16), that β is roughly proportional to 1/P * , β b/P * , we obtain a Mohr-Coulomb relation, Eq. 1, for the stresses in the critical state: with the same value of internal friction as in the dry case, µ * 1 = µ ∞ 0 and a macroscopic cohesion given by

c = bµ ∞ 0 F 0 a 2 . ( 18
)
Fig. 18 is a plot of σ 12 versus σ 22 -the yield locus -in which the predictions of relation 17, both with the measured coefficient β (Fig. 16), and with the one predicted as (z + z c )Φ/(2πP * ) from ( 16), are confronted with the numerical results. The admittedly crude prediction of relation (17) appears surprisingly close to the numerical results on this plot. The relative error in the prediction for stress ratio µ * 0 , with the measured value of β, is actually about 5% at P * = 10, increasing to 20% at P * = 0.436, and the value of µ * 0 for P * = 0.1 ( 1.6) from the measurements at I = 0.01 is largely overestimated, at 2.7.

One may directly test for the validity of a Mohr-Coulomb relation to the data by fitting a linear form for the data of Fig. 18. Given the error bars (which are small and do not appear on the graph), an attempted straight line fit through all 5 data points with P * ≥ 0.43 in Fig. 18 is unambiguously rejected by the standard likelihood criterion. A linear fit is (barely) acceptable upon ignoring the value P * = 0.436, yielding µ * 1 = 0.340 ± 0.001 and a 2 c/F 0 = 0.267 ± 0.005 for the Mohr-Coulomb parameters. From (17) the predicted apparent macroscopic cohesion is above 0.3F 0 /a 2 , and varies according to which data are used to identify b in (18). The result µ * 1.6 for P * = 0.1 (corresponding to a 2 σ 12 /F 0 = µ * P * = 0.16) is thus in contradiction with the Mohr-Coulomb model, which becomes increasingly inadequate for smaller P * , as apparent in the insert in Fig. 18.

The performance of the simple effective pressure prediction for the P * dependence of µ * 0 is better visualized in Fig. 19, which, unlike Fig. 18, is not sensitive to stress scale. The global increase of µ * 0 is predicted, yet overes-FIG. 19: (Color online) Apparent quasistatic friction coefficient µ * 0 versus 1/P * -showing the value of µ ∞ 0 for 1/P * = 0. Measurements and predictions of (17), with exact and estimated coefficient β, same symbols as in Fig. 18.

timated for the smallest P * values.

One aspect that is not captured by this approach is the dependence of µ * 0 on meniscus volume (Tab. IV): the variations of coefficient β (from 1.8 to about 2.2 as V increases from 10 -6 a 3 to 5.10 -3 a 3 ) are insufficient to account for the increase of the friction coefficient.

There are quite a few reasons for the effective pressure approach to fail: while the mechanical properties are supposed to be the same once stresses are corrected, the density of the material, for one thing, is different in the dry and the wet case (with Φ varying between 0.525 and 0.595 as P * grows from 0.46 to infinity); capillary forces also contribute to shear stress, the force network is bound to be different, etc. Nevertheless, although admittedly crude, the prediction based on (17) proves apt to capture the trend of the change of µ * 0 with P * , although it overestimates its growth at small P * . As to the Mohr-Coulomb representation of yield stresses, it might be used as an approximation for P * ≥ 1, but the observations clearly preclude the definition of unique values of macroscopic cohesion and friction coefficient according to (1) for smaller pressures.

Pierrat et al. [52] report on a laboratory study of quasistatic yield loci (σ 12 versus σ 22 at the onset of plastic yielding and flow) of various kinds of wet granular assemblies in the pendular regime, including glass beads, which offer a suitable experimental comparison to our results. It proposes (under the name 'shift theory'), exactly the same effective pressure approach as the one we have attempted here, and concludes that it provides a good approximation, by which the yield condition of wet materials is deduced from the one of the dry grains. Interestingly, the investigated P * values in this study range from about 0.2 to ∼ 2.5, and on the yield locus the increase of σ 12 with σ 22 is slightly sublinear, as in our numerical results. Measured values of µ * are similar to our results (with e.g., µ * 0.7 for P * = 1), and little change is obtained by increasing saturation by a factor of 3. Some possible differences between those experiments and our simulations could result from the different state of the material: the experiments are not necessarily carried out in steady state quasistatic shear flow, and could depend on the initial assembling process. The intergranular friction coefficient might also differ. Thus the internal friction coefficient of dry grains, µ * 0 , appears to be larger (over 0.4) in the experiments of Ref. [52]. However, the interesting semi-quantitative agreement between our results and those published data is to be noted. Unlike Pierrat et al., Richefeu et al., in the experimental part of their 2006 paper [19] explicitly assume a Mohr-Coulomb form for the yield criterion. They obtain, with P * values of a few units, a macroscopic cohesion c agreeing with a theoretical formula which coincides with our estimate c = zΦµ ∞ 0 F 0 πa 2 up to a factor of about 1.5. Given the uncertainties on the measurements of c, µ * and coordination number z in that study, our results are quite compatible with their experimental (and numerical) data as well.

In the following sections, for a better assessment of the rheophysical effect of attractive capillary forces, microscopic and microstructural aspects of force networks are investigated in greater detail.

VI. FORCE DISTRIBUTION

The distribution of intergranular force values in a granular material in equilibrium [START_REF] Agnolin | [END_REF][53][54][55][56], or in inertial flow [1,5] has received a lot of attention in the recent literature. While the probability distribution function of force values in cohesionless systems tends to decrease exponentially, on a scale given by the average F N , in cohesive granular assemblies, characterized by the contact tensile strength F 0 , the equilibrium force distribution evolves, as P * decreases to low values, towards a roughly symmetric distribution about zero, with values of both signs of order F 0 F N [13,15]. As compared to the two-dimensional results of Refs. [13,15], the present 3D numerical study of wet spherical grain assemblies does not investigate very small P * states, but involves longerranged distant interactions. The positive force wing of the p.d.f. of normal forces near the quasistatic limit is shown in Fig. 20 the cohesionless distribution shape, and the transition to a cohesion-dominated force network with values of order F 0 , ratio F N /F 0 being approximately proportional to P * as discussed in Sec. V B. At low reduced pressure, as for P * = 0.436, it is more appropriate to normalize the distribution by F 0 , as in Fig. 21. This plot shows the influence of inertia parameter I, which is, for large positive values, qualitatively similar to the one observed with dry grains: the distribution widens, large forces being associated with collisions between grains or groups of grains. Another effect of increasing the inertial number is, as expected from the results of Fig. 12, a depletion of the population of contacts, compensated by a greater number of distant grains joined by a meniscus. To understand better the distribution shape for negative values, Fig. 22 distinguishes the distributions of contact and distant (attractive) forces. Contact force distributions exhibit a maximum in zero, with negative values becoming more frequent as I increases. The larger value of the pdf near -F 0 signals then the opening of more contacts. The distant interactions are responsible for the non-monotonic part of the pdf. On the one hand, the sharp maximum near -F 0 signals a large population of grain pairs at close distance, in agreement with the fast increase of z(h) at small h visible in Fig. 14. On the other hand, the increase near the minimum attractive force at rupture distance D 0 merely reflects the slow variation of function F Cap (h) (Fig. 2).

The "effective pressure" concept relies on the assump- tion that the effect of attractive capillary forces are similar to that of a larger applied isotropic stress. One way to test such an idea at the microscopic scale is to compare the distributions of normal elastic forces: if normalized by the average elastic force, related to the effective pressure, those should be independent on P * and similar to the force distribution in a cohesionless system. Fig. 23 compares the distributions of elastic normal forces, normalized by their mean value, for small I and different values of P * . Those distributions are roughly similar, but show, as expected, notable discrepancies for values of order F 0 . The decay for large values is faster in cohesive systems, reflecting a difference in force networks.

VII. AGGLOMERATION

The aggregation of cohesive grains is observed and reported in many numerical and experimental studies, and is exploited in industrial processes [16,57]. It was directly observed in flow of cohesive granular assemblies, both in numerical model materials [16], and in experiments with wet powders [58]. A numerical study of steady state chute flow [59] reports an increase of the number of long-lasting contacts in the presence of cohesive forces. Weber et al. in [60], carried out a detailed study of the effect of capillary forces on agglomerate duration and size. The agglomeration phenomena in steady shear flow is studied here, first, by measuring contact ages and meniscus ages, depending on state parameters. Then, the age-dependent size of clusters is measured, depending on P * and I. These clustering properties are related to the material rheology.

A. Age of contacts and of distant interactions

The distribution of the age τ c of contacts for I = 10 -1 and different values of P * is shown in Fig. 24. P (τ c γ) is the probability distribution of contact ages τ c , expressed as a strain τ c γ. The decrease of P (τ c γ) is slower for smaller P * , showing that for the stronger cohesive forces the contacts survive over larger strain intervals [59,60]. For large enough strains, τ c γ > 0.5, these probability density functions decay with an exponential form, P (τ c γ) ∝ e -τ c /τ0 . Values of decay times τ 0 , given in Tab. VI, increase as we decrease P * . Average contact ages, τ c avg , also provided in the table, show the same behavior (τ c avg is smaller than τ c 0 because the distribution is not exponential for short times -see inset on the figure). Fig. 25 shows the evolution of the pdf with I for two different values of P * , revealing, as expected, that contact ages (in units of 1/ γ) decrease in faster flows. For I ≤ 10 -2 , curves appear to coincide, showing nearly quasistatic behavior. The probability distribution function of the age of interactions P (τ i γ) (i.e., the age of liquid bridges) is also shown in Fig. 26 for different values of P * . Liquid bridges survive for quite large strain intervals, reaching several units with a probability of order 0.1, which increase as P * decreases. Initially, most liquid bridges survive at least for strains of order 0.1. Beyond unit strain curves might be fitted by an exponential function too, defining a decay time τ i 0 . P * -dependent values of τ i 0 and of the average meniscus age τ i avg are listed in Tab. VI. Remarkably, the curves do not present any notable difference for different values of I: the pairs may lose their contacts in faster flows, but they are still bonded with liquid bridges. The age of contacts and of distant interactions thus reveal the formation of aggre- gates in the presence of capillary forces. These clusters are transported by the flow for some distance before they are broken or restructured. They may survive for strain intervals of a few units. , and for all interactions, τ i 0 , obtained by an exponential fit to the data of Fig. 24 and of Fig. 26); average contact age τ c avg and interaction age τ i avg , for different values of P * and I = 0.1. All four times are normalized by shearing time 1/ γ.

B. Clusters

Clusters are defined as sets of grains connected by liquid bonds for a minimum time, τ cl , and the clustering tendency might be appreciated on recording the τ cldependent mass-averaged cluster size:

S cl m = i S 2 i i S i , (19)

VIII. FABRIC ANISOTROPY

The capacity of granular assemblies to form anisotropic force networks is the only origin of shear strength with frictionless grains [5,56,61], and is known to play a central role in the shear strength of frictional grains as well. To understand how contact and distant capillary forces contribute to the shear stress, it is instructive to study the distribution of contact orientations (normal vectors n on the unit sphere Σ), E(n), which, to lowest order, is characterized by the fabric tensor:

F αβ = n α n β = Σ E(n)n α n β d 2 n (20)
The connection between fabric and normal force contri-bution to stresses, or σ N , is quite direct, as one has:

σ N αβ = zΦ 3πa 2 Σ E(n) F N n n α n β d 2 n, (21)
an integral over the unit sphere in which F N n denotes the average normal force carried by the pairs with orientation n. As reported in Sec. V C, the contribution of the normal forces, σ N , amounts to more than 80% of the shear stress. The contribution of fabric parameters F 12 to shear stress σ 12 might be visualized in Fig. 28. On average, if pairs are preferentially oriented with the normal vector within a compression quadrant in the shear flow, then F 12 < 0 will tend to increase the absolute value of σ 12 if forces are positive, and to decrease it they are negative. On the other hand, negative forces will increase the absolute value of σ 12 if preferentially oriented in the extension quadrants in the shear flow. The evolutions of fabric parameters F c 12 and F d 12 , pertaining respectively to contact and distant normal forces, versus I, for different P * values, are displayed in Fig. 29.

F c

12 is negative, signaling the contribution of normal contact forces to shear strength (as F N,c is the dominant contribution to F N , related to P by ( 12)). The largest value, and the largest variation of F c 12 with I, is obtained in the dry system (P * = ∞). The decrease of this anisotropy parameter for smaller P * may be understood in reference to the clustering phenomena and to the larger duration of contacts evidenced in Sec. VII. Longer-lived contacts rotate in the shear flow, and are less favorably oriented in the compression quadrant. Shear flow carries agglomerates for some distance before they break and thus their random tumbling motion increases the isotropy of the contact orientations. In faster flows ( 15. The different rule of meniscus formation (at contact) and breakage (at distance D 0 ) explains, in part, this large fabric anisotropy of attractive forces: approaching particles are not attracted to each other, whence a small number of distant interacting pairs in the compressive quadrant, with a negative contribution to F 12 ; as particles get separated, receding pairs are still attracted to each other, whence a positive contribution to F 12 from the extension quadrant. In the model without meniscus hysteresis, assuming capillary attraction to appear as soon as grains approach within distance D 0 , F d 12 strongly decreases, from 0.14 to about 0.07 at P * = 0.436 and small I.

IX. SUMMARY AND DISCUSSION

The rheological properties of unsaturated granular materials, in which a small amount of wetting liquid, forming liquid bridges and transmitting attractive capillary forces between particles, generalize, in many respects, previous observations on cohesive granular materials, with macroscopic properties exhibiting similar dependences on I and P * . Thus, compared to dry materials, the apparent internal friction coefficient µ * = σ 12 /σ 22 is enhanced (from 0.33 to more than 1 in the explored range P * ≥ 0.1); looser structures are stabilized, even in the quasistatic limit (Φ 0.52, for P * = 0.436 is below all packing densities with cohesionless grains), even though contact coordination numbers, due to the absence of rattlers, may be larger. Our results describe those effects in quantitative form, in the range P * ≥ 0.1 and 10 -4 ≤ I ≤ 0.56, and specify the dependence on various features of the model. We only predict quite a small dependence of the rheology on saturation within the pendular range (up to 5-10%), in agreement with experimental observations [19,52].

More accurate models of the capillary force dependence on intergranular distance, or of the distribution of liquid between menisci with varying volumes, would hardly change the results. Interestingly, though, some variants of the model, although arguably not realistic, have notably different rheological properties. Thus, reducing the meniscus volume to very small values would have quite a notable effect on internal friction and density -but, in practice, menisci are unlikely to form with such small liquid contents. Assuming menisci form as soon as grains approach to their maximum extension distance (range of capillary force) would also strongly affect macroscopic properties.

Shear localization systematically affects shear flows at low P * , and we could not measure the constitutive behavior at P * = 0.1 except for some intermediate I values of order 0.01. Localized states are characterized by velocity profiles with gradients concentrated within narrow bands, where the solid fraction is well below its bulk value. The band thickness lies in the range of 5 to 10 grain diameters a at small I, but might be as small as about 1.5a in faster flows (I ≥ 0.1).

We also record normal stress differences, which are larger than for dry grains, and tend to grow with decreasing P * . The second normal stress difference, in particular, reaches 20% of the imposed normal stress σ 22 in the quasistatic limit for small P * .

The effective pressure approach to the yield criterion of wet grains ignores such sophistications, as well as density or microstructural changes due to capillary forces. It assumes critical states to be in correspondence for different values of P * , as though the introduction of capillary forces, pushing grains against their neighbors, were equivalent to the application of a larger confining pressure. Such a crude approach is in fact surprising successful as a rough approximation, to predict the increase of µ * for decreasing, but not too small, P * , say P * ≥ 1 (below P * = 1, the increase of µ * is overestimated). The effective pressure might be evaluated upon adding, to the applied pressure, the capillary contribution to the average normal stress. This contribution might itself be estimated from density and coordination numbers, which leads to a Mohr-Coulomb form (1) for the variation of shear strength µ * 0 σ 22 with normal stress σ 22 (if Φ and z do not vary too much). Such a form of the critical state plasticity criterion proves however inadequate to describe the whole range of reduced pressures P * : data are incompatible with a P * -independent macroscopic cohesion c. At low P * , the important differences in micromorphology and force networks are such that one cannot simply relate the properties of wet granular assemblies to those of dry ones. The flow configurations become looser, yet better coordinated, and force distributions are dominated by the characteristic scale of adhesion forces.

Meanwhile, the fabric anisotropy of the contact network become smaller, so that the stress anisotropy has to stem, to a larger extent, from anisotropically distributed force values.

Another remarkable feature of the measured rheology, compared to dry granular materials, is the slow variation with I of both macroscopic rheological parameters such as µ * , and microscopic data such as coordination numbers. Although we could not observe a direct correlation, a slower increase of function µ * (I) could signal a shear banding tendency.

Many of those rheophysical features are explained by, or, at least, related to, the strong clustering tendency emerging as attractive forces gradually become the dominant ones, upon decreasing P * . As contacts are stabilized by attractive forces, they do not so easily open as the network is being sheared. When they do (which happens preferentially in the extension direction within the average shear flow, whence a fabric anisotropy of distant interactions contributing to shear strength), the network of grains bonded by liquid bridges might still be connected, forming enduring connected clusters. The survival of such clusters over quite notable strain intervals (reaching several unities with sizable probability) should limit the dilating tendency of faster flows. It also maintains a network in which the capillary forces act in closer similarity to an effective pressure. In broad, qualitative terms, the capillary forces, which have a finite range, maintain some kind of effective compression, thereby limiting the disruptive effects of collisions on the contact network that are experienced by dry granular assemblies.

Admittedly, this is still a descriptive rheophysical scenario. More quantitative studies should be carried out to better characterize the deformation mechanisms of the grain clusters. The shear banding phenomenon certainly deserves detailed investigations, in which sample size and shape effects should be systematically assessed, and partly localized velocity and density fields analyzed and related to a stability analysis.

Our results in the quasistatic limit are in agreement with the limited available experimental results, as regards the enhancement of shear resistance brought about by capillary forces [19,52]. More laboratory data should however be used, with, if possible information both on rheology and on micromorphology and liquid distribution, for better confrontations of numerical simulations and experimental data. In particular, experiments would be especially valuable to guide the design of complete numerical models, capable of dealing with saturations exceeding the limited pendular range, and of describing the liquid motion. Such numerical models could be applied to the mixing process of the grains with the liquid as well as to the rheology of the mixture. The Lattice Boltzmann method for a diphasic interstitial fluid medium, coupled to a DEM description of grain motion, is a promising perspective [62,63].

FIG. 1 :

 1 FIG.1:(Color online) A meniscus between two spherical grains of diameter a = 2R, with distance h between solid surfaces, filling angle ϕ, contact angle θ.

3 FIG. 2 :

 32 FIG. 2: (Color online) Force law F Cap (h), for two different meniscus volumes, according to the Maugis model and to the Soulié formula.

FIG. 3 :

 3 FIG. 3: (Color online) Solid fraction Φ versus shear strain γ. Time series is obtained with P * = 1, I = 0.1 and N = 4000 when the rupture distance is D0 = 0.1.

FIG. 4 :

 4 FIG. 4: (Color online) Shear stress |σ 12 | (lower curve, red, left axis) and normal stress σ 22 (upper curve, blue, right axis) versus shear strain γ. Note the different scales on left and right axes. Time series obtained with P * = 1, I = 0.1 and N = 4000 when the rupture distance is D0 = 0.1.

FIG. 5 :

 5 FIG. 5: (Color online) Velocity profile for P * = 0.1, I = 0.178 (a) and for P * = 0.436 and I = 10 -3 (b) at different shear strain values.

  value ∆ = 1 in the case of a perfect localization within a plane,

FIG. 6 :

 6 FIG. 6: (Color online) Velocity profiles (lower axis, (red) dotted curve) and local densities of grains (upper axis, (blue) continuous curve) for two configurations with shear strain γ = 20 (a) and γ = 352 (b) when P * = 0.1 and I = 0.178. Average solid fraction, Φ s , is shown as vertical (blue) dashed line, with value of 0.44 in (a) and 0.49 in (b).

Fig. 7 FIG. 7 :

 77 Fig.7is the plot of ∆ as a function of strain γ corresponding to the same simulation as in Fig.5a. It initially shows small fluctuations near zero, and suddenly increases near γ = 250 when the velocity gradient localizes in a shear band.

FIG. 8 :

 8 FIG. 8: (Color online) Macroscopic friction coefficient µ * (a) and solid fraction Φ (b) versus inertial number I for different values of reduced pressure P * (D0 = 0.1a = V 1/3 ).

2 IFIG. 9 :

 29 FIG. 9: (Color online) First (a) and second (b) normal stress differences as functions of I for different P * values. The same symbols and color codes apply to both figures. Both N1 and N2 decrease as P * increases. N1, for large P * , changes sign as a function of I (N1 = 0 is visualized by the thin horizontal line in graph (a)).

6 FIG. 10 :

 610 FIG. 10: (Color online) Macroscopic friction coefficient µ * (a) and solid fraction Φ (b) versus inertial number I for different values of P * and meniscus volume V (with D0 = V 1/3 ).

FIG. 11 :

 11 FIG. 11: (Color online) Macroscopic friction coefficient µ * (a) and solid fraction Φ (b) versus inertial number I for P * = 0.436 and D0 = 0.1, with and without capillary hysteresis.

FIG. 12 :

 12 FIG. 12: (Color online) Coordination numbers: (a) of contacts, zc; (b) of distant interactions, z d .

FIG. 13 :

 13 FIG. 13: (Color online) Total coordination number z = zc + z d .

FIG. 14 :

 14 FIG. 14: (Color online) Coordination number of neighbor grains versus interparticle distance h, for different P * values and for I = 10 -3 . Inset shows detail at small h.

1 FIG. 15 :

 115 FIG. 15: (Color online) Contribution of distant interactions to shear stress σ 12 (a) and to normal stress σ 22 (b).

1 FIG. 16 :

 116 FIG. 16: (Color online) Contribution of capillary forces to stress σ 22 .

1 FIG. 17 :

 117 FIG. 17: (Color online) Contributions of tangential (a) and capillary (b) interactions to total shear stress σ 12 .

FIG. 18 :

 18 FIG.18:(Color online) σ 12 versus σ 22 , in quasistatic flow, in units of F0/a 2 . Square dots: numerical results (error bars are smaller); red crosses: predictions of (17), with exact coefficient β; blue circles: same with estimated β. Insert: detail of numerical data for small P * , including additional point at P * = 0.1.

FIG. 20 :

 20 FIG. 20: (Color online) Distributions of normal forces normalized by average normal force F N for I = 10 -3 and different values of P * .

3 FIG. 21 :FIG. 22 :

 32122 FIG. 21: (Color online) Distributions of normal forces for P * = 0.436 and different values of I, normalized with F0.

FIG. 23 :

 23 FIG. 23: (Color online) Distributions of normal elastic forces at small I for different P * values.

2 FIG. 24 : 3 FIG. 25 :

 224325 FIG. 24: (Color online) Distribution of the age of contacts for different values of P * and I = 0.1. Inset shows the same graph in a shorter range of τ c γ.

5 FIG. 26 :

 526 FIG. 26: (Color online) Distribution of the age of menisci for different values of P * (same for all values of I). Inset shows the same graph in a shorter range of τ i γ.

  in which the summations run over all clusters i containing S i grains. The values of S cl m for the different values of P * and I are shown in Fig. 27, as functions of τ cl γ. For small values of τ cl γ almost all particles are gathered in a

FIG. 28 :FIG. 27 :

 2827 FIG. 28: (Color online) Sketch of approaching (ϕ near 3π/4) and receding pairs (ϕ near π/4) in macroscopic shear flow.

FIG. 29 :

 29 FIG. 29: (Color online) Fabric parameter F12 for contacting pairs (a) and distant interactions (b) versus I for different P * .

TABLE II :

 II Parameters of the fit of function µ

	P *	µ * 0	α	c
	0.436	0.867 ± 0.003	0.70 ± 0.05	0.30 ± 0.01
	1	0.607 ± 0.003	0.76 ± 0.05	0.37 ± 0.02
	2	0.473 ± 0.007	0.72 ± 0.06	0.42 ± 0.02
	5	0.387 ± 0.006	0.70 ± 0.05	0.46 ± 0.02
	10	0.366 ± 0.004	0.74 ± 0.04	0.48 ± 0.02
	∞	0.332 ± 0.004	0.71 ± 0.03	0.50 ± 0.01
	P *	Φ0	ν	e
	0.4360 0.5243 ± 2.10 -4 1.73 ± 0.05 0.497 ± 0.017
	1	0.5559 ± 10 -4 1.34 ± 0.012 0.512 ± 0.005
	2	0.5726 ± 10 -4	1.21 ± 0.01 0.547 ± 0.003
	5	0.5851 ± 10 -4	1.12 ± 0.01 0.580 ± 0.003
	10	0.5900 ± 10 -4	1.09 ± 0.01 0.594 ± 0.004
	∞	0.5970 ± 10 -4 0.96 ± 0.015 0.562 ± 0.008

1, measured in reasonably homogeneous flows ∆ ≤ 0.1), indicate µ * 1.6 for I = 10 -2 . Meanwhile, the material becomes looser, with Φ reaching values that cannot be observed without cohesion in quasistatic conditions. * (I) by Eq. 10, for different values of P * .

TABLE III :

 III Parameters of the fit of function Φ(I) by Eq. 10, for different values of P

* .

TABLE IV :

 IV ). Effect of meniscus volume or saturation level on different parameters for I = 10 -2 and P

	V /a 3	Sw	z	Φ	µ *
	10 -2	7.137 × 10 -2	6.863	0.520	1.071
	5 × 10 -3	3.418 × 10 -2	6.556	0.522	1.003
	10 -3	6.305 × 10 -3	5.970	0.524	0.875
	2 × 10 -4	1.075 × 10 -3	5.534	0.525	0.787
	10 -6	5.539 × 10 -6	4.836	0.530	0.661

* = 0.436.

  for larger I), while contacts tend to open in cohesionless materials, enhancing fabric anisotropy, cohesive contacts can resist flow agitation and inertial effects better, whence a smaller I influence; if they open, they transform into attractive distant interactions, and the anisotropy of distant interactions also decreases. Those distant capillary forces are characterized by a comparatively large anisotropy, about three times as large as |F c 12 |. As F d 12 is positive, those distant attractive forces contribute to increase the internal friction coefficient. Unlike |F c 12 |, F d 12 increases for smaller P * values, which corresponds to the growing contribution of distant interactions to shear strength shown in Fig.


B Chareyre

Saeed Khamseh

Jean-Noël Roux

François Chevoir

About "Flow of wet granular materials: a numerical study". Response to a Comment by

Keywords: numbers: 83.80.Fg, 47.57.Qk Keywords: granular flow, capillary cohesion, Mohr-Coulomb criterion

, is bound to fail in strongly cohesion-dominated material states.

In our 2015 Physical Review E paper [START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF] we investigate, by discrete numerical simulations, the constitutive laws of wet spherical grains in steady shear flow, in relation to micromorphology and micromechanics. Assemblies of beads of diameter a, made of a material with mass density ρ, are subjected to a slow shear flow, with controlled shear rate γ = ∂v 1 ∂x 2 and normal stress σ 22 .

(index 1 labelling the flow direction and index 2 the velocity gradient direction). We characterize the material behavior through the apparent friction coefficient, µ * , conventionally defined, as in cohesionless materials, as the ratio of the shear stress to the normal stress: µ * = σ 12 /σ 22 . One major control parameter is the reduced pressure P * comparing the controlled normal stress σ 22 to the capillary attractive force F 0 between contacting grains: P * = a 2 σ 22 /F 0 . Together with inertial number I = γa ρ/σ 22 , reduced pressure P * determines the internal state of the material in steady uniform shear flow. The issue discussed in the Comment is the definition of effective stresses σ eff , such that, in the quasistatic limit of I → 0, the flow condition (yield criterion) takes the form σ eff 12 = µ * ∞ 0 σ eff 22 , with the static internal friction coefficient µ * ∞ 0 of the dry material (corresponding to the limits of P * → ∞ and I → 0).

In our paper, we argue that the effect of capillary forces, contributing σ cap to the total stress, could be likened to a shift in stresses transmitted through the contact network, and we adopt these shifted stress values as effective stresses. We choose nevertheless to ignore the shift in shear stress, and indeed our simplifying approach amounts to writing σ eff 12 = σ 12 and σ eff 22 = σ 22 -σ cap 22 . This omission is criticized by Chareyre in his Comment, * Electronic address: jean-noel.roux@ifsttar.fr as entailing poorer predictions of shear strength. He points out that on accounting for all terms of the capillary stresses, the quasistatic criterion written as in (his) Eq. 1 very closely fits all simulation data points (Fig. 1) pertaining to homogeneous quasistatic flows. As, on ignoring σ cap 12 , we obtain a relative error growing to 25% at small P * , he therefore concludes that we are unduly pessimistic in our assessment of the performance of shear strength predictions based on the effective stress approach. We find this conclusion well founded for the main data series of Ref. [START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF], we agree that the success of the assumption that the contact stresses (defined as σ cont = σ-σ cap ) are effective stresses is remarkable, and that this observation nicely agrees with other results [START_REF] Scholtès | Micromechanics of granular materials with capillary effects[END_REF]. Thus, we regard the Comment as very relevant and interesting.

The perspectives of our work [START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF] (see also [START_REF] Than | Basic Mechanical Properties of Wet Granular Materials: A DEM Study[END_REF]), in which many other aspects of wet granular material rheology are addressed, were actually slightly different. Our main motivation on invoking effective stress ideas was to provide a physical explanation for the strong increase of friction coefficient µ * 0 , the quasistatic limit of µ * , as P * decreases through the investigated range. In this respect, the simplified approach ignoring σ cap 12 proves quite successful, as it provides a semi-quantitative prediction of the shear strength with simple means. As an additional advantage, it justifies the existence of a cohesion in the sense of Mohr-Coulomb, with a value that relates in satisfactory approximation to F 0 , solid fraction Φ and coordination number z as

c = zΦµ * ∞ 0 F 0 πa 2 . (1)
This lends itself to interesting comparisons with the literature [START_REF] Pierrat | Effect of moisture on the yield locus of granular materials: theory of shift[END_REF][START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF] and correctly assesses cohesive contributions to strength for P * ≥ 1. Using the full contribution of capillary forces to the stresses is more accurate, but σ cap 12 is not as easily related to the simplest state variables Φ and z. In recent simulations with a slightly different sys-tem [6], we also observed σ cap 12 to vary significantly with P * .

A more complete exploration of the range of validity of postulate σ eff = σ cont is still needed, as its physical origin remains mysterious. The basic objections put forth in [START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF], about the microstructure differences induced by cohesion, should bear some relevance. Cohesive systems might form stable structures with very low densities [START_REF] Than | Basic Mechanical Properties of Wet Granular Materials: A DEM Study[END_REF] (e.g., Φ = 0.25) inaccessible to cohesionless ones, in which case the comparison with dry materials subject to the effective stress makes no sense. Such configurations are however associated with low P * values for which homogeneous shear flow were not observed. Similar numerical simulations of wet bead assemblies in shear flow were recently carried out in our group [6], with a different value of the intergranular friction coefficient, µ = 0.09, motivated by quantitative comparisons with an experimental system. The effective stress postulate (accounting for σ cap 12 ) is observed to predict the shear resistance with an accuracy better than 5%, except for P * < 1 and small liquid contents (V m < 10 -4 a 3 with the notations of [START_REF] Khamseh | Flow of wet granular materials: a numerical study[END_REF]), in which case the error grows to order 10%.