L Goren

E Aharonov

D Sparks

R Toussaint

The Mechanical Coupling of Fluid-Filled Granular Material Under Shear

The coupled mechanics of fluid-filled granular media controls the physics of many Earth systems such as saturated soils, fault gouge, and landslide shear zones. It is well established that when the pore fluid pressure rises, the shear resistance of fluid-filled granular systems decreases, and as a result catastrophic events such as soil liquefaction, earthquakes, and accelerating landslides may be triggered. Alternatively, when the pore pressure drops, the shear resistance of these geosystems increases. Despite the great importance of the coupled mechanics of grains-fluid systems, the basic physics that controls this coupling is far from understood. Fundamental questions that need to be addressed are what are the processes that control pore fluid pressurization and depressurization in response to deformation of the granular skeleton? and how do variations of pore pressure affect the mechanical strength of the grains skeleton? To answer these questions, a formulation for the pore fluid pressure and flow is developed from mass and momentum conservation, and is coupled with a granular dynamics algorithm that solves the grain dynamics, to form a fully coupled model. The pore fluid formulation reveals that the evolution of pore pressure obeys a viscoelastic rheology in response to pore space variations. Elastic-like behavior dominates with undrained conditions and leads to a linear relation between pore pressure and overall volumetric strain. Viscous-like behavior dominates under well drained conditions and leads to a linear relation between pore pressure and volumetric strain rate. Numerical simulations reveal the possibility of liquefaction under drained and initially over-compacted conditions, which were often believed to be resistant to liquefaction. Under such conditions liquefaction occurs during short compactive 1 phases that punctuate the overall dilative trend. In addition, the more established generation of elevated pore pressure under undrained compactive conditions is observed. Simulations also show that during liquefaction events stress chains are detached, the external load becomes completely supported by the pressurized pore fluid, and shear resistance vanishes.

Introduction

Fluid-filled granular media are ubiquitous in the Earth, mostly in the upper crust. Soils, fault gouge and landslide shear zones located below the water table are geosystems that are best described as fluid-filled granular media. Geometrically, such materials are composed of a 3D skeleton built out of contacting grains, whose exact configuration defines the pore space where fluid may reside. The mechanical strength of such systems is a function of both phases: the pore fluid and the grains. Already at the beginning of the 20th century [e.g. [START_REF] Terzaghi | Theoretical Soil Mechanics[END_REF], Terzaghi understood that it is not the stress which controls the solid-fluid system strength, but instead a quantity termed the 'effective stress':

σ ′ ij = σ ij -δ ij P, (1)
where σ ij is the applied stress tensor, P is the pressure experienced by the fluid within the pores of a granular or porous material, δ ij is Kronecker's delta, and σ ′ ij is the effective stress tensor. Later on, the effective stress was found to depend also on the properties of the bulk material composing the grains, the properties of the granular skeleton, and the properties of the pore fluid. These dependencies were formulated using an effective stress coefficient, 0 < α ≤ 1, that multiplies P , where α was found to be different for different physical quantities [START_REF] Wang | Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology[END_REF][START_REF] Pride | Relationships between seismic and hydrological properties[END_REF]. Still, it was shown that generally α is very close to one and when the material composing the solid matrix is incompressible relative to the pore fluid, α = 1, and Terzaghi's formulation is valid [START_REF] Nur | An exact effective stress law for elastic deformation of rock with fluids[END_REF][START_REF] Robin | Note on effective pressure[END_REF][START_REF] Wang | Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology[END_REF][START_REF] Pride | Relationships between seismic and hydrological properties[END_REF].

The most important consequence of the law of effective stress, equation (1), is that the shear stress, τ , required to shear the system is not a function of the normal stress as in Coulomb's law, but a function of the effective stress instead:

τ = µ(σ n -P ), (2)
where σ n is the total applied stress normal to a shear surface, µ is the surface friction coefficient, and the cohesion is neglected. It is immediately observed that if P increases to be equal to σ n then the system completely loses its shear strength, while if P decreases the system has higher resistance to shear. Therefore, the pore fluid pressure is of critical importance in the mechanics of fluid-filled granular (and porous) systems undergoing shear. What is less clear is (1) what is the physics behind the pore pressure control over the shear strength (equation ( 1) and ( 2))?, and ( 2) what are the mechanical processes that control the evolution of pore pressure? These two questions are at the heart of this work, and their answer lays the foundations for predicting the coupled mechanics of grains and pore fluid. Characterization and understanding of shear deformation is of particular importance in geosystems: Any differential forcing or gradients in material properties may lead to shear deformation. For example: the passage of seismic shear waves through a soil column induces shear deformation of the soil; tectonic loading accumulated in a fault zone will eventually lead to shear sliding of the fault; and gravitational forces may lead to landslides that shear along a confined zone at their base. In these examples, the presence of pore fluid changes the onset and dynamics of shear deformation in response to forcing, because the fluid pressure affects the resistance to shear, in accordance with equation (2). Next, we review in details the role of pore pressure in soil liquefaction and in motion along fault zones and landslide shear zones.

Soil liquefaction. In the process of earthquake induced soil liquefaction, the passage of seismic waves deforms the granular matrix and the fluid in such a way as to lead to pore pressure rise [START_REF] Das | Principles of Soil Mechanics[END_REF][START_REF] Kramer | Geotechnical earthquake engineering[END_REF]. The consequent reduction of shear resistance causes the granular system, which under normal conditions behaves like a solid that resists shear, to flow as a fluid. Once liquefied, soils can no longer support the infrastructure and a catastrophic collapse of buildings, roads, bridges and other structures may take place (e.g., damage during earthquakes at Niigata, 1964, [Kawakami and[START_REF] Kawakami | Damage to the ground and earth structures by the Niigata earthquake of june 16,1964[END_REF], or Izmit, 1999, [START_REF] Cetin | Seismically induced landslide at Degirmendere Nose, Izmit bay during Kocaeli (Izmit)-Turkey earthquake[END_REF]).

The coupled physics controlling soil liquefaction is not completely understood. The classical approach suggests that cyclic loading (such as the passage of shear waves during an earthquake) leads to irreversible collapse of initially under-consolidated pore volume. When drainage is poor, 'the tendency for volume reduction' of the loose granular skeleton may lead to pore fluid pressurization and to liquefaction [START_REF] Sawicki | Developments in modeling liquefaction of granular soils, caused by cyclic loads[END_REF]. This basic understanding guides most engineering practices, yet the classical approach still leaves open questions: What is the role of fluid compressibility in the pressurization process [START_REF] Garga | Volume changes in undrained triaxial tests on sands[END_REF]?, What are the relevant drainage conditions -can liquefaction occur also when fluid inflow and outflow to and from the system are allowed [Seed et al., 1976, El Shamy and[START_REF] Shamy | A micro-mechanical investigation of the dynamic response and liquefaction of saturated granular soils[END_REF]? And what is the role of the initial packing -can densely packed layers still liquefy [Soga, 1998, Gabet and[START_REF] Gabet | The mobilization of debris flows from shallow landslides[END_REF]?

The importance of these questions is demonstrated when comparing two numerical models of grain-fluid systems that study soil liquefaction. [START_REF] Okada | Coupling pore-water pressure with distinct element method and steady state strengths in numerical triaxial tests under undrained conditions[END_REF] model an undrained system (with impermeable boundaries) with a highly compressible pore fluid, forced by a compressive constant strain rate. The results of [START_REF] Okada | Coupling pore-water pressure with distinct element method and steady state strengths in numerical triaxial tests under undrained conditions[END_REF] may be interpreted as an example of the classical mechanism of soil liquefaction, as they observe pore fluid pressurization when compacting an initially loosely packed layer subjected to undrained conditions, a situation that was observed also experimentally to lead to pore pressure rise and liquefaction [e.g. [START_REF] Seed | Liquefaction of saturated sand during cyclic loading[END_REF][START_REF] Peacock | Sand liquefaction under cyclic loading simple shear conditions[END_REF][START_REF] Seed | Soil liquefaction and cyclic mobility evaluation for level ground during earthquakes[END_REF]. Unlike this classical approach, El [START_REF] Shamy | A micro-mechanical investigation of the dynamic response and liquefaction of saturated granular soils[END_REF] model a drained system (where the fluid is allowed to flow freely across the top boundary), with completely incompressible pore fluid, an assumption that follows many engineering interpretations of experiments [e.g. [START_REF] Garga | Volume changes in undrained triaxial tests on sands[END_REF]Zhang, 1997, Kozlov et al., 1998]. El [START_REF] Shamy | A micro-mechanical investigation of the dynamic response and liquefaction of saturated granular soils[END_REF] report on significant pore fluid pressurization and liquefaction when forcing their system with periodic shear acceleration at its base. These results are somewhat unexpected because the drained boundary conditions contradict the classical view of liquefaction, which requires poor drainage.

Fault zones and landslides. In fault zones, deformation often localizes along a fault plane filled with fault-gouge. Fault-gouge is a granular layer formed from fragments that are the product of wear during shear between the fault-walls. In gouge layers and in similarlyformed landslide's shear zones, pressurization and depressurization may occur as a result of irreversible rearrangement of the granular skeleton during continuous shearing. In these cases, drainage conditions and porosity evolution were shown to affect the evolution of pore pressure and the strength of the shear zone. Porosity within gouge and shear zones is a function of shearing velocity [START_REF] Marone | Frictional behavior and constitutive modeling of simulated fault gouge[END_REF] and stress conditions [Aharonov andSparks, 1999, Iverson et al., 2000]. When the fault is sealed, porosity increase (i.e. dilation) during shear is often suggested to prohibit unstable sliding via pore pressure reduction leading to strain hardening following equation (2) [START_REF] Scholz | Earthquake prediction -physical basis[END_REF][START_REF] Scholz | Velocity anomalies in dilatant rock[END_REF][START_REF] Rudnicki | Stabilization of rapid frictional slip on a weakening fault by dilatant hardening[END_REF][START_REF] Segall | Dilatancy, compaction, and slip instability of fluid-infiltrated fault[END_REF][START_REF] Moore | Slow episodic shear of granular materials regulated by dilatant strengthening[END_REF][START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF][START_REF] Samuelson | Shear induced dilatancy of fluid saturated faults: experiment and theory[END_REF], while compaction of under-compacted gouge was shown experimentally to lead to extreme weakening and unstable sliding [START_REF] Blanpied | An earthquake mechanism based on rapid sealing of faults[END_REF].

However, similarly to soil liquefaction, the mechanisms that are responsible for pore pressure evolution in shear zones and the effect of pore pressure variation on the mechanics of grains-fluid shear zones are not completely clear. A basic question that is still debated is the possibility of significant pressurization, and as a consequence reduction of shear strength, in an initially densely packed shear zones. This is an important question since despite our knowledge that natural shear zones are in most cases initially over-consolidated and thus, according to conventional thought, resistant to liquefaction [e.g. [START_REF] Iverson | Acute sensitivity of landslide rates to initial soil porosity[END_REF], ample observations point to significant pore fluid pressurization during earthquake and landslides: hydrofractures and liquefied injection are reported in sesmically active fault zones [Boullier et al., 2009, Sagy and[START_REF] Sagy | Geometric and rheological asperities in an exposed fault zone[END_REF] and along landslide shear zones [START_REF] Anders | Stratified granular media beneath large slide blocks: Implications for mode of emplacement[END_REF], and transient liquefaction was observed in experimental landslides that were constructed with an initially densely packed grain and rod layers [START_REF] Iverson | Dynamic pore-pressure fluctuations in rapidly shearing granular-material[END_REF]].

Overview of existing research approaches

To study the mechanics of pore pressure evolution in deforming granular skeleton, and its implications to shear strength, there is a need for a fully coupled theory for the mechanics of fluid-filled granular systems. Such full coupling should include two-directional mechanics:

(1) the effect of granular matrix deformation on the pore fluid pressure and flow, and (2) the effect of fluid flow and pressure gradients on the deformation of the granular matrix. A continuum theory for the first direction (the solid effect on the pore fluid) is available and reviewed below. However, the second direction requires a continuum description for the general dynamics of a collection of grains. Despite a recent advance in this field in the form of constitutive relations for the flow of dry granular material [START_REF] Pouliquen | Flow of dense granular material: towards simple constitutive laws[END_REF], Jop et al., 2006], it is not clear if these relations may apply also to the case of frictional dominated dense suspension of grains immersed in fluid [Rondon et al.], which is the system that we study here. Therefore a well-founded continuum theory for the second direction of our system (the fluid effect on the grains and the resultant grain dynamics) is probably still missing. In the absence of such a theory, alternative approaches are used. One approach that is adopted in many engineering applications for the study of the coupled mechanics of grains and pore fluid is the use of phenomenological models that are based on continuum mixture theory formulations [START_REF] Zienkiewicz | Computational Geomechanics with Special Reference to Earthquake Engineering[END_REF], and include many parameters that require calibration.

A second approach combines the continuum theory for the fluid with a discrete elements numerical method for the solid grains dynamics. Indeed, a common physical method, developed over the last 30 years for dry granular systems, is to use simple interaction laws between individual grains, with few parameters to characterize their interactions [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. The granular rheology then arises from their collective behavior. The combination of the continuum and discrete constitutes allows solving the first direction of the solid effects on the fluid with the continuum component, and the second direction of the fluid effect on the grain dynamics with the discrete component. Such an approach was used in the modeling of instabilities in the flow of granular media and gas fluids [START_REF] Mcnamara | Grains and gas flow: Molecular dynamics with hydrodynamic interactions[END_REF], Vinningland et al., 2007a,b, Johnsen et al., 2006, 2007[START_REF] Johnsen | Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or a gas[END_REF][START_REF] Vinningland | Influence of particle size in rayleigh taylor granular flow instability[END_REF] and liquids of various compressibilities and viscosities [START_REF] Niebling | Sedimentation instabilities: impact of the fluid compressibility and viscosity[END_REF], for modeling hydrofractures [START_REF] Flekkøy | Modeling hydrofracture[END_REF], and for the study of soil liquefaction [START_REF] Shamy | A micro-mechanical investigation of the dynamic response and liquefaction of saturated granular soils[END_REF]Zeghal, 2007, Okada and[START_REF] Okada | Coupling pore-water pressure with distinct element method and steady state strengths in numerical triaxial tests under undrained conditions[END_REF]. Such an approach is adopted in the current work.

The continuum component is developed in Goren et al. [2010] for compressible pore fluid pressurization and flow in response to general (reversible and irreversible) granular matrix deformation (the first direction of the full coupling). This paper couples this continuum component for the fluid with a discrete element granular dynamics algorithm following the scheme presented in [START_REF] Mcnamara | Grains and gas flow: Molecular dynamics with hydrodynamic interactions[END_REF], to form a fully coupled model that may be used in the study of any general granular matrix deformation, and any form of drainage boundary conditions. The analysis of the continuum component by itself already supplies interesting results. The equation describing the evolution of pore pressure in response to grain matrix deformation is shown here to lead to a viscoelastic type of behavior, where the pore pressure (stress) depends on both the pore volume change (overall volumetric strain) and the rate of porosity change (volumetric strain rate). Two end-member behaviors for the evolution of pore pressure emerge from the pore pressure equation: Viscous-like response, in which the pore pressure is linearly dependent on the volumetric strain rate, dominates when the shear zone is well drained. Elastic-like response, in which the change of pore pressure is proportional to the volumetric strain and is inversely proportional to the fluid compressibility, dominates when the shear zone is effectively undrained. This basic recognition of viscoelastic behavior of the pore pressure was predicted in Goren et al. [2010] assuming the pore fluid responds to grain rearrangement, but grains deformation is not affected by the fluid (infinitely stiff approximation). An important goal of this current work is to check whether the two end-member behaviors also characterize the fully coupled two-way model.

Overview of current research

The first goal of the current paper is to develop the full coupling between the pore pressure model of Goren et al. [2010] and a granular dynamics model, and to validate it. Such a coupled model is presented in section 2, and its validation is presented in section 3. Since the pore pressure evolution and its two end-member behaviors were predicted to control the coupled grain-fluid response, section 4 is dedicated to reviewing in detail the pore pressure evolution that emerges from the formulation of the pore fluid, and the conditions that control the pore pressure response to granular deformation.

The rest of the paper is dedicated to simulations with the fully coupled model, their analysis, and comparison with theoretical prediction: The simulations (presented in section 5) investigate shear of a closely-packed fluid-filled granular layers, under constant shear velocity and normal stress. During the simulations we monitor the evolution of pore pressure, porosity, and apparent friction, µ a , which is the friction the layer appears to have (the shear stress that is required for shearing the layer at a constant velocity divided by the applied normal stress). The importance of the apparent friction is a result of our inability to define a single value of pore pressure that may be assigned in equation ( 2). Thus, equation ( 2) is often substituted by τ = µ a σ n .

(3)

We have performed two types of simulations using the fully coupled grains and pore fluid model that differ in their boundary conditions: undrained and drained. Simulation results indeed agree with the two end-member behaviors predicted by the simple pore pressure model in section 4 and in Goren et al. [2010]: when the boundaries are undrained, the pore pressure response is elastic-like. Dilation with respect to the initial configuration causes pore pressure reduction and system hardening (increase of µ a ). Conversely, when shear leads to compaction, the elastic-like behavior causes rapid liquefaction with pore pressure that becomes equal to the applied normal stress throughout the system, and to a steady-state loss of shear resistance. This behavior corresponds to the classical mechanism of liquefaction.

Less intuitive results arise when simulations are conducted with well-drained boundary conditions. Here, viscous-like evolution of pore pressure may lead to high pore pressure values even if the granular system is undergoing long term dilation. The pressurization occurs during short compactive phases that always punctuate the dilative trend. In some cases, instantaneous events of localized high pore pressure may overcome the applied normal stress and lead to transient liquefaction with a complete loss of shear resistance, (µ a ≤ 0). Thus, we predict here that liquefaction in an initially densely packed system is possible, provided that drainage is good. This surprising result agrees with the simplified model of Goren et al. [2010], and suggest a new look at liquefaction potential of natural systems. Section 6 discusses the implications of our model and simulations to natural systems of grains and fluids. We show that the two end-member behaviors we observe are seen both in other simulations and were seen previously in experiments and natural conditions, but because they were not understood well, they were not dwelled upon. In this section we also introduce a new measurable parameter, the liquefaction potential, that is derived from our model and may be used for evaluating the potential for liquefaction of various natural systems with various scales and boundary conditions.

Coupled grains and fluid model

The numerical model for the coupled dynamics of grains and pore fluid is constructed as a two-phase two-scale model. The granular phase makes up the finer scale where each grain is modeled as a discrete lagrangian element. The fluid phase composes the coarser scale and is modeled on an eulerian grid that is super-imposed on the granular level. The fluid grid spacing is approximately the size of two grain diameters. This choice of grid spacing ensures that each cell is larger than a Darcy representative volume element. Similar models reported in [START_REF] Johnsen | Pattern formation during air injection into granular materials confined in a circular Hele-Shaw cell[END_REF] and Vinningland et al. [2007a] have shown that simulation results are mostly insensitive to the exact choice of grid spacing when it is between 2-10 grain diameters. The fluid does not see the detailed pore space geometry imposed by the grains, but instead an average field of porosity and permeability as explained below. Furthermore, although the current model is 2D, the porosity is assumed to be always connected via the third dimension to allow simultaneous percolation of both the grains and the fluid.

Granular phase

To simulate the mechanics of a collection of grains we use a two dimensional discrete elements granular dynamics algorithm [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. Each individual grain is treated as an inelastic soft disc. Grain interactions, body forces and the force induced by the interstitial fluid lead to linear and rotational acceleration of the grains. Interaction force between two grains i and j is resolved when the distance between the centers of the two grains, r ij , is less than the sum of their radii, R i + R j . The grain overlap is expressed as

ξ ij = R i + R j -r ij .
Interaction force on a contact has a normal component, F n , and a shear component, F s , that are resolved with Hertz-Mindlin contact model (figure 1):

F n ij (t) = kn ξ ij -γm ij (ṙ ij • nij ) nij , (4)
where the first term on the right hand side of equation ( 4) is a nonlinear repulsive force and the second term is a damping force that depends on the damping coefficient, γ, the harmonic mean of the grains mass, m ij , and the relative velocity between the grains along the direction of the contact ṙij • nij , where r ij is a vector connecting the grains centers and ṙij is the relative grain velocity. nij = (r ij • x, r ij • ŷ)/r ij is a unit vector normal to the contact. The coefficient of the normal repulsive force is the nonlinear normal stiffness [START_REF] Schäfer | Force Schemes In Simulations If Granular Materials[END_REF]:

kn = √ 2E 3(1 -ν 2 ) (R ij ξ ij ) 1/2 (5)
where E and ν are the grains bulk modulus and Poisson's ratio, respectively, and R ij is the harmonic mean of the grains radii. The shear force is determine using an elastic/friction law:

F s ij (t) = -min( ks ∆s, µF n ij ) ŝij , ( 6
)
where ∆s is the shear displacement since the formation of the contact, µ is the surface friction coefficient, and ŝij = (r ij • ŷ, -r ij • x)/r ij is a unit vector tangent to the contact. The coefficient of the tangent repulsive force is the nonlinear tangent stiffness:

ks = 2 √ 2E (2 -ν)(1 + ν) (R ij ξ) 1/2 . (7)
Equation ( 6) implies that the shear force opposes the sense of the relative tangential motion between the grains. It is important to note that as long as F s ij < µF n ij no sliding occurs along the contact, but there might still be relative tangential displacement between the centers of grains i and j. When Coulomb failure criteria along the contact is met, i.e. ks ∆s ≥ µF n ij , the contact slides with a constant shear force, µF n ij , (figure 1). The motion of each particle is determined by monitoring collisions between grains, resolving the induced contact forces and torques and using them in the momentum equations:

m i ui = m i g + j F ij - ∇P • V i 1 -Φ , (8)
I i ẇi = j R i nij × F ij , (9)
where u i and w i are the translational and rotational velocity vectors of grain i (a superposed dot indicates time derivative). m i is the grain mass, g is the gravitational acceleration, I i is the grain moment of inertia, and F ij refers to inter-grain force at the contact between grain i and grain j. The arm of the force in the torque balance, equation (9), is expressed as R i nij because all grains are perfect discs. The last term on the right-hand side of equation ( 8) refers to the force exerted on grain i by the pressure gradient, ∇P , of the fluid surrounding it, normalized by the solid fraction, (1 -Φ), in its vicinity, where Φ is the porosity and V i is the volume of grain i [START_REF] Mcnamara | Grains and gas flow: Molecular dynamics with hydrodynamic interactions[END_REF].

Fluid phase

The formulation for the physics of the pore fluid is developed in Goren et al. [2010]. For clarity, we briefly review it here. First, mass conservation equations are written for the grains and for the fluid:

∂[(1 -Φ)ρ s ] ∂t + ∇ • [(1 -Φ)ρ s u s ] = 0, (10)
∂[Φρ f ] ∂t + ∇ • [Φρ f u f ] = 0, ( 11
)
where t is the time, ρ s and ρ f are the densities of the solid grains and fluid, respectively, and u s and u f are the solid and fluid velocity fields, respectively. These velocities are defined for mesoscopic volumes containing at least a few grains, where Darcy's law is applicable. In that sense u s is an average of u i over spatially close grains.

The full fluid momentum equation includes inertial terms and forcing terms of pressure gradient and viscous drag. Here we choose to neglect fluid inertia to allow a theoretical analysis of the pore fluid equation. Goren et al. [2010] show that for the deformation field and parameter range that are used here, fluid inertia is mostly negligible, but it is important to note that such an approximation is more suitable when the permeability and grain size are relatively small. The full fluid momentum equation is treated for example in El [START_REF] Shamy | A micro-mechanical investigation of the dynamic response and liquefaction of saturated granular soils[END_REF], and although it is hard to compare the two models, the overall observed behavior when drainage is good (the scenario treated in El [START_REF] Shamy | A micro-mechanical investigation of the dynamic response and liquefaction of saturated granular soils[END_REF]) is similar. When neglecting inertial terms the fluid momentum equation is approximated by Darcy's law:

Φ(u f -u s ) = - k i η ∇P, (12)
where k i is the permeability, η is the fluid viscosity and P is the excess (over hydrostatic) fluid pressure. In the small system we will consider, we take hydrostatic pressure to be uniform throughout the system. The fluid density is given by the fluid state equation:

ρ f = ρ 0 (1 + βP ), ( 13
)
where ρ 0 is the fluid density at hydrostatic pressure level, and β is the adiabatic fluid com-pressibility. We assume that the compressibility of a grain is negligible relative to the fluid compressibility, so that ρ s can be approximated as constant. We also assume that βP ≪ 1 [Goren et al., 2010]. Equations ( 10) to (13) then lead to:

βΦ ∂P ∂t -∇ • [ k η ∇P ] + ∇ • u s + βΦu s • ∇P = 0. ( 14
)
By further assuming that the length scale of pore pressure diffusion is always larger than the diameter of a single grain [Goren et al., 2010], the last term of equation ( 14) may be neglected.

This assumption is revisited in section 4 where a non-dimensional analysis of equation ( 14) is presented. As a result, the pore pressure evolution equation becomes a three terms equation:

∂P ∂t - 1 βΦη ∇ • [k∇P ] + 1 βΦ ∇ • u s = 0. ( 15
)
The first term of equation ( 15) expresses the temporal evolution of pore pressure, the second term expresses pore pressure diffusion, and the third term is the forcing to the pore pressure, which arises due to divergence in the grains velocity. When this divergence is negative, the pore volume collapses, and the fluid pressurizes and flows away from the collapsing pores.

When the divergence is positive, the pore volume expands, and the pore fluid depressurizes and flows into the expanding pore volume.

It is sometimes convenient to express the forcing term as a function of the porosity evolution rather than the divergence of the solid grains velocity. From the grains mass conservation, equation (10), it can be shown that

(1 -Φ)∇ • u s - ∂Φ ∂t -u s • ∇Φ = 0. ( 16
)
The grains velocity divergence term scales with local compaction and dilation (local grains rearrangement), while the porosity gradient term that describes the advection of porosity scales with the imposed shear velocity over the whole layer thickness. Because the length scale associated with the former term is expected to be much smaller than the whole layer thickness (which is associated with the later term), the advection of potosity is neglected, so that

(1 -Φ)∇ • u s ≈ ∂Φ/∂t. (17)
(A similar conclusion is adopted in [START_REF] Walder | Porosity reduction and crustal pore pressure development[END_REF] and [START_REF] Snieder | The liquefaction cycle and the role of drainage in liquefaction[END_REF]). Under this approximation, equation ( 15) may be rewritten as

∂P ∂t = 1 βΦη ∇ • [k∇P ] - 1 βΦ(1 -Φ) ∂Φ ∂t . (18)
Equations ( 15) and (18) were shown in Goren et al. [2010] to be a general form of previous formulations in [START_REF] Biot | General theory for three-dimensional consolidation[END_REF], [START_REF] Wang | Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology[END_REF] and [START_REF] Bachrach | Liquefaction and dynamic poroelasticity in soft sediments[END_REF] that assume only elastic deformation of the grains skeleton, and [START_REF] Walder | Porosity reduction and crustal pore pressure development[END_REF] that assume a specific law for the evolution of porosity. Similar formulations also appear in [START_REF] Iverson | Differential-equations governing slip-induced pore-pressure fluctuations in a water-saturated granular medium[END_REF], [START_REF] Rudnicki | Stabilization of rapid frictional slip on a weakening fault by dilatant hardening[END_REF], [START_REF] Miller | Permeability as a toggle switch in fluid-controlled crustal processes[END_REF], [START_REF] Snieder | The liquefaction cycle and the role of drainage in liquefaction[END_REF] and [START_REF] Samuelson | Shear induced dilatancy of fluid saturated faults: experiment and theory[END_REF].

Coupling between the grains and the fluid phases

In the current model, we couple a 2D granular dynamics algorithm with a continuous solver of the pore fluid pressure and flow. To allow the coupling, information must be transferred between the two phases of the model. Fluid pressure gradients are needed in order to solve the grain force balance, equation ( 8), and the divergence of grains velocity and the porosity are required for the solution of the fluid pressure, equation (15). To achieve this full coupling we use a 2D linear interpolation scheme between the two scales of the model. For a grain whose center is located in coordinate x i and for a fluid grid point located in coordinate x the interpolation scheme is represented by the function s:

s(x i -x) = 1 -|x i -x| lx 1 -|y i -y| ly , |x i -x| < l x , |y i -y| < l y 0 otherwise (19)
where l x and l y are the horizontal and vertical grid spacing. Each grain contributes its sweighted area and momentum to the grid points surrounding it. The porosity along a grid point is calculated as

Φ(x) = 1 - A s (x) A , ( 20
)
where A is the area of a grid cell, and

A s (x) = N i=1 s(x i -x)A i , (21)
A i being the area of grain i, and N is the number of grains. The solid velocity field is defined as the ratio of granular momentum and granular mass. For equal density grains the mass dependency is reduced to an area dependency, and the solid velocity field may be evaluated as:

u s (x) = N i=1 s(x i -x)A i u i N i=1 s(x i -x)A i . ( 22
)
In the simulations presented here, the grain size distribution is close to being monodispersed and equation ( 22) is reduced to:

u s (x) = N i=1 s(x i -x)u i N i=1 s(x i -x) , (23)
where N i=1 s(x i -x) is the on-site mass density [START_REF] Mcnamara | Grains and gas flow: Molecular dynamics with hydrodynamic interactions[END_REF]. Interpolated granular velocities are calculated on a staggered grid with respect to the porosity, so that the velocity divergence is defined exactly on the porosity grid.

The ratio of pressure gradient to the solid fraction, ∇P/(1 -Φ), that is calculated on the fluid grid by solving equation ( 15), is interpolated back from the fluid grid to the grains surrounding this grid using the same interpolation function s, equation ( 19).

The permeability is calculated with a Carman-Kozeny-like relation. However, Carman-Kozeny relation gives the permeability as a function of the volume fraction of spheres, while the porosity in our model is computed with area fraction of discs. Consequently, we transform the area fraction in the simulations (1 -Φ) (2D) to an equivalent volume fraction in 3D. The simplest map of 2D to 3D solid fraction, which ensures that the pure fluid state and the random close packing state correspond between the two dimensionalities, is (1 et al., 2000]. This mapping results in the following relation between porosity and permeability:

-Φ) (3D) = (2/3)(1 -Φ) (2D) [McNamara
k i = k c (1 + 2Φ) 3 (1 -Φ) 2 . ( 24
)
where k c is a prefactor with unit of m 2 , and Φ is the 2D porosity.

Stability and accuracy requirements force us to take a time-step small enough to resolve the evolution of forces during collision of grains. In each time step, equations ( 8) and ( 9) are solved to find the new location, velocity and acceleration of each grain. Then, the granular velocity and the porosity are interpolated from the granular level to the fluid grid. In the next stage, an alternating-direction-implicit (ADI) algorithm is used in the solution of the fluid pressure, equation ( 15), and the pressure gradients are interpolated back to the granular level and assigned in the last term of equation ( 8) at the next time step.

Model validation

In order to validate the coupled grains-fluid model, we perform three tests in which we compare simulation results to analytical predictions.

Effective stress validation test. The first test verifies that the model reproduces correctly the law of effective stress. [START_REF] Nur | An exact effective stress law for elastic deformation of rock with fluids[END_REF] develop an effective stress law for volumetric strain of fluid-filled porous material: σ ′ ij = σ ij -αδ ij P , where α, the effective stress coefficient, is a function of the compressibility of the solid grains and of the matrix. For the case of incompressible grains α = 1. [START_REF] Nur | An exact effective stress law for elastic deformation of rock with fluids[END_REF] assert their law by a series of experiments showing that there is no correlation between the applied stress and the resulting measured volumetric strain, while they find a linear relation between the effective stress and the measured volumetric strain, with a stress-strain curve similar to dry samples [START_REF] Nur | An exact effective stress law for elastic deformation of rock with fluids[END_REF], their figure 2]. In our test we reproduce numerically the experimental series of [START_REF] Nur | An exact effective stress law for elastic deformation of rock with fluids[END_REF]. We perform a series of numerical simulations where in each simulation a system of variable size grains is packed under confining stress. The systems are periodic in direction x so that a normal stress, σ n , applied to the top and bottom walls corresponds to a uniform confining stress. The upper and lower boundaries are composed of half grains that are glued together to form rough walls. In each simulation, a fluid pressure, P , is introduced and maintained constant during the simulation as if the pore fluid in the granular system is connected to a big reservoir. The top and bottom boundaries are undrained and fluid may not flow out or into the system across these boundaries. Simultaneously with the introduction of pressure, the applied unidirectional stress is increased by ∆σ n . Thus, each experiment is characterized by a couple, ∆σ n and P . We measure the volumetric (vertical due to periodicity) strain, ∆ǫ, that results from the extra loading, ∆σ n , under constant pore pressure, P . The setup of the numerical simulations is depicted in figure 2.

Figures 3a and3b show ∆ǫ, the volumetric strain, as a function of ∆σ n and ∆σ n -αP , respectively, similarly to Nur and Byerlee [1971, figure 2]. A set of dry simulations, with no pore fluid, serves as a reference and is depicted by 'x'. Wet simulations are depicted by 'o'. Figure 3a shows that ∆ǫ and ∆σ n in the wet simulations are poorly correlated. In contrast, figure 3b shows a linear relation between ∆ǫ and ∆σ n -αP , for α = 1 (as expected for incompressible solid grains), with the same slope as the dry simulations. Simulation results show that the model successfully reproduces the effective stress behavior that is observed experimentally in fluid-filled porous and granular material.

Since pore pressure is maintained constant within the system, ∇P = 0. Therefore, the fluid does not exert a force on internal grains, and the last term of equation ( 8) vanishes. However, since the system is undrained and the pore pressure outside the system is assumed zero, boundary grains do feel the effect of fluid pressure because it is exerted on them only from one side (from the bottom on the top boundary and from the top on the bottom boundary). Thus, a pressure gradient force acts on the boundary grains in the outward normal direction to the boundaries, and opposes the external force induced by σ n . We note here an important insight regarding the effective stress: in models, one could be tempted to apply the effective In each experiment, a unidirectional stress, ∆σ n , is applied, and a constant pore pressure, P , is maintained. The vertical strain, ∆ǫ, is measured. The white pore area in between the grains is filled with fluid. stress law at each grain contact, but that would lead to an incorrect formulation of the forces, as only pressure gradients exert a net force on grains. The effective stress law may therefore be viewed as the macroscopic manifestation of microscopic gradients of pore pressure.

Sedimentation validation test. The second test compares the sedimentation velocity under gravity of grains suspended in fluid, with the theoretical prediction that assumes that particles fall without acceleration (following McNamara et al. [2000]). The volume of grains transported downward must be compensated by an equal volume of fluid upflow. Thus (1 -Φ)u s = -Φu f . Combining this equation with Darcy law, equation ( 12), results in u s = k(Φ)∇P/η. The pressure force must balance the weight of the grains, so that ∇P = ρ s g(1-Φ). Assigning equation ( 24) for k(Φ), with k c = r 2 /540 [START_REF] Mcnamara | Grains and gas flow: Molecular dynamics with hydrodynamic interactions[END_REF], r being the average grain radius, leads to:

u s u sed = - [1 -2 3 (1 -Φ)] 3 (1 -Φ) ẑ, (25)
where u sed = r 2 ρ s g/20η, and ẑ is a unit vector in the opposite direction of gravity. Each sedimentation simulation starts with a different uniform porosity. Then, every 50,000 time steps we average the porosity and the granular velocity over the grid points. In order to x 10

-3 α = 1 (a) (b) ∆σ n [MPa] ∆σ n -α P [MPa]
∆ε ∆ε exclude material that was already settled on the bottom, and the clear region above the settling grains, we do not include in the average grid points with granular velocity smaller than 0.5 max(u s ), and grid points with (1 -Φ) < 0.25. Finally, each simulation is temporally averaged to obtain a space and time average of u s and (1 -Φ) for the simulation. These couples are depicted by 'o' in figure 4 and show good fit to the theoretical prediction of equation ( 25).

Fluidization validation test. The third test is of fluidization of a granular layer. Initially the grains rest at the bottom of the system after sedimentation. Then, a constant fluid pressure gradient is applied between the top and the bottom of the system. It is predicted that when ∇P = -(ρ s -ρ f )(1 -Φ)g, the force exerted by ∇P will exactly balance the weight of the grains that rest at the base of the system [START_REF] Richardson | Incipient fluidization and particulate system[END_REF]. For larger ∇P the grains will be lifted, while for smaller ∇P the grains will not move. Figure 5 shows the granular velocity averaged in space over the grains and in time over the first million time steps, ūs , plotted vs. ∇P for several simulations. The grains move only when ∇P is larger than the predicted critical value, while for smaller values of ∇P , ūs = 0.

4 Pore pressure evolution mechanisms

After presenting the fully coupled model and verifying it, we turn to discuss the two mechanisms that are responsible for the evolution of pore pressure, and that arise from the pore fluid formulation presented in section 2.2. These mechanisms, which are generic and independent of the specific deformation of grain dynamics, are discussed at length in Goren et al. [2010]. The two mechanisms depend on the system boundary conditions, and on the relative magnitude of the different terms in equations ( 15) and ( 18). In order to estimate the magnitude of these terms we perform a non-dimensional analysis of equation ( 14). The characteristic scales of the variables in the system are defined as: P = P /β, u s = ûs u 0 , k = kk 0 , and t = tt 0 , where the ˆnotation denotes non-dimensional variables, and u 0 , k 0 , and t 0 are the velocity, permeability, and time scale factors, respectively. The divergence operators resulting from grain mass conservation, equation ( 10), represent grain-scale rearrangements, and therefore are scaled by d -1 , where d is a characteristic grain diameter. The gradient operators that result from Darcy law, equation ( 12), and act on the pressure, P , are scaled by l -1 , where l = min(ζ, D i /u 0 ) is a larger length scale that corresponds to the distance that is reached by the pore pressure signal. ζ is the system half thickness, which is also the distance between the system's center and its boundaries. D i = k i /βηΦ is the internal diffusion coefficient of the system that depends on the internal permeability, k i , with accordance to equation (24) (while ignoring the boundary permeability), and D i /u 0 is the internal diffusion length. D i /u 0 may be cast as a (more common) diffusion length scale of the form

0 10 -2 10 0 10 -1 u s (m/s) ∇P [-(ρ s -ρ f )(1-φ)g]
√ D i t 1 , with t 1 = D i /u 2 0 .
To understand the meaning of t 1 note that when grains move at velocity u 0 , the pressure next to where the motion takes place changes first due to diffusion and later due to advection. t 1 is the time at which the diffusive and advective influences balance. According to the definition of l, if the system is relatively small or the internal permeability is relatively large, then ζ < D i /u 0 and a pore pressure signal will interact with the boundaries since the system is well connected with diffusion, leading to l = ζ. If ζ > D i /u 0 , then diffusion does not have the necessary time to level out the pressure changes advected with the matrix at velocity u 0 during shear, and l = D i /u 0 . The time scale factor, t 0 , is defined as the time scale of deformation:

t 0 = d u 0 . ( 26
)
The permeability scale factor, k 0 , requires further discussion. When, ζ > D i /u 0 and l = D i /u 0 , then the boundaries are not expected to interact with a pore pressure signal originating at the system interior, and therefore, k 0 = k i . However, when ζ < D i /u 0 and l = ζ, k 0 should represent the effective permeability that accounts both for the value of the internal permeability, k i , over a layer with thickness 2ζ, and for the value of the boundary permeability, k b , over a thin (but finite) boundary layer with thickness δ. As the permeability is proportional to Darcy velocity, k 0 is estimated as harmonic mean:

k 0 = 2ζk i k b /(δk i + 2ζk b ), (27)
where 2ζ + δ ≈ 2ζ. The harmonic mean is for permeabilities transverse to the boundaries and so it gives larger weight to the smaller permeability. As a result (28)

k 0 =          k i for
Assigning the non-dimensional variables in equation ( 14) results in:

∂ P ∂ t -De -1 ∇ • ( k ∇ P ) + 1 Φ ∇ • ûs + d l ûs • ∇ P = 0. ( 29
)
The last term of equation ( 29) is negligible with respect to the first and third terms because for any natural system, and in particular for our choice parameters, the diffusion length is significantly larger than the diameter of a single grain, i.e. d ≪ l [Goren et al., 2010]. Equation (29) then becomes:

∂ P ∂ t -De -1 ∇ • ( k ∇ P ) + 1 Φ ∇ • ûs = 0, (30)
which is the non-dimensional form of equation ( 15). The coefficient of the second term in equation ( 30) is a function of the Deborah number,

De = t d t 0 , (31)
that is defined as the ratio of relaxation time scale and a characteristic process time scale [Osswald, 1998, p. 54]. Here, the relaxation time scale, t d = ld/D, is the time scale for pore pressure diffusion, where D = k 0 /βηΦ is the system diffusion coefficient (note that D ≤ D i ).

The characteristic process time scale, t 0 , is simply the time scale of deformation (equation ( 26)). In terms of the system parameters, De may be expressed as:

De = lu 0 D = lu 0 βηΦ k 0 . ( 32
)
Deborah numbers are normally used in the characterization of viscoelastic materials for estimating the relative importance of the viscous and elastic rheologies. In the following we show that viscoelasticity is a good analog for the description of the pore fluid pressure evolution.

4.1 Small system (ζ < D i /u 0 )

When the system is relatively small compared to the diffusion length, then we take l = ζ in the definition of De,equation (32). If the boundaries are undrained, k 0 is zero and De -1 = 0. As a result, the diffusion term in equations ( 15) and ( 18) should be ignored. The non-dimensional equation (30) then becomes

∂ P ∂ t + 1 Φ ∇ • ûs = 0. ( 33
)
Under such conditions, the dimensional equation ( 18) can be solved for the spatially averaged pore pressure P :

∆ P = - ∆Φ βΦ(1 -Φ) , (34)
where Φ(1 -Φ) is approximated as a constant, ∆ P = P (t ′′ ) -P (t ′ ), for any t ′′ > t ′ , and ∆Φ is defined in a similar manner. Here, changes of average pore pressure, ∆ P , are linearly related to the overall change of porosity, ∆Φ, with a proportionality factor that depends on the fluid bulk modulus, β -1 , so that pore pressure responds 'elastically' to pore strain. Pore fluid that is trapped within a shrinking pore volume (∆Φ < 0) is pressurized, while pore fluid trapped in expanding pore volume (∆Φ > 0) is depressurized.

When the boundaries are well drained (k 0 is approximated as k i ), then always (by the definition of De and by the condition ζ < D i /u 0 ) De < 1. If De ≪ 1, for example when ζ is small or k 0 is large, then the non-dimensional equation (30) reveals that the first time dependent term is negligible relative to the second diffusion term. The reason the forcing term is not negligible with respect to the diffusion term is that the forcing is the source of pore pressure variations and therefore it is considered the pivot of the equation [Goren et al., 2010]. The non-dimensional pore fluid equation, (30), then takes the form:

De -1 ∇ • ( k ∇ P ) = 1 Φ ∇ • ûs . ( 35
)
Appendix A presents the solution of the dimensional equation ( 18) for the pore pressure, while neglecting the time dependent term. Appendix A shows that the pore pressure within a system that is characterized by De ≪ 1 may be approximated as:

P (z, t) = - η 2k 0 d Φ(ζ, t) dt ζ 2 -z 2 , ( 36
)
where Φ(ζ, t) is the overall or average porosity in the zone between the center of the system and its boundary, that is located at distance ζ from the center (see figure 7). In the derivation of equation ( 36) it was assumed that the top and bottom boundaries are completely drained, i.e. P (ζ, t) = P (-ζ, t) = 0, and that the rate of change of the average porosity is uniform in space. Note that since (ζ 2 -z 2 ) ≥ 0, compaction with d Φ(ζ, t) /dt < 0 leads to pressurization, while dilation with d Φ(ζ, t) /dt > 0 leads to depressurization. Here, the evolution of pore pressure is controlled by the strain rate of porosity d Φ(ζ, t) /dt, with a coefficient that depends on the fluid viscosity, η, and is inversely proportional to the permeability. Therefore, the pore pressure responds 'viscously' to the deformation of the matrix.

When De 1, which may occur for larger ζ or smaller k 0 , the diffusion term is not strictly larger than the time dependent term, and both terms are expected to contribute to the evolution of pore pressure. Still equation (36) should give a rough estimation of the pore pressure evolution.

Large system (ζ

> D i /u 0 )
When the system is large compared to the diffusion length, then we take l = D i /u 0 and k 0 = k i in equation ( 32), which always results in De = 1 (and t d = t 0 ). In this case, the diffusion and time dependent terms are of the same order. The non-dimensional equation (30) then becomes

∂ P ∂ t -∇ • ( k ∇ P ) + 1 Φ ∇ • ûs = 0, ( 37
)
and the dimensional evolution of the pore pressure is governed by the three terms equations ( 15) and ( 18).

The behavior of the pore pressure evolution

This analysis shows that the pore pressure equations ( 15) and ( 18) express a viscoelastic-like rheology. The two end-members of this rheology, elastic and viscous, lead to two mechanisms that control the evolution of pore pressure, and operate under different drainage conditions. A schematic representation of the two mechanisms is depicted in figure 6. Elastic-like pore pressure evolution dominates when the system is effectively undrained (De ≫ 1). Under such conditions, pore fluid that cannot escape and is trapped within a shrinking pore volume is pressurized (figure 6a), while pore fluid that is trapped in an expanding pore volume is depressurized. The evolution of the average pore pressure will follow equation ( 34), and the magnitude of pressurization and depressurization depends on the inverse of the fluid compressibility, β -1 , and is controlled by the overall change of porosity, ∆Φ. In that sense, the elastic end-member holds memory of the initial state of porosity.

Viscous-like pore pressure evolution dominates when the system is effectively drained (De ≪ 1), and is a less intuitive mechanism. Here, because of mass conservation, convergence (or divergence) of grains causes the pore fluid that resides between the grains to flow out of (or into) this region (figure 6b). Because of fluid momentum conservation (here Darcy law), pressure gradients must arise between the location of converging (or diverging) grains and the surrounding region, to generate these flows. Here, pore pressure evolves from the arising pressure gradients, and is governed by equation ( 36). The evolving pore pressure depends linearly on the fluid viscosity, η, and inversely on the permeability, k 0 . The magnitude of pore pressure is controlled by the instantaneous strain rate of porosity, d Φ(ζ, t) /dt. Therefore, this mechanism holds no memory of previous states of porosity. Note that when De ≈ 1, an intermediate behavior is expected with some short term memory.

The viscous end-member has normally not been offered as a mechanism for liquefaction, although it may lead to significant pressurization. Moreover, because of its 'lack of memory', this mechanism may lead to generation of high pore pressure even when an initially dense granular matrix is sheared. Indeed, upon shearing an over-compacted layer, it will first dilate (figure 6b, left to center), and then oscillate around its critical porosity [Aharonov andSparks, 2002, Gabet and[START_REF] Gabet | The mobilization of debris flows from shallow landslides[END_REF]]. In the oscillatory stage, any compaction phase, with d Φ(ζ, t) /dt < 0, will lead to pressurization despite the fact that the instantaneous porosity may be significantly larger than the initial porosity (figure 6b, center to right).

It is of interest to note that the viscous-like rheology arises when De ≪ 1, and the elasticlike rheology arises when De ≫ 1. Indeed, a Deborah number of zero represents a viscous fluid and an infinite Deborah number represents an elastic solid [START_REF] Osswald | Polymer Processing Fundamental[END_REF]]. In the following, we present our simulation results in light of the pressurization mechanisms that are reviewed here. In order to facilitate this flow, pressure gradient must arise with low pressure within the system interior. Upon compaction (center to right), fluid escapes from the system, and an opposite pore pressure gradient arises with higher pore pressure in the system interior. In this mechanism the evolving pore pressure depends linearly on fluid viscosity, η, and inversely on the permeability, k 0 .

Simulations of shearing granular layers

To study the coupled mechanics of granular matrix deformation and pore fluid pressurization and flow we perform simulations of shearing a fluid-filled granular layer at a constant shear velocity. The simulations are performed in a rectangular system with about 1680 (24 × 70) or 864 (24×36) grains. Grain diameters are drawn randomly from a Gaussian distribution with an average d, and a standard deviation d, clipped at plus/minus 0.2d. Although there is no gravity in the simulations, we define the vertical and horizontal directions for convenience.

The top and bottom walls are composed of half grains of variable size that are glued together along their centerlines to form rigid rough walls. In the horizontal direction the system is periodic, and is thus analogous to a rotary shear apparatus. Each simulation is initiated by compacting a system of loosely packed grains under some normal stress, σ n , until the porosity equilibrated. Then, we assume the pore space is filled with fluid at zero excess fluid pressure. Variations of pore pressure are measured with respect to the initial zero value that corresponds to hydrostatic conditions. For that reason, σ n is interpreted as the applied external stress minus hydrostatic pore pressure, i.e. the initial effective stress. Finally, a constant shear velocity, V sh , is applied on the top wall. During a simulation, σ n and V sh are maintained constant, and we follow the systems' compaction and dilation, the shear stress that is required in order to shear the top wall at a constant velocity, and the evolution of pore pressure. Dividing the shear stress by σ n results in the apparent friction coefficient, µ a , equation (3). As the grains themselves are considered incompressible, compaction and dilation are accommodated by change in pore volume. Figure 7 summarizes the settings of the simulations.

In the simulations, we take the pore fluid to be water, and as a consequence fluid compressibility and viscosity are β = 4.5 × 10 -10 Pa -1 and η = 10 -3 Pa s, respectively. The granular phase is assumed quartz grains with a bulk modulus, E = 8 × 10 10 Pa, Poisson's ratio, ν = 0.25, density of 2640 kg• m -3 , surface friction coefficient µ = 0.5, and damping coefficient (in equation ( 4)), γ = 0.8. We further assume that the grains have an average diameter d = 10 -3 m. As we are interested in studying the role of the permeability on the evolution of pore pressure, and since the permeability varies by orders of magnitude between different natural terrains [e.g. Saar and Manga, 2004], we vary k c in our simulations, where k c is the permeability prefactor that appears in equation ( 24). For the original 2D Carman-Kozeny relation with average grain diameter of d = 10 -3 m, k c = 4.6 × 10 -10 m 2 , but in the present simulations we use a range of k c = 4.6 × 10 -10 -4.6 × 10 -15 m 2 . The lower value of k c may alternatively be obtained by choosing a smaller grain diameter of d ∼ 10 -6 m, or by taking a heterogeneous grain size distribution. We choose to work with relatively large grains in order to allow simulations of thick layers (up to ∼ 0.07 m) during 0.5 second, in a b We report here the maximal vertical half thickness of the system during a simulation. ζ fluctuates by as much as 3% c k 0 is defined in equation ( 28). When k 0 = k i , equation ( 24) is used with Φ = 0.2. When k 0 = k b = 0, the value of k i is reported in parenthesis.

d In the calculation of the Deborah number, Φ = 0.2, and the velocity scale factor, u 0 , is estimated as V sh . In simulations D10, l = D i /u 0 , for the rest of the simulations l = ζ e Liquefaction potential. For simulations D1 -D9 and D11 -D12 equation ( 40) is used. For simulation D10 equation (??) is used. For simulation D14 equation ( 41) is used.

f Simulations U13 and U14 differ in their initial porosity.

reasonable computation time, but we still desire to study the role of the permeability. The applied shear velocity, V sh , is either 0.76 m/s or 7.6 m/s. We use applied normal stresses of σ n = 24 or 2.4 MPa, corresponding to depths of about 2 km and 200 m, respectively. Two types of simulations are performed that differ in their drainage boundary conditions: drained and undrained. Table 1 summarizes the system parameters in the simulations: normal stress, σ n , half thickness, ζ, permeability scale factor, k 0 , and shear velocity, V sh . The dimensionless Deborah number, De, whose significance is explained is sections 4, and the dimensionless liquefaction potential, LP, that is based on section 6.3, are also presented in table 1.

Drained systems

Completely drained systems are simulated by setting the pore pressure to be zero along the top and bottom walls. All drained simulations start with a well compacted system that dilates in the initial stages of shear deformation. After the dilation stage, the porosity oscillates around some mean critical value. Figure 8 compares the temporal evolution of three spatially-averaged quantities (porosity, average pore pressure, and apparent friction, µ a ) between simulations D3 and D9 that differ only in their permeability. Two differences are observed for the average pore pressure in these simulations: First, the peaks of the average pore pressure are smaller by an order of magnitude in the high-permeability simulation D3 (figure 8c) with respect to the low-permeability simulation D9 (figure 8d). This results from the dependency of the pore pressure on the inverse of the permeability for drained systems with De ≪ 1. Indeed, equation ( 36) predicts that when the permeability is smaller, as is the case for simulation D9, the pore pressure will be higher. Second, the average pore pressure is negative during the first half of the low-permeability simulation D9, while for the high-permeability simulation D3, the average pore pressure oscillates around zero from the onset of the simulation. This is due to the competition between the rate of deformation and the rate of pore fluid flow, and will be discussed further in section 6.1. The differences in the evolution of pore pressure between simulations D3 and D9 causes their apparent frictions to differ. The apparent friction starts with a larger value and exhibits larger peaks in the low-permeability simulation D9 (figure 8f). The higher initial value results from the initially negative pore pressure that increases the effective normal stress and the shear resistance following equation (2). Similarly, the larger peaks in the apparent friction result from a larger magnitude of negative values in the pore pressure.

In the drained simulations D1 -D7 (with De ≪ 1) the average pore pressure, P , is well correlated with the temporal derivative of the porosity, dΦ/dt, but is not correlated with the actual value of the porosity, Φ. An example of these relations is depicted in figure 9 for the drained simulations D7. Indeed, equation (36) predicts that when the boundaries are well-drained and De ≪ 1, the evolution of pore pressure is controlled by the temporal derivative of the porosity. Although equation ( 36) was developed based only on the pore pressure formulation, the good correlation that is depicted in figure 9a indicates that viscouslike evolution of pore pressure occurs also when the grain and pore fluid are fully coupled.

In the drained simulations D8, D9, and D11 (with De 1) and in simulation D10 (with De = 1), it is less clear which of the two correlations ( P vs. dΦ/dt and P vs. Φ) is more dominant. In fact, the average pore pressure seems to be influenced both by the temporal derivative of the porosity, dΦ/dt, and by the actual value of the porosity, Φ, as can be seen in figure 10 for simulation D9. Indeed, according to the mechanistic analysis presented in section 4, when De 1 or De = 1 both dependencies are expected. The drained simulation D12 (with De = 6.41 × 10 -2 ) shows a fairly good correlation between P and -dΦ/dt (that is slightly less distinct than this correlation for simulations D1 -D7), and also a minor influence of -Φ on the evolution of P .

Simulations D1 -D10 are conducted under normal stress of σ n = 24 MPa. The average pore pressure in these simulations is significantly lower than the applied normal stress, with F igure 8: Time evolution of space-averaged quantities for drained simulations D3 and D9 that differ in their assigned permeability. Evolution of porosity (a, b), average pore pressure (c, d), and apparent friction (e, f). See table 1 for simulation parameters. peaks that are one or two orders of magnitude smaller than σ n . Still, in runs D8 -D10 that are characterized by small permeabilities and De 1 or De = 1, there were short events where P > σ n . These events, being also very localized in space, did not seem to affect the system's strength. For that reason, we next consider simulations where the applied normal stress is reduced to σ n = 2.4 MPa. It is observed that although σ n does not enter the formulation of the pore fluid, its reduction leads to decrease of pore pressure. However, an order of magnitude reduction in σ n did not cause the pore pressure peaks to decrease by an order of magnitude, but only by a factor of 2-3. Therefore, the average pore pressure peaks may become of the order of the applied stress and lead to a significant reduction of the system strength. Here, the system strength is expressed as the apparent friction coefficient, µ a , where small µ a results from small shear resistance. If the apparent friction becomes zero or even negative, then the shear resistance of the system is completely lost, and it is liquefied.

Figure 11 shows the evolution of µ a and of the average pore pressure for the drained simulation D11. A brief liquefaction event with a reduction of the apparent friction below zero (circled) is observed. This short event is correlated with system average pore pressure of P = 1.9 MPa. Although P < σ n , this liquefaction event is characterized by horizontal layers that experience high pore pressure of P ≃ σ n , and few localized zones with P > 4 MPa (figure 12e).

Figure 12 shows three snapshots of the grain system configuration and the corresponding pore pressure map, before, during and after the liquefaction event circled in figure 11. In frames 12a -12c the thickness of the lines connecting grain centers indicate force on grain contacts. Stress chains are observed as connected force lines that percolate from the top of the system to its base. The color code of the grains corresponds to the overall normal stress that a Figure 11: Apparent friction (blue) and average pore pressure (green) for drained simulation D11 are plotted as functions of time, and show good correlation. Note that the pore pressure axis increases downwards. The apparent friction reduces below zero (circled) at time 0.32 s, when the average pore pressure becomes high, P ∼1.9 MPa, which is of the order of magnitude of the applied normal stress, σ n = 2.4 MPa. Negative apparent friction is defined here as liquefaction. Framed letters mark the time when the snapshots in figure 12 are taken.

grain sustains due to contact forces, with warmer colors for higher stress. In each of the frames 12a -12c, the grain color code is relative to the frame itself, where the grain that supports the maximum load is red and the grains with minimal load are light blue. Therefore, grains that are connected by thick lines and participate in the major stress chains that support the external loading are more reddish. It is clearly observed that the stress chains that support the external load before the liquefaction event, in figure 12a, disappear during the liquefaction event in figure 12b, and the grains at the top and bottom of the system become light blue, indicating that the compressive stress on them arising from granular contacts is minimal.

Percolating stress chains reappear after the event, in figure 12c. The pore pressure before the liquefaction event in figure 12d is mostly negative (lower than hydrostatic). During the event the pore pressure at the top and bottom of the system becomes highly positive (figure 12e), and exceeds σ n . Zones of very high pore pressure in figure 12e correlate to zones with no stress chains in figure 12b, and indicate that the high pore pressure replaces the stress chains in supporting the external load. This high pore pressure also allows for the apparent friction to become negative. Shear localization is observed along the highly pressurized layers, but It is not clear if the localization preceded (and triggered) pressurization or vise versa.

Simultaneously with the reappearance of stress chains in figure 12c, the pore pressure reduces (figure 12f), and the system strengthens.

To study the effect of system size, simulation D12 is conducted with the same parameters as simulation D11, but with a layer half as thick. Figure 13 shows a liquefaction event (circled) during simulation D12 that occurs simultaneously with an increase of average pore pressure, P . The average pore pressure during this event, P ≃ 2.4 MPa, is very close to the applied normal stress. Figure 14 shows snapshots of the granular configuration and the corresponding pore pressure map before, during, and after this liquefaction event. While in simulation D11, high pore pressure is localized close to the top and bottom boundaries (figure 12e), here, when the layer is thinner, stress chains disappear and pore pressure rises throughout the whole system interior, as seen in figures 14b and 14e. It should be noted that several such liquefaction events occurred during simulations D11 and D12.

Undrained systems

Undrained systems are simulated by assigning zero fluid flux across the top and bottom boundaries, and as a result, fluid mass is conserved within the grains layer. In the undrained simulations we assign large internal permeability that allows for rapid fluid flow within the system interior, and as a consequence rapid homogenization of pore pressure. Two simulations are performed, both with σ n = 2.4 MPa and V sh = 0.76 m/s, but with different initial porosities. Simulation U13 starts with the same initial configuration as the drained simulation D11, with an initial 2D porosity of 0.1719. This simulation is referred to as dense. In simulation U14, before the onset of shear, all grains that have no contacting neighbors and thus are not participating in stress chains are removed to form high initial 2D porosity of 0.2385. Simulation U14 is referred to as loose. Dense simulation, U13, dilates when sheared and the pore pressure decreases and becomes negative (i.e., smaller than hydrostatic). Due to the large internal permeability, the pore pressure is uniform throughout the system so that the average pore pressure is very close to the pore pressure at any point within the system. Unlike the drained simulations, here, the evolution of pore pressure is correlated with the evolution of porosity, Φ, (figure 15b), and not with the temporal derivative of the porosity, dΦ/dt, (figure 15a). The dependency between the average pore pressure and the porosity when the boundaries are completely undrained and De = INF is predicted by equations (34). Therefore, although equation ( 34) is developed based only on the fluid formulation, simulation results indicate that an elastic-like evolution of the pore pressure is observed also when the two phases of the system, pore fluid and grains, are fully coupled. The apparent friction coefficient, depicted in figure 15c, increases with decreasing pore pressure. High apparent friction results from high shear resistance due to increasing effective normal stress, following equations ( 1) and ( 2). An alternative view might be that the reduced pore pressure within the system interior exerts a suction force that hardens the system. Simulation U13 thus exhibits 'dilatancy hardening'.

Loose simulation, U14, exhibits a short transient dilation phase with pore pressure decrease upon the onset of shear, but then compaction occurs very fast. Upon compaction, pore pressure increases to the value of the normal stress, P = 2.4 MPa, and the system liquefies (figure 16a). Here, the liquefaction event is a steady state and the pore pressure remains equal to σ n until the end of the simulation. During this long liquefaction event, grains detach, stress chains disappear, and the external load is completely balanced by the pore pressure (figures 17). With the exception of very few contacts arising due to small collisions, the force balance on the grains is zero and they do not touch. Since grains do not slide past each other (and since viscous resistance to sliding is not accounted for in our model), there is no resistance to shear and the apparent friction drops to zero, as depicted in figure 16b.

Discussion

This paper uses a continuum equation (( 15) or ( 18)) to describe the pore pressure response to granular media deformation, coupled with a discrete description of the grain dynamics. First in the discussion we address the pore pressure behavior under any general deformation of the granular skeleton, as it turns out that this behavior by itself is very rich (as already suggested by a simplified model [Goren et al., 2010]). Section 4 reveals that the evolution of the pore pressure in response to granular skeleton deformation obeys a viscoelastic-like rheology. Indeed, equation ( 18) resembles the Maxwell model of viscoelasticity: dǫ/dt = η -1 σ + βdσ/dt, where the strain rate, dǫ/dt, is a linear combination of the stress, σ, and of the time derivative of the stress, dσ/dt. In equation ( 18), the time derivative of the porosity, dΦ/dt, stands for the strain rate, and the pore pressure, P , stands for the stress. The pore pressure diffusion term in equation ( 18) acts as the term η -1 σ in Maxwell model (where the spatial derivative can be approximated as P/ld). Section 4 also shows that the system parameters, as expressed by the Deborah number, De, lead to the emergence of two end-member mechanisms for the evolution of pore pressure that are encapsulated in this rheology -elastic and viscous. Goren et al. [2010] showed that these two end-member mechanisms describe well the evolution of pore pressure when the loading is assumed infinitely stiff, i.e., when granular deformation is externally prescribed and pore pressure gradients do not enter the grains momentum balance. Here, our simulation results show that the two end-member mechanisms also apply for the fully coupled system, i.e., when the skeleton deformation affects the evolution of pore pressure, and pore pressure gradients add up to the force balance on the grains and contribute to the deformation of the skeleton. Therefore, we claim that a viscoelastic analog for describing the evolution of pore pressure applies generally to fluid-filled granular (and porous) systems, and it is independent of the exact details of grain dynamics.

Following these insights, one would like to understand how to estimate De, as it plays a crucial role in determining the pore pressure behavior in real systems. Some parameter groups are of particular importance in the estimation of De. For example, when the system is large with respect to the internal diffusion length, i.e. ζ > D i /u 0 , (where ζ is the system half thickness, D i is the internal diffusion coefficient in the system that ignores the drainage boundary conditions, and u 0 is the rate of deformation), then pore pressure signals originating close to the layer center are unaffected by the boundaries. This situation leads to De = 1, which means that the pore pressure responds both to volumetric strain rate and to overall volumetric strain, resulting in a combined viscoelastic-like behavior. When the layer is relatively thin, as is often the case for experimental systems and field configurations (see for example the cases discussed in section 6.3), then ζ < D i /u 0 , De = 1, and the drainage boundary conditions play an important role in the evolution of pore pressure. These drainage conditions are accounted for in the permeability scale factor, k 0 . Here, we have considered only completely drained and completely undrained systems. However, our choice to express k 0 as the harmonic mean of interior and boundary permeabilities is general: It may be used to evaluate De for relatively complex systems, composed of several layers with different permeability, and to evaluate the expected effect of more complex drainage boundary conditions that assume a linear combination of pressure and pressure gradient leading to an intermediate situation of partial drainage.

Although the pore pressure equations ( 15) and ( 18) are rich and highly predictive of the classes of behavior observed in the simulations, there are some non-linear effects that arise due to the coupled response with the grain dynamics, which cannot be predicted analytically. Such an effect is the relation between the pore pressure and the applied normal stress, σ n . The pore fluid pressure formulation, equations ( 15) and ( 18), predicts that the evolution of pore pressure is independent of σ n . Yet, it is observed that when σ n decreases, the pore pressure is generally smaller. This indirect relation between the applied stress and the pore pressure arises due to coupling with grain contact forces that transmit stresses of the order of σ n . In response to grains convergence, pore fluid pressurization resists the converging grains by exerting pressure gradients across them. Grains then rearrange so that the skeleton forcing on the fluid is relaxed and so is the pore pressure. If σ n is large, large pressure gradients (and large pressures) are needed to push the grains aside and overcome the contact forces. If σ n is smaller, then smaller pressure gradients are sufficient to overcome the granular contact forces, push aside converging grains, and relax the fluid forcing.

The evolution of pore pressure with drained conditions

Next, we turn to analyze the behavior of the pore fluid pressure in simulations with drained boundary conditions. Figure 9a shows the good correlation between the spatially averaged pore pressure and the temporal derivative of the porosity for a representative drained simulation with De ≪ 1, following the prediction of equation ( 36). In order to further validate this correlation, we compare the slopes of the graphs of P vs. -dΦ/dt between the drained simulations D1 -D10 and equation (36). For each time step in a simulation, we plotted the global extremum of the pore pressure and the instantaneous -dΦ/dt. We then take the slope of the linear regression line of P vs. -dΦ/dt and plot it as a function of the permeability, k 0 , in figure 18, as 'o'. In the figure, the slope predicted by equation (36), η(ζ 2 -z 2 )/2k 0 , is depicted by a solid line, where the center of the grains layer is assigned for z, i.e. z = 0. A good fit is found between simulation results and the analytical prediction. It is of interest to note that also simulations D8 and D9 that have De 1, and simulation D10 with De = 1 show good fit to the analytic prediction of equation ( 36) that is developed under the assumption of De ≪ 1.

Although the overall behavior of pore pressure with drained conditions follows equation (36), there are some differences in between the drained simulations that arise from the different parameters that are used. To study these differences we use again the value of the Deborah number. The definition of De in equation ( 32) may also be viewed as a velocity ratio between the velocity of deformation, u 0 , and the velocity of diffusion D/l [START_REF] Samuelson | Shear induced dilatancy of fluid saturated faults: experiment and theory[END_REF]. When De ≪ 1, the deformation is slow enough to allow for a pore pressure front originating at any depth in the layer to reach the drained boundaries in the time scale of deformation. When De approaches 1, the deformation becomes more rapid, or the boundaries are further away from the center of the grains layer, and pore pressure cannot always diffuse across the whole system during the time scale of deformation. In simulation D2, De ≪ 1. Therefore, pore pressure diffusion originating at any depth may reach the drained boundaries at the time scale of deformation. As a result, the pore pressure everywhere in the system is expected to follow equation ( 36), resulting in a parabolic profile, as depicted in figure 19a. In contrast, in simulation D8, De 1. Therefore, a pore pressure front originating at the center of the layer may not always reach the drained boundaries during the deformation time scale. For that reason, non-parabolic profile is observed along the center of the layer (far from the boundaries), as depicted in figure 19b. Deviation from the parabolic profile is probably the result of transient elastic effects arising from the relation between P and Φ when De ≈ 1.

The Deborah number may be used to explain also the differences in the evolution of the average pore pressure between the high-permeability simulations D3 (with De ≪ 1) and the low-permeability simulation D9 (with De 1), which are observed in figure 8. Upon the onset of the simulations, dilation occurs that increases pore space. In simulation D3, the small Deborah number allows for an immediate compensation of the newly generated pore space by fluid inflow from the drained boundaries, because the rate of PP diffusion across the system is fast with respect to the rate at which pore volume expansion occurs, due to the high permeability. For that reason, a transient effect of PP reduction is not observed (figure 8c). In simulation D9, the larger Deborah number causes a delay in fluid inflow that compensates for the expanding pore space, leading to a transient negative pore pressure (figure 8d).

An additional result of the drained simulations is the observed correlation between the average pore pressure and the apparent friction, as depicted in figues 11 and 13. Such correlation is predicted by equation ( 2). However, it is important to note that the law of effective stress, equation (1), and the relation between pore pressure and shear resistance as appears in equation ( 2), apply either only locally, or when the pore pressure is uniform within a granular or porous layer. Therefore, application of equations ( 1) and ( 2) to a meso-scale implicitly assumes uniform pore pressure, a situation that only arises in quasi-static and undrained conditions. In the drained simulations, the pore pressure is not uniform within the layer, but it varies dynamically in space and in time during shear. Still, the correlations between the average pore pressure and the shear resistance in figures 11 and 13 are good, indicating that although equation (2) may not be valid quantitatively for drained conditions, the concept of a meso-scale effective stress is still useful when the boundaries are drained, albeit, it merits further investigation.

Simulations of liquefaction events with drained conditions

After studying the general evolution of pore pressure and shear resistance in the drained simulations, we turn to discuss special events that are characterized by high pore pressure and low shear resistance. In the drained simulations D11 and D12 the average pore pressure has peaks that are of the order of magnitude of σ n , with zones that experience pore pressure that significantly exceeds σ n . During such events the effective stress vanishes, and the shear resistance becomes negative (figures 11 and 13). We define these events of µ a ≤ 0 as liquefaction. As long as P < σ n throughout most of the grains layer, granular stress chains support the external normal load by transmitting stress from top to bottom (figures 12a and 14a ). During liquefaction events, stress chains disappear in zones that have P > σ n (figures 12b and 14b), and the external load can no longer be transmitted through the granular phase. Instead, it is supported by the highly pressurized pore fluid (figures 12e and 14e).

The localization of high pore pressure close to the drained boundaries during the liquefaction event in simulation D11, (figures 12b and 12e), may be explained by its close to 1 De. In simulation D11 high pore pressure is not generated within the center of the layer since the De 1 limits the distance of pore fluid flow during the time scale of deformation. As a result, communication by pore pressure diffusion between the drained boundaries and the layer's center is uncommon. In simulation D12, the system half size, ζ, is smaller, and as a result, the Deborah number is smaller. For that reason, the pore pressure is high throughout the whole layer during the liquefaction event, and is not localized close to the boundaries (figures (14b) and (14e)).

Such transient liquefaction events may have significant implications to natural systems that are driven by dynamic forcing, such as tectonic loading. The simulations presented here are driven kinematically, i.e. a constant shear velocity is imposed on the top wall, and the systems cannot accelerate. It is speculated that if the systems would have been driven dynamically, by shear stress boundary conditions, then the short liquefaction events would have generated instabilities leading to acceleration. Following this initial acceleration, the absolute value of the porosity strain rate is expected to increase, potentially leading to even larger pressurization and further acceleration.

Previously, the only mechanism that was suggested to induce liquefaction was compaction of loosely packed grains under undrained conditions [START_REF] Sawicki | Developments in modeling liquefaction of granular soils, caused by cyclic loads[END_REF], leading to an elastic-like pore pressure evolution, (as modeled in our undrained loose simulation, U14). Therefore, the occurrence of liquefaction events with densely packed drained conditions, leading to viscous-like evolution of pore pressure, is a surprising result of our model (that was also obtained in our approximate model [Goren et al., 2010]). Such conditions are traditionally believed to be unfavorable for liquefaction [START_REF] Seed | Pore-water pressure changes during soil liquefaction[END_REF] despite field and experimental evidence that suggest otherwise. Such evidence are reviewed next.

Field and experimental evidence of liquefaction events with drained conditions

A first field evidence for liquefaction of initially dense but well drained layers comes from the report of [START_REF] Gabet | The mobilization of debris flows from shallow landslides[END_REF] on debris flow mobilization from dense soils. [START_REF] Gabet | The mobilization of debris flows from shallow landslides[END_REF] find correlation between mobilization and fines/sand ratio, where soils with a small ratio are mobilized. Assuming that a small ratio of fines contributes to good drainage, then the emergence of debris flows may be attributed to viscous-like pressurization during shear deformation induced by gravity. In such a case, the short compactive stages that followed the initial dilative phase, a scenario that is reviewed in [START_REF] Gabet | The mobilization of debris flows from shallow landslides[END_REF], may have led to pore pressure rise, liquefaction, and mobilization to debris flow. A second example comes from the famous liquefaction event in Kobe, Japan, that followed the 1995 Great Hanshin Earthquake (M=6.9). [START_REF] Soga | Soil liquefaction effects observed in the Kobe earthquake of 1995[END_REF] reviewed the damage in the port facilities that were built on reclaimed islands. It was found that soils that were vibro-compacted were still liquefied, although they sustained less deformation. Such soils are expected to be highly dilative, and according to previous conventional wisdom, liquefaction resistance. It is suggested that the moderate liquefaction that was observed may be attributed to viscous-like pore pressure rise during instantaneous and short compactive phases that punctuated the overall dilative path of the vibro-compacted fill material.

Another evidence for liquefaction with drained boundary conditions comes from two sets of experiments conducted by [START_REF] Iverson | Dynamic pore-pressure fluctuations in rapidly shearing granular-material[END_REF]. In the first set, a layer of initially compacted uniform rods immersed in water was sheared at a constant shear velocity, while fluid was allowed to drain out of the top boundary. During a significant portion of the experiment, the 'rods above the slip surface lost contact with the underlying rods and glided on a cushion of water', i.e. 'the pressure gradients between the shearing layer and the top boundaries were high enough to support the rods weight' [START_REF] Iverson | Dynamic pore-pressure fluctuations in rapidly shearing granular-material[END_REF]]. In the second set of experiments, natural fluid-filled sand was laid on a tilted table and slid under gravity. At the onset of slide, the pore pressure declined, presumably due to dilation, but then, when motion became steady, pore pressure rose and started fluctuating, with pressure gradient transiently supporting the grains layer overburden. This description is very similar to our observation of the pore pressure evolution during drained simulations D11 and D12 that showed transient liquefaction events.

The evolution of pore pressure with undrained conditions

Here we discuss the pore pressure response to granular skeleton deformation with undrained boundary conditions. Figures 15b and16a show good correlations between the evolution of pore pressure and the evolution of porosity for our undrained simulations. Indeed, equation (34) predicts that with undrained conditions (when De -1 = 0) pore pressure evolution should be elastic-like, with ∆P that depends on ∆Φ. To further validate this relation, figures 15d and 16c compare between undrained simulation (U13 and U14) results and the analytical prediction of equation ( 34), for the relation between ∆P and -∆Φ, and show good fit when assigning the initial porosity for Φ in equation (34).

Simulations of liquefaction and hardening events with undrained conditions

Our simulations show that when the boundaries are undrained, the response of shear resistance to shear deformation depends on the initial packing. When the system is initially densely packed, then dilative shear causes the pore pressure to decrease with respect to its initial value (figure 15b). Momentary compaction events will only slightly increase the pore pressure, but if the porosity remains higher than its initial value, the pore pressure cannot increase above its initial value. The decreasing pore pressure causes the effective normal stress and shear resistance to rise. Such 'dilatancy hardening' may have important implications for nucleation of earthquakes along fault gouge [START_REF] Scholz | Earthquake prediction -physical basis[END_REF][START_REF] Scholz | Velocity anomalies in dilatant rock[END_REF][START_REF] Lockner | Dilatancy in hydraulically isolated faults and the suppression of instability[END_REF][START_REF] Samuelson | Shear induced dilatancy of fluid saturated faults: experiment and theory[END_REF][START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF], possibly retarding the onset of an earthquake instability. However, when the accumulating tectonic load will eventually reach the system shear resistivity, then the slip may potentially be more rapid. This may occur for example when the very initial slip damages the sealed boundaries, which allows for fluid flow into the gouge layer. As a result, the effective stress will decrease abruptly, and the tectonic shear stress will be far greater than the system resistivity, and runaway accelerating slip may develop. This is a plausible mechanism for dynamic weakening by fluid inflow. When the layer is initially loosely packed, then compaction occurs with shear, and pore pressure quickly rises. Equation (34) predicts that for the pore pressure to increase to the value of the normal stress, a relatively small change of porosity of ∆Φ = -2 × 10 -4 is required (for pore water). The volumetric strain in simulation U14 that results from such a small reduction of porosity is 0.027%. Such a small strain may not be measurable in the lab, and it may appear that liquefaction occurs without any volumetric strain. Simulation U14 shows that from the onset of liquefaction and onward, the pore pressure remains at a constant value that exactly balances the external load, P = σ n . This steady state is achieved because a normal force balance on the top and bottom walls is achieved, and further volumetric strain is prevented. In simulation U14, the stress chains that detach at the onset of liquefaction (figure 17b) do not rejoin. As a result, shear deformation is completely accommodated within the fluid phase, and the system loses its shear resistance as depicted by the zero apparent friction in figure 16b.

Experiments of liquefaction and hardening with undrained conditions

According to conventional understanding, poor drainage is believed to be a favorable condition for liquefaction. For that reason most engineering studies of the process of soil liquefaction used undrained boundary conditions. Empirical studies of both cyclic loading (simulating the passage of shear waves during earthquakes) [e.g. [START_REF] Seed | Liquefaction of saturated sand during cyclic loading[END_REF][START_REF] Peacock | Sand liquefaction under cyclic loading simple shear conditions[END_REF][START_REF] Finn | Sand liquefaction in triaxial and simple shear tests[END_REF], and continuous loading (simulating mass sliding under gravity) [e.g. [START_REF] Castro | Liquefaction and cyclic mobility of saturated sand[END_REF] have confirmed that loosely packed systems are more prone to liquefaction with undrained conditions. Furthermore, when cyclic loading was stress-controlled, sand layers showed an abrupt increase in their strain amplitude at the point of liquefaction. Following simulation U14, the increase of strain amplitude may be attributed to the detachment of stress chains during liquefaction, which transfer the support of the external load from the grains to the pore fluid, which exerts significantly smaller resistance to shear.

As for the mechanism of liquefaction, it was attributed to the tendency of loose soil to compact with drained conditions [START_REF] Sawicki | Developments in modeling liquefaction of granular soils, caused by cyclic loads[END_REF], a tendency that was observed in dry and completely drained experiments [e.g. [START_REF] Youd | Compaction of sands by repeated shear straining[END_REF]. However, to the best of our knowledge, volume reduction was never reported during undrained experiments, and the relation between pore volume reduction and pore fluid pressurization was not acknowledge with undrained conditions. The reason is probably the limitations in measuring the tiny strains associated with compaction under undrained conditions, which may be as small as tenth of per mill. Unlike experiments, even very small pore volume change may be easily measured during simulations. Indeed, pore volume reduction is measured in our loose undrained simulation U14, and was shown to lead to liquefaction with accordance to equation (34).

While undrained conditions have shown to cause liquefaction when loading a loose specimen, they are also believed to cause hardening by pore volume increase (dilation) and pore pressure decrease when loading a densely packed layer [START_REF] Rudnicki | Stabilization of rapid frictional slip on a weakening fault by dilatant hardening[END_REF]Chen, 1988, Scholz, 2002]. [START_REF] Moore | Slow episodic shear of granular materials regulated by dilatant strengthening[END_REF] performed stress-controlled shear experiments of dilative saturated granular layers. They report that shearing of fine-grained sediments produced smaller deformation velocity than coarse-grained sediments, presumably because the fine sediments contributed to poor drainage leading to pore pressure reduction and hardening (higher shear resistance) upon dilation. When [START_REF] Samuelson | Shear induced dilatancy of fluid saturated faults: experiment and theory[END_REF] performed a double direct shear experiment of well drained grain layers, they observed no hardening upon dilation due to the good drainage that allowed for an immediate pore fluid inflow into the newly generated pore volume, and prevented pore pressure reduction and layer strengthening. Indeed, simulation U13 confirms that when an undrained initially dense system is sheared, then dilation will cause pore pressure reduction and layer hardening (with increased aparent friction, figure 15c). When an initially dense but drained system is sheared, then hardening may occur transiently (as in the onset of simulation D9, figures 8d and 8f) if the internal permeability is relatively low (leading to De 1). If the internal permeability is high and De ≪ 1, hardening will not be observed at all (figures 8c and 8e).

Estimating liquefaction potential

Following the analysis presented so far we may attempt to estimate the potential for liquefaction with various boundary conditions and different parameters. A precondition for liquefaction is the occurrence of compaction. When drainage is poor the overall compaction matters, and when drainage is good the rate of compaction matters. Indeed, many engineering analyses of liquefaction focus on the compaction potential as a function of the initial packing [e.g. [START_REF] Castro | Liquefaction and cyclic mobility of saturated sand[END_REF] and of the applied cyclic strain [e.g. [START_REF] Vucetic | Cyclic threshold shear strains in soils[END_REF], where the applied cyclic loading presumably induces progressive compaction of a loosely or unevenly packed system [START_REF] Youd | Compaction of sands by repeated shear straining[END_REF]. In the current analysis of the pore pressure evolution mechanisms, parameters such as compaction potential or number of loading cycles are not accounted for explicitly. For that reason, when we estimate the liquefaction potential in the following, we do not account explicitly for the initial packing and for the amount of imposed shear strain. We do however introduce a statistical factor that accounts for the chances of getting high enough pore pressure in large enough zone to completely detach a layer of stress chains by the shear strain applied in the simulation. This statistical factor may be thought of as the statistical equivalent to empirical measures such as the number of shear cycles needed for liquefaction.

When estimating the potential for liquefaction, one should first calculate the system Deborah number in order to decide which of the pore pressure evolution regimes will dominant. Figure 20 summarizes the various options. If De ≪ 1 and viscous-like evolution of pore pressure is expected to dominate, a first order approximation of the liquefaction potential is possible via the use of equation ( 36). Equation ( 36) gives the expected pore pressure as a function of the temporal derivative of the porosity and the system parameters. Approximating d Φ(ζ, t) /dt during compaction as -V sh /ζ, and accounting for the pore pressure at the center of the system, z = 0, equation (36) becomes:

P (z, t) ≈ ηV sh ζ 2k 0 . (38)
Dividing equation ( 38) by σ n gives a non-dimensional pressure ratio that expresses the liquefaction potential:

LP d = ηV sh ζ 2k 0 σ n . ( 39
)
When LP d ≪ 1, liquefaction is not expected because the generated pore pressure due to compaction is significantly smaller than the applied external stress, σ n . When LP d ≥ 1 liquefaction is possible, because the evolving pore pressure may reach and even surpass the value of σ n . Calculating LP d from equation ( 39) for the simulations it is found that in simulations D1 -D5, LP d < 1 and indeed no liquefaction is observed. In simulations D6 and D7, LP d ≈ 1 but still the pore pressure maxima are observed to be smaller than σ n . For simulations D8 -D9, LP d > 10, and pore pressure maxima surpass σ n . However, these maxima occur in highly localized zones, a situation that is shown to be insufficient for causing a complete loss of shear strength and liquefaction. In simulations D11 -D12, LP d > 100, and the pore pressure is shown to surpass σ n in relatively large volumes and to cause liquefaction (figures 12e and 14e). Note that simulation D10 is not accounted for here because it has De = 1. It is concluded that LP d > 1 is not a sufficient condition for liquefaction, however it is still expected that larger LP d will lead to larger chances of liquefaction. Our simulations show that when LP d > 100, (and when De < 1, where equations ( 36) and (38) are valid) liquefaction occurs (Simulations D11 and D12). It is proposed that the threshold of LP d > 100 comes about as a statistical measure for the chances of getting fast enough compaction in large enough area during the course of shear strain application. In that sense it is possible to rewrite equation (39) as

LP d = λ ηV sh ζ 2k 0 σ n , (40)
where λ ≪ 1 is a newly introduced statistical factor. The measure of LP d in accordance to equation ( 40) with λ = 0.01 exactly distinguishes between simulations that do not generate liquefaction and are characterized by LP d < 1, and simulations that show liquefaction and are characterized by LP d > 1 (see table 1). Yet it is not completely clear wether and how λ scales with system size and system dimensionality.

To estimate the liquefaction potential of field cases consider for example a layer of saturated soil buried at depth 10 m (this is also the distance to the boundary, ζ), with permeability of 10 -10 m 2 . For liquefaction to occur, the excess pore pressure should reach the initial effective stress at depth of 10 m, which is ∼0.15 MPa. Following equation ( 40) with λ = 0.01, for LP d > 1, the peak ground velocity (PGV) should be V sh > 0.3 m/s. Indeed, [START_REF] Kostadinov | Assessment of liquefaction-inducing peak ground velocity and frequency of horizontal ground shaking at onset of phenomenon[END_REF] estimated that the minimal PGV that may generate soil liquefaction is 0.1 m/s. This observation suggests that indeed also in the field λ may be approximated as 10 -2

When De ≫ 1, the evolution of the average pore pressure in the system is elastic-like and it follows equation ( 34). Here, liquefaction potential may be estimated with the non-dimensional pressure ratio:

LP u = - ∆Φ βΦ(1 -Φ)σ n . ( 41
)
Liquefaction is expected when LP u ≥ 1. For our undrained initially loose simulation U14, figure 16a shows that a porosity reduction of ∆Φ = -2 × 10 -4 has led to liquefaction. Such a porosity reduction gives LP u 1 in equation ( 41). Note that a statistical factor is not required here because equation ( 34) is for the average pore pressure, and since in simulation U14 the pore pressure homogenizes fast within the system, LP u ≥ 1 means that the pore pressure is of the order of the applied normal stress throughout the whole system.

To estimate the liquefaction potential of a field case with De ≫ 1, consider for example a thin gouge layer within a fault zone that is buried at depth of 1 km, and is bounded by undrained blocks. Still it is assumed that initially the pore pressure within the gouge is hydrostatic. For liquefaction of the gouge, the pore pressure should reach a value of σ n = 15 MPa. Assuming the initial porosity of the gouge is 0.1, then according to equation (41), in order to achieve LP u = 1, the reduction of porosity should be ∆Φ = -6 × 10 -4 , which results in a very small volumetric strain of less than 1%.

When De ≈ 1, the evolution of pore pressure is expected to be controlled both by viscouslike and elastic-like behaviors. Following the approximated model presented in Goren et al. [2010, their Appendix B], it is suggested that under such conditions the pore pressure in these mixed mode systems evolves as:

P (z, t) ≈ ηΦV sh d πβk 0 . (42)
For our simulation D10 with De = 1, equation ( 42) predicts P = 233 MPa. This prediction gives an order of magnitude approximation under the assumption that deformation is localized at the system's center [Goren et al., 2010, their Appendix B]. In simulation D10 the deformation zone is more widespread and the maximum pore pressure that was generated is 80MPa. To estimate liquefaction potential, equation ( 42) is divided by σ n and the statistical nature of the pore pressure evolution is accounted for using the prefactor λ, similarly to equation 40:

LP m = λ ηΦV sh d πβk 0 σ 2 n . (43)
For simulation D10, and using again λ = 0.01, equation ( 43) gives LP m ≈ 0.1. Indeed, high pore pressure that surpasses σ n is generated, but only in localized zones, and a complete loss of shear strength is not observed in this simulation.

Conclusions

In this work, a fully coupled model for the mechanics of fluid-filled granular media is developed from two components: A continuum formulation that describes the evolution of pore pressure in response to granular matrix deformation, and a granular dynamics algorithm, that solves the grain dynamics. The resulting fully coupled model is extremely general, as it is capable of simulating various loading scenarios leading to both reversible and irreversible granular matrix deformation, with various drainage conditions for the pore fluid. Furthermore, in the coupled model we do not explicitly implement the law of effective stress on the micro (granular) level, but apply on the grains forces that arise from pressure gradient in the pore fluid . Yet, we have validated that the effective stress law arises macroscopically. To the best of our knowledge this is the first time that the micro origin of the macro effective stress law was demonstrated. Analysis of the pore fluid formulation together with simulation results reveals that the evolution of pore pressure may be described as having two end-member behavior, 'viscous-like' and 'elastic-like'. These two behaviors control pore fluid pressurization and depressurization and the dominant terms in the pore pressure equation. The choice of which mechanism dominates depends on the Deborah number, De, which determines whether the system is effectively drained or undrained. When drainage is good (De ≪ 1), pore pressure evolution is viscous-like as it is a function of the volumetric strain rate (pore volume strain rate), and it depends on the fluid viscosity and the inverse of the permeability. When drainage is poor (De ≫ 1), pore pressure evolution is elastic-like, with pore pressure variations being a function of the overall volumetric strain (pore volume strain). Here, the pore pressure is a function of the fluid compressibility. Depending on the systems parameters and the boundary conditions, pore pressure evolution may follow one of these end-members, or be a mixture of the two.

Simulations of fluid-filled granular layers under constant normal stress and constant shear velocity reveal that pressurization and liquefaction may occur also in initially densely-packed layers as long as the boundaries are drained. Such conditions were previously often believed to be resistant to liquefaction [e.g. [START_REF] Seed | Pore-water pressure changes during soil liquefaction[END_REF]]. Here we show that liquefaction events can occur under such conditions, since viscous-like pore pressure evolution (that arises when some drainage exists) is a function of the instantaneous rate of change of porosity, and it has 'no memory' of the initial void ratio of the layer. Simulations with initially denselypacked undrained boundaries show 'dilatancy hardening' with pore pressure reduction and an increase of the layer resistance to shear. Shear of loose initial-packing under undrained conditions leads to a steady-state liquefaction upon very small volumetric strain (that may not be measurable in the lab).

Finally, we conclude the manuscript by addressing the two questions that were posed in the introduction: (1) what is the physics behind the pore pressure control over the shear strength? and (2) what processes alter the pore pressure? To answer the first question we have seen that when pore pressure rises to the value of the applied normal stress, then the force exerted by pressure gradients across the grains may be large enough to counter-balance the solid stresses, and thus acts to detach stress chains and separate previously contacting grains. When a large enough region experiences large pressure gradient, frictional resistance Integrating equation ( 46) between the center of the layer at z = 0 and some distance z from the center (see figure 7) results in: where Φ(z, t) is the average porosity between the system's center and distance z from it. Equation ( 47) then leads to the relation

u sz (z, t) = u sz (0, t) + ∂ Φ(z, t) ∂t z. (48)
Assigning equation (48) in equation ( 45) results in:

∂P (z, t) ∂z = η k 0 ∂ Φ(z, t) ∂t z + C 1 (t). ( 49
)
Integrating equation ( 49) between the layer's center and distance z leads to:

P (z, t) = P (0, t) + η k 0 d Φ(z, t) dt z 2 2 + C 1 (t)z, ( 50
)
where the rate of change of the average porosity, d Φ(z, t) /dt, is approximated as uniform in space. Requiring complete drainage across the boundaries, i.e. P (ζ, t) = P (-ζ, t) = 0, equation (50) leads to:

P (z, t) = - η 2k 0 d Φ(ζ, t) dt ζ 2 -z 2 .
(51)

Figure 1 :

 1 Figure

Figure 2 :

 2 Figure 2: Setup of a numerical experiment designed to verify that the grains-fluid model reproduces correctly the law of effective stress.In each experiment, a unidirectional stress, ∆σ n , is applied, and a constant pore pressure, P , is maintained. The vertical strain, ∆ǫ, is measured. The white pore area in between the grains is filled with fluid.

Figure 3 :Figure 4 :

 34 Figure 3: Effective stress validation test. Dry simulations are marked by 'x' and wet simulations by 'o'. (a) Volumetric strain, ∆ǫ, is plotted against the change of confining stress, ∆σ n , and show no correlation for the wet simulations. (b) Volumetric strain is plotted against the effective stress showing linear relation with the same slope as dry simulations.

Figure 5 :

 5 Figure5: Fluidization validation test. Comparison between theory and simulations for the minimal pressure gradient, ∇P , required to fluidize a layer of grains with porosity Φ under gravity g. ūs is the average grains velocity. ρ s and ρ f are the densities of the solid grains and of the fluid, respectively.

  well drained systems with k i ≪ k b 2ζk b /δ for poorly drained systems with k i ≫ k b 0 for completely undrained systems with k b = 0.

Figure 6 :

 6 Figure6: Two mechanisms control the evolution of pore pressure in equation (18). (a) When the boundaries are undrained (marked by double solid lines) leading to De ≫ 1, pore pressure responds elastically to any strain of pore volume. Compaction will lead to pressurization and dilation will lead to depressurization. The magnitude of pore pressure change depends on the overall porosity change, ∆Φ, and on the inverse of fluid compressibility, β -1 . (b) When the boundaries are drained (marked by dashed lines) and De ≪ 1, pore pressure evolves viscously in response to instantaneous strain rate of porosity, dΦ/dt. Upon dilation (left to center) fluid will flow into the system. In order to facilitate this flow, pressure gradient must arise with low pressure within the system interior. Upon compaction (center to right), fluid escapes from the system, and an opposite pore pressure gradient arises with higher pore pressure in the system interior. In this mechanism the evolving pore pressure depends linearly on fluid viscosity, η, and inversely on the permeability, k 0 .

  Applied drainage B.C. drained: P = 0 k b = INF undrained: ∂P/∂z = 0 k b = 0 Wrap around B.C.

  ,x) k i (t,x) P(t,x) μ a (t,z = ζ )

Figure 7 :

 7 Figure 7: Grains-fluid simulations setup. In each simulation a collection of grains are packed within a rectangular box with wrap around boundary conditions along the horizontal direction. The equation for the pore pressure, equation (15), is solved on a super-imposed grid. Normal stress, σ n , and shear velocity, V sh , are applied and maintained constant. The spatial and temporal evolution of porosity, Φ, and of pore pressure, P , and the temporal evolution of the apparent friction µ a are measure. k i is the internal permeability that is set by the local porosity with accordance to equation (24). k b is the boundary permeability that expresses the drainage boundary conditions. When the top and bottom boundaries are drained, P = 0 on the boundaries as if k b = INF. When the top and bottom boundaries are undrained, no pressure flux arises accros the boundaries as if k b = 0.

Figure 9 :Figure 10 :

 910 Figure 9: Results of drained simulation D7. (a) dΦ/dt (blue) and average pore pressure (green) for simulation D7 are plotted as a function of time to show good correlation as predicted by equation (36). (b) Φ (blue) and average pore pressure (green) show no correlation. Note that P increases downwards.

Figure 12 :Figure 13 :Figure 14 :

 121314 Figure 12: Snapshots of the granular configuration (top) and corresponding pore fluid pressure map (bottom) for drained simulation D11, before, during, and after the liquefaction event shown in figure 11. The times when the snapshots are taken are indicated by the framed letters in figure 11. Grains color code corresponds to the overall contact-induced compressive force they sustain. Warm colors indicate high force and cold colors low force. Contact forces are depicted by lines that connect contacting grains. The width of a line correlates to the magnitude of the normal force along the contact.

Figure 15 :Figure 16 :

 1516 Figure 15: Results of undrained initially dense simulation U13. (a) No correlation is found between dΦ/dt (blue) and average pore pressure (green). (b) Good correlation is found between porosity (blue) and average pore pressure (green). (c) Correlation is shown between apparent friction (blue) and average pore pressure (green). (d) Correlation between -∆Φ and ∆P in the simulation (circles) are well matched by the linear relation predicted by equation (34) (red line). Note that in (a), (b), and (c) the average pore pressure increases downwards.

Figure 17 :

 17 Figure 17: Snapshots of the granular configuration (top) and corresponding pore fluid pressure map (bottom) for undrained initially loose simulation U14, before, and during the liquefaction event depicted in figure 16. The times when the snapshots are taken are indicated by the framed letters in figure 16. For the color code of the grains see the caption of figure 12.

F

  Figure 18: Log-log plot of the slope of the relation P vs. -dΦ/dt as a function of the permeability scale factor, k 0 , for drained simulations. Theoretical prediction according to equation (36) appears as a solid line and simulations result as circles.

Figure 20 :

 20 Figure 20: Diagram suggesting the path to estimate liquefaction potential for grains-fluid systems under shear.

Table 1 :

 1 Numerical simulations. # a σ n [MPa] ζ [m] b k 0 [m 2 ] c V sh [m/s] De d

	LP e

a In the simulations numbering, D stands for drained and U stands for undrained.

Acknowledgements

RT acknowledges the support of the CNRS INSU program, the regional REALISE program, the ANR SISCA program and the european SAFELAND program.

to sliding of the layer is suppressed and shear is accommodated within the pressurized fluid phase. The answer to the second question is that grain compaction causes the pore pressure to rise and grain divergence causes the pore pressure to decrease. The magnitude of pore pressure change depends on both the volumetric strain rate under well drained conditions and on the absolute volumetric strain under undrained conditions.

Appendix A Pore fluid pressure evolution for De ≪ 1

In this section, the evolution of pore pressure is studied for drained boundaries with De ≪ 1. Under such conditions, the time dependent term in equations ( 15) and ( 18) becomes negligible with respect to the diffusion term because De ≪ 1 in the non-dimensional equation (30). Equation ( 15) then becomes

Formulation similar to equation ( 44) is developed by [START_REF] Iverson | Differential-equations governing slip-induced pore-pressure fluctuations in a water-saturated granular medium[END_REF] for drained conditions. For the 1D case, after integration, equation ( 44) becomes:

where C(t) is an integration factor, k(z, t) is approximated as the permeability scale factor, k 0 , and u sz is the horizontally averaged z component of the solid grains velocity. In order to express the pressure as a function of the temporal derivative of the porosity, ∂Φ/∂t as in equation ( 18), we use the 1D form of equation ( 17):