# Microstructural study of liquefaction in highly polydisperse granular media

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

► To cite this version:

Carolina Castro-Malaver, Manuel Cárdenas-Barrantes, David Cantor, Mathieu Renouf, Carlos Ovalle, et al.. Microstructural study of liquefaction in highly polydisperse granular media. Powders & Grains 2025 - 10th International Conference on Micromechanics on Granular Media, Dec 2025, Goa, India. pp.07012, 10.1051/epjconf/202534007012. hal-05466090

HAL Id: hal-05466090

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

Submitted on 19 Jan 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
# Microstructural study of liquefaction in highly polydisperse granular media

Carolina Castro-Malaver<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,*, Manuel Cárdenas-Barrantes<fnref n="1" />,<fnref n="2" />, David Cantor<fnref n="1" />,<fnref n="2" />, Mathieu Renouf<fnref n="3" />, Carlos Ovalle<fnref n="1" />,<fnref n="2" />,**, and Emilien Azéma<fnref n="1" />,<fnref n="3" />,<fnref n="4" />

**Abstract.** During earthquakes, rapid loading on loose, water-saturated silty sands can lead to undrained (constant volume) conditions that induce high pore water pressures. This process, known as liquefaction in geotechnical engineering, involves a loss of stress in the solid phase (effective) and can result in structural failures, such as frequent mine tailings dam collapses. Understanding the particle-scale mechanisms behind liquefaction is crucial for predictive modeling. However, this aspect remains poorly explored due to experimental limitations. In this study, we use discrete element method (DEM) simulations on one highly polydisperse granular material to investigate liquefaction. Samples of varying density are prepared by removing different amounts of floating particles (rattlers) after consolidation. The samples are then sheared under constant volume to the critical state. The results show that loose samples lose all strength, medium-loose samples temporarily liquefy but regain strength at large strains, and denser samples do not liquefy and exhibit continued shear strain hardening. At the micro-mechanical scale, permanent liquefaction is linked to heterogeneous solid fraction distributions (macro-pores), while samples with uniformly distributed local solid fraction either resist liquefaction or recover from it.

## 1 Introduction

Mine tailings are usually generated as slurries composed of silt and sand, which are deposited into a Tailings Storage Facility (TSF) through hydraulic filling. In this state, tailings stay saturated and unconsolidated, potentially leading to disastrous mudflows. Numerous TSF failures were documented throughout the 20th century and continue to occur today [1, 2].

TSF failures occur primarily by liquefaction after rapid loading, such as during earthquakes. Saturated soil liquefaction is defined as an increase in pore water pressure that leads to vanishing effective stresses and, therefore, a loss of shear strength [3, 4]. This phenomenon occurs in loose soils that exhibit a contractive response when exposed to shear. Under rapid loads, the material behaves undrained (constant volume) and the contractive tendency is restricted, leading to increased pore pressure and ultimately triggering liquefaction [5].

With the aim of better understanding the particle-scale mechanisms that trigger soil liquefaction, in this study we use Discrete Element Method (DEM) simulations of highly polydisperse granular materials. We shear samples of varied density to identify loose liquefiable behavior and dense non-liquefiable cases, along with the transition between them. The analyses are presented in terms of macro-mechanical behavior and micro-structural descriptors.

## 2 Numerical method

We simulated static liquefaction in a granular medium using strain-controlled constant volume shear tests. The simulations were carried out with the Contact Dynamics (CD) method, as implemented in the open-source platform LMGC90 [6].

In its standard formulation, the CD method treats particles as perfectly rigid and enforces non-overlapping contact conditions. However, under constant-volume shear, large stress concentrations can locally violate the non-overlap constraint, leading to numerical instabilities. To address this, we employed a hybrid contact model that introduces a linear elastic repulsion force, allowing for minimal controlled overlaps. This elastic extension of CD follows the formalism introduced by [7], which integrates particle elasticity into the rigid-body framework without compromising its core assumptions. The repulsion stiffness is chosen to be significantly higher than the stress levels involved, thereby preserving the rigid-particle approximation while ensuring numerical robustness. The numerical procedure comprises three stages: (1) sample generation, (2) uniaxial consolidation, and (3) constant-volume shearing.

## 3 Sample preparation and numerical test

We built two-dimensional samples of disks with a high size-polydispersity. The size span parameter $S = (d_{max} –

<footnote n="1">¹Department of Civil, Geological and Mining Engineering, Polytechnique Montréal, Québec, Canada</footnote>

<footnote n="2">²Research Institute on Mines and Environment (RIME) UQAT-Polytechnique, Québec, Canada</footnote>

<footnote n="3">³LMGC, Université de Montpellier, Montpellier, France</footnote>

<footnote n="4">⁴Institut Universitaire de France, France</footnote>

<footnote>*e-mail: lady-carolina.castro-malaver@polymtl.ca</footnote>

<footnote>**e-mail: carlos.ovalle@polymtl.ca</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Initial arrangement of particles for loose ($\beta$ = 0.0) and dense sample ($\beta$ = 1.0) after uniaxial consolidation at $P_0$. Colors refers to number of contacts per particle ($N_c$); purple color indicates floating particles ($N_c$ ∈ [0, 1]). $v_x$ is the velocity of the shear-strain.</figcaption></figure>

$d_{min}$) / ($d_{max}$ + $d_{min}$) accounts for the range variability of the maximum ($d_{max}$) and minimum particle size ($d_{min}$). We set $S$ = 0.7. A total of 17000 particles were uniformly distributed in volume within this range of sizes. The rigid disk particles were deposited into a square container of height $h$ = 48$d_{max}$, following a potential-based protocol [8]. The inter-particle friction coefficient was set to $\mu$ = 0.4, and was kept constant throughout all stages.

Once the particles were deposited, an uniaxial consolidation pressure of $P$ = 10kPa was applied at the bottom and top boundaries of the sample. Consolidation continued until the density stabilized, defined as the ratio between two consecutive values of packing fraction lower than $\Delta\phi$ = $10^{-4}$, with the packing fraction $\phi$ = $V_{grains}/V_{total}$; $V_{grains}$ is the volume occupied by the particles and $V_{total}$ is the total volume of the sample.

To generate samples of varied solid fraction or solid density, a fraction of floating particles were randomly removed after the consolidation stage. Floating particles refer to particles with fewer than two contacts and, therefore, do not belong to the force network. Thus, the mechanical stability of the sample reached at the consolidation stage is unaffected by removing them. The parameter $\beta$ refers to the ratio of floating particles remaining in the sample over the original floating particles. $\beta$ = 0.0 means that no floating grains remain in the sample, representing the loosest sample. In contrast, $\beta$ = 1.0 is the densest one. We created a set of eleven samples by varying $\beta$ uniformly from 0 to 1. Each sample had the same initial stress-engaged microstructure but different packing fractions and initial number of floating grains. The packing fraction after consolidation was $\phi_0$ = 0.865, which corresponds to the sample with $\beta$ = 1.0. The densities reached after the removal process varied up to $\phi$ = 0.806 for the loosest samples ($\beta$ = 0.0). Figure 1 shows the loosest and densest samples created for this study.

Before shearing the consolidated samples, the lateral walls of the container were replaced with periodic boundaries, meaning that any particle reaching one boundary reappears at the opposite side, allowing large shear strain to ensure critical states. Additionally, the top and bot-

tom walls were roughened by fixing some of the previously deposited particles, ensuring that deformation did not localize at these boundaries but was instead homogeneously distributed throughout the sample. The constant volume tests were performed by fixing the position of the upper and bottom plates along the y-axis [9]. The samples were strain-controlled sheared by moving these fixed plates along the x-axis. The strain rate was set to $\gamma$ = $4.5 \times 10^{-5}$, which ensures an inertial number of $I$ = $2.5 \times 10^{-5}$ meaning quasi-static flow.

## 4 Macroscopic response

The macroscopic response of the granular assemblies is described in terms of the deviatoric stress, $q$ = ($\sigma'_1$ - $\sigma'_2$) / 2, and the effective mean stress, $p'$ = ($\sigma'_1$ + $\sigma'_2$) / 2. $\sigma'_1$ and $\sigma'_2$ are the principal stresses of the granular stress tensor, $\sigma'_{ij}$, defined as:

$$ \sigma'_{ij} = \frac{1}{V} \sum_{\forall c} f_i^c l_j^c, \qquad (1) $$

where $i$ and $j$ run over the $x$ and $y$ components, $f_i^c$ is the $i$-th component of the contact force at contact $c$, $l_j^c$ is the $j$-th component of the branch vector and $V$ is the sample volume. The branch vector, $l^c$, is the vector joining the particles' mass centers interacting at contact $c$ [10].

Figure 2 shows the macroscopic response for all of the samples tested. These results are presented as a function of the shear strain $\gamma$ = $\Delta x/h$, where $\Delta x$ is the deformation in the x direction. The granular assemblies exhibited three types of behaviors. First, a continuous increase of stress up to a critical strength (i.e. non-liquefaction) for samples with $\beta$ ≥ 0.7. Second, a drop followed by an increment of the stresses (i.e. temporary liquefaction) in sample $\beta$ = 0.6. And third, vanishing of $q$ and $p'$ (i.e. liquefaction) showed in samples with $\beta$ ≤ 0.5. This complete loss of strength under monotonic loading is known as static liquefaction.

The logarithmic scale in Figure 2 allows evidencing the early occurrence of liquefaction in loose samples ($\beta$ ≤ 0.5). In fact, $\beta$ ≤ 0.5 corresponds with a transition point on density, where the system is not able to recover or gain resistance and liquefy. Further increases in density ($\beta$) delayed the occurrence of liquefaction. Besides, it can be observed that the higher the density, the greater the strength at critical state for the cases that do not exhibit liquefaction ($\beta$ ≥ 0.6). All of the samples have reached their critical strength value at $\gamma$ ≥ 0.18 regardless of the density.

## 5 Microstructure and connectivity

The mechanical coordination number ($z_M$), as defined in Eq. 2, quantifies the connectivity in granular assemblies by averaging the number of contacts per engaged particles (i.e.,those in the force network). Thus, $z_M$ excludes floating particles, particles with 1 and 0 contacts ($N_{1,2}$). Two-dimensional assemblies of rigid and frictional particles achieve mechanical stability at $3 \le z_M \le 4$; values of $z_M$ below this range indicate mechanical instability [11].
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Deviatoric *q* and (b) mean stress *p'* normalized by the initial compaction pressure *P*₀ against the shear strain *γ* for all initial proportion of floating particles *β*.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Mechanical coordination number *z*<sub>M</sub> as function of shear strain *γ* in logarithmic scale for *S* = 0.7</figcaption></figure>

Therefore, static liquefaction may be assessed by analyzing the evolution of z<sub>M</sub> as a function of shear strain, as shown in Figure 3.

$$z_M = \frac{2N_c - N_1}{N_p - N_0 - N_1} \qquad (2)$$

Each sample began with the same z<sub>M</sub>, which remained nearly constant in the early shear stage, z<sub>M</sub>(γ ≤ 10<sup>-3</sup>) ≈ 3.5. This indicates their initial mechanical stability. Subsequently, z<sub>M</sub> gradually increased in non-liquefaction cases, whereas it suddenly dropped with the liquefaction cases. In the latter cases, z<sub>M</sub> remained below 2.5. The increase in β delayed the abrupt decline below z<sub>M</sub> ≤ 3.3. The temporary liquefaction cases exhibited a sudden decrease to z<sub>M</sub> ≤ 3.4 followed by a gradual recovery to values above 3.4. Note that the stress recovery of β = 0.5 sample did not extend enough to exceed z<sub>M</sub> ≤ 3.4, resulting in its remaining liquefied state. The instability point at z<sub>M</sub> ≈ 3.4 seems to represent the liquefaction threshold.

## 6 Local density

Despite the strength and connectivity loss at liquefaction, the mechanism of recovery from liquefaction remains unclear [12]. Therefore, we perform a spatial analysis to find any variations in the sample's density along the shear deformation. We discretized the sample space with a mesh

to compute the density in each of the cells. The mesh has Δx = Δy = 1.04d<sub>max</sub> cell sizes, considering that Δx ≥ d<sub>max</sub>. The local cell density (φ<sub>L</sub>) was computed for each time step to assess its evolution through the deformation. Figure 4 shows the distribution of the normalized local density (φ<sub>L</sub>/⟨φ<sub>L</sub>⟩) and its variability through the strain γ. Three cases were analyzed: liquefaction (β = 0.0), temporary liquefaction (β = 0.6), and non-liquefaction (β = 1.0).

In the non-liquefied case at the beginning of the shear, the local density distribution ranges around 0.8 ≤ φ<sub>L</sub>/⟨φ<sub>L</sub>⟩ ≤ 1.2. The distribution of φ<sub>L</sub> is similar over the strain, indicating an homogeneous density distribution during shearing. Similar behavior is exhibited by the temporal liquefied sample. In contrast, the distribution of φ<sub>L</sub> in the liquefied case shifted towards lower values, indicating the occurrence of lower density cells in the sample. The three cases follow a normal distribution function, with p-values=according to a Shapiro-Wilk test

Every case was further analyzed by plotting the standard deviation of the local density (STDφ<sub>L</sub>/⟨φ<sub>L</sub>⟩) against shear strain γ (Figure 4b)). The STD in the non-liquefied and temporary liquefied samples slightly varied through the shear, accounting for the density homogeneity in the sample. This differs from the liquefied case in which STD fluctuated, suggesting the density heterogeneity. In fact, the change in the STD for the liquefaction case is consistent with the decrease in the local density below the mean value ⟨φ<sub>L</sub>⟩. The above results indicate that the homogeneity preserved in the local density is crucial to recover from the liquefied state.

## 7 Conclusions

In this study, we systematically investigated the influence of material density on static liquefaction using two-dimensional numerical simulations based on an extended Contact Dynamics (CD) method. This adapted framework incorporates contact elasticity following the approach of [7], enabling us to model highly polydisperse granular assemblies and to apply constant-volume shear conditions that emulate undrained loading. Eleven samples were created with the same initial micro-structure and different initial packing fractions, which varied as a function of the remaining proportion of floating particles (β) in each sample. The samples with the same particle size distribution (high size-polydispersity) were sheared exceeding 30% of shear strain.

The samples exhibited three behaviors: liquefaction, temporary liquefaction or non-liquefaction. Those samples that did not liquefy exhibited high values of the deviatoric stress at large deformations; the higher the amount of initial floating particles, the higher the critical deviatoric stress. On the other hand, low initial proportion of floating grains (β ≤ 0.5) resulted in liquefaction. Temporary liquefaction took place in the sample with 60% of the initial floating grains.

By conducting a multi-scale analysis we observed that those samples underwent liquefaction exhibited an average number of contacts below 3, regardless of the initial
<figure><img src="image_5.png" /><figcaption>**Figure 4.** a) Probability distribution functions and b) standard deviation of local density $\phi_L$ normalized by the average local density $\langle\phi_L\rangle$ throughout shear-strain $\gamma$ for three cases: liquefaction ($\beta = 0.0$), temporary liquefaction ($\beta = 0.6$), and non-liquefaction ($\beta = 1.0$).</figcaption></figure>

floating grains. The liquefied cases displayed the average connectivity continuously decreasing reaching their ultimate state without particles connected, whereas the non-liquefied cases' connectivity remained constant. This suggests that under undrained conditions the inter-particle connectivity loss caused liquefaction.

We present a possible microstructural explanation of the macromechanical observations. Since the overall solid fraction within a sample remains constant during our shear tests, the causes of liquefaction should be tracked at the local density. As expected, the analyses show that non-liquefiable materials maintain a homogeneous local density throughout the volume. On the other hand, liquefied samples exhibit heterogeneous local densities, with the appearance of local macropores that hinder structural reorganization between particles, resulting in vanishing stresses. Surprisingly, if a liquefied sample can maintain a homogeneous local density distribution, further shear deformation will prompt structural reorganization followed by stress increasing, allowing recovery from the liquefied state. Although computationally expensive, future work could be addressed for polydispersity in size by systematically considering the S parameter, polydispersity in shape, along with tridimensional grains in order to compare the results with experimental work on real soils.

## References

[1] J.C. Santamarina, L.A. Torres-Cruz, R.C. Bachus, Why coal ash and tailings dam disasters occur, Science **364**, 526–528 (2019). 10.1126/science.aax1927

[2] R. Dobry, L. Álvarez, Seismic failures of chilean tailings dams, Journal of the Soil Mechanics and Foundations Division **93**, 237 (1967). 10.1061/JS-FEAQ.0001054

[3] R. Verdugo, K. Ishihara, The steady state of sandy soils, Soils and Foundations **36**, 81 (1996). https://doi.org/10.3208/sandf.36.2_81

[4] M. Jefferies, K. Been, Soil Liquefaction: A Critical State Approach, Second Edition (CRC Press, 2015), ISBN 9780429153914, http://dx.doi.org/10.1201/b19114

[5] A. Casagrande, Liquefaction and cyclic deformation of sands. a critical review., in Proceedings of the Fifth Pan-American Conference on Soil Mechanics and Foundation Engineering (Harvard Soil Mechanics Series, Buenos Aires, 1975), 88, p. 27

[6] F. Dubois, V. Acary, M. Jean, The contact dynamics method: A nonsmooth story, Comptes Rendus Mécanique **346**, 247 (2018), the legacy of Jean-Jacques Moreau in mechanics / L'héritage de Jean-Jacques Moreau en mécanique. https://doi.org/10.1016/j.crme.2017.12.009

[7] K. Krabbenhoft, J. Huang, M.V. da Silva, A.V. Lyamin, Granular contact dynamics with particle elasticity, Granular Matter **14**, 607 (2012). 10.1007/s10035-012-0360-1

[8] C. Voivret, F. Radjaï, J.Y. Delenne, M.S. El Youssefi, Space-filling properties of polydisperse granular media, Phys. Rev. E **76**, 021301 (2007). 10.1103/PhysRevE.76.021301

[9] M. Cárdenas-Barrantes, C. Ovalle, Multiscale insights into sliding surface liquefaction through dem simulations, Comput. Geotech. **183**, 107191 (2025). 10.1016/j.compgeo.2025.107191

[10] R.J. Bathurst, L. Rothenburg, Micromechanical aspects of isotropic granular assemblies with linear contact interactions, Journal of Applied Mechanics **55**, 17 (1988). 10.1115/1.3173626

[11] B. Andreotti, Y. Forterre, O. Pouliquen, Granular Media: Between Fluid and Solid (Cambridge University Press, 2013)

[12] R. Verdugo, Static liquefaction in the context of steady state/critical state and its application in the stability of tailings dams, Soil Dynamics and Earthquake Engineering **176**, 108270 (2024). https://doi.org/10.1016/j.soildyn.2023.108270