Drainage explains soil liquefaction beyond the earthquake near-field

Domain: earth_science
Published: 2023-09-27
Editorial outcome: major_revision

Abstract
Earthquake-induced soil liquefaction is a devastating phenomenon associated with loss of soil rigidity due to seismic shaking, resulting in catastrophic liquid-like soil deformation. Traditionally, liquefaction is viewed as an effectively undrained process. However, since undrained liquefaction only initiates under high energy density, most earthquake liquefaction events remain unexplained, since they initiate far from the earthquake epicenter, under low energy density. Here we show that liquefaction can occur under drained conditions at remarkably low seismic-energy density, offering a general explanation for earthquake far-field liquefaction. Drained conditions promote interstitial fluid flow across the soil during earthquakes, leading to excess pore pressure gradients and loss of soil strength. Drained liquefaction is triggered rapidly and controlled by a propagating compaction front, whose velocity depends on the seismic-energy injection rate. Our findings highlight the importance of considering soil liquefaction under a spectrum of drainage conditions, with critical implications for liquefaction potential assessments and hazards.

Introduction
Seismically induced soil liquefaction is a natural hazard that commonly occurs during earthquakes1. During liquefaction, a soil that initially possessed an elasto-plastic rheology and was capable of supporting the load of infrastructure, loses its strength and stiffness in response to earthquake shaking, consequently exhibiting fluid-like rheology. Earthquake-induced soil liquefaction results in buildings and infrastructures sinking2, floating and tilting2, ground lateral spreading2, settlement3, and landsliding4. Liquefaction damage often leads to extensive human casualties4,5, destruction of lifelines2,6, and economic losses6–9, that may result in complete abandonment of formerly inhabited areas7, posing a significant challenge to community resilience10.

The classical mechanism explaining seismically induced soil liquefaction1 considers the soil as an effectively undrained medium. Upon cyclic shear, an initially loosely-packed soil tends to reduce its pore volume, as readily documented under dry and drained conditions11. If the pore fluid flow rate is slow compared to the rate of porosity reduction, as is expected in an undrained soil response, the pore fluid is trapped within the contracting pores and its pressure increases. If the pore pressure builds up to the level of the overburden stress (commonly lithostatic values), then the effective stress reduces to zero12, the soil loses its shear strength and stiffness and is said to be liquefied8,11. Undrained lab experiments8,11,13,14 showed that during continuous shaking, and depending on the initial soil density and the applied shear stress, the pore pressure builds up gradually and reaches lithostatic values after several to tens of shear cycles.

Despite the overall success of the undrained perspective in describing the conditions leading to pore pressure rise and soil strength and stiffness loss during earthquakes, it struggles to explain soil liquefaction beyond the near-field, far from the earthquake’s epicenter, where the seismic energy density input is small. Empirical inferences established a lower bound of 30 J m−3 for the seismic energy density required to induce liquefaction by undrained consolidation15,16. Consequently, as the seismic energy decays away from the earthquake’s epicenter, liquefaction events beyond the near-field should become uncommon. Nonetheless, the majority of the events in an extensive soil liquefaction compilation16,17 were triggered beyond the earthquake near field, at a distance greater than one fault rupture length, where the seismic energy density is well below 30 J m−3, and as low as 0.1 J m−3 (Fig. 1). The discrepancy between the leading theory and field observations of soil liquefaction indicates that our understanding of the conditions and processes associated with earthquake-induced soil liquefaction is incomplete.

The main attempt to reconcile theory and observations16 invoked seismically induced permeability enhancement7,18 between deep pressurized aquifers and the liquefied layer. However, dynamic permeability increase and the availability of buried high fluid-pressure sources represent unique geometric and hydrologic conditions, likely precluding it from being a general mechanism for liquefaction beyond the near field.

In a rare movie capturing soil liquefaction19 at the Makuhari Seaside Park in Chiba (Japan), during the Tohoku earthquake (2011)20, the photographer commented while pointing to the lawn: “... there was this water just coming up right there, on the ground... and the ground is just swaying right now”. This testimony suggests that fluid drainage toward the surface during an earthquake could play an important role in the process of soil liquefaction, which agrees with recent theoretical and experimental studies proposing that rapid fluid flow could be instrumental in initiating liquefaction21–28. This implies that the undrained liquefaction initiation mechanism does not necessarily cover the full spectrum of conditions leading to soil liquefaction.

Drained liquefaction initiation implies that the timescale of fluid flow is shorter than the timescale associated with earthquake-induced soil deformation. In this scenario, porous fluid flow toward a drained boundary is accompanied by pore pressure gradients that exert seepage forces on the soil grains, supporting their weight, weakening grain contacts, and reducing soil strength. The notion that pressure gradients and seepage forces could fully support grains is, in itself, not novel. Static vertical pressure gradients supporting a layer of grains is known as quicksand conditions29,30. Similarly, interstitial fluid ejection leading to ground settlement is a known post-liquefaction failure mechanism31–33.

The evolution of the pore pressure in a deformable saturated granular media can be described by a diffusion equation with a source term related to the granular skeleton deformation21,22,27,34:

∂P0 / ∂t = ∇ · [κ∇P0] + 1/β ηϕ ∇ · u_s = 0 (1)

where P0 is the dynamic pore pressure deviation from hydrostatic value (P0 = P − Phyd), β and η are the fluid compressibility and viscosity, respectively, κ is the permeability, t is time and ∇ is a spatial derivative. The second term in Eq. (1) is a diffusion term arising from Darcy flow, while the third term describes the internal source for dynamic pore pressure, due to divergence of solid grain velocities (u_s). This term can be approximated22,27 as the rate of pore space compaction and dilation, ∇ · u_s ≈ 1/(1−ϕ) ∂ϕ/∂t, where ϕ is the porosity.

Since fluid drainage within the soil granular media is expected to obey Darcy’s flux law, the characteristic velocity scale in Eq. (1) is identified with u_0 = (κ_0 / ϕη) · (σ_h / h), where σ_h / h describes the initial effective lithostatic stress gradient, which is also the pressure gradient during liquefaction. σ_h is the initial effective lithostatic stress at depth h, and κ_0 is the characteristic permeability. There are two length scales characterizing the system: Stress and pressure change gradually over the layer depth, h, yet grain divergence and convergence could occur over a different length scale, l, which could be as small as several grains wide21,27. Subsequently, non-dimensional parameters (represented by caret symbols, ^) can be defined as follows: ∇ = ∇^ / l where ∇ appears as a gradient operator, u^ = u / u_0, t^ = t / t_0, P^ = P_0 / σ_h, and κ^ = κ / κ_0. Eq. (1) can then be re-written as:

∂P^ / ∂t^ = [l^2 / (t_0 D)] ∇^ · [κ^∇^ P^] + [ϕ_0 l^2 / (t_0 D)] ∇^ · u^ = 0 (2)

The non-dimensional coefficient in front of the first term in Eq. (2) is known as the Deborah number (De)21,22,27,35,36:

De = t_d / t_T = β_f ηϕ_0 l^2 / t_T κ_0, (3)

which expresses the ratio between the timescale for pressure diffusion t_d = l^2 / D, where D = β_f κ_0 / ηϕ_0 is the pore pressure diffusion coefficient, and the timescale of deformation imposed by the shaking period, t_T = T. The De number provides a metric for evaluating the system’s drainage conditions. When De ≪ 1 and t_d ≪ T, pore pressure diffusion is sufficiently rapid, so as to allow drainage during shaking. In this case, which we term “drained”, the first term in Eq. (2) becomes negligible, and the diffusion (second term) balances the source term (third term). When De ≫ 1, the layer is “undrained”, the diffusion term is negligible, and the source term is balanced by the temporal derivative of the dynamic pressure21. Notably, contrary to previous studies that used the term ’drained’ to describe an end-member with no change in fluid pressure8,37, in the current study, ’drained’ implies De ≪ 1, and pore pressure gradients could emerge in response to skeleton deformation. For a discussion on drainage-related terminology and its relation to Eq. (2), see Supplementary Note 1.

To apply this general formulation to the deformation of a shallow soil column, we consider a layer of saturated cohesionless grains with a free and drained surface (where the pressure is maintained at a constant value, P = 0) and where no internal permeability barriers are present (Fig. 1). When such a soil layer is relatively loosely packed, with initial porosity ϕ_0, and when it is subjected to horizontal shaking, an upward propagating compaction front (also referred to as a “solidification front”) develops27,31,38,39. The propagating compaction front separates two regions within the layer27 (Fig. 1a, b): a lower region, in which grains have compacted to porosity ϕ_c < ϕ_0 and are approximately stationary in the vertical direction. In this lower region, the pore pressure gradient is nearly hydrostatic, although the pore pressure itself is elevated to the value of the pressure at the front (Fig. 1c). In the region above the front, grains continuously settle at a uniform velocity while maintaining their initial porosity, ϕ_0. The settling grains exchange place with the upward flowing pore fluid. The pore fluid pressure gradient above the front can, and normally will, become as high as lithostatic (Fig. 1c and Supplementary Note 2). This fluid pressure gradient is the source of upward directed seepage forces that support the settling grains. Under these conditions, the upward front velocity, u_front, is theoretically predicted to be27:

u_front = ϕ_0 / (ϕ_0 − ϕ_c) · (κ_0 / η) · (dσ_0 / dz), (4)

where dσ_0 / dz is the effective static normal stress gradient (Fig. 1). The uniform porosities above (ϕ_0) and below (ϕ_c) the front imply that active compaction occurs only across the relatively narrow front, which acts as an in-situ, migrating, pressure source, continuously forcing upward fluid drainage27 (Fig. 1).

Here we show that the drained liquefaction mechanism, with its associated upward migrating compaction front, predicts liquefaction events triggered under low seismic-energy density. Thereby, the drained mechanism provides a general model for explaining the vast number of liquefaction events triggered beyond the earthquake near-field. This, in turn, has critical implications for the physics of liquefaction, the conditions for liquefaction triggering, and consequently, the evaluation of liquefaction potential and associated hazards.

Results
We performed grain-scale simulations and experiments of horizontally shaken layers of water-saturated cohesionless grains with a free surface. The simulations used a coupled Discrete Element Method (DEM)-Computational Fluid Dynamics (CFD) model40, and the experiments were conducted in a transparent box, allowing inferences of grain motion and measurements of pore pressure by using an array of pressure transducers (see Methods). The simulations and experiments show the dynamics predicted by the compaction front model.

Evaluating the drainage conditions
The drainage conditions in the simulations and experiments were evaluated by estimating the De number following equation (3). The length-scale controlling pressure gradients (and thus fluid fluxes) is conservatively chosen here as h = H, where H is the layer height41. In situations of homogeneous compaction and dilation41, H would also control the divergence of grain motion. However, when a compaction front is present, grain compaction and dilation are localized36 at the front. The natural length scale that thus emerges for the divergence of grain motion is l = w < H, where w (Fig. 1) is the width of the compaction front. The simulations and experiments show that w spans several tens of grain diameters (≈20). As a consequence, the maximal value of De is ~10−2 in the simulations and ~10−4 in the experiments (see Table 1 for simulation and experiment parameters). This analysis indicates that the behavior we observed in the experiments and simulations arises from drained layer dynamics.

Liquefaction indicators in drained layers
Simulations and experiments determined to be controlled by drained dynamics show four indicators that are widely associated with soil liquefaction in the field and the lab: pore pressure rise, soil settlement, attenuation of shear waves, and degradation of shear modulus.

The dynamic pore pressure rises quickly in response to the onset of horizontal shaking and reaches approximately the value of the initial effective vertical solid stress (Fig. 2a). The duration at which the pore pressure remains elevated is a function of depth27 and is set by the compaction front arrival. Once the front passes a certain depth, the pore pressure starts to decrease, so the closer a point is to the surface, the longer the pressure remain elevated at that point. The event ends at a time, t_e, which corresponds to the time it takes to initiate liquefaction, t_i, plus the time it takes the compaction front to propagate a distance L ≤ H, from its initiation depth (Supplementary Fig. 1) to the surface, t_e = t_i + L / u_front. The initiation time, t_i, is found here to be exceedingly short, with a conservative median value of 0.25 s in simulations and 2.5 s in experiments (see also Section “Evaluating the compaction front velocity (u_front), the duration of liquefaction event (t_e) and the surface settlement (ΔH)”).

Concurrently with the pore pressure rise, the excited soil layer compacts continuously and linearly (Fig. 2b). Despite continued shaking, we find that the soil stops settling and reaches a new equilibrium configuration after time t_e.

Grains are shaken horizontally by the shear waves propagating from the excited layer base. Figure 2c presents the mean horizontal grain velocity time series at different depths. Shortly after the onset of shaking, shear wave amplitudes become strongly attenuated throughout the layer, as expected from a fluid-like medium. At any given depth, attenuation persists until the compaction front arrives, after which shear-waves resume the amplitude of the input shear. The black lines in Fig. 2c follow the positions where the velocity amplitude increases back to >60% of the shaking velocity imposed at the bottom boundary, chosen here to depict the front position. The observed trend indicates that, similar to the pore pressure dynamics, attenuation lasts longer closer to the surface and overall continues up to ~t_e.

Figure 2d presents the relation between the shear stress and the shear strain at the mid-depth of a simulation layer. The mean slope of the stress-strain curve, known as the shear modulus, is used as a metric for the shear strength of a material42,43. We observe that the stress-strain curve flattens soon after the application of shaking, over less than two shear cycles, indicating that the saturated soil layer has dynamically lost its shear strength. The soil progressively regains its strength as the front progresses upwards, displaying a finite stress-strain slope.

Drained liquefaction beyond the near field in simulations and experiments
The simulations and experiments were forced with a range of shaking amplitudes (A) and angular frequencies (ω) (Tables 2 and 3), leading to an energy density range of 0.07–7 J m−3. The average seismic energy density in one shear cycle is calculated as16,44 e = (ρ_s / 4) PGV^2, where PGV = Aω is the amplitude of the imposed harmonic cyclic velocity. Thus, the four liquefaction indicators described above emerged although the input energy density corresponded to low, far-field, values and was smaller than the previously established liquefaction triggering threshold of 30 J m−3 15,16.

Analysis of simulation results further shows that the change in porosity across the compaction front, Δϕ = ϕ_0 − ϕ_c, correlates with the applied seismic energy density (Fig. 3a), and has an even better correlation with the rate of the seismic energy density input (Fig. 3b), which can be evaluated as the seismic energy density over one period of shaking, e / T (or as the seismic power proportional to PGV · PGA in monofrequency harmonic oscillations).

Consistent with the prediction of Eq. (4), we further find that a larger porosity change is associated with a slower propagating compaction front (Fig. 3c) and a longer liquefaction event, t_e ∝ 1 / u_front.

Combining the above dependencies (Fig. 3a, c), a power-law relation emerges between the front velocity scaled by the permeability and the energy density (log10 [e] = 9.42 − 0.6 log10 [u_front (κ_0 (1 − ϕ_0))−1]; Fig. 3d). Forcing the system with a large energy density, yet lower than the previously predicted liquefaction triggering threshold based on undrained consolidation, generates more compaction, a slower front velocity, and a longer liquefaction event. In contrast, a small energy density input induces only a small change in porosity across the front, leading to a rapid front propagation and a short-lived liquefaction event.

The co-seismic dependency of the compaction front velocity and the amount of compaction on the shaking characteristics (seismic power), differ from the often observed “solidification front”31,38, where a front controls post-seismic consolidation, following undrained liquefaction. In contrast to such post-seismic, post-liquefaction, compaction, in drained liquefaction, the very migration of the co-seismic compaction front is the source of the drained liquefaction above it.

The dynamics of drained liquefaction
The four liquefaction indicators observed in the low De number simulations and experiments demonstrate that liquefaction can initiate under drained conditions. In such cases, efficient drainage is key in facilitating the dynamic rheological change of the soil layer. The upward fluid flow between the compaction front and the free surface generated lithostatic pore pressure gradients (Fig. 2a) that supported the weight of the grains, so that granular contact forces vanished. The loss of grain contacts caused the shear modulus to drop (Fig. 2d) and the shear waves to attenuate (Fig. 2c). Ongoing evacuation of fluid from the compaction front facilitated homogeneous continuous layer settlement (Fig. 2b).

Relying on the inferred drained conditions (De ≪ 1), Eq. (2), the description of the pore pressure evolution can be simplified by neglecting the first term relative to the second and third terms. Consequently, it reduces to a two term equation27:

u_sz = (κ_0 / η) (∂P0 / ∂z), (5)

where u_sz is the downward solid grains velocity. Equation (5) describes a compaction-pressurization feedback whereby the pressure gradient responds to instantaneous grain velocity and holds no memory of the previous pressure state22. In accordance with the prediction of the compaction front model, below the compaction front, where u_sz ≈ 0, no dynamic pressure gradient develops and the total pressure gradient is approximately hydrostatic. Above the front, the grains settle at a uniform velocity, leading to a uniform pressure gradient. The compaction front coincides with the location where u_sz changes from finite to zero (Fig. 1b and insets).

Similar to the undrained end-member, drained liquefaction is triggered by shaking-induced destabilization of the granular skeleton through sliding and rolling over grain contacts. At the lowest position of failure, compaction occurs relative to the stable grains below, potentially prescribing the initiation depth of the compaction front, z_front(t) (see Supplementary Note 3 for more details). The pressure gradient and seepage forces that develop in response to this initial compaction only partially support the weight of the settling grains. As long as the pressure gradient remains smaller than lithostatic, the force balance on the settling grains promotes downward acceleration and faster settlement, leading to greater pressure gradients. Once the pressure gradient reaches lithostatic values, it fully supports the weight of the grains. The force balance over the settling grains is then zero, and the grains continue to settle at a constant, terminal velocity27:

u_szC = (κ_0 / η) (dσ_0 / dz), (6)

where subscript C stands for a constant velocity. This terminal constant velocity dictates the linear compaction trend observed in Fig. 2b.

The timescale associated with pore pressure rise to lithostatic values in the simulations and experiments is short, and likely related to a rapid downward propagating liquefaction front. Such a behavior was previously identified in experiments as an unloading front38,45,46 (see Methods Section “Evaluating the compaction front velocity (u_front), the duration of liquefaction event (t_e) and the surface settlement (ΔH)” and Supplementary Note 6). We observed that the deepest location to which the down-going unloading front reaches correlates with the imposed shaking frequency (Supplementary Note 3), and that the unloading front reaches this deepest location in less than two shear cycles, consistent with previous experiments conducted under drained conditions31,32. Notably, our theory does not predict a depth, or normal stress limit for the depth, at which drained liquefaction could be triggered, and similar dynamics at much greater depths25,39 can be identified (see Supplementary Note 5 for details). Another timescale operating in the system is the time required for an isolated, fully immersed, grain to reach its terminal downward velocity. However, since this timescale is exceedingly small, 10−3–10−8 s (see Supplementary Note 2), the acceleration of a single grain is not a rate-limiting process for triggering drained liquefaction. Recent cyclic triaxial experiments28 found that the number of cycles required to initiate liquefaction under drained conditions is smaller than under undrained conditions, supporting the hypothesis that a pressurization time of the order of a few cycles could be indicative of drained liquefaction initiation. Such a consideration might apply to a recent ground motion analysis showing that, in some cases, the time for liquefaction triggering is as short as ~1.7 s from the onset of recorded earthquake shaking47,48.

Different approaches have been proposed to evaluate soil liquefaction potential. Among these, the shear stress or the earthquake peak ground acceleration (PGA)26 forms the theoretical basis for the widely used “simplified procedure for evaluating soil liquefaction potential”49–51. Other approaches emphasize the shear strain14 or the seismic energy15,52–54 in identifying liquefaction triggering thresholds. Importantly, although the three approaches are mechanically linked, their predictions could differ55. The current numerical liquefaction events show a good correlation between the seismic energy density (which follows PGV2) and settlement magnitude (Fig. 3a). Furthermore, both numerical and experimental results show a good correlation between the seismic energy density and the front velocity (Fig. 3d). An even better correlation is found with the rate of seismic energy density input (seismic power, Fig. 3b, d). Both robust correlations emerged despite inherent differences in the boundary conditions, geometrical setup, and particle shapes between the simulations and experiments, suggesting that, within the framework of drained liquefaction triggering, the seismic energy density, and possibly a new measure, the rate of seismic energy density input (seismic power), can be considered as controlling parameters on the magnitude and duration of liquefaction events.

A leading energy-based approach for evaluating soil liquefaction potential uses the earthquake’s Arias intensity54. While the Arias intensity is a cumulative measure that accounts for the amplitude and frequency content throughout the duration of the earthquake, the rate of seismic energy density input (e / T) considered here, can be interpreted as a quasi-instantaneous Arias intensity or an average power of ground shaking over one shear cycle. The excellent performance of the seismic power in explaining the amount of compaction (Fig. 3b) and the front velocity (Fig. 3d), and in defining the clearest threshold between liquefied and non-liquefied simulations (Fig. 3b) suggests that the drained liquefaction dynamics depends on the momentary power rather than on the cumulative power. This is likely in contrast to undrained liquefaction, which is a cumulative process by nature (the volumetric strain required to initiate liquefaction is accumulated over many shear cycles13), hence it might depend on a cumulative energy measure like Arias intensity.

Drained liquefaction beyond the near field in nature
Our simulations and experiments show that drained liquefaction (with De ≪ 1) is triggered when forced with an energy density < 30 J m−3 and as small as 0.1 J m−3. In natural settings, as well, the De number (Eq. (3)) can be evaluated to be smaller than one. For example, using representative values of a 5 m deep soil layer, comprising 1 mm diameter grains, and assuming a 20 grains thick compaction front, gives De = 10−4–10−1, when the permeability range is κ_0 = 10−9–10−12 m2. Consequently, drained liquefaction initiation can be invoked as a general mechanism to explain field observations of liquefaction beyond the earthquake near field, accounting for the previously puzzling 61% of the events reviewed in refs. 16,17 (Fig. 4a).

We propose that the compaction-pressurization feedback, inherent to the drained compaction front dynamics27, is a pivotal player in neutralizing the energy density threshold. With this feedback, small compaction induced by low energy density56 (or more precisely, low rate of energy density input, e / T), presumably facilitated by failure of the weakest grain contacts57,58, generates the initial pressure gradient. The associated seepage forces partially support the weight of the surrounding grains, weakening their contacts and promoting further sliding between grains, compaction, and pressurization, until a lithostatic pressure gradient is achieved and complete liquefaction occurs.

The data of field liquefaction events16,17 show that the number of recorded events decays relatively rapidly below e = 1 J m−3 and no events are recorded when e < 0.1 J m−3 (Fig. 4a). Others59 observed a similar trend regarding the decay of field liquefaction events with low PGV (proportional to the square root of the seismic energy density44), where no liquefaction was observed below PGV = 0.03 m s−1 (e ≈ 0.5 J m−3). The control of the input energy density and the layer permeability, on the compaction front velocity, can explain these observations: With a lower energy density or a larger permeability, the front velocity increases (Eq. (4)), producing a short-lived liquefaction event (Fig. 4b), that is less likely to be observed or recorded. Furthermore, since under a lower energy density, compaction across the front, ϕ_0 − ϕ_c, is smaller (Fig. 3a), less ground settlement occurs (Fig. 4c), so that the associated hazard in such low-energy cases is potentially reduced. The decay of the number of events, and the limit on documented liquefaction event16,17, could therefore be explained within the drained liquefaction triggering framework as a combination of an asymptotically shorter-lived (Fig. 4b) and smaller-settlement events (Fig. 4c), and a seismic energy threshold below which liquefaction does not occur (Fig. 3a, b). Alternatively, it is possible to define a PGA threshold26,27,60 for liquefaction triggering in the current simulation set, which we find to have an exceptionally low value of Aω / g ≈ 0.05 (see Table 2). The numerically-identified low liquefaction triggering threshold might be partially related to model parameters, including the use of perfectly smooth and spherical grains. More realistic grain shapes could potentially necessitate a higher triggering threshold61. Nevertheless, our experiments with natural sand grains showed liquefaction triggering at PGA values as low as Aω / g ≈ 0.1 (the precise threshold was not investigated). This suggests that the intrinsic characteristics of natural grains are unlikely to alter the outcomes significantly, permitting low PGA liquefaction triggering under real-world conditions. Furthermore, after the triggering phase, once granular contacts are minimized, the dynamics of the liquefied layer are expected to be independent of grain shape, as evidenced by the congruence of front velocities across simulations, experimental results, and theory (depicted in Fig. 3d).

The analysis above does not account for the finite duration of earthquakes, t_EQ, which can be shorter or longer than the event termination time, t_e (dictated by the front velocity, as discussed above). After time t_EQ, the seismic energy input drops to zero. The post-seismic evacuation of the excess pore pressure occurs at a rate that depends mainly on the permeability and, in our simulations and experiments, is accompanied by negligible residual compaction. The real duration of a drained liquefaction event is thus the minimum between t_e and t_EQ, where t_EQ is typically of the order of tens of seconds for moderate to large earthquakes (Fig. 4b, c). When the permeability is high and the seismic energy low, t_e < t_EQ, the total compaction is dictated by the front passing through the whole layer, converting the porosity from ϕ_0 to ϕ_c, where Δϕ is controlled by seismic energy density (sloping gray line in Fig. 4c). However, when the permeability is low, and the seismic energy is high, t_e > t_EQ, the front does not have sufficient time to sweep the whole layer during the earthquake. The total compaction, in this case, is independent of seismic energy density, and instead depends on permeability. The permeability controls the settling velocity of grains in the liquefied sub-layer (Eq. (6)), and the total compaction is then the integral of the settling velocity over time t_EQ (horizontal lines in Fig. 4c).

Finally we discuss briefly two common scenarios of more complicated geometries: 1) when the water table does not coincide with the surface, and 2) a soil with layered permeability. In our simulations and experiments, the water table was taken to coincide with the grain layer free surface. In such settings, fluid expulsion out of the soil layer will start concurrently with drained liquefaction initiation. However, if the water table lies much below the surface, fluid expulsion out of the soil during drained liquefaction initiation could be delayed. In cases where the water table is sufficiently deep, and the energy density is relatively low (inducing only a small Δϕ), the water table might not reach the surface during liquefaction, and no fluid expulsion out of the surface would be observed. Coseismic settling of the ground surface, on the other hand, will take place even in the absence of fluid expulsion. Delayed fluid expulsion should therefore, not be a-priori considered as an indication for undrained liquefaction initiation followed by a breach of low permeability barrier, or as an indicator for liquefaction by pore pressure diffusion from a distant source7,16,18. Instead, it could also be consistent with a drained initiation if the water table was originally relatively deep (see Supplementary Note 4).

In our simulations and experiments, we did not include sub-layers with variable permeability. Centrifuge experiments62 and numerical simulations63 show that a water-film may form below a low permeability seam. If the fluid flow upwards across the seam is slow and the seam is not broken yet, the water-film may change its volume to preserve a constant lithostatic water pressure boundary condition for the layer below the seam (by “pushing” the seam upwards62). This promotes a behavior very similar to the presented drained compaction front, where the higher permeability sub-layer below the seam is compacting and the fluid drains toward the water film (rather than directly to the surface). A further investigation of such a setting is needed to examine the initial pressurization process and its sensitivity to the seismic energy density.

In conclusion, we summarize that theory, simulations, and experiments demonstrate that drained liquefaction triggering could be invoked to explain ubiquitous and previously puzzling observations of liquefaction beyond the earthquake near-field16, where the seismic energy density is lower than an empirical threshold inferred for the onset of undrained liquefaction.

Effectively drained conditions are associated with a pore pressure diffusion timescale shorter than the grain skeleton deformation timescale, such that the non-dimensional Deborah number, De ≪ 1. Such conditions represent combinations of site and event properties, where high permeability throughout the soil column and low shaking frequency contribute to the drained conditions. The latter is also associated with smaller energy density input. This suggests that the drained initiation end-member likely dominates far away from the earthquake epicenter (in the intermediate to far-field).

The dynamics of drained liquefaction are controlled by a co-seismic, upward-migrating compaction front that induces co-seismic, spatially heterogeneous compaction. Theoretical predictions, simulations, and experiments further reveal that the compaction front velocity inversely correlates with the seismic energy density, and shows an even better correlation with the rate of seismic energy density input, a new seismic-intensity-based control parameter. Consequently, exceedingly small events with energy density ≲ 0.1 J m−3 will be characterized by a rapid compaction front and small co-seismic compaction, potentially hindering field detection of liquefaction and explaining the lack of documented liquefaction events at very low energy density16.

A critical implication of this study is that liquefaction potential and risk evaluation should account for the possibility of drained liquefaction triggering, with its general mechanistic model for liquefaction of well-drained soils, i.e., young fluvial and beach sediments and reclaimed lands, under exceedingly small seismic energy density input.

Methods
A general description of the grain-fluid system
We study the coupled grain-fluid dynamics of a fully saturated granular layer subjected to 1D horizontal harmonic shear displacement. The top of the layer is unconfined. Horizontal cyclic shear is applied to the base of the layer, which acts as a no-flow boundary for the fluid. This geometry represents a shallow soil layer overlaying a bedrock that is agitated by an upward traveling horizontally polarized seismic shear wave.

Numerical simulations
The numerical approach is described in 22,27. Here, we repeat its main details. We use a two-phase coupled model. The grains are modeled using the discrete element method64, and the interstitial pore fluid is modeled as a continuum on a superimposed Eulerian grid34,40,65–68. Grain velocity and position are resolved by time integration of the linear and rotational momentum conservation equations27:

m_i u_s,i = m_i g + V_imm,i ρ_f g + Σ_j F_ij − ∇ P_0 ϕ V_i (7)

I_i ω_s,i = Σ R_i n^ × F_ij (8)

where u_s,i and ω_s,i are the translational and rotational accelerations of grain i (dot notation indicates time derivative) and m_i and I_i are the mass and moment of inertia of grain i. R_i is the radius of grain i and n^ is a unit vector along the direction connecting the centers of grains i and j. In Eq. (7), the first term on the right-hand side is the gravitational force, where g is the gravitational acceleration. The second term is the buoyancy force, whose magnitude depends on the grain immersed volume V_imm,i and the fluid density ρ_f26. The third term is the sum of contact forces (F_ij) over all grains j that are in contact with grain i, calculated with a linear contact model64. The fourth term represents the seepage force exerted by the gradient of the dynamic pore pressure, ∇P_0, where V_i is the volume of grain i.

The evolution of the interstitial fluid pressure is represented by Eq. (1)22, which is solved by using an implicit scheme over a square grid, with a grid spacing of two average grain diameters21,22,34,66. No a-priori assumption is made regarding the value of the De number (Eq. (3)), and the full three term equation is solved.

The two-way coupling between the grains and the fluid is implemented as follows. The fourth term on the right-hand side of Eq. (7) is evaluated via a bilinear interpolation of ∇P_0 / (1 − ϕ) from the surrounding grid nodes to grain i. The second and third terms of Eq. (1) are evaluated by defining smooth fields of grain velocity and porosity over the grid through a bi-linear interpolation of grain radius and velocity from individual grains surrounding each grid node. The permeability, κ, in Eq. (1) is calculated based on a three dimensional Kozeny–Carman relation69:

κ = κ_0(x, y, t) = α r^2 ϕ^3 / (1 − ϕ)^2, (9)

where r^2 is the bi-linearly interpolated squared grain radii in the surroundings. κ_0 = α<r>^2 is a constant prefactor, and r_0^2 f(ϕ) captures permeability variations in space and time. <r> is the mean grain radius in the system and r_0 is the local deviation from it, such that r = <r> r_0. In the original Kozeny–Carman relation, α = 1/4569 is a geometrical prefactor for spheres. In our simulations, we vary α to directly control the order of magnitude of the permeability independent of the grain size22,27.

The geometry of the numerical layer (Fig. 5a) is a Hele-Shaw cell comprising spherical grains with grain radii between 0.8–1.2 cm, drawn from a normal distribution with a mean of 1 cm and a standard deviation of 1 cm. The system’s horizontal dimension is 0.4 m. The layer is prepared as follows: First, a target height is specified. Then, grains are sedimented under gravity onto the bottom wall in a fluid-free environment. Next, to slightly compact the layer, a short horizontal shaking phase is applied over 0.62 s with f = 12 Hz and amplitude of A = 0.0431 cm, followed by 0.13 s relaxation, where no external forces aside from gravity are applied. Finally, the fluid is added so its height approximately coincides with the top of the grain layer and the layer is relaxed again. In the simulations presented here, the initial layer height, following the preparation stage, is H ≈ 1.44 m.

The bottom wall of the numerical Hele-Shaw cell is made of half grains glued together. The boundary condition for the bottom wall is zero velocity in the vertical direction (u_sz(z = bottom, t) = 0) and sinusoidal displacement in the horizontal direction, x(z = bottom, t) = A(1 − cos(ωt)), where A and ω are the shearing amplitude and angular frequency, respectively. At the top boundary, there are no normal or shear stresses. The boundary conditions for the fluid phase are no flow boundary at the bottom (∂P_0 / ∂z(z = bottom, t) = 0) and constant pressure boundary at the top (P_0(z = top, t) = 0). The water level is maintained at its initial height throughout the simulation. The side boundaries are periodic for the grains and pore fluid, mimicking a laterally infinitely long layer.

Table 1 summarizes the simulations’ parameters. Table 2 lists the simulations presented here with their applied shear amplitude and frequency. The pressure signal in Fig. 2a is smoothed over a window of two cycles. The compaction in Fig. 2b is calculated as the time integral of the mean vertical velocity of grains in the topmost sub-layer.

Experiments
The experiments (Fig. 5b–e) comprise a 12 × 12 × 12 cm3 transparent box. The box is attached to a horizontal shaker (Tira® S51120) fed with a harmonic signal from a signal generator (Agilent® 33220A) through an amplifier (BAA500). The box’s face perpendicular to the shaking direction is filmed by a high-speed camera (Photron® SA5) at a rate of 250 frames per second. The frames are analyzed by using MATLAB® image processing toolbox and PIVlab70,71, an open-source MATLAB® toolbox, to identify changes in the layer’s height and define instantaneous grain velocity. An array of three pressure transducers (Honeywell 24PC) is mounted vertically on the opposite parallel face of the box at depths of 1, 3.6, and 6.1 cm above the box base.

Before the experiment starts, the experimental box is filled with tap water, and the pressure transducers are calibrated under hydrostatic conditions. Then, sand grains with a mean diameter of 200 micrometers (SIFRACON 34) and density ρ_s = 2650 kg/m3 are gradually poured into the box. We aim for a situation where the water table and the top grains approximately coincide. Horizontal shaking is applied for 30 s with a displacement amplitude of at least one mean grain diameter. The pore pressure is measured at a frequency of 104 Hz from 30 s before the application of shaking and until 240 s after shaking stops. The pressure at the top is atmospheric, such that the top boundary is fully drained. All the other box faces exert no flow conditions.

Table 1 summarizes the experiments’ parameters. Table 3 lists the experiments presented here with their applied shear amplitude and frequency. The initial porosity, ϕ_0, presented in Table 3 is evaluated as follows: First, the pore water volume is evaluated as the difference between the water volume used in the experiment and the volume of the thin water film above the grains. Then, the pore water volume is divided by the total volume of the saturated grains based on the height of the grain layer as recorded by the first high-speed camera image. The mean value of the pre-shaking pressure measurements is used to determine the hydrostatic pressure reference. The pressure signal in Fig. 2a is filtered using a low pass filter with a cutoff frequency of 1 Hz. The shaded red area in Fig. 2a, representing the uncertainty on the pressure, follow the 95% confidence bounds on the parameters of the linear regression between voltage and pressure based on the pressure transducers calibration stage.

The calculation of the normalized compaction in Fig. 2b is based on an edge detection algorithm that identifies the top boundary (edge) of the grain layer in individual images. The algorithm was executed several times while varying the top boundary of the search frame within which the algorithm searches for the edge. From search frame 211 (corresponding to a scaled time of t / t_e = 0.084) and on, the edge becomes independent of the search frame height, and thus Fig. 2b (red curve) shows the compaction trend only from frame 221. The red curve represents the normalized mean edge topography within the search frame (smoothed by a moving average window of 0.8 s). The shaded red band represents the uncertainty on the normalized edge height based on the standard error of the edge topography. The mean and standard error of the edge topography in the first search frame are based on a manual edge extraction.

The instantaneous grain velocity field is measured in every frame by using PIVlab70,71, which relies on sub-frame correlation between timely-adjacent frames. The vertical velocity is then averaged over sub-layers, yielding the vertical velocity of grains as a function of depth and time (see Section “Evaluating the compaction front velocity (u_front), the duration of liquefaction event (t_e) and the surface settlement (ΔH)”). To minimize boundary effects from the box’s walls, the averaging is done only close to the box’s center (approximately in the middle 2/4 of the box’s total width).

The energy, e, in Fig. 4c, is based on an estimation of the imposed PGV in the experiments. The input shaking frequency was accurately controlled by setting the frequency of the shaker. The shaking amplitude was estimated based on four markers placed close to the corners of the experiment box. Markers’ position was traced across frames. The temporal mean of the markers position was subtracted from the position time series of each marker, and the four position time series were averaged. Then, the peaks of the combined, averaged time series were extracted, and the shaking amplitude was estimated as the average over the absolute value of the peaks through time t = [0, t_e]. The uncertainty in evaluating PGV is related to the standard error of the absolute value of the peaks time series. The error propagated to log10 [e] is smaller than the symbol size in Fig. 4c.

The permeability in the experiments was evaluated based on five static permeability tests. A constant head was applied across a saturated sand layer in each test, prepared similarly and with the same geometry as the shaking experiments. The outlet point that was located 1.4 cm above the base of the box imposed a 3D porous flow field in the box. The cumulative outflow was measured through time, and its time derivative was used as the discharge (with units of m3 s−1) in a 1D Darcy’s law to determine the permeability. A correction factor from a true 1D porous flow to the specific 3D flow structure in these tests was derived by simulating the two geometries in COMSOL Multiphysics. For the same material permeability, the discharge in the 3D geometry was smaller by a factor of 10 with respect to its 1D counterpart. The permeability of each experiment was then estimated as being larger by a factor of 10 with respect to the measured quantity. The hydraulic head was varied between the five experiments, and the permeability used in Fig. 4c is the mean over the five measurements. The uncertainty on the permeability is evaluated as the standard error over the five permeability measurements. When propagated to the y-axis of Fig. 4c, log10 [u_front / (κ_0 (1 − ϕ_0))], the uncertainty is smaller than the symbol size.

Evaluating the compaction front velocity (u_front), the duration of liquefaction event (t_e) and the surface settlement (ΔH)
The compaction front velocity is defined based on the ratio between the horizontally averaged vertical grain velocity and the grain terminal velocity defined by Eq. (6). Averaging is performed over sub-layers of two average grain diameter thickness. The averaged velocity is smoothed in time using a running average window of ~0.67 s in the simulations and ~1.2 s in the experiments. In the simulations, we further smooth the vertical dimension using a running average window of 10 cm. Plotting the averaged, smoothed and normalized velocity as a function of depth and time, u_norm_sz(z, t) results in a map that highlights settling vs. stagnant grains (Fig. 6).

The compaction front depth at each time, t, is extracted by scanning u_norm_sz(z, t) from the bottom upward and identifying the first depth where u_norm_sz(z, t) ≥ 0.5 in simulations and ≥ 0.01 in experiments. This depth is defined as the front location at time t, z_front(t). Finally, we manually pick the time when the front starts migrating upward continuously (t_i) and the time of front arrival to the top of the layer (t_f). In some experiments and simulations, z_front(t) loses its coherent slope at some stage, and we choose t_f to be the last time step showing a coherent slope. The average front velocity is calculated as the average slope of the front depth-time curve between these two times.

The time of liquefaction event termination, t_e, used as the timescale factor in Fig. 2, is determined as the time at which the soil compaction slows down significantly (see Fig. 2b), as an approximation of t_e = t_i + L / u_front, where L is the distance to the surface from z_front(t = t_i). In most cases, t_i is negligible in comparison to L / u_front since the downward moving unloading front (Fig. 6 inset) is very fast.

The theoretical time of a liquefaction event in Fig. 4b (semi-transparent symbols) was extrapolated to a 5-m deep soil layer, based on the relationship between the measured compaction front velocity in simulations and experiments, and the input seismic energy density. The synthetic high permeability data (κ = 10−9 m2) was calculated based on the power-law prediction derived from simulations (Fig. 3d and section “Drained liquefaction beyond the near field in simulations and experiments”), by substituting the numerical characteristic permeability with the synthetic permeability. The actual time of liquefaction (opaque symbols in Fig. 4b) was calculated as the minimum value between the theoretical time (t_e) and earthquake duration of t_EQ = 30 s. The sloping blue line is the best fit for the data points which satisfy t_e < t_EQ. The surface settlement in Fig. 4c was calculated as the integral of the settlement velocity (Eq. (6)), between t = 0 and t = min(t_e, t_EQ). The blue and red horizontal trend lines represent the mean value of data points associated with t_e > t_EQ. The gray sloping line represents the amount of settlement resulting from the synthetically adapted power-law prediction, as described above.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Seismic liquefaction of sand is a quite complex phenomenon. Its generation, development and\ndissipation are mainly restricted by the physical properties, stress state and boundary conditions of\nsoil. There are many influencing factors, such as soil density, soil structure, clay content,\nsaturation, gradation, water permeability, initial stress state and dynamic load characteristics.\nGenerally speaking, the mechanism of earthquake liquefaction can be expressed as follows：\nThe shear wave propagates upward from the underlying rock and soil layer during the earthquake,\nand causes alternating stress in the soil, resulting in vibration pore water pressure, which is the\nmain reason for the liquefaction of saturated sand. Under the action of alternating stress, new\nstress will be generated at the contact point of soil particles. When this stress reaches a certain\nvalue, the original connection and structural state between soil particles will be destroyed, and the\nsand particles will be separated from each other. At this time, the stress originally transmitted by\nsand particles through the contact point will be transmitted to the water in the pores to bear,\nresulting in the increase of pore water pressure. With the increase of the number of stress cycles,\nthe pore water pressure increases due to the gradual accumulation. On the one hand, pore water\ntries to discharge upward under a certain vibration pore water pressure; On the other hand, soil\nparticles try to sink down under the action of self weight, so that at the moment of structural\ndamage or within a certain period of time, the downward sinking of soil particles is hindered by the\nupward discharge of pore water, so that the soil particles are partially or completely suspended,\nand the shear strength disappears locally or completely, which means that the soil has varying\ndegrees of deformation or complete liquefaction.\nIn other words, the seismic liquefaction of sand is related to the ground motion excitation and the\nproperties of soil. The ground motion excitation here mainly refers to the peak acceleration and\nthe duration of earthquake. For the earthquake liquefaction phenomenon in the region far away\nfrom the epicenter mentioned by the authors, being far away from the epicenter does not mean\nthat liquefaction will not occur, and the peak acceleration in the region far away from the epicenter\nmay not be low. Studies have found that liquefaction occurs when the earthquake intensity is as\nsmall as moment magnitude 4.5 (Green and Bommer, 2020). The existing mechanism of seismic\nliquefaction mentioned by the authors is based on undrained conditions. In fact, drainage and\nundrained conditions are not unnecessary conditions for seismic liquefaction of sand. Undrained\nconditions are only to explain that pore water pressure offsets gravity. Even under drainage\nconditions, pore water pressure may dissipate slowly, but sand liquefaction is slow. Soil structure\nis also an important factor to control the liquefaction of sand. After the structure is damaged by\nearthquake vibration, the vibration pore water pressure gradually increases with the increase of the\nnumber of stress cycles, and its size is finally enough to make the saturated sand lose its\nshear capacity locally or completely. That is, the seismic liquefaction of sand is affected by the\ncompactness of soil. Even if there is a new mechanism of sand seismic liquefaction, it should be\nbased on particle mechanics, such as the so-called shear jamming mechanism. The innovation of\nthis study is limited.\nIn general, according to the current status of this study, the reviewer regrets that, it is not\nrecommended to be published in this professional journal.", "reviewer_2": "The manuscript titled “Drainage explains soil liquefaction beyond the earthquake near-field”\nexamines the phenomenon of earthquake-induced liquefaction, focusing on the importance of co-\nseismic drainage, with implications on the conventional view of liquefaction triggering.\nMethodologically, it combines analytical with numerical and experimental work. It is a very\nwelcome addition to the growing body of work that challenges the traditional view of liquefaction\nas an undrained phenomenon. The results reported are noteworthy and I expect them to be of\nsignificance to the field. The work presented supports the conclusions and claims. I would\nrecommend the manuscript for publication subject to the following points being addressed:\n1. Line 109: the description of an upwards moving “compaction” front – referred to in literature\nalso as a “solidification” front – which separates the layer it two parts, a compacted, or solidified,\nbelow the front and an uncompacted, or unsolidified, agrees with literature (also with references\n28, 33, 34), though it is often referring to post-liquefaction reconsolidation, following an undrained\nassumption for excess pore water pressure generation. Following your compaction front analysis,\nhow can the initial, downwards propagation of a front separating liquefied from non-liquefied soil\nbe explained (e.g. [33])? Is this an aspect of response that can be captured with your analytical\napproach? This downwards propagation of a “liquefaction” front precedes the reconsolidation\nprocess, where an upwards moving “solidification” (or here “compaction”) front is observed – this\nreconsolidation process seems to be explained well with your compaction front analysis and indeed\nexperimental observation shows that reconsolidation typically starts during the earthquake and not\nnecessarily after its end, agreeing with your view of liquefaction as a “drained” phenomenon.\n2. Some further explanation on when and where the compaction front is formed when applying\nyour analytical (or theoretical) approach would be welcome. In some cases (e.g. Fig2c, Experiment\ne2) only the top part of a deposit liquefies. Is it possible to apply your analytical approach to\npredict the greatest depth where soil is liquefied, and when this happens, or is this taken as a\ngiven and the compaction front analysis starts from there?\n3. Line 112: The assumption of stationary particles below the compaction front is essentially the\nassumption of Florin and Ivanov [33] but experimental data from centrifuge testing (e.g. [28],\n[29]) has shown that the pore water gradient below the front is not actually hydrostatic, so that\nfurther consolidation does take place below the “compaction” front, meaning that compaction\nand fluid flow would still take place in region 1 of Figure 1. In any case, even if the pore water\ngradient were hydrostatic, some compaction is inevitable as the effective stress below the compaction front\ngradually increases with the upwards propagation of the front. Though it is understood that some\nsimplification is necessary, a comment to clarify to the reader the potential limitations of the\nassumption made would be welcome.\n4. Figure 3b: is the equation plotted with a black line eq3 or eq4?\n5. References [15] and [16] are cited regarding the lower bound of energy density required to\ncause liquefaction. Cases below 30 J m^{-3} are given the term “unexplained”. Based on this\nassessment, 61% of events examined in figure 4 are unexplained – an alarmingly high number.\nHowever, in practice the most common methodology to assess liquefaction triggering is the stress\nbased simplified approach, introduced by Seed & Idriss (1967). See Boulanger & Idriss (2016) for\na recent iteration of the method. This method gives a factor of safety for liquefaction occurrence.\nWould 61% of the events you examined remain “unexplained” if the simplified method was used?\nIt would be very useful for the reader to know the number of “unexplained” events based on the\nmethodology that is typically used in practice.\nSeed, H. B. & Idriss, L. M. (1967). Analysis of Soil Liquefaction: Niigata Earthquake. Journal of the\nSoil Mechanics and Foundations Division 93, No. SM3, 83–108.\nBoulanger, R. W. & Idriss, I. M. (2016). CPT-Based Liquefaction Triggering Procedure. Journal of\nGeotechnical and Geoenvironmental Engineering 142, No. 2, 04015065\n6. Line 253-254: Some clarification is needed in this phrase to avoid confusion. You state that\nno further significant compaction is accumulated after the end of the earthquake. You should\nclarify if you are referring to the whole liquefiable deposit and if so, how you explain significant\npost-earthquake settlements observed during post-earthquake reconsolidation e.g. [33], [28], [29].\n7. Line 357: Typically in liquefaction experiments, laminar boxes are used, where each horizontal\nlamination can move relative to the others. This would be closer to the periodic boundary used for\nthe numerical simulation. Given the practically rigid box used and the small dimensions compared\nto similar experiments in literature (as already cited in this manuscript), a comment on potential\nboundary effects is necessary. Moreover, it is implied that the input motion was harmonic but it\nwould be better if this was clearly stated in the text.", "reviewer_3": "This is an interesting paper that aims to explain the soil liquefaction in the far-field where the\nearthquake shaking intensity is small. It is true that such events have been observed in previous\nearthquake events and reported in case histories. Several theories exist for that such as soil\namplification in some regions of the far-field that cause liquefaction while attenuation of\naccelerations in other regions may not lead to full liquefaction.\nThe following comments are offered for the authors to consider.\n1. It may be better to distinguish between ‘liquefaction’ and ‘fluidisation’ throughout the\nmanuscript. I would use ‘liquefaction’ for excess pore pressures generated due to seismic shaking\nand ‘fluidisation’ for loss of shear strength of soil due to pore fluid migration (from lower layers\ntowards the soil surface) and can happen into regions where shaking intensity is not high.\nAuthors may want to clarify that when they said ‘drained liquefaction’ if they mean ‘fluidisation’.\n2. In the same vein, it may be better to distinguish clearly between the ‘drained’ and ‘undrained’\nterms. In liquefaction research (see Refs 22, 24 and 29) the term ‘partially drained’ was\nintroduced to demarcate this. In this paper the authors used the Debroah number (De) to indicate\npartially drained events.\n3. The authors may want to clarify what are the time scales for the liquefaction and fluidisation\nevents. One would expect liquefaction to happen (for loose sands) in the first few cycles of\nearthquake shaking. Fluidisation can happen and continue into the post seismic period.\n4. The concept of ‘liquefaction front’ moving from soil surface towards the bedrock as the shaking\nstarts and the ‘densification front’ moving from the bedrock towards the soil surface have been\npreviously identified by many researchers (e.g. ref 33). There is a large body of dynamic\ncentrifuge test data to show these fronts experimentally occurring in layers of ~ 10m deep. The\nauthors may want to include some of these references into the manuscript.\n5. There is considerable body of research on liquefaction that establishes cyclic shear strain\namplitude threshold i.e. earthquakes that cause shear strains above this threshold cause excess\npore pressure generation while those below won’t. I feel this is more physical as shear strains are\nphysically linked to the movement of individual soil grains and the length parameter. Incidentally I\nagree with the authors argument about length scale in Sec. 3.1. This almost corresponds to\nTatsuoka’s observations of formation of shear bands and their thickness in sands (~20 particle\ndiameters). However, I do not see why the energy argument is used in this paper. Measures like\nArias intensity that depend on the PGV are global measures, where as liquefaction (in my view)\ndepends on the relative displacement between particles. I am not talking about fluidisation here\n(as distinguished in point 1 above). The authors may want to argue why PGV is a good measure\nfor liquefaction or otherwise.\n6. In Sec. 4.1 the authors argue that the time it takes to reach ‘full drained liquefaction’ is\nexceedingly small i.e. 10-3 ~ 10-8 seconds. I presume that the authors mean that this is different\nfrom ‘seismic liquefaction’ and are saying fluidisation can occur much more quickly. Firstly, for\nseismic liquefaction, the excess pore pressure generation can take anything from 1 to 2 cycles to\nmany cycles (~8 to 10 cycles) depending on the relative density of the sand. Loose sands liquefy\nmore easily than dense sands. Secondly, the time for fluidisation surely depends on the\npermeability, porosity of the soil strata and pressure gradient applied which will yield much longer\ntime scales, as pore fluid has to migrate from one region to the other?\n7. I have several objections on the experimental program reported in Sec. 10. Also not sure why\nthe authors do not use data from dynamic centrifuge tests on liquefaction.\na) The experiments are in a very small-scale model container 120 mm (cube). The boundary\neffects will be significant in such a small model container.\nb) What is the relative density of the soil specimens? This should be included in Table 2. Porosity is\nnot convenient for this as RD is measure of the soil state with respect to the loosest and densest\npossible states.\nc) The frequency of shaking is very high (10 Hz, most earthquakes are in the 1 to 5 Hz range)\nd) The amplitude of shaking is very small 0.8 mm is the largest amplitude used. How does this\ncompare to the mean particle size (D20 or D50) for this sand?\ne) Referring to Fig.2a, the excess pore pressures do not show any suction cycles, cyclic oscillations\n(to indicate the crossing of phase transformation line (following Ishihara) nor reach full liquefaction\nof 1. The authors may want to explain why.\nf) The time for generation of excess pore pressure in Fig. 2a is much more slower than in the\nsimulations and magnitudes are smaller.\ng) In fig. 2d, the shear strains are biased towards the positive side? Why?\nh) In Fig.2d, are the experimental results and simulations at the same frequency? The\nrecoupling and increase in oscillations in the experiment are not clear."}, "author_response": "Dear reviewers,\nWe wish to thank you for putting the time and effort in reviewing our manuscript. The comments\nyou made are significant and helped us notice some previously hidden aspects of our work, as\nwell as highlight several “blind spots” that we had when communicating our study.\nWe believe that the changes we have introduced to address your comments have greatly\nimproved the manuscript, both scientifically and in terms of science communication. We hope that\nyou’ll see those changes similarly to how we perceive them.\nSincerely,\nShahar Ben-Zeev on behalf of the co-authors\n\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nSeismic liquefaction of sand is a quite complex phenomenon. Its generation, development and\ndissipation are mainly restricted by the physical properties, stress state and boundary conditions\nof soil. There are many influencing factors, such as soil density, soil structure, clay content,\nsaturation, gradation, water permeability, initial stress state and dynamic load characteristics.\nGenerally speaking, the mechanism of earthquake liquefaction can be expressed as follows：\nThe shear wave propagates upward from the underlying rock and soil layer during the\nearthquake, and causes alternating stress in the soil, resulting in vibration pore water pressure,\nwhich is the main reason for the liquefaction of saturated sand. Under the action of alternating\nstress, new stress will be generated at the contact point of soil particles. When this stress\nreaches a certain value, the original connection and structural state between soil particles will\nbe destroyed, and the sand particles will be separated from each other. At this time, the stress\noriginally transmitted by sand particles through the contact point will be transmitted to the water\nin the pores to bear, resulting in the increase of pore water pressure. With the increase of the\nnumber of stress cycles, the pore water pressure increases due to the gradual accumulation.\nOn the one hand, pore water tries to discharge upward under a certain vibration pore water\npressure; On the other hand, soil particlestry to sink down under the action of self weight, so\nthat at the moment of structural damage or within a certain period of time, the downward sinking\nof soil particles is hindered by the upward discharge of pore water, so that the soil particles are\npartially or completely suspended, and the shear strength disappears locally or completely,\nwhich means that the soil has varying degrees of deformation or complete liquefaction.\nIn other words, the seismic liquefaction of sand is related to the ground motion excitation and\nthe properties of soil. The ground motion excitation here mainly refers to the peak acceleration\nand the duration of earthquake. For the earthquake liquefaction phenomenon in the region far\naway from the epicenter mentioned by the authors, being far away from the epicenter does not\nmean that liquefaction will not occur, and the peak acceleration in the region far away from the\nepicenter may not be low. Studies have found that liquefaction occurs when the earthquake\nintensity is as small as moment magnitude 4.5 (Green and Bommer, 2020). The existing\nmechanism of seismic liquefaction mentioned by the authors is based on undrained conditions.\nIn fact, drainage and undrained conditions are not unnecessary conditions for seismic\nliquefaction of sand. Undrained conditions are only to explain that pore water pressure offsets\ngravity. Even under drainage conditions, pore water pressure may dissipate slowly, but sand\nliquefaction is slow. Soil structure is also an important factor to control the liquefaction of sand.\nAfter the structure is damaged by earthquake vibration, the vibration pore water pressure\ngradually increases with the increase of the number of stress cycles, and its size is finally\nenough to make the saturated sand lose its shear capacity locally or completely. That is, the\nseismic liquefaction of sand is affected by the compactness of soil. Even if there is a new\nmechanism of sand seismic liquefaction, it should be based on particle mechanics, such as the\nso-called shear jamming mechanism. The innovation of this study is limited.\nIn general, according to the current status of this study, the reviewer regrets that, it is not\nrecommended to be published in this professional journal.\nWe agree with many of the reviewer’s comments describing soil liquefaction and\nthe processes associated with it. Yet, the undrained perspective, which has been\nperceived as an axiom in the liquefaction research community, has been recently\nquestioned (e.g., Goren et al 2010,2011; Madabhushi & Haigh 2012; Lakeland 2014,\nOdamidis & Madabhushi 2016,2018). Here, we propose the non-dimensional\nDeborah (De) number framework to estimate the soil drainage conditions and\nexpected response, based on basic physics. This allows questioning the validity of\nthe undrained assumption in specific setups, geometries, and boundary conditions.\nWe further show that pressurization can be achieved even under well-drained\nconditions when the fluid flow itself produces the pressure gradients and pore\npressure values needed to liquefy the soil. This should be distinguished from\nliquefaction initiation under undrained conditions that precede a post-seismic\ndissipation of pore pressure.\nThe current work addresses the puzzling nature of liquefaction beyond the\nearthquake near-field, aiming to explain and predict the many observed liquefaction\nevents which occur under a very low seismic energy input (Wang, 2007), much\nbelow the minimum input required to liquefy soils in undrained lab tests (Green &\nMitchell 2004). We show that drained liquefaction is not limited by the undrained\nseismic energy threshold and thus can be invoked to explain the occurrence of\nliquefaction under very low seismic energy input (with values typical to the\nearthquake’s intermediate and far fields).\nWe think that the community will benefit from studies like this, that promote\nrethinking of basic assumptions like the undrained approach. Especially, since as\ncited above, there is already a non-negligible amount of research pointing in that\ndirection. Our analytical approach is laid out and supported by numerical\nsimulations and experiments. It should be noted, that we do not aim to replace the\nundrained mechanism, we simply demonstrate that liquefaction can be triggered\nunder a full spectrum of drainage conditions, and that drained liquefaction can be\nthe missing piece in liquefaction related puzzles (e.g., far-field liquefaction).\nReviewer #2 (Remarks to the Author):\nThe manuscript titled “Drainage explains soil liquefaction beyond the earthquake near-field”\nexamines the phenomenon of earthquake-induced liquefaction, focusing on the importance of\nco-seismic drainage, with implications on the conventional view of liquefaction triggering.\nMethodologically, it combines analytical with numerical and experimental work. It is a very\nwelcome addition to the growing body of work that challenges the traditional view of liquefaction\nas an undrained phenomenon. The results reported are noteworthy and I expect them to be of\nsignificance to the field. The work presented supports the conclusions and claims. I would\nrecommend the manuscript for publication subject to the following points being addressed:\n1. Line 109: the description of an upwards moving “compaction” front – referred to in literature\nalso as a “solidification” front – which separates the layer it two parts, a compacted, or solidified,\nbelow the front and an uncompacted, or unsolidified, agrees with literature (also with references\n28, 33, 34), though it is often referring to post-liquefaction reconsolidation, following an\nundrained assumption for excess pore water pressure generation.\nWe agree with this comment. To clarify this point, we added a paragraph at the end\nof section 2:\n“This compaction front is in many aspects similar to previously\ndescribed post-liquefaction \"solidification fronts\", which are preceded by an\nundrained liquefaction initiation and result in consolidation (Florin &\nIvanov, 1961; Wang et al. 2013). In the drained end-member described here, the\ncompaction front forms and migrates co-seismically, and its very migration\ncauses liquefaction above it, as demonstrated next.”\nFollowing your compaction front analysis, how can the initial, downwards propagation of a front\nseparating liquefied from non-liquefied soil be explained (e.g. [33])? Is this an aspect of\nresponse that can be captured with your analytical approach? This downwards propagation of a\n“liquefaction” front precedes the reconsolidation process, where an upwards moving\n“solidification” (or here “compaction”) front is observed – this reconsolidation process seems to\nbe explained well with your compaction front\nanalysis and indeed experimental observation shows that reconsolidation typically starts during\nthe earthquake and not necessarily after its end, agreeing with your view of liquefaction as a\n“drained” phenomenon.\nAt this stage, our theory cannot predict analytically the downward progressing\nliquefaction front. However, following this suggestion, we looked for and identified\nthe liquefaction front in our liquefied simulations. A sample simulation exhibiting\nthe very rapid liquefaction down-going front is now shown in Fig 6, inset. Our\nlimited numerical resolution prevents us from studying it thoroughly. We do see\nthat this liquefaction front controls the liquefaction onset, since the upwards\ncompaction fronts initiates from the deepest point that the liquefaction front\nreaches. Similar to our simulations, in the experiments presented in Wang et al.\n(2013) and Adamidis & Madabhushi (2016) the liquefaction front is very fast and\nalmost undetectable. Following this suggestion of the referee, we added a\ndiscussion related to the existence of the downward liquefaction front as part of\nsec 4.1 which discussed the initiation timescale of liquefaction:\n“The timescale associated with pore pressure rise to lithostatic\nvalues in the simulations and experiments is short, and likely related to a\nrapid downward propagating liquefaction front. Such a behavior, viewed as\n\"unloading\" fronts (Florin & Ivanov, 1961), was previously identified in\nexperiments (Florin & Ivanov, 1961; Sharp & Dobry, 2002; Gonzalez et al.\n2005) (see Methods section 10.4 and section A in the supplementary\nmaterial). We observed that the deepest location to which the down-going\nliquefaction front reaches correlates with the imposed shaking frequency\n(Appendix C), and that the liquefaction front reaches this deepest location in\nless than two shear cycles, consistent with previous experiments conducted\nunder drained conditions (Wang et al., 2013; Adamidis & Madabhushi,\n2016).”\nIn the revised version, in addition to demonstrating the presence of the downward\nliquefaction front in Fig. 6 (inset), we also analyze the depth to which the\nliquefaction front progresses in a newly added Appendix C.\nFurthermore, a new Supplementary Material section is introduced, which shows the\ngrains vertical velocity over depth and time for all of our simulations. Some of these\nfigures allow identification of the down-going liquefaction front, despite its high\nspeed.\n2. Some further explanation on when and where the compaction front is formed when applying\nyour analytical (or theoretical) approach would be welcome. In some cases (e.g. Fig2c,\nExperiment e2) only the top part of a deposit liquefies. Is it possible to apply your analytical\napproach to predict the greatest depth where soil is liquefied, and when this happens, or is this\ntaken as a given and the compaction front analysis starts from there?\nThe initiation depth of the front is an important topic, and although it was not the\nfocus of this study, following this comment, we added appendix C which shows that\nthe initiation depth in simulations is proportional to the shaking frequency. In\naddition, as explained in our answer to the previous comment, we report a finding\nthat the compaction front initiation depth is the deepest point that the liquefaction\nfront reaches. This is expressed in line 516:\nOur results indicate that 𝒕 may be controlled by the descending\n𝒊\n\"unloading\" liquefaction front that precedes the upward compaction front.\nLiquefaction initiation is only complete once this unloading front reaches the\ncompaction front initiation depth, as seen in Fig. 6 and in section A in the\nsupplementary materials.\nA possible implication for the observed frequency dependence of the initiation\ndepth for the ascending compaction front, is that under some shaking frequencies,\nthinner layers will not liquefy because their thickness is smaller than the\ncompaction front initiation height, for this frequency. This topic however merits\nfurther study, beyond the scope of this study.\n3. Line 112: The assumption of stationary particles below the compaction front is essentially the\nassumption of Florin and Ivanov [33] but experimental data from centrifuge testing (e.g. [28],\n[29]) has shown that the pore water gradient below the front is not actually hydrostatic, so that\nfurther consolidation does take place below the “compaction” front, meaning that compaction\nand fluid flow would still take place in region 1 of Figure 1. In any case, even if the pore water\ngradient were hydrostatic, some compaction is inevitable as the effective stress below the compaction\nfront gradually increases with the upwards propagation of the front. Though it is understood that some\nsimplification is necessary, a comment to clarify to the reader the potential limitations of the\nassumption made would be welcome.\nWe agree that some compaction does occur below the front. Following this\ncomment, we rephrased the relevant sentence using the adverbs “nearly” and\n“approximately” when referring to the grain layer below the (in section 2 and section\n4.1). Lines 127-132 now reads:\n“The propagating compaction front separates two regions within the\nlayer (Ben-Zeev et al. 2020) (Fig. 1a-b): a lower region, in which grains have\ncompacted to porosity 𝝓 < 𝝓 and are approximately stationary in the\n𝒄 𝟎\nvertical direction. In this lower region, the pore pressure gradient is nearly\nhydrostatic, although the pore pressure itself is elevated to the value of the\npressure at the front (Fig. 1c).”\nIt is worth noting that for the physics described in the current study, stagnant grains\nbelow the front are a good approximation, as reflected by the agreement of the\nnumerical and experimental results with the analytic predictions in figure 2a. This\napproximation is also validated by the newly added supplementary material section\n(grain velocity maps).\n4. Figure 3b: is the equation plotted with a black line eq3 or eq4?\nEq. 4. Thank you for the catching the typo.\n5. References [15] and [16] are cited regarding the lower bound of energy density required to\ncause liquefaction. Cases below 30 J m^{-3} are given the term “unexplained”. Based on this\nassessment, 61% of events examined in figure 4 are unexplained – an alarmingly high number.\nHowever, in practice the most common methodology to assess liquefaction triggering is the\nstress based simplified approach, introduced by Seed & Idriss (1967). See Boulanger & Idriss\n(2016) for a recent iteration of the method. This method gives a factor of safety for liquefaction\noccurrence. Would 61% of the events you examined remain “unexplained” if the simplified\nmethod was used? It would be very useful for the reader to know the number of “unexplained”\nevents based on the methodology that is typically used in practice.\nSeed, H. B. & Idriss, L. M. (1967). Analysis of Soil Liquefaction: Niigata Earthquake. Journal of\nthe Soil Mechanics and Foundations Division 93, No. SM3, 83–108.\nBoulanger, R. W. & Idriss, I. M. (2016). CPT-Based Liquefaction Triggering Procedure. Journal\nof Geotechnical and Geoenvironmental Engineering 142, No. 2, 04015065\nThe factor of safety in the stress-based simplified procedure is structured as the\nratio between the CSR (Cyclic Stress Ratio), which represents the seismic loading\non the soil in terms of acceleration/stress, and the CRR (Cyclic Resistance Ratio)\nwhich represents the soil resistance. The CRR is evaluated empirically by plotting\nthe CSR during liquefaction\non-liquefaction vs. in situ soil resistance index (SPT-\nStandard Penetration Test / CPT-Cone Penetration Test/ other) from case history\ndatabase, and drawing the most probable line between liquefaction to non-\nliquefaction. Almost by definition, events that were part of the original database\nfrom which the CRR line was constructed, are “explained” by the simplified\nprocedure.\nIn that sense, it is of interest to look for new events that mismatch the original model\nand can be called “unexplained”. Interestingly, one of the motivations for updating\nthe “stress based simplified approach” (Seed & Idriss 1967), mentioned by the\nreviewer (Boulanger et al. 2016), is the accumulation of new data over the course of\ntime elapsed between the introduction of the approach and present time. Such\nperiodic updating of predictions is an inherent feature of data-driven or semi-data-\ndriven approaches like the “stress based simplified approach”. In contrast, a\nphysics-based approach (as suggested in this study for example), provides an\nindependent prediction, that does not need change with accumulation of data,\npointing to the importance of physically-based studies.\nThe current study aims to provide a physical model to fill an existing gap of\nknowledge, and not to perform a rigorous database study or databases\ncomparisons. The observations we use are from Wang (2007), who demonstrated\nthat many events fall below the expected energetic threshold for undrained\nliquefaction, and by this identified a gap in knowledge. Although not popular as the\nstress-based approach, the energy approach is not new, and pursuing it as an\nalternative was recommended by many experts’ workshops and reviews (e.g., Youd\net al. 2001, National Academies of Sciences, 2016). The energetic threshold marked\nby Wang (2007, after Green & Mitchel, 2004) arises from undrained experiments,\nindependently from the field data.\nThe fact that many real liquefaction events mismatch what is expected from\nundrained prediction (by orders of magnitude) is an excellent demonstration of the\nfailure of the undrained approach to explain all soil liquefaction occurrences,\nmaking it a compelling observation for researchers pursuing an alternative\nmechanism.\nThat being said, we do acknowledge that the use of ‘unexplained’ might be\nmisleading, and following this comment, we change the wording in lines 15&302 to\n‘puzzling’.\n6. Line 253-254: Some clarification is needed in this phrase to avoid confusion. You state that\nno further significant compaction is accumulated after the end of the earthquake. You should\nclarify if you are referring to the whole liquefiable deposit and if so, how you explain significant\npost-earthquake settlements observed during post-earthquake reconsolidation e.g. [33], [28],\n[29].\nWe agree with the reviewer that this important point needs clarification.\nWithin the parameter range of our tests, our drained simulations show (to be\npublished as part of a work in progress about repeated liquefaction) that the amount\nof post-seismic compaction is negligible compared to the co-seismic\ncompaction/settlement. A similar observation regarding the settlement rate was\nalready made by (Madabhushi & Haigh, 2012; Adamidis & Madabhushi, 2016) based\non centrifuge tests.\nThe residual compaction and pressure diffusion after the cessation of seismicity is\nnot a focus of the current study, but we can propose some scenarios where\nsignificant soil settlements in the post-seismic phase or at least long duration post-\nseismic settlement can occur:\nThe first scenario is due to a low permeability or a very thick liquified layer- After\nthe earthquake stops, the rate in which the excess pore pressure dissipates and\nfluid flows out depends mainly on the permeability of the soil. Hence, long pressure\ndissipation is expected where the permeability is low, or the liquefied layer is very\ndeep. Both contributing to a high De number and would tend to push the event\ntoward the undrained limit. This high De is outside the regime of our experiments\nand simulations. Furthermore, in a very thick layer, it is possible that minor local\nchange in porosity, while the residual pore pressure is dissipating, will be\nwitnessed as larger soil settlement of the total layer (by integrating small changes\nover the layer and accumulating a total volumetric strain).\nAnother possibility for witnessing a significant post-seismic settlement, is where\nthere is a low permeability seam (which is not the case in [28,29] but is common in\nthe field, as in [Holzer et al. 1989]). This again poses an undrained scenario, outside\nthe scope of the present study, Even if the seam does not yield, a slow fluid flow\nacross it will act as a pressure source for the layer above it, that can be sustained\nover longer durations (similar to fluidization experiment, or coastal quicksand\ndriven by up-flow). In this scenario, the reduction of fluid volume accumulated\nbelow the seam could contribute to observations of post-seismic settlement.\nFollowing this comment, we modified the text in section 4.2 as follows:\n“An intrinsic timescale for drained liquefaction is the event\ntermination time, 𝒕 , which depends on the front velocity. However, in natural\n𝒆\nevents, the earthquake shaking duration, 𝒕 , is finite, and after time 𝒕 the\n𝑬𝑸 𝑬𝑸\nseismic energy input drops to zero. The post-seismic evacuation of the\nexcess pore pressure occurs at a rate that depends mainly on the\npermeability and in our simulations and experiments is accompanied by\nnegligible residual compaction.”\n7. Line 357: Typically in liquefaction experiments, laminar boxes are used, where each\nhorizontal lamination can move relative to the others. This would be closer to the periodic\nboundary used for the numerical simulation. Given the practically rigid box used and the small\ndimensions compared to similar experiments in literature (as already cited in this manuscript),\na comment on potential boundary effects is necessary.\nThe primary purpose of the experiments was to examine the grain dynamics and\nfluid pressurization visually. To achieve that, we chose to use a completely\ntransparent box to allow filming of at least one face of the box. We succeeded in\nachieving this goal with the current experimental setup, as we were able to validate\nvisually the formation of the compaction front and its velocity relation to the\nimposed energy density (and rate of energy density input). We are aware of the\nlimitations of this experimental setup, but it allowed us to validate the trends that\nemerged from theory and simulations in a way that centrifuge tests or laminar box\nsetup would not allow. Even more so, the observation that despite their limitations\n(small size, boundary effects, rigid), the experiments showed the same dynamics\nas the simulations (the formation of the compaction front and its velocity relation\nwith the imposed seismic energy density and rate of seismic energy input, see\nfigs.2a-d & 3a-b,d), and also agreed perfectly with the theoretical prediction, is a\nstrong support for both the robustness of the compaction front dynamics theory,\nand the suitability of the current experimental setup for studying it.\nAs there are limitations to any experimental choice, the experimental setup here did\nserve its designated purpose.\nFollowing this comment, we added the following text to briefly express these views\n(section 4.1):\n“…Both robust correlations emerged despite inherent differences in\nthe boundary conditions and geometrical setup between the simulations and\nexperiments, suggesting that, within the framework of drained liquefaction\ntriggering, the seismic energy density, and possibly a new measure, the\nseismic energy density rate, can be considered as controlling parameters on\nthe magnitude and duration of liquefaction events.”\nWe also note that to minimize the wall’s effect we limit the grain velocity\nmeasurements to a region close to the box center. Thank you for pointing the\nabsence of this explanation from the manuscript. We added the following text to\nsection 10.3:\n“The instantaneous grain velocity field is measured in every frame by\nusing PIVlab (Thielicke, & Stamhuis 2014; Thielicke & Sonntag, 2021), which\nrelies on sub-frame correlation between timely-adjacent frames. The vertical\nvelocity is then averaged over sub-layers, yielding the vertical velocity of\ngrains as a function of depth and time (see section 10.4). To minimize\nboundary effects from the box's walls, the averaging is done only close to\nthe box's center (approximately in the middle 2/4 of the box's total width).“\nMoreover, it is implied that the input motion was harmonic but it would be better if this\nwas clearly stated in the text.\nThank you. The revised text clarifies this point in section 10.3:\n“The box is attached to a horizontal shaker (Tira S51120) that is fed with harmonic\nsignal from a signal generator (Agilent 33220A) through an amplifier (BAA500).”\nReviewer #3 (Remarks to the Author):\nThis is an interesting paper that aims to explain the soil liquefaction in the far-field where the\nearthquake shaking intensity is small. It is true that such events have been observed in previous\nearthquake events and reported in case histories. Several theories exist for that such as soil\namplification in some regions of the far-field that cause liquefaction while attenuation of\naccelerations in other regions may not lead to full liquefaction.\nThe following comments are offered for the authors to consider.\n1. It may be better to distinguish between ‘liquefaction’ and ‘fluidisation’ throughout the\nmanuscript. I would use ‘liquefaction’ for excess pore pressures generated due to seismic\nshaking and ‘fluidisation’ for loss of shear strength of soil due to pore fluid migration (from lower\nlayers towards the soil surface) and can happen into regions where shaking intensity is not high.\nAuthors may want to clarify that when they said ‘drained liquefaction’ if they mean ‘fluidisation’.\nThe term “liquefaction” or “seismically induced liquefaction” is used throughout\nthe manuscript to describe the loss of shear strength of the soil due to seismic\nshaking-induced pressurization. One of our major findings is that this\npressurization can be achieved via drained liquefaction, which adds to the familiar\nundrained end member. Distinguishing the high De # liquefaction (undrained) from\nthe low De# (drained) liquefaction, by calling each a different name (the 1st\n“liquefaction” and the 2nd “fluidization”) will artificially differentiate between end-\nmember behaviors that are in fcat a part of a continuum and present similar field\nmenifistations (ground settlement, wave attenuation, rigidity loss, high fluid\npressure). Such an artificial division may only decrease the ability to understand\nand analyze liquefaction and pressurization processes as part of a continuum, as\nrevealed by our new work.\nFurthermore, the term ‘fluidization’, which usually refers to grain suspension via an\nexternally imposed fluid flow, could be misleading in this context, since there is no\nimposed flow through the boundary (e.g., fluid injection at the bottom), rather the\nflow is the outcome of in-situ compaction-pressurization feedback across the\ncompaction front.\nTo clarify the terminology confusion, we included a new appendix (A) that\ndiscusses terminology issues (also relevant to the following comment).\n2. In the same vein, it may be better to distinguish clearly between the ‘drained’ and ‘undrained’\nterms. In liquefaction research (see Refs 22, 24 and 29) the term ‘partially drained’ was\nintroduced to demarcate this. In this paper the authors used the Debroah number (De) to\nindicate partially drained events.\nThe current study does not refer to \"drained\" and \"undrained\" in the context\nof pore pressure value, but rather by comparing the time scale of pore pressure\ndiffusion and the time scale of the forcing inducing pore space deformation.\nThe perspective we attempt to bring forward as part of this work is that the non-\ndimensional De number can be used to adequately define the drainage conditions.\nThe De number and the pore fluid equation (eq. 1) describe two end-members. When\nDe >> 1 the layer response is undrained, fluid cannot dissipate in the timescale of\ngranular deformation and the diffusion term is negligible. The other end-member,\nwhere De << 1, fluid can easily flow within the timescale of pore space change, and\nthe first term of eq. 2 is negligible, is referred to as drained.\nWe believe that previous uses of \"partially drained\" to describe scenarios of De <<\n1 reflect the degree to which the undrained hypothesis has become an axiom in the\nliquefaction research community. A prefix “partially” was felt needed to explain\nobservations of fluid flow during liquefaction. We suggest the contrary: Two\ntheoretical end members should be considered, “drained” and “undrained”. Once\na clear definition of them is agreed, it is easier to discuss the in-between cases\n(which are probably common in nature and characterized by De ~ 1), and they\nshould be referred to as \"partially drained\" or \"partially undrained\".\nTo extend our point, we note that Darcy’s law describes fluid flow (drainage) in\nrelation to non-hydrostatic pressure gradients. I.e., one cannot categorically\nexclude excess pore pressure in drained conditions.\nFollowing this comment, we added a discussion of terminology related issues in\nnew appendix (A).\n3. The authors may want to clarify what are the time scales for the liquefaction and fluidisation\nevents. One would expect liquefaction to happen (for loose sands) in the first few cycles of\nearthquake shaking. Fluidisation can happen and continue into the post seismic period.\nPlease refer to our reply to your question 1, and reply 6 to reviewer #2. We\nmodified Appendix B - Liquefaction initiation timescale to discuss also the time\nscale of the downward liquefaction front.\nGenerally, the time scale of liquefaction initiation is longer under undrained\nconditions (De>>1) and shorter under drained conditions (De<<1). This was\npreviously shown (independent from this study) by Adamidis & Anastasopulis\n(2022).\nWithin the parameter range of our tests, our drained simulations show (to be\npublished as part of a work in progress about repeated liquefaction) that the amount\nof post-seismic compaction is negligible compared to the co-seismic\ncompaction/settlement. Consequently, the timescale of pore pressure dissipation\nis also short (given the characteristic permeability). Based on centrifuge tests, a\nsimilar observation regarding the post-seismic settlement rate was previously\nmade by (Madabhushi & Haigh, 2012; Adamidis & Madabhushi, 2016).\n4. The concept of ‘liquefaction front’ moving from soil surface towards the bedrock as the\nshaking starts and the ‘densification front’ moving from the bedrock towards the soil surface\nhave been previously identified by many researchers (e.g. ref 33). There is a large body of\ndynamic centrifuge test data to show these fronts experimentally occurring in layers of ~ 10m\ndeep. The authors may want to include some of these references into the manuscript.\nWe agree with this comment. Please see our reply to reviewer #2 1st question and\nthe related text changes we introduced following these comments.\n5. There is considerable body of research on liquefaction that establishes cyclic shear strain\namplitude threshold i.e. earthquakes that cause shear strains above this threshold cause\nexcess pore pressure generation while those below won’t. I feel this is more physical as shear\nstrains are physically linked to the movement of individual soil grains and the length parameter.\nIncidentally I agree with the authors argument about length scale in Sec. 3.1. This almost\ncorresponds to Tatsuoka’s observations of formation of shear bands and their thickness in\nsands (~20 particle diameters). However, I do not see why the energy argument is used in this\npaper. Measures like Arias intensity that depend on the PGV are global measures, where as\nliquefaction (in my view) depends on the relative displacement between particles. I am not\ntalking about fluidisation here (as distinguished in point 1 above). The authors may want to\nargue why PGV is a good measure for liquefaction or otherwise.\nWe agree that liquefaction needs volumetric strain to occur. The relative motion of\ngrains beyond elastic limits is the response of the soil to external excitation, which\nin this paper, we found to be related to the seismic power and the seismic energy\ndensity (scales with 𝑷𝑮𝑽𝟐) This is in agreement with Goren et al. (2010,2011) who\nshowed that unlike undrained conditions, where the volumetric strain control\nliquefaction, the drained response depends on the strain rate (which scales with\nPGV).\nFurthermore, an energy-based approach for the evaluation of liquefaction potential\nis not new, and pursuing it as an alternative to stress-based and strain-based\napproaches was recommended by experts’ workshops and reviews (e.g., Youd et\nal. 2001, NAS 2016).\nSpecifically, the energy-based measure Arias Intensity was previously shown to\ncorrelate with liquefaction potential (Kayen & Mitchell, 1997). We believe that it is\nlikely a good predictor for the undrained end-member. We found that for the drained\nend-member studied here, the best controlling parameter is similar to Arias\nIntensity per cycle when considering a harmonic signal (The seismic power). This\nview in now communicated in the revised version (section 4.1 line 242):\n“A leading energy-based approach for evaluating soil liquefaction potential\nuses the earthquake's Arias intensity (Kayen & Mitchell, 1997). While the Arias\nintensity is a cumulative measure that accounts for the amplitude and frequency\ncontent throughout the duration of the earthquake, the rate of seismic energy\ndensity input (𝒆/𝑻), which we consider here, can be interpreted as a quasi-\ninstantaneous Arias intensity or an average power of ground shaking over one\nshear cycle. The excellent performance of the rate of seismic energy density input\nin explaining the amount of compaction (Fig. 3b) and the front velocity (Fig. 3d\ninset), and in defining the clearest threshold between liquefied and non-liquefied\nsimulations (Fig. 3b) suggests that a dynamic process as drained liquefaction\ndepends on the momentary power rather than on the cumulative power. This is\nlikely in contrast to undrained liquefaction, which is a cumulative process in nature\n(the volumetric strain required to initiate liquefaction is accumulated over many\nshear cycles (de Alba et al.,1976), hence it might depend on a cumulative energy\nmeasure like Arias intensity.”\n6. In Sec. 4.1 the authors argue that the time it takes to reach ‘full drained liquefaction’ is\nexceedingly small i.e. 10-3 ~ 10-8 seconds. I presume that the authors mean that this is\ndifferent from ‘seismic liquefaction’ and are saying fluidization can occur much more quickly.\nFirstly, for seismic liquefaction, the excess pore pressure generation can take anything from 1 to\n2 cycles to many cycles (~8 to 10 cycles) depending on the relative density of the sand. Loose\nsands liquefy more easily than dense sands. Secondly, the time for fluidisation surely depends\non the permeability, porosity of the soil strata and pressure gradient applied which will yield\nmuch longer time scales, as pore fluid has to migrate from one region to the other?\nFollowing this comment, we edited the mentioned paragraph as follows:\n“The timescale associated with pore pressure rise to lithostatic\nvalues in the simulations and experiments is short, and likely related to a\nrapid downward propagating liquefaction front. Such a behavior, viewed as\n\"unloading\" fronts (Florin & Ivanov, 1961), was previously identified in\nexperiments (Florin & Ivanov, 1961; Sharp & Dobry, 2002; Gonzalez et al.\n2005) (see Methods section 10.4 and section A in the supplementary\nmaterial). We observed that the deepest location to which the down-going\nliquefaction front reaches correlates with the imposed shaking frequency\n(Appendix C), and that the liquefaction front reaches this deepest location in\nless than two shear cycles, consistent with previous experiments conducted\nunder drained conditions (Wang et al., 2013; Adamidis & Madabhushi,\n2016).\nAnother timescale operating in the system is the time required for an\nisolated, fully immersed, grain (equation 5, Appendix B) to reach its terminal\ndownward velocity. However, since this timescale is exceedingly small,\n𝟏𝟎−𝟑−𝟏𝟎−𝟖 seconds, the acceleration of a single grain is not a rate-limiting\nprocess for triggering drained liquefaction. Recent cyclic triaxial\nexperiments (Adamidis & Anastasopulis, 2022) found that the number of\ncycles required to initiate liquefaction under drained conditions is smaller\nthan under undrained conditions, supporting the hypothesis that a\npressurization time of the order of a few cycles could be indicative of drained\nliquefaction initiation. Such a consideration might apply to a recent ground\nmotion analysis showing that, in some cases, the time for liquefaction\ntriggering is as short as ∼𝟏.𝟕 seconds from the onset of recorded shaking\n(Greenfield, 2017; Ozener et al., 2020).”\nWe stress that in our view, the term “fluidization” does not adequately describe our\nresults, since no pressure gradient is imposed externally. Although we did not\nexplore a wide range of permeabilities in the current study (the simulations and\nexperiments had different permeabilities), a previously published study (Ben-Zeev\net al., 2020) found similar dynamics over a wide range of permeabilities.\n7. I have several objections on the experimental program reported in Sec. 10. Also not sure why\nthe authors do not use data from dynamic centrifuge tests on liquefaction.\nThe main purpose of the experiments was to examine visually the grain dynamics\nand fluid pressurization through the experiments. To achieve that, we chose to use\na completely transparent box to allow filming at least one face of the box. This goal\nwas achieved as we were able to validate visually the formation of the compaction\nfront and its velocity in relation with the imposed power and energy. Despite the\nlimitations of this experimental setup (small size, boundary effects, rigid), it allowed\nus to validate the trends that emerged from theory and simulations in a way that\ncentrifuge tests or laminar box setup would not allow. The experimental results\nagree very well with simulations (the formation of the compaction front and its\nvelocity relation with the imposed energy and power, see figs.2a-d & 3c). The fact\nthat the same dynamics emerged despite the difference in setup between\nsimulations and experiments, and the fact that both experiments and simualtions\nagree with theory, shows that the compaction front dynamics are robust, and the\nthe experimental setup, in spite of its limitations, suffices to capture it.\nTo convey this massage, we added in the discussion (section 4.1) the following\nsentence:\n“Both robust correlations emerged despite inherent differences in the\nboundary conditions and geometrical setup between the simulations and\nexperiments, suggesting that, within the framework of drained liquefaction\ntriggering, the seismic energy density, and possibly a new measure, the\nseismic energy density rate, can be considered as controlling parameters on\nthe magnitude and duration of liquefaction events.”\na) The experiments are in a very small-scale model container 120 mm (cube). The boundary\neffects will be significant in such a small model container.\nWe calculate velocities in the middle part of the box’s face. We note that despite\nthe box’s small dimensions, the same dynamics emerged as in simulations, and\nboth experiments and simulations show remarkable agreement with our\nliquefaction theory, suggesting the experimental setup is indeed sufficient to\nproduce the liquefaction process.\nFollowing this comment, we added information in section 10.3:\n“The instantaneous grain velocity field is measured in every frame by\nusing PIVlab (Thielicke, & Stamhuis 2014; Thielicke & Sonntag, 2021), which\nrelies on sub-frame correlation between timely-adjacent frames. The vertical\nvelocity is then averaged over sub-layers, yielding the vertical velocity of\ngrains as a function of depth and time (see section 10.4). To minimize\nboundary effects from the box's walls, the averaging is done only close to\nthe box's center (approximately in the middle 2/4 of the box's total width).“\nb) What is the relative density of the soil specimens? This should be included in Table 2.\nPorosity is not convenient for this as RD is measure of the soil state with respect to the loosest\nand densest possible states.\nWe are aware of the costume to report the relative density of the soil specimens in\nthe geotechnical community and consequently we would have like to follow the\nreviewer recommendation. However, this very study (see fig. 3a) has revealed to us\nthat the densest state of a soil specimen is a function of the shaking energy and the\nrate of seismic energy input (see fig. 3b), and therefore, the relative density does\nnot have a unique value. On the other hand, porosity is a uniquely defined physical\nquantity, that appears in the theory we developed (from basic physics).\nc) The frequency of shaking is very high (10 Hz, most earthquakes are in the 1 to 5 Hz range)\ncertainly not the maximum frequency, see for example (Souriau) 2006 and (Tosi et\nal. (2012). In order to experimentally explore seismic energy density values in a\nrange similar in magnitude to the seismic energy density in simulations, and to\nensure the low seismic energy density range representative to the far field, we had\nto choose a high frequency for the experiments to overcome limitations on the\nmaximum allowed displacement of our shaker. The agreement of our results with\ntheory reassured that the high frequency did not effect the physics considerably.\nYet, following this comment we did indeed probe more deeply the effect of\nfrequency on our results. We found that one effect of frequency is on the initiation\ndepth of the compaction front. We thus added a new Appendix C showing the\ninitiation depth as function of frequency.\nd) The amplitude of shaking is very small 0.8 mm is the largest amplitude used. How does this\ncompare to the mean particle size (D20 or D50) for this sand?\nthe mean grain diameter in the experiments (reported in Table 1) is 0.02 cm and\nthe minimum shaking amplitude is 0.024 cm (reported in Table 3), i.e., the shaking\namplitude is at least of the order of the mean grain size.\nFollowing this comment, we expanded the method section (10.3) to address this\npoint:\n“Horizontal shaking is applied for 30 seconds with displacement\namplitude of at least one mean grain diameter.”\ne) Referring to Fig.2a, the excess pore pressures do not show any suction cycles, cyclic\noscillations (to indicate the crossing of phase transformation line (following Ishihara) nor reach\nfull liquefaction of 1. The authors may want to explain why.\nThe pressure signal is smoothed in order to highlight the most important\ntrend (and it is a very noisy signal). Furthermore, we see in others non-jacketed\ntests (e.g., centrifuge) that suction cycles are not observed (for example in Wang\net al. (2013) and Adamidis & Madabhushi (2016)\nThe fact that the pressure ratio does not reach 1 can be explained by a number of\nreasons. First, the errorbar (pink shaded area in fig.2a) is calculated based on a\nstatic calibration procedure during the pre-shaking stage. It is possible that the\nerrorbars should be larger during the dynamic stage of shaking. I.e., the error\nshown is the minimal one.\nAlternatively, the pressure in fig. 2a is normalized by the initial static solid stress\nwhich is evaluated analytically from porosity, grain density and depth below the\nsurface. The pressure sensors remain in place during the experiment (they are\nattached to the box’s face), but since during drained liquefaction the grains settle,\nthe overburden at the sensor depth is expected to decrease (relative to its initial\nvalue). This will bring the curve closer to 1 (dividing the pressure by a smaller\nnumber). Another possible effect is that the sidewalls (which are absent in\nsimulations) carry some constant amount of the soil load, instead of the fluid.\nTogether with the above suggestions, the fact that the pressure value reaches a\nplateau before dropping back to zero might indicate that the pore pressure ratio\nreached the value of 1 during the liquefied state. A stable pressure value is easier\nto maintain at a complete liquefaction. Despite these suggestions, to remain\nconservative, we prefer to keep the current form of the pressure curve.\nf) The time for generation of excess pore pressure in Fig. 2a is much more slower than in the\nsimulations and magnitudes are smaller.\nThe time axis in figure 2a is normalized by the duration of the liquefaction event,\nwhich is not identical in the simulation and the experiment presented. Please see\nthe figure below, showing that the number of cycles until complete liquefaction is\nN=2 in the simulation and N=7 in the experiment. The time until complete\nliquefaction is 0.34 and 0.65 seconds, respectively. This difference is attributed to\nthe delay of the shaker before reaching the designated amplitude (due to inertia and\nfriction), unlike in simulations where the imposed shaking BC is immediately\napplied.\ng) In fig. 2d, the shear strains are biased towards the positive side? Why?\nThis is due to a complication arising when defining shear strain in an unconfined\nlayer. If the layer is confined between two sheared plates, then it is easy to define\nthe shear strain in the layer as the relative motion between the two plates (where\nthe shear is imposed) divided by the layer’s thickness. Here, with the free surface\nat the top, we define the shear strain in each sub-layer as the difference of the\nhorizontal velocity of grains integrated over time, divided by the vertical thickness\nof that sub-layer. The data presented in figure 2d is from the middle of the\nsimulation box. The strain axis being not symmetric around zero is a remnant of the\ninitial effect where the grains were moving in the positive direction before\nliquefying. Also, the imposed motion is of (1-cos) which is not centered around zero.\nh) In Fig.2d, are the experimental results and simulations at the same frequency? The\nrecoupling and increase in oscillations in the experiment are not clear.\nWe assume you refer to fig 2c. The simulation is of ~6.7 Hz and experiment of 10\nHz. The time axis is normalized by the total duration of the liquefaction event, which\nis different between the presented simulation and experiment. Indeed, the change\nin amplitude of motion of the experiment grains is less dramatic than in the\nsimulations. This is likely because in the experiments, the sidewalls impose\nhorizontal motion even to the liquefied layer, unlike simulations where the motion\nis imposed only at the bottom. Despite this behavior, there is a clear trend in the\nexperiments relating to regaining 60 % of the imposed shear amplitude. Following\nthis comment, to achieve clarity, we added a magnification inset around the front in\nthe experiment panel of fig.2c."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have meaningfully addressed all comments from the reviewing process. The\nresubmitted manuscript is improved as a result and I would support its publication.", "reviewer_3": "The authors have addressed many of the comments I have raised. I would still standby my\ncomments that it would have perhaps been better to compare the authors theory against dynamic\ncentrifuge test data with PIV analyses, rather than their small scale experiments. Not withstanding\nthat, I am happy with the changes they have made.", "reviewer_4": "This paper deals with the important topic of liquefaction under drained conditions. Unfortunately,\nits contribution seems unclear and untimely. It is not clear from the manuscript what exactly is the\ncontribution here. For example, the authors claim:\n\"The discovery that high permeability, well-drained soils should not be apriori assumed liquefaction\nresistant\" That is not a \"discovery\" since that is well known in the literature. There have been at\nleast one decade worth of studies on static liquefaction and drained instabilities showing that the\nstatement above is correct.\nAlso, the authors rely on a coupled solid-fluid simulation that is based on spherical DEM. On the\nother hand, the experiments are done using a transparent tank filled with sand and water. The\nability to compare between simulations and experiments here is at best qualitative. Of course, a\nloose packing of spheres or sand under unconfined conditions can display instabilities, especially\ntowards the top of the column where effective stresses approach zero anyway. Any excitation that\nwould increase excess pore pressure would tend to provoke instabilities. Again, this is not deep or\nnew.\nUnfortunately this study as is cannot be accepted for publication given its inability to clearly\ncommunicate its contribution and its timelines."}, "author_response": "Dear reviewers,\nWe wish to thank you for putting the time and effort in reviewing our manuscript.\nThe manuscript’s main objective is to solve the formerly puzzling observations of liquefaction in\nthe earthquake’s intermediate and far fields, under low seismic energy density. To achieve this\nobjective, we harness a mechanism for liquefaction under drained conditions, which is not limited\nby the same seismic energy threshold as undrained liquefaction. We further show that the unique\ncharacteristics of drained liquefaction, i.e., co-seismic heterogeneous compaction controlled by\nthe compaction front model (whose magnitude and velocity are related to the seismic excitation\nmagnitude), can be invoked to explain the decay in the number of liquefaction observations with\ndistance from earthquakes epicenters.\nTo address the main issue raised in the last round of review, the relatively small size of our\nexperimental setup, we wrote a new appendix that compares our findings to previously published\nlarge-scale simulations and experiments. We also rewrote the entire conclusions section to better\nemphasize the “take home massage” of the manuscript.\nSincerely,\nShahar Ben-Zeev on behalf of the co-authors\n\nPoint by point response:\nReviewer #2 (Remarks to the Author):\nThe authors have meaningfully addressed all comments from the reviewing process. The\nresubmitted manuscript is improved as a result and I would support its publication.\nWe are happy that our revisions were found to be satisfactory for publication, and\nthank Rev#2 for meaningful comments.\nReviewer #3 (Remarks to the Author):\nThe authors have addressed many of the comments I have raised. I would still standby my\ncomments that it would have perhaps been better to compare the authors theory against\ndynamic centrifuge test data with PIV analyses, rather than their small scale experiments. Not\nwithstanding that, I am happy with the changes they have made.\nWe are happy that our revisions were found to be satisfactory for publication, and\nthank Rev#3 for meaningful comments.\nRegarding the comment that it is desirable to compare our results to large scale\ncentrifuge experiments. Our physics-based theory, which is scale-less and fits well\nour experimental and computer modeling results, serves to validate the\napplicability of our results to all scales, including large scale. Yet the reviewers’\nrequest to compare our theory and simulations also to large scale experiments is\nunderstandable. To this end we have searched previously published experimental\nliterature, and were able to identify two manuscripts reporting large-scale\ncentrifuge experiments (but without PIV) and simulations under enhanced gravity\n(simulating greater stresses and depths). These previous studies reported\nobservables that allow us to partly compare the dynamics they observed to our\ntheory and observations.\nWe added a new appendix E (called “High gravity experiments and simulations\nsupport depth independent triggering of drained liquefaction”), in which we\ncompare different aspects of our predictions and results to the centrifuge\nexperiment (50 g) of Adamidis & Madabhushi (2018) and to the high gravity\nsimulations (30 g) of El Shamy & Zeghal (2007). For the latter, we show the existence\nof a co-seismic compaction front (triggered at prototype depth of 5.2 meters) with a\nfront propagation velocity that fits very well our prediction for the relation between\nfront velocity and both the imposed seismic energy density and rate of seismic\nenergy density. Overall, these two studies support depth-independence (and\nconsequently stress-independence) of the model proposed in our manuscript, and\nconfirm the applicability of our model, experiments, and theory to the large-scale.\nReviewer #4 (Remarks to the Author):\n1) This paper deals with the important topic of liquefaction under drained conditions.\nUnfortunately, its contribution seems unclear and untimely. It is not clear from the manuscript\nwhat exactly is the contribution here. For example, the authors claim:\n\"The discovery that high permeability, well-drained soils should not be apriori assumed\nliquefaction resistant\" That is not a \"discovery\" since that is well known in the literature. There\nhave been at least one decade worth of studies on static liquefaction and drained instabilities\nshowing that the statement above is correct.\nWe thank reviewer #4 for his review and comments.\nWe agree that the above phrasing (from the Conclusions Section) was misleading,\nand should not have appeared in the text. The manuscript actually is not centered\naround the “discovery” of drained liquefaction: Already in the introduction, we cite\nstudies that discuss the possibility of drained liquefaction (e.g., Madabhushi &\nHaigh, 2012; Lakeland et al., 2014; Adamidis & Madabhushi, 2018; Adamidis et al.,\n2022), as well as studies that we ourselves have published starting from 2011 on\ndrained liquefaction. The reviewer’s comment made us realize that the entire\nconclusions section should be improved, and we rewrote it entirely, to better reflect\nthe content of the paper and its title. The main take home message of this\nmanuscript is not the mere existence of drained liquefaction, but its ability to\nexplain previously puzzling far-field liquefaction events.\n2) Also, the authors rely on a coupled solid-fluid simulation that is based on spherical DEM. On\nthe other hand, the experiments are done using a transparent tank filled with sand and water.\nThe ability to compare between simulations and experiments here is at best qualitative.\nHundreds of previously published papers quantitatively compare physical granular\nexperiments (with and without water) to DEMs, and we fail to see the problem in\nthat. We explicitly make a quantitative comparison between simulations and\nexperiments in Fig 3d, showing the seismic energy density range over which the\nsimulations and experiments were conducted (far-field energies) and their general\nagreement regarding the dependency between the compaction front velocity and\nthe imposed seismic energy density. The general quantitative agreement between\nsimulations and experiments, and also the excellent fit between these 2 scientific\nmethodologies and our physics-based theory, demonstrate that our findings are\nrobust, and can be identified even if the experimental setup is not identical in all\nminute details to the simulation’s setup. In fact, the agreement between simulations\nand experiments performed with different materials and boundary conditions points to the\nuniversality of the process - a clear advantage and strength of the model.\n3) Of course, a loose packing of spheres or sand under unconfined conditions can display\ninstabilities, especially towards the top of the column where effective stresses approach zero\nanymore. Any excitation that would increase excess pore pressure would tend to provoke\ninstabilities. Again, this is not deep or new.\nThis comment suggests that our observed liquefaction is due to the low stresses\nthat characterize small-scale experiments. First, we note that we demonstrated in\nthe manuscript the four main seismically induced liquefaction indicators in the field\nand lab. The liquefaction starts at the base of the grain column, ~1 m below the\nsurface in the simulations (not so shallow). Second, we note again that our physics-\nbased theory predicting drained liquefaction is scale-less, and thus does not\nrequire or assume small stresses. Finally, as detailed in our answer to Reviewer 3,\nfollowing this concern we searched the literature and were able to identify two\npreviously published manuscripts that conducted similar experiments (but without\nPIV) and simulations under enhanced gravity. These studies reported some\nobservables that allow us to partially compare the dynamics they observed to our\ntheory and observations.\nWe added a new appendix E, called “High gravity experiments and simulations\nsupport depth independent triggering of drained liquefaction”, in which we\ncompare different aspects of our predictions and results to the centrifuge\nexperiment (50 g) of Adamidis & Madabhushi (2018) and to high gravity simulations\n(30 g) of El Shamy & Zeghal (2007). For the latter, we show the existence of a co-\nseismic compaction front (triggered at prototype depth of 5.2 meters) with a velocity\nthat fits very well our prediction for the relation between front velocity and the\nimposed seismic energy density and the input rate of seismic energy density.\nOverall, these two studies indicate that the model proposed in our manuscript\noperates also at greater soil depths and at high stress levels.\n4) Unfortunately this study as is cannot be accepted for publication given its inability to clearly\ncommunicate its contribution and its timelines.\nThe processes that lead to soil liquefaction are not yet completely understood,\nwhile liquefaction risk is evident and expected to increase (due to current and\nplanned massive sea reclamation projects around the world). In this manuscript we\nattempted to provide an answer to how soil liquefaction is initiated at the far-field,\ndespite a very low seismic energy density excitation (Wang & Manga, 2021). We\nshowed that if liquefaction occurs as a drained phenomenon, it possesses unique\ncharacteristics that provide a possible answer to this question. Hence, we think that\nthis study is timely and can impact liquefaction studies as well as hazard\nassessment practices. In addition, following this comment, we rewrote the entire\nConclusions Section, to better communicate our contributions."}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have carried out comparison with previously published high gravity centrifuge test\ndata and included in the Appendix E. I am happy with the changes they have made.", "reviewer_4": "The authors have significantly improved the manuscript, including its main contribution. However,\nthere is one item that still needs to be solved. It has been made worse. The authors \"fail to see\nwhat the problem is\" to use classic spherical DEM to compare with experiments using sands.\nThe effect of particle shape has been clearly demonstrated. DEM has remained qualitative until very\nrecently. At the very least, the authors must acknowledge this explicitly in the manuscript. Current\nspherical simulations are qualitative and require tuning of ad-hoc parameters like damping and\nrotational resistance to capture the effect of particle shape. This has been documented in recent\npublications such as Kawamoto et al. 2018 (All you need is shape. JMPS).\nOnce this issue is resolved by inserting the appropriate remark in the manuscript, the paper\ncould be accepted for publication."}, "author_response": "Point-by-point response to the reviewers’ comments:\nManuscript title: “Drainage explains soil liquefaction beyond the earthquake near-field”.\nReviewer #3 (Remarks to the Author):\nThe authors have carried out comparison with previously published high gravity centrifuge test\ndata and included in the Appendix E. I am happy with the changes they have made.\nWe thank you for helping improve our manuscript.\nReviewer #4 (Remarks to the Author):\nThe authors have significantly improved the manuscript, including its main contribution.\nHowever, there is one item that still needs to be solved. It has been made worse. The authors\n\"fail to see what the problem is\" to use classic spherical DEM to compare with experiments\nusing sands. The effect of particle shape has been clearly demonstrated. DEM has remained\nqualitative until very recently. At the very least, the authors must acknowledge this explicitly in\nthe manuscript. Current spherical simulations are qualitative and require tuning of ad-hoc\nparameters like damping and rotational resistance to capture the effect of particle shape. This\nhas been documented in recent publications such as Kawamoto et al. 2018 (All you need is\nshape. JMPS).\nOnce this issue is resolved by inserting the appropriate remark in the manuscript, the paper\ncould be accepted for publication.\nWe thank you for helping improve our manuscript. We edited the following\nparagraph in section 2.5:\n“The numerically-identified low liquefaction triggering threshold might be partially\nrelated to model parameters, including the use of perfectly smooth and spherical grains.\nMore realistic grain shapes could potentially necessitate a higher triggering threshold61.\nNevertheless, our experiments with natural sand grains showed liquefaction triggering at\nPGA values as low as Aω/g≈0.1 (the precise threshold was not investigated). This\nsuggests that the intrinsic characteristics of natural grains are unlikely to alter the\noutcomes significantly, permitting low PGA liquefaction triggering under real-world\nconditions. Furthermore, after the triggering phase, once granular contacts are minimized,\nthe dynamics of the liquefied layer are expected to be independent of grain shape, as\nevidenced by the congruence of front velocities across simulations, experimental results,\nand theory (depicted in Fig. 3d).”"}]

Decision letter