S Athani

B Metzger

Y Forterre

R Mari

Transient flows and migration in granular suspensions: Key role of Reynolds-like dilatancy

, which prescribes the dilation rate of the suspension over a strain scale γ 0 , quantitatively captures the suspension dilation/compaction over the whole range of parameters investigated. Together with the Darcy flow induced by the pore pressure gradient during dilation or compaction, this Reynolds-like dilatancy implies that the early stress response of the suspension is nonlocal, with a nonlocal length scale ℓ which scales with the particle size and diverges algebraically at jamming. In regions affected by ℓ, the stress level is fixed, not by the steady-state rheology, but by the Darcy fluid pressure gradient resulting from the dilation/compaction rate. Our results extend the validity of the Reynolds-like dilatancy flow rule, initially proposed for jammed suspensions, to flowing suspension below ϕ c , thereby providing a unified framework to describe dilation and shear-induced migration. They pave the way for understanding more complex unsteady flows of dense suspensions, such as impacts, transient avalanches or the impulsive response of shearthickening suspensions.

Introduction

Granular suspensions consisting of rigid particles and fluids are ubiquitous in nature (landslides including debris flow, mud flow and submarine avalanches) and industrial applications (among others manufacturing and handling of concrete, drilling slurries or molten chocolate) [START_REF] Guazzelli | Rheology of dense granular suspensions[END_REF]. In many flow situations, the particle volume fraction of the suspension is not constant but can evolve both in space and time. This is observed for instance in volume-imposed configurations with neutrally-buoyant particles, when shear-induced migration occurs as a result of an inhomogeneous flow field, such as in pipe or large gap Couette flows [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF]. It is also observed in pressure controlled configurations, for instance in gravity driven flows with non-buoyant particles, when an immersed granular avalanche dilates in order to flow [START_REF] Iverson | Acute Sensitivity of Landslide Rates to Initial Soil Porosity[END_REF]. For both of these situations, the flow involves local changes of the volume fraction.

Upon dilation, pore space is created in between the particles, which has to be filled by the incompressible fluid. The resulting relative motion between the fluid and solid phases creates a drag on the particles which transiently increases the particle pressure. A similar scenario, but of opposite sign, occurs upon compaction: this time fluid is expelled out of the contracting pores, which transiently decreases the particle pressure.

Depending on the scientific community and the flow regime considered, this two-phase flow coupling has been comprehended in two very different ways. In the suspension community, transient changes of the particle volume fraction below ϕ c (the critical volume fraction at which the suspension jams), are usually described using the so-called Suspension Balance Model (SBM) [START_REF] Nott | Pressure-driven flow of suspensions: Simulation and theory[END_REF][START_REF] Morris | Pressure-driven flow of a suspension: Buoyancy effects[END_REF][START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF][START_REF] Nott | The suspension balance model revisited[END_REF]. In this framework, particle migration is driven by particle stress gradients arising in inhomogeneous flow field [START_REF] Deboeuf | Particle Pressure in a Sheared Suspension: A Bridge from Osmosis to Granular Dilatancy[END_REF]. SBM is built on phase averaged momentum and mass conservation laws [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid[END_REF] and two closures: an interphase coupling and a constitutive law for the particle phase stress. In the Stokes regime, the interphase coupling is usually written as a drag proportional to the relative phase velocities. The standard choice for the constitutive law is the steady-state rheological flow rules for the shear and normal particle stresses, which have now been well characterized and tested (Denn & Morris 2014;[START_REF] Guazzelli | Rheology of dense granular suspensions[END_REF][START_REF] Ness | The Physics of Dense Suspensions[END_REF]. These constitutive laws can be equivalently expressed either as "volume-imposed" or "pressure-imposed" flow rules. In volume-imposed configurations, dimensional analysis requires that, in steady state, stresses are viscous and proportional to η f γ, where η f is the fluid viscosity and γ the shear rate. The constitutive laws thus resume to the volume fraction dependence of the shear η s (ϕ) and normal η n (ϕ) viscosities, where ϕ is the particle volume fraction. In pressure-imposed configurations, the dimensionless number controlling the flow is the viscous number J ≡ η f γ/P p . Here P p is the socalled "particle pressure", or more exactly the opposite of the particle normal stress along the y direction. It is the part of the (opposite of the) normal stress which is due to contact and lubrication forces between adjacent grains, although in practice for the large volume fractions considered here the particle pressure is dominated by contact forces [START_REF] Gallier | Rheology of sheared suspensions of rough frictional particles[END_REF]. The constitutive laws are then provided by the relations µ(J) and ϕ SS (J), where µ = τ /P p is the suspension friction coefficient, with τ the applied shear stress. So far, SBM has been tested and found to be in reasonable agreement with experimental results on particle migration in wide-gap Couette [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF][START_REF] Sarabian | Fully developed and transient concentration profiles of particulate suspensions sheared in a cylindrical couette cell[END_REF], in pipe flows of concentrated suspensions (Snook et al. 2016), or during resuspension [START_REF] Acrivos | Shear-induced resuspension in a couette device[END_REF]d'Ambrosio et al. 2021;[START_REF] Saint-Michel | X-ray radiography of viscous resuspension[END_REF].

In the granular and soil mechanics communities, volumetric changes of the particle phase are usually studied for packings initially at rest, and prepared very close to, or even beyond ϕ c . A central concept in this field is the Reynolds dilatancy [START_REF] Reynolds | on the dilatancy of media composed of rigid particles in contact. with experimental illustrations[END_REF][START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF]), which stipulates that, in order to deform, an initial dense (resp. loose) packing must dilate (resp. compact) toward its critical state of volume fraction ϕ c in the quasistatic regime. In the presence of an interstitial fluid, this change in volume fraction leads to Darcy back-flow (or negative/positive pore pressure) whose drag on the particles sets the particle stress. This Darcy-Reynolds coupling between granular dilatancy and Darcy back-flow is essential for predicting the behavior of soils under drained or undrained conditions [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF]. It was also shown to have major consequences for the transient flow of granular suspensions close to jamming, such as for the onset of debris flows and submarine avalanches [START_REF] Iverson | Acute Sensitivity of Landslide Rates to Initial Soil Porosity[END_REF][START_REF] Pailha | Initiation of underwater granular avalanches: Influence of the initial volume fraction[END_REF][START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF][START_REF] Topin | Collapse dynamics and runout of dense granular materials in a fluid[END_REF][START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF][START_REF] Montellà | A two-fluid model for immersed granular avalanches with dilatancy effects[END_REF], during impacts [START_REF] Jerome | Unifying impacts in granular matter from quicksand to cornstarch[END_REF][START_REF] Jørgensen | Deformation upon impact of a concentrated suspension drop[END_REF], silo discharge [START_REF] Kulkarni | Particle-pressureinduced self-filtration in concentrated suspensions[END_REF], or when shearing dense clouds of particles [START_REF] Metzger | Clouds of particles in a periodic shear flow[END_REF]. In the context of initiations of avalanches from initially loose sedimented packings, [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] modeled this feedback by introducing a dilatancy dynamics for ϕ, which is essentially a relaxation of ϕ towards its steady-state law ϕ SS (J) on a typical strain scale γ 0 ∼ O(1), extending previous Reynolds dilatancy laws proposed for dry granular media in the quasi-static regime [START_REF] Roux | Texture-Dependent Rigid-Plastic Behavior[END_REF].

Interestingly, the theory based on Reynolds dilatancy and the Suspension Balance Model describe the same type of interplay between dilation (compaction) of the granular skeleton and fluid interphase drag forces. However, this coupling is treated in very different ways. In the Suspension Balance Model, the particle stress field is set by the steady-state rheological flow rules and the migration rate adapts to satisfy the force balance between this stress field and the interphase drag force. Conversely, in the theory based on Reynolds dilatancy, the rate of dilation is geometrically imposed by a "dilatancy angle" (the distance between the actual volume fraction to ϕ SS (J)). In this case, it is the particle stress field that adapts to this transient kinematically-constrained evolution. This fundamental difference raises important questions. A particularly significant issue is the origin of the stress levels observed during a sudden change of flow condition, as depicted in the examples given in figure 1(a). Are the typical pore and particle pressures set by the steady-state rheology, or by the Darcy back-flow induced by a geometrically constrained dilation (compaction) of the granular phase? Another question is whether the concepts of Reynolds dilatancy and shear-induced migration could be reconciled and described in a unified framework? Clearly, the Suspension Balance Model cannot be used for very compact granular layers, as above ϕ c , the steady-state rheological flow rules are not defined. However, whether the Reynolds dilatancy concept applies to describe migration for systems below ϕ c has never been tested.

In this article, we address the above questions by performing Discrete Element Method (DEM) simulations in the canonical setup shown in Figure 1(b): a neutrally buoyant sheared suspension in a pressure-imposed rheology configuration, subject to a sudden change in shear rate or external pressure (i.e. a step change in J ext = η f γ/P ext ). This simple configuration allows us a detailed characterization of the dilation (compaction) dynamics, while varying all the relevant control parameters (initial and final J ext , ratio of system size to particle radius H/a). We then compare these numerical results to the two different continuum models described above: the Suspension Balance Model based on the steady-state rheological flow rules (from now on simply called steady-state rheology model ), and the two phase model coupling the Darcy drag force and the Reynolds-like dilatancy closure of [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] (from now on called Darcy-Reynolds model, following [START_REF] Jerome | Unifying impacts in granular matter from quicksand to cornstarch[END_REF]), which has so far only been validated when the initial state is a loose sedimented packing, as opposed to a flowing steady state.

We find that, while at large strains both continuum models give similar results in agreement with DEM, at small strains, only the Darcy-Reynolds model can quantitatively capture the DEM results over the whole range of control parameters. The steady-state approximation is not only quantitatively poorer, it also provides qualitatively wrong predictions. Specifically, we show that the Darcy-Reynolds model induces non-locality in the stress upon dilation (compaction), on a typical length scale ℓ proportional to the grain size and diverging at the jamming transition, which the steady-state rheology model is completely oblivious to. This nonlocality implies that the transient excess (respectively deficit) of macroscopic stresses observed during dilation (respectively compaction) is sizedependent. Systems whose size H is smaller than ℓ show smaller (respectively larger) level of stresses than expected from the steady-state rheology. In this case, stresses are set by the Darcy back-flow induced by the volumetric strain. Conversely, systems much larger than ℓ show larger (respectively smaller) stresses, well predicted by the steady constitutive law. Overall, our results highlights that during transients, stress levels in granular suspensions, even below ϕ c , are controlled by a Reynolds-like dilation (compaction) of the granular phase and not by the steady-state rheological flow rules.

The manuscript is organized as follows: After presenting the numerical set-up in § 2, the DEM results are reported in § 3. We then detail both continuum models in § 4 and compare them to DEM results in § 5. Conclusions are drawn in § 6. The system is composed of mobile particles (in grey) which are immersed in a suspending Newtonian fluid and placed between the top and bottom walls constructed with frozen particles (in blue). The control parameters are the external stress Pext on the top-wall and the shear rate γ imposed to the entire system. In this pressure-imposed rheology configuration, the top wall can move vertically, allowing the volume fraction to adjust freely owing to the value of the imposed Jext = η f γ/Pext. The fluid is not explicitly modeled but accounted for through lubrication and drag forces. b) The drag force fD exerted by the fluid on particle i is proportional to the difference of between the local fluid and particle velocities. This force ensures that grains are advected by a homogeneous background shear flow in the x-direction, and models the pore pressure feedback (or interphase drag coupling) occurring during dilation (compaction) of the granular skeleton in the y-direction.

Numerical Setup

Pressure-imposed rheology DEM simulations are developed based on an established DEM code for volume-imposed simulations presented in detail in [START_REF] Mari | Shear thickening, frictionless and frictional rheologies in non-brownian suspensions[END_REF]. As depicted in Fig. 2 (a), we consider a monolayer of neutrally buoyant non-Brownian hard disks (in gray) immersed in a Newtonian fluid of viscosity η f , placed between two rigid top and bottom walls built out of frozen particles arranged in a disordered configuration (in blue). The averaged vertical position of the particles in the bottom wall defines the vertical coordinate y = 0, while the averaged vertical position of the particles in the top wall is at y = H. The bottom wall is not permeable to the fluid nor to the particles, while the top wall is permeable to the fluid, but not to the particles. Furthermore, the bottom wall is fixed, whereas the top wall is free to move vertically, subject to an externally applied normal stress P ext , and follows a prescribed horizontal velocity γH when the suspension is sheared at an imposed shear rate γ. This corresponds to an applied viscous number J ext ≡ η f γ/P ext . This setup is, in spirit, the two-dimensional equivalent of the pressure-imposed rheometer of [START_REF] Boyer | Unifying suspension and granular rheology[END_REF]. We choose a monolayer setup over a three-dimensional one to be able to reach reasonably large linear extensions, without simulating a prohibitively large number of particles. The particle size distribution is bidisperse, with a size ratio of 1.4, mixed in equal volume. We use periodic boundary conditions along the horizontal x-direction, and call l x the system length in this direction. The solid fraction of the system is then ϕ ≡ N πa 2 /(l x H), with a the mean particle radius. We assume that the system is in the limit of vanishing Reynolds number, and neglect inertial effects. In this limit, we thus follow the Stokesian Dynamics approach and do not explicitly simulate the fluid phase, which only comes in as hydrodynamic interactions acting on the particles [START_REF] Brady | Stokesian Dynamics[END_REF]. The dynamics being overdamped, the equation of motion consists of mechanical force and torque balance on every particle.

The forces acting on bulk particles are contact forces, with a Coulomb frictional model with friction coefficient µ p = 0.5, implemented in a standard Cundall-Strack manner [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], and hydrodynamic forces. The particle stiffness is chosen such that the particle overlap is always smaller than 2 % of their diameter. We neglect long-range hydrodynamic interactions, which are screened in a dense system, and consider only the short-ranged pairwise lubrication forces [START_REF] Mari | Shear thickening, frictionless and frictional rheologies in non-brownian suspensions[END_REF]. Finally, a Stokes drag is applied to particles, and is assumed proportional to the relative velocity between a particle and the background fluid velocity evaluated at the center of the particle, u ∞ (y) = ( γy, 0) if the particle is at height y: for particle i with radius a i and velocity u i , the drag force is f D = -6πη f a i (u i -u ∞ (y i )), see Figure 2 (b). This force plays two important roles: first it ensures the granular suspension is advected by an homogeneous background shear flow in the x-direction. Second, it is used to model, in the simplest possible way, the pore pressure feedback (or interphase drag coupling) occurring during dilation (compaction) of the granular phase in the y-direction. Note that two important simplifications are made. First, we neglect the vertical fluid velocity component that should arise to satisfy the suspension incompressibility during dilation (compaction). It is important to note that, for consistency, the same approximation will be made in the two-phase flow modeling in Section 4, where we show that this approximation does not change the physics of two-phase coupling and only results in an underestimation of the hydrodynamic resistance by a factor of less than 2. Second, we do not explicitly simulate the fluid flow between pores, nor do we use standard mesoscopic models as the empirical hindered settling function proposed by [START_REF] Richardson | Sedimentation and fluidisation: Part i[END_REF], which scales as (1-ϕ) α , with α > 0 [START_REF] Davis | Sedimentation of Noncolloidal Particles at Low Reynolds Numbers[END_REF][START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF]Snook et al. 2016).

This choice of interaction between the solid and fluid phases, through a simple Stokes drag force, retains the minimal physics needed to describe the transient dynamics considered here: a drag force proportional to the difference of vertical velocities between fluid and particle phases. Such a choice will prove trivial to model when it comes to the comparison of our DEM results with continuum modeling. This is an important aspect of our approach, as it will allow a much finer testing of the constitutive laws. Indeed, instead of having to attribute discrepancies between the continuum model and the numerical simulations to one of the two closures (the interphase drag and the constitutive law), which can rapidly prove impossible to disentangle, we know here by construction that the observed discrepancies are coming from inaccurate constitutive modeling only. The downside to this modeling choice is that we will not achieve quantitative agreement with experimental data, but this was already prevented by our choice of a monolayer.

Forces acting on the top wall are the imposed external force P ext l x , the contact and lubrication forces coming from interaction with the bulk particles, and a viscous drag proportional to the relative motion of the wall with respect to u ∞ (H), which reads -κl x (u wall -u ∞ (H)), with u wall = ( γH, ∂ t H) the velocity of the top wall and κ the wall hydraulic resistance. Except when noted, we take κ as if a Stokes drag was acting on each particle composing the wall, that is κ = 6πη f N wall a/l x , with N wall the number of particles in the top wall. We also performed simulations with κ = 0 which yields similar results. The vertical motion of the top wall is obtained by writing that the wall is force-free and its equation of motion is derived in Appendix A.

The simulation protocol is as follows: the suspension is first pre-sheared under constant pressure P ext and constant shear rate γ so that it reaches an initial steady state characterized by the initial viscous number J i . We then suddenly impose a step change in J ext . This can equivalently be done by keeping P ext constant while changing the imposed shear rate γ (as shown in Figure 1), or by keeping γ constant while changing the imposed external pressure P ext . Before the system eventually reaches a final steady state characterized by the viscous number J f , we systematically investigate the transient evolution of the system by monitoring the evolution of the top wall vertical position H, the depth-averaged shear stress τ , and the volume fraction ϕ and particle vertical normal stress P p profiles within the granular layer along the y-direction. Simulations are performed for both dilation (J f > J i ) and compaction (J f < J i ), over a wide range of viscous number J ext ∈ [10 -4 -10 -1 ], and for three systems sizes differing by their particle number (N = 1500, 1000 and 600) or corresponding vertical normalized height (H/a = 140, 68 and 33, respectively, when J ext = 10 -4 ).

DEM simulation results

Macroscopic phenomenology

Figure 3 shows the evolution of the granular layer subject to a step change in J ext for two cases of dilation (J f > J i , blue line) and compaction (J f < J i , red line). After the step, we observe that when J f > J i (resp. J f < J i ), the thickness of the granular layer H increases (resp. decreases) as a result of the dilation (resp. compaction) of the granular skeleton. This behavior is expected since the steady-state rheological rule ϕ SS (J) prescribes that the volume fraction of the suspension must decrease with its viscous number. Interestingly, the strain scale to reach the new steady state at J f strongly differs for dilation and compaction. As seen in Fig. 3a, the volumetric change occurs on a strain γ of order 1 for compaction, while it takes more than 10 strain units to complete for dilation. This simply comes from the fact that dilation and compaction are asymmetric as they occur under different external values of J ext .

Remarkably, during the transient, both the mean particle pressure Pp /P ext (averaged over the whole layer thickness) and the depth-averaged shear stress τ /τ SS exhibit a sudden and large increase (resp. decrease) for dilation (resp. compaction), before relaxing towards their steady-state value. Understanding what sets the stress levels during these transients is the main goal of the present study. We can readily compare the magnitude of these stress jumps, occurring right after the step, to the prediction of the steady-state rheology model. At γ = 0 + , the particle packing fraction profile must be the same as that of the initial steady-state preparation. It is therefore homogeneous along the y-direction such that ϕ| γ=0 + = ϕ| γ=0 -= ϕ SS (J i ). Assuming the steady-state flow rule ϕ SS (J) applies at all times, the viscous number inside the granular layer just after the jump is thus also the same as that just before the jump, i.e. J| γ=0 + = J| γ=0 -= J i , even though the external imposed viscous number has changed to J ext = J f . Since the shear rate is imposed to be homogeneous, the particle pressure profile must also remain homogenous and equal

to P p | γ=0 + = (η f γ/J)| γ=0 + = η f γ| γ=0 + /J i yielding P p P ext γ=0 + = J f J i . (3.1) Similarly, one obtains (τ /τ SS )| γ=0 + = µ(J i )J f /µ(J f )J i .
The steady-state rheology model thus predicts that, right after the step, the magnitude of the particle stress jump is given by the ratio of applied J ext , after and before the step. This early strain prediction of the steady-state rheology is highlighted by the orange stars in Figure 3. Strikingly, we find that the steady-state rheology overestimates by more than one order of magnitude the particle pressure under dilation obtained from the DEM simulations, and underestimates it by a similar amount under compaction. The same discrepancy is observed for the shear stress overshoots, τ /τ SS , which are largely overestimated under dilation, and underestimated under compaction by the steady-state rheology. These initial observations show that steady state rheology fails to capture the instantaneous stress response of a suspension subjected to a sudden disturbance. This discrepancy, although occurring on a short strain scale of order 0.1, is a serious problem when describing the configurations shown in the figure 1(a). For instance, during an impact, knowing the stress level after the impact is crucial since it sets the deceleration time scale of the drop.

Transient evolution inside the bulk

To further investigate the transient response of the granular layer, we now turn to local measurements of the particle pressure P p (y) and solid fraction ϕ(y) profiles inside the bulk of the granular layer for a step increase in J ext (dilation case), see Figure 4. After the step, the particle pressure P p (y) and the volume fraction ϕ(y) progressively evolve from their initial (red dashed line) to final (black dashed lines) steady-state profiles. The striking observation is that, even at a strain as small as γ = 0.01 after the step, the particle profile obtained from the DEM simulations (dark blue solid line) significantly departs from the prediction of the steady-state rheology model evaluated at γ = 0 + : P p (γ = 0 + , y)/P ext = constant = J f /J i (grey solid line). The particle pressure exhibits smaller stress levels than expected and a large gradient over the cell height. The volume fraction profile also shows a gradient, but it becomes significant only after a strain of order 1. At γ = 0.01 after the step, the pressure is roughly three times smaller at the top than at the bottom of the system, but the volume fraction is nearly uniform. These observations further evidence the inability of the steady-state rheology to describe the early stress state of the granular layer.

To characterize more systematically these differences between the prediction of the steady-state rheology and DEM results, we perform simulations for many combinations of the three control parameters of the system: J i , J f and H/a. In Figure 5, we report the normalized particle pressure at the bottom of the cell P p (y = 0)J i /P ext J f , just after the step change in J ext , versus J f /J i . In this representation, the steady-state rheology prediction at γ = 0 + (Eq. 3.1) corresponds to the horizontal line P p (y = 0)J i /P ext J f = 1 (thick grey line). We find that many simulation data are in good agreement with this prediction, but equally many deviate from it, sometimes by more than an order of magnitude. We can also observe that the largest discrepancy occurs for small systems (purple and green squares in Figure 5), while in large systems the pressure tends to be well predicted by the steady-state rheology model (brown and blue diamonds). This points to the existence of a characteristic length scale in the suspension, which, in the case shown in Figure 5, would be of order of a few tens of particle diameters. For systems larger than this length scale, the steady rheology accounts well for the pressure level at the bottom, while for smaller systems, finite-size effects affect the pressure levels. Interestingly, when we vary J f at a fixed J i value, for small systems the compaction and dilation cases show asymptotically different behaviors. Under abrupt dilation (J f /J i ≫ 1), P p (y = 0) saturates to a J f -independent value different from P i , while under abrupt compaction (J f /J i ≪ 1), P p (y = 0) roughly scales as P ext .

The above results suggest that the steady-state rheology is in many cases unable to predict the early transient evolution of the stress levels in the granular suspension. Note that in practice, the earliest strain at which DEM results are reported in Fig. 4 and 5 is γ = 0.01, since before that, DEM results may be dependent on the stiffness of the particles. One could thus argue that the discrepancy between DEM results and the steady-state rheology predictions arises from the fact they are not evaluated exactly at the same strains (γ = 0.01 versus γ = 0 + , respectively). However, in what follows, by deriving the full transient evolution of the granular layer using both the steady-state and Darcy-Reynolds models, we show that this discrepancy arises from a more fundamental mechanism.

Continuum models

In this section, we recall the usual two-phase continuum description for suspensions [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid[END_REF][START_REF] Jackson | The dynamics of fluidized particles[END_REF]. Importantly, the derivation is done by following the same simplifications that we adopted in the DEM simulations. First, we use a simplified expression for the interaction between the particle and fluid phase based on a Stokes drag. This drag force is used to simply account for the pore pressure feedback induced by the dilation (compaction) of the particle phase. Second, the vertical fluid velocity resulting from the mass conservation (fluid counterflow) is neglected, which, as described in Appendix B, only marginally affects the results. We then present two competing constitutive models to close this continuum description: one assuming steady rheology at all time (steady rheology model), the other taking into account a Reynolds-like dilatancy equation for the transient dynamics of the particle packing fraction (Darcy-Reynolds model).

Conservation laws

The continuum model aims at describing the time evolution of the particle phase stress profile P p (y) in the y-direction and solid fraction profile ϕ(y), and ultimately of the height H of the top wall. It relies on the mass and momentum conservations and a constitutive relation relating the volume fraction ϕ with the particle stress P p and the shear rate γ. As our dilation problem is translation invariant in the x-horizontal direction, it is essentially a one-dimensional problem in the y-direction. We call u p (y) the vertical velocity of the particle phase and we neglect the vertical fluid velocity (see Appendix B for the full equations taking into account the fluid velocity). Mass conservation for the particle phase reads

∂ t ϕ + ∂ y (u p ϕ) = 0. (4.1)
Momentum conservation, in the approximation of the Suspension Balance Model (SBM) [START_REF] Nott | Pressure-driven flow of suspensions: Simulation and theory[END_REF][START_REF] Morris | Pressure-driven flow of a suspension: Buoyancy effects[END_REF][START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF][START_REF] Nott | The suspension balance model revisited[END_REF], balances the pressure gradient ∂ y P p (y) with the local interphase drag, which is directly proportional to u p as we recall that the vertical fluid velocity is neglected:

∂ y P p (y) = - η f a 2 ϕR(ϕ)u p , (4.2)
with R(ϕ) the hydrodynamic resistance of the particle matrix, which is dimensionless in 3D. This expression is adapted to our 2D DEM simulations with a simple Stokes drag using R(ϕ) = 6a, where the particle pressure P p is now a force per unit length. Importantly, the full two-phase equations including the fluid vertical velocity also leads to a particle pressure gradient that only depends on u p , but which differs from Eq. 4.2 by a factor (1 -ϕ) -1 , see Eq. B 5 in Appendix B. Therefore, neglecting the fluid velocity only amounts to underestimating the interphase resistance by a factor (1 -ϕ) -1 ≈ 2 (with the same approximation in the DEM), without impacting the qualitative features of the model. The boundary conditions for the problem are as follows. We require vanishing vertical velocity at the bottom wall u p (0) = 0, which from (4.2) yields ∂ y P p (y = 0) = 0. Moreover, force balance on the top wall prescribes

P ext + κ∂ t H -P p (H) = 0, (4.3)
with κ the resistance of the top wall, and we impose that the velocity of the particle phase follows the top wall velocity in y = H, u p (H) = ∂ t H. Similarly, to match our 2D DEM setup, we use κ = 6πη f N wall a/l x , with N wall the number of particles in the top wall.

We can non-dimensionalize this model with a the unit length, 1/ γ the unit time, and η f γ (η f γa in 2D) the unit stress. We denote by X the adimensionalized X, for any physical quantity X. Inserting momentum balance in the mass balance equation to eliminate u p , we obtain

∂ tϕ -∂ ŷ R(ϕ) -1 ∂ ŷ Pp = 0, (4.4)
with boundary conditions

0 = ∂ ŷ Pp 0 (4.5) 0 = J -1 ext - κ ϕ( Ĥ) R(ϕ( Ĥ)) ∂ ŷ Pp H -Pp ( Ĥ), (4.6)
with κ = κa/η f and R = R(ϕ) in 3D, and κ = κ/η f and R = R(ϕ)/a in 2D. Note that in practice the granular layer is always much thicker than the upper wall thickness, such that the hydraulic resistance of the upper wall is negligible relative to the bulk resistance of the suspension.

As such, the model is not closed, as we need a further relation between the normal stress and the solid fraction, that is, a constitutive model. In the following, we will consider two of them, both reducing to the usual [µ(J), ϕ SS (J)] rheology in steady state, but differing in their transient behaviors.

Steady-state rheology model

In the first model, we assume, as it is usually done when solving a shear-induced migration problem [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF]Snook et al. 2016;[START_REF] Sarabian | Fully developed and transient concentration profiles of particulate suspensions sheared in a cylindrical couette cell[END_REF], that the steady-state rheology is valid at all times and can describe the state of the suspension even during the transient migration of the particle phase. This simply sets ϕ = ϕ SS (J), (4.7)

with J = 1/ Pp and ϕ SS (J) behaving close to jamming as

ϕ SS (J) ≈ ϕ c -KJ β , (4.8)
which is equivalent to η n (ϕ) = K 1/β (ϕ c -ϕ) -1/β in a volume imposed formulation [START_REF] Degiuli | Unified theory of inertial granular flows and non-brownian suspensions[END_REF]. Inserting Eq. 4.7 in Eq. 4.4, we obtain that the volume fraction obeys the following diffusion equation

∂ tϕ -∂ ŷ [D(ϕ)∂ ŷ ϕ] = 0, (4.9) with a ϕ-dependent diffusion coefficient D(ϕ) = K 1/β /( R(ϕ)βJ 1+β ).
Just like for the temperature governed by the heat equation, after a step change of external pressure on the upper wall or after a global change of shear rate, the particle stress remains strictly homogeneous at γ = 0 + and exhibit a discontinuity at the upper wall boundary, see Fig. 6 (a). The dilation/compaction process starts from the upper wall and gradually diffuses within the bulk of the granular layer, simply as heat would diffuse within the layer after a step change in temperature of the upper wall, with a ϕ-dependent diffusion coefficient.

Darcy-Reynolds model

In the second model, we use the constitutive relation proposed by [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF], which assumes that during transients the volume fraction ϕ relaxes towards its steady-state value ϕ SS (J) in a finite strain scale γ 0 as

∂ tϕ + ûp ∂ ŷ ϕ = - 1 γ 0 [ϕ -ϕ SS (J)] . (4.10)
This simple closure was built as an extension at finite J of previous Reynolds dilatancy laws proposed for dry granular material [START_REF] Roux | Texture-Dependent Rigid-Plastic Behavior[END_REF][START_REF] Roux | Statistical Approach to the Mechanical Behavior of Granular Media[END_REF], and to recover the ϕ SS (J) rheology in steady state. It prescribes that the rate of dilation of the suspension is kinematically imposed by a "dilatancy angle" proportional to the difference between the actual volume fraction ϕ and ϕ SS (J). The relaxation strain scale γ 0 could be ϕ dependent, but for the sake of simplicity we will consider it a material constant, as a fit parameter.

Crucially, this finite relaxation strain scale for dilation completely modifies the way the stress profile is set up. In contrast to the diffusive dynamics for the steady-state rheology model (Eq. 4.9), inserting (4.10) in (4.4), yields a second-order nonlinear ODE for the stress (4.11) where the time derivative has dropped out. This fundamental difference has an important consequence regarding the early transient response during a change of J ext (or more generally to a change in boundary conditions). In the ODE (Eq. 4.11), information can travel infinitely fast and instantaneous changes of the particle stress can occur at γ = 0 + a finite distance from the wall. Said otherwise, in contrast with the steady-state rheology model, the Darcy-Reynolds model induces a non-local dynamics for the particle stress †.

∂ ŷ 1 ϕ R(ϕ) ∂ ŷ Pp + 1 γ 0 ϕ [ϕ -ϕ SS (J)] = 0,
For infinitesimal changes of the boundary conditions, i.e for δJ ≡ |J f -J i | ≪ 1, we can work out the "non-locality" length scale l over which Pp is affected by the change in J ext . Linearizing Eq. 4.11 at γ = 0 + for which ϕ = ϕ i , we obtain that δ Pp = Pp -Pp i satisfies

∂ 2 ŷ δ Pp +
δ Pp l2 = 0, (4.12) † While Eq. 4.11 formally looks similar to the stress dynamics postulated in non-local constitutive models [START_REF] Goyon | Spatial cooperativity in soft glassy flows[END_REF][START_REF] Kamrin | Nonlocal Constitutive Relation for Steady Granular Flow[END_REF], the similarity is only superficial, as here the transient constitutive law is local. The non-locality in Eq. 4.11 is induced by the combination of two processes, the local dilation coupled to the pore flow resisting it.

with l = -γ 0 RJ 2 i ϕ ′ SS (J i ) = γ 0 Kβ R(ϕ i )J β+1 i . (4.13)
The solution of Eq. 4.12, given the boundary conditions provided in Eq. 4.6, leads to the particle stress profile

Pp (ŷ) = 1 J i - 1 J i - 1 J f cosh(ŷ/ l) l-1 ϕ -1 i R(ϕ i ) -1 κ sinh( Ĥi / l) + cosh( Ĥi / l) , (4.14)
with Ĥi the thickness of the granular layer at γ = 0 + . A step change in J ext here leads to an instantaneous change in particle pressure on a depth of order l below the top wall, while deeper regions remain unchanged relative to the initial steady state at J i , see Fig. 6 (b). In particular, the particle pressure at the bottom of the layer P p (y = 0) follows:

P p (y = 0) P ext = J f J i + J f J i -1 1 l-1 ϕ -1 i R(ϕ i ) -1 κ sinh( Ĥi / l) + cosh( Ĥi / l) . (4.15)
The Darcy-Reynolds model thus predicts that, when Ĥi / l ⩽ 1, a step change in J ext instantly affects the particle pressure all the way to the bottom of the layer. Conversely, in the limit Ĥi / l → ∞, the steady-state rheology prediction P p (y = 0)/P ext = J f /J i is recovered. Moreover, we see that l may diverge at the jamming transition if we assume that γ 0 remains finite at ϕ c , with an asymptotic behavior l ∝ J .6 , using β ≈ 0.44 [START_REF] Athani | Transients in pressure-imposed shearing of dense granular suspensions[END_REF]. Investigating the particle pressure profile just after the step in small systems which are close to jamming, i.e. for which Ĥi / l ⩽ 1, should thus provide an unambiguous way to discriminate the two models.

-(1+β)/2 i ≈ J -0.7 i , which corresponds to l ∝ (ϕ c -ϕ) -(1+β)/2β ≈ (ϕ c -ϕ) -1
Note that in the limit case where J i = 0, i.e. for a granular layer initially at rest, we have ℓ → ∞. In such a case, the stress profiles is instantaneously affected over the entire height of the suspension, whatever the system size. We compute the stress profile in this limit case in Appendix B.

Model testing

In this section, we compare the results of the DEM simulations to the solution of the two models presented above ( §4). The predictions of the steady-state rheology model are fit-free, since the values of β = 0.44 and K = 0.67 are set from the steady-state rheological law ϕ SS (J) obtained from the DEM simulations [START_REF] Athani | Transients in pressure-imposed shearing of dense granular suspensions[END_REF]). In the Darcy-Reynolds model, we evaluated γ 0 , with γ 0 = 0.48 for dilation, and γ 0 = 0.28 for compaction separately.

We start by comparing the two models to the macroscopic quantities presented in Fig. 3. Overall, the long time evolution of the thickness H of the granular layer, the average particle stress Pp and the depth-averaged shear stress τ (see Fig. 3 b, c andd, respectively) are fairly well predicted by both the steady-state and Darcy-Reynolds models. There are some discrepancies on the shape of the curves, however, the agreement is pretty satisfactory given that there is no fitting parameters besides γ 0 . In particular, both models capture well the typical amplitude and the strain scale over which dilation/compaction occur. While from these results one could conclude that the refinement of the Darcy-Reynolds model over the steady-state rheology model is in practice unnecessary, the value of the particle stress at short times after the step reveals that the prediction of the steady-state rheology model is off by more than one order of magnitude. Conversely, the Darcy-Reynolds model achieves good quantitative predictions even at early strains after the step change of J ext , as highlighted in the inset of Fig. 3. Both models converge and start providing similar results only after a strain scale of order γ 0 . This result is expected since the strain scale required for stress diffusion to affect the stress profile on a length scale l is D(ϕ)/ l2 = γ 0 . It is also consistent with the dilatancy law Eq. 4.10, which states that locally, the steady-state rheology is recovered on a strain scale of order γ 0 . These macroscopic observations provide a first illustration of the fundamentally different dynamics for the particle stress in the two models. In the dilation case shown in Fig. 3, the nonlocal length scale given by the Darcy-Reynolds model is much smaller than the system height ( Ĥi / l ≈ 0.08). The change of boundary conditions at the top wall is thus felt instantaneously within the granular layer, an indication that Darcy-Reynolds coupling is key to capture the early transient dynamics of the granular layer.

To evidence this further, we compare in Fig. 5 the two continuum models to the particle pressure obtained with the DEM simulations at the bottom of the layer and just after the step for many combinations of J i , J f , and several system sizes. We have already seen that DEM results can differ from the steady-state rheology prediction at γ = 0 + (thick horizontal grey line), but this difference could be attributed to the finite strain at which DEM results are reported. Importantly, both models and DEM results are now compared at the same strain γ = 0.01. We find that in the steady-state rheology model, stress diffusion from the upper wall can lead to a finite decrease of the particle pressure at γ = 0.01. However, significant discrepancies remain (e.g. purple dashed lines). By contrast, the Darcy-Reynolds model (e.g. purple solid lines) provides surprisingly good quantitative agreement with the DEM simulations over the whole range of parameters.

In the insets of Fig. 5, we also show the full pressure profiles at γ = 0.01 for two examples. In the top profile, the step is performed far from jamming (i.e. starting from J i = 0.01). The "non-locality" length scale ℓ is thus much smaller than the system size and as a result, the steady-state and Darcy-Reynolds predictions are undistinguishable at the bottom of the cell. Conversely, in the bottom profile, the step is performed closer to jamming (i.e. starting from J i = 0.001). This time the "non-locality" length scale is larger and as a result the bottom particle pressure is instantaneously modified after the step, as predicted by the Darcy-Reynolds model. In both cases however, the Darcy-Reynolds model makes much better predictions for the pressure than the steady-state model near the top wall. Note that these particle stress profiles (along with similar plots at other strain values γ and other J f values) are used to fit the value of γ 0 , both for compaction (J f = 1 × 10 -4 ), yielding γ 0 = 0.28, and dilation (J f = 0.1) yielding γ 0 = 0.48. Those values are then kept constant when evaluating the Darcy-Reynolds model predictions.

Our simulation data unambiguously show that the instantaneous response in particle pressure can either be set by the fluid Darcy flow through the particle phase associated with dilation or compaction, or by the steady-state rheology, depending on the value of the parameters involved in the step in J ext . What decides which of these two scenarii dominates? We showed that in the case of the Darcy-Reynolds model, the linearized pressure ODE (4.12) predicts that the Darcy flow is associated with a length scale l corresponding to the typical length over which a dilation or compaction (and the associated Darcy flow) is instantaneously initiated when the boundary conditions on the top wall are modified. A natural expectation is thus that, even in the non-linear regime, the pressure response is dominated by the Darcy flow when the initial height of the system Ĥi is of order l. On the contrary, when Ĥi / l ≫ 1, we may expect that the steady rheology sets the pressure level in most of the system except a thin layer of height l below the top wall.

In Fig. 7 we therefore show as a function of Ĥi / l the rescaled bottom pressure ∆P = (P p (y = 0)/P ext -1) / (1 -J i /J f ) + 1. This rescaling is useful as for any step in J ext , ∆P is bounded such that 0 < ∆P < 1, with ∆P = 1 corresponding to the prediction of the steady rheology at γ = 0 + . We find that for all the combination of J i and J f investigated, data collapse on a master curve, where ∆P monotonically increases when plotted versus Ĥi / l, and saturates at 1 for Ĥi / l ≳ 10. This confirms that the boundary between the steady-rheology dominated and Darcy-Reynolds dominated regimes is set by l. In Fig. 7, we also show the prediction of the Darcy-Reynolds model (solid lines). The agreement with DEM simulation data is good, over the whole range of parameter investigated, that is for Ĥi / l varying over three orders of magnitude and J i varying over four orders of magnitude. The fact that the best fits of the particle stress profiles are obtained for two distinct values of γ 0 in dilation (γ 0 = 0.48) and compaction (γ 0 = 0.28) implies a J f -dependence of l which is ignored in the linearized problem, but probably cannot be in the cases considered here, for which J i -J f is of the same order of or even larger than J i .

Discussion & Conclusion

Through DEM simulations, we investigated the transient rheological behavior of a neutrally buoyant suspension in pressure imposed rheology configurations, subject to a sudden change in shear rate or external pressure. We compare these simulations with two competing continuum two-phase models: the standard Suspension Balance Model (SBM) which assumes the steady-state rheology to be valid at all times, and a "Darcy-Reynolds" model in which the volume fraction locally relaxes towards its steady-state value on a strain scale γ 0 , by analogy with Reynolds dilatancy in soil mechanics. This study shows that the early stress response of the suspension is not set by the steadystate rheological flow rules, but instead arises from the Darcy back-flow resulting from the geometrically imposed dilation rate of the granular phase (4.8). Before discussing the consequences of these results, let us recall that they were obtained for a dynamics which is essentially one-dimensional. This is because we considered a model setup where dilation occurs in only one direction and shear is homogeneous in the sample. Some results we obtain may be quantitatively modified in a less idealized setup, as for instance a non-uniform shear rate could locally modulate the Reynolds dilatancy. Moreover, in our configuration dilation occurs in the direction of the gradient of velocity. In other situations, dilation may also occur orthogonally to the velocity gradient direction. A model experimental realization of this phenomenon is the resuspension observed in a vertical cylindrical Couette cell [START_REF] Acrivos | Shear-induced resuspension in a couette device[END_REF]. It is however well-known from triaxial tests in soil mechanics that dilation is qualitatively similar in all directions [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF]. We also use a very simplified model for the fluid/particles interactions, through a constant (independent of ϕ) Darcy permeability coefficient embedding the complex flow between the particles. It would be interesting in future works to refine this description by performing more realistic three-dimensional hydrodynamic simulations which could also tackle more complex flow configurations (drop impact, avalanches). However, we believe that the minimal description we propose here captures the main physics involves during transients and that the qualitative picture which emerges from our results should remain. Indeed, our results have several important consequences and implications.

First, our study extends the domain of application of the Reynolds-like dilatancy law (4.10) proposed by [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF], which was introduced to describe the transient dilation/compaction of the avalanche of an initially dense/loose sediment under gravity. Here, we provide evidence that this law also applies for continuously sheared suspensions below ϕ c , when the flow parameters (external pressure, shear rate) are suddenly changed. The transient migration and stresses are quantitatively captured over a wide range of control parameter, J ∈ [10 -4 -10 -1 ], and system sizes, H/a ∈ [33 -140]. Our study thereby shows that the concepts of Reynolds dilatancy and shear-induced migration can be described within a unique framework provided by the Darcy-Reynolds model (4.11).

Second, our study reveals that after a sudden change of flow parameters, the stress levels inside the suspension can be very different than that predicted by the steady-state rheology (SBM model). More precisely, the Darcy-Reynolds model gives rise to a nonlocal length scale ℓ, which scales with the particle size and diverges algebraically at jamming (4.13). This length scale corresponds to the distance from the free boundary over which a finite rate of dilation/compaction occurs after a step change of flow conditions. In this region, the stress level is thus fixed, not by the steady-state rheology, but by the Darcy fluid pressure gradient resulting from this dilation/compaction rate. In practice, taking for instance a 3D granular layer with H = 20a, we estimate directly from our model (Eq. 4.13) that the Reynolds-dilatancy framework will dominate as soon as the distance to jamming (ϕ c -ϕ)/ϕ c is smaller than 7 % -a range often reached in dense suspensions (here we use for the hydrodynamic resistance the Richardson-Zaki expression for 3D spheres: R(ϕ) = (9/2)(1 -ϕ) 1-n with n = 5 and ϕ c = 0.58). Note that in our study, we consider steps in shear rate (viscous number J), where the SBM fails to predict the early stress response of the granular layer. Equally large discrepancies are expected under steps of applied shear stress (e.g. a rapid change of slope for submarine avalanches or a sudden change of pressure gradients in a pipe flow), for the temporal response of the suspension.

Accounting for this mechanism is key to predict the stress levels at early strains (γ ⩽ γ 0 ), which we quantify by the difference between the external imposed pressure and the particle pressure in the bulk, ∆P p . As summarized in the flow regime diagram sketched in Fig. 8, for a suspension of volume fraction ϕ, when the system size H is much larger than ℓ, ∆P p is well predicted by the steady-state rheology of the suspension. In a situation where the shear rate is suddenly modified from γ to γ + ∆ γ, the quasi-Newtonian steady-state rheology predicts ∆P steady-state p ∼ η f η n (ϕ)∆ γ. On the contrary, when H ⩽ ℓ, the transient pressure is set by the drag of the fluid on the particle phase generated by the Darcy flow coming from dilation or compaction, which spans the whole system. It therefore depends on the particle and system sizes. As the pressure gradient across the system scales as ∆P p /H, from Eq. 4.11 we can evaluate the Darcy-Reynolds scaling ∆P Darcy-Reynolds p ∼ η f γH 2 ∆ϕ/a 2 γ 0 , where ∆ϕ = |ϕ -ϕ SS (J f )| is the amount of dilation/compaction between the initial and final states. The boundary between the steady-rheology dominated and Darcy-Reynolds dominated regimes is thus set by the length scale ratio H/ℓ = O(1). This boundary also corresponds to the conditions for which the steady-state and Darcy-Reynolds prediction for the stress are equal, i.e. ∆P Darcy-Reynolds p /∆P steady-state p = O(1). Indeed, for a small flow perturbation ∆ γ ≈ γ∆J/J 2 while, ∆ϕ ≈ KβJ β-1 ∆J, which gives ∆P Darcy-Reynolds p /∆P steady-state p = H/ℓ. Our demonstration of the validity of the Reynolds-like dilatancy law (4.10) for flow below ϕ c and the identification of the nonlocal length scale ℓ for estimating the particle stress should help to better understand the short time dilation and migration dynamics of suspensions, as observed during impacts [START_REF] Nicolas | Spreading of a drop of neutrally buoyant suspension[END_REF][START_REF] Peters | Splashing onset in dense suspension droplets[END_REF][START_REF] Grishaev | Complex drop impact morphology[END_REF][START_REF] Schaarsberg | From splashing to bouncing: The influence of viscosity on the impact of suspension droplets on a solid surface[END_REF][START_REF] Boyer | Drop impact of shear thickening liquids[END_REF], submarine avalanches [START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF][START_REF] Topin | Collapse dynamics and runout of dense granular materials in a fluid[END_REF][START_REF] Iverson | Elementary theory of bed-sediment entrainment by debris flows and avalanches[END_REF][START_REF] Bougouin | Granular collapse in a fluid: Different flow regimes for an initially dense-packing[END_REF][START_REF] Montellà | A two-fluid model for immersed granular avalanches with dilatancy effects[END_REF], or unsteady two-phase flows in general [START_REF] Kulkarni | Particle-pressureinduced self-filtration in concentrated suspensions[END_REF]Snook et al. 2016;[START_REF] Saint-Michel | X-ray radiography of viscous resuspension[END_REF]d'Ambrosio et al. 2021). For instance, the Darcy-Reynolds model was successfully used to described the impact of a sphere on a suspension initially prepared above ϕ c [START_REF] Jerome | Unifying impacts in granular matter from quicksand to cornstarch[END_REF]. However, the behavior of the system below ϕ c could not be described by the quasi-static approach used in this study. Similarly, when prepared just below ϕ c , suspension droplets impacting a rigid plane observe large spreadings that cannot be captured using the steady-state viscosity of the suspension [START_REF] Jørgensen | Deformation upon impact of a concentrated suspension drop[END_REF]. For such small systems close to jamming, it is likely that H/ℓ < 1, for which the Darcy-Reynolds scaling applies. Note that these configurations are apparently under volume imposed conditions. However, the presence of a free surface allows a slight dilation of the granular network after the impact. We thus anticipate that in these configurations too, the level of stress right after the impact should originate, not from the steadystate rheology, but from the transient geometrical dilation of the granular phase and the associated pore pressure feedback effect.

Perhaps unexpectedly, Darcy-Reynolds physics could also play a role in uniform flows under constant volume. Indeed, even though in these flows the volume fraction is constant and uniform on average, it nonetheless shows small spatio-temporal fluctuations. The dynamics of these volume fraction fluctuations involve constant local dilation and compaction phenomena. We can then speculate that the typical extent of these fluctuations could be affected by the Darcy-Reynolds length scale ℓ.

Our study could also have implication to rationalize the transient dynamics of shear thickening suspensions. Shear thickening suspensions are characterized by a stressdependent critical packing fraction [START_REF] Seto | Discontinuous shear thickening of frictional hard-sphere suspensions[END_REF][START_REF] Wyart | Discontinuous shear thickening without inertia in dense non-brownian suspensions[END_REF]. Upon a sudden change of boundary conditions, an initially unjammed suspension may be driven above its maximum packing fraction and shear-jam. This shear-jammed regime exhibits fascinating transient features such as impact-activated solidification and traveling "jamming front" [START_REF] Waitukaitis | Impact-activated solidification of dense suspensions via dynamic jamming fronts[END_REF][START_REF] Han | High-speed ultrasound imaging in dense suspensions reveals impact-activated solidification due to dynamic shear jamming[END_REF]Han et al. , 2018a)), that cannot be described using steady-state rheological flow rules. Predicting the resistive stress in this shear-jammed regime is an important issue which is still largely unresolved. Most existing models assume a minimal description in which the medium jams in a finite strain [START_REF] Waitukaitis | Impact-activated solidification of dense suspensions via dynamic jamming fronts[END_REF][START_REF] Han | High-speed ultrasound imaging in dense suspensions reveals impact-activated solidification due to dynamic shear jamming[END_REF]Han et al. , 2018a) ) without considering two-phase flow coupling (see however [START_REF] Jerome | Unifying impacts in granular matter from quicksand to cornstarch[END_REF] and [START_REF] Brassard | Viscous-like forces control the impact response of shearthickening dense suspensions[END_REF]). However, in such configurations standing above ϕ c where the nonlocal length scale ℓ diverges, one expects the Darcy-Reynolds coupling to play a major role.

From the definition of the Darcy-Reynolds length scale in Eq. 4.13, one might expect that ℓ decreases with an increase of the hydrodynamic resistance R, see Eq. 4.13. However, for hard spheres this expectation is not met, as the increase of R (which comes from an increase of ϕ) comes together with a rapid increase of the viscosity, which has the opposite and actually dominant effect of increasing ℓ. In the case of deformable spheres the situation might be different, and increasing the volume fraction may or may not lead to a weakened Darcy-Reynolds physics. Separately, particle deformability also affects the constitutive law, which is then augmented with an elastic component. It is well-known that along with Darcy flow, an elastic constitutive law leads to a stress diffusion characteristic of poroelasticity (in the context of soil mechanics, this is known as Terzaghi's theory of soil consolidation [START_REF] Terzaghi | Theoretical Soil Mechanics[END_REF]). Thus, we can speculate that the effects of particle deformability both on the constitutive law and the Darcy resistance would lead to a reduction of the non-locality effects, and would instead favor a diffusive behavior of the stress. In our study however, we made sure that particle finite deformability does not affect our results, and that the particles can be assumed as infinitely rigid.

An important open question concerns the microscopic origin of the Reynolds-like dilatancy law (4.10) proposed by [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]. In the quasistatic regime, Reynolds dilatancy is usually explained from geometric/kinematic arguments through the introduction of a dilatancy angle [START_REF] Reynolds | on the dilatancy of media composed of rigid particles in contact. with experimental illustrations[END_REF][START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF]). More recently, this law was reinterpreted as a normal stress relaxation law in a medium of finite compressibility [START_REF] Bouchut | A two-phase two-layer model for fluidized granular flows with dilatancy effects[END_REF][START_REF] Lee | Two-phase modelling of submarine granular flows with shear-induced volume change and pore-pressure feedback[END_REF][START_REF] Montellà | A two-fluid model for immersed granular avalanches with dilatancy effects[END_REF], in line with the "Reynolds pressure" concept of [START_REF] Ren | Reynolds pressure and relaxation in a sheared granular system[END_REF]. Compared to the steady-state rheology model, the Reynolds-like dilatancy law requires the microstructure to reorganize in a strain scale γ 0 of order 1. Capturing the microstructure dynamics in a linear relaxation of a scalar variable is certainly an over simplification, as one probably needs to generally consider a tensorial descriptor of the microstructure (Chacko et al. 2018b;[START_REF] Gillissen | Modeling sphere suspension microstructure and stress[END_REF]. However, in the case of a simple shear in a direction constant in time, a scalar relaxation has already proven useful to describe transient responses [START_REF] Mari | Nonmonotonic flow curves of shear thickening suspensions[END_REF]Chacko et al. 2018a;Han et al. 2018b).

It is tempting to relate the value of γ 0 to the typical strain scale γ trans for transients in imposed-volume rheology, observed e.g. in startup flows or shear reversal experiments Gadala-Maria & [START_REF] Gadala-Maria | Shear-Induced Structure in a Concentrated Suspension of Solid Spheres[END_REF], or to the strain scale γ rev of the reversibleirreversible transition [START_REF] Pine | Chaos and threshold for irreversibility in sheared suspensions[END_REF]. These strain scales are all of order 1, provided that they are measured sufficiently far from jamming. A careful examination however concludes that the relation, if it exists, is subtle. Indeed, both γ trans and γ rev decrease upon increasing the volume fraction, and most probably vanish at jamming (see [START_REF] Peters | Rheology of non-Brownian suspensions of rough frictional particles under shear reversal: A numerical study[END_REF] for γ trans and [START_REF] Corté | Self-Organized Criticality in Sheared Suspensions[END_REF]; [START_REF] Pham | Origin of critical strain amplitude in periodically sheared suspensions[END_REF] for γ rev ). This is not the case for γ 0 , which we measure independent of the volume fraction. Moreover, γ 0 is by nature tied to volumetric change, and as such is not observable under volume-imposed conditions. This is obvious when one considers the same protocol than studied here but under imposed-volume conditions: the response to a step increase (or decrease) of the shear rate by some proportion is an immediate increase (or decrease) of the shear stress by the same proportion, which is fundamentally a consequence of the linearity of Stokes equation. This is in contrast to shear under imposed pressure, for which we show that there is a transient in the stress response occuring over a finite strain scale γ 0 .

It is worth noting that in [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF], for dilation experiments starting from above the jamming point (ϕ > ϕ c ), Reynolds dilatancy follows γ 0 ≈ 0.3. This value is intriguingly close to the values of γ 0 we find for our ϕ < ϕ c conditions. It is thus tempting to conjecture that the value for γ 0 is independent of ϕ, or at least only weakly dependent on ϕ, around the jamming transition. This would mean that the classical Reynolds dilatancy of the critical state theory is just a limiting case of a more general dilatancy law valid over a wide range of volume fractions. forces and the externally applied force for wall particles. Contact forces are modelled using a system of springs and dashpots [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. Using the notation F C ≡ (f C,1 , . . . , f C,N , t C,1 , . . . , t C,N ) for the N d(d + 1)/2 vector of contact forces and torques on each particle (in spatial dimension d), we have

F C = -R C FU • U + F C,S (A 1)
with U (resp. Ω) the vector of velocities (resp. angular velocities) for each particle, R C FU the resistance matrix associated to dashpots, and F C,S the part of contact forces coming from springs. For hydrodynamic forces and torques F H , we have [START_REF] Jeffrey | The calculation of the low Reynolds number resistance functions for two unequal spheres[END_REF])

F H = -R H FU • U ′ + R FE : E ∞ (A 2)
with U ′ = U -U ∞ the vector of non-affine velocities and angular velocities for each particle, where U ∞ are the background velocities and angular velocities evaluated at the particles centers, and E ∞ the symmetric part of the imposed velocity gradient ∇U ∞ . The resistance matrices R H FU and R FE contain Stokes drag and regularized lubrication forces at leading order in particle separation [START_REF] Mari | Shear thickening, frictionless and frictional rheologies in non-brownian suspensions[END_REF].

The equation of motion is thus

F H + F C + F Ext = 0 , ( A 3)
with F Ext the external force/torque applied on each particle. Of course, F Ext takes nonzero values only for the wall particles. As the simulated system (see Fig. 2) comprises of frozen wall and bulk particles suspended in the fluid, we separate the total resistance matrix as

R C FU + R H FU ≡ R FU = R bb FU R bw FU R wb FU R ww FU (A 4)
The matrices R bb FU , R bw FU , R wb FU and R ww FU indicate the hydrodynamic resistance matrices including bulk-bulk, bulk-wall, wall-bulk and wall-wall interactions respectively. Similarly, the non-affine velocities have also been separated into bulk particle non-affine velocities U b -U ∞ and wall particle non-affine velocities U w -U ∞ ,

U ′ = U ′b U ′w . (A 5)
Note that the wall particles do not follow the affine flow U ∞ for two reasons. First, their velocity have a vertical component when the system dilates or contracts. Second, even the horizontal component of the velocity of a wall particle, say particle i, located at height y i , does not match U ∞ (y i ) because the wall moves rigidly and is made of particles located at slightly different heights. We can then solve the equations of motion, Eq. A 3 for the bulk velocities, to get

U ′b = R bb FU -1 • K b -R bw FU • U ′w (A 6) with K b K w = R FE : E ∞ + F C,S -R C FU • U ∞ (A 7)
We further decompose the wall non-affine velocity in horizontal and vertical components

U ′w = U ′w h + v y Y w
, where Y w is the vector corresponding to a unit non-affine vertical velocity for particles belonging to the upper wall, and vanishing non-affine velocity for particles of the lower wall. Thus the scalar v y = ∂ t H is the upper wall vertical speed.

Using ∂ ŷ Pp as an integrating factor (that is, multiplying by ∂ ŷ Pp and integrating on ŷ), we get

1 2 ∂ ŷ Pp 2 = - R(ϕ c )K γ 0 (1 -β) P -β+1 p + C 2 , ( C 2)
with C an integration constant. Remembering that from a jammed configuration we expect dilation, that is ∂ ŷ Pp < 0, leads to

∂ ŷ Pp = --2 R(ϕ c )K γ 0 (1 -β) P -β+1 p + C. ( C 3)
To go further we assume β = 1/2, close to the actual measured value, and separating variables we can integrate to get

ŷ = √ 2γ 2 0 3 R(ϕ c ) 2 K 2 R(ϕ c )K γ 0 Pp + C -2 R(ϕ c )K γ 0 Pp + C + C ′ , ( C 4)
with a new constant of integration C ′ . Because of the no-flux boundary condition at the bottom wall, ∂ ŷ Pp (0) = 0, which with Eq. C 3 imposes that -2 R(ϕc)K γ0

Pp + C = 0 for ŷ = 0. We therefore have C ′ = 0, as well as C = 2 R(ϕ c )K Pp (0)/γ 0 > 0. Introducing r = Pp and A = R(ϕ c )K/γ 0 , we then get a cubic equation for r,

0 = 2A 3 r 3 + Cr 2 - C 3 3A 2 + 3ŷ 2 2 ≡ Q ŷ (r). ( C 5)
For each height ŷ, one solution r(ŷ) of this polynomial gives the square root of the particle pressure. The discriminant of the polynomial is

∆(ŷ) = - 16A 2 3 ŷ2 - 2C 3 9A 2 + ŷ2 . (C 6)
This discriminant defines two different ranges for ŷ: for ŷ ∈]0, ŷ0 [, with ŷ0 = √ 2C 3 /3A, ∆(ŷ) > 0, while for ŷ ∈]ŷ 0 , +∞[, ∆(ŷ) < 0. At two points, ŷ = 0 and ŷ = ŷ0 , the discriminant vanishes. Let us focus on the ŷ ∈]0, ŷ0 [ range. A strictly positive discriminant corresponds to having three real solutions for this region. Furthermore, the derivative Q ′ ŷ (r) = 2r(Ar + C) vanishes for r = 0 and r = -C/A < 0. This implies that for ŷ ∈]0, ŷ0 [, one of the three solutions must be ∈ [-∞, -C/A], one must be ∈ [-C/A, 0] and the last one must be in [0, ∞]. As r = P p must be strictly positive, the latter is the correct solution. For ŷ = ŷ0 , and Q ŷ becomes

Q ŷ=ŷ0 (r) = r 2 2A 3 r + C , ( C 7)
so that there are only two roots (as expected from the fact that ∆(ŷ 0 ) = 0), a degenerate one in r = 0, and another one in r = -3C/2A < 0. Again, r > 0 up to the top wall, which implies that Ĥ < ŷ0 . As, ŷ0 = √ 2C 3 /3A, this gives us the condition C > (3A Ĥ/ √ 2) 2/3 . The three roots of Q ŷ are, with k = 0, 1, 2

r k = - C 2A 1 + 2 cos 1 3 φ 3A C 3/2 ŷ + (2k + 1)π , ( C 8)
with φ such that e iφ(x) = 1 -x 2 -ix √ 2 -x 2 . We can identify the correct root with the . Particle pressure at the bottom wall Pp(y = 0) at γ = 0 + scaled by the imposed pressure at the top wall Pext, as a function of the system height Ĥ = H/a, when the initial configuration is jammed (Ji = 0), from Eq. C 9 . The pressure scales as Ĥ4/3 in the large Ĥ limit.

case ŷ = 0, for which we know that there are two degenerate negative solutions and the positive one we are interested in. Indeed, as φ(0) = 0, we have that r 0 = r 2 = -C A for ŷ = 0, while r 1 = C 2A , which implies that r 1 is the root we are looking for. As a result, replacing A by its value in terms of the problem parameters, we have the stress profile at γ = 0

Pp (ŷ) = γ 2 0 C 2 4 R(ϕ c ) 2 K 2 1 -2 cos 1 3 φ 3 R(ϕ c )K γ 0 C 3/2 ŷ 2 . ( C 9)
We can then pick C to satisfy the top wall boundary condition defined in Eq. 4.6, but this must be done numerically as there is no mathematically closed expression for C. A few important things can however be inferred from Eq. C 9. First, although a cosine appears in the r.h.s., it of course does not mean that the stress profile is non-monotonic. Indeed, as noted earlier C is bounded from below by (3A Ĥ/ √ 2) 2/3 , which means that the argument of φ in Eq. C 9 is satisfying 3 R(ϕc)K γ0C 3/2 ŷ < √ 2. It is easy to show that φ(0) = 0 and φ( √ 2) = π, so that the cosine in r.h.s of Eq. C 9 is decreasing from the bottom wall to the top wall, its argument spanning at most the interval [0, π/3]. Moreover, the stress at the top of the cell must be different from the stress at the bottom, and from the argument of the cosine that it implies that C must scale as Ĥ2/3 for large Ĥ.

Second, we expect from the divergence of l when ϕ i = ϕ c that the stress profile is instantaneously modified at all points in the system at γ = 0. The fact that a branch of cosine taking ŷ in its argument appears in Eq. C 9 is pointing towards such a scenario, but because of the nonlinear behavior of φ it is not immediately apparent that the profile is not flat in the bottom of the cell. We know that ∂ ŷ Pp (0) = 0 by construction (which can easily be verified on Eq. C 9), so we have to evaluate the second derivative at the bottom wall to be informed about the flatness of the profile there. A quick calculation, using that φ ′ (0) = -√ 2, shows that the curvature of the stress profile at the bottom wall is

∂ 2 ŷ Pp (0) = - 2 C ∝ Ĥ-2/3 , ( C 10)
which is to be compared with the ∝ exp(-Ĥ/ l) for ϕ i < ϕ c . This weak scaling implies that the stress at the bottom will be height dependent, which is easy to verify, as (C 11) here again in stark contrast to the case J i > 0, for which the stress at the bottom is initially converging to -1/J i exponentially in Ĥ/ l. This approach to the asymptotic behavior is shown in Fig. 9.

Pp (0) = γ 2 0 C 2 4 R(ϕ c ) 2 K 2 ∝ Ĥ4/3 ,

C.2. Flowing initial state

If the initial state is a steady state under viscous number J i , we again can derive the stress profile at γ = 0 for a rheology given by Eq. 4.8 in the case β = 1/2. Injecting Eq. 4.8 in Eq. 4.11, with a uniform solid fraction ϕ i = ϕ SS (J i ), we get

∂ 2 ŷ Pp + γ -1 0 R(ϕ i )K   -J i + 1 Pp   = 0. (C 12)
Using ∂ ŷ Pp as an integrating factor, we get

1 2 ∂ ŷ Pp 2 = γ -1 0 R(ϕ i )K 2 Pp (0) -Pp + J i Pp -Pp (0) , ( C 13)
where we used the bottom wall boundary condition ∂ ŷ Pp (0) = 0 and we recall that Pp (0) is the particle pressure at the bottom wall. This leads to ∂ ŷ Pp = ± 2γ -1 0 R(ϕ i )K 2 Pp (0) -Pp + J i Pp -Pp (0) .

(C 14)

Here the choice of the r.h.s. sign depends on the change of J ext . A dilation, with J f > J i , implies ∂ ŷ Pp ⩽ 0, whereas a compaction, with J f < J i , implies ∂ ŷ Pp ⩾ 0. Introducing Finally, Pp (0) is set by the top wall boundary condition, Eq. 4.6, which makes Eq. C 16 an implicit but complete solution for the stress profile.

Figure 1 .

 1 Figure 1. (a) Examples of configurations involving transient dilation/migration in granular suspensions. From left to right: a submarine avalanche with a slope-break; pipe flow with a sudden change of section; impact of a solid object on a grain/fluid mixture; impact of a drop of concentrated suspension. (b) Canonical configuration investigated in the present study. A suspension of neutrally buoyant particles is homogeneously sheared under an imposed external pressure Pext and shear rate γ. A step change in Jext = η f γ/Pext is then applied by suddenly varying Pext or γ. Upon the step increase in Jext shown in the Figure, the particle skeleton dilates and yields: a transient fluid Darcy back-flow (blue arrows), a transient decrease of pore pressure P f and a transient increase of particle stress Pp. The red segments highlight the increase of the particle stress.

Figure 2 .

 2 Figure2. a) Pressure-imposed rheology DEM simulations. The system is composed of mobile particles (in grey) which are immersed in a suspending Newtonian fluid and placed between the top and bottom walls constructed with frozen particles (in blue). The control parameters are the external stress Pext on the top-wall and the shear rate γ imposed to the entire system. In this pressure-imposed rheology configuration, the top wall can move vertically, allowing the volume fraction to adjust freely owing to the value of the imposed Jext = η f γ/Pext. The fluid is not explicitly modeled but accounted for through lubrication and drag forces. b) The drag force fD exerted by the fluid on particle i is proportional to the difference of between the local fluid and particle velocities. This force ensures that grains are advected by a homogeneous background shear flow in the x-direction, and models the pore pressure feedback (or interphase drag coupling) occurring during dilation (compaction) of the granular skeleton in the y-direction.

Figure 3 .

 3 Figure 3. Macroscopic response of the granular layer subject to a step change in Jext applied at γ = 0. In blue, a step increase from Ji = 1 × 10 -4 to J f = 1 × 10 -1 , corresponding to a dilation from ϕ ≈ 0.8 to ϕ ≈ 0.57, and in red, a step decrease from Ji = 8 × 10 -2 to J f = 1 × 10 -4 , corresponding to a compaction from ϕ ≈ 0.61 to ϕ ≈ 0.8. Evolution with strain γ of (a) the imposed Jext, (b) the layer thickness height H/a, (c) the normalized depth-averaged particle pressure Pp/Pext, (d) the normalized depth-averaged shear stress τ /τSS at the top wall, where τSS = µ(Jext)Pext when the system is at steady state, and (e) of the early transient particle stress response in log-log representation in the compaction case. All results are averaged over 10 different realizations and were obtained for the system size H/a = 33 at Jext = 1×10 -4 . Red and blue lines: DEM simulations, Dashed orange lines: model based on steady-state rheology (orange star: prediction at γ = 0 + ), Solid green lines: Darcy-Reynolds model (green star: prediction at γ = 0 + ).

Figure 4 .

 4 Figure 4. Successive normalized particle pressure Pp/Pext and volume fraction ϕ profiles along the y-direction after a step change in Jext, from Ji = 1.6 × 10 -3 to J f = 1 × 10 -1 , corresponding to a dilation from ϕ ≈ 0.77 to ϕ ≈ 0.57 (respectively red and black dashed lines).

Figure 5 .

 5 Figure 5. Normalized particle pressure at the bottom of the cell Pp(y = 0)Ji/PextJ f just after the step change in Jext, for various combinations of Ji and J f values corresponding to values of ϕ in the range 0.56 -0.8, and different system sizes H/a. Markers: DEM simulations evaluated at γ = 0.01. Thick grey line: prediction of the steady-state rheology at γ = 0 + , Dashed lines: prediction of the steady-state rheology at γ = 0.01 (note that they are falling on top of the predictions of the steady-state rheology at γ = 0 + under compaction, for J f /Ji < 1), Solid lines: prediction of the Darcy-Reynolds model at γ = 0.01. DEM simulation results are averaged over at least 5 different realizations when Ji is fixed, and at least 40 different realizations when J f is fixed. Insets: two examples of particle pressure profiles along the system height.

Figure 6 .

 6 Figure 6. Expected normalized particle pressure profile at γ = 0 + for a) the steady-state rheology model and b) the Darcy-Reynolds model for different values of H/ℓ.

Figure 7 .

 7 Figure 7. The rescaled bottom particle pressure ∆P = (Pp(y = 0)/Pext -1) / (1 -Ji/J f ) + 1 plotted versus Ĥi/ l at γ = 0.01 shows that the boundary between the steady-state rheology and the Darcy-Reynolds dominated regime is set by the nonlocal length scale l, and that the steady-state rheology only applies when Ĥi/ l ≫ 1. Symbols: DEM simulations, Solid lines: Darcy-Reynolds model, Grey dashed line: steady-state rheology model at γ = 0 + . Steps in Jext considered here are from Ji values in the range 1×10 -4 to 1×10 -1 (from ϕ ≈ 0.8 to ϕ ≈ 0.57) to J f = 1 × 10 -1 , corresponding to ϕ ≈ 0.57, for the red data, and to J f = 1 × 10 -3 , corresponding to ϕ ≈ 0.78, for the blue data.

Figure 8 .

 8 Figure 8. Domain of application of the steady-state rheology model and adequate scaling for the particle stress for the early response of a suspension of initial volume fraction ϕ and typical size H, subject to a sudden change of flow condition ∆J. Here ∆ γ is the shear rate jump and ∆ϕ = |ϕ -ϕSS(J f )|.

Figure 9

 9 Figure9. Particle pressure at the bottom wall Pp(y = 0) at γ = 0 + scaled by the imposed pressure at the top wall Pext, as a function of the system height Ĥ = H/a, when the initial configuration is jammed (Ji = 0), from Eq. C 9 . The pressure scales as Ĥ4/3 in the large Ĥ limit.

  g( Pp ) = 2 Pp (0) -Pp + J i Pp -Pp (0) , (C 15)we can integrate Eq. C 14 by separation of variables,

Acknowledgements

This work was funded by ANR ScienceFriction (ANR-18-CE30-0024).

Declaration of interests

The authors report no conflict of interest.

Appendix A. Numerical method

As described in Sec. 2, the granular suspension considered in this study is neutrally buoyant and we work in vanishing Stokes number and Reynolds number regime; consequently, the equation of motion is just the force balance between hydrodynamic forces and contact forces for bulk particles, and between hydrodynamic forces, contact Injecting Eq. A 6 in the wall part of Eq. A 3, and using Eq. A 5, we get

Taking the dot product with Y w , and using Y w • F Ext = P ext l x we get the vertical wall velocity

Imposing that the vertical component of the upper wall velocity is v y ensures that the total force on the upper wall vanishes.

Appendix B. Two-phases equations: effect of the vertical fluid velocity

In this appendix, we show that taking into account the vertical fluid velocity in the two-phases equations [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid[END_REF][START_REF] Jackson | The dynamics of fluidized particles[END_REF] only changes the Darcy pressure term by a factor (1 -ϕ) -1 . Mass conservation for the whole suspension and the particle phase reads

while the vertical momentum balance for the whole suspension and the particle phase reads

where P f is the opposite of the normal fluid stress and η f a 2 ϕR(ϕ) (u f -u p ) is the interphase drag force with R(ϕ) the hydrodynamic resistance. Since the bottom boundary is fixed and impermeable, the whole suspension vertical velocity is zero, i.e. u = ϕu p +(1-ϕ)u f = 0, yielding u f -u p = -u p /(1 -ϕ), which enables to express the particle pressure gradient (Eq. B 4) only in terms of u p :

This equation, which includes the effect of the vertical fluid velocity, differs from Eq. 4.2 only by a factor (1 -ϕ) -1 .

Appendix C. Stress profile from uniform state

C.1. Jammed initial state

As noted in Sec. 4, within the Darcy-Reynolds model upon change of J ext at γ = 0, the viscous number is instantaneously changing in the part of the suspension just below the top wall, on a length scale l. When the initial state is at J i = 0, l however diverges, and the response is qualitatively different from the one exposed in Eq. 4.14. We can illustrate this as Eq. 4.11 has an analytically tractable solution for an initial profile ϕ i = ϕ c when β = 1/2.

From Eq. 4.11, first with an initially arbitrary β, for ϕ i = ϕ c with Eq. 4.8 the stress satisfies


domain: info math sdv.bibs

text:
Coralie Muller
email: coralie.muller@inria.frarie.wortsman-zurich@inria.fr

Arie Wortsman

Pablo Andres

Ugaldes Salas

Clémence Frioux

Simon Labarthe

metabolites demand for inputs of metabolites export of metabolites between compartments We extracted and cleaned reactions occurring in colonocytes (epithelial cell

Keywords: curation, harmonization, optimization

from a metabolic network of the whole human body [2].

We want to model the dynamics of a host-microbiota-pathogen system through the dFBA method. Dynamic FBA (dFBA) integrates the fluxes calculated in FBA into a system of ordinary differential equations (ODEs), enabling us to model the temporal dimension of metabolic simulations. In order to speed up the computation of the FBA models needed for a thorough numerical exploration, the metabolic models will be approximated by a metamodel method [1].

Dynamic simulation fit better to the biological reality

Bioinformatics Mathematics

We model the host-microbiota-pathogen system by gathering three metabolic models : A long manual curation process was performed in order to clean up and make the model usable. It consists of five main steps:

Step 1 Step 2

Step 5 In order to make the simulation with the 3 models, we have harmonized the networks between them on two points:

Identify blocked reactions

Different Ids for the same metabolite and differenciate name for compartments between networks For the simulation, we consider that the system works by maximizing biomass production. Also, we assume that the system is in a stationary state. This way, can model the system by solving the linear optimization problem:

We assume all three of the system's compartment's follow the FBA assumptions. We also define the function that maps constraints to the optimal flux according to FBA [5] by: Our model assumes that FBA holds for each instant [5]. Then, for simulating it for a period of time, it is necessary to solve FBA on each compartment. After this timestep, an ODE system models the fluxes between different compartments, until the next step. . . How does the equations of the ODE system look like?. For the concentration of a component, we have the following form :

Luminal compartment:

Epithelial compartment:

Each time iteration of the model implies an optimization step in a high dimensional space, which means that the simulation has a very high computational cost.

There are still reactions in the colonocyte network that do no behave as expected. Mainly, the consumption of butyrate in the network has some problems we've been looking at. Solving this is a work in progress. Modelling the dynamics of Salmonella infection in the gut at the bacterial and host levels

Blood Port Compartment

[bpC]

  Figure 1:

Fix

  Are there ways to solve this?Is there a big trade-off between computation cost and precision?We want to simulate two main cases to validate our model of but by commensal bacteria -Accumulation of O2 and transfert to [luLI] Øyås and Jörg Stelling. Genome-scale metabolic networks in time and space. Current Opinion in Systems Biology, 8:51-8, 2018.


Coralie Muller
email: coralie.muller@inria.fr

Arie Wortsman

Pablo Andres

Ugaldes Salas

Clémence Frioux

Simon Labarthe

SK sink for outputs of metabolites demand for inputs of metabolites export of metabolites between compartments

Colonocyte metabolic network curation

Gut microbiota -Metabolic network -Biological system -System of ODE -Metamodelling -Numerical metabolic modelling

References:

Three main challenges: curation, harmonization, optimization

We want to model the dynamics of a host-microbiota-pathogen system using the Dynamic Flux Balance Analysis method. dFBA integrates the fluxes calculated in FBA into a system of ordinary differential equations (ODEs), enabling us to model the temporal dimension of metabolic simulations. In order to speed up the computation of the FBA models needed for a thorough numerical exploration, the metabolic models will be approximated by a metamodel method [ 1]. We model the host-microbiota-pathogen system by gathering three metabolic models : Faecalibacterium prausnitzii [3] as a representative butyrate producing commensal bacterium, Salmonella Enterica serovar Thiphimurium [4] for the enteric pathogen , these models previously used in CEMRACS[1] project, and a host genome scale metabolic network of human [5] called iHsa, to simulate the epithelial cell.

Dynamic simulation fit better to the biological reality

A manual curation process was performed in order to clean up five main steps:

Step 1 Step 2

Step 5

Step 3

Metabolic network of colonocyte

In order to make the simulation with the 3 models, we have harmonized the networks between them on two points:

Different Ids for the same metabolite and differenciate name for compartments between networks Fix fluxes for essential reactions in the right compartments for the exchanges

butyrate[c] but[cFprau] but[luLI] BUty[c] but[cSalmo]

Harmonization of the models

Next steps

Simulation: FBA -dFBA

For the simulation, we consider that the system works by maximizing biomass production. Also, we assume that the system is in a stationary state. This way, can model the system by solving the linear optimization problem What does ODE system look like?

Each time iteration of the model implies an optimization step in a high dimensional space, which means that the simulation has a very high computational cost. Simulations involving the uptake of O2 and butyrate on ATP synthase activity are encouraging. The next steps of the project will be:

Our model assumes that FBA holds for each instant [2]. Then, for simulating it for a period of time, it is necessary to solve FBA on each compartment. After this timestep, an ODE system models the fluxes between different compartments, until the next step...

Beta oxydation

Butyrate

Ove Øyås and Jörg Satellin . Genome-scale metabolic networks in time and space. Current Opinion in Systems Biology, 8:51-8, 2018.

FBA on each compartment

FBA on each compartment

Simulation of growth in the human

Effect of but and O 2 on energy production

We performed simulations with a custom medium to observe biomass growth as a function of O2 and butyrate uptake.

We observed a dependence of biomass growth on O2 uptake but not on butyrate uptake.

The ATP synthase is the main energy production reaction of the mitochondria and reflects the energy production [6] of the cell. With this activity, we want to monitor the Salmonella infection.

Lumen Compartment [luC]

Test several human models iHsa model ? -Enterocyte model ?

The electron transport chain (ETC) is the source of energy production of the cell, it is located in the mitochondria and allows to produce at the end of the chain ATP. To function, ETC needs in part O2 and FADH2 to enable electron transport. FADH2 is a product of the Krebs cycle reactions (TCA) whose main reactant is butyrate in enterocytes. Butyrate is produced by the beta oxydation of fatty acids in the cell.

TCA FADH2 Electron transport chain ATP O2

These results are not surprinsing because biomass includes many are not directed related to the energy metabolism and therefore do not depend on butyrate. with dFBA to validate the reduced of butyrate uptake lead to O2 accumulation in the lumen [luLI].

Learn the behaviour of the system build a metamodel able to efficiently approximate the dynamics.

  Figure 1:

[ 3 ]

 3 El-Semman et al. "Genome-scale metabolic reconstructions of Bifidobacterium adolescentis L2-32 and Faecalibacterium prausnitzii A2-165 and their interaction." BMC systems biology 8.1 (2014): 1-11

  We want to simulate two main cases to validate our model but

  Modelling the dynamics of Salmonella infection in the gut at the bacterial and host levels Frioux et al. Accelerating metabolic models evaluation with statistical metamodels: applicationto Salmonella infection models. working paper or preprint, April 2022.

  transcriptome based reconstruction of salmonella human integrated metabolic network to screen potential drug targets. PLOS ONE, 17, e0268889 https://doi.org/10.1371/JOURNAL.PONE.0268889

1The butyrate/O2 uptake has an impact on the ATP synthase activity, when the other uptake is stable.The effect of O2 uptake has a greater impact on ATP synthase activity than butyrate uptake.