When loosely packed water-saturated granular soils, for example sands, are subjected to strong earthquake shaking, they may liquefy, causing large deformations with great destructive power. The phenomenon is quite general and occurs in any fluid-saturated granular material and is a consequence of the transfer of stress from inter-grain contacts to water pressure. In modern geotechnical practice, soil liquefaction is commonly considered to be an ‘undrained’ phenomenon; pressure is thought to be generated because earthquake deformations are too quick to allow fluid flow, which enables water depressurization. Here, we show via a first principles analysis that liquefaction is not a strictly undrained process; and, in fact, it is the interplay between grain rearrangement, fluid migration and changes in permeability, which causes the loss of strength observed in so many destructive earthquake events around the world. The results call into question many of the common laboratory and field methods of evaluating the liquefaction resistance of soil and indicate new directions for the field, laboratory and scale-model study of this important phenomenon.
The general phenomenon of the liquefaction of granular materials is applicable to a wide variety of grains and fluids. Here, we focus on one of the most important aspects of this phenomenon, namely the role of water flow in earthquake-induced liquefaction of sand near the ground surface. Our analysis uses a non-dimensional formulation and asymptotic analysis, which is adapted to the first few tens of metres of soil where the interplay of gravity, grain rearrangement, water compressibility and water flow combine to cause these destructive events.
In the geotechnical engineering field, soil liquefaction is commonly understood as a consequence of water pressure build-up owing to rapid squeezing of pore space, without sufficient time for water to flow through the grains and drain the pressure, e.g. [1–3]. When grains are loosely packed, during earthquake motion the tendency is for soil grains to move closer together, squeezing the water and rapidly increasing the pressure owing to the high bulk modulus of water.
Taken across a thin horizontal section of the soil, the total vertical force is the sum of the contact forces between grains and the water pressure times the cross-sectional area. This ‘total stress’ minus the water pressure is the so-called ‘effective stress’, commonly used in geotechnical analysis , which is a measure of the contribution of grain–grain interactions in the soil. If the contact forces drop to zero, then the water pressure carries the entire normal stress, and shear stresses induce flow in the manner of a viscous fluid.
The goal of liquefaction assessment up to now has been to determine how the water pressure in the soil will change during cyclic loading. The methodology employed has primarily been to build models based on the results of laboratory triaxial, hollow cylindrical and simple shear tests, as well as more extensive physical centrifuge models. Laboratory tests use samples of sand typically approximately 10–20 cm in characteristic size. The sand is surrounded by an impermeable, flexible membrane to trap the pore water, and the entire sample is contained in a pressurized vessel to simulate the overall pressure conditions in the ground. Cyclic loading of various forms is applied and the total and effective stress states, and evolution of pore water pressure are tracked through the cycling. These experiments have been extensively performed. To complement these small laboratory-scale tests, extensive physical modelling in various types of geotechnical centrifuge apparatus has been carried out. These tests use centrifugal acceleration to model the stresses induced in deposits of soil that are 10–100 times deeper than the model scale and allow full three-dimensional geometries to be simulated without the ambiguity of numerical simulations of complex soil materials. Combined with these physical observations, engineers have employed correlations with observed field conditions in post-seismic investigations. An overview of the current state of liquefaction research may be found in .
Implicit in the use of small laboratory samples with impermeable membranes and explicit in many textbook definitions of liquefaction is the assumption that during the earthquake there is no significant water flow or change in water volume, which is known as the ‘undrained’ condition (cf. [2,3]). The assumption that soil liquefaction occurs under undrained conditions has led to extensive research into undrained tabletop experiments, such as triaxial and simple shear tests, with several methods suggested to overcome the small volume changes allowed by the compliance of the rubber membrane . Although centrifuge and laboratory test data have long shown that water migration can occur [5,6], these situations have been treated as if they were exceptions to the normal situation of undrained pore pressure increase. However, given the thermodynamics of water pressure, the undrained assumption can lead to physically incorrect predictions. If the water is not allowed to change volume whatsoever, then the only mechanism for pressurization is heating. This shows that some care is required to determine properly the role of water flow, water compressibility and thermal expansion, which was the initial impetus for this research.
With the advent of large computing power, more recent studies of liquefaction phenomena have used discrete element models (DEMs) which operate at the grain scale and calculate the equations of motion for thousands of individually tracked cylindrical or spherical grains. In , a two-dimensional DEM model was coupled to a continuum model for fluid flow, and the interactions of grains and fluid were calculated for a sample of a few thousand grains. Their conditions are typical of 200 m- to 2 km-deep thin deposits sheared at ord(1) to ord(10) m/s1 which are relevant for fault gouge conditions. Holding boundary conditions of their sample at either zero fluid mass flow, or constant fluid pressure, they were able to observe liquefaction for both dense and loose sheared assemblies under drained and undrained conditions. While their method is, in principle, applicable to a wide variety of situations, computing the interactions of large deposits of sand over metres or tens of metres would be computationally prohibitive. Their model does show, however, that our understanding of even tabletop-sized experiments may be flawed, as they observe liquefaction under all conditions in assemblies whose bulk density was either relatively loose or dense. In an earlier paper , a non-dimensional continuum equation for dynamic fluid flow was derived which is valid for mesoscopic scales and based on conservation of mass together with Darcy's law. Their equation provides much of what is required to analyse a realistic soil deposit, though they explicitly neglect certain aspects, for example thermal heating, and they do not extensively analyse the equation in the context of typical near-surface liquefaction conditions.
In this paper, we analyse the liquefaction of saturated sandy and silty soils in the first ord(10) metres below the ground where initial water pressures are within a few atmospheres of the total vertical stress. Our method is to derive a one-dimensional equation for fluid flow in the vertical direction based on the assumptions of mass conservation and Darcy's law for fluid flow through porous media in a manner very similar to . Although the same derivation can be trivially extended to three dimensions mathematically, it provides no additional qualitative understanding for the main point of this study, which is that fluid flow and soil inhomogeneity are of critical importance in the liquefaction phenomenon. In the development of full predictive models for realistic heterogeneous three-dimensional soil deposits however, lateral water flow must be accounted for, as it will influence the water pressure as well.
By treating a comprehensive set of dependencies, including temperature effects and the compressibility of water potentially containing some small quantity of gas, such as that arising from organic processes or dissolved gasses, we arrive at a very general result, which leaves no further degrees of freedom that could be expected to contribute significantly to water pressurization in these soils. By non-dimensionalizing the model using a distinguished limit , we focus attention on the relative size of the various effects that can cause fluid pressurization and transfer of stress from soil particles to the water under the relevant circumstances. By showing that heating and compressibility of water are both perturbation effects, our results clarify that both grains collapsing into denser arrangements, and fluid migration via pressure diffusion are the dominant processes for earthquake liquefaction. Indeed, we show that for loose sands, diffusion is equally important to densification on timescales significantly shorter than the duration of an earthquake and even shorter than a single-loading cycle.
For layered deposits or samples with spatially varying permeability, gradients in the permeability field are critical in determining where liquefaction onset might occur. Though this has been shown in a variety of experiments, it has never been explicitly analysed mathematically to determine how important the effect is in general. Our model explains the dominant effects in a range of conditions from very low permeability soils typical of silts, through sands, and even in the case of high permeability gravels where diffusion is more dominant. One value of the model is that it predicts certain new features, such as large settlements together with bulk fluid flow to the surface in loose high permeability deposits, and also suggests a new mechanism of liquefaction potentially mediated by heating within silts even if the silt's initial density is such that sustained contractive loading does not occur. We compute solutions for several example scenarios that mimic the observed dynamics in tabletop [6,10] and centrifuge experiments .
Based on the results of these analyses, we believe that liquefaction evaluation approaches should explicitly acknowledge the importance of fluid flow. Future efforts should be placed on evaluating permeability, porosity and the energetic potential for earthquake-induced porosity changes over entire deposits of soil.
In addition to a description of the fluid flow, the prediction of the grain motions is of great interest, as it is a driving force for the fluid pressurization and migration. At present the authors are unaware of any continuum model, which could reliably substitute for the results of fully coupled DEM and yet be computationally tractable for large scales of spatially variable soil deposits. Because of this limitation, we apply an assumed grain deformation, usually a constant in space and time, and investigate the predicted results in terms of the induced fluid flow.
2. Model construction
We assume a horizontally layered deposit of soil and a vertical coordinate axis pointing upward with zero reference point at a deep impermeable layer below the region of interest. We choose this coordinate system rather than the one oriented downward from the ground surface because the ground surface may settle downward. We are interested in the first few tens of metres below the ground surface, as this is where small pore volume strains can produce changes in water pressure that can easily produce a zero effective stress condition. We consider a thin horizontal slice of this material with cross-sectional area A which is very large compared with the typical size of a grain, and thickness dz which is of the order of several grain diameters. For the saturated sand we consider, the fluid in our differential volume at pressure P and density ρ fills the inter-grain spaces whose volume is Vf=ϕA dz, where ϕ is the bulk porosity (the fraction of total volume occupied by pore space). Since the bulk modulus of the grains is typically an order of magnitude higher than that of the fluid, which is already quite high, we assume that the individual grains do not compress by any appreciable amount. However, grain rearrangement can cause changes in porosity.
While the usual assumption in liquefaction research is that fluid flow does not occur on the time scale of the earthquake, we instead model the rate of change of fluid mass within the volume of a thin slice by examining the difference in the fluid flux out of the top surface and into the bottom surface of the slice. We can estimate the fluid volumetric flux rate Q using Darcy's law: Q/A=−(kd/μ)(∂P/∂z). This expression defines kd mathematically, and the scientific content of Darcy's law is that kd is constant for a given grain structure over a wide variety of pressure gradient conditions. It has been frequently validated in experiments typically performed when the grains are stationary in a laboratory reference frame. To maintain Galilean invariance of the expression with respect to an arbitrary velocity of an observer, we add ϕvg, where vg is the apparent velocity of the grains and ϕ is the fraction of the cross section that contains water. This gives the correct fluid flux when a camera is passed by a static tube of soil and water at any non-relativistic velocity. Although Darcy's law represents an approximation of the full fluid flow equation when only viscous effects and pressure gradients dominate and a fuller equation would be required when grains are accelerating and moving rapidly relative to each other, it nevertheless can be assumed that this simplified expression also remains asymptotically valid when vg is sufficiently slowly spatially varying such as during settlement and liquefaction. This assumption will be borne out later by qualitative comparison with experimental results.
Assuming that the horizontal slice can be made thin enough relative to the length scale over which fluid velocity varies, the resulting differential equation for change in fluid mass within the differential slice is 2.1
Here, the left-hand side expresses the change in a short time dt of the total fluid mass in the slice. The right-hand side expresses the difference in fluid mass flux Q.
Because liquefaction of soils takes place over a very small range of fluid density and temperature conditions, we can eliminate the density ρ from the equation by solving for it in terms of pressure and temperature using the Taylor expanded state equation: 2.2where P0, ρ0 and T0 are, respectively, reference pressure, density and temperature levels for the water, KB is the bulk modulus of water (−V (∂P/∂VT)), α is the coefficient of volumetric thermal expansion ((1/V)(∂V /∂TP)), and ρ and T are the actual density and temperature, respectively. The temperature dependence can be derived from allowing a volume of water to be heated and expand under constant pressure, and then be squeezed back to its original volume at constant temperature. For our purposes, we neglect the very small difference between the adiabatic and isothermal bulk modulus of water.
The term ΔMstatic is a correction necessary to account for the fact that there is no change in the slice's water mass when static gravitational conditions apply. For example, owing to gravity there is a static pressure gradient. Also, the permeability may change with depth and water density changes minutely with depth. Setting the left-hand side of the equation equal to zero when static gravitational conditions apply and solving for the correction gives 2.3
In equations (2.1) and (2.3), for a vertical cartesian axis z, A, dz and dt are constant, kd depends on location z, μ depends on temperature T, ρ is determined from P and T using the state equation, and P, T and ϕ depend on time t and position z.
Note that we may compensate for the presence of a small fraction of the voids containing gas by adjusting KB. If a volume of water contains small bubbles occupying a small fraction ϵ of the volume, then we can approximate the compressibility (1/KB) using Vtot=VW+Vg and 1/KBmix=−(1/Vtot)(∂Vtot/∂P) together with the state equation Vg=KT/P for an ideal gas and the Taylor series state equation for the liquid Vw=V0(1−(P−P0)/KB). Assuming isothermal gas compression owing to the thermal equilibrium with the water, Vg=V0P0/P. In a mixture with ϵ fraction of gas, this leads to the equation: evaluating at P0, we get This is a strong function of ϵ as P0/KB is small for atmospheric-scale pressures. This may be relevant in situations involving gas bubbles formed by biological organisms, recent rains bearing an excess of dissolved gas or in tidal areas where fluxes of water and gas bubbles are large.
As we used the state equation involving P and T to eliminate ρ, our equation now includes the temperature T and in expanded form will include both ∂T/∂t and ∂T/∂z. We use the following heat equation for a flowing fluid moving at velocity vf to relate these quantities 2.4
In (2.4) cv is a bulk averaged specific heat capacity, ρavg is the bulk average density, kT is the effective thermal conductivity and ∂S/∂tQ is the change in entropy density owing to heat generation.
In addition, the conservation of grain mass was used to replace ∂vg/∂z in equation (2.1) with a term involving ∂ϕ/∂t and a convective term involving vg, which will become a perturbation owing to the small settlement velocities. 2.5
With these substitutions in place, equation (2.1) is not yet useful for analysis and calculation, as within each derivative there are multiple variables whose values are dependent on time or position. Our next task will be to expand the fluid equation (2.1) through the comprehensive application of the chain rule, so we can arrive at an equation for the rate of change of pressure ∂P/∂t and then make this equation dimensionless to help to evaluate the relative importance of various terms. Although it is tempting at this point to try to simplify this equation by making assumptions about the negligibility of certain effects, such as compressibility of water, the heat generation rate and similar effects, in order to reduce the complexity of the expanded equation, we believe that such early simplifying assumptions have led the field astray in general, and we instead choose to expand the equation fully and deal with the negligible terms only after a full dimensionless analysis. Details of this analysis are available in the electronic supplementary material.
We employed the freely available Maxima computer algebra system  to calculate the expanded version of equation (2.1). By identifying and dropping negligible terms, we arrived at a simplified ‘full’ equation, which captured all of the dynamics of interest including first-order perturbations, however, in order to arrive at our simplified equation, we needed to understand the dominant balances of effects and identify the negligible terms. To do this, a proper non-dimensionalization was required. The choice of the characteristic scales for non-dimensionalization is one of the critical aspects of our analysis, because it focuses the equation on the specific situations of interest.
To focus our understanding on geotechnically relevant scales, we choose scales for P, kd and ϕ by reference to the situation in the first several metres of sandy soil, where P/P0=O(1) throughout. We introduce kd0 which represents a typical order of magnitude for permeability of moderately coarse loose sands, such that the permeability range for most sands falls within about a factor of 4 of this value. We also introduce ϕ0, which represents a typical value for porosity in a loose sand. General sand is between about 70 and 110% of this ϕ0 value. These values were chosen as follows by reference to tables given in . 2.6
We wrote equations for the dimensional variables in terms of non-dimensional ‘primed’ or ‘hatted’ variables as follows: 2.7Here, Cv0≈1594 J kg−1 K−1 represents a mass weighted average of the specific heat capacity of grains and water. We took numerical values for the heat capacities of water and grains as Cvw=4156.7 J kg−1 K−1 and Cvgr=835 J kg−1 K−1, respectively. We then substituted these expressions into the expanded dimensional form of equation (2.1) and solved for the dimensionless ∂P′/∂t′.
During liquefaction, owing to the high bulk modulus of water, ϕ does not need to change very much in order to induce a unit of dimensionless pressure. We have introduced the variable ϕ′ which changes O(1) when ϕ changes by this tiny amount required for pressurization. On the other hand, represents a normalized value of ϕ as a fraction of its typical size, a value that can be considered very nearly constant in time during the event and that changes from place to place on this scale. The scaling of dimensional velocity vg was chosen so that 1 unit of vg′ represents a ground settlement of about 10 cm during a 30 s earthquake after computing Δϕ, as described below.
The remaining scale constants Δz, Δt, ΔS and Δϕ can be chosen arbitrarily, however the situation with the maximal interplay of different effects is the so-called distinguished limit , where all of the relevant terms in the equation have coefficient equal to 1 under typical conditions so that one unit of change in any variable or derivative induces the same size of effect compared with changes of any other variable or derivative. By calculating our scale factors to force the equation into this form, we identify the regime in which maximal interplay occurs. We can then examine the numerical size of the scale constants and determine whether the actual typical conditions are within this fully dynamic regime.
To calculate the scale constants, we evaluate the coefficients under the conditions kd′=1 and ϕ=ϕ0 as representative of typical conditions, and then require the important terms ∂2P′/∂z′2, ∂ϕ′/∂t′ and (∂P′/∂z′)(∂kd′/∂z′), which represent the effect of pressure diffusion, grain skeleton contraction, and spatial variability in permeability, respectively, to have coefficients equal to 1. In addition, we require the initial dimensionless gravitational static pressure gradient to be 1. These four equations are sufficient to constrain the values of Δz, Δt, ΔS and Δϕ, and we solve for them as follows: 2.8
Clearly, our assumption that diffusion, contraction and permeability variation are all important in realistic flow regimes is validated by the magnitude of these scale constants. The length scale Δz is O(1) relative to the size of deposits of interest. The value of Δt, which is a diffusion timescale related to (Δz)2/D, where D is the pressure diffusivity, is small compared with an earthquake duration of around 30 s, and even to a single cycle of loading of approximately 1 s. The heating scale ΔS remains to be considered below, and Δϕ is a reasonably achievable fluctuation in ϕ. The fact that Δt is small compared to the earthquake duration or the single-loading cycle means that a near-steady-state approximation of the flow should be valid for predicting the pressure field at the end of earthquake shaking. This is true even though the forcing may change significantly over a loading cycle as the diffusion term ∂2P′/∂z′2 implies a Gaussian weighted averaging over one dimensionless unit of length during one dimensionless unit of time. Therefore, all points in the deposit are influencing all other points and no localized pressure changes can last more than O(0.14) s without both rapid and sustained grain skeleton contraction in a single direction to balance the diffusion. We expect that drag forces on particles prevent the sustained contraction in a single direction from becoming too rapid (as will be shown in later examples), and because the typical earthquake loading fluctuates in both directions, the fast fluctuations in loading are averaged out and the pressure field remains near-steady state when time averaged over at least a cycle.
Note in the above scales we have left γS as a non-dimensional parameter which measures the fraction of the ‘distinguished limit’ level of heating which actually occurs in practice. If γS=1, then this implies a characteristic temperature change rate of 1.05 K s−1 which for a 30 s earthquake implies 32 K temperature change, a clearly untenable quantity: the sand would be scalding hot, a phenomenon not reported previously. The distinguished limit with heating is therefore more relevant for some mechanical apparatus in which external heating is applied, or perhaps a geothermal or fault gouge condition or within rapidly shearing bands where large quantities of energy are being dissipated in small spatial regions, such as in landslides or thin partially liquefied layers. However, as pointed out by  and verified in field observations by , dissipated mechanical energy is highly correlated with water pressure development. Clearly, the wave energy is used by the grains as activation energy for rearrangement and this process is irreversible and must produce some heat. To estimate the actual scale of heating within the soil during the early stages of typical liquefaction events, and thereby constrain γS, we made reference to the entropy generated by the flow of water through a characteristic soil under a characteristic pressure gradient. We calculate this entropy scale by assuming that the water flows horizontally for Δt seconds through a constant pressure gradient P0/Δz at a velocity given by Darcy's law, and the work done by the pressure field (the change in enthalpy) is completely converted to entropy at temperature T0, via drag from the particles 2.9
This is a lower bound for the relative entropy generation in the soil as it does not include any effect of the grain motion. To bound the range of reasonable values for γS, we assumed that during a 30-s earthquake soil temperature should not change more than ord(1) K as an order of magnitude more would not be ignorable and temperature must change at least ord(10−4) K owing to the flow entropy calculated above. As a representative value, we therefore chose γS≈5×10−3 under the assumption that the grain motion contributes significantly to entropy production but that the temperature change is still a fraction of a Kelvin. The implied characteristic change in temperature for the soil overall is 0.24 K in 30 s. This would seem to be a large fraction of the amount that could occur while still being neglected in most field or laboratory experiments. Even if the true value is an order of magnitude larger than this, our treatment of the effect as a small perturbation will remain valid.
After non-dimensionalizing, we replaced each of four important non-dimensional groups with a single non-dimensional symbol for each group as follows: 2.10 2.11 2.12 and 2.13Above, γvol expresses the compressibility of water in terms of the volumetric strain required to generate a characteristic unit of pressure, γμT represents fractional change in viscosity for a non-dimensional unit of temperature change, and Ggr is the relative density of grains to water, typically about 2.6–2.7; here, we have used 2.655. β is the local non-dimensional bulk modulus of the mixture. It gives the dimensionless pressure induced by a unit change of ϕ′ in the absence of flow, and is a function of the local value ranging in practice between about 1 and 2 for realistic values of . For important constants related to water, we used the values of density, heat capacity, viscosity and other important properties available online from the NIST database of the properties of water for a reference pressure of P=100 kPa and temperature T=293.15 K (20°C).
By construction, for typical conditions of soil liquefaction all the variables and their derivatives are O(1) so that the relative size of each term is determined from its coefficient. In certain conditions, these assumptions are violated, for example near-sudden transitions in permeability. To eliminate the negligible terms and identify the leading order terms of the perturbation we performed an asymptotic analysis of the equation detailed in the electronic supplementary material. Retained perturbation terms are relevant only near large gradients in permeability, large gradients in pressure or extremely low permeability regions. The final equation is 2.14
Equation (2.14) describes all relevant features of our model including small first-order perturbations owing to water compressibility and thermal effects. In this equation, is the non-dimensional relative heat capacity as a fraction of Cv0 which changes owing to spatial changes in porosity, and is the non-dimensional bulk modulus of the mixture defined above (equation (2.13)). In equation (2.14), the eliminated terms all had coefficients smaller than about 10−9.
It is at this point that the assumption of ‘fluid incompressibility,’ which is often invoked in soil mechanics, can be applied in a rational manner by taking a limit as . This assumption is a helpful one and justified for most purposes involving soil liquefaction. The primary conditions in which compressibility can become significant are high-pressure gradients in low permeability soils, such as at the centre of thin layers of silt, the spatial regions where these conditions occur are necessarily small, as discussed in the oscillating permeability example below. Thus the term multiplied by γvol is dropped, and the resulting simplified equation is 2.15
3. Qualitative results
Some basic simplified models for the qualitative behaviour expected from these equations are now considered. Although researchers such as Zienkiewicz et al.  have previously investigated the role of drainage and shown that drainage is usually important, these results have not generally been appreciated. Here, we examine several key experiments performed by others in the context of the predictions from our equations. First, recall that the characteristic time scale Δt≈0.14 s for our liquefiable sand system is smaller than the time scale for a typical earthquake loading period of about 1 s and much smaller than a typical earthquake duration of around 10–60 s. This means that the behaviour of the water pressure under realistic loading is always nearly in steady state. The diffusion term in the equation implies Gaussian averaging of the pressure field over the entire deposit in time Δt∼0.14 s so that any localized pressure feature will spread over the whole deposit in this short time via pressure diffusion induced by imbalances in flow described by the right-hand side of equation (2.1). Any ord(1) imbalance in any of the terms of equation (2.15) would induce ord(1) bars of pressure change over 0.14 s, returning the system rapidly to near-steady state. Even for fine sand where kd∼kd0/20 the time scale would be approximately the same order as a single-loading cycle and significantly shorter than typical earthquake duration. Only when permeability is similar to the permeability of silts kd0/1000 we can assume that the pressure diffusion time is long compared to earthquake duration. Under those conditions, the pressure changes are determined by the tiny volumetric changes caused by squeezing the nearly incompressible water, and the ‘undrained’ approximation might be reasonable. Even in regions of silty soils, as we shall see in examples, when a silt layer is small and in contact with a sandy layer, very high curvature of the pressure field and large gradients generated near these interfaces can cause the silt to reach steady-state pressures on the same time scale as the sand. Clearly, flow cannot be neglected in earthquake liquefaction of sand. In fact, for sands quite the opposite case occurs, the flow required to induce steady-state pressures happens extremely quickly.
To support this analysis with physical data, we can reference, for example, the centrifuge model B of  involving a fine sand overlain by silt. Shaking starts at approximately 2 s, and lasts until approximately 12 s. By 4 s, the pore pressure in the sand has come to a locally time-averaged steady state with a fast oscillation superimposed owing to continued shaking, and the silt has a similar if slightly slower response, achieving a locally time-averaged steady state at approximately 10 s. Both of these steady states are achieved prior to the end of shaking at 12 s (times given in prototype dimension). A steady-state solution of equation (2.15) based on time-averaged rate of change of porosity could be expected to be quite accurate at 10 s, even before the shaking stopped, if the loading were taken to be the time average of the actual loading over a few seconds.
Because of this near-steady state, we can approximate the pressure field at liquefaction by solving for the steady state of equation (2.15) numerically to find the average forcing ∂ϕ′/∂t′ required to make the water pressure field tangent to the total vertical stress. Once we have this solution we can find the location of tangency which is where the onset of liquefaction occurs. We use ∂ϕ′/∂t′ constant in space except where noted below. As heating was earlier determined to be a small perturbation for sands, we ignore it here for simplicity except in the example involving thermal liquefaction in silts. Also, we set and ignore the small effect of spatial variability in ϕ except as it influences kd which is spatially varying as shown in each figure. By adjusting the characteristics of the soil, especially the permeability field and the spatial variation in the forcing, we can determine how the permeability and permeability gradient affect the onset of liquefaction. The model is only valid during the onset of liquefaction. After the soil has fully liquefied, we may have significant changes in the fluid dynamics: for example, Darcy's law with a constant permeability will not hold in a fluidized bed. Several authors have pointed out that the permeability required to satisfy Darcy's law changes dramatically when effective stresses drop near zero [17,18].
Example 3.1 (Uniform Loose Hydraulic Fill)
A model of a uniform loose hydraulic fill, such as those found at reclaimed land in bays and port facilities worldwide would be a uniform deposit with kd′=1 of depth z′≈1. We assume for simplicity that the water table is at the surface of the ground, and we apply a 0.12 non-dimensional pressure unit load on the ground designed to model 1/2 m of concrete pad to ensure that the effective stress is not zero at the surface (figure 1). In this situation, owing to the uniform permeability, the term in equation (2.15) involving the spatial derivative of kd disappears. If we assume the contractive forcing ∂ϕ′/∂t′ is constant in space and we are solving for steady-state conditions, then the curvature of the pressure field is constant and we have a parabolic pressure curve.
In figure 1, we show the results of this simulation and we find that the soil just initiates liquefaction via the loss of grain contacts at z′≈0.61 when ∂ϕ′/∂t′=−1.53 corresponding to a ground surface settlement velocity of (∂ϕ′/∂t′)ΔϕΔz/Δt= −0.08 cm s−1. Liquefaction occurs at this location not because of special properties of the sand at this point as all sand is treated equally, but rather because of the flow and its interaction with the boundary conditions as will be explained later. As our model is of the initiation of liquefaction, we can assume further deformations occur after initiation. For an event that takes 3 s to initiate over a 10 m deposit, this implies a vertical settlement of about 2.4 cm in the first 3 s. For a 30 s earthquake with continued deformation at similar rates, total settlement would be of the order of 24 cm. Normally, if the ground drops several centimeters we would see serious consequences for pipelines or structures on spread footings.
Example 3.2 (Fine-scale variability in permeability around an average value)
Owing to the importance of permeability variation, we consider a model of a loose hydraulic fill sand whose layering structure causes a fine-scale oscillation in permeability of 30% about a typical value . This oscillation occurs on a small length scale O(1/40) of the total length, and is a simple model for variability in hydraulic fill material and the natural segregation of sandy soil by varying sinking rates owing to varying grain sizes, as pointed out by . Note that the water pressure curve for this case is very similar to that of figure 1 (see the electronic supplementary material).
This example allows us to highlight the role of the variability in permeability. The following argument shows that the diffusive term dominates over the permeability gradient-related terms except in regions where the permeability is quite small which must necessarily be thin regions.
If the permeability field is O(1) throughout, but has a fine-scale oscillation of dimensionless length scale λ≪1, then the term ∂kd′/∂z′ becomes O(1/λ), which can be quite large. To balance this in the steady state, we can either generate curvature of the pressure field to induce diffusion, or if diffusion is quenched by kd(∂2P′/∂z′2)∼0 then (∂P′/∂z′+1) must be O(λ) to keep (∂kd′/∂z′)(∂P′/∂z′+1)= O(1) implying a very small flow rate.
Suppose that a fine-scale variation in P′ develops so that P′=F(z′/λ)+G(z′) where F is an O(1) function representing the fine-scale variation and G is O(1) representing the large-scale variation. Then, |∂2P′/∂z′2|=O(1/λ2) with a factor 1/λ coming from each derivative of F with respect to z′. In order for diffusion to be quenched, we must have kd′≪λ2 so that kd′(∂2P′/∂z′2)≈0 in equation (2.15). As in real materials 0<kd′<O(1), necessarily the region where both kd′≪λ2 and |∂kd′/∂z′|=O(1/λ) is extremely small. To see this, consider that if kd′ is small, and the derivative of permeability is negative, then it must become non-negative quickly to keep kd′>0; or if the derivative is positive, kd′ increases rapidly, so it does not remain small for long.
If the overall size of kd′ is small everywhere, then we can rescale the kd′ variable everywhere, and by rescaling the t′ variable the equation remains symmetric though now all processes take longer in dimensional time. This shows that for lower permeability soils an overall lower level of forcing is required for liquefaction as the same change in ϕ′ would liquefy the soil if imposed over a longer dimensional time. In general, the conclusion is that large stable gradients in pressure occur primarily in small regions of very low permeability soils surrounded by higher permeability ones as shown in a later example with a thin silty layer.
By using a Fourier series approach to approximate the analytical solution we find that the oscillations in permeability are completely dominated by the tendency for pressure diffusion in the steady-state case so that the pressure field is nearly indistinguishable from the earlier constant permeability field (figure 1), and does not develop an important fine-scale oscillation as the permeability never becomes small enough; in this example, ∂ϕ′/∂t′=−1.45 and z′=0.599 compared with the constant permeability model where ∂ϕ′/∂t′=−1.53 and z≈0.61. As the figure for this example is very similar to figure 1, it is reproduced in the electronic supplementary material.
Example 3.3 (Silty Inter-Layer)
Another case of interest is a simple model for the experiment of . In this situation, we have loose sandy soil with kd′=1 and a thin layer of much lower permeability kd′∼0.001 with a large gradient in permeability around this minimum (figure 2). Because of the large permeability gradients, there can be fast-changing boundary layers in which the pressure gradients are high and develop over small regions so that the curvature of pressure is quite large. If we solve the steady-state equation we find that with ∂ϕ′/∂t′=−1.23 the confining of pressure by the inter-layer forces a localized loss of grain contacts at z′=0.497 just below the centre of the layer. In the experiment of  there is the formation of a thin water film below the silt inter-layer in a manner consistent with these results.
A follow up to the above experiment considered a silt inter-layer in a fine sand, in which a long-lasting stable layer of water below the silt layer was observed . The water layer will be long lasting if the upper soil cannot fall quickly through the water owing to drag. This can be explained in general by the ratio of drag force to weight for small particles. If we consider a single spherical soil particle of radius R, the well-known formula for Stokes drag is fd=−6πμRv , whereas the weight for a spherical grain is w=(4/3)πR3ρg. If we take the fluid velocity in the fluid layer to be the apparent velocity in the porous medium immediately below the layer v=(kd/μ)(∂P/∂z) the characteristic value of fd/w for one non-dimensional unit of each variable is therefore 9kd0/2GgrR2. For sand grains with R∼0.5 mm this ratio is 0.001, whereas for a silt grain with R∼0.01 mm the ratio is 2.5 indicating that silt particles which fall down into the water layer are easily pushed back up by the drag caused by fluid flow, whereas the sand will easily sink under these conditions. In addition to this basic interpretation,  and  show that the drag on the particle increases dramatically as it integrates into the porous medium at the water-silt boundary so that once a layer is formed, the particles above are unlikely to detach and sink. The effect of grain sizes is discussed more comprehensively in an additional section below.
In another experiment by Kokusho & Kojima , which involved a fine sand overlying a coarse sand, a transitory turbulent water layer was observed at the interface between the sands. The differences between this coarse- and fine-sand case compared with the fine sand and silt case can be understood in terms of two factors: the lack of a long-term strong gradient in pressure which would induce a large flow velocity and the larger grain sizes involved in the coarse- and fine-sand experiment. In further cases, involving a fine sand inter-layer in a coarse sand, the results were substantially similar to the sand with silt inter-layer experiment with an altered timescale owing to the different permeabilities involved . Finally, the authors alter their fine sand overlying a coarse-sand experiment where only transient turbulent water films formed at the interface by reducing the height of the water table. In this case, the water layer at the interface is stabilized . The lack of water above the water table means that the effective size of the upper fine sand layer is only 7 cm compared to 90 cm in the original experiment, because the unsaturated sand above the water level does not affect the flow and the fluid pressure drops to atmospheric at the fluid/air boundary. This effectively thinner upper soil layer means a larger pressure gradient can form and higher fluid velocities develop, thus higher drag develops on the upper sand particles which stabilizes the water film.
Example 3.4 (Liquefiable Sand Inter-Layer)
A typical warning sign for geotechnical engineers is the existence of a loose sand layer between layers of denser and especially lower permeability material. This is the inverse of the silty inter-layer example given above. When we make permeability kd′=1/8 outside of a region of width about 20 cm (10 times as wide as the silty inter-layer) where permeability is 1 we find that the contraction of the deposit pushes water upwards out of the coarse sandy layer and the deposit liquefies in the fine sand above at z′=0.65. Also ∂ϕ′/∂t′=−0.198 which is much lower loading than in either the uniform or silty layer case. In this case, the location and manner of liquefaction is less important than the magnitude of the loading required (figure reproduced in the electronic supplementary material for completeness). Owing to the lower permeability of the main deposit, the diffusion is less dominant, and with less drainage, it is easier to cause liquefaction. However, since the gradients in pressure remain O(1) the perturbation terms do not become important here. This highlights the phenomenon that ‘drainage’ is a non-local phenomenon which is affected by the entire permeability field. Sand itself may be locally well drained, but can be confined by low permeability boundaries.
Example 3.5 (Sand overlain by silt)
This case (figure 3) is modelled after  in which a silt overlies a loose sand and is subjected to shaking in a geotechnical centrifuge. In these experiments, pore water pressure was recorded in several locations within the sand and silt. For our model of this case, we have rescaled the z-axis to correspond to the depth of the experiment and correspondingly rescaled P′, ϕ′ and t′. Contractive forcing ∂ϕ′/∂t′ is spatially varying and appreciable only in the lower sand region; the silt region, whose permeability is ϕ′=0.001 is assumed not to contract. Clearly, in figure 3 water is squeezed out of the lower sand and pressure diffuses into the upper silt, liquefying it at the base. The equilibriation time in the silt is much longer, and only the sand region and transition region would be likely to equilibriate during a typical earthquake duration. The experimental results of  qualitatively agree with this steady-state solution. Note that the consequences of liquefaction, a failure via a sand-boil, were observed to occur much later, presumably near a localized weakness within the silt layer. In our example, liquefaction occurs at z′=0.695, which is slightly above the centre of the permeability transition, which is located at z=0.629. The required loading is ∂ϕ′/∂t′=−0.011, which ultimately is not rescaled (both ϕ′ and t′ are rescaled by the same factor, hence cancelling). This is significantly smaller than in the uniform sand example (figure 1), where ∂ϕ′/∂t′=−1.53. The relatively thick silty layer makes drainage difficult, as expected and this contributes to the observed reduction in the necessary loading.
Example 3.6 (Theoretical process of thermal liquefaction of silt)
By substituting in the incompressible equation (2.15), we can consider the case when typical permeability is small relative to sand, such as ϵ<0.001. When this is the case, diffusion does not act on the same time scale, and in fact the diffusion time scale Δt/ϵ>140 s is larger than the typical duration of the earthquake. This rescaling process also makes the γvol term more dominant, and the full compressibility effects might become important though we ignore them in this example as our purpose is simply to highlight a new phenomenon, which could provide an avenue for further study. If a silt is somehow precluded from sustained grain skeleton contraction ∂ϕ′/∂t′, perhaps owing to its initial porosity being near the equilibrium state, it could still theoretically liquefy by absorption of wave energy into thermal energy. The simplest case is a single homogeneous layer of silt similar to the homogeneous layer of sand above. Instead of contractive forcing, we apply ∂S′/∂t′ and attempt to liquefy the silt via heating. The equation is similar to uniform contraction, and produces the same pressure field as figure 1. We find that the homogenous layer liquefies at z′≈0.61 when ∂S′/∂t′=1.53ϵ/γS. The question arises then as to whether ϵ/γS=ord(1) in some regimes of shaking and grain size implying that silt might thermally liquefy by absorbing wave energy to heat. This is an interesting question for which a simple centrifuge experiment using a well-insulated bucket and measuring before and after temperatures could give insight, or perhaps borehole measurements of thermal changes during an earthquake might be available to provide insight.
4. The role of grain sizes and drag forces for liquefaction consequences
As we are considering the liquefaction of silt, it should be mentioned that pure silts and gravels are generally thought to be lower risk for liquefaction than sands, though a more nuanced view of silt and gravel liquefaction is developing (cf. [21,22]). Examining the cases for different grain sizes within the context of simple steady-state water flow described here can be instructive.
In order for sustained contraction of the grain skeleton to occur we need gravitational and overburden forces on the particles to be balanced by drag forces of the water so that the grains travel downwards at a slow terminal velocity while the water is pressurized and travels upwards towards the surface. Because of this balance of forces, and the rapidity of pressure diffusion, the details of grain motion at each point in the loading cycle are less important than the overall average contractive trend.
In the examples by , the sand particles descend rapidly through the water while the silt particles sit on top of the water layer. These settlement velocities are related to the drag forces on the particles. However, when all particles are silt sized, even if all particle contacts are lost, settlement will be extremely slow and particles will not move far from their initial positions during a loading cycle owing to high drag forces relative to gravitational settlement forces. Consequences of such an event on level ground may be significantly less destructive than that in sand where drag forces do not affect mobility as strongly and settlement is more rapid. For gravels, mobility is even more rapid than sand, but contraction of the gravel does not induce strong pressurization because the pressure dissipates extremely rapidly owing to the high permeability, unless surrounding lower permeability soils impede flow. This implies that we should expect fast settlements but only transient liquefaction in contractive gravels.
Although in the vertical direction we expect rapid settlement from gravel particles and slow settlement from silt, in terms of shear displacements the roles may be reversed. The rapidity with which gravel particles descend may allow them to regain contact more rapidly thereby reducing the time during which shear strength is lost, whereas silt particles will take a long time to regain their shear-transferring ability through settlement. Soils with a wide range of particle sizes will be more complex. The development of a soil behaviour model which takes the fluid drag forces into account explicitly is a worthy goal. Whatever the case, this analysis highlights the gradient in particle size, which is also related to the gradient in permeability, as playing an important role in the liquefaction process especially with regard to the development of fluid lubrication layers, as seen in [6,10].
5. The problems with small-scale empirical tests
To understand the role of water flow on small-scale samples, we can consider our equation rescaled so that the typical length scale is similar to that of a tabletop triaxial sample. The resulting analysis will show that water pressure within a triaxial sample is at all points and all times in equilibrium with the pressure applied to the water at the boundary of the sample via the flexible membrane. Hence, results of water pressure measurements from triaxial samples are exclusively informative about how the rubber membrane with a constant external pressure bath stretches and interacts with the grains. Furthermore, by analysis of one of our simplest examples (example 3.1), we will show that the location of initiation of liquefaction in that example is entirely determined by the diffusion and flow of water towards the drained surface boundary conditions. As the rate of pore space collapse ∂ϕ′/∂t′ was spatially constant in the example, the onset of liquefaction was not determined by any special soil-mechanical characteristics of the soil at that point.
To begin, consider example 3.1 in figure 1 where all the sand is uniformly similar in initial porosity and permeability conditions. If we take a slice of dimensionless length 0.01 from this example, it corresponds to a sample similar in size to a hypothetical triaxial sample of about 10 cm characteristic length. Recall that all points in the model of example 3.1 were squeezed equally as ∂ϕ′/∂t′ was a constant in space. The rate of squeezing in example 3.1 was 1.53 units of Δϕ per 1 unit of dimensionless time, or a volumetric strain of about 1% in 130 s, where 130 s is a reasonable timescale for the duration of a tabletop undrained liquefaction experiment. In the absence of drainage, at this volumetric strain rate, we would generate in our hypothetical sample 100 kPa of additional pressure every 0.6 s, and at the end of 130 s we would have about 22.5 MPa of water pressure. Instead, in example 3.1, we have about 100 kPa of excess pressure (1 dimensionless unit over atmospheric) at infinite time, as shown in figure 1 at z=0.6 where liquefaction occurs. This discrepancy is owing to drainage in steady-state conditions. Note that the pressure field in example 3.1 is influenced by three factors: the rate of squeezing (∂ϕ′/∂t′=−1.53), the constant permeability, and the boundary conditions (impermeable at z=0, constant atmospheric pressure at z=1 and zero lateral flow). This is in contrast to the triaxial test where drainage is prevented (impermeable at top and bottom), and the lateral boundary conditions are much closer, and more complex, owing to the surrounding rubber membrane and constant pressure bath.
In example 3.1, the reason the point at z=0.6 has liquefied is not that it has been squeezed harder or that its grain skeleton was more susceptible than other points, it is that the fluid is flowing continuously through the sample, and owing to a combination of squeezing, and imbalance in fluid flow, the density of water at this point is just slightly higher than it was before we began squeezing the sample. The role of water flow is highly non-trivial even in this extremely simple example.
By contrast, consider taking our slice of dimensionless length 0.01, removing it from the ground, and placing it into a triaxial machine. Just as in the full-scale equation, the water pressure in our impermeable-membrane-wrapped triaxial sample is dependent on three things: the dynamics of the loading, the uniform permeability which remains unchanged and the pressure induced by the boundary conditions.
Now, the boundary conditions have been brought 100 times closer to the sample (about 10 cm instead of 10 m) and have been changed. Although a triaxial test is not really a one-dimensional system, we can begin to qualitatively analyse the effect of investigating the small-scale samples using the same equation (2.15) by choosing new scale constants to interpret the results.
We leave the characteristic pressure at P0=100 KPa because this is the order of magnitude of pressure we must induce to liquefy the sample. As Δϕ depends only on the pressure scale, it remains unchanged. We ignore heating for these purposes. We are left with re-evaluating Δz and Δt. The characteristic size of our sample is Δz=0.1 m, and the new value of Δt is determined from Δz and equal to Δt=(Δz)2ϕ0μ/kd0KB≈1.3×10−5 s or 13 μs. Next we consider the loading, ∂ϕ′/∂t′. It is safe to say that this is nearly zero on our new scales, as the loading is such that the order of magnitude of pressure changes that these experiments induce is enormously less than P0=100 kPa in 13 μs.
The gradients in pressure are tiny as we will never induce anywhere near 100 kPa fluid pressure differences across the 10 cm sample; in fact, gradients will be only of the order of 1 kPa across the vertical distance of the sample owing to gravity. In the absence of any significant loading or flow, the only remaining factor that influences the water pressure is the equilibrium with the boundary conditions. It is no wonder that experimentalists have found an interest in accounting for the ‘membrane penetration effect’ (e.g. ), as this is the only effect at work on the water pressure in these experiments. Thus, we should be highly skeptical of any results from tabletop triaxial or similar specimens insofar as they are supposed to represent how the water pressure in the soil behaves in the ground. As we have shown at the beginning of this section, under the ground, water flow is extremely important and tabletop experiments, as currently performed, simply cannot model this process.
By deriving a comprehensive equation for the water pressure development within layered soils and analysing some critical examples we have shown that fluid flow driven by porosity and permeability changes is the main factor that causes pore pressure change during liquefaction. By comparison with experimental results [5,6,10], we show that sandy soil is well drained enough that this soil comes to steady state on the same time scale as the earthquake. Also, owing to our model of this drainage, we can identify several modes of liquefaction, from fluid-flow-induced ‘fast settlement’ to localized water film formation at permeability boundaries, to flow-induced liquefaction of low permeability layers that overly a high permeability layer, and even potentially thermally induced liquefaction when permeability is very low. In each of these example cases, we do not attempt to couple the behaviour of the grains to that of the water, or provide time-varying dynamic solutions, as no large-scale continuum model is available for grain–fluid interaction behaviour comparable in quality to the DEM model of . Our intention in this paper has been to highlight the importance of fluid flow and spur an understanding of the need to develop soil models that can account for flow. The need to account for flow has also been determined based on empirical grounds by authors, for example .
Modelling a 1/4 m square cross section 10 m long with 1 mm diameter typical grain sizes would require perhaps 108 grains and be computationally prohibitive. However, for qualitative understanding of the liquefaction process, it is sufficient to perform a careful analytical analysis of a continuum model of water flow through soil driven by changes in pore volume. Our analysis provides tremendous insight without supercomputing requirements, from which we make several important conclusions:
— Soil liquefaction in water saturated sands is a process usually involving fluid flow owing to gravitational pumping, not purely fluid stagnation in an ‘undrained’ event. Thin regions of low permeability silt do not prevent drainage, as large pressure gradients will develop to force fluid through these layers. However, their presence significantly influences the location and ease of onset of liquefaction. Large regions of low permeability soil can significantly reduce drainage, leading to liquefaction at significantly lower rates of grain contraction.
— Owing to water pressure diffusion, we cannot examine the liquefaction process by simply considering local properties of soil in a small sample on a tabletop apparatus, given the characteristic length scales of ord(10) metres and corresponding characteristic time scales of tenths of seconds. In a tabletop apparatus with samples on the order of 10 cm the timescale for diffusion of pressure to the impermeable rubber boundary and hence to the constant pressure boundary conditions is on the order of 10 μs so that the pressure throughout the sample is equal to the boundary pressure at the membrane. Therefore, geotechnical centrifuge experiments are the only meaningful physical experiments for this phenomenon.
— In lower permeability materials, such as thick regions of fine sand or silt, the time scales for the water flow process can be longer. These long characteristic times occur when the typical permeability throughout the material is O(10−13) m2, which is well below the sandy range . Owing to grain size effects, fluid flow in silt causes high drag on the particles, limiting their mobility. Liquefaction in extensive silt deposits is an interesting and different flow regime.
— Impedance to flow caused by local sharp reductions in permeability is the main driving force in localized liquefaction events and can cause the formation of water films that lubricate layers of soil allowing large lateral spreading. For level ground, failure is expected via localized sand boils at fissures or the liquefaction of overlying layers as fluid flow stagnates. This phenomenon can only be explained by explicitly acknowledging the flow of water.
— ‘Liquefaction’ as the term is commonly used is not necessarily associated with the total loss of effective stress, as is already known. Large settlements could occur without the formation of a fully liquefied layer if loose soil settles while water flows to the surface. This may especially occur in contractive gravel or coarse sand fills.
— In silts, where permeability is small and diffusion cannot act on the earthquake duration timescale, extensive shaking might concievably cause liquefaction through heating effects, a previously unexplored phenomenon.
↵1 Note on our use of asymptotic notation: x=O(y) means that |x/y|<C for some positive constant C, which is expected to be not extremely large, x=ord(y) means both |x/y| and |y/x| should be treated as close to 1 and x=o(y) means |x/y| is negligibly small.
- Received July 10, 2013.
- Accepted January 23, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.