# Static liquefaction: The role of grain size polydispersity from a micro-structural perspective

Carolina Castro-Malaver, Manuel Cárdenas-Barrantes, David Cantor, Mathieu Renouf, Emilien Azéma, Carlos Ovalle, Arcesio Lizcano

► To cite this version:

Carolina Castro-Malaver, Manuel Cárdenas-Barrantes, David Cantor, Mathieu Renouf, Emilien Azéma, et al.. Static liquefaction: The role of grain size polydispersity from a micro-structural perspective. Computers and Geotechnics, 2026, 193, pp.107996. 10.1016/j.compgeo.2026.107996. hal-05519553

HAL Id: hal-05519553

https://hal.science/hal-05519553v1

Submitted on 19 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Research paper

# Static liquefaction: The role of grain size polydispersity from a micro-structural perspective

Carolina Castro-Malaver<sup>a,b,c ID</sup>, Manuel Cárdenas-Barrantes<sup>a</sup>, David Cantor<sup>d ID</sup>, Mathieu Renouf<sup>c</sup>, Emilien Azéma<sup>a,c,e ID</sup>, Carlos Ovalle<sup>a,b ID</sup>,<sup>*</sup> Arcesio Lizcano<sup>f ID</sup>

## ARTICLE INFO

**Keywords:**

Static liquefaction

Grading

DEM

Micro-structure

## ABSTRACT

Static liquefaction, defined as the loss of strength after strain softening under undrained shear, is often the cause of catastrophic failures in loose earthfills, waste dumps and tailings storage facilities. The number of failures does not cease to increase yearly, revealing operational faults, construction defects, and a lack of knowledge regarding liquefaction triggering. While liquefaction is well understood at the scale of a representative elementary volume of soil, the triggering mechanisms at the particle scale are seldom studied. This study aims to analyze the links between the micro- and macro-mechanisms controlling static liquefaction by means of numerical simulations using the Discrete Element Method. We consider mono to highly size polydisperse 2D samples of discs, prepared at varied densities from loose to dense packings. As expected, loose samples liquefy, while dense ones continue to exhibit resistance under undrained shearing. Some medium-dense cases liquefy only temporarily, recovering their shear strength at larger strains. We reveal a dual mechanism for liquefaction through multi-scale analysis: a collapse of the contact network—marked by the coordination number dropping below 3—and the emergence of low-density regions. Temporary liquefaction involves transient connectivity loss with minor density fluctuations, enabling stress recovery. In contrast, full liquefaction combines both effects, leading to an irreversibly disconnected and heterogeneous microstructure.

## 1. Introduction

Waste storage structures—such as tailings dams and waste dumps— are particularly vulnerable to static liquefaction. This phenomenon refers to a sudden loss of strength in loose, saturated soils subjected to undrained (constant volume) monotonic loading. As shearing progresses, pore water pressure builds up, reducing effective stress to near zero and resulting in flow-like deformation. Numerous examples of failures due to liquefaction in tailings storage facilities have been documented throughout the 20th century (Dobry and Álvarez, 1967; Fourie et al., 2001; Verdugo and González, 2015; Santamarina et al., 2019) and continue to occur in recent years (Morgenstern et al., 2015, 2016; Robertson et al., 2019). Similarly, several failures have been reported in dumps composed of discarded construction materials (Zhao et al., 2024), municipal solid waste (Blight and Fourie, 2005), and coal mine waste (Dawson et al., 1998).

Significant research efforts have been dedicated to understand the mechanisms of static liquefaction in loose granular soils (Castro, 1969; Verdugo and Ishihara, 1996; Lade and Yamamuro, 2011). It is well known that liquefaction susceptibility decreases with increasing the state parameter (Been and Jefferies, 1985) (i.e., decreasing the initial relative density and/or increasing the consolidation mean stress), and is sensitive to the granular fabric (Li and Dafalias, 2012; Lü et al., 2017). Moreover, empirical evidence shows that intrinsic material parameters, such as particle size distribution (PSD), could significantly influence the liquefaction susceptibility (Troncoso and Verdugo, 1985; Karim and Alam, 2014; Kwa and Airey, 2017; Pires-Sturm and DeJong, 2022). Undrained strength at a given state parameter increases with particle size polydispersity, as a result of enhanced dilative tendencies (Ahmed et al., 2023). For instance, Yamamuro and Lade (1997) studied the

<footnote><sup>a</sup> Department of Civil, Geological and Mining Engineering, Polytechnique Montréal, 2500 Chem. de Polytechnique, Montréal, H3T 0A3, Québec, Canada</footnote>

<footnote><sup>b</sup> Research Institute on Mines and Environment (RIME) UQAT-Polytechnique, 2500 Chem. de Polytechnique, Montréal, H3T 0A3, Québec, Canada</footnote>

<footnote><sup>c</sup> LMGC, Université de Montpellier, 860 Rue de St - Priest, Montpellier, 34090, France</footnote>

<footnote><sup>d</sup> Department of Civil and Environmental Engineering, University of Washington, 3760 E. Stevens Way NE, Seattle, 98195, WA, USA</footnote>

<footnote><sup>e</sup> Institut Universitaire de France, 1 Rue Descartes, Paris, 75005, France</footnote>

<footnote><sup>f</sup> SRK Consulting, 320 Granville, Vancouver, V6C 1S9, British Columbia, Canada</footnote>

<footnote>* Corresponding author at: Department of Civil, Geological and Mining Engineering, Polytechnique Montréal, 2500 Chem. de Polytechnique, Montréal, H3T 0A3, Québec, Canada.
E-mail address: carlos.ovalle@polymtl.ca (C. Ovalle).</footnote>
influence of non-plastic fine content (FC) on fabric and undrained strength of loose silty sands. They conceptually proposed that, at low FC, fines partially occupy the voids between sandy grains without significantly affecting the load-bearing granular structure; i.e., mainly acting as floating particles that do not bear interparticle forces. This condition is called underfilled and is characterized by a soil skeleton primarily formed by the coarse fraction of the soil (Thevanayagam et al., 2002). On the other hand, at high FC the fines generate an overfilled condition, in which they dominate the force chain network and carry most of the stresses. Between these two extremes lies a transitional state, where granular chain forces include both fines and coarse grains, creating unstable force chains that lead to a compressible fabric highly susceptible to liquefaction. As a result, although the transitional packing corresponds to the densest state, it exhibits the lowest liquefaction resistance (Thevanayagam et al., 2002; Lade et al., 2009; Huang et al., 2023). Beyond this threshold, further addition of fines promotes a more stable overfilled fabric, thereby increasing the liquefaction resistance (Sadrekarimi, 2014). However, while recent empirical studies exist (Rahmani and Naeini, 2020; Gobbi et al., 2022; Yazdani et al., 2022), there are not measurable microstructural parameters to support this conceptual framework. In this context, a gap persists in establishing robust connections between micro-scale mechanisms and macro-scale responses associated with static liquefaction in grain size polydisperse materials. Advancing our understanding of this phenomenon could enhance the reliability of predictive models. However, this requires a comprehensive multi-scale approach that captures the fundamental processes governing material behavior across scales.

To the best of the authors' knowledge, only a limited number of particle-scale physical experiments on liquefiable soils have been reported, likely due to the inherent complexity such tests would involve. For instance, Sukhumkitcharoen et al. (2024) performed biaxial tests and applied image-based analysis to track the evolution of local density and void anisotropy in liquefiable 2D granular samples. They observed a strongly heterogeneous evolution of these quantities, with localized changes preceding macroscopic instability. However, the void-based anisotropy measures remain primarily descriptive and cannot be directly related to force transmission or contact-network organization. To address this gap, several researchers have proposed numerical simulations based on the Discrete Element Method (DEM). In DEM, undrained conditions can be modeled either by enforcing constant volume or by coupling the solid phase with fluid dynamics. The latter approach, which involves solving fluid flow equations and integrating them with DEM to capture fluid-solid interactions (Polanía et al., 2022; Coppin et al., 2023; Wang et al., 2024; Zeghal and El Shamy, 2004), demands substantial computational resources. Alternatively, when fluid transport processes such as erosion or particle dragging are neglected, constant volume conditions offer a practical and widely accepted approximation for simulating undrained behavior in granular materials (Dyvik et al., 1987; Shafipour and Soroush, 2008; Liu et al., 2015; Otsubo et al., 2022).

DEM has opened a window into the grain-scale mechanisms that drive failure. At the microscopic level, static liquefaction reflects a fundamental breakdown in the internal structure of granular materials (Gong et al., 2012). As shearing continues under monotonic undrained conditions on loose specimens, the network of contacts that transmits forces between particles gradually weakens (Huang et al., 2018). This process can happen abruptly: the number of force-carrying contacts drops, contact chains collapse and the internal fabric loses its ability to resist shear (Cárdenas-Barrantes and Ovalle, 2025). The majority of these DEM-based investigations have focused on cyclic strain-induced liquefaction (Soroush and Ferdowsi, 2011; Zhang and Evans, 2020; Yang et al., 2022; Rahman et al., 2021), emphasizing the role of the granular fabric in the ability to rebuild the force network during cyclic mobility (Evans and Zhang, 2019; Wang et al., 2016; Wang and Wei, 2016; Wei et al., 2018; Yang et al., 2021). Moreover, recent DEM studies that systematically control the initial granular

fabric have demonstrated its significant influence on undrained stress paths (Salimi et al., 2025; Irani et al., 2024), providing microstructural insights into how loading orientation and sample preparation methods affect liquefaction susceptibility (Yang and Wu, 2017; Zhang et al., 2023; Yang and Taiebat, 2024). For instance, using DEM numerical simulations, Otsubo et al. (2022) demonstrated the important role of inherent fabric anisotropy and local void orientation in controlling both static and cyclic liquefaction resistance. Interestingly, Martin et al. (2020) observed that fabric evolution under monotonic loading differs from that under cyclic loading. In the latter, contact loss is uniformly distributed throughout the specimen, and the main contact network does not collapse upon cycling. In contrast, during monotonic undrained loading of loose packings, the force network breaks down rapidly following the peak deviatoric stress (Gong et al., 2012).

Few micro-structural DEM studies on the specific effect of grading on liquefaction have been reported (Zuo et al., 2023; Basson et al., 2024b; Banerjee et al., 2023). It has been shown that polydispersity increases the undrained strength due to higher particle connectivity, and the force network during constant volume shearing exhibits increased anisotropy of strong contacts, thus improving the undrained strength (Basson et al., 2024a). Yet, important questions remain open. Although expanded grain size polydispersity enhances packing density, it has no influence on the normalized critical shear strength (i.e., internal critical friction angle) (Yang and Luo, 2018; Amirpour Harehdasht et al., 2018; Polanía et al., 2023; Girumugisha et al., 2024). This decoupling between solid fraction and strength raises complementary issues for liquefaction. Whether a broader PSD enhances micro-structural stability by facilitating denser packings, or instead promotes local instabilities due to phenomena such as stress concentration or local density heterogeneity, remains an open question. More fundamentally, the relationship between the evolving contact network and the macroscopic softening that culminates in failure at large strains is still not fully understood. Addressing these issues requires a detailed investigation of the material's micro-structure—a challenge for which DEM may excel.

With the aim of gaining deeper insights into the particle-scale mechanisms responsible for triggering soil liquefaction, this study employs DEM simulations of highly grain size polydisperse assemblies. We simulate undrained shearing conditions by imposing constant volume boundary conditions on samples with varying PSD and initial densities. The density of the specimens is controlled by randomly removing a defined number of floating particles after initial consolidation. As a result, before shearing the specimens share an identical initial fabric among force-carrying grains, while their overall densities vary. We capture the response of loose specimens that exhibit pronounced liquefiable behavior, dense specimens that remain stable and non-liquefiable, and intermediate cases that lie within the transition zone between these two extremes. The results are analyzed from both a macro-mechanical perspective and a micro-mechanical viewpoint, using descriptors such as coordination number, amount of floating particles, and local solid fraction distributions. By bridging the micro- and macro-scales, this approach allows for a more comprehensive understanding of the conditions that govern the onset and evolution of static liquefaction in granular soils.

This paper is organized as follows: Section 2 describes the numerical approach used to prepare 2D granular assemblies of discs at different PSD and initial densities, and to shear them under constant volume at large strain; Section 3 presents the macro-mechanical behavior observed in the numerical tests and identifies the cases of liquefaction, temporary liquefaction and no-liquefaction, depending on both the initial density and PSD; Section 4 presents and discusses the micro-structure and connectivity, the analysis of the evolution of local density, and the bridge between micro- to macro-mechanical descriptors; finally, Section 5 summarizes our findings, conclusions and perspectives.
## 2. Numerical method

### 2.1. Contact dynamics with linear elasticity

To simulate undrained shearing of granular materials under constant volume conditions, a numerical method must allow for internal rearrangements and local contact deformations. In most cases, this is achieved using soft-particle DEM (Cundall and Strack, 1979; Weinhart et al., 2019), where particles are assumed to slightly overlap at contacts. This local deformation is modeled through force-displacement relationships (e.g., linear spring or Hertz-Mindlin contact law), often coupled with damping to account for energy dissipation. Therefore, in the soft-DEM framework the contact force arises from the overlap between particles, and Newton's equations of motion are integrated explicitly with small time steps to resolve the dynamic interactions. Such methods inherently allow for constant volume shearing via contact compliance.

This study employs a different approach to simulate constant volume shear tests, using Contact Dynamics (CD) with regularized contact laws inspired by soft-DEM techniques. In classical CD, particles are treated as perfectly rigid and contact interactions are governed by non-penetration and Coulomb friction constraints, resolved implicitly at each time step $\delta t$. This makes CD particularly efficient and numerically robust for dense granular systems, especially when dealing with highly size polydisperse assemblies, as it allows for significantly larger time steps compared to soft-DEM. The equations of motion are integrated over $\delta t$, and contact interactions are enforced through kinematic constraints formulated as complementarity conditions. In the normal particle contact direction, the Signorini condition ensures unilateral contact, where the normal force $f_n$ and the normal relative velocity $u_n$ must satisfy a non-penetration condition (see Fig. 1b). Coulomb's friction law governs the tangential interaction, where the tangential force $f_t$ must satisfy $|f_t| \le \mu_s f_n$; sliding is activated when $u_t \ne 0$ (see Fig. 1c). Energy dissipation is governed by the friction coefficient $\mu_s$, along with restitution coefficients in the normal and tangential directions. A global solution for contact forces and particle velocities is obtained using a nonlinear, iterative Gauss-Seidel-type algorithm, which can be parallelized efficiently (Renouf et al., 2004). Further details on the algorithm and its mathematical formulation are provided by Radjaï and Dubois (2011).

In constant volume shearing, strains must therefore be accommodated through local rearrangements and slight elastic deformations at contacts. This physical behavior can be incorporated into the CD framework without modifying its fundamental structure. As shown by Krabbenhoft et al. (2012), a linear elastic contact law can be reformulated as a Signorini-type condition via a simple change of variables. Specifically, considering a linear normal force-displacement relation $f_n = k_n \delta$, where $k_n$ is the contact stiffness and $\delta$ the normal overlap between two contacting particles, one can define a shifted gap variable $\tilde{\delta} = \delta + f_n/k_n$. Note that in the present study, only circular grains in a 2D framework are considered. Under this assumption, a linear normal contact law is used and the same stiffness $k_n$ is assigned to all particles (da Cruz et al., 2005). Therefore, the elastic law is mapped into a form compatible with the classical Signorini condition, in this case expressed in terms of $\tilde{\delta}$ and $f_n$. As shown in Fig. 1a, this is equivalent to replacing the normal velocity $u_n$ by a modified velocity $\tilde{u}_n = \tilde{\delta}/\delta t$ in the usual contact graph representations. This reformulation allows elastic compliance to be included directly within the CD algorithm, still satisfying its implicit, non-smooth resolution scheme. The method thus retains its key advantages — robust contact handling, strict volume conservation, and large stable time steps — while being capable of capturing the small elastic deformations necessary to simulate constant volume shearing. See Krabbenhoft et al. (2012) for more details about the numerical implementation.

All simulations presented in this work were performed using the open-source platform LMGC90, a computational framework developed at the University of Montpellier for modeling the dynamics of collections of deformable or rigid particles with various contact laws (Dubois et al., 2011, 2018).

<figure><img src="image_3.png" /><figcaption>Fig. 1. Graphic representation of (a) the standard linear elasticity contact law, (b) the regularized linear elasticity at non-smooth contact dynamics for normal forces and (c) the tangential part.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Particle size distributions (PSDs) used in this study for each size span S. The horizontal axis represents normalized particle diameter ($d/d_{max}$).</figcaption></figure>

### 2.2. Sample preparation

We generated two-dimensional granular samples made of discs with controlled size polydispersity, characterized by the size span parameter S defined as:

$$S = \frac{d_{\max} - d_{\min}}{d_{\max} + d_{\min}}, \qquad (1)$$

where $d_{\max}$ and $d_{\min}$ denote the maximum and minimum particle diameters, respectively. We considered eight values of S ranging from 0.1 (quasi-monodisperse) to 0.8 (highly polydisperse) in increments of 0.1. For each value of S, Fig. 2 shows the particle size distribution (PSD) constructed by dividing the size range [$d_{\min}, d_{\max}$] into $N_c = 12$ discrete classes. The total particle volume was uniformly distributed among these classes, ensuring a controlled range of particle sizes.

Particles were deposited into a square container of dimensions $l_x$ and $l_y$, using a potential-based deposition algorithm (Taboada et al., 2005; Voivret et al., 2007). To ensure that the number of grains exceeds a minimum representative elementary volume (REV) (Quiroz Rojo et al., 2025; Cantor and Ovalle, 2025), the total number of particles increased linearly from 10,000 for S = 0.1 to 18,750 for S = 0.8. Therefore, $l_x$ depends on S to ensure a REV; for instance: $l_x \sim 87d_{max}$ for S = 0.1, and $l_x \sim 36d_{max}$ for S = 0.8.

### 2.3. Consolidation stage

Once the particles were deposited, each sample was uniaxially compressed by applying a vertical pressure $P_0 = 10$ kPa, as schematically presented in Fig. 3a, until the packing fraction $\phi = V_s/V$ stabilized ($V_s$ is the total volume of solid particles and V is the sample volume). A sample was considered stable when the relative variation of $\phi$ between two successive time steps was less than $\Delta\phi = 10^{-4}$. In all simulations, gravity was set to 0 to avoid a stress gradient, and the interparticle
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Schematic of the numerical protocol: (a) consolidation stage with rigid boundaries, and (b) simple shear configuration with periodic boundaries along the horizontal axis—shown by red dashed lines—, white discs are particles glued to the top and bottom walls to generate a rough surface.</figcaption></figure>

<table><caption>**Table 1**
Summary of the main model parameters.</caption><thead><tr><td>Parameter</td><td>Symbol</td><td>Value</td></tr></thead><tbody><tr><td>Maximum particle diameter</td><td>d<sub>max</sub></td><td>1.0 m</td></tr><tr><td>Interparticle friction coefficient</td><td>μ</td><td>0.3</td></tr><tr><td>Particle density</td><td>ρ</td><td>2.7 g/cm<sup>3</sup></td></tr><tr><td>Consolidation pressure</td><td>P<sub>0</sub></td><td>10 kPa</td></tr><tr><td>Normal contact stiffness</td><td>k<sub>n</sub></td><td>10<sup>5</sup> P<sub>0</sub></td></tr><tr><td>Top and bottom walls speed</td><td>v<sub>x</sub></td><td>4 · 10<sup>-3</sup> d<sub>max</sub> s<sup>-1</sup></td></tr></tbody></table>

friction coefficient was set to $μ = 0.3$. The normal contact stiffness was fixed at $k_n = 10^3 P_0$. Table 1 summarizes the parameters used in the DEM model.

To enable large deformations and minimize boundary effects — particularly during shearing —, periodic boundary conditions were imposed along the horizontal direction (i.e., at the vertical walls). Under this condition, any particle leaving one side of the container immediately re-enters from the opposite side, simulating an infinite horizontal domain (see Fig. 3a). Under periodic horizontal boundary conditions and zero gravity, the lateral stress is not prescribed during consolidation and emerges from the internal equilibrium of the contact network. As a result, the initial stress state is not strictly isotropic. This is analogous to simple shear conditions, where only the vertical consolidation stress is controlled. Fig. 4 shows the horizontal to vertical stress ratio ($\sigma_{xx}/\sigma_{yy}$) during the consolidation phase for all samples as a function of the normalized time $t/t_c$, where $t_c$ is the time at the end of this stage. The stress state is nearly isotropic for quasi-monodisperse samples (i.e., $S = 0.0 - 0.5$) and deviates to anisotropic consolidation of $\sigma_{xx}/\sigma_{yy} \sim 70\%$ for $S = 0.7$. All samples are prepared under the same imposed vertical pressure $P_0$, providing a consistent reference state for comparing the subsequent simple shear response across all $S$. Thus, the periodic boundary configurations are particularly well-suited for simple shear simulations, as they enable reaching large strains and critical state conditions.

To generate assemblies spanning a wide range of initial packing fractions within a given polydispersity class, we introduce a numerical control procedure based on the removal of floating particles. This procedure is not intended to mimic a physical densification or loosening mechanism, but rather to systematically vary the packing fraction while preserving the same load-bearing contact network after consolidation. The number of removed floating is given by the parameter $\beta$, defined as the fraction of floating particles $N_f$ remaining in the sample after this process. Clearly, this number depends on $\beta$ and the polydispersity: $N_f(\beta, S)$. Since floating particles do not contribute to force transmission or particle interlocking, removing them at the end of the consolidation stage does not modify the topology of the load-bearing contact network of the consolidated sample, nor the force anisotropy within a given

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Stress ratio $\sigma_{xx}/\sigma_{yy}$ as a function of consolidation time for all samples.</figcaption></figure>

polydispersity set. It is worth-noting that the floating-particle removal procedure does not correspond to a laboratory preparation method and is not intended to reproduce a specific densification or loosening process. It is used here as a numerical tool to generate initial samples spanning a wide range of packing fractions within a given polydispersity, while preserving the load-bearing contact network at the end of the consolidation stage.

For each value of $S$, we generated a set of eleven samples by uniformly varying $\beta$ within [0, 1]. This approach allows us to produce samples with identical contact network micro-structure after consolidation, but different packing fractions. Fig. 5(a) shows that for $S = 0.1$, $\phi$ varies from 0.78 to 0.83 across $\beta = 0$ to 1, respectively; similarly, the range of $\phi$ for $S = 0.8$ over all $\beta$ values is from 0.82 to 0.88. A total of 88 samples were prepared in this manner. Fig. 6 presents examples of initial configurations for two samples with $S = 0.8$: corresponding to (a) a loose state with $\phi = 0.82$ ($\beta = 0.0$), and (b) a medium-dense configuration $\phi = 0.86$ ($\beta = 0.7$). Within a given $S$, samples with the lowest packing fraction are referred to as loose, those with the highest packing fraction as dense, and intermediate cases as medium-dense.

The removal of floating particles modifies the initial PSD, with the extent of this modification increasing with $S$ and decreasing with $\beta$. However, the resulting changes in grading remain relatively minor, as illustrated in Fig. 5(b), which shows that the removal process primarily affects the finer fractions. Since the coarsest particles are consistently engaged in force chains, the number of grains in the coarser size classes remains unchanged. Conversely, even when the maximum number of floating particles $N_f(0, S)$ is removed, the inlet in Fig. 5(b) shows that the proportion of grains of size $d_{min}$ decreases by less than 3% across all $S$.

To assess the robustness of the floating-particle removal procedure with respect to its inherent randomness, additional preparation tests were conducted. For three representative polydispersity levels ($S = 0.1, 0.5$, and $0.8$), the entire preparation protocol — including deposition, consolidation, and floating-particle removal — was repeated using three independent random seeds differing only in the stochastic aspects of particle deposition; thus exhibiting distinct initial configurations. This resulted in a total of 154 preparation runs. For each $S$, the standard deviation of the solid fraction across realizations remained below 0.023, while the standard deviation of the coordination number did not exceed 0.042. These variations are small compared to the corresponding mean values and to the intrinsic fluctuations expected in granular assemblies. Consequently, these additional results demonstrate that the preparation protocol yields statistically robust initial states and that the trends reported in the following sections are not sensitive to the specific realization of the random removal process.

To further assess whether the removal of floating particles affects the initial void topology, we computed the probability distribution
<figure><img src="image_7.png" /><figcaption>Fig. 5. (a) Packing fraction and (b) PSD of $S = 0.8$ after removing floating particles across all $\beta$; inset shows the change in the percent finer than $d_{min}$ for $\beta = 0$ against $S$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Snapshots of the initial state for two samples with $S = 0.8$: (a) a loose sample with $\beta = 0.0$, and (b) a medium-dense sample with $\beta = 0.7$. Colors indicate the number of contacts per particle.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Probability distribution functions of the normalized local density $\phi_L/\langle\phi_L\rangle$ at the end of the compression stage for representative polydispersity cases ($S \in [0.1, 0.4, 0.8]$). Four cases within each $S$ case are presented: $\beta \in [0.0, 0.5, 0.7, 1.0]$.</figcaption></figure>

functions of the normalized local packing fraction $\phi_L/\langle\phi_L\rangle$ at the end of the consolidation stage. To this end, we computed the local solid fraction $\phi_L$ by discretizing the sample into a regular mesh of cell size $\Delta x = \Delta y = d_{max}$. This resolution was chosen after a sensitivity analysis to balance spatial detail and computational cost. At each time step, $\phi_L$ was measured as the ratio of grain area to cell area, allowing us to track how density is redistributed as $S$ and $\beta$ change. Fig. 7 presents the distributions for representative values of the $S \in \{0.1, 0.4, 0.8\}$ and all values of $\beta$. For all polydispersities, the distributions are close to Gaussian and centered around the mean value. Increasing $\beta$ mainly leads to a sharpening of the distribution peak, reflecting a progressive reduction of local density fluctuations. This effect is barely noticeable for quasi-monodisperse samples and becomes increasingly pronounced as the size polydispersity increases, while the overall shape of the distributions remains unchanged. These observations indicate that the floating-particle removal procedure does not significantly alter the initial void topology, but rather provides a controlled way to adjust density heterogeneity prior to shearing.

## 2.4. Simple shear test

After the consolidation stage, each sample was subjected to a strain-controlled simple shear test under constant volume (isochoric) conditions. To enforce volume conservation, the top and bottom walls of the box were kept fixed in the vertical direction. As shown in Fig. 3b, shearing was applied by imposing a constant horizontal velocity $v_x$ on the upper and lower walls, in opposite directions. To minimize slippage at the boundaries and ensure that shear deformation was transmitted uniformly into the granular bulk, a layer of particles was glued to both the top and bottom walls. These glued grains, selected among those deposited during the preparation phase and later consolidated, act as rough elements and are shown in white in Fig. 3b.

The intensity of shear can be characterized by the inertial number $I = \dot{\gamma}d\sqrt{\rho/p}$, where $\dot{\gamma}$ is the shear rate, $d$ is the mean particle diameter, $\rho$ is the particle density, and $p$ is the mean pressure. Note that, since the system evolves under constant volume, $p$ is not fixed but evolves during shearing. Therefore, to ensure a consistent and quasi-static regime at
<figure><img src="image_10.png" /><figcaption>Fig. 8. Stress paths of all the numerical constant volume shear tests; dashed red lines correspond to the critical state line $q/p = 0.3$; red points indicate the average stress values at critical state.</figcaption></figure>

the beginning of shearing, a horizontal velocity $v_x$ was imposed such that the initial inertial number satisfies $I < 10^{-3}$ for a reference pressure of $p = P_0$. This procedure results in shear rates ranging from $\dot{\gamma} = 4.6 \times 10^{-5} s^{-1}$ for the quasi-monodisperse samples $S = 0.1$, to $\dot{\gamma} = 1.1 \times 10^{-4} s^{-1}$ for the most polydisperse ones ($S = 0.8$), ensuring all tests are initiated under quasi-static flow conditions.

### 3. Macroscopic response

The macroscopic mechanical response of the granular assembly is characterized using the stress tensor, defined as:

$$ \sigma_{ij} = \frac{1}{V} \sum_{c} f_{i}^{c} f_{j}^{c}, \qquad (2) $$

where indices $i$ and $j$ run over the Cartesian coordinates $x$ and $y$, the term $f_i^c$ denotes the $i$th component of the contact force at contact $c$, and $f_j^c$ the $j$th component of the branch vector connecting the centers of mass of the two particles in contact.

Since we consider a dry granular medium, all stress components represent effective stresses. We define the mean stress $p$ and the deviatoric stress $q$ as:

$$ p = \frac{\sigma_1 + \sigma_2}{2}, \quad q = \frac{\sigma_1 - \sigma_2}{2}, \qquad (3) $$

where $\sigma_1$ and $\sigma_2$ are the principal components of the stress tensor.

Fig. 8 shows the stress states reached by each sample in the $(p, q)$ plane for all simulations. Each set corresponds to a fixed polydispersity index $S$, with varying values of $\beta$, from loose ($\beta = 0$) to dense initial states ($\beta = 1$). The data at large shear strains generally collapse onto a linear trend corresponding to the critical state line. A linear fitting of all the tests at critical state results in $q/p = 0.3$, with a coefficient of determination $R^2 = 0.99$. As discussed in the introduction, this is consistent with previous results suggesting that the internal friction angle is not affected by the PSD.

As expected, in denser samples ($\beta \to 1$) the tendency to dilate develops relatively high stresses along the critical state line. Moreover, higher values of both $q$ and $p$ are attained in samples with higher $S$, suggesting that polydispersity promotes stronger packing configurations capable of sustaining greater loads. In contrast, stresses tend to vanish in loose samples ($\beta \to 0$), which is indicative of mechanical instability

and is consistent with the phenomenon of static liquefaction. These initial observations qualitatively suggest that a substantial portion of the loose specimens undergoes a liquefaction-like transition during shearing.

To better distinguish the samples that liquefy from those remaining stable, Fig. 9 presents the stress-strain responses of all numerical tests: (a) $q/P_0$ and (b) $p/P_0$ are plotted versus $\gamma$ up to 0.3. Note that Fig. 9 is not intended to provide a quantitative comparison of shear strength across different values of $S$, since polydispersity and initial packing fraction are intrinsically correlated. Instead, this figure is used to identify and classify the different mechanical responses as a function of $S$ and $\beta$. Since the peak stresses occur at very low strains of $\gamma \le 0.01$, both plots employ logarithmic scaling for the strain axis to enhance the visualization. Three distinct mechanical responses are typically observed. First, liquefaction is characterized by a rapid and irreversible loss of both $q$ and $p$ after reaching a peak. This behavior, observed in loose samples with low $\beta$, reflects a transition to mechanical instability and is known as static liquefaction. Second, in dense samples with large $\beta$, both $q$ and $p$ increase with strain and reach a plateau at large $\gamma$, corresponding to a steady critical state. These samples exhibit strain hardening — i.e. increase in shear strength upon undrained shear deformation — and maintain mechanical integrity throughout loading. Third, an intermediate behavior that we call “temporary liquefaction”, occurs in samples with intermediate values of $\beta$. These samples exhibit an initial increase in $q$ and $p$, followed by a sharp drop and subsequent recovery toward a critical state at large strains. This transient collapse followed by re-strengthening indicates a temporary loss of resistance that is recovered upon continued shearing. These behaviors are consistently found across different values of polydispersity $S$, but the $\beta$ ranges associated with each regime depend strongly on $S$. For instance, in the case of $S = 0.8$, samples with $\beta \le 0.6$ undergo full liquefaction, while $\beta = 0.7$ shows temporary liquefaction, and $\beta \ge 0.8$ does not exhibit liquefaction and reaches a critical strength plateau. Conversely, for $S = 0.1$, all samples undergo static liquefaction except $\beta = 1$, which exhibits temporary liquefaction. Furthermore, temporary liquefaction appears across different combinations of $\beta$ and $S$—specifically for ($S = 0.1, \beta = 1.0$), ($S = 0.7, \beta = 0.6$), and ($S = 0.8, \beta = 0.8$).

In the following section, we shift our focus toward a micro-structural description of the assemblies, aiming to identify the underlying mechanisms responsible for these macroscopic behaviors. In particular, we
<figure><img src="image_11.png" /><figcaption>Fig. 9. Stress–strain response for all tests: (a) deviatoric stress q and (b) mean stress p as a function of the shear strain γ, for samples from monodisperse (S = 0.1) to highly polydisperse (S = 0.8), and β ∈ [0, 1].</figcaption></figure>

investigate how particle connectivity and packing homogeneity evolve during shear, and how these structural features correlate with the observed mechanical responses.

## 4. Microscopic analysis

As defined in Eq. (2), the granular stress tensor is constructed from quantities associated with the contact and force network. This expression can be conveniently rewritten as an average over all contacts as:

$$ \sigma_{ij} = n_c \langle f_i^c l_j^c \rangle_{\forall c \in V}, \quad (4) $$

where $n_c = N_c/V$ is the contact density, with $N_c$ the total number of active contacts (i.e., contact carrying a non-zero force) within the assembly. Then, it is possible to demonstrate that $n_c$ is itself related to two key micro-structural descriptors: $\phi$ and the mechanical coordination number $z_M$, defined as the average number of active contacts per particle, excluding grains with 1 and 0 contacts. Indeed, for a two-dimensional packing of discs with mean grain diameter $\langle d \rangle$, one obtains $n_c = 2z_M\phi/(\pi\langle d \rangle^2)$ by assuming $\langle d^2 \rangle \approx \langle d \rangle^2$. While this

approximation slightly underestimates $\langle d^2 \rangle$ in the most polydisperse cases (by up to about 20% for S = 0.8), it remains accurate enough for the present analysis and is consistent with the trends observed in our simulations across all S and β. Consequently, the stress tensor can be expressed in terms of simple micro-structural parameters — solid fraction and coordination number — along with the average contact force and branch length, using the following form (Agnolin and Roux, 2007; Khalili et al., 2017):

$$ \sigma_{ij} \approx 2 \frac{z_M \phi}{\pi \langle d \rangle^2} \langle f_i^c l_j^c \rangle_{\forall c \in V}, \quad (5) $$

This highlights that variations in connectivity and local density directly control the macroscopic shear response. In the following, we analyze these two descriptors separately: first the evolution of the contact network, then the heterogeneity of local packing density.

### 4.1. Connectivity of the contact network

For two-dimensional packings of rigid frictional grains, mechanical stability is typically achieved for $z_M$ lying between 3 and 4 (Andreotti
<figure><img src="image_12.png" /><figcaption>Fig. 10. Mechanical coordination number $z_M$ as function of shear strain $\gamma$ for three sets $S = 0.1$, $S = 0.7$ and $S = 0.8$.</figcaption></figure>

et al., 2013; van Hecke, 2009); values below this range indicate an unstable state. Fig. 10 shows the evolution of $z_M$ with $\gamma$ for three representative polydispersity levels. Initially, all samples start with essentially the same initial connectivity ($z_M \approx 3.5$ at $\gamma \le 10^{-3}$), demonstrating that the initial micro-structure is comparable across sets and that removing floating grains during sample preparation does not alter it. In non-liquefying samples, $z_M$ gradually increases as the contact network strengthens. In fully liquefied samples, by contrast, $z_M$ experiences an abrupt drop to values below 2.5 and remains low throughout the test, signaling a breakdown of the contact network that coincides with the macroscopic loss of shear strength. Temporary liquefaction exhibits an intermediate behavior: after a sharp decrease to around 3.3–3.4, $z_M$ partially recovers toward its initial value. However, in some cases this recovery remains insufficient to exceed the apparent threshold of $z_m \sim 3.4$, leaving the material in a mechanically unstable state. The magnitude and timing of this drop depend strongly on both density and $S$: higher polydispersity softens the abruptness of the drop, while increasing density delays its onset. These observations suggest that, in our samples, a critical coordination number of about 3.4 serves as microscopic indicator of liquefaction. Similar findings were reported by Yang et al. (2022) and Zhang et al. (2023). It should be emphasized that the critical coordination number $z_M \approx 3.4$ identified here is not intended as a universal threshold, as its absolute value may depend on dimensionality and particle properties. However, within the present framework of 2D assemblies of rigid frictional grains prepared with identical force-bearing micro-structures at the end of the consolidation stage, this threshold emerges consistently across a wide range of initial densities and size polydispersities (Kruyt and Rothenburg, 2014; Azéma et al., 2007; Taghizadeh et al., 2017).

While $z_M$ provides a global average measure of connectivity, a complementary perspective is obtained by examining the probability distribution of contacts $P_c$, which quantifies the fraction of particles with $c$ contacts. This descriptor highlights how different contact numbers contribute to the overall network structure (Nguyen et al., 2014). Fig. 11 shows the evolution of $P_c$ for $S = 0.8$ (polydisperse) and $S = 0.1$ (quasi-monodisperse), in different mechanical behaviors (liquefaction, temporary liquefaction, and hardening). In the non-liquefying case at $S = 0.8$, the distributions remain virtually unchanged throughout shear, with $P_3$ and $P_4$ dominating the connectivity, while $P_0 \sim 0.35$. Thus, roughly 90% of the particles have 0 to 4 contacts, and fewer than 10% of the grains form more than 4 contacts. By contrast, liquefying samples show a drastic reorganization of the contact network. For the highly polydisperse case ($S, \beta = (0.8, 0.0)$, $P_3$ collapses from 50% to about 10% at $\gamma \sim 0.0025$, while $P_0$ surges to nearly 100%, indicating a fully disconnected state. A similar trend occurs in the quasi-monodisperse liquefied sample ($S, \beta = (0.1, 0)$, where $P_3$ and $P_4$ initially account for about 80% of the connectivity before dropping sharply; note that these samples initially have no floaters (i.e., $P_0 = P_1 = 0$). Temporary liquefaction cases ($S, \beta = (0.8, 0.7)$ and ($S, \beta = (0.1, 1.0)$ follow a similar initial trajectory as liquefying cases, with a sudden rise of $P_0$ and collapse of $P_3$. However, after recovery, $P_0$ decreases and $P_3$ rises back toward

their pre-liquefaction values, restoring partial connectivity. This re-covery mechanism distinguishes temporary liquefaction from complete liquefaction and reflects the reformation of a load-bearing network.

These analyses demonstrate that liquefaction corresponds to a sharp and nearly complete breakdown of contact connectivity, whereas temporary liquefaction represents a transient loss followed by partial re-covery. However, connectivity alone does not fully explain why some samples recover and others do not. To address this point, we now examine local density variations during shear, which provide complementary insights into the mechanisms behind recovery or persistent failure.

## 4.2. Local density

The connectivity analysis revealed that the collapse of the contact network is a key marker of liquefaction. Moreover, the stress tensor formulation (Eq. (5)) shows that both connectivity and packing fraction contribute to the macroscopic response. Potential differences between liquefaction and temporary liquefaction should arise from local variations in density. Capturing these heterogeneities is therefore essential to explain why some samples recover their strength while others do not. The local packing fraction is computed following the method described in Section 2. For the most polydisperse set $S = 0.8$ and three values of $\beta$, Fig. 12(a) presents the probability distribution functions (PDF) of the normalized local density $\phi_L/\langle\phi_L\rangle$ over $\gamma$, where $\langle\phi_L\rangle$ is the mean local density of the sample. At the beginning of the test, the liquefying sample already exhibits slightly looser local regions than the temporary and non-liquefying cases. As shear progresses, its PDF shifts markedly toward lower local densities and develops long tails, revealing the emergence of extended low-density zones. In contrast, the PDFs for the temporary and non-liquefying cases remain narrow and largely stable, indicating limited local rearrangements despite ongoing shear.

To quantify the temporal variations of the local density, Fig. 12(b) displays the standard deviation (STD) of the normalized local density during shear. This metric provides a simple measure of how heterogeneous the sample becomes upon constant volume shearing. Across different levels of polydispersity $S \in \{0.1, 0.7, 0.8\}$, non-liquefying and temporary-liquefying samples maintained low, practically constant STD values. In contrast, fully liquefied samples exhibit a sharp increase in STD that fluctuates with shear strain, coinciding with the collapse of connectivity reported earlier. This correlation indicates that liquefaction involves not only contact loss but also the development of strong local density fluctuations, whereas recovery is possible only when these heterogeneities remain moderate.

Although statistical indicators (PDFs and STD) reveal the presence of density fluctuations, they do not provide spatial information regarding the location of these heterogeneities. To clarify it, we examine the spatial distribution of $\phi_L$ during shear to identify whether liquefaction is associated with localized voids or a more diffuse rearrangement. Fig. 13 displays $\phi_L$ maps for (a) a case of liquefaction and (b) a temporary liquefaction case, at key deformation stages: initial state,
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Probability of having *c* contact *P*<sub>c</sub> as function of *γ* for two sets: (a) polydisperse *S* = 0.8 and (b) quasi-monodisperse *S* = 0.1. The polydisperse case (a) shows 3 cases: liquefaction *β* = 0.0, temporary liquefaction *β* = 0.7, and Hardening *β* = 1.0. The quasi-monodisperse case shows 2 cases: liquefaction *β* = 0.0 and temporary liquefaction *β* = 1.0.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** (a) Probability distribution functions (PDF) of the local density $\phi_L$ normalized by the average local density $\langle \phi_L \rangle$ throughout shear-strain for three cases of most polydisperse case $S = 0.8$: liquefaction ($\beta = 0.0$), temporary liquefaction ($\beta = 0.7$) and hardening ($\beta = 1.0$). (b) Standard deviation (STD) of the distribution of $\phi_L/\langle \phi_L \rangle$ as a function of $\gamma$, for samples that exhibit liquefaction, temporary liquefaction and Hardening.</figcaption></figure>

peak stress, minimum stress (post-liquefaction), and critical state. In the liquefied sample, large macro-voids appear after the peak, near the sample boundaries, and progressively migrate inward, persisting until the end of the test. Conversely, the temporary liquefaction case exhibits transient low-density zones that dissipate during continued shear, leading to the recovery of connectivity and macroscopic strength. Similar conclusions were reached experimentally by Sukhumkitcharoen et al. (2024), where liquefaction resistance was shown to be strongly

controlled by the spatial organization of local voids. The present numerical results provide complementary evidence that the development of extended low-density regions is a key micro-structural mechanism underlying static liquefaction.

The heterogeneity of the local density is also reflected in the normal contact force network. Fig. 14 illustrates the evolution of force chains at different shear strains for two representative cases: (a) fully liquefying sample (S, β) = (0.1, 0.0) and (b) temporarily liquefying sample
<figure><img src="image_15.png" /><figcaption>Fig. 13. Color map of the local density $\phi_L$ for different shear strain values for (a) liquefaction case in a quasi-monodisperse sample ($S, \beta$) = (0.1, 0), and (b) temporary liquefaction in a polydisperse sample ($S, \beta$) = (0.8, 0.7).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Normal contact force networks for different shear strain levels: (a) liquefaction case ($S, \beta$) = (0.1, 0.0) and (b) temporary liquefaction case ($S, \beta$) = (0.8, 0.7); red lines represent the strong force network ($\phi > \phi_n$) and black lines are the weak network ($\phi \le \phi_n$).</figcaption></figure>

($S, \beta$) = (0.8, 0.7). The network is represented by strong forces (greater than the mean normal force $\langle f_n \rangle$) and weak forces ($\le \langle f_n \rangle$). In the initial state, force chains percolate throughout the assembly in both cases, indicating a mechanically stable load-bearing structure. Upon shearing, however, markedly different behaviors emerge. In the fully liquefying case, force-chain percolation progressively degrades, and force-free regions first develop near the top and bottom boundaries. By contrast, in the temporarily liquefying case, the network experiences only a transient loss of connectivity, followed by reorganization and re-percolation, without the formation of low-force or void regions adjacent to the boundaries. Overall, these observations indicate that

boundary-adjacent regions are the first to respond to the collapse of the force-chain network. Nevertheless, the subsequent development of low-density zones is governed by the ability of the contact network to reorganize, rather than by an artificial boundary effect.

### 4.3. Connecting micro-structural stresses to macroscopic response

At the macroscopic scale, the onset of full or temporary liquefaction was systematically associated with a sharp drop in mean stress $p$ (see Fig. 9b). To understand this behavior, we can combine the definition
<figure><img src="image_17.png" /><figcaption>Fig. 15. Comparison of the mean stress obtained from the granular stress tensor (p) and from the micro-parameter-based expression (p<sub>micro</sub>) for two cases: S = 0.1 (left) and S = 0.8 (right), each one for two samples β = 0 and 1.</figcaption></figure>

of the mean stress p (Eq. (3)) together with the micro-mechanical expression of the granular stress tensor in Eq. (5). Assuming that the branch vectors and contact forces are statistically uncorrelated, a condition well verified in our simulations, p can be written in terms of micro-structural descriptors:

$$p_{\text{micro}} \sim \frac{z_M \phi \langle f \rangle}{\pi \langle d \rangle}, \qquad (6)$$

where ⟨−⟩ denotes a mean value, thus ⟨f⟩/⟨d⟩ is a measure of the mean local stress. This micro-to-macro relation is classical in granular mechanics and has been widely used in various contexts (Agnolin and Roux, 2007), such as identifying mechanism of compaction in soft granular materials (Cantor et al., 2020; Cárdenas-Barrantes et al., 2021) or rupture in wet granular media (Khamseh et al., 2015). Note that Eq. (6) does not rely on a small-strain assumption and remains formally valid at large deformation, as it directly derives from the microscopic definition of the stress tensor. However, its simplified form assumes weak correlations between contact forces and branch vectors, as well as negligible cross-effects between force anisotropy and fabric anisotropy. In highly polydisperse samples, these assumptions may be less strictly applicable than in monodisperse samples owing to enhanced geometric heterogeneity, which explains the slightly larger discrepancies observed in the most polydisperse cases. In the present study, we test whether Eq. (6) can capture the precursors of liquefaction. Specifically, we compare p obtained directly from the stress tensor (Eq. (3)), with p<sub>micro</sub> from Eq. (2). Fig. 15 shows this comparison for the quasi-monodisperse S = 0.1 and most polydisperse S = 0.8 samples. Despite small quantitative differences, p<sub>micro</sub> closely tracks p over the entire shearing path, including the abrupt drop associated with liquefaction. The small discrepancy arises from the assumption that the branch vectors and contact forces in Eq. (4) are uncorrelated across all polydispersities. Nevertheless, Polanía et al. (2025) reported that this correlation increases slightly with higher polydispersity. This close match demonstrates that the evolution of p, and thus the macroscopic strength, is primarily governed by micro-structural variables, namely coordination number and packing fraction.

Linking micro- to macro-scale mechanisms reveals that liquefaction emerges from a dual process: the collapse of the contact network combined with the formation of extended low-density regions. In temporary liquefaction, this is manifested as a transient loss of connectivity with limited density heterogeneity, allowing recovery. In contrast, full liquefaction involves a severely disconnected and highly heterogeneous micro-structure, preventing any structural recovery.

## 5. Conclusions

Static liquefaction remains one of the most critical failure modes in tailings storage facilities and granular waste deposits. Although its macroscopic features are well established, the particle-scale mechanisms that trigger instability remain elusive. In this work, we addressed

this gap using two-dimensional numerical simulations based on the Discrete Element Method in its Contact Dynamics formulation. We extended the method to include contact elasticity, allowing us to capture small force fluctuations and network rearrangements during instability. We systematically varied the initial density of the assemblies, while preserving a constant force-bearing micro-structure. To accomplish this, we prepared samples by selectively removing a controlled fraction of floating particles after consolidation. This approach allowed us to isolate the role of density from other micro-structural factors and explore the combined effect of density and size polydispersity on liquefaction.

Across eight particle size distributions, ranging from quasi-monodisperse to polydisperse samples, each spanning densities from loose to dense states, and sheared to critical states at large strain of 0.3, three distinct mechanical responses emerged: full liquefaction, temporary liquefaction and hardening. As expected, loose samples systematically liquefied, while dense samples did not. Intermediate cases displayed temporary liquefaction, characterized by a transient drop in shear resistance, followed by recovery at large strains. Increasing size polydispersity and initial density enhanced resistance to liquefaction and, at large strains, led to higher ultimate shear strengths. Our results are consistent with trends reported in both experimental and theoretical studies.

An important contribution of this study lies in the multi-scale analysis linking micro-structural evolution to macroscopic response. The onset of liquefaction coincided with the collapse of the contact network, captured by the mechanical coordination number dropping below ~3, regardless of size polydispersity or initial density. Simultaneously, local packing fraction maps revealed the formation of extended low-density regions in fully liquefied cases, while temporarily liquefied samples maintained comparatively homogeneous density fields. By combining these insights, we identify a dual mechanism for liquefaction: loss of connectivity at the contact scale, coupled with the emergence of heterogeneous density domains. Temporary liquefaction corresponds to transient connectivity loss with limited density heterogeneity, enabling recovery. Full liquefaction involves both effects, producing a disconnected and highly heterogeneous micro-structure from which recovery is impossible.

The perspective of this study is to apply multi-scale approaches to investigate the role of the material fabric on liquefaction behavior. For instance, in anisotropic structured soils the undrained effective stress path strongly depends on the principal stress directions, and might significantly increase the static liquefaction susceptibility along a given orientation. The evolution of the micro-structure is crucial in that case. Likewise, the influence of non-plastic fines in undrained behavior could be tracked down to the interactions between fines and coarse grains forming unstable force-chains. Also of major relevance, this study should be extended to 3D systems, larger assemblies, and to non-spherical or angular particles, where dimensionality, finite-size
effects, and shape-induced interlocking may modify the absolute value of liquefaction thresholds. Nevertheless, the onset of liquefaction is anticipated to remain associated with a collapse of the contact network, even if the critical coordination number shifts in 3D. Also exploring coupled hydro-mechanical effects and validating these findings against physical model experiments would provide further insights into the relevance of these micro-scale mechanisms for real liquefiable tailings and natural granular soils.

From an engineering perspective, the mechanisms identified in this study provide a particle-scale interpretation of liquefaction triggering that complements classical state-based approaches used in practice. The collapse of the contact network offers a micro-mechanical explanation for the sudden loss of effective stress observed during undrained loading, while the development of extended low-density regions sheds light on the irreversible nature of full liquefaction in loose granular deposits. The existence of temporary liquefaction highlights that transient instability does not necessarily imply permanent loss of strength, which may help explain post-instability strength recovery observed in some tailings storage facilities and granular embankments (Wang et al., 2002; Verdugo, 2024). Finally, the stabilizing role of size polydispersity provides micro-structural support to common engineering strategies favoring well-graded materials and densification to mitigate liquefaction susceptibility. Together, these findings bridge particle-scale mechanisms and engineering observations, and contribute to a more mechanistic understanding of liquefaction risk in granular soils and tailings.

### CRediT authorship contribution statement

**Carolina Castro-Malaver:** Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis. **Manuel Cárdenas-Barrantes:** Writing – review & editing, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. **David Cantor:** Writing – review & editing, Validation, Software, Methodology, Formal analysis. **Mathieu Renouf:** Writing – review & editing, Validation, Supervision, Formal analysis, Conceptualization. **Emilien Azéma:** Writing – review & editing, Supervision, Funding acquisition, Formal analysis, Conceptualization. **Carlos Ovalle:** Writing – review & editing, Supervision, Resources, Funding acquisition, Formal analysis, Conceptualization. **Arcesio Lizcano:** Writing – review & editing, Funding acquisition, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors gratefully acknowledge the financial support provided by MITACS [project Accelerate IT30341], SRK Consulting Vancouver, and the industrial partners of the Research Institute on Mines and Environment UQAT-Polytechnique (irme.ca/en).

### Data availability

Data will be made available on request.

### References

Agnolin, I., Roux, J.-N., 2007. Internal states of model isotropic granular packings. I. Assembling process, geometry, and contact networks. Phys. Rev. E 76, 061302. http://dx.doi.org/10.1103/PhysRevE.76.061302.

Ahmed, S.S., Martinez, A., DeJong, J.T., 2023. Effect of gradation on the strength and stress-dilation behavior of coarse-grained soils in drained and undrained triaxial compression. J. Geotech. Geoenvironmental Eng. 149 (5), 04023019. http://dx.doi.org/10.1061/JGGEFK.GTENG-10972.

Amirpour Harehdasht, S., Hussien, M.N., Karray, M., Roubtsova, V., Chekired, M., 2018. Influence of particle size and gradation on shear strength-dilation relation of granular materials. Can. Geotech. J. 56 (2), 208-227. http://dx.doi.org/10.1139/cgj-2017-0468.

Andreotti, B., Forterre, Y., Pouliquen, O., 2013. Granular Media: Between Fluid and Solid. Cambridge University Press.

Azéma, E., Radjai, F., Peyroux, R., Saussine, G., 2007. Force transmission in a packing of pentagonal particles. Phys. Rev. E 76, 011301. http://dx.doi.org/10.1103/PhysRevE.76.011301, URL https://link.aps.org/doi/10.1103/PhysRevE.76.011301.

Banerjee, S.K., Yang, M., Taiebat, M., 2023. Effect of coefficient of uniformity on cyclic liquefaction resistance of granular materials. Comput. Geotech. 155, 105232. http://dx.doi.org/10.1016/j.compgeo.2022.105232.

Basson, M.S., Martinez, A., DeJong, J.T., 2024a. DEM investigation of the effect of gradation on the strength, dilatancy, and fabric evolution of coarse-grained soils. J. Geotech. Geoenvironmental Eng. 150 (8), 04024060. http://dx.doi.org/10.1061/JGGEFK.GTENG-12310.

Basson, M.S., Martinez, A., DeJong, J.T., 2024b. DEM simulations of the liquefaction resistance and post-liquefaction strain accumulation of coarse-grained soils with varying gradations. Comput. Geotech. 174, 106649. http://dx.doi.org/10.1016/j.compgeo.2024.106649.

Been, K., Jefferies, M.G., 1985. A state parameter for sands. "Géotechnique" 35 (2), 99-112. http://dx.doi.org/10.1680/geot.1985.35.2.99.

Blight, G.E., Fourie, A., 2005. Catastrophe revisited – disastrous flow failures of mine and municipal solid waste. Geotech. Geol. Eng. 23, 219-248. http://dx.doi.org/10.1007/s10706-004-7067-y.

Cantor, D., Cárdenas-Barrantes, M., Preechawuttipong, I., Renouf, M., Azéma, E., 2020. Compaction model for highly deformable particle assemblies. Phys. Rev. Lett. 124, 208003. http://dx.doi.org/10.1103/PhysRevLett.124.208003.

Cantor, D., Ovalle, C., 2025. Sample size effects on the critical state shear strength of granular materials with varied gradation and the role of column-like local structures. Géotechnique 75 (1), 29-40. http://dx.doi.org/10.1680/geot.23.00032.

Cárdenas-Barrantes, M., Cantor, D., Bares, J., Renouf, M., Azéma, E., 2021. Micromechanical description of the compaction of soft pentagon assemblies. Phys. Rev. E. 103, 062902. http://dx.doi.org/10.1103/PhysRevE.103.062902.

Cárdenas-Barrantes, M., Ovalle, C., 2025. Multiscale insights into sliding surface liquefaction through DEM simulations. Comput. Geotech. 183, 107191. http://dx.doi.org/10.1016/j.compgeo.2025.107191.

Castro, G., 1969. Liquefaction of Sands (Ph.D. thesis). Harvard University, Cambridge, MA, p. 112, reprinted as Harvard Soil Mechanics series No. 81.

Coppin, N., Henry, M., Cabrera, M., Azéma, E., Dubois, F., Legat, V., Lambrechts, J., 2023. Collapse dynamics of two-dimensional dry and immersed granular columns of elongated grains. Phys. Rev. Fluids 8, 094303. http://dx.doi.org/10.1103/PhysRevFluids.8.094303.

da Cruz, F., Emam, S., Prochnow, M., Roux, J.-N., Chevoir, F., 2005. Rheophysics of dense granular materials: Discrete simulation of plane shear flows. Phys. Rev. E 72, 021309. http://dx.doi.org/10.1103/PhysRevE.72.021309.

Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Géotechnique 29 (1), 47-65. http://dx.doi.org/10.1680/geot.1979.29.1.47.

Dawson, R., Morgenstern, N., Stokes, A., 1998. Liquefaction flowslides in rocky mountain coal mine waste dumps. Can. Geotech. J. 35 (2), 328-343. http://dx.doi.org/10.1139/t98-009.

Dobry, R., Álvarez, L., 1967. Seismic failures of Chilean tailings dams. J. Soil Mech. 237-269. http://dx.doi.org/10.1061/JSFEAQ.0001054.

Dubois, F., Acary, V., Jean, M., 2018. The contact dynamics method: A nonsmooth story. Comptes Rendus Mécanique 346 (3), 247-262. http://dx.doi.org/10.1016/j.crme.2017.12.009.

Dubois, F., Jean, M., Renouf, M., Mozul, R., Martin, A., Bagnéris, M., 2011. LMGC90. In: 10e colloque national en calcul des structures. Giens, France.

Dyvik, R., Berre, T., Lacasse, S., Raadim, B., 1987. Comparison of truly undrained and constant volume direct simple shear tests. Géotechnique 37 (1), 3-10. http://dx.doi.org/10.1680/geot.1987.37.1.3.

Evans, T.M., Zhang, L., 2019. A numerical study of particle friction and initial state effects on the liquefaction of granular assemblies. Soil Dyn. Earthq. Eng. 126, 105773. http://dx.doi.org/10.1016/j.soildyn.2019.105773.

Fourie, A.B., Blight, G.E., Papageorgiou, G., 2001. Static liquefaction as a possible explanation for the merriesfruit tailings dam failure. Can. Geotech. J. 38 (4), 707-719. http://dx.doi.org/10.1139/t00-112.

Girumugisha, G., Ovalle, C., Ouellet, S., 2024. Grading scalping and sample size effects on critical shear strength of mine waste rock through laboratory and in-situ testing. Int. J. Rock Mech. Min. Sci. 183, 105915. http://dx.doi.org/10.1016/j.ijrmms.2024.105915.
Gobbi, S., Reiffsteck, P., Lenti, L., d'Avila, M.P.S., Semblat, J.-F., 2022. Liquefaction triggering in silty sands: effects of non-plastic fines and mixture-packing conditions. *Acta Geotech.* 17 (2), 391-410.

Gong, G., Thornton, C., Chan, A.H., 2012. DEM simulations of undrained triaxial behavior of granular material. J. Eng. Mech. 138 (6), 560-566. http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000366.

van Hecke, M., 2009. Jamming of soft particles: geometry, mechanics, scaling and isostaticity. J. Phys.: Condens. Matter. 22 (3), 033101. http://dx.doi.org/10.1088/0953-8984/22/3/033101.

Huang, X., Kwok, C., Hanley, K., Wu, W., 2018. DEM analysis of the onset of flow deformation of sands: linking monotonic and cyclic undrained behaviours. *Acta Geotech.* 13 (5), 1061-1074. http://dx.doi.org/10.1007/s11440-018-0664-3.

Huang, C., Xu, S., Shen, Z., Zhang, L., Wang, L., Liu, C., 2023. Effects of non-plastic fines on liquefaction properties of saturated silt using discrete element modeling. *Eng. Geol.* 317, 107091. http://dx.doi.org/10.1016/j.enggeo.2023.107091.

Irani, N., Salimi, M., Golestaneh, P., Tafili, M., Wichtmann, T., Lederer, J., 2024. Deep learning-based analysis of true triaxial DEM simulations: Role of fabric and particle aspect ratio. Comput. Geotech. 173, 106529. http://dx.doi.org/10.1016/j.compgeo.2024.106529.

Karim, M.E., Alam, M.J., 2014. Effect of non-plastic silt content on the liquefaction behavior of sand-silt mixture. Soil Dyn. Earthq. Eng. 65, 142-150. http://dx.doi.org/10.1016/j.soildyn.2014.06.010.

Khalili, M.H., Roux, J.N., Pereira, J.M., Brisard, S., Bornert, M., 2017. Numerical study of one-dimensional compression of granular materials. I. Stress-strain behavior, microstructure, and irreversibility. Phys. Rev. E 95, 032907. http://dx.doi.org/10.1103/PhysRevE.95.032907.

Khamseh, S., Roux, J.N., Chevoir, F., 2015. Flow of wet granular materials: A numerical study. Phys. Rev. E 92, 022201. http://dx.doi.org/10.1103/PhysRevE.92.022201.

Krabbenhoft, K., Huang, J., da Silva, M.V., Löwe, T., 2012. Granular contact dynamics with particle elasticity. Granul. Matter 14, 607-619. http://dx.doi.org/10.1007/s10035-012-0360-1.

Kruyt, N., Rothenburg, L., 2014. On micromechanical characteristics of the critical state of two-dimensional granular materials. Acta Mech. 225, 2301-2318. http://dx.doi.org/10.1007/s00707-014-1128-y.

Kwa, K.A., Airey, D.W., 2017. Effects of fines on liquefaction behaviour in well-graded materials. Can. Geotech. J. 54 (10), 1460-1471. http://dx.doi.org/10.1139/cgj-2017-0016.

Lade, P.V., Yamamuro, J.A., 2011. Evaluation of static liquefaction potential of silty sand slopes. Can. Geotech. J. 48 (2), 247-264. http://dx.doi.org/10.1139/T10-063.

Lade, P.V., Yamamuro, J.A., Liggio Jr., C.D., 2009. Effects of fines content on void ratio, compressibility, and static liquefaction of silty sand. Geomech. Eng. 1 (1), 1-15. http://dx.doi.org/10.12989/gae.2009.1.1.001.

Li, X.S., Dafalias, Y.F., 2012. Anisotropic critical state theory: Role of fabric. J. Eng. Mech. 138 (3), 263-275. http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000324.

Liu, G., Rong, G., Peng, J., Zhou, C., 2015. Numerical simulation on undrained triaxial behavior of saturated soil by a fluid coupled-DEM model. Eng. Geol. 193, 256-266. http://dx.doi.org/10.1016/j.enggeo.2015.04.019.

Lü, X., Huang, M., Andrade, J.E., 2017. Predicting the initiation of static liquefaction of cross-anisotropic sands under multiaxial stress conditions. Int. J. Numer. Anal. Methods Geomech. 41 (17), 1724-1740. http://dx.doi.org/10.1002/nag.2697.

Martin, E.L., Thornton, C., Utili, S., 2020. Micromechanical investigation of liquefaction of granular media by cyclic 3D DEM tests. *Géotechnique* 70 (10), 906-915. http://dx.doi.org/10.1680/igeot.18.P.267.

Morgenstern, N.R., Vick, S.G., Van Zyl, D., 2015. Independent Expert Engineering Investigation and Review Panel Report on Mount Polley Tailings Storage Facility Breach. Review Report, British Columbia, Canada.

Morgenstern, N.R., Vick, S.G., Viotti, C., Watts, B., 2016. Report on the Immediate Causes of the Failure of the Fundão Dam, Fundão Tailings Dam Review Panel. Technical Report.

Nguyen, D.H., Azéma, E., Radjaï, F., Sornay, P., 2014. Effect of size polydispersity versus particle shape in dense granular media. Phys. Rev. E 90, 012202. http://dx.doi.org/10.1103/PhysRevE.90.012202.

Subuo, M., Chitrawel, S., Kuwano, R., Hanley, K.J., Kyokawa, H., Koseki, J., 2022. Linking inherent anisotropy with liquefaction phenomena of granular materials by means of DEM analysis. Soils Found. 62 (5), 101202. http://dx.doi.org/10.1016/j.sandf.2022.101202.

Pires-Sturm, A.P., DeJong, J.T., 2022. Influence of particle size and gradation on liquefaction potential and dynamic response. J. Geotech. Geoenvironmental Eng. 148 (6), 04022045. http://dx.doi.org/10.1061/(ASCE)GT.1943-5606.0002799.

Polanía, O., Cabrera, M., Renouf, M., Azéma, E., 2022. Collapse of dry and immersed polydisperse granular columns: A unified runout description. Phys. Rev. Fluids 7, 084304. http://dx.doi.org/10.1103/PhysRevFluids.7.084304.

Polanía, O., Cabrera, M., Renouf, M., Azéma, E., Estrada, N., 2023. Grain size distribution does not affect the residual shear strength of granular materials: An experimental proof. Phys. Rev. E 107, L052901. http://dx.doi.org/10.1103/PhysRevE.107.L052901.

Polanía, O., Renouf, M., Cabrera, M., Estrada, N., Azéma, E., 2025. Monodisperse behavior of polydisperse flows. Phys. Rev. E 11, L043401. http://dx.doi.org/10.1103/PhysRevE.111.L043401.

Quiroz Rojo, P., Cantor, D., Renouf, M., Ovalle, C., Azéma, E., 2025. Rev assessment of granular materials with varied grading based on macro- and micro-mechanical statistical data. Acta Geotech. 20, 1585-1598. http://dx.doi.org/10.1007/s11440-024-02498-3.

Radjaï, F., Dubois, F., 2011. Discrete-Element Modeling of Granular Materials. Wiley-ISTE, p. 425, HAL: hal-00691805, version v1.

Rahman, M.M., Nguyen, H.B.K., Fourie, A.B., Kuhn, M.R., 2021. Critical state soil mechanics for cyclic liquefaction and postliquefaction behavior: DEM study. J. Geotech. Geoenvironmental Eng. 147 (2), 04020166. http://dx.doi.org/10.1061/(ASCE)GT.1943-5606.0002453.

Rahmani, H., Naeini, S.A., 2020. Influence of non-plastic fine on static liquefaction and undrained monotonic behavior of sandy gravel. Eng. Geol. 275, 105729.

Renouf, M., Dubois, F., Alart, P., 2004. A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media. J. Comput. Appl. Math. 168 (1), 375-382. http://dx.doi.org/10.1016/j.cam.2003.05.019.

Robertson, P.K., de Melo, L., Williams, D., Wilson, G., 2019. Report of the Expert Panel on the Technical Causes of the Failure of Feijão Dam I. Technical Report.

Sadrekarimi, A., 2014. Effect of the mode of shear on static liquefaction analysis. J. Geotech. Geoenvironmental Eng. 140 (12), 04041069. http://dx.doi.org/10.1061/(ASCE)GT.1943-5606.0001182.

Salimi, M., Irani, N., Tafili, M., Wichtmann, T., 2025. Micromechanical Influence of Fabric Anisotropy and Stress Path Dependency on Liquefaction Susceptibility in Granular Soils. Open Geomech. 6, 1-18. http://dx.doi.org/10.5802/ogeo.23.

Santamarina, J.C., Torres-Cruz, L.A., Bachus, R.C., 2019. Why coal ash and tailings dam disasters occur. Science 364 (6440), 526-528. http://dx.doi.org/10.1126/science.aax1927.

Shafipour, R., Soroush, A., 2008. Fluid coupled-DEM modelling of undrained behavior of granular media. Comput. Geotech. 35 (5), 673-685. http://dx.doi.org/10.1016/j.compgeo.2007.12.003.

Soroush, A., Ferdowsi, B., 2011. Three dimensional discrete element modeling of granular media under cyclic constant volume loading: A micromechanical perspective. Powder Technol. 212 (1), 1-16. http://dx.doi.org/10.1016/j.powtec.2011.04.007.

Sukhunkitcharoen, R., Koseki, J., Otsubo, M., 2024. Relationship between void characteristics and re-liquefaction resistance: An image analysis study. Soils Found. 64 (6), 101506. http://dx.doi.org/10.1016/j.sandf.2024.101506.

Taboada, A., Chang, K.J., Radjaï, F., Bouchette, F., 2005. Rheology, force transmission, and shear instabilities in frictional granular media from biaxial numerical tests using the contact dynamics method. J. Geophys. Res.: Solid Earth 110 (B9), http://dx.doi.org/10.1029/2003JB002955.

Taghizadeh, K., Gaël, C., Luding, S., 2017. ALERT doctoral school 2017: Discrete element modeling.

Thevanayagam, S., Shenthan, T., Mohan, S., Liang, J., 2002. Undrained fragility of clean sands, silty sands, and sandy silt s. J. Geotech. Geoenvironmental Eng. 128 (10), 849-859. http://dx.doi.org/10.1061/(ASCE)1090-0241(2002)128:10(849).

Troncoso, J., Verdugo, R., 1985. Silt content and dynamic behavior of tailings sands. In: XI Int. Conf. on Soil Mechanics and Foundation Engineering, vol. 3, Buenos Aires, pp. 1311-1314.

Verdugo, R., 2024. Static liquefaction in the context of steady state/critical state and its application in the stability of tailings dams. Soil Dyn. Earthq. Eng. 176, 108270. http://dx.doi.org/10.1016/j.soildyn.2023.108270.

Verdugo, R., González, J., 2015. Liquefaction-induced ground damages during the 2010 Chile earthquake. Soil Dyn. Earthq. Eng. 79, 280-295. http://dx.doi.org/10.1016/j.soildyn.2015.04.016.

Verdugo, R., Ishihara, K., 1996. The steady state of sandy soils. Soils Found. 36 (2), 81-91. http://dx.doi.org/10.3208/sandf.36.2.81.

Voivret, C., Radjaï, F., Delenne, J.Y., El Youssoufi, M.S., 2007. Space-filling properties of polydisperse granular media. Phys. Rev. E 76, 021301. http://dx.doi.org/10.1103/PhysRevE.76.021301.

Wang, R., Fu, P., Zhang, J.M., Dafalias, Y.F., 2016. DEM study of fabric features governing undrained post-liquefaction shear deformation of sand. Acta Geotech. 11, 1321-1337. http://dx.doi.org/10.1007/s11440-016-0499-8.

Wang, F., Li, W., Yu, P., He, C., 2024. Research on the static liquefaction failure of upstream tailings dams under continuous discharge conditions. Comput. Geotech. 176, 106785. http://dx.doi.org/10.1016/j.compgeo.2024.106785.

Wang, F., Sassa, K., Wang, G., 2002. Mechanism of a long-runout landslide triggered by the August 1998 heavy rainfall in Fukushima prefecture, Japan. Eng. Geol. 63 (1), 169-185. http://dx.doi.org/10.1016/S0013-7952(01)00080-1.

Wang, G., Wei, J., 2016. Microstructure evolution of granular soils in cyclic mobility and post-liquefaction process. Granul. Matter 18, 51. http://dx.doi.org/10.1007/s10035-016-0621-5.

Wei, J., Huang, D., Wang, G., 2018. Microscale descriptors for particle-void distribution and jamming transition in pre- and post-liquefaction of granular soils. J. Eng. Mech. 144 (8), 04018067. http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0001482.

Weinhart, T., Orefice, L., Post, S., van Schrojenstein Lantman, M.P., Denissen, I.F.C., Tunuguntla, D.R.V., Tsang, J., Cheng, X., Shaheen, A., Shi, H., Rapino, S., Grannonio, R., Losacco, G., Barbosa, R., Jing, L., Alvarez Naranjo, C., Roy, A., den Otter, W.K., Thornton, A.R., 2019. Fast, flexible particle simulations — An introduction to MercuryDPM. Comput. Phys. Comm. 249, 107129. http://dx.doi.org/10.1016/j.cpc.2019.107129.
Yamamuro, J.A., Lade, P.V., 1997. Static liquefaction of very loose sands. Can. Geotech. J. 34 (6), 905-917. http://dx.doi.org/10.1139/t97-057.

Yang, J., Luo, X., 2018. The critical state friction angle of granular materials: does it depend on grading? Acta Geotech. 13 (3), 535-547. http://dx.doi.org/10.1007/s11440-017-0581-x.

Yang, M., Taiebat, M., 2024. Exploring the role of fabric anisotropy in cyclic liquefaction resistance under non-hydrostatic consolidation: Insights from DEM analysis. E3S Web Conf. 544, 10011. http://dx.doi.org/10.1051/e3sconf/202454410011.

Yang, M., Taiebat, M., Mutabaruka, P., Radjaï, F., 2021. Evolution of granular materials under isochoric cyclic simple shearing. Phys. Rev. E 103, 032904. http://dx.doi.org/10.1103/PhysRevE.103.032904.

Yang, M., Taiebat, M., Radjaï, F., 2022. Liquefaction of granular materials in constant-volume cyclic shearing: Transition between solid-like and fluid-like states. Comput. Geotech. 148, 104800. http://dx.doi.org/10.1016/j.compgeo.2022.104800.

Yang, Z.X., Wu, Y., 2017. Critical state for anisotropic granular materials: A discrete element perspective. Int. J. Geomech. 17 (2), 04016054. http://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0000720.

Yazdani, E., Nguyen, A., Evans, T.M., 2022. Shear-induced instability of sand containing fines: Using the equivalent intergranular void ratio as a state variable. Int. J. Geomech. 22 (8), 04022121. http://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0002486.

Zeghal, M., El Shamy, U., 2004. A continuum-discrete hydromechanical analysis of granular deposit liquefaction. Int. J. Numer. Anal. Methods Geomech. 28 (14), 1361-1383. http://dx.doi.org/10.1002/nag.390.

Zhang, L., Evans, T.M., 2020. Investigation of initial static shear stress effects on liquefaction resistance using discrete element method simulations. Int. J. Geomech. 20 (7), 04020087. http://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0001720.

Zhang, A., Jiang, M., Wang, D., 2023. Effect of fabric anisotropy on the cyclic liquefaction of sands: Insight from DEM simulations. Comput. Geotech. 155, 105188. http://dx.doi.org/10.1016/j.compgeo.2022.105188.

Zhao, Y., Liu, Z., Liang, T., He, F., Zhan, L., Chen, Y., Ling, D., Wang, J., 2024. Soil fluidisation induced by fine particles migration: Insights from the shenzhen 2015 landfill landslide. Eng. Geol. 343, 107783. http://dx.doi.org/10.1016/j.enggeo.2024.107783.

Zuo, K., Gu, X., Zhang, J., Wang, R., 2023. Exploring packing density, critical state, and liquefaction resistance of sand-fines mixture using DEM. Comput. Geotech. 156, 105278. http://dx.doi.org/10.1016/j.compgeo.2023.105278.