Jérôme Duriez

Richard Wan

Stress in wet granular media with interfaces via homogenization and discrete element approaches

The nature of the stress tensor for an unsaturated pendular-state granular medium is investigated following two micromechanical approaches. Firstly, a stress tensor is analytically derived through stress-homogenization of the medium with internal surfaces being explicitly incorporated in addition to the solid, liquid and gaseous volumes. As such, the derivation identifies a surface stress tensor associated with the liquid-gas interface endowed with distributed surface tension forces. Secondly, numerical simulations of unsaturated conditions are pursued within the Discrete Element Method (DEM) which can only consider resultant point forces, while actual internal forces are indeed distributed in nature, e.g., the liquid pressure which acts over the wetted surfaces. Despite this shortcoming, stress descriptions provided by the above two fundamentally distinct approaches are found to be equivalent for unsaturated media subjected to mechanical and hydraulic loading in the pendular regime.

Moreover, both approaches indicate that the capillary stress, interpreted as the part of the total stress representing the combined effects of the liquid and gas phases and interfaces, is driven by the microstructure and is thus generally non-spherical.

INTRODUCTION

Granular materials form a special class of porous media that are ubiquitously encountered in various types of engineering applications. In the realm of civil engineering, granular soils such as sand and gravel are widely present, whereas powders of varied nature are equally relevant in materials science, as well as in food and chemical industries. In the most general 1 Dept. of Civil Engrg., Univ. of Calgary. E-mail: jerome.duriez@ucalgary.ca 2 Dept. of Civil Engrg., Univ. of Calgary. E-mail: wan@ucalgary.ca 1 case, such granular materials are not dry, but their pore spaces may be partially saturated with a liquid while the remaining space is occupied by a gas. For instance, granular soils are subjected to such unsaturated conditions above the groundwater table in the so-called vadose zone where water is present in the smallest pores due to capillary rise, whereas air occupies the bigger pores. Such unsaturated conditions have a drastic influence on the mechanical properties of wet materials, as well-illustrated by the change in consistency of moist food or wet sand from the dry state. This practical engineering issue has taken on a more theoretical importance recently since a thorough understanding of how the liquid-gas mixture in the pore spaces changes the mechanical properties of granular materials is still lacking. One of the long-standing controversies in geotechnical engineering surrounds the question of whether the effective stress concept [START_REF] Terzaghi | Soil Mechanics in Engineering Practice Third Edition[END_REF]) can be theoretically applied to unsaturated conditions so that the mechanical behavior of an unsaturated soil can be predicted from its properties in dry conditions [START_REF] Nuth | Effective stress concept in unsaturated soils: Clarification and validation of a unified framework[END_REF].

Given the inherent microstructure of granular media, the above question can be readily addressed by adapting a multiscale approach that has clearly been proven successful in dry conditions [START_REF] Bathrust | Observations on stress-force-fabric relationships in idealized granular materials[END_REF]Bagi 1996;[START_REF] Wan | Anisotropy of a tensorial Bishop's coefficient for wetted granular materials[END_REF][START_REF] Duriez | Constitutive modelling of cohesionless soils and interfaces with various internal states: An elasto-plastic approach[END_REF].

Generally speaking, multiscale modeling approaches can be either numerical or analytical in nature. For instance, Chateau andDormieux (1995, 2002), [START_REF] Gray | Analysis of the solid phase stress tensor in multiphase porous media[END_REF], [START_REF] Nikooee | Effective stress in unsaturated soils: A thermodynamic approach based on the interfacial energy and hydromechanical coupling[END_REF], [START_REF] Madeo | A continuum model for deformable, second gradient porous media partially saturated with compressible fluids[END_REF], [START_REF] Wan | Micromechanical analysis of force transport in wet granular soils[END_REF] developed analytically based formulations for the multiscale description of unsaturated porous media, whereas [START_REF] Gili | Microstructural deformation mechanisms of unsaturated granular soils[END_REF], [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF][START_REF] Wan | Anisotropy of a tensorial Bishop's coefficient for wetted granular materials[END_REF][START_REF] Wan | Anisotropy of a tensorial Bishop's coefficient for wetted granular materials[END_REF], to name a few, treated the condition of unsaturation numerically using the Discrete Element Method (DEM). According to the standard DEM calculation paradigm, unsaturated conditions in a wet granular medium are idealized by reducing all internal forces arising from fluid pressures and surface tension phenomena into resultant interaction capillary forces that are conveniently applied to a pair of interacting particles or so-called discrete elements (DE). This subtle difference bears no consequence for what concerns the computations of displacements of rigid solid particles or is that total stresses as computed in DEM are necessarily expressed in terms of these resultant point forces solely [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF], suggesting a possible breakdown of the model since resultant forces and internal stresses are distinct mechanical concepts. By contrast, the distributed nature of actual internal forces, e.g. the liquid pressure that acts along wetted surfaces, can be readily incorporated into analytical derivations of the total stress tensor based on homogenization techniques; see [START_REF] Wan | Micromechanical analysis of force transport in wet granular soils[END_REF], or even earlier, Chateau andDormieux (1995, 2002). The central point is to distinguish any model prediction difference between the two calculation paradigms that both provide the total stress for an assemblage of wet particles. In this connection, Wan et al. (2015) identified stress calculation discrepancies between the DEM and analytical homogenization approaches when the liquid volumes become non-negligible. Such result appears to be logical since liquid pressure imparts mechanical actions on the particles that clearly deviate from the point force nature assumed in DEM calculations for significant liquid contents.

It appears that the brief comparison presented by Wan et al. (2015) rests on an incomplete homogenization formula with regard to the exact inclusion of liquid-gas interfaces through the so-called contractile skin [START_REF] Fredlund | Stress state variables for unsaturated soils[END_REF]) that is a dominant characteristic of wet granular media. For this reason, we aim in this paper to revisit the comparison of stress descriptions for wet media provided by DEM and analytical homogenization. A second objective is to underscore the liquid-gas interface contribution to the total stress of a wet granular soil through a rigorous derivation of an interface term in the form of a surface stress tensor. It turns out that the final total stress expression, herein based on volume averaging, is equivalent to the one presented by [START_REF] Chateau | Homogenization of a non-saturated porous medium: Hill's lemma and applications[END_REF] via the principle of virtual work. Although the two expression are the same, our formulation has been derived independently and an extended presentation is herein proposed. Then, the DEM modeling of wet granular media is presented, before addressing the comparison of the two stress descriptions. Once the stress descriptions have been carefully examined, the final section provides some micromechanical insights on the stress state of granular soils in unsaturated conditions.

It is noted that most DEM models used in unsaturated conditions apply to the pendular regime, e.g. [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF], [START_REF] Wan | Anisotropy of a tensorial Bishop's coefficient for wetted granular materials[END_REF], and this saturation regime will be the scope of most sections, except the opening section which is written for the most general case. Herein, the pendular regime is defined by low liquid content, such that the liquid phase forms distinct menisci or bridges between solid particle pairs.

STRESS DESCRIPTION VIA A HOMOGENIZATION APPROACH

We consider a representative elementary volume (REV) V of unsaturated granular soil as a ternary mixture of solid (s), liquid (l) and gaseous (g) phases; see Fig. dition to the partitioned solid, liquid and gaseous volumes V α (α " s, l, g), we further distinguish within V the liquid-gas interface S lg associated with surface tension forces and specific energies as evidenced by the spherical shape that bubbles or droplets adopt under ideal conditions. Such a consideration dating back to the early works of [START_REF] Morrow | Physics and thermodynamics of capillary action in porous media[END_REF] and [START_REF] Fredlund | Stress state variables for unsaturated soils[END_REF] is also present in formulations presented by Chateau andDormieux (1995, 2002), [START_REF] Gray | Analysis of the solid phase stress tensor in multiphase porous media[END_REF], [START_REF] Nikooee | Effective stress in unsaturated soils: A thermodynamic approach based on the interfacial energy and hydromechanical coupling[END_REF], [START_REF] Madeo | A continuum model for deformable, second gradient porous media partially saturated with compressible fluids[END_REF], among others. Thus, we also endeavour to apply the same treatment to our previous works [START_REF] Wan | Micromechanical analysis of force transport in wet granular soils[END_REF]Wan et al. 2015) by explicitly accounting for such liquid-gas interfaces in the relevant average calculations. As a prelude and to provide the background of subsequent developments at the REV level, we will first propose an extended presentation of the key role of interfaces in the micromechanics of mixtures.

Interface description

Considering a system of two immiscible fluids α and β, the corresponding interface S αβ includes material points belonging to both fluids; see Fig. 2. Approaching the interface from the α-side of S αβ , material points experience mechanical actions exerted by the internal forces within α arising from σ α " u α δ, external tractions from the β-side, as well as surface tension forces proportional to γ αβ . A similar argument can be extended to the case where we approach the interface from the β-side, without changing the final result.

Enforcing equilibrium conditions to the interface S αβ leads to:

ż S αβ pσ α ´σβ q n dS " ż C γ αβ ν dl " ż C γ αβ n ˆ dl (1)
In Eq. ( 1), it is considered that the unit normal vector n points from α to β with the surface contour C " BS αβ , while ν is an inward conormal and dl " || dl|| the arc length, as illustrated in Fig. 2.

ν n C S αβ α β dl FIG. 2. Interface S αβ
Assuming surface tension γ αβ to be homogeneous, the right hand side of Eq. ( 1) is shown to be equal to the surface integral, along S αβ , of γ αβ div p nq n (see, e.g., Wan et al. 2015). A boundary condition valid along S αβ and including capillarity effects then follows:

pσ α ´σβ q n " γ αβ div p nq n (2)
which is a statement of stress balance at the interface. The r.h.s. of Eq. ( 2) would be zero in the absence of capillarity, thus coinciding with the classical stress continuity condition.

However, for the case at hand where α and β are immiscible fluids exerting pressures u α and u β respectively, Eq. ( 2) indeed reveals that, in crossing the interface, the fluid pressure undergoes a jump by an amount that depends on the curvature div p nq and the surface tension γ αβ , corresponding to the Young-Laplace equation, i.e.

u α ´uβ " γ αβ div p nq (3)
Considering the interface as a distinct medium or phase, a stress-like tensor π αβ can be invoked to describe the surface tension force field as internal forces for the interface. Adopting soil mechanics sign convention with compressive stresses being positive, and considering an elementary surface S Ă S αβ with the inward conormal ν (see Fig. 2), π αβ necessarily obeys the following equation to describe the tensile state associated with surface tension:

π αβ ν " ´γαβ ν @ x P BSzBS αβ (4)
so that an adequate expression for π αβ takes the form:

π αβ " γ αβ p n b n ´δq (5)
Eq. ( 5) basically expresses π αβ in terms of the projection tensor onto S αβ , i.e. pδ n b nq. Thermodynamic justifications for such an expression have been presented by e.g., [START_REF] Chateau | Homogenization of a non-saturated porous medium: Hill's lemma and applications[END_REF], [START_REF] Gray | Analysis of the solid phase stress tensor in multiphase porous media[END_REF], demonstrating that π αβ is workconjuguate to an interface kinematics field.

Stress homogenization

Let us now turn back to the REV level where the unsaturated mixture includes the solid, liquid and gaseous phases V α (α " s, l, g) as well as the liquid-gas interface phase S lg as a new addition. The macroscopic stress tensor for the REV, Σ, is obtained by an averaging of 'microscopic' stresses existing in all of the above-mentioned four phases:

Σ " 1 V ˜ÿ α"s,l,g ż Vα σ α dV `żS lg π lg dS ¸(6)
The stress tensors in the fluid phases are readily expressed in terms of their respective uniform pressure, i.e. σ α " u α δ @ x P V α (α " l, g). As for the solid phase, we may classically transform the volume integral into surface integral from tensor calculus, since for any continuum in equilibrium under no body force:

ż V σ dV " ż S pσ nq b x dS (7)
Here, Eq. ( 7) is applied successively to all particles p that constitute the solid phase of granular media. The particles surfaces S p include:

• the contact lines Γ p where the three phases intersect and liquid-gas surface tension imparts on solid particles,

• the wetted and non-wetted surfaces (resp. S p,sl and S p,sg ) over which fluid pressures (resp. u l and u g ) act, and

• finally, contact surfaces between solid particles giving rise to contact forces. Considering rigid solid particles, such contact surfaces reduce to contact points x c .

Then, after some algebraic manipulations, we finally obtain:

Σ´u g δ " 1 V « ÿ c f c b l ´s ˆżS sl n b x dS `Vl δ ˙´ż Γ γ lg b x dl ´żS lg γ lg pδ ´ n b nq dS ff (8)
The first r.h.s. term of Eq. ( 8) considers all contacts c between two solid particles 1-2, with f c the contact force exerted by 1 on 2, and l the so-called branch vector linking the centre of 1 to that of 2. This indeed corresponds to the classical Love-Weber formula that expresses Σ in dry conditions. For easier reference, this stress contribution is denoted as the contact stress tensor σ cont throughout the manuscript. The contact stress tensor is specific to the solid phase and appears as a natural choice for an effective stress that governs the behavior of the granular skeleton in unsaturated conditions, as considered by, e.g., [START_REF] Lu | Suction stress characteristic curve for unsaturated soil[END_REF]. As a matter of fact, it has been shown to unify the failure description in dry and unsaturated conditions from numerical data [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF][START_REF] Wan | Micromechanical analysis of force transport in wet granular soils[END_REF]Wan et al. 2015) or experimental ones [START_REF] Lu | Suction stress characteristic curve for unsaturated soil[END_REF]. However, a comprehensive description of the stress-strain behavior of the unsaturated REV would require accounting for the coupling between the various phases in addition to the only knowledge of σ cont [START_REF] Madeo | A continuum model for deformable, second gradient porous media partially saturated with compressible fluids[END_REF]Wan et al. 2015).

The second r.h.s. term of Eq. ( 8) describes mainly the action of the suction, or capillary pressure, s " pu g ´ul q ą 0 along the wetted surfaces S sl " Ť p S p,sl . Interestingly, this term is not spherical in the general case as it depends on the microstructure of the fluid phase distribution. Thus, it cannot be normally associated with an averaged fluid pressure that would act equally in any direction, as assumed in the classical Bishop's expression of an effective stress [START_REF] Bishop | The principle of effective stress[END_REF][START_REF] Bishop | Some aspects of effective stress in saturated and partly saturated soils[END_REF]. This point will be developed in more details in the next sections.

As for the third r.h.s. term, Γ is the set of contact lines where the three phases intersect:

Γ " Ť p Γ p , and γ lg dl " γ lg ν dl is an infinitesimal surface tension force as experienced by the solid particles along Γ p . This third tensorial term depends again on the fluid phase distribution, and also on the wettability of the soil particles expressed by the contact angle.

The fourth and last r.h.s. term represents the surface tension contribution to the total stresses from the liquid-gas interfaces as a membrane-like stress Π lg {V :

1 V Π lg " 1 V ż S lg π lg dS (9)
We denote as capillary stresses σ cap [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF]Wan et al. 2015) the part of the total stresses that is due to the fluid mixture contribution, i.e.

σ cap " Σ ´ug δ ´σcont " ´1 V « s ˆVl δ `żS sl n b x dS ˙`γ lg ˜żS lg pδ ´ n b nq dS `żΓ ν b x dl ¸ff (10)
Such capillary stress terminology that we use in this paper corresponds to the suction stress defined by [START_REF] Lu | Suction stress characteristic curve for unsaturated soil[END_REF] and [START_REF] Nikooee | Effective stress in unsaturated soils: A thermodynamic approach based on the interfacial energy and hydromechanical coupling[END_REF]. As it will be emphasized in the next sections, the capillary stresses not only directly depend on both the suction s and the surface tension γ lg , but also on the microstructure of the fluid distribution (V l , S sl , S lg , Γ).

Note finally that Eq. ( 8) is consistent with Terzaghi's equation when conditions of full saturation are considered. In this case S lg " Γ " H and S sl " S s (there are no contact surfaces since we consider rigid particles) and we have furthermore:

ż Ss n i x j dS " ż Vs Bx j Bx i dV " ż Vs δ ij dV " V s δ ij (11)
such that σ cap " ´s δ and Σ " σ cont `ul δ is finally obtained, as in Terzaghi's equation.

Stresses in an idealized granular material

We specialize, from now on, the derivations to idealized assemblies of spherical particles p which can be of different radii R p so that x " R p n @ x P S p , and the capillary stress in Eq.

(10) takes on the following form:

σ cap " ´1 V « s ˜Vl δ `ÿ p R p ż S p,sl n b n dS ¸`γ lg ˜żS lg pδ ´ n b nq dS `ÿ p R p ż Γp ν b n dl ¸ff (12)
We next focus on the isotropic part of Eq. ( 8) to identify the mean stress, i.e.

p ´ug " p cont `pcap

p cont " 1 3 trpσ cont q p cap " 1 3 trpσ cap q " ´s ˜n S r `1 3 V ÿ p R p S p,sl ¸´γ lg 3 V ˜2 S lg `ÿ p R p Γ p sin θ ¸(13)
where n is the porosity: n " pV l `Vg q{V , S r the saturation ratio: S r " V l {pV l `Vg q and θ is the contact angle describing the wettability of the solid particles with the liquid (see next Fig. 3); so that we find p cont from Eq. ( 13) as:

p cont " p ´ug `s ˜n S r `1 V ÿ p R p S p,sl ¸`γ lg 1 3 V ˜2 S lg `ÿ p R p Γ p sin θ ¸(14)
Keeping in mind our previous comments on the effective nature of σ cont , it is still useful to draw some parallels between the relation ( 14) providing the mean contact stress in an idealized spherical-particle assembly in unsaturated conditions, and Bishop's equation (1959,1963) classically used for unsaturated soils in the form of p 1 " p´u g `s χ , with p " 1{3 trpΣq, p 1 the mean effective stress, and the effective stress parameter χ P r0; 1s. It becomes evident that the scope of Eq. ( 14) goes beyond Bishop's equation by providing a micromechanical interpretation of χ in terms of degree of saturation and fluid distribution details, and a second term proportional to γ lg that accounts for interfacial tension associated with the contractile skin. The absence of these details in the initial Bishop's equation explains the difficulties encountered in validating the former from experiments (see, e.g., the discussion by [START_REF] Nuth | Effective stress concept in unsaturated soils: Clarification and validation of a unified framework[END_REF]). Yet, another shortcoming of Bishop's expression is the spherical nature it confers to the capillary stress pΣ ´ug δ ´σ1 q, whereas our expression, Eq. ( 10) or ( 12), is more general in nature and admits a deviatoric component. The last part of the paper investigates such aspect, based on a microstructural description of unsaturated soil which allows the explicit calculation of capillary stresses using Eq. ( 12).

DEM DESCRIPTION OF UNSATURATED SOILS

DEM model formulation

For the purpose of this paper, we carried out a DEM implementation of unsaturated condition that extends the version of [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF] to incorporate non-zero contact angles. In the numerical model, an assemblage of spherical DE as an idealized granular system under unsaturated conditions in the pendular regime is simulated considering two types of particle interactions.

The first type of particle interaction is elementary as the contact force between pairwise contacting particles f c is derived from relative displacements between them according to a linear elastic-plastic contact law. As such, there are only three basic parameters involved: Y and P that govern the normal and tangential contact stiffnesses, and ϕ as the local friction angle that restricts the tangential contact force through Coulombic friction law. More details have been presented by Wan et al. (2015).

To accommodate for unsaturated conditions in the pendular regime as a dispersion of liquid bridges and associated capillary forces in between the discrete particles, a second type of particle interaction is necessary. In particular, for a given suction value specified to the particle assemblage, potential distinct liquid bridges, as shown in Fig. 3, are determined for pairs of both contacting and distant particles as follows. First, a comprehensive data

R 2 ≥ R 1 R 1 θ θ δ 2 δ 1 z ρ(z) d FIG.
3. Liquid bridge (half-)geometry set of constant-curvature axisymmetric bridge profiles ρpzq is generated in accordance with Laplace-Young equation. The associated numerical procedure is based on a Taylor's expansion of ρpzq from assumed various boundary conditions and satisfying the Laplace-Young equation through relations between ρ and its derivatives ρ 1 , ρ 2 [START_REF] Lian | A theoretical study of the liquid bridge forces between two rigid spherical bodies[END_REF]. Then, throughout the DEM simulations, potential liquid bridges between pairs of particles are searched interpolating from this data set depending on the contact angle θ, the imposed suction, the interparticle distance d ě 0, as well as the radii R 1 and R 2 ě R 1 . The liquid bridge distribution that results from such a calculation procedure conforms with uniform suction conditions obtained at thermodynamic equilibrium, see [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF] and [START_REF] Wan | Anisotropy of a tensorial Bishop's coefficient for wetted granular materials[END_REF]. In this instance and contrary to [START_REF] Gili | Microstructural deformation mechanisms of unsaturated granular soils[END_REF], no pore flow responsible of liquid transfers is described by the model. Associated capillary forces f cap as a byproduct of Laplace-Young's solution are then readily applied to the two particles, once the liquid bridges are determined:

f cap " π R 1 sin δ 1 ps R 1 sin δ 1 `2 γ lg sinpθ `δ1 qq z " π R 2 sin δ 2 ps R 2 sin δ 2 `2 γ lg sinpθ `δ2 qq z (15)
As previously mentioned, the DEM model is restricted to the pendular regime, and the numerical procedure expounded in the above applies to isolated menisci bridging not more than two particles as depicted in Fig. 3.

All DEM simulations are performed considering a numerical 3-D sample composed of 20,000 particles with a mean diameter D 50 « 0.05 mm (Wan et al. 2015). Under an isotropic pressure of 1 kPa, the porosity of the sample comes out as n « 0.36. As such, the sample displays a typically dense behavior for the range of confining pressures considered here (tenths of kPa). Note that the behavior of the model is particle size-dependent in unsaturated conditions, in line with experimental evidences.

All the parameters of the model are summarized in Table 1. The surface tension value retained in our study refers to air-water interface at 20 0 C. For the purpose of a broader study, different contact angle values are considered in the range r0 0 ; 60 0 s.

DEM description of the liquid phase

Once all liquid bridges between pairs of particles are determined in the DEM model, the detailed liquid phase distribution at both the pore and macro scales, and hence liquid (SWCC) is readily computed in the case of perfect wetting (θ " 0 0 ) for the sample being subjected to hydraulic loading at a constant isotropic stresses of 10 kPa. Since the liquid bridge distribution depends on the particle distribution, i.e. packing, the proposed SWCC is specific for this mechanical state. As liquid condensates primarily along contacting solid surfaces, a pseudo-primary wetting path is simulated disregarding liquid bridges between distant particles. On the other hand, the simulation of a pseudo-primary drying path considers all possible liquid bridges between both contacting and distant particles, as long as a solution for the Laplace-Young equation can be found. We caution that other existing physical mechanisms affecting the liquid pore transfer, such as the contact angle hysteresis, are neglected. This explains the limited hysteresis obtained in the DEM simulations as shown in Fig. 4. However, crudely speaking, the numerical SWCCs are comparable with the experimental curve characteristic of Ottawa sand.

DEM description of the contractile skin

The liquid bridge calculations outlined in the preceding discussion provide other important information such as the liquid bridge profile ρpzq and the associated interface surface, as well as their statistical distibutions throughout the DEM sample. In particular, and of utmost interest, is the explicit computation of the membrane stress tensor Π lg as given in Eq. ( 9). Due to the oriented nature of the meniscus, which does not conform with a spherical shape, the deviatoric component pΠ zz m ´Πxx m q is non-zero. To gain some physical insights in the nature of this membrane stress tensor, examples of the mean and deviatoric components of Π m , in a non-dimensionalized form, as a function of the dimensionless distance d ˚" d{R 2 are illustrated in Fig. 5.

When transferring the calculations to the sample scale, it is nevertheless observed that the deviatoric components of all local membrane stresses associated with each liquid bridge of the REV approximatively cancel each other for both isotropic and anisotropic microstructures with preferred menisci orientations. Under such an instance, the membrane stress is almost spherical, i.e.

Π lg V " 1 V ż S lg π lg dS « ´2 3V γ lg S lg δ (16) 0 0.05 0.1 0.15 0.2 -0.5 -0.4 -0.3 -0.2 -0.1 0 d* (-) Π m / (γ lg R 2 2
) (-) p: s* = 2 q: s*=2 p: s* = 20 q: s* = 20 FIG. 5. Mean (p) and deviatoric (q) components of the membrane stress for selected menisci (R 2 {R 1 " 2; θ " 20 0 and distinct dimensionless suctions s ˚" s R 2 {γ lg )

The specific interfacial area S lg {V as computed by the model upon primary wetting and drying for an isotropic sample (under a constant isotropic stress of 20 kPa) and different contact angles is given Fig. 6. Independently of the hydraulic path and the saturation, lower

DEM stress description

As signalled in the opening section of this paper, the DEM calculation paradigm only admits resultant point forces as internal forces within the particle assemblage, but a total stress tensor can still be readily calculated from the interaction forces applying the celebrated Love-Weber formula to DE, considering quasi-static conditions and no body forces. Thus,

Σ " 1 V ÿ DE ż V σ dV " 1 V ÿ DE ż S σ n b x dS " ´1 V ÿ DE ÿ f f b x (17)
Since we are considering unsaturated conditions both capillary forces and contact forces must be advocated, so that the total stress tensor in Eq. ( 17) finally gives [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF]):

Σ " 1 V ˜ÿ c f c b l `ÿ m f cap b l ¸" σ cont `σcap (18)
where σ cont is exactly the same contact stress tensor than the one identified in the homogenization approach, see Eq. ( 8). The second term of Eq. ( 18) has the same meaning as the capillary stresses σ cap defined in the homogenization approach, Eq. ( 10), both being equal to Σ ´σcont , if gas pressure is neglected. However, the expressions of these two capillary stresses are quite different. On the one hand, Eq. ( 18) expresses the "DEM capillary stresses" directly from the capillary resultant forces ( f cap being here the capillary force acting on the particle 2 of a liquid-bonded 1-2 pair). On the other hand, Eq. ( 10) obtained through the homogenization approach considers the actual distributed nature of the internal forces, e.g.

the fluid pressures and surface tensions acting over the solid surface.

Obviously, the DEM capillary forces f cap , Eq. ( 15), do correspond to the integral of these fluid pressures and surface tensions over the solid surfaces. However, as far as stress computations are concerned where dyadic products of forces are involved, it is different to consider resultant forces, as in DEM, than the physically distributed tractions, as in the homogenization approach. This is because f b x " p ş S σ n dSq b x ‰ ş S σ n b x dS in the general case. For instance, a uniform fluid pressure does not induce any resultant force on solid particles in saturated conditions, whereas it does engender an internal stress. The DEM capillary stress as calculated from Eq. ( 18) is then zero for such saturated conditions, whereas the homogenized capillary stress according to Eq. ( 10) is not.

From the above discussion, it appears that the calculation approach adopted in DEM may be inadequate in describing stresses of unsaturated soils in the case wetted surfaces are significant. As a matter of fact p ş S t dSqb x " ş S tb x dS is recovered for infinitesimal surfaces S, thus guaranteeing consistent stress descriptions for both DEM and the homogenization approaches at very low liquid content. For higher liquid contents remaining still within the pendular regime, discrepancies have been identified by Wan et al. (2015). Because these previous comparisons considered an homogenization approach that omitted the membrane stress given in Eq. ( 9), the next sections readdress this contentious issue.

DEM VS HOMOGENIZATION STRESS DESCRIPTIONS

Isotropic case

As a first case, we consider a sample with isotropic distributions of contact normals and fluid phases. This results into an isotropic (spherical) capillary stress σ cap " p cap δ. For comparison purposes, Eq. ( 19) recalls the expression of p cap as derived from the homogenization approach, whereas Eq. ( 20) gives p cap according to the DEM stress description:

p cap hom " ´1 3 V " s ˆ3 V l `ř p R p S p,sl ˙`γ lg ˆ2 S lg `ř p R p Γ p sin θ ˙ « ´1 3 V ř m s `3 V m `4πR m 3 p1 ´cos δ m q ˘`γ lg `2 S m `4πR m 2 sin δ m sin θ ˘(19) p cap DEM " 1 3 V tr ˆř m f cap b l ˙" 1 3 V ř m tr ´ f cap b l ¯" 1 3 V ř m f cap . l " ´1 3 V ř m πR 1 sin δ 1 ps R 1 sin δ 1 `2 γ lg sinpθ `δ1 qq pR 1 `R2 `dq « ´1 3 V ř m πR m sin δ m ps R m sin δ m `2 γ lg sinpθ `δm qq p2 R m `dq (20)
In order to facilitate the comparison between p cap hom and p cap DEM , both Eqs. ( 19) and ( 20) are approximated, considering that all menisci m connect a pair of spheres of mean radius

R m « R 1 « R 2 .
Under such assumption, filling angles on both spheres are also equal, i.e.

δ 1 « δ 2 « δ m .
The two mean capillary stress expressions are different, with for instance the menisci volume V m and surface S m entering Eq. ( 19), but not (20). Nevertheless, the drying path under constant isotropic stresses (p " 20 kPa) shows a remarkable agreement between the DEM and the homogenization calculations as shown in Fig. 7. Such a comparison suggests 7. Mean capillary pressure as expressed from DEM or the homogenization approach: drying path on an isotropic sample that both expressions ( 19) and ( 20) are in fact equivalent. Due to intricacies in suction and surface tension terms together with the liquid bridge geometry -through Laplace-Young's equation -in Eqs. ( 19) and ( 20), the equivalence of the two equations cannot be ruled out based simply on the apparent differences between them.

General case

We now conduct a more general comparison, considering simple shear loading as shown in Fig. 8 where principal stresses rotate and capillary stress tensors involve full components.

Starting with an isotropic distribution of contact normals in a DEM sample, an initial wetting x y z FIG. 8. Simple shear loadings: Bv x {By " cst, Σ yy " Σ zz " Σ lat " cst " 20 kPa is applied so that the fluid phase disperses through it in the form of liquid bridges computed only at contacting particles, as previously discussed. Upon mechanical loading, whenever particles get too far away after contact is lost, liquid bridges will inevitably rupture. This is captured in the DEM simulations when no physical solution to Laplace-Young equation for a meniscus between a pair of particles can be found, while considering a constant and uniform capillary pressure throughout the sample. On the other hand, new menisci form at every new contact as dictated by the Laplace-Young equation. These changes in liquid phase distribution occur under fully drained conditions, without any consideration of pore flow within the sample, as opposed to [START_REF] Gili | Microstructural deformation mechanisms of unsaturated granular soils[END_REF], for instance. Then, the induced anisotropy of contact normals that builds up during loading is similarly reflected in the liquid bridge distribution. Accordingly, such anisotropy in the fluid phase distribution necessarily leads to a non-zero deviatoric part for the capillary stresses, s cap " σ cap ´pcap δ. The capillary stress tensors as obtained from homogenization and DEM calculations -Eq. ( 12) or (18) -are then compared in terms of mean pressure p cap , deviatoric stress |q cap | " a 3{2 ||s cap ||, and Lode angle ϑ cap P r0; 60 0 s. The orientation ψ cap of the minor principal capillary stress (the greatest in absolute value) with respect to the x-axis as depicted in Fig. 8 is also examined. Principal directions orientation (the curves are identical for θ " 40 0 and not plotted for clarity purposes) FIG. 9. Capillary stresses as expressed from DEM or the homogenization approach during simple shear loading. For s = 10 kPa; and θ " 0 0 (S r = 10 ˘1 %) or θ " 40 0 (S r = 3.7 ˘0.5 %)

Uniqueness of micromechanical description of the pendular regime

The various comparisons of the stress computations based on DEM and homogenization approaches presented in the above all point to the same stress description at the REV scale, neglecting the gas pressure, so that:

Σ " 1 V ÿ DE ż V DE σ dV " 1 V ˜żVs σ s dV `żV l σ l dV `żVg σ g dV `żS lg π lg dS ¸(21)
the equivalence between the two calculation approaches is remarkable since the DEM model does not explicitly include the internal forces arising from fluid phases V l , V g , and the liquidgas interface S lg . Instead, the distributed actions of these phases enter the DEM model only through the resultant forces they exert on each DE. With such a simplification, it turns out that this procedure is still appropriate for describing properly the mechanical contributions from the fluid phases and interfaces when deriving stresses in the DEM model.

Indeed, stresses computed from resultant point forces applied to the each DE turn out to coincide with correct stresses expected in an unsaturated medium where all phases are explicitly accounted for. It seems that a particular DE in an unsaturated DEM model must be considered as an entity made up of a solid particle and half of a meniscus. Indeed, we recall that the stress associated to one discrete element ş V DE σ dV cannot be equal to the stress calculated for one solid particle ş Vp σ s dV . This is because the sole consideration of resultant point forces is not appropriate to analyze stresses. We finally note that if the net resultant forces on every DE is zero such as for the saturated case, the equivalence between the DEM and a homogenization approach for stress description would be lost.

This comparison of two modeling approaches does have practical implications in the realm of laboratory experimental testing. Similar to the numerical analysis described in this paper, the capillary stresses could be alternatively determined from actual physical samples in which the fluid phase distribution is measured experimentally using imaging techniques [START_REF] Willson | Quantification of grain, pore, and fluid microstructure of unsaturated sand from x-ray computed tomography images[END_REF]. Because of the equivalence between the DEM and the homogenization approach, it would then be possible to avoid measuring comprehensive details of the fluid phase structure as needed in the homogenization approach, and judiciously use fewer parameters of the wetted surfaces as evoked in the DEM. For instance, the measurements of filling angles δ and contact angle θ would be enough to compute the mean capillary pressure considering Eq. ( 20), without the need of additional measures such as the volume and interface area of the liquid phase, V m and S m , if Eq. ( 19) were to be considered.

Capillary stress descriptions

The capillary stresses obtained from micromechanics -either the DEM or the homogenization approach since they are equivalent -are now quantitatively compared with classical capillary stresses descriptions from Bishop's equation with χ " S r :

σ cap " ´s S r δ (22)
Even though [START_REF] Bishop | Some aspects of effective stress in saturated and partly saturated soils[END_REF] himself identified the difficulties in interpreting experimental data with Eq. ( 22), this equation has been thermodynamically verified by several authors [START_REF] Houlsby | The work input to an unsaturated granular material[END_REF][START_REF] Nuth | Effective stress concept in unsaturated soils: Clarification and validation of a unified framework[END_REF][START_REF] Lu | A closed-form equation for effective stress in unsaturated soil[END_REF][START_REF] Nikooee | Effective stress in unsaturated soils: A thermodynamic approach based on the interfacial energy and hydromechanical coupling[END_REF]. Contrary to the micromechanics here presented, such thermodynamic considerations adopt an isotropic mechanical description with scalar stresses and strains. Also, additional hypotheses are usually involved. In particular, using the "suction stress" designation, [START_REF] Lu | A closed-form equation for effective stress in unsaturated soil[END_REF] neglected the energy contribution due to the interface to demonstrate Eq. ( 22). It is noteworthy that [START_REF] Lu | A closed-form equation for effective stress in unsaturated soil[END_REF] considered for S r the effective saturation that disregards the bound residual water layers, as it is done here. Considering the same drying paths as in Fig. 7, Eq. ( 22) is then compared to the capillary stresses obtained by micromechanics (Fig. 10).

Here, Eq. ( 22) underestimates drastically the capillary stresses obtained from the micromechanical approaches. Moreover, qualitative differences clearly exist at very low water contents. Indeed, the micromechanics-based capillary stress has a non-zero finite limit when S r tends asymptotically towards zero, contrary to the capillary stress obtained from Eq. ( 22). While such Bishop-like capillary stress expression could be used to interpret experimental results for various soils with significant water content [START_REF] Lu | A closed-form equation for effective stress in unsaturated soil[END_REF], the underlying omission of the interfaces invalidates the application of Eq. ( 22) to the pendular regime. [START_REF] Chateau | Micromechanics of saturated and unsaturated porous media[END_REF] also concluded through theory the invalidity of Bishop's stress in cases where interfaces are significant.

Finally, we again highlight the deviatoric nature of the capillary stresses in the general 10. Mean capillary pressure from micromechanics or from p cap " ´s S r . Both approaches predict zero capillary pressure for strictly zero saturation (dry conditions) S r " 0, which is not included in the data case, with non zero |q cap | (Fig. 9). Such deviatoric nature cannot be accounted for in Eq. ( 22).

For instance, when θ " 40 0 , |q cap {p cap | reaches 0.25 due to the induced anisotropy during the simple shear loading considered in Fig. 9. As a crude interpretation of such value, it is to be noted that a sand with a hypothetical friction angle of 7 0 would fail under a triaxial loading with q{p " 0.25.

CONCLUSIONS

We highlight the importance of a membrane stress as a micromechanical feature governing the stress in wet granular media. In this connection, two calculation approaches based on DEM and an analytical homogenization method have been presented to derive the same continuum stress models of an unsaturated granular soil in the pendular regime. The consistency between the DEM-and homogenization-stress calculations is remarkable, considering that DEM does not explicitly account for the fluid and interface phases of the unsaturated soil. Rather, the actions of the latter phases are transmitted solely through resultant forces to the solid particles. While resultant point forces involved in DEM calculations are adequately tied to micro-kinematics such as displacements of solid particles and hence strains in a granular soil, they may not properly describe the statics, i.e. stresses. In spite of this, it appears that the actions of all the concerned phases (solid, liquid, gaseous and interface) can be considered, one way or another, by the sole application of resultant forces on particles to give the correct stress description in the discrete element modeling of unsaturated conditions.

As seen in the preceding discussions, the above interesting numerical result has experimental implications such as the facilitation of the stress state measurement in unsaturated granular soils based on fluid phase distribution parameters determination by imaging techniques.

Turning to capillary stresses that describe the actions of the fluid mixture, we argue that a precise description of such capillary stresses in the pendular regime requires a comprehensive knowledge of the microstructural details. In particular, because of possible preferred orientations of wetted surfaces, the capillary stresses are deviatoric (anisotropic) in the general case. Using an averaged fluid pressure that is isotropic in nature to describe these capillary stresses is then conceivably inadequate.

Further extensions of this work are possible by considering the capillarity along solid-fluid interfaces to broaden the analysis of interfacial tension. Also, further analytical demonstration of the equivalence between the two models is currently being worked out.

  FIG. 1. Schematized REV of a pendular unsaturated soil

  FIG. 4. SWCCs of the DEM model in case of perfect wetting (θ " 0 0 ). Ottawa sand data from[START_REF] Willson | Quantification of grain, pore, and fluid microstructure of unsaturated sand from x-ray computed tomography images[END_REF]

  FIG. 6. Liquid-gas interfacial specific area

Fig. 9

 9 Fig.9shows a remarkable agreement between the capillary stress tensors obtained respectively from the DEM and homogenization calculations. Thus, this numerical result suggests that the two approaches are definitely equivalent. The next section delves into a more detailed discussion of these results.

TABLE 1 .

 1 DEM model parameters

	Parameter	Value
	Y (MPa)	50
	P (-)	0.5
	ϕ ( 0 )	30
	D max D min D 50 (mm)	3 0.05
	γ lg (N/m)	0.073
	θ ( 0 )	P r0; 60s
	volumes, follow immediately. For illustrative purposes, the Soil Water Characteristic Curve

ACKNOWLEDGEMENTS

This work was funded by the Natural Sciences and Engineering Research Council of Canada and Foundation CMG. The authors gratefully acknowledge rich discussions with Rakulan Sivanesapillai (Institute of Mechanics, Ruhr-Universität Bochum, Germany) and Félix Darve (3SR, Grenoble Universités, France). They also extend a warm thank-you to Mahdad Eghbalian, a Doctoral student within the authors' research group, for his attempts to demonstrate analytically the equivalence between Eqs. (19) and (20) in the framework of the toroidal approximation.


A Shah

J Coste

J J Lemaire

E Taub

R Guzman

E Schkommodau

S Hemm-Ode

Optimization of deep brain stimulation surgery for Parkinson's disease with quantitative rigidity evaluation

Background

Deep brain stimulation (DBS) is now a widely accepted surgical treatment for Parkinson's disease (PD). Electrodes are implanted in the patient's brain after intraoperative test stimulation. Changes in parkinsonian rigidity during test stimulation are detected by an evaluator, usually a neurologist, by identifying changes in the resistance of the patient's arm to a passive movement. We hypothesised that at the moment of reduction in rigidity, the speed with which the evaluator moves the patient's arms increases and that this change and its amplitude can be detected with an acceleration sensor. The aim of the present study was to test this hypothesis by collecting data during DBS surgery. Furthermore, to know more about the optimal stimulation target, these quantitative data were categorized based on the anatomical location of the electrode during test stimulation.

Discussion

 The additional acceleration measurements during the surgery did not increase operation time or the patient's discomfort.  Higher sensitivity when using the accelerometer recording system; effective stimulation amplitudes were found for 33 additional test stimulations.  Conventionally targeted STN requires the lowest stimulation amplitude to reduce rigidity, but has significantly higher chances of side effect occurrence. The Fields of Forel have slightly higher stimulation amplitudes but have much lower change of causing side effects.  Sufficient baseline data is necessary for proper identification of BQAs.  There is an inherent subjective component in the acceleration analysis because the evaluation is done by the neurologist.

Conclusion

• Changes in rigidity of PD patients can be quantified during passive movements by measuring data from the evaluator.

• Acceleration measurements confirm the subjective evaluation, but they seem to be more sensitive (Fig 4).

• STN may not be the most efficacious target structure. The patient may benefit from an electrode placed closer to the Fields of Forel.

Acknowledgements

This research has been supported by the Swiss National Science Foundation (SNSF) and the Germaine de Stael program of the Swiss Academy of Engineering Sciences (SATW).

  Standard Deviation Signal Energy Peak Frequency Amplitude

Figure 2 .Figure 3 .Figure 1 .

 231 Figure2. A 3-axis accelerometer housed in a 3D printed plastic case is mounted on the evaluators wrist and connected to a laptop with in-house developed recording software. This software is also connected to electrophysiology system which records test stimulation amplitudes.

Figure 4 :

 4 Figure 4: Box plot showing comparison between BCA and BQA

  Out of the 188 test stimulations, 138 evaluations were used for comparison between BCA and BQA. For 14 evaluations none of the thresholds were found, for 30 evaluations no BCA were found and for 6 evaluations no BQA were found.  Results of Wilcoxon 2 sided rank test showed that BQAs were significantly lower than BCAs(p<0.001, Fig4).  The 138 evaluations were distributed in 5 structures: 2 in Substantia Nigra, 3 in Thalamus, 27 in FF, 26 in ZI and 80 in STN.  STN had the lower average values for BCA and BQA, but highest occurrences of side-effects. The comparison with FF and ZI can be seen in Figure5. Effective stimulation amplitudes inducing a reduction in rigidity were identified (best quantitative amplitude, BQA, Fig3).  BQA and BCA were compared using Wilcoxon 2 sided signed rank test.  Data were grouped based on the anatomical structure in which the electrode was present during test stimulation for a particular position.  Average BCA, Average BQA and number of side effect occurrences were compared for the STN, zona incerta (ZI) and Fields of Forel (FF).

	Results  Clinical study (University Hospital in Clermont-Ferrand): 9 rigidity patients undergoing DBS surgery  Preoperative manual outlining of sub-thalamic nucleus (STN) and its anatomic neighbors was done using iPlan (Brainlab, Feldkirchen, Germany, Fig 1A and B) and the target point was selected.  Intraoperative microelectrode recording (MER) and test stimulations (Fig 1C and Fig 1D) were performed.  Maximum reduction in rigidity during passive movements and the corresponding amplitude (best clinical amplitude, BCA) were noted for all positions.  188 test stimulations in total.  One anatomical structure was attributed to every test stimulation position [1].  Acceleration data were recorded and evaluated for all test stimulations (Fig 2).  Data filtering and statistical feature extraction preceded normalization of features to baseline recordings.  Methods 	A	B	C		D
					A
					B
				Dr. Simone Hemm-Ode	Ashesh Shah	Dr. Jerome Coste
			Contact:	T +41 61 4674-796 F +41 61 4674-701	T +41 61 4674-413 ashesh.shah@fhnw.ch	T +33 4 73 751 002 jerome.coste@udamail.fr
				simone.hemm@fhnw.ch