William Chèvremont

Hugues Bodiguel
email: hugues.bodiguel@univ-grenoble-alpes.fr

Bruno Chareyre

Lubricated contact model for numerical simulations of suspensions

Keywords: suspension, granular, contact, lubrication, friction

Discrete granular models are a natural choice when simulating dense suspensions, when the small distances between the particles lead to dominant contributions by lubrication and contact forces. In such case one can get rid of the costly resolution of Navier-Stokes equations, using closed form expressions for lubrication terms. However, those terms diverge when two hard spheres approach contact, and there are issues when integrating them directly with the finite precision of floating point calculations. In this paper, we introduce a visco-elasto-plastic interaction model for suspended spheres, which combines lubrication and elastic-frictional contact behaviour depending on surface roughness. An integration scheme is proposed for that model. Unlike earlier methods, the scheme enables an unconditionally stable time-integration of the interactions.

The case of perfectly smooth spheres (null roughness), namely, is integrated correctly.

The theoretical results are well reproduced in benchmark tests on two-sphere systems:

one sphere sedimenting on one other and two spheres in a shear flow. From these benchmark tests, we propose phase diagrams showing the interplay between viscosity, roughness and stiffness. The second test case highlights the origin of non-reversibility particle trajectories. It is controlled by the particle roughness for rigid particles, and by the particle deformation when the capillary number is higher than the relative roughness.

Introduction

Suspensions of solid particles dispersed in a viscous fluid are ubiquitous in natural, industrial as well as biological flows. For dense suspensions, the resistance to flow depends on a combination of frictional contacts between the particles and viscous interactions mediated by the interstitial fluid. Those combined effects might result in non-Newtonian behaviour even for simple cases such as suspensions of monodisperse spheres in a Newtonian fluid [START_REF] Stickel | Fluid mechanics and rheology of dense suspensions[END_REF][START_REF] Morris | A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow[END_REF]. Particle-scale simulation with discrete element methods (DEM) is a way to investigate the intricate micro-scale processes at play in such flow. There are a few caveats in the concrete implementation of a lubricated contact model in DEM though, and coverage of this point in the literature is scarse whereas many other types of interaction model have been analysed in details (see, namely, [START_REF] Thornton | A theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres[END_REF][START_REF] Thornton | Granular dynamics, contact mechanics and particle system simulations[END_REF] and other papers by the same author). This paper contributes a robust method.

Several simulation frameworks have yet been developed in relation to suspensions. Direct Numerical Simulations (DNS) using some Navier-Stokes solvers are feasible though computationally demanding. They require a high spatial resolution to resolve the flow between close particles. Regardless of spatial resolution, though, some singularities will not be captured. Those singularities are associated to the divergence at the approach of contact of the so-called lubrication forces, i.e. the viscous resistance to relative motion between two immersed particles. The change of distance between two rigid spheres, namely, leads to a resisting normal viscous force which scales like un u -1 n if u n is the gap between the solid surfaces. Since getting this resistance from DNS would require to shrink the spatial resolution to virtually 0 when u n → 0, the models generally incorporate lubrication corrections to include what's beyhond mesh resolution. The approaches coupling DNS with suspended particles include the Force-Coupling Method (FCM) [START_REF] Peters | Rheology of non-brownian suspensions of rough frictional particles under shear reversal: A numerical study[END_REF][START_REF] Yeo | Simulation of concentrated suspensions using the forcecoupling method[END_REF], Fictious Domain method (FD) [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF] or Smooth Particle Hydrodynamics (SPH) [START_REF] Vázquez-Quesada | Rheology and microstructure of non-colloidal suspensions under shear studied with smoothed particle hydrodynamics[END_REF][START_REF] Tanaka | Simulation method of colloidal suspensions with hydrodynamic interactions: Fluid particle dynamics[END_REF].

While the aforementioned methods spend most of the CPU time in resolving the fluid, multiple authors found that the steady-state flow of suspensions could be simulated accurately by leaving the fluid unresolved and introducing only the lubrication terms through closed form expressions. Down to solid fractions of the order of φ = 0.2, analytical solutions for regular arrays of spheres [START_REF] Van Den Brule | Modeling of concentrated suspensions[END_REF] as well as numerical solutions in three-dimensional disordered systems [START_REF] Ball | A simulation technique for many spheres in quasi-static motion under frame-invariant pair drag and brownian forces[END_REF][START_REF] Mari | Shear thickening, frictionless and frictional rheologies in non-brownian suspensions[END_REF][START_REF] Marzougui | Microscopic origin of shear stress in dense fluid-grain mixtures[END_REF][START_REF] Gallier | Simulations of sheared dense noncolloidal suspensions: Evaluation of the role of long-range hydrodynamics[END_REF][START_REF] Chèvremont | Quantitative study of the rheology of frictional suspensions: Influence of friction coefficient in a large range of viscous numbers[END_REF], provide rather accurate predictions of bulk viscosity by taking lubrication as the only form of viscous dissipation.

As lubrication forces decrease sharply with distance the pairs farther than one particle diameter, typically, are insignificant and can be safely ignored. It of course leads to much greater computational efficiency. Worth noting, good agreement with experimental data needs to account not only for the repulsion-attraction effects but also for the viscous response to contact shearing [START_REF] Marzougui | Microscopic origin of shear stress in dense fluid-grain mixtures[END_REF] and for contact friction [START_REF] Chèvremont | Quantitative study of the rheology of frictional suspensions: Influence of friction coefficient in a large range of viscous numbers[END_REF].

The lubrication effects being pair-wise interactions they fit well in conventional discrete element methods (DEM), which track the motion of interacting particles with an explicit integration scheme [START_REF] Marzougui | Microscopic origin of shear stress in dense fluid-grain mixtures[END_REF], and which was adopted for this work.

The aforementioned models share the same closed-form solutions for the lubrication terms overall. However they differ in the way lubrication and solid contact behaviour are combined. It must be noted that in most DEM models a contact is defined when u n < 0. This situation, which corresponds to a small overlap between the geometrical spheres, reflects a deformation near the contact region, following Hertzian models or some linearized form of them. u n < 0 comes in contradiction with the fact that the lubrication forces diverge when u n = 0, which would in principle prevent contacts. A classical argument to resolve the contradiction, supported by empirical facts [START_REF] Zhao | Interaction of two touching spheres in a viscous fluid[END_REF], is that asperities of the solid surfaces can be in contact when there is still a fluid film of finite thickness between the surfaces. Along this line a surface roughness εa (length scale of the asperities, a being the mean sphere radius and ε relative roughness) is introduced in the models such that the normal lubrication force scale like (u n + εa) -1 [START_REF] Brändle De Motta | Numerical modelling of finite-size particle collisions in a viscous fluid[END_REF][START_REF] Seto | Discontinuous Shear Thickening of Frictional Hard-Sphere Suspensions[END_REF][START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF][START_REF] Marzougui | Microscopic origin of shear stress in dense fluid-grain mixtures[END_REF],

and thus it takes a finite value when a contact is created (i.e. when u n = 0). Various implementations of this idea have been proposed. Some authors proposed that whenever u n ≤ 0 the viscous multiplier be kept constant [START_REF] Brändle De Motta | Numerical modelling of finite-size particle collisions in a viscous fluid[END_REF], or that all lubrication terms be dropped [START_REF] Brändle De Motta | Numerical modelling of finite-size particle collisions in a viscous fluid[END_REF][START_REF] Seto | Discontinuous Shear Thickening of Frictional Hard-Sphere Suspensions[END_REF][START_REF] Gallier | Simulations of sheared dense noncolloidal suspensions: Evaluation of the role of long-range hydrodynamics[END_REF] and the interaction model replaced by a linear springdashpot model (LSD [START_REF] Walton | Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks[END_REF], conventional in dry granular models), tuned in order to produce a specific coefficient of restitution. As pointed out by some authors, though, the restitution coefficient as well as the distance at which lubrication is turned off are additional parameters and the results depend on tuning them relatively arbitrarily [START_REF] Gallier | Simulations of sheared dense noncolloidal suspensions: Evaluation of the role of long-range hydrodynamics[END_REF].

Finally, some authors considered that solid contact and lubrication could act simultaneously and be combined in a visco-elasto-plastic (VEP) model. Therein interaction force at a lubricated contact combines an elastic force in response to some surface deflection and a viscous force in response to changes of gap distance. A simple method to define changes of the gap distance is to set it equal to the total deflection (u n + εa) even if u n < 0 [START_REF] Trulsson | Transition from the viscous to inertial regime in dense suspensions[END_REF], as if soft asperities would deform at the surface of otherwise rigid spheres. Another method [START_REF] Marzougui | Microscopic origin of shear stress in dense fluid-grain mixtures[END_REF][START_REF] Chèvremont | Quantitative study of the rheology of frictional suspensions: Influence of friction coefficient in a large range of viscous numbers[END_REF], which we follow overall, is to let the total contact deflection be split in two parts. One part is the change of gap, coming with the flattening of asperities and producing lubrication effects. The other part is an elastic length u e corresponding to change of shape around the contact (Hertzian scale) under the combined action of contact forces through asperities and lubricating pressure in the gap.

The gap is then u n + εau e . This model recovers others as special cases if spheres and asperities have distinct stiffness: a low sphere to asperity ratio corresponds to the hard sphere limit, a high ratio corresponds to the hard asperity limit (constant gap).

A potential downside of the VEP models is that the smooth limit (ε → 0 + ) may be difficult to approach. In the hard sphere limit the interaction forces diverge when (u n + εa) → 0 + (with u n < 0 at contacts), which sets an implicit lower bound to ε. In addition, viscosity dominated regimes when (u n + εa) is small though strictly positive, may cause instabilities with explicit time integration schemes. This is partly regularized by introducing u e , but then the stiff problem becomes the time integration of u e in itself. As a matter of fact ε < 10 -2 remained untouched overall with VEP models while other methods were applied successfully to ε = 10 -3 [START_REF] Seto | Discontinuous Shear Thickening of Frictional Hard-Sphere Suspensions[END_REF] or ε = 0 [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF].

In this paper a model and a time integration method are proposed to alleviate the limitations of previous VEP approaches; namely • the physical model is free of arbitrary tricks in combining lubrication and solid contact;

• the limits ε → 0 and/or (u n + εa) → 0 is approached robustly.

The proposed model is based on the assumption that the elastic deformations comes from both the compliance of individual contact asperities and the global compliance of the solid surface. Robustness is obtained using an implicit time integration of the interaction forces in an otherwise explicit time-stepping algorithm.

The governing equations and the time integration scheme are given in the first part of the paper. In the second part two test cases are presented to assess the robustness of the approach and to identify different regimes by parametric analysis. The first test case is a suspended sphere subjected to gravitational load and bouncing on another, fixed, sphere; the other involves the interaction of two spheres in a sheared fluid. In both cases theoretical solutions exist and are recovered by our model.

Governing equations and time integration

The interaction between two immersed spheres is defined as the sum of two contributions from, respectively, direct (repulsive) contact between the solid surfaces and lubrication by the suspending fluid. It is assumed that they both contribute to deform the spheres elastically and that the deformation occurs through, either, a flattening of the surface asperities leading to a change of the interstitial gap u, or a more general deflection of the surface associated to a local change of shape, as if the sphere radii where modified locally by a small distance u e . For simplicity linear elasticity is assumed for both modes of deformation in what follows. It will lead to a visco-elastic model similar to the standard linear solid model, yet with variable viscosity. The approach could be extended to non-linear contact models, such as Hertzian ones, without major difficulty.

In the following, we detail the governing equations of the model. Since time integration only requires the calculations of interaction forces and torques between the two spheres at a given time step and is furthermore not the focus of the present work, we start by presenting these forces and torques. Since they involves intermediate quantities, and in particular the interstitial gap u between the (deformed) surface, their evaluation is not straightforward. In section 2.1, we first detail the normal component of the forces and explain how to efficiently compute the interstitial gap using an implicit scheme. Then, we detail the tangential components in section 2.2. These are more direct as there is no direct coupling with u contrary to the case of normal ones. However, several cases need to be distinguished, depending on the type of contact (no contact, sticking contact and slipping contact). Time integration is finally described briefly in section 2.3, but is beyond the scope of this paper. It is assumed that a finite repulsive contact force arises whenever the interstitial gap is less than the size of asperities. For two particles of mean radius a = (r 1 + r 2 )/2 the normal component of the contact force is taken proportional to the deflection of the asperities, following

F c n = -k n max(0, εa -u)n (1)
where εa defines the characteristic size of the asperities, and ε is termed relative roughness. k n is the normal stiffness of the asperities, and u is the distance between the solid surfaces (not counting the size of asperities). n is the unit normal of the contact (see figure 1). Here and in the following, subscripts n and s denote normal and shear components of the interactions, respectively. Superscript c and l will be used to distinguish contact and lubrication contributions.

In the elastic regime the tangential component of the contact force is an incrementally linear function of the tangential displacement. The contact model includes a threshold on the magnitude of the shear force following Coulombian friction, as in Cundall and Strack [START_REF] Cundall | Discrete numerical-model for granular assemblies[END_REF]. The shear force thus has to satisfy the inequality:

F c s ≤ µ m F c n , (2)
where µ m is the coefficient of contact friction.

Lubrication

The lubrication terms are defined as in Marzougui et al. [START_REF] Marzougui | Microscopic origin of shear stress in dense fluid-grain mixtures[END_REF]. They define the viscous resistance to relative motion between particles.

F l n = 3 2 πη f a 2 u u n (3)
F l s = πη f 2 -2a + (2a + u) ln 2a + u u v (4)
T l r = πη f a 3 3 2 + 63 500 u a ln a u ω × n (5)
T l t = πη f a 2 u ln a u (ω • n)n (6)
where F l n , F l s , T l r and T l t are respectively normal lubrication force, tangential lubrication force, rolling lubrication torque and twisting lubrication torque. v is the tangential displacement, ω the spin velocity, and η f the fluid velocity. The expressions are based on Frankel and Acrivos [START_REF] Frankel | On the viscosity of a concentrated suspension of solid spheres[END_REF] for the forces and Jeffrey and Onishi [START_REF] Jeffrey | The forces and couples acting on two nearly touching spheres in low-reynolds-number flow[END_REF] for the torques.

Although the above expressions define the lubrication terms for virtually any distance u, they are exact only asymptotically in the limit u/a → 0. In practice, they are computed only for particle pairs within a distance range. The torques (equations [START_REF] Peters | Rheology of non-brownian suspensions of rough frictional particles under shear reversal: A numerical study[END_REF][START_REF] Yeo | Simulation of concentrated suspensions using the forcecoupling method[END_REF] are computed as long as u ≤ a. For u > a the direction of the torques switches to the same direction as the relative rotation and thus physical consistency is lost -which justifies the upper bound unambiguously. There normal and shear forces (equations 3-4) are also cut off beyond a certain distance yet there is no simple argument to fix the maximum distance in that case. The cut-off distance is left as a model parameter for now. It will be further examined on the basis of the second test case.

For writing simplicity, prefactors ν n , ν t (u) are introduced in the expression of forces, such that: u n depends only on the particles positions which, in concrete cases, will result from integrating the equations of motion for the particles (which is done independently of the present derivation). A change of the center-to-center distance relative to the stress free configuration implies that either the distance between the surfaces is changing, or the surfaces are deforming, or both. In all cases u n accumulates the other two displacements and the following equality holds:

ν n = 3 2 πη f (7) ν t (u) = πη f 2 -2a + (2a + u) ln 2a + u u (8)
F l n = ν n a 2 u u n (9)
F l s = ν t (u)v (10
u n = u + u e (11)
It is assumed that the lubricating pressure and the contact forces through asperities both contribute to deflect the surfaces by a distance u e in the normal direction, and that the compliance is defined by a single coefficient of stiffness k b such that the total normal force F n = F l n + F c n and u e are linearly dependent:

F n = k b u e n (12)
Combining equations 1, 3 and 12 and projecting in the normal direction n yields the scalar form

k b u e = -k n max(0, εa -u) + ν n a 2 u u (13)
Elimination of u e using equation 11 yields a differential equation governing the visco-elastic coupling, where u is the only unknown (owing to the fact that u n results from motion integration). Normalization by the length a and the characteristic visco-elastic time τ = ηa/k b leads to the dimensionless form

(1 + α)u * -u * n -αε + u * u * = 0 (14) α(u * ) =            0 when u * > ε (no contact) k n k b when u * <= ε (contact) (15)
Here and in the following " * " is used to denote the dimensionless form of a variable (e.g. u * = u/a). In this paper, k n /k b is set equal to 1 at contacts; in other words the stiffness of asperities and the stiffness of the surfaces are taken equal. The numerical scheme is compatible with all positive values.

[Figure 1 about here.]

Tracking the evolution of the interaction with time asks for u integration using equation 14 while continuously updating u n to reflect particles motion. For this purpose, and in order to insure unconditional stability, an implicit backward Euler method is used. The time derivative of u * is approximated by

u * + = u * + -u * - ∆t * + O(∆t * ), ( 16
)
where ∆t * is the time step normalized by τ, and the exponents " -" and " + " refer to the times at which a quantity is evaluated (the start and the end of the time interval). Using equation 14, we obtain

u * + 2 (1 + α + ) -u * + u * n + α + ε - 1 ∆t * - u * - ∆t * = 0, (17)
which is a second order polynomial equation, of positive discriminant. The smallest solution is always negative, hence rejected since the gap u is positive. Finally the updated gap and the associated force are given by

u * + = 1 2(1 + α + )          u * n + α + ε - 1 ∆t * + u * n + α + ε - 1 ∆t * 2 + 4(1 + α + ) u * - ∆t *          (18)
F n = k b a u * n -u * + . (19)
Interestingly, equation 18 remains applicable for every positive value of u and it does not suffer from the singular evolution of lubrication as u → 0 + . This is a key property of the proposed approach.

The function α that is used to write the general solution introduces discontinuous derivatives at the transition between contact and no-contact solutions, and it has to be handled with care. Setting u * + = ε in equation 17 gives the time delay after which the transition occurs, i.e.

t * c = ε -u * - ε(ε -u n ) . ( 20
)
Where t * c is the dimensionless critical collision time. If t * c is positive and the dimensionless time step ∆t * > t * c , a contact transition occurs within the time step. In such a case, we integrate on a fraction of the nominal time-step, noted ∆t * , and starting from an intermediate solution (variables subscripted with i) which correspond to the transition:

                 u * - i = ε ∆t * i = ∆t * -t * c α + i = 1 -α + (21)
Implementing equation 18 directly would have detrimental effects in terms of accuracy when both u * and u * n are small. This is because with such settings the terms 1/∆t * dominate all other terms in equation 18, which leads to a substraction of two nearly equal terms, beyond the accuracy of floating point operations. The algorithm would effectively give u * ≤ 0 after some time if a contact pair is compressed permanently. The fact that u * appears often in a logarithm suggests the change of variable ξ = log(u * ) to circumvent this accuracy issue. The following equation on ξ now has to be solved:

(1 + α + ) exp(ξ + ) -u n + α + ε - 1 ∆t * - exp(ξ --ξ + ) ∆t * = 0 ( 22
)
This equation can be solved by conventional non-linear solvers, initialized at ξ + = ξ -.

The bisections method is used in the current implementation.

Tangential components

The tangential forces are combined in a visco-elastic model as shown in figure 2.

Like for the normal part, lubrication and contact act in parallel, and they both contribute v s is always known before computing the tangential interaction. Just like u n for the normal component it results from motion integration. The gap u can be considered known, too, since it was derived in the previous section independently of the tangential components.

[Figure 2 about here.]

The evolution of the shear components is governed by the system of equations

           v s = v e + v F s = -k t v e = min ||F c s ||, µ m ||F c n || F c s ||F c s || -ν t (u)v. ( 23
)
The equations express, respectively, the additivity of elastic and plastic displace-

The three cases are detailed below.

No contact

If the normal contact force is null, the shear contact force is also null as a consequence of Coulomb condition. Equations 23 degenerate to a Maxwell-type viscoelastic equation given by

-k t v e = -ν t (u)(v s -ve ). ( 25
)
Using equation 24 and rearranging yields to

v + e = -ν t (u)(v s ∆t + v - e ) ν t (u) + k t ∆t (26)
and to

F + s = ν t (u)(F - s -k t vs ∆t) ν t (u) + k t ∆t . (27)

Sticking contact

If the contact is sticking over a time interval, the elastic displacement rate is equal to total displacement rate: ve = vs (no slip). The incremental change of v e can be obtained by direct integration of vs and in this purely elastic regime the force changes linearly with respect to this change, i.e.

v + e = vs ∆t + v - e , (28)
F e+ s = F - s -k t vs ∆t.. (29)
The force obtained herein by assuming no-slip is a trial elastic force, which may or may not satisfy Coulomb condition. If it does, then we set simply F + s = F e+ s . If it does not then we proceed to the third case (slipping contact) to find a visco-elasto-plastic solution.

Slipping contact

The trial elastic force from Eqs. 29 defines the direction of the contact force and of the plastic slip if it occurs. An equivalent form of equations 23 is then

-k t v e = µ m ||F c n || F e+ s ||F e+ s || -ν t (u)(v s -ve ). ( 30
)
Using equation. 24, we obtain

v + e = µ m ||F c n || F e+ s ||F e+ s || ∆t -ν t (u) vs ∆t + v - e ν t (u) + k t ∆t , (31)
and

F + s = k t µ m ||F c n || F e+ s ||F e+ s || ∆t + ν t (u) F - s -k t vs ∆t ν t (u) + k t ∆t . ( 32
)
Further substitution and factorization by F e+ s gives a more compact expression:

F + s = F e+ s ν t (u) + k t ∆tµ m ||F c n || ||F e+ s || ν t (u) + k t ∆t . ( 33
)

Note on motion integration in DEM

The classical motion integration in DEM requires a single evaluation of the interaction forces at time t to compute positions and velocities at time t + ∆t. This is done classically with a symplectic integrator using a centered, second order accurate, finite difference approximation of the translational and rotational accelerations. The details of this integration is beyond the scope of the paper yet an important feature is that it sets an upper bound to ∆t for numerical stability, which has a critical influence on the total computation time (a detailed derivation can be found in [START_REF] Hosn | Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction[END_REF]).

In brief the computational time-step for rate-independent (non-viscous) interactions depends on the incremental stiffness of each interaction, i.e. the partial derivatives of F with respect the particle positions.

For rate-dependent interaction models such as the LSD however, an additional stability condition appears in the viscosity dominated regimes, since the forces depend not only on positions but also on velocities. To avoid this constraint many authors tune the interaction viscosity in order to stick to the under-damped, elasticity dominated, regime. Obviously such a trick is not acceptable in a lubrication model, where fluid viscosity has to be taken for what it is.

A key feature of the visco-elastic interaction defined in previous sections is that it leaves the stability of explicit motion integration independent of viscosity ν, and even of ν/u. Indeed the stiffness g is an upper bound of the incremental stiffness regardless of other parameters (the upper bound being approached when viscosity dominates).

In practice, it is thus enough to determine ∆t using the same expressions as for rateindependent interactions, using g as the contact stiffness. It does not only make selecting the computational time-step straightforward, but also, and more importantly, it lets limit cases be approached with no impact on the computational cost (e.g. viscosity dominated regimes, perfectly smooth particles, vanishing gap distances...).

Test cases

In this section, two 2-spheres simulations are performed using the computational scheme introduced previously and the DEM code Yade-DEM. In this code, motion is integrated by a conventional explicit, 2nd order, central finite difference scheme. A link to the algorithm used to generate these test cases is provided in reference 25.

Falling sphere

In this first test case two immersed spheres undergo normal motion only. The simulation consists in one fixed sphere and one free sphere subjected to gravitational acceleration. The free sphere is initially static at position (0, u 0 , 0). The fixed sphere remains at position (0, 0, 0). This test shows the robustness and numerical stability of the model down to very small gaps. Various regimes are expected depending on input parameters.

When the viscosity is high enough, the trajectory is dominated by viscous effects and the gap distance tends to decay exponentially with time. If roughness is finite, the gap closure converges to a final value corresponding to static equilibrium, else it converges to zero. For less viscous situations, the trajectory may appear nearly parabolic in the early stage. However, viscosity will necessarily come into play when the spheres are close enough. Finally, the free sphere may bounce, oscillate around the static position, or approach it in over-damped mode.

[Figure 3 about here.]

The other relevant number for this test case is the Stokes number, generally defined as St = aρu/η f . Since the characteristic velocity for this problem is the undamped free-fall velocity, u2 = 2u 0 g, the Stokes number reads

St = aρ 2u 0 g η f (39)
When the Stokes number is small enough, the trajectory is driven by viscous effects, as in over-damped oscillators. In contrast, if this number is high, the sphere should collide as in a nearly undamped, or under-damped, oscillator.

[Figure 4 about here.]

Figure 4 shows several typical solutions of this test case. The initial trajectory is always dominated by gravity since viscous forces are initially null (no initial velocity), which results in a quadratic evolution of position with time just after the start. Then, for rigid particles or smooth particles at low Stokes numbers, the movement is damped (a), and the falling sphere reaches static equilibrium smoothly. With soft particles, the trend is for the free sphere to be entrapped by lubrication (b), as if the surfaces where sticking.

The oscillations in u are so small that viscous dissipation becomes negligible, hence the center of mass keeps oscillating in a nearly undamped mode. For high Stokes number as in (c), the rebound energy is sufficient for the sphere to detach after the collisions and we recover damped bouncing. Note that the restitution coefficient progressively decreases at each bouncing, and the number of bouncing is finite. Situation (a) can be seen as the last rebound of a series starting like (c). Finally, (d) is intermediate between (b) and (c), the stiffness is small enough to enable oscillations but still high enough for bouncing.

[Figure 5 When the Stokes number is small enough, all trajectories are completely driven by viscous forces, and the solution calculated in (36) is recovered until the two spheres create a contact. Figure 6a shows it is indeed the case down to very small gaps for smooth particles. For rough particles, this solution is also obtained before reaching contact. Rough particles with negative equilibrium distances u eq exhibit an exponential behaviour with a longer characteristic time associated to the apparent weight discussed above. As shown in Figure 6c an excellent agreement is found between the gap reached at the end of the simulation and the calculated equilibrium one.

[Figure 6 about here.]

To conclude, this first test case enables to verify the implementation of the method by comparing the numercial solution to the analytical one at low Stokes number for both rough and smooth spheres. It also highlights that for the latter, it is robust down to arbitrary small gaps.

Two spheres in a shear flow

[Figure 7 about here.]

[Figure 8 about here.]

The second test case concerns the motion of two neutrally buoyant spheres suspended in a sheared fluid. This problem has been studied extensively by numerical methods and compared to closed-form solutions in earlier works [START_REF] Metzger | Irreversibility and chaos: Role of lubrication interactions in sheared suspensions[END_REF][START_REF] Yeo | Simulation of concentrated suspensions using the forcecoupling method[END_REF][START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF][START_REF] Kromkamp | Lattice boltzmann simulation of 2d and 3d non-brownian suspensions in couette flow[END_REF][START_REF] Gallier | Rheology of sheared suspensions of rough frictional particles[END_REF]. The closed-form solutions come from analytical integration of Stokes' equations in linear flow fields [START_REF] Batchelor | The hydrodynamic interaction of two small freelymoving spheres in a linear flow field[END_REF][START_REF] Van De Ven | The microrheology of colloidal dispersions: Iv. pairs of interacting spheres in shear flow[END_REF], of which shear flow is a particular case [START_REF] Da Cunha | Shear-induced dispersion in a dilute suspension of rough spheres[END_REF][START_REF] Lin | Slow motion of two spheres in a shear field[END_REF].

The spheres are subjected to entrainment by the suspending fluid, and lubrication forces as they approach each other. The entrainment is modelized by Stokes drag, for which force F d and torque T d are given by

F d = -6πη f a ṙ -w f (40)
T d = -8πη f a 3 ω - 1 2 ∇ × w f (41)
where w f is the background fluid velocity given by (γy; 0; 0), ṙ and ω are the particle velocity and spin velocity, respectively. The initial particle positions are (0, 0, 0) and (-10a, d 0 sin θ, d 0 cos θ), with initial velocity and spin corresponding to the force and torque free case in the above equations. All distances are normalized by a and noted with " * " for the comparison with analytical results. Besides θ and d 0 /a, the relevant dimensionless numbers of the problem are the roughness and the ratio of viscous stress to stiffness: η f γ/E. This last number controls the deformation of the sphere induced by viscous forces. Examples of soft systems for which this coupling has been studied include bubbles, drops, vesicles or elastic capsules, for which the stiffness is an interfacial property (surface tension or surface modulus). By analogy to these cases and similarly to other authors (see e.g. reference 33), let us call the above number the Capillary Number, Ca = η f γ/E.

Whereas suspended smooth hard spheres follow reversible paths, both roughness (i.e. the possibility of solid contact) and deformability of the particles are known to induce non reversible, the particles being driven to different streamlines as a result of the interaction [START_REF] Da Cunha | Shear-induced dispersion in a dilute suspension of rough spheres[END_REF][START_REF] Lac | Hydrodynamic interaction between two identical capsules in simple shear flow[END_REF]. Figure 8 illustrates this three typical cases: smooth hard spheres, rough hard shperes, deformable particles.

[Figure 9 about here.]

The case of smooth hard spheres is particularly interesting since it can be directly compared to the analytical solution of Da Cunha and Hinch [START_REF] Da Cunha | Shear-induced dispersion in a dilute suspension of rough spheres[END_REF]. As it can be seen on figure 8a, the obtained trajectories follow the analytical solutions from Da Cunha and Hinch [START_REF] Da Cunha | Shear-induced dispersion in a dilute suspension of rough spheres[END_REF] for smooth hard spheres with a maximum deviation of 2%. However, looking into details to the gap between the spheres, a discrepancy could be detected between the numerical and the analytical one. Although it has no consequences on the trajectories, it is interesting to discuss its origin. In the numerical model, the best match to the analytical solution of Da Cunha and Hinch [START_REF] Da Cunha | Shear-induced dispersion in a dilute suspension of rough spheres[END_REF] being obtained at 4a. This cut-off value was thus retained since there was no other objective way to define it, and the minimal distance was somehow sensitive to it. That the same cutoff should be used in all cases is no implied since the two-sphere system is not representative of typical applications. The influence of the cut-off distance is much smaller when simulating dense or, even, semi-dilute suspensions. The dominant lubrication terms are indeed due to the nearest neighbours in general, and the nearest neighbours are generally much closer than 4a. Any cut-off value beyond one particle diameter will lead to quantitatively similar results.

The very good agreement between with the analytical solution for smooth hard sphere is interesting as it clearly shows that lubrication forces and torques are sufficient to describe the viscous interaction between two spheres. This conclusion differs to the one which is mentioned in Metzger et al. [START_REF] Metzger | Irreversibility and chaos: Role of lubrication interactions in sheared suspensions[END_REF], where a significant discrepancy is found between the numerical solution obtained with lubrication interaction and the analytical solution. This has been interpreted as the consequence of neglecting the long-range interaction which are present the Stokes equation. In fact, in reference 26, only the normal component of the lubrication force was considered; additionally no drag torque was considered. By turning off tangential forces or sphere rotation, we indeed obtain a similar discrepancy with the analytical solution , as shown in figure 10.

This consideration allows us to highlight that tangential lubrication forces are crucial to get quantitative results -they are too often neglected in numerical models. The cost of computing the shear component in the present model is small compared to the normal component and it does not introduce additional constraints in terms of numerical stability. Therefore there is no substantial computational gain to expect by neglecting these terms -besides the fact that ignoring them would make the model quantitatively wrong.

[Figure 10 about here.]

A set of typical trajectories for rough stiff spheres is presented in figure 8b. As expected, roughness breaks the symmetry as soon as contact occurs. The role of roughness is apparent in the plots of u * , which look as if the smallest distances were truncated. The sphere then moves to a streamline corresponding to min(u * ) ≈ ε irrespective of the initial position of the particle. As seen in figure 9 for smooth spheres (ε = 0), u * reaches a finite minimum value when d * 0 → 0, and this distance is about 3 10 -5 . Consistently, the results are not modified by a finite roughness as long as ε < 3 10 -5 .

Another set of trajectories is presented for soft spheres in figure 8c. Up to our knowledge, there is no analytical solution in that case. However, changes of streamline have been reported in the literature, based on experiments and on simulations [START_REF] Lac | Hydrodynamic interaction between two identical capsules in simple shear flow[END_REF]. The trajectories with large capillary number are qualitatively similar to those found in the literature. A noticeable effect is that the closure of the gap is smaller as compared to stiff spheres, which implies -for a given value of roughness -a smaller normal contact force and therefore a lower resistance to contact sliding.

[Figure 11 about here.]

We have performed a systematic analysis of the streamline deviation, varying both roughness and capillary numbers. The results are summarized in figure 11, which reports the magnitude of the streamline deviation for almost coaxial initial distances (d * whereas when Ca ε, it is controlled by the particle deformation. Strikingly, when Ca and ε are of the same order of magnitude, the streamline deviation is smaller than the case of rigid particles of similar roughness or than the case of smooth deformable particles. This non-monotonic effect is probably due to the fact that the particle stiffness also controls the deformation of the roughness. In case of intermediate roughness and capillary numbers, the apparent roughness is smaller than at low Ca.

To conclude, this test case shows that the model reproduces very well analytical results in the case of rigid particles. Extension to slightly deformable particles allow to predict the existence of reversible trajectories at low Ca and for rather smooth particles.

It also lead to a phase diagram of the origin of irreversibility.

Conclusion

A complete model of lubricated contact and an unconditionally stable integration scheme for the interaction forces have been introduced in section 2. The model can be readily implemented as an interaction model in conventional granular simulations.

The scheme handles the case ε = 0 (perfectly smooth spheres) gracefully, and regularization by elasticity removes the singularity of lubrication when u → 0 + . Furthermore the visco-elastic coupling is such that the stability of explicit integrator for particle dynamics is independent of viscous terms, whereas viscosity dominated regimes were challenging most explicit schemes until now.

This work contributed an implementation in the open source discrete element code

Yade-DEM.org [35] (used for running the test cases of last section). The lubricated contact model presented here is available as part of that code and our results can be reproduced independently. Two test cases have been investigated. The first one evidenced the robustness of the method. In the viscous regime, the decreasing exponential trajectory has been recovered. In the regime where inertia is not negligible, some behaviours like damped bouncing (close to pure-elastic limit) and vibration after collision (soft limit) have been reproduced by the model. Two dimensionless numbers are shown to control the response in that simple two-sphere system, and a phase diagram is proposed, showing the interplay between stiffness, roughness and viscosity. These dimensionless numbers might bear some relevance to shear flows of granular suspensions.

The second test case, two spheres suspended in a sheared liquid, has been extensively used to compare models and theories. We found a good agreement with theoretical solutions for stiff spheres, and we reproduced qualitatively some features observed previously with soft particles. A phase diagram has been determined concerning the deviation from reversible trajectories in shear flows, which originates either from particle roughness, either from particle deformation.

n k b ε k n ν n u e u u n

) 2 . 1 . 3 . 1 .

 2131 Visco-elastic coupling Contact and lubrication are combined in a rheological model summarized in figure Three different length appear in the model: the gap between the solid surfaces u, the deflection of surfaces u e , and finally the change of center-to-center distance with respect to the stress free configuration: u n = |r 1 -r 2 | -2a, with r 1,2 the positions vectors and 2a the sum of radii. The gap u is always positive, u n and u e can be positive or negative (u e > 0 when lubrication produces traction between separating bodies).

  to distord (or to shear) the solid surfaces. Unlike the normal behaviour, described by a single scalar equation, the quantities v s , v e , and v, defining the tangential displacements, are vector quantities. They need to be handled as such to reflect changes in the direction of shearing and rotations of the particle pair. v s is the geometrical shear displacement, obtained by integrating shear velocity over time. v e is the elastic part of v s produced by shearing the particles around a sticking contact area. v is the accumulated plastic slip, i.e. v = v s -v e . The lubrication forces are activated by plastic slip, assuming that viscous stresses in the interstitial fluid are negligible when the contact sticks.

  ment, and the force evaluated alternatively in the spring component then in the viscoplastic component. A Coulombian slider accounts for the frictional solid contact (components of force F c s and F c n ), and a Newtonian dashpot with u-dependent viscosity reflects lubrication by the fluid. There is no need to integrate the total plastic slip vector v over time in practice since its rate of change suffice to integrate the above equations. The problem is now solved for the tree possible cases: no contact (u > εa), sticking elastic contact (v = 0), and slipping contact (||F c s || = µ m ||F c n ||). If there is contact, the force is first evaluated by assuming a sticking regime. Then, if the magnitude of the trial force F c * s exceeds the frictional threshold, the Coulomb condition is taken into account. Similarly as the normal part, we derive an implicit scheme based on backward Euler approximation for the time derivative, i.e. ve = v + e -v - e ∆t + O(∆t).

  about here.] The three dimensionless numbers (ε, K, and St) have been investigated in an extensive set of simulations in order to obtain a complete phase diagram of the type of response (damping, bouncing, oscillating). The complete phase diagram is a volume and some slices for St from 10 to 10 4 are presented in figure5. As showed in figure4d, simulations close to the transitions exhibits both behaviours. Then, the transition between areas are smooth rather than sharp. For the highest Stokes simulations (a), the trajectories exhibit mainly either oscillating or bouncing according to the particles stiffness. The softest ones (K < 0.7) are oscillating after first collision whereas the stiffest ones bounces. Even at this Stokes number, there are still some parameters where the ideal exponential damping predicted by equation 36 is recovered. Reducing the Stokes number (from (a) to (d)) let more and more configurations being damped, but only for stiff cases. The oscillating area is being affected for lower Stokes number (e,f).

  the fluid is accounted by a combination of Stokes drag and lubrication force and torque. As stated in the model description, the lubrication analytical expressions are only valid for small gaps, and are set to zero in the model when the gap exceeds an arbitrary distance. We have varied this lubrication cut off distance from 2a to 7a. The results are displayed in Fig 9, in terms of the minimal gap reached during the simulation. Clearly, this cut off has a rather strong influence on this minimal gap -although the trajectories remain nearly unchanged overall. A larger cur-off distance tends to increases the minimal gap,

Figure 1 :

 1 Figure 1: Rheological scheme for normal component.

Figure 2 :

 2 Figure 2: Rheological scheme the tangential component.

Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 345678 Figure 3: Representation of falling sphere test case

Figure 10 : 2 *Figure 11 :

 10211 Figure 10: Comparison of trajectories and distances obtained with all components of lubrication, with normal lubrication only and compared to the analytical solution of Da Cunha and Hinch [31].

=

-2 ) For this kind of initial separation distance, the streamline is very close from the closed-loop streamline border for smooth hard spheres. Therefore, if any streamline deviation occurs, the far deflection gives the height of an area that is not reachable anymore after spheres interaction. The red area highlights the cases where this distance is strictly 0. It approximately corresponds to cases where ε < 3 × 10 -5 and Ca < 10 -4 . From this figure, we can deduce the main origin of irreversibility above these limits. When ε Ca, the streamline deviaton is controlled by the roughness,

Acknowledgements

The authors acknowledge Univ. Grenoble Alpes for financial support. 3SR and LRP are members of Labex Tec21, and LRP of Institut Carnot Polynat. Computations presented in this paper were performed using the CIMENT infrastructure (https:// ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT

The equation of motion for this system gives

mü n = -mgk b u e mü n = -mgk n max(0, εau) -3 2 πη f a 2 u u u n = u + u e [START_REF] Lac | Hydrodynamic interaction between two identical capsules in simple shear flow[END_REF] This system is non-linear, and as far as we know, it cannot be solved analytically. However, it is still possible to obtain partial solutions with some simplifications. For quasi-static motion controlled by viscous damping and hard spheres (u e → 0), this system can be reduced to:

which has the solution

where τ c is the characteristic sedimentation time:

At static equilibrium, there must be a persistent contact through asperities (for finite roughness at least), with the normal force balancing weight, and all velocities must be null. The elastic displacement at equilibrium is thus given by

Note that this equilibrium does not exist for heavy yet soft particles, when a < mg/k n .

In this case, the gap distance still decreases exponentially but it does so with an apparent weight reduced by the contact force corresponding to u = 0: F c n (u = 0) = -k n εan. The problem involves two dimensionless numbers, in addition to the relative roughness . Let us define the "contact number" K, which balances the contact force and the particle weight at equilibrium, i.e.

This number can be generalized in the case of suspensions, as K = Eε/P p , where E is the particle elastic modulus and P p the confining pressure applied on the particle phase. For the present test case, if the contact number goes below 1, there is no static equilibrium.


Quang Tuan

MD Rémy Nguyen

MD Alban Gravier

MD Constance Lesoil

MD Alexandre Bedet

MD Camille Petit-Hoang

MD Matthieu Mahevas

MD, PhD Armand Mekontso-Dessap

MD, PhD Jérôme Hodel

Anne-Catherine Bachoud-Lévi

MD, PhD Laurent Cleret De Langavant
email: laurent.cleret@gbhi.org

MD, PhD Anne-Catherine Bachoud Levi

Affiliations

Acute Hippocampal Encephalopathy in Heavy Cannabis Users: About 2 Cases

Keywords: Cannabis, Memory, Encephalopathy, Hippocampus, Addiction Running head Cannabis-related hippocampal encephalopathy

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Abstract:

Background:

Cannabis use is increasing worldwide despite the various health effects of this substance.

Methods:

We report two cases of acute hippocampal encephalopathy in heavy cannabis users (>10 joints/day).

Results:

In both male patients, acute encephalitis was suspected. Brain MRI diffusion-weighted sequences showed bilateral high signal abnormalities in hippocampal regions. Patients had renal dysfunction, rhabdomyolysis and inflammatory syndrome. Investigations showed no evidence of infectious or autoimmune encephalitides. Repeated electroencephalograms revealed no epileptic activity. Clinical, biological and MRI acute abnormalities improved within weeks. New exposition to cannabis yielded a new episode of encephalopathy. In both patients, severe long-lasting episodic memory impairment associated with hippocampal atrophy were observed several months later.

Conclusions:

Health professionals should be aware of this cannabis-related syndrome given its severe and long-lasting effects.

Introduction:

The number of cannabis users is increasing along with the progressive legalisation of both recreational and medical cannabis worldwide. Acute cardiovascular, respiratory, cognitive, psychiatric and public health consequences have been reported in cannabis users, as well as chronic cardiovascular and respiratory effects. [START_REF] Sachs | Safety and Toxicology of Cannabinoids[END_REF][START_REF] Volkow | Adverse health effects of marijuana use[END_REF] The cannabinoid hyperemesis syndrome [START_REF] Allen | Cannabinoid hyperemesis: cyclical hyperemesis in association with chronic cannabis abuse[END_REF] associates cyclical nausea, vomiting and abdominal pain relieved by hot showers.

Here, we report two cases of acute hippocampal encephalopathy in heavy cannabis users, with long-lasting memory impairment and delayed hippocampal atrophy.

Case 1

In October 2016, a 32-year-old man was admitted for delirium associated with inflammatory syndrome, rhabdomyolysis, and acute kidney injury.

He was a heavy user of cannabis (10 joints/day since his twenties), tobacco and alcohol (30-45 units/week). He had generalized epilepsy diagnosed nine years earlier and treated with Levetiracetam 1000 mg/day and Lamotrigine 200 mg/day. He had not travelled abroad in the preceding year.

The patient was found at home, disoriented and calm. The first clinical examination revealed delirium and agitation. The first blood tests found elevated leukocytes (30000/ mm 3 , 90% neutrophils), subnormal creatinine (117 μmol/L) and normal creatine kinase (CK = 250 IU/L). Within 12 hours, the patient's clinical condition worsened with severe delirium, agitation, and fever (38.5°C). A second blood test demonstrated rhabdomyolysis (CK = 1600 IU/L), renal dysfunction (serum creatinine 558 μmol/L at day 3) and inflammatory syndrome (33.900 leukocytes/mm 3 , C-Reactive Protein CRP = 61.5 mg/L, procalcitonin = 2.96 μg/L). Troponin T (20ng/L) and NT-proBNP (8579 IU) raised with no electrocardiogram changes. Serum lactate and carboxyhemoglobin were normal. The patient was admitted in intensive care unit for suspicion of acute encephalitis.

The cerebrospinal fluid analysis found no evidence of meningitis (<10 leukocytes /mm 3 , no germs, normal protein and glucose levels). Further investigations showed no evidence of antibodies related to autoimmune encephalitides. The toxicological test was positive for cannabis but negative for other drugs and alcohol. The brain MRI achieved 2 days later showed bilateral hippocampal abnormalities (Figure 1) with no gadolinium enhancement. The repeated electroencephalograms (EEG) showed no epileptic activity.

The patient was treated by intravenous acyclovir despite two negative Herpes Simplex Virus

Polymerase Chain Reaction (PCR) in the cerebro-spinal fluid. Delirium and agitation improved after one week, but episodic memory impairment remained severe (Table 1).

Biological abnormalities were normalized within five days. One month later, the patient returned home and discontinued his alcohol and cannabis intoxications. In May 2017, memory performance remained severely impaired despite slight improvement (Table 1).

In June 2017, the patient was admitted for delirium, rhabdomyolysis (CK = 12984 IU/L), acute kidney injury (creatinine = 1158 μmol/L), and hippocampal abnormalities on brain MRI. The evolution was similar to the first episode, without any treatment with acyclovir.

Investigations revealed heavy and recent consumption of cannabis.

Brain MRI revealed hippocampal atrophy in September 2017 (Figure 1). Memory performance remained impaired in July 2018, despite discontinuation of cannabis and alcohol intoxications. The patient resumed work under supervision.

Case 2

In July 2018, a 38-year-old man was admitted for delirium, fever, rhabdomyolysis and renal dysfunction.

He was a heavy user of cannabis (10 to 15 joints/day for 10 years) but not of tobacco or alcohol. He had a seizure associated with transient cognitive disorders, renal dysfunction and rhabdomyolysis in October 2016 and had been treated since with Sodium Valproate 1500 mg/day. He had not travelled abroad in the preceding year.

In July 2018, the patient was found at home sitting on the toilet and stained with stool and urine. The clinical examination confirmed a delirium without fever. Blood tests revealed severe rhabdomyolysis (CK = 64815 IU/L), renal dysfunction (creatinine 124 μmol/L) and inflammatory syndrome (neutrocytes 37.000/mm 3 , but CRP=3.1mg/L). Serum lactate and carboxyhemoglobin were normal. Within hours, a fever (38.2°C) appeared and the patient's consciousness deteriorated. He was intubated and hospitalized in intensive care unit for suspicion of acute encephalitis.

Cerebro-spinal fluid analysis found no evidence of meningitis (<10 leukocytes, no germs, normal glucose level but elevated protein 0.7 g/L). Blood and cerebro-spinal fluid tests showed neither evidence of antibodies suggestive of autoimmune encephalitides.

Toxicological testing was positive for cannabis and negative for synthetic cannabinoid or other drugs. The repeated EEG did not reveal any epileptic activity. Despite raised troponin T (297 ng/L), there was no sign of acute coronary syndrome on electrocardiogram or echocardiography. The brain MRI on day 2 after admission showed bilateral hippocampal abnormalities, with no gadolinium enhancement (Figure 1).

The patient was treated with intravenous acyclovir which was halted after two negative Herpes Simplex Virus PCR in the cerebro-spinal fluid. Under the suspicion of autoimmune limbic encephalitis, the patient received intravenous immunoglobulins (2g/kg). Consciousness improved allowing extubation (day 2) and biological disorders normalized within a week.

Episodic memory disorders remained severe (Table 1). Thirteen days after admission, brain MRI found persistent hippocampal abnormalities.

Eight months later, episodic memory remained severely impaired with a significant impact in daily life and incapacity to return to work (Table 1). At that time, the brain MRI revealed hippocampal atrophy (Figure 1).

Discussion:

We reported two cases of acute hippocampal encephalopathy in heavy cannabis users. Both patients had a history of epilepsy. The biological abnormalities included renal dysfunction, rhabdomyolysis, and inflammatory syndrome. The clinical, biological and brain MRI abnormalities improved within days to weeks with the exception of long-lasting episodic memory impairment and hippocampal atrophy. Both patients had two episodes following heavy consumption of natural cannabis.

We found no evidence of conditions known to cause hippocampal lesions: ischemic stroke, transient global amnesia, status epilepticus, herpetic encephalitis, autoimmune limbic encephalitis, anoxic, hypoglycemic or carbon monoxide encephalopathy were discarded. [START_REF] Förster | Diffusion-weighted imaging for the differential diagnosis of disorders affecting the hippocampus[END_REF] Rather, these abnormalities could correspond to the effects of two psychoactive components of cannabis, ∆-9-tetrahydrocannabinol (∆-9-THC) and cannabidiol (CBD), that have distinct affinity for cannabinoid type 1 (CB1R) and type 2 receptors (CB2R). The CB1R are found in neuron terminals of several brain regions including the hippocampus [START_REF] Iversen | Cannabis and the brain[END_REF] as well as in enteric nervous system, striated myocardial and skeletal muscles [START_REF] Mendizabal-Zubiaga | Cannabinoid CB1 Receptors Are Localized in Striated Muscle Mitochondria and Regulate Mitochondrial Respiration[END_REF] and kidneys. [START_REF] Larrinaga | Expression of cannabinoid receptors in human kidney[END_REF] The cannabinoid type 2 receptors (CB2R) are found in immune cells. [START_REF] Pertwee | Cannabinoid pharmacology: the first 66 years[END_REF] Although CBD is described as anticonvulsant, 9 ∆-9-THC might be proconvulsant. [START_REF] Malyshevskaya | Natural (∆9-THC) and synthetic (JWH-018) cannabinoids induce seizures by acting through the cannabinoid CB1 receptor[END_REF] In both cases here, cannabis consumption preceded epilepsy occurrence suggesting a causal role of cannabis. Acute kidney injuries were reported in cannabis users and explained by acute tubular necrosis secondary to dehydration in hyperemesis syndrome or to rhabdomyolysis, or hypothetically caused by direct toxicity of synthetic cannabinoid. [START_REF] Park | Cannabinoids and the kidney: effects in health and disease[END_REF] Noteworthily, our patients used natural cannabis and did not exhibit hyperemesis syndrome. The spontaneous improvement of their renal function suggested acute tubular necrosis. Rhabdomyolysis was reported after synthetic cannabis consumption. [START_REF] Tournebize | Acute effects of synthetic cannabinoids: Update 2015[END_REF] Cannabis might interfere with CB1R in striated muscle mitochondria impacting oxidative activity and glucose metabolism. [START_REF] Mendizabal-Zubiaga | Cannabinoid CB1 Receptors Are Localized in Striated Muscle Mitochondria and Regulate Mitochondrial Respiration[END_REF] A long-lasting episodic memory disorder followed the acute hippocampal encephalopathy episodes in both patients. The tropism of cannabis neurotoxicity for hippocampus may be explained by hippocampal high concentration of CB1R. Long term memory deficits are reported in chronic cannabis users. [START_REF] Rocchetti | Is cannabis neurotoxic for the healthy brain? A meta-analytical review of structural brain alterations in non-psychotic users: Cannabis effects in non-psychotic users[END_REF] Reduced hippocampal grey matter volume is associated to cannabis use. [START_REF] Yücel | Hippocampal harms, protection and recovery following regular cannabis use[END_REF] These long-lasting effects may be caused by ∆-9-THC 15 while CBD may be neuroprotective. [START_REF] Lorenzetti | The Role of Cannabinoids in Neuroanatomic Alterations in Cannabis Users[END_REF] Acute hippocampal encephalopathy can occur in heavy cannabis users, but both reported cases could have genetic predisposition to this syndrome. [START_REF] Schacht | Associations between cannabinoid receptor-1 (CNR1) variation and hippocampus and amygdala volumes in heavy cannabis users[END_REF] Furthermore, cannabis withdrawal after admission might have worsened patients' condition. Depriving the CB1R might cause a cascade of cerebral, renal, muscular, and cardiac disorders, while the CB2R liberation might cause an inflammatory syndrome. Future work is needed to ascertain the role of cannabis components and to assess the therapeutic value of using CBD in this syndrome.

Figure legend Figure 1 .

 legend1 Figure legend

Table legend Table 1 . Cognitive performance at acute and chronic phases.

 legend1 Episodic memory is severely impaired in both patients, at both acute and chronic phases. The patient case 2 additionally shows naming deficits partly explained by lower proficiency in French language. *** means that patients perform below 3 standard deviations from the mean score of a group of controls of the same age and educational level. FCSRT: Free and Cued Selective Reminding Test (Buschke, 1984). ROCF: Rey-Osterrieth Complex Figure (Osterrieth, 1944). Trail Making Test (Reitan, 1955). Stroop test (Stroop, 1935). Language naming (Deloche, 1996). Pyramids and Palm Trees test (Howards and Patterson, 1992). Mini Mental Status (Folstein, 1975).

	Table 1
		Case 1	Case 2
	Cognitive Assessment Date	Nov 2016	May 2017	Aug 2018	Feb 2019
	Episodic Memory
	FCSR Identification /16	16	16	16	16
	FCSRT Immediate Cued Recall /16	12***	11***	8***	9***
	FCSRT Total Free Recall /48	3***	14***	7***	11***
	FCSRT Total Recall 1 /16	4***	8***	4***	6***
	FCSRT Total Recall 2 /16	0***	13***	6***	7***
	FCSRT Total Recall 3 /16	0***	11***	8***	12***
	FCSRT Delayed Free Recall /16	0***	6***	0***	0***
	FCSRT Delayed Total Recall /16	0***	14***	6***	5***
	ROCF Recall	0***	0***	3***	9***
	Attention and Executive Function
	Forward Digit Span	6	6	6	6
	Backward Digit Span	4	4	4	5
	Trail Making Test A (sec)	24	27	37	42
	Trail Making Test B (sec)	47	53	87	102
	Stroop Colors	57	52	55	56
	Stroop Words	47	44	42	46
	Stroop Interference	90	93	107	99
	Categorical fluency (animals 120 sec)	27	26	23	17
	Lexical fluency (letter P 120 sec)	21	18	8	10
	Instrumental Function
	Language naming /80	79	80	60***	69***
	Pyramid and Palm Tree Test /52	49	50	-	-
	Calculation /28	28	28	28	28
	Gestures (ideomotor, motor) /40	40	40	40	40
	ROCF Copy	36 (type I)	36 (type I)	36 (type IV)	36 (type I)
	Global efficiency
	Mini Mental Status Examination /30	24	26	22

Acknowledgements

The authors thank both patients who agreed for reporting their history.

Funding source

LCL and ACBL are supported by ANR-17-EURE-0017.

Conflicts of Interest

The authors report no conflict of interest.