A General Criterion for Liquefaction in Granular Layers with Heterogeneous Pore Pressure

A General Criterion for Liquefaction in Granular Layers with Heterogeneous Pore Pressure

Liran Goren, Renaud Toussaint, Einat Aharonov, David W. Sparks, and Eirik Flekkøy

ETH Zurich, Switzerland; email: liran.goren@erdw.ethz.ch

Institut de Physique du Globe de Strasbourg (IPGS), UMR 7516, Université de Strasbourg/EOST, CNRS; 5 rue Descartes, F-67084 Strasbourg Cedex, France; email: renaud.toussaint@unistra.fr

Hebrew University, Israel; email: einatah@cc.huji.ac.il

Texas A&M University, TX, USA; email: sparks@geo.tamu.edu

University of Oslo, Norway; email: e.g.flekkoy@fys.uio.no


Fluid-saturated granular and porous layers can undergo liquefaction and lose their shear resistance when subjected to shear forcing. In geosystems, such a process can lead to severe natural hazards of soil liquefaction, accelerating slope failure, and large earthquakes. Terzaghi’s principle of effective stress predicts that liquefaction occurs when the pore pressure within the layer becomes equal to the applied normal stress on the layer. However, under dynamic loading and when the internal permeability is relatively small the pore pressure is spatially heterogeneous and it is not clear what measurement of pore pressure should be used in Terzaghi’s principle. Here, we show theoretically and demonstrate using numerical simulations a general criterion for liquefaction that applies also for the cases in which the pore pressure is spatially heterogeneous. The general criterion demands that the average pore pressure along a continuous surface within the fluid-saturated granular or porous layer is equal to the applied normal stress.


Dynamic loading can induce liquefaction when applied to fluid-saturated granular or porous media. From the macro-mechanical perspective, liquefaction means that the medium transitions from a solid-like material that can sustain some shear stresses without irreversible deformation to a fluid-like material that has no shear resistance. As a consequence, even the smallest shear stresses are relaxed by irreversible flow deformation. In the upper crust of the Earth, liquefaction can occur in response to shear loading, and it is the cause of many natural hazards. Soil liquefaction can take place in response to the passage of earthquake induced seismic waves through a soil column. It causes the soil to lose its shear resistance [Das, 1993, Kramer, 1996], and as a result, the soil can no longer support the foundations of infrastructures and catastrophic collapse might take place. When slope material collapses under gravity, liquefaction can lead to further acceleration of the sliding mass since the gravitational forces are no longer opposed by the frictional resistance of the sliding body [Iverson and Lahusen, 1989]. When fault zones slide during an earthquake, motion is localized along a shear zone that is in many cases filled with granular gouge. Liquefaction of the gouge layer can further enhance the earthquake instability, promoting a longer and faster event [Blanpied et al., 1992].

Terzaghi’s effective stress principle [Terzaghi, 1943] is commonly invoked to explain the process of liquefaction in such Earth systems. Terzaghi [1943] understood that it is not the applied stress that controls the strength of a solid-fluid system, but instead a quantity termed the effective stress:


where is the applied stress tensor, is the pressure experienced by the fluid within the pores of a granular or porous material, is Kronecker’s delta, and is the effective stress tensor. According to this principle, the shear resistance, , of a fluid-filled granular or porous layer is:


where is the total applied stress normal to the layer, is a friction coefficient, and it is assumed that the material is cohesionless. As a consequence, when , , the layer has no shear resistance, and is will deform irreversibly in response to shear stresses. Such a granular layer is said to be liquefied.

An important assumption, that is commonly not stated explicitly, is that the pore pressure in the layer is uniform. However, under dynamic conditions, such as simple shear forcing, and when the internal permeability of the layer is relatively small, the pore pressure is not necessarily spatially uniform. In such a case, it is not clear how Terzaghi’s effective stress principle, equation (1), the shear resistance of the layer, equation (2), and the criterion for liquefaction should be defined.

In this work we develop a general theoretical criterion for liquefaction that depends on the pore pressure and applies also when the pore pressure is spatially heterogeneous. We further demonstrate the validity of this criterion using simulations of saturated granular material under forcing of shear at a constant velocity. The general criterion states that a sufficient and a necessary condition for liquefaction of a granular layer that is finite in at least one dimension is that there exist a surface within this layer that separates its boundaries, and along this surface the average pore pressure is equal to the applied normal stress on the layer’s boundaries.


For most applications, it is sufficient to consider a granular layer that is finite in at least one dimension, i.e. it is bounded in two opposite directions, and the normal stress on its boundaries is uniform. For simplicity, we assume that the two boundaries are perpendicular to the direction, along which gravity operates. The layer is further assumed to be very long or even periodic (and thus infinite) along the horizontal dimensions, and . We start by writing again the total stress tensor in a fluid-filled granular medium over a representative element volume (REV), with negligible fraction of the surface transmitting the solid stress [Biot, 1956]:


where is the solid stress tensor associated with the forces transmitted between neighboring grains through solid contacts, and is the total stress. Note that and can be identified with the effective stress tensor and the stress tensor from equation (1), respectively. We use here the notation of total and solid stresses to emphasize the different physical components in the layer. The conservation of momentum statement over a small REV is:


where , is the bulk density of the solid grains, is the density of the interstitial fluid, and is the porosity. is the gravitational acceleration, and is the total acceleration. Combining equations (3) and (4) leads to:


We further decompose the pressure, , into a hydrostatic component, , and a non-hydrostatic component, . Equation (5) then becomes:


A general liquefaction criterion

Next, we use the force balance, equation (6), together with some assumptions in order to derive a simple and general criterion for liquefaction. The first assumption is that the inertia term, the right hand side of equation (6), is relatively small with respect to the terms on the left hand slide. This assumption is mostly valid for slow granular flow, and we revisit it when we analyze the numerical simulations. Equation (6) then reduces to a balance of forces:


To motivate the next stages in the derivation, we note that during liquefaction, shear is not opposed by resisting forces that arise from grain contacts. Therefore, when shear stress or strain is applied to one (or both) of the layer’s boundaries, the boundary can slide with no shear resistance only if there exists a continuous surface that separates the two boundaries along which there are no forces arising from grain contacts, i.e. stress chains do not percolate from one boundary to the other. The need for such a continuous surface arises from the non-local nature of the rigidity of granular media, which means that a local region with no solid contact forces can be bypassed by stress chains around it that support the system and maintain resistivity to shear.

In the following, we will consider a planar surface parallel to the boundaries instead of a general surface. This assumption simplifies the derivation, and it is mostly valid, as we observe in numerical simulations that the separating surfaces responsible for liquefaction tend to be subhorizontal. Again, we discuss the validity of this assumption and its implications in the section dealing with the numerical simulations, below. Since only planes that are perpendicular to the direction are considered, we focus on the vertical force balance of equation (7):


Then, we average equation (8) along a horizontal plane:


where is the average pore pressure along a plane, , whose area is , and is defined in a similar way. in equation (9) refers to the porosity above the considered plane. Note, that upon integration, shear stresses on vertical planes are assumed to average out, so the second and third terms in equation (8) become zero. This is the result of the very long or periodic nature of the layer in the and directions. Equation (9) is then integrated along the direction to give:


Where is an integration constant.

The final assumption that we adopt here is that the gravity term, the third term on the left hand side of equation (10), is relatively small with respect to the total stress. Such an assumption is valid as long as the layer is thin with respect to the overburden that causes the total stress. Using this assumption, , where the normal stress on the boundaries of the layer. Equation (10) then reduces to:


where the prime symbol is dropped from the as we neglected the hydrostatic contribution to the pore pressure. Note that equation (11) is similar to equation (3), but is valid for quantities averaged over a horizontal plane rather than for a REV. Since solid stresses are always positive, the argument presented before, that liquefaction occurs when there are no solid contact forces across some continuous plane, can be expressed as along some plane at depth . Therefore, equation (11) leads to a necessary and sufficient conditions for liquefaction in terms of the pore pressure:


which means that along some horizontal plane located at depth the average pore pressure is equal to the applied normal stress over the boundaries.

Numerical Approach - a Coupled Grains and Fluid Model

Numerical simulations are used to better illustrate the conditions that are related to liquefaction and in order to validate the theoretical general liquefaction criterion that we develop herein. Simulations are performed with a coupled model that combines discrete element method for the granular phase and a continuum fluid solver for the interstitial fluid phase. The details of the coupled approach are discussed in [Goren et al., 2010, 2011]. For clarity, we briefly review it here. The discrete element component for the solid granular phase implements a 2D granular dynamics (GD) algorithm, in which each individual grain is treated as dynamically tracked particle, subjected to a variety of external forces. Grain interactions, body forces, and the force induced by the interstitial fluid lead to linear and rotational acceleration of the grains:


where and are the translational and rotational velocity vectors of grain (a superposed dot indicates time derivative). is the grain mass, is the gravitational acceleration, is the grain moment of inertia, is the radius of the grain, and is a unit vector normal to the contact between grains and . The last term on the right-hand side of equation (13) refers to the force exerted on grain by the pressure gradient, , of the fluid surrounding it, normalized by the solid fraction, , in its vicinity, and by its volume, [McNamara et al., 2000]. refers to the inter-granular contact force between grains and [Goren et al., 2011]. Equations (13) and (14), which are solved for each individual grain, give rise to the emergent dynamics of the granular layer.

The formulation for the fluid component is based on mass conservation statements for the solid and fluid phases, Darcy’s law, which is a reduced form of the fluid conservation of momentum equation under the assumption of negligible fluid inertia, and a state equation for a compressible fluid. Combined together, a simple three terms equation can be written:


where and are the compressibility and viscosity of the pore fluid, respectively. The permeability, , is calculated with a Carman-Kozeny relationship, modified for a 2D volume fraction formulation [McNamara et al., 2000], and is the solid velocity field averaged at a scale over which Darcy’s law applied. Equation (15) is derived under the assumptions that (1) the compressibility of a grain is negligible relative to the fluid compressibility, (2) the pore fluid pressure is not too large, i.e. , and (3) the length scale of pore pressure diffusion is always larger than the diameter of a single grain. The fluid equation (15) is solved on a grid that is super-imposed over the granular layer, with grid spacing of about 2 grain diameters. The fluid and solid components interact by mutual interpolation, such that granular packing and grain motion induce interstitial fluid flow and pressure changes, and pore fluid pressure gradients affect the force balance of individual grains. This type of coupling was shown to reproduce pattern formation during sedimentation of grains in liquid [Niebling et al., 2010a, b] or gas [Vinningland et al., 2007a, b, Niebling et al., 2010b], and mechanisms of fracture due to fluid injection [Johnsen et al., 2006, 2007, Niebling et al., 2012].

We use the coupled grains-fluid model to perform simulations of shearing a saturated granular layer at a constant shear velocity. Simulations are performed in a rectangular domain with approximately 1680 (24 70) grains. Grain diameters are drawn randomly from a Gaussian distribution with an average 1 mm, and a standard deviation 1 mm, clipped below 0.8 mm and above 1.2 mm. The simulations are performed at a constant and relatively high total normal stress, =2.4 MPa, such that the gravity in the thin layer is negligible. Although there is no gravity, we still define a vertical direction, and is applied to the two boundaries that are perpendicular to this vertical direction. The layer is periodic in the horizontal direction, and thus analogous to a rotary shear apparatus.

Each simulation is initiated by compacting a system of loosely packed grains under , until the porosity equilibrates. We then assume that the pore space is filled with water ( Pa s, Pa) at zero excess fluid pressure. Variations of pore pressure are measured relative to the initial zero value that corresponds to hydrostatic conditions. For this reason, is interpreted as the initial effective stress, under hydrostatic conditions. Finally, a constant shear velocity, m/s, is applied to the top wall. During a simulation, and are maintained constant, and we measure compaction and dilation in the layer, the shear stress required to shear the top wall at constant velocity, and the evolution of pore pressure. The macro-scale friction coefficient of the system, , is found by dividing the measured shear stress by the applied . We report here on two simulations that differ in their drainage boundary conditions: an undrained simulation, in which there is no fluid flux across the top and bottom boundaries, and a drained simulation, in which the pressure at the top and bottom boundaries is set to zero.

Liquefaction in Constant Shear Velocity Simulations

Liquefaction under undrained conditions

First we study an undrained granular layer with an initial high porosity. Such a system undergoes compaction during shear. The internal permeability within the layer is of the order of m. The time scale related to pore pressure diffusion throughout the layer is defined as , where 70 mm is the height of the granular layer. Due to the high internal permeability, sec, which is much smaller than the time scale related to shear deformation that is defined as sec. Therefore, the pore pressure equilibrates rapidly within the layer and can be considered uniform. Simulation results follow the theory developed in [Goren et al., 2010, 2011], and show that as the pore volume reduces with respect to the initial porosity, the pore pressure rises. Figure 1 shows that when the pore pressure becomes equal to the applied normal stress , the solid contact forces disappear, volumetric strain stops growing, the global friction, , drops to zero, and the fluid saturated granular system has liquefied. Note that here, because of the uniformity of the pore pressure, the classical interpretation of Terzaghi’s effective stress principle is valid, and the pore pressure at any point within the system or the average pore pressure within the layer could be used in equations (1) and (2), and in Terzaghi’s criterion for liquefaction, i.e. .


Figure 1: Pore pressure (solid line) and global friction (dotted line) for an undrained saturated layer of grains that compacts during shear. Figures to the left and right show the grain configurations and stresses at the beginning and end of this time period, respectively.

Liquefaction under drained conditions

Next, we turn to a drained simulation that starts from a smaller initial porosity and it has a smaller internal permeability of the order of m. The initial well-compacted configuration leads to an overall dilative trend upon shear, and the smaller internal permeability results in a longer pore pressure diffusion time scale, sec. As a consequence, the pore pressure cannot equilibrate within the granular layer during the time scale of shear deformation, which is 0.05 sec, and the pore pressure is heterogeneous (see Figure 2f). It might be expected that the dilative trend should prevent generation of high pore pressure, but as shown in Goren et al. [2010] and Goren et al. [2011], in drained conditions, , so that any small, but rapid, compaction that punctuates the overall dilative trend will lead to positive excess pore pressure. Indeed, despite the dilative trend, the macro friction coefficient, , reduces to zero or to a small negative value several times during the deformation, in correlation with rapid compaction events. Times with are transient liquefaction events. Figure 2d, shows the evolution of through time. Arrows mark the times when which we regard here as liquefaction. A finite value of 0.01 is used because the record of is not continuous during the simulation.

During shear deformation, we cannot define a macro-scale Terzaghi effective stress, and we cannot simply relate a value of to the liquefaction events because the pore pressure is highly heterogeneous. Instead, we test three different metrics for the spatially variable pore pressure in order to define a criterion for liquefaction. Goren et al. [2011] showed a good correlation between the average pore pressure within the system and the macro-scale friction coefficient, . However, as can be seen in Figure 2a, pore pressure averaged over teh whole layer, , never exceeds , even during liquefaction events; therefore is not a sufficient criterion for liquefaction. Alternatively, we can consider the maximum value of pore pressure within the system at any given time during shear. Figure 2b shows that the maximum value of pore pressure exceeds for almost 70% of the simulation duration, and therefore a criterion of the form is not a necessary condition for liquefaction.

As a third metric, we calculate along all horizontal planes within the system through time. Figure 2c shows the locations along the layer’s height where during the simulation. Arrows mark the times when such a condition is true and the friction curve drops to a value smaller than 0.01. The correlation between the occurrences of and is good, and supports the validity of our general liquefaction criterion, equation (12). For the times in which the friction is low but for all (light-color arrows in Figure 2d), we speculate that a non-horizontal surface, , for which we do not account here satisfies . For the times in which but the friction do not drop to low enough value, we offer that inertia is non-negligible, and therefore the general criterion is not entirely valid. Note that also when , the inertia is not negligible, but this does not affect the predictive nature of a general liquefaction criterion of the form .


Figure 2: (a)-(d) are time series of pore pressure observables and global friction for a drained grains-fluid simulation. (a) Average pore pressure. (b) Maximum pore pressure. (c) General liquefaction criterion. Locations in which are marked by plack line. (d) Global friction coefficient. Dark shading arrows show transient liquefaction events that are explained by the general liquefaction criterion. Light shading arrows are transient liquefaction events not explained by the general principle. (e) Granular packing and contact forces, and (f) Pore pressure map during the transient liquefaction event that is marked by the long arrow. Note that the pore pressure is very heterogeneous and that the solid contact forces disappear in localities where the pore pressure is high.


Terzaghi’s principle of effective stress can be cast as a criterion for liquefaction along a fluid-filled granular layer when, . However, this criterion is valid as long as the pore pressure is spatially uniform. Under dynamic loading and when the permeability is relatively small, the pore pressure is spatially heterogeneous and therefore the classical criterion is not useful because it doesn’t define where the pore pressure should be measured.

We develop here a general criterion that is suitable also for systems with heterogeneous distribution of pore pressure. The criterion states that a necessary and sufficient condition for liquefaction is the existence of a continuous surface along which the average pore pressure is equal or greater than the applied normal stress on the boundaries of the layer. According to this criterion, the value of pore pressure at local points or the average value of pore pressure as measured in experimental or natural settings are not sufficient for relating the pore pressure to loss of shear resistance. Instead, a full pore pressure map is required and an analysis of the general criterion is needed.


  • Biot [1956] M. A. Biot. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range. J. Acoust. Soc. Am., 28:168–178, 1956.
  • Blanpied et al. [1992] M. L. Blanpied, D. A. Lockner, and J. D. Byerlee. An earthquake mechanism based on rapid sealing of faults. Nature, 358(6387):574–576, 1992.
  • Das [1993] M. B. Das. Principles of Soil Mechanics. PWS-Kent, Boston, Mass., 1993.
  • Goren et al. [2010] L. Goren, E. Aharonov, D. Sparks, and R. Toussaint. Pore pressure evolution in deforming granular material: A general formulation and the infinitely stiff approximation. J. Geophys. Res., 115(B09216), 2010. doi: 10.1029/2009JB007191.
  • Goren et al. [2011] L. Goren, E. Aharonov, D. Sparks, and R. Toussaint. The mechanical coupling of fluid-filled granular material under shear. Pure and Applied Geophysics, 168(12), 2011. doi: 10.1007/s00024-011-0320-4.
  • Iverson and Lahusen [1989] R. M. Iverson and R. G. Lahusen. Dynamic pore-pressure fluctuations in rapidly shearing granular-material. Science, 246(4931):796–799, 1989.
  • Johnsen et al. [2006] Ø. Johnsen, R. Toussaint, K. L. Måløy, and E. G. Flekkøy. Pattern formation during air injection into granular materials confined in a circular Hele-Shaw cell. Phys. Rev. E, 74(1), 2006. doi: 10.1103/PhysRevE.74.011301.
  • Johnsen et al. [2007] Ø. Johnsen, R. Toussaint, K. J. Måløy, E. G. Flekkøy, and J. Schmittbuhl. Coupled air/granular flow in a linear Hele-Shaw cell. Phys. Rev. E, 77(1), 2007. doi: 10.1103/PhysRevE.77.011301.
  • Kramer [1996] S. L. Kramer. Geotechnical earthquake engineering. Prentice Hall, Inc., Upper Saddle River, New Jersey, 1996.
  • McNamara et al. [2000] S. McNamara, E. G. Flekkøy, and K. J. Måløy. Grains and gas flow: Molecular dynamics with hydrodynamic interactions. Phys. Rev. E, 61(4):4054 – 4059, 2000.
  • Niebling et al. [2010a] M. J. Niebling, E. G. Flekkøy, K. J. Måløy, and R. Toussaint. Sedimentation instabilities: impact of the fluid compressibility and viscosity. Phys. Rev. E, 82, 2010a. 10.1103/PhysRevE.82.051302.
  • Niebling et al. [2010b] M. J. Niebling, E. G. Flekkøy, K. J. Måløy, and R. Toussaint. Mixing of a granular layer falling through a fluid. Phys. Rev. E, 82, 2010b. 10.1103/PhysRevE.82.011301.
  • Niebling et al. [2012] M. J. Niebling, R. Toussaint, E. G. Flekkøy, and K. J. Måløy. Dynamic aerofracture of dense granular packings. Phys. Rev. E, 86, 2012. 10.1103/PhysRevE.86.061315.
  • Terzaghi [1943] K. Terzaghi. Theoretical Soil Mechanics. John Wiley, New York, 1943.
  • Vinningland et al. [2007a] J. L. Vinningland, Ø. Johnsen, E. G. Flekkøy, R. Toussaint, and K. J. Måløy. Granular Rayleigh-Taylor instability: Experiments and simulations. Phys. Rev. Lett., 99(4), 2007a. doi: 10.1103/PhysRevLett.99.048001.
  • Vinningland et al. [2007b] J. L. Vinningland, Ø. Johnsen, E. G. Flekkøy, R. Toussaint, and K. J. Måløy. Experiments and simulations of a gravitational granular flow instability. Phys. Rev. E, 76, 2007b. doi: 10.1103/PhysRevE.76.051306.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description