# A MICROMECHANICS-BASED CLASSIFICATION OF THE REGIMES DELINEATING THE BEHAVIOUR OF GAP-GRADED SOILS

Peter Adesina, Antoine Wautier, Nadia Benahmed

► To cite this version:

Peter Adesina, Antoine Wautier, Nadia Benahmed. A MICROMECHANICS-BASED CLASSIFICATION OF THE REGIMES DELINEATING THE BEHAVIOUR OF GAP-GRADED SOILS. Computers and Geotechnics, 2024, 168, pp.106165. 10.1016/j.compgeo.2024.106165. hal-04464980

HAL Id: hal-04464980

https://hal.inrae.fr/hal-04464980

Submitted on 19 Feb 2024

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 Attribution 4.0 International License
A MICROMECHANICS-BASED CLASSIFICATION OF THE REGIMES
DELINEATING THE BEHAVIOUR OF GAP-GRADED SOILS

Peter Adesina<fnref n="1" />*, Antoine Wautier<fnref n="2" />, Nadia Benahmed<fnref n="3" />

¹,<fnref n="2" />,³INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182

Aix-en-Provence, France

*Corresponding author. E-mail addresses: peter.adesina@inrae.fr, greatpetson121@gmail.com

(P. Adesina); antoine.wautier@inrae.fr (Antoine Wautier); nadia.benahmed@inrae.fr (Nadia Benahmed)

**Abstract**

This study presents a micromechanical evaluation of the regimes delineating the behaviour of gap-graded granular assemblies, using discrete element simulations. Dense and loose bimodal assemblies of different fines content were prepared and subjected to drained triaxial compression until the critical state was reached. The regimes delineating the behaviour of the assemblies were evaluated, characterised and their significance discussed. While two regimes demarcated by the threshold fines content were identified based on the analysis of the macroscale characteristics of the assemblies, up to four regimes were identified based on the contributions of the particle size fractions and contact types to the total mean stress. Contrary to previous studies according to which fines control the mechanical behaviour of gap-graded assemblies from the threshold fines content, $f_c^{th}$, we found that the fines do not play a primary role in stress transmission until beyond a significantly larger fines content, $f_c^{eq}$ (the equivalent fines content), which depends on density and stress state. Based on the correlation found between the critical state strength and the stress-based skeleton void ratio proposed in this study, we conclude that stress-based skeleton void ratio can be useful in understanding the mechanical response of gap-graded materials at the critical state.

**Keywords: Gap-graded assemblies, regimes, density, triaxial compression, DEM**

<footnote n="1">¹Post doctoral researcher, INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence, France. Email: peter.adesina@inrae.fr; greatpetson121@gmail.com</footnote>

<footnote n="2">²Research Fellow, INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence, France. Email: antoine.wautier@inrae.fr</footnote>

<footnote n="3">³Research Fellow, INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence, France. Email: nadia.benahmed@inrae.fr</footnote>
# A MICROMECHANICS-BASED CLASSIFICATION OF THE REGIMES DELINEATING THE BEHAVIOUR OF GAP-GRADED SOILS

Peter Adesina<fnref n="1" />*, Antoine Wautier<fnref n="2" />, Nadia Benahmed<fnref n="3" />

¹,<fnref n="2" />,³INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182

Aix-en-Provence, France

*Corresponding author. E-mail addresses: peter.adesina@inrae.fr, greatpetson121@gmail.com

(P. Adesina); antoine.wautier@inrae.fr (Antoine Wautier); nadia.benahmed@inrae.fr (Nadia Benahmed)

## Abstract

This study presents a micromechanical evaluation of the regimes delineating the behaviour of gap-graded granular assemblies, using discrete element simulations. Dense and loose bimodal assemblies of different fines content were prepared and subjected to drained triaxial compression until the critical state was reached. The regimes delineating the behaviour of the assemblies were evaluated, characterised and their significance discussed. While two regimes demarcated by the threshold fines content were identified based on the analysis of the macroscale characteristics of the assemblies, up to four regimes were identified based on the contributions of the particle size fractions and contact types to the total mean stress. Contrary to previous studies according to which fines control the mechanical behaviour of gap-graded assemblies from the threshold fines content, $f_c^{th}$, we found that the fines do not play a primary role in stress transmission until beyond a significantly larger fines content, $f_c^{eq}$ (the equivalent fines content), which depends on density and stress state. Based on the correlation found between the critical state strength and the stress-based skeleton void ratio proposed in this study, we conclude that stress-based skeleton void ratio can be useful in understanding the mechanical response of gap-graded materials at the critical state.

**Keywords: Gap-graded assemblies, regimes, density, triaxial compression, DEM**

<footnote n="1">¹Post doctoral researcher, INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence, France. Email: peter.adesina@inrae.fr; greatpetson121@gmail.com</footnote>

<footnote n="2">²Research Fellow, INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence, France. Email: antoine.wautier@inrae.fr</footnote>

<footnote n="3">³Research Fellow, INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence, France. Email: nadia.benahmed@inrae.fr</footnote>
31

# 1. Introduction

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

Soil mixtures containing a coarser and a finer fraction (as in silty sands, sandy gravels, etc.) are abundant in nature. It is well established that the mechanical behaviour of these mixtures depends on the proportions of the coarser and the finer fraction within the mixtures. However, there is no consensus in the literature on specific fine contents delineating different behaviour for gap-graded granular mixtures. For example, as a result of the conceptual analysis presented in experimental studies, it has been suggested that the fines play a primary role in the shear behaviour of gap-graded assemblies beyond a threshold fines content (Thevanayagam et al., 2002); and control the behaviour beyond a limiting fines content (Lade & Yamamuro, 1997; Salgado et al., 2000; Skempton & Brogan, 1994; Vallejo, 2001). The threshold fines content is usually related to geometric properties, and in particular the measure of the minimum void ratio (Cubrinovski & Ishihara, 2002; Yang et al., 2005; Zuo & Baudet, 2015). From the threshold fines content, the fines start to disperse the coarse particles. The limiting fines content corresponds to the situation where the fines completely disperse the coarse particles (floating grains in a fine matrix) (Skempton & Brogan, 1994; Thevanayagam et al., 2002). Since quantifying the stress or the dispersion of individual particles within a matrix is difficult in experiments, it is almost impossible to ascertain the delineating fines content of the conceptual models without the use of a micromechanical framework.

In this regard, the discrete element method (DEM) proposed by (Cundall & Strack, 1979) has been found an effective numerical tool for probing the micromechanics of granular materials. DEM offers the opportunity to evaluate the assertions from the conceptual models in earlier experimental studies where the stress transmitted by individual particles cannot be determined. In the earlier studies where DEM has been used to assess the behaviour of gap-graded assemblies, the approach adopted involved determining the proportion of the stress transmitted by each particle size fraction and contact types. However, no consensus has been reached on the specific fines content delineating different behaviour; in fact, the existence of some of the specific fines contents suggested by experimental studies is sometimes questioned (as in Sufian et al., (2021)) or not identified.

In this study, we conducted four analyses organised into four sections to evaluate the regimes<fnref n="4" /> delineating the behaviour of gap-graded sand-silt mixtures. In Section 3, we determined the macroscale and micromechanical behaviour of the assemblies of sand-silt mixtures with fines

<footnote n="4">⁴ The word “regime” in this paper refers to zones of unique behaviour within the range of fines content in gap-graded assemblies.</footnote>
# 1. Introduction

Soil mixtures containing a coarser and a finer fraction (as in silty sands, sandy gravels, etc.) are abundant in nature. It is well established that the mechanical behaviour of these mixtures depends on the proportions of the coarser and the finer fraction within the mixtures. However, there is no consensus in the literature on specific fine contents delineating different behaviour for gap-graded granular mixtures. For example, as a result of the conceptual analysis presented in experimental studies, it has been suggested that the fines play a primary role in the shear behaviour of gap-graded assemblies beyond a threshold fines content (Thevanayagam et al., 2002); and control the behaviour beyond a limiting fines content (Lade & Yamamuro, 1997; Salgado et al., 2000; Skempton & Brogan, 1994; Vallejo, 2001). The threshold fines content is usually related to geometric properties, and in particular the measure of the minimum void ratio (Cubrinovski & Ishihara, 2002; Yang et al., 2005; Zuo & Baudet, 2015). From the threshold fines content, the fines start to disperse the coarse particles. The limiting fines content corresponds to the situation where the fines completely disperse the coarse particles (floating grains in a fine matrix) (Skempton & Brogan, 1994; Thevanayagam et al., 2002). Since quantifying the stress or the dispersion of individual particles within a matrix is difficult in experiments, it is almost impossible to ascertain the delineating fines content of the conceptual models without the use of a micromechanical framework.

In this regard, the discrete element method (DEM) proposed by (Cundall & Strack, 1979) has been found an effective numerical tool for probing the micromechanics of granular materials. DEM offers the opportunity to evaluate the assertions from the conceptual models in earlier experimental studies where the stress transmitted by individual particles cannot be determined. In the earlier studies where DEM has been used to assess the behaviour of gap-graded assemblies, the approach adopted involved determining the proportion of the stress transmitted by each particle size fraction and contact types. However, no consensus has been reached on the specific fines content delineating different behaviour; in fact, the existence of some of the specific fines contents suggested by experimental studies is sometimes questioned (as in Sufian et al., (2021)) or not identified.

In this study, we conducted four analyses organised into four sections to evaluate the regimes<fnref n="4" /> delineating the behaviour of gap-graded sand-silt mixtures. In Section 3, we determined the macroscale and micromechanical behaviour of the assemblies of sand-silt mixtures with fines

<footnote n="4">⁴ The word “regime” in this paper refers to zones of unique behaviour within the range of fines content in gap-graded assemblies.</footnote>
content, $f_c$, within 10%$\le f_c \le$70%; and assessed the existence of regimes delineating the
characteristics of the assemblies. The macroscale characteristics considered are the void ratio,
strength and dilatancy, while the micromechanical characteristic considered is the coordination
number. We then conducted in Section 4, a particle-scale analysis of the contribution of the
finer and the coarser fractions to the total mean stress transmitted by an assembly. This enabled
a particle stress-based evaluation of the regimes delineating the behaviour of the assemblies.
Following this, in Section 5, we conducted a contact-scale analysis of the contribution of each
contact type to the total mean stress, again to identify existing regimes. Finally, in Section 6,
we present a stress-based skeleton void ratio as an alternative void ratio index to interpret the
strength properties of the granular mixtures studied. We assessed the performance of alternative
void ratio indexes such as the mechanical and the void ratio (global), in interpreting the strength
exhibited by the granular mixtures studied, at both the peak and the critical state. An assessment
of the internal instability of the gap-graded assemblies (i.e. their susceptibility to internal
erosion) is beyond the scope of this study. However, the implications of the findings with
respect to internal erosion are discussed in the conclusion.

77

## 2. Numerical Simulation approach

The DEM simulations in this study were conducted using the open-source code YADE
(Šmilauer et al, 2021). All assemblies were generated to have a bimodal grading of fine
particles with diameter, $D_f$=0.375 mm and coarse particles with diameter $D_c$=3.075 mm such
that the size ratio, $\lambda = D_c/D_f$ = 8.2. The bimodal grading is considered the simplest type of gap-
graded material and was employed in this study to ensure the results obtained are solely from
the interaction between the coarser and the finer fractions without the interference of the
grading effect from within either of the fractions. The size ratio value was selected sufficiently
large to allow for fine migration in the pore space (should water flow be considered) (Lade et
al., 1998; Rahman et al., 2008; Shire et al., 2016; Thevanayagam et al., 2002).

After a parametric study to determine the representative element volume (REV), each assembly
contains 100,000 particles (Fig.1). Increasing the sample size to 200,000 particles for the $f_c$ =
30% case at both the dense and loose states did not significantly influence the stress-strain
responses (see Appendix A). The slight sensitivity to sample size observed in the volumetric
strains was also reported in the DEM study on REV for granular materials by Adesina et al.,
(2022) and can be attributed to the small variation in the initial void ratio of the assemblies.
content, $f_c$, within 10%$\le f_c \le$70%; and assessed the existence of regimes delineating the characteristics of the assemblies. The macroscale characteristics considered are the void ratio, strength and dilatancy, while the micromechanical characteristic considered is the coordination number. We then conducted in Section 4, a particle-scale analysis of the contribution of the finer and the coarser fractions to the total mean stress transmitted by an assembly. This enabled a particle stress-based evaluation of the regimes delineating the behaviour of the assemblies. Following this, in Section 5, we conducted a contact-scale analysis of the contribution of each contact type to the total mean stress, again to identify existing regimes. Finally, in Section 6, we present a stress-based skeleton void ratio as an alternative void ratio index to interpret the strength properties of the granular mixtures studied. We assessed the performance of alternative void ratio indexes such as the mechanical and the void ratio (global), in interpreting the strength exhibited by the granular mixtures studied, at both the peak and the critical state. An assessment of the internal instability of the gap-graded assemblies (i.e. their susceptibility to internal erosion) is beyond the scope of this study. However, the implications of the findings with respect to internal erosion are discussed in the conclusion.

## 2. Numerical Simulation approach

The DEM simulations in this study were conducted using the open-source code YADE (Šmilauer et al., 2021). All assemblies were generated to have a bimodal grading of fine particles with diameter, $D_f$=0.375 mm and coarse particles with diameter $D_c$=3.075 mm such that the size ratio, $\lambda = D_c/D_f$ = 8.2. The bimodal grading is considered the simplest type of gap-graded material and was employed in this study to ensure the results obtained are solely from the interaction between the coarser and the finer fractions without the interference of the grading effect from within either of the fractions. The size ratio value was selected sufficiently large to allow for fine migration in the pore space (should water flow be considered) (Lade et al., 1998; Rahman et al., 2008; Shire et al., 2016; Thevanayagam et al., 2002).

After a parametric study to determine the representative element volume (REV), each assembly contains 100,000 particles (Fig.1). Increasing the sample size to 200,000 particles for the $f_c$ = 30% case at both the dense and loose states did not significantly influence the stress-strain responses (see Appendix A). The slight sensitivity to sample size observed in the volumetric strains was also reported in the DEM study on REV for granular materials by Adesina et al., (2022) and can be attributed to the small variation in the initial void ratio of the assemblies.
Table 2 shows the number of the fines and the coarse particles for each percentage of fines content, $f_c$ (%).

The assemblies studied here were prepared by subjecting them to isotropic compression at 100 kPa within six frictionless walls. The standard linear elasto-plastic model was employed to simulate the interactions between particles. Table 1 shows the simulation parameters employed in the simulations. Following Jiang et al., (2018), the material density in the simulations was scaled to 1000 times the original value in order to increase the critical time step, $\Delta t_{cr} \propto r\sqrt{\rho/E}$, where r is the particle radius, $\rho$ is the density and E is the Young's modulus. A fraction (0.9$\Delta t_{cr}$) of this critical time step was then applied to ensure numerical stability. While this density scaling approach was adopted for computational efficiency, it does not compromise the quasi-static conditions of the simulations as long as the inertial number remains sufficiently small (Thornton, 2000; Thornton & Anthony, 1998). Based on the standard practice adopted in prior DEM studies including (Thornton, 2000) for assembly preparation, here, assemblies of different initial densities were generated by using a friction coefficient $\mu=0.03$ for the dense assemblies and $\mu=0.5$ for the loose assemblies, during isotropic compression. The void ratio of the dense assemblies is denoted $e_{min}$, while the void ratio of the loose assemblies is denoted as $e_{max}$. It is important to note that these $e_{min}$ and $e_{max}$ values cannot be directly mapped to the values obtained using standard procedures employed in experiments.

After isotropic compression, the assemblies were subjected to triaxial compression under a constant confining pressure of $\sigma_{xx}=\sigma_{yy}=100 \text{ kPa}^5$ in both the x-direction and y-direction and a strain rate of $0.01 \text{ s}^{-1}$ in the z-direction. A numerical damping coefficient of 0.05 is used. These conditions ensured that for our simulations, the inertial number for the coarse particles, $I_{dmax} = 1.60 \times 10^{-4}$, and that for the fines, $I_{dmin} = 1.95 \times 10^{-5}$, therefore the simulations were deemed quasi-static (da Cruz et al., 2005). Prior to shearing, additional number of cycles at $\mu=0.5$ was applied to equilibrate the assemblies. Also, the unbalanced force ratio$^6$ was checked to remain below a limit of 0.01 to guarantee static equilibrium. All assemblies were sheared at $\mu=0.5$ regardless of the friction coefficient used during the assembly preparation stage. A total of 22 shearing simulations were conducted in this study (11 fines contents are considered for 2 relative densities) (Table 2).

<footnote>$^5$Throughout this paper, soil mechanics conventions are used with compression being positive and extension being negative.</footnote>

<footnote>$^6$The maximum resultant force on the different grains divided by the mean contact force.</footnote>
Table 2 shows the number of the fines and the coarse particles for each percentage of fines content, $f_c$ (%).

The assemblies studied here were prepared by subjecting them to isotropic compression at 100 kPa within six frictionless walls. The standard linear elasto-plastic model was employed to simulate the interactions between particles. Table 1 shows the simulation parameters employed in the simulations. Following Jiang et al., (2018), the material density in the simulations was scaled to 1000 times the original value in order to increase the critical time step, $\Delta t_{cr} \propto r\sqrt{\rho/E}$, where r is the particle radius, $\rho$ is the density and E is the Young's modulus. A fraction (0.9$\Delta t_{cr}$) of this critical time step was then applied to ensure numerical stability. While this density scaling approach was adopted for computational efficiency, it does not compromise the quasi-static conditions of the simulations as long as the inertial number remains sufficiently small (Thornton, 2000; Thornton & Anthony, 1998). Based on the standard practice adopted in prior DEM studies including (Thornton, 2000) for assembly preparation, here, assemblies of different initial densities were generated by using a friction coefficient $\mu=0.03$ for the dense assemblies and $\mu=0.5$ for the loose assemblies, during isotropic compression. The void ratio of the dense assemblies is denoted $e_{min}$, while the void ratio of the loose assemblies is denoted as $e_{max}$. It is important to note that these $e_{min}$ and $e_{max}$ values cannot be directly mapped to the values obtained using standard procedures employed in experiments.

After isotropic compression, the assemblies were subjected to triaxial compression under a constant confining pressure of $\sigma_{xx}=\sigma_{yy}=100 \text{ kPa}^5$ in both the x-direction and y-direction and a strain rate of $0.01 \text{ s}^{-1}$ in the z-direction. A numerical damping coefficient of 0.05 is used. These conditions ensured that for our simulations, the inertial number for the coarse particles, $I_{dmax} = 1.60 \times 10^{-4}$, and that for the fines, $I_{dmin} = 1.95 \times 10^{-5}$, therefore the simulations were deemed quasi-static (da Cruz et al., 2005). Prior to shearing, additional number of cycles at $\mu=0.5$ was applied to equilibrate the assemblies. Also, the unbalanced force ratio$^6$ was checked to remain below a limit of 0.01 to guarantee static equilibrium. All assemblies were sheared at $\mu=0.5$ regardless of the friction coefficient used during the assembly preparation stage. A total of 22 shearing simulations were conducted in this study (11 fines contents are considered for 2 relative densities) (Table 2).

<footnote>$^5$Throughout this paper, soil mechanics conventions are used with compression being positive and extension being negative.</footnote>

<footnote>$^6$The maximum resultant force on the different grains divided by the mean contact force.</footnote>
<table><caption>125 **Table 1** Simulation parameters</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Diameter of coarse particles, D<sub>c</sub>, mm</td><td>3.075</td></tr><tr><td>Diameter of fine particles, D<sub>f</sub>, mm</td><td>0.375</td></tr><tr><td>Particle size ratio D<sub>c</sub> / D<sub>f</sub></td><td>8.2</td></tr><tr><td>Particle density, ρ, kg/m<sup>3</sup></td><td>2700 × 10<sup>3</sup></td></tr><tr><td>Contact law</td><td>Elasto-frictional</td></tr><tr><td>Inter-particle friction coefficient during isotropic compression</td><td>0.03 (dense), 0.5 (loose)</td></tr><tr><td>Inter-particle friction coefficient during shearing</td><td>0.5</td></tr><tr><td>Wall-particle friction coefficient</td><td>0</td></tr><tr><td>Particle normal stiffness, k<sub>n</sub>, N/m<sup>2</sup></td><td>3.56 × 10<sup>8</sup></td></tr><tr><td>Stiffness ratio</td><td>1.0</td></tr></tbody></table>

126

<table><caption>127 **Table 2** Number of particles for each assembly</caption><thead><tr><th rowspan="2">Specimens</th><th rowspan="2"></th><th colspan="10">Percentage of fine particles by weight</th></tr><tr><th>10%</th><th>15%</th><th>20%</th><th>25%</th><th>30%</th><th>35%</th><th>40%</th><th>45%</th><th>50%</th><th>60%</th><th>70%</th></tr></thead><tbody><tr><td rowspan="2">Number of particles</td><td>Coarse</td><td>1,612</td><td>1,021</td><td>723</td><td>543</td><td>423</td><td>337</td><td>272</td><td>222</td><td>182</td><td>121</td><td>78</td></tr><tr><td>Fine</td><td>98,388</td><td>98,979</td><td>99,277</td><td>99,457</td><td>99,577</td><td>99,663</td><td>99,728</td><td>99,778</td><td>99,818</td><td>99,879</td><td>99,922</td></tr><tr><td colspan="2">Total</td><td></td><td></td><td></td><td></td><td></td><td></td><td>100,000</td><td></td><td></td><td></td><td></td></tr></tbody></table>

128

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

129

130 **Fig. 1.** Dense specimens with different percentages of fines content: (a) f<sub>c</sub> = 10 % (b) f<sub>c</sub> = 30 %

131 (c) f<sub>c</sub> = 50 % (d) f<sub>c</sub> = 70 %

132
<table><caption>125 **Table 1** Simulation parameters</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Diameter of coarse particles, D<sub>c</sub>, mm</td><td>3.075</td></tr><tr><td>Diameter of fine particles, D<sub>f</sub>, mm</td><td>0.375</td></tr><tr><td>Particle size ratio D<sub>c</sub> / D<sub>f</sub></td><td>8.2</td></tr><tr><td>Particle density, ρ, kg/m<sup>3</sup></td><td>2700 × 10<sup>3</sup></td></tr><tr><td>Contact law</td><td>Elasto-frictional</td></tr><tr><td>Inter-particle friction coefficient during isotropic compression</td><td>0.03 (dense), 0.5 (loose)</td></tr><tr><td>Inter-particle friction coefficient during shearing</td><td>0.5</td></tr><tr><td>Wall-particle friction coefficient</td><td>0</td></tr><tr><td>Particle normal stiffness, k<sub>n</sub>, N/m<sup>2</sup></td><td>3.56 × 10<sup>8</sup></td></tr><tr><td>Stiffness ratio</td><td>1.0</td></tr></tbody></table>

126

<table><caption>127 **Table 2** Number of particles for each assembly</caption><thead><tr><th rowspan="2">Specimens</th><th rowspan="2"></th><th colspan="10">Percentage of fine particles by weight</th></tr><tr><th>10%</th><th>15%</th><th>20%</th><th>25%</th><th>30%</th><th>35%</th><th>40%</th><th>45%</th><th>50%</th><th>60%</th><th>70%</th></tr></thead><tbody><tr><td rowspan="2">Number of particles</td><td>Coarse</td><td>1,612</td><td>1,021</td><td>723</td><td>543</td><td>423</td><td>337</td><td>272</td><td>222</td><td>182</td><td>121</td><td>78</td></tr><tr><td>Fine</td><td>98,388</td><td>98,979</td><td>99,277</td><td>99,457</td><td>99,577</td><td>99,663</td><td>99,728</td><td>99,778</td><td>99,818</td><td>99,879</td><td>99,922</td></tr><tr><td colspan="2">Total</td><td></td><td></td><td></td><td></td><td></td><td></td><td>100,000</td><td></td><td></td><td></td><td></td></tr></tbody></table>

128

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

129

130 **Fig. 1.** Dense specimens with different percentages of fines content: (a) f<sub>c</sub> = 10 % (b) f<sub>c</sub> = 30 %

131 (c) f<sub>c</sub> = 50 % (d) f<sub>c</sub> = 70 %

132
### 3. Macro-scale and micromechanical characteristics of the binary mixtures for varying fines content.

The initial void ratios of the gap-graded assemblies with different fines content are plotted in Fig. 2a. As earlier mentioned, the void ratio of the dense and loose assemblies studied are denoted as $e_{min}$ and $e_{max}$, respectively. In both the dense and the loose assemblies, two regimes of distinct trends are observed. The first regime is characterised by a reduction in the void ratio as $f_c$ increases. This is as a result of the fines progressively filling the voids within the coarse particles. At a certain fines content referred to as the critical fines content (Skempton & Brogan, 1994), the threshold fines content, $f_c^{th}$ (Lade et al., 1998; Thevanayagam et al., 2002), or the transitional fines content (Yang et al., 2005), the voids within the coarse particles are fully filled by the fines. The second regime corresponds to $f_c>f_c^{th}$ and it is characterised by an increase in the void ratio as the $f_c$ increases. In this regime, the fines progressively disperse the coarse particles, thereby causing an increase in the void ratio. These two regimes have been reported in previous studies involving gap graded materials (Cubrinovski & Ishihara, 2002; Kuerbis, 1989; Lade & Yamamuro, 1997; W. Li et al., 2023; Y. Li et al., 2022; Minh et al., 2014; Sufian et al., 2021; Vallejo, 2001; Zuo & Baudet, 2015). The two identified regimes are referred to as the underfilled and the overfilled categories by dam engineers (ICOLD, 2013; Shire et al., 2014), and are demarcated by the threshold fines content (i.e. the filled state).

One particular feature that is less frequently highlighted in the literature is the fact that the $f_c^{th}$ slightly depends on the relative density of the mixture. The $f_c^{th}$ here was attained at $f_c = 0.30$ for the dense assemblies whereas $f_c^{th} = 0.35$ for the loose assemblies. Indeed, in looser assemblies, the pores are larger and more fine grains are required to fill the voids before the coarse grains are dispersed. This is in agreement with the 3D DEM study of binary mixtures by Minh et al., (2014). In the literature, different limits within which the threshold fines content is attainable have been reported. For example, Skempton & Brogan (1994) suggested that in practice, $f_c^{th}$ is unlikely to occur beyond 24% and 29% for dense and loose packings of sandy-gravels, respectively. Lade & Yamamuro (1997) suggested a limit of 20%-30% fines content for sand-silt mixtures of different gradations. Other studies show that the $f_c^{th}$ can occur outside Skempton and Brogan's and Lade & Yamamuro's limits, depending on the grain size distribution and the particle shape considered (Evans & Zhou, 1995; Sarkar et al., 2020; Shire et al., 2014; Sufian et al., 2021; Wang et al., 2022; Zuo & Baudet, 2015).
### 3. Macro-scale and micromechanical characteristics of the binary mixtures for varying fines content.

The initial void ratios of the gap-graded assemblies with different fines content are plotted in Fig. 2a. As earlier mentioned, the void ratio of the dense and loose assemblies studied are denoted as $e_{min}$ and $e_{max}$, respectively. In both the dense and the loose assemblies, two regimes of distinct trends are observed. The first regime is characterised by a reduction in the void ratio as $f_c$ increases. This is as a result of the fines progressively filling the voids within the coarse particles. At a certain fines content referred to as the critical fines content (Skempton & Brogan, 1994), the threshold fines content, $f_c^{th}$ (Lade et al., 1998; Thevanayagam et al., 2002), or the transitional fines content (Yang et al., 2005), the voids within the coarse particles are fully filled by the fines. The second regime corresponds to $f_c>f_c^{th}$ and it is characterised by an increase in the void ratio as the $f_c$ increases. In this regime, the fines progressively disperse the coarse particles, thereby causing an increase in the void ratio. These two regimes have been reported in previous studies involving gap graded materials (Cubrinovski & Ishihara, 2002; Kuerbis, 1989; Lade & Yamamuro, 1997; W. Li et al., 2023; Y. Li et al., 2022; Minh et al., 2014; Sufian et al., 2021; Vallejo, 2001; Zuo & Baudet, 2015). The two identified regimes are referred to as the underfilled and the overfilled categories by dam engineers (ICOLD, 2013; Shire et al., 2014), and are demarcated by the threshold fines content (i.e. the filled state).

One particular feature that is less frequently highlighted in the literature is the fact that the $f_c^{th}$ slightly depends on the relative density of the mixture. The $f_c^{th}$ here was attained at $f_c = 0.30$ for the dense assemblies whereas $f_c^{th} = 0.35$ for the loose assemblies. Indeed, in looser assemblies, the pores are larger and more fine grains are required to fill the voids before the coarse grains are dispersed. This is in agreement with the 3D DEM study of binary mixtures by Minh et al., (2014). In the literature, different limits within which the threshold fines content is attainable have been reported. For example, Skempton & Brogan (1994) suggested that in practice, $f_c^{th}$ is unlikely to occur beyond 24% and 29% for dense and loose packings of sandy-gravels, respectively. Lade & Yamamuro (1997) suggested a limit of 20%-30% fines content for sand-silt mixtures of different gradations. Other studies show that the $f_c^{th}$ can occur outside Skempton and Brogan's and Lade & Yamamuro's limits, depending on the grain size distribution and the particle shape considered (Evans & Zhou, 1995; Sarkar et al., 2020; Shire et al., 2014; Sufian et al., 2021; Wang et al., 2022; Zuo & Baudet, 2015).
The conceptual distinction between underfilled and overfilled assemblies led Yin et al., (2014) to propose the following equation for determining the void ratio, e, for sand-silt mixtures at different fines content $f_c$:

$$e = [e_{hc}(1 - f_c) + a f_c] \frac{1 - \tanh[\xi(f_c - f_{th})]}{2} + e_{hf} \left(f_c + \frac{1 - f_c}{(R_d)^m}\right) \frac{1 + \tanh[\xi(f_c - f_{th})]}{2} \quad (1)$$

where $e_{hc}$ is the void ratio of the pure sand, $e_{hf}$ is the void ratio of the pure silt, a is a material constant depending on the fabric structure of the soil mixture, $\xi$ is a material constant controlling the transition from a coarse grain matrix to a fine grain matrix, $f_{th}$ is the threshold fines content at which the coarse and fine grains contribute equally to the global void ratio, $R_d$ is the ratio of the mean size of the coarse grains, $D_{50}$, to the mean size of the fine grains $d_{50}$ and m ($0 < m < 1$) is a coefficient that depends on grain characteristics and fine grain packing. The fitting parameters for the dense assemblies are $e_{hc}=0.67; a = -0.64; \xi = 16.6; f_{th} = 0.28; e_{hf}=0.62; m=0.73$ while the parameters for the loose assemblies are $e_{hc}=0.71; a = 0.03; \xi = 12.6; f_{th} = 0.26; e_{hf}=0.76; m=0.99$. It is clear, as shown in Fig. 2a, that this phenomenological equation fits our data.

<figure><img src="image_2.png" /><figcaption>Fig. 2. Effect of fines content, $f_c$, on (a) the initial void ratio, e, (b) the range of attainable void ratios, $e_{max} - e_{min}$ for the assemblies of binary mixture studied</figcaption></figure>

Fig. 2b shows the range of attainable void ratio, $e_{max} - e_{min}$, for each fines content. $e_{max} - e_{min}$ for the underfilled assemblies (i.e. $f_c<f_c^{th} = 35\%$) is generally higher in comparison to the overfilled assemblies. In the underfilled assemblies, we observed an initial increase in $e_{max} - e_{min}$ with the fine content which indicates that the fines do not simply fill the voids in the loose case. Indeed, they enable the formation of larger pore structures.
The conceptual distinction between underfilled and overfilled assemblies led Yin et al., (2014) to propose the following equation for determining the void ratio, e, for sand-silt mixtures at different fines content $f_c$:

$$e = [e_{hc}(1 - f_c) + a f_c] \frac{1 - \tanh[\xi(f_c - f_{th})]}{2} + e_{hf} \left(f_c + \frac{1 - f_c}{(R_d)^m}\right) \frac{1 + \tanh[\xi(f_c - f_{th})]}{2} \quad (1)$$

where $e_{hc}$ is the void ratio of the pure sand, $e_{hf}$ is the void ratio of the pure silt, a is a material constant depending on the fabric structure of the soil mixture, $\xi$ is a material constant controlling the transition from a coarse grain matrix to a fine grain matrix, $f_{th}$ is the threshold fines content at which the coarse and fine grains contribute equally to the global void ratio, $R_d$ is the ratio of the mean size of the coarse grains, $D_{50}$, to the mean size of the fine grains $d_{50}$ and m ($0 < m < 1$) is a coefficient that depends on grain characteristics and fine grain packing. The fitting parameters for the dense assemblies are $e_{hc}=0.67;a = -0.64; \xi = 16.6; f_{th} = 0.28$; $e_{hf}=0.62;m=0.73$ while the parameters for the loose assemblies are $e_{hc}=0.71;a = 0.03; \xi = 12.6; f_{th} = 0.26; e_{hf}=0.76;m=0.99$. It is clear, as shown in Fig. 2a, that this phenomenological equation fits our data.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Effect of fines content, $f_c$, on (a) the initial void ratio, e, (b) the range of attainable void ratios, $e_{max} - e_{min}$ for the assemblies of binary mixture studied</figcaption></figure>

Fig. 2b shows the range of attainable void ratio, $e_{max} - e_{min}$, for each fines content. $e_{max} - e_{min}$ for the underfilled assemblies (i.e. $f_c<f_c^{th} = 35\%$) is generally higher in comparison to the overfilled assemblies. In the underfilled assemblies, we observed an initial increase in $e_{max} - e_{min}$ with the fine content which indicates that the fines do not simply fill the voids in the loose case. Indeed, they enable the formation of larger pore structures.
Fig. 3 shows the evolution of the stress ratio, *q*/*p*, and the volumetric strain (ε<sub>vol</sub>), against the axial strain (ε<sub>a</sub>), for the dense assemblies of binary mixtures with different fines content. The corresponding data for the loose assemblies are presented in Fig. 4. The assemblies were sheared until ε<sub>a</sub> ≈ 0.5 in order to reach the critical state. Here, the deviatoric stress, *q* = σ<sub>zz</sub> − (σ<sub>yy</sub> + σ<sub>xx</sub>)/2 and the mean stress, *p* = (σ<sub>zz</sub> + σ<sub>yy</sub> + σ<sub>xx</sub>)/3 (axisymmetric conditions). In a typical fashion to sand behaviour, the *q*/*p* for the initially dense assemblies increased until reaching a peak and thereafter softened to the critical state where the *q*/*p* fluctuates around a mean value, with more pronounced fluctuations observed for *f*<sub>c</sub> < 45% (Fig. 3a &b). The stiffness of the packings as influenced by the fines content is shown in the inset of Fig. 3a-c. These dense assemblies exhibited a dilative volumetric response (ε<sub>vol</sub> < 0) during shearing (Fig. 3d-f). At large strains (ε<sub>a</sub> > 35%), ε<sub>vol</sub> tend to increase from *f*<sub>c</sub> = 10% to *f*<sub>c</sub> = 35% (indicating a less dilative behaviour), and then decreased with further increase in the *f*<sub>c</sub>. In contrast, at small strains (ε<sub>a</sub> < 2%) (inset of Fig. 3d-f), a more dilative response was observed from *f*<sub>c</sub> = 10% to *f*<sub>c</sub><sup>th</sup>=30% as ε<sub>vol</sub> decreases. The stress ratio *q*/*p* for the initially loose assemblies (Fig. 4a-c) increased monotonically from the beginning of shearing to the critical state without significant softening observed, as expected. These loose assemblies exhibited a contractive response (ε<sub>vol</sub> > 0) during shearing (Fig. 4d-f).

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Effect of fines content, *f*<sub>c</sub>, on the stress-strain and the volumetric strain responses for dense assemblies (a,d) *f*<sub>c</sub> = 10%-25% (b,e) *f*<sub>c</sub> = 30%-40% (c,f) *f*<sub>c</sub> = 45%-70%. Inset graphs correspond to the start of the loading.</figcaption></figure>
Fig. 3 shows the evolution of the stress ratio, *q*/*p*, and the volumetric strain (ε<sub>vol</sub>), against the axial strain (ε<sub>a</sub>), for the dense assemblies of binary mixtures with different fines content. The corresponding data for the loose assemblies are presented in Fig. 4. The assemblies were sheared until ε<sub>a</sub> ≈ 0.5 in order to reach the critical state. Here, the deviatoric stress, *q* = σ<sub>zz</sub> - (σ<sub>yy</sub> + σ<sub>xx</sub>)/2 and the mean stress, *p* = (σ<sub>zz</sub> + σ<sub>yy</sub> + σ<sub>xx</sub>)/3 (axisymmetric conditions). In a typical fashion to sand behaviour, the *q*/*p* for the initially dense assemblies increased until reaching a peak and thereafter softened to the critical state where the *q*/*p* fluctuates around a mean value, with more pronounced fluctuations observed for *f*<sub>c</sub> < 45% (Fig. 3a &b). The stiffness of the packings as influenced by the fines content is shown in the inset of Fig. 3a-c. These dense assemblies exhibited a dilative volumetric response (ε<sub>vol</sub> < 0) during shearing (Fig. 3d-f). At large strains (ε<sub>a</sub> > 35%), ε<sub>vol</sub> tend to increase from *f*<sub>c</sub> = 10% to *f*<sub>c</sub> = 35% (indicating a less dilative behaviour), and then decreased with further increase in the *f*<sub>c</sub>. In contrast, at small strains (ε<sub>a</sub> < 2%) (inset of Fig. 3d-f), a more dilative response was observed from *f*<sub>c</sub> = 10% to *f*<sub>c</sub><sup>th</sup>=30% as ε<sub>vol</sub> decreases. The stress ratio *q*/*p* for the initially loose assemblies (Fig. 4a-c) increased monotonically from the beginning of shearing to the critical state without significant softening observed, as expected. These loose assemblies exhibited a contractive response (ε<sub>vol</sub> > 0) during shearing (Fig. 4d-f).

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Effect of fines content, *f*<sub>c</sub>, on the stress-strain and the volumetric strain responses for dense assemblies (a,d) *f*<sub>c</sub> = 10%-25% (b,e) *f*<sub>c</sub> = 30%-40% (c,f) *f*<sub>c</sub> = 45%-70%. Inset graphs correspond to the start of the loading.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 4. Effect of fines content, $f_c$, on the stress-strain and the volumetric strain responses for loose assemblies (a,d) $f_c$ = 10%-25% (b,e) $f_c$ = 30%-40% (c,f) $f_c$ = 45%-70%. Inset graphs correspond to the start of the loading.</figcaption></figure>

In Fig. 5a, the stress ratio $q/p$ at the peak, $(q/p)_{peak}$, is plotted against $f_c$ alongside the void ratio at the peak, $e_{peak}$, for comparison. Similarly, Fig. 5b shows the plot of the critical state stress ratio, $(q/p)_{crit}$, against $f_c$ alongside the critical state void ratio, $e_{crit}$. Here, the $(q/p)_{peak}$ is the maximum $q/p$ attained between the start of shearing and the critical state; $(q/p)_{crit}$ and $e_{crit}$ are mean values determined from the start of the critical state to the end of shearing. The critical state marks the state at which there is no noticeable change in the void ratio or the volume of a sample during shearing. This generally occurred at $\epsilon_a > 0.35$ for the studied assemblies. As the fines progressively fill the voids within the coarse particles in the underfilled regime where $f_c < f_c^{th}$, the assemblies become denser such that $(q/p)_{peak}$ increases with $f_c$ until the threshold fines content where the highest $(q/p)_{peak}$ is attained (Fig. 5a). For $f_c>f_c^{th}$, the fines disperse the coarse particles leaving more voids within the assemblies thereby causing a monotonic decrease in $(q/p)_{peak}$ as $f_c$ increased.

In agreement with earlier experimental studies on Ottawa sand and Carmague silty sand with $f_c \le 20\%$ under drained shearing (Benahmed et al., 2015; Chang & Yin, 2011; Salgado et al., 2000), a marginal increase in $(q/p)_{crit}$ was observed from $f_c$ = 10% to $f_c$ = 20%. This is followed by a significant decrease in $(q/p)_{crit}$ until $f_c$ = 40%, and finally a monotonic increase
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Effect of fines content, $f_c$, on the stress-strain and the volumetric strain responses for loose assemblies (a,d) $f_c$ = 10%-25% (b,e) $f_c$ = 30%-40% (c,f) $f_c$ = 45%-70%. Inset graphs correspond to the start of the loading.</figcaption></figure>

In Fig. 5a, the stress ratio $q/p$ at the peak, $(q/p)_{peak}$, is plotted against $f_c$ alongside the void ratio at the peak, $e_{peak}$, for comparison. Similarly, Fig. 5b shows the plot of the critical state stress ratio, $(q/p)_{crit}$, against $f_c$ alongside the critical state void ratio, $e_{crit}$. Here, the $(q/p)_{peak}$ is the maximum $q/p$ attained between the start of shearing and the critical state; $(q/p)_{crit}$ and $e_{crit}$ are mean values determined from the start of the critical state to the end of shearing. The critical state marks the state at which there is no noticeable change in the void ratio or the volume of a sample during shearing. This generally occurred at $\epsilon_a > 0.35$ for the studied assemblies. As the fines progressively fill the voids within the coarse particles in the underfilled regime where $f_c < f_c^{th}$, the assemblies become denser such that $(q/p)_{peak}$ increases with $f_c$ until the threshold fines content where the highest $(q/p)_{peak}$ is attained (Fig. 5a). For $f_c>f_c^{th}$, the fines disperse the coarse particles leaving more voids within the assemblies thereby causing a monotonic decrease in $(q/p)_{peak}$ as $f_c$ increased.

In agreement with earlier experimental studies on Ottawa sand and Carmague silty sand with $f_c \le 20\%$ under drained shearing (Benahmed et al., 2015; Chang & Yin, 2011; Salgado et al., 2000), a marginal increase in $(q/p)_{crit}$ was observed from $f_c$ = 10% to $f_c$ = 20%. This is followed by a significant decrease in $(q/p)_{crit}$ until $f_c$ = 40%, and finally a monotonic increase
229

until $f_c = 70\%$ (Fig. 5b). The $e_{crit}$ follows a similar trend as the $e_{peak}$, although the threshold fines content shifted from $f_c^{th}=30\%$ at the initial state (Fig. 5a) to $f_c^{th}=35\%$ at the critical state (Fig. 5b). While $e_{peak}$ could be used to explain the trend observed for $(q/p)_{peak}$, $(q/p)_{crit}$ does not correlate with $e_{crit}$. In Section 6, we provide a discussion on the use of an alternative void index, the stress-based skeleton void ratio, to explain the trend observed at the critical state. It is worthy of note that having a lot of fine grains filling the void space increases the peak strength but decreases the critical state strength. This might be interpreted in the capacity of fine grains to i) provide lateral support to force chains if the initial state is dense, but at the same time ii) act as ball bearings at critical state and ease shearing with limited deviatoric stress.

230

231

232

233

234

235

236

237

238

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

239

240

241

242

243

<caption>**Fig. 5** Effect of fines content on (a) stress ratio at peak, $(q/p)_{peak}$, for the dense assemblies, and peak void ratio, $e_{peak}$ (b) critical state stress ratio, $(q/p)_{crit}$, (unique value for dense and loose assemblies), and critical state void ratio, $e_{crit}$.</caption>

244

245

246

247

248

249

250

251

252

253

254

Fig. 6a shows the maximum dilatancy rate, $(-d\epsilon_v/d\epsilon_a)_{max}$ against $f_c$. In the underfilled regime, $(-d\epsilon_v/d\epsilon_a)_{max}$ for the dense assemblies increased from $f_c = 10\%$ until the threshold fines content $f_c^{th}=30\%$, as the fines progressively fill the void space (Fig. 6a). For $f_c>f_c^{th}$, $(-d\epsilon_v/d\epsilon_a)_{max}$ decreased monotonically with an increase in $f_c$. The increase in dilatancy observed within the underfilled regime was also reported for Ottawa sand having $f_c \le 20\%$, under drained triaxial shearing (Chang & Yin, 2011; Salgado et al., 2000); and was attributed to increased particle interlock as the fines occupy the voids within the coarse particles (Kuerbis, 1989; Salgado et al., 2000). Indeed, for the assemblies of binary mixture studied here, there is a correlation between the maximum dilatancy angle, $\psi_{max}$ and the initial void ratio, $e_{min}$, (Fig. 6b), where $\psi_{max} = \sin^{-1}[((-d\epsilon_v/d\epsilon_a)_{max})/(2 + (-d\epsilon_v/d\epsilon_a)_{max})]$ following (Xiao et al., 2017). A linear relationship was also established between the peak friction angle, $\phi_p$ and the
until $f_c = 70\%$ (Fig. 5b). The $e_{crit}$ follows a similar trend as the $e_{peak}$, although the threshold fines content shifted from $f_c^{th}=30\%$ at the initial state (Fig. 5a) to $f_c^{th}=35\%$ at the critical state (Fig. 5b). While $e_{peak}$ could be used to explain the trend observed for $(q/p)_{peak}$, $(q/p)_{crit}$ does not correlate with $e_{crit}$. In Section 6, we provide a discussion on the use of an alternative void index, the stress-based skeleton void ratio, to explain the trend observed at the critical state. It is worthy of note that having a lot of fine grains filling the void space increases the peak strength but decreases the critical state strength. This might be interpreted in the capacity of fine grains to i) provide lateral support to force chains if the initial state is dense, but at the same time ii) act as ball bearings at critical state and ease shearing with limited deviatoric stress.

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Effect of fines content on (a) stress ratio at peak, $(q/p)_{peak}$, for the dense assemblies, and peak void ratio, $e_{peak}$ (b) critical state stress ratio, $(q/p)_{crit}$, (unique value for dense and loose assemblies), and critical state void ratio, $e_{crit}$.</figcaption></figure>

Fig. 6a shows the maximum dilatancy rate, $(-d\epsilon_v/d\epsilon_a)_{max}$ against $f_c$. In the underfilled regime, $(-d\epsilon_v/d\epsilon_a)_{max}$ for the dense assemblies increased from $f_c = 10\%$ until the threshold fines content $f_c^{th}=30\%$, as the fines progressively fill the void space (Fig. 6a). For $f_c>f_c^{th}$, $(-d\epsilon_v/d\epsilon_a)_{max}$ decreased monotonically with an increase in $f_c$. The increase in dilatancy observed within the underfilled regime was also reported for Ottawa sand having $f_c \le 20\%$, under drained triaxial shearing (Chang & Yin, 2011; Salgado et al., 2000); and was attributed to increased particle interlock as the fines occupy the voids within the coarse particles (Kuerbis, 1989; Salgado et al., 2000). Indeed, for the assemblies of binary mixture studied here, there is a correlation between the maximum dilatancy angle, $\psi_{max}$ and the initial void ratio, $e_{min}$, (Fig. 6b), where $\psi_{max} = \sin^{-1}[((-d\epsilon_v/d\epsilon_a)_{max})/(2 + (-d\epsilon_v/d\epsilon_a)_{max})]$ following (Xiao et al., 2017). A linear relationship was also established between the peak friction angle, $\phi_p$ and the
$\psi_{max}$ (Fig. 6c), where the friction angle, $\phi = \sin^{-1}[(\sigma_{zz} - \sigma_{yy})/(\sigma_{zz} + \sigma_{yy})]$. This relationship has been found in earlier studies for clean sands (Adesina et al., 2024; Bolton, 1986; Vaid & Sasitharan, 1992), and was also reported for sand-silt mixtures having $f_c \le 20\%$ in the experimental study on sand with non-plastic fines (Xiao et al., 2017).

<figure><img src="image_6.png" /><figcaption>Fig. 6 (a) Effect of fines content, $f_c$, on the maximum dilatancy rate, $(-d\epsilon_v/d\epsilon_a)_{max}$. Relationship between (b) maximum dilatancy angle, $\psi_{max}$ and initial void ratio, $e_{min}$ (c) peak friction angle, $\phi_p$, and maximum dilatancy angle, $\psi_{max}$ for dense assemblies of binary mixture.</figcaption></figure>

The coordination number, i.e. the average number of contacts per particle in a granular system is the simplest measure of the connectivity of the contact network which is related to the structural stability of the system (Thornton, 2015). Here, we employ the mechanical coordination number as a micromechanical parameter to identify and characterise the regimes delineating the behaviour of gap-graded materials. The mechanical coordination number, $Z_m$, is defined as the average number of contacts per particle when the particles which do not contribute to force transmission (i.e., particles with less than two contacts which are referred to as rattlers) are excluded<fnref n="7" />. Thornton (2000) defined the $Z_m$ as:

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

where $N_c$ is the total number of contacts, $N_p$ is total number of particles and $N_0$ and $N_1$ are the number of particles with zero and one contact in the granular system, respectively. Figs. 7a-c show $Z_m$ for all particles, $Z_{mc}$ for the coarse particles only, and $Z_{mf}$ for the fine particles only, respectively, alongside the void ratio, $e$, at the initial, the peak and the critical states. While $Z_m$ considers all particles, $Z_{mc}$ and $Z_{mf}$ considers only the contacts made by the coarse particles and the fine particles, respectively<fnref n="8" />. Figs. 7d-e show the number of fines and coarse particles

<footnote n="7">⁷ Note that the coordination number of particles with more than two contacts accounts for the contacts from particles with only one contact in the standard definition of $Z_m$.</footnote>

<footnote n="8">⁸ The contacts made by the coarse particles include both the coarse-coarse and the coarse-fine contacts while the contacts made by the fine particles include both the fine-fine and the coarse-fine contacts.</footnote>
$\psi_{max}$ (Fig. 6c), where the friction angle, $\phi = \sin^{-1}[(\sigma_{zz} - \sigma_{yy})/(\sigma_{zz} + \sigma_{yy})]$. This relationship has been found in earlier studies for clean sands (Adesina et al., 2024; Bolton, 1986; Vaid & Sasitharan, 1992), and was also reported for sand-silt mixtures having $f_c \le 20\%$ in the experimental study on sand with non-plastic fines (Xiao et al., 2017).

<figure><img src="image_7.png" /><figcaption>Fig. 6 (a) Effect of fines content, $f_c$, on the maximum dilatancy rate, $(-d\epsilon_v/d\epsilon_a)_{max}$. Relationship between (b) maximum dilatancy angle, $\psi_{max}$ and initial void ratio, $e_{min}$ (c) peak friction angle, $\phi_p$, and maximum dilatancy angle, $\psi_{max}$ for dense assemblies of binary mixture.</figcaption></figure>

The coordination number, i.e. the average number of contacts per particle in a granular system is the simplest measure of the connectivity of the contact network which is related to the structural stability of the system (Thornton, 2015). Here, we employ the mechanical coordination number as a micromechanical parameter to identify and characterise the regimes delineating the behaviour of gap-graded materials. The mechanical coordination number, $Z_m$, is defined as the average number of contacts per particle when the particles which do not contribute to force transmission (i.e., particles with less than two contacts which are referred to as rattlers) are excluded<fnref n="7" />. Thornton (2000) defined the $Z_m$ as:

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

where $N_c$ is the total number of contacts, $N_p$ is total number of particles and $N_0$ and $N_1$ are the number of particles with zero and one contact in the granular system, respectively. Figs. 7a-c show $Z_m$ for all particles, $Z_{mc}$ for the coarse particles only, and $Z_{mf}$ for the fine particles only, respectively, alongside the void ratio, $e$, at the initial, the peak and the critical states. While $Z_m$ considers all particles, $Z_{mc}$ and $Z_{mf}$ considers only the contacts made by the coarse particles and the fine particles, respectively<fnref n="8" />. Figs. 7d-e show the number of fines and coarse particles

<footnote n="7">⁷ Note that the coordination number of particles with more than two contacts accounts for the contacts from particles with only one contact in the standard definition of $Z_m$.</footnote>

<footnote n="8">⁸ The contacts made by the coarse particles include both the coarse-coarse and the coarse-fine contacts while the contacts made by the fine particles include both the fine-fine and the coarse-fine contacts.</footnote>
that do not contribute to force transmission (i.e. the fines-rattlers and coarse-rattlers) at the initial, the peak and the critical states. Generally, two trends demarcated by the $f_c^{th}$ can be observed in Fig. 7a; a decrease in the $Z_m$ for $f_c<f_c^{th}$ (although the decrease is less significant for the dense assemblies at the initial state); and a monotonic increase for $f_c>f_c^{th}$. The trends observed here are similar to the data shown in the DEM study by Liu et al., (2022) for the $Z_m$ of sand-silt mixtures having particle size ratio, $\lambda = 8.4$ and 18.1, and at the initial states. Generally, an initial increase in the $f_c$ resulted in no change in the $Z_{mc}$ and the $Z_{mf}$ until $f_c^{th}$, after which a monotonic increase followed (Figs. 7b & c). The monotonic decrease in the number of fines-rattlers from $f_c^{th}$ (Fig. 7d) indicates an increase in the number of fines mobilised in force transmission as $f_c$ increases. As expected, the number of coarse-rattlers decreased with an increase in $f_c$ until $f_c^{th}$, beyond which no rattler exists among the coarse particles (Fig. 7e).

<figure><img src="image_7.png" /><figcaption>**Fig. 7** Effect of fines content on the (a) mechanical coordination number, $Z_m$, for all particles (b) $Z_{mc}$ for coarse particles (c) $Z_{mf}$ for fine particles (d) number of fines-rattlers (e) number of coarse-rattlers within assemblies of binary mixture.</figcaption></figure>

Fig. 8a shows the contact lost or gained during shearing, $\Delta Z$, determined here as $Z_{ini} - Z_{crit}$, where $Z_{ini}$ and $Z_{crit}$ are the coordination number at the initial states and at the critical state, respectively. Fig. 8b shows the excess friction angle, which indicates the strain softening exhibited by the assemblies. From these figures, we observe the emergence of three groups of
that do not contribute to force transmission (i.e. the fines-rattlers and coarse-rattlers) at the initial, the peak and the critical states. Generally, two trends demarcated by the $f_c^{th}$ can be observed in Fig. 7a; a decrease in the $Z_m$ for $f_c<f_c^{th}$ (although the decrease is less significant for the dense assemblies at the initial state); and a monotonic increase for $f_c>f_c^{th}$. The trends observed here are similar to the data shown in the DEM study by Liu et al., (2022) for the $Z_m$ of sand-silt mixtures having particle size ratio, $\lambda = 8.4$ and 18.1, and at the initial states. Generally, an initial increase in the $f_c$ resulted in no change in the $Z_{mc}$ and the $Z_{mf}$ until $f_c^{th}$, after which a monotonic increase followed (Figs. 7b & c). The monotonic decrease in the number of fines-rattlers from $f_c^{th}$ (Fig. 7d) indicates an increase in the number of fines mobilised in force transmission as $f_c$ increases. As expected, the number of coarse-rattlers decreased with an increase in $f_c$ until $f_c^{th}$, beyond which no rattler exists among the coarse particles (Fig. 7e).

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Effect of fines content on the (a) mechanical coordination number, $Z_m$, for all particles (b) $Z_{mc}$ for coarse particles (c) $Z_{mf}$ for fine particles (d) number of fines-rattlers (e) number of coarse-rattlers within assemblies of binary mixture.</figcaption></figure>

Fig. 8a shows the contact lost or gained during shearing, $\Delta Z$, determined here as $Z_{ini} - Z_{crit}$, where $Z_{ini}$ and $Z_{crit}$ are the coordination number at the initial states and at the critical state, respectively. Fig. 8b shows the excess friction angle, which indicates the strain softening exhibited by the assemblies. From these figures, we observe the emergence of three groups of
f_c values with similar ΔZ and φ_p-φ_c values. This suggests the existence of three regimes, the first being f_c < f_c^th where ΔZ ≈ 0 and strain softening is relatively minimal. In the second regime delineated by f_c^th ≤ f_c < 45%, ΔZ is highest, suggesting that the highest amount of contacts were lost within this regime during shearing. This is accompanied by significant softening occurring within the regime, in agreement with the DEM study by Adesina et al., (2024), where a linear relationship is established between the degree of softening and the contact loss during shearing by linearly-graded granular assemblies. The significant loss of contacts exhibited by the assemblies in this regime can be linked to the large fluctuations in their shearing responses at the critical state (see Figs. 3 & 4). It has been shown in the DEM study by Adesina et al., (2022) that the magnitude of fluctuations in the shearing responses of granular assemblies at the critical state increases with a decrease in the total number of contacts in the assemblies. This observation can also be related to the size of the REV, which depends on the fine contents. While the condition is largely met for small and large fine contents with 100,000 grains, this is less the case for intermediate fine contents in which local microstructure rearrangements reflect more in the macroscopic response. This is observed in Fig. 16 through the evolution of the boundary term in the stress decomposition. In Appendix A, we showed that our sample size of 100,000 particles can be considered as a REV.

In the third regime (i.e. $f_c \ge 45\%$), relatively minimal ΔZ is observed; the strain softening within this regime transitions from moderate to minimal. ΔZ values in the last two regimes are higher for dense assemblies in comparison to the loose assemblies, for the same $f_c$ value; ΔZ<0 for the loose assemblies in regime 3 indicating an ultimate contact gain. While Skempton & Brogan (1994) and Thevanayagam et al., (2002) provided a conceptual delineation of the regimes based on void ratio, here, we provide a mechanistic delineation of the regimes with our analysis of contact evolution and strain softening. In Table 3, we present a summary of the unique trends observed for each of the three regimes identified based on our analysis of the macromechanical and the micromechanical characteristics of the studied gap-graded assemblies.
f_c values with similar ΔZ and φ_p-φ_c values. This suggests the existence of three regimes, the first being f_c < f_c^th where ΔZ ≈ 0 and strain softening is relatively minimal. In the second regime delineated by f_c^th ≤ f_c < 45%, ΔZ is highest, suggesting that the highest amount of contacts were lost within this regime during shearing. This is accompanied by significant softening occurring within the regime, in agreement with the DEM study by Adesina et al., (2024), where a linear relationship is established between the degree of softening and the contact loss during shearing by linearly-graded granular assemblies. The significant loss of contacts exhibited by the assemblies in this regime can be linked to the large fluctuations in their shearing responses at the critical state (see Figs. 3 & 4). It has been shown in the DEM study by Adesina et al., (2022) that the magnitude of fluctuations in the shearing responses of granular assemblies at the critical state increases with a decrease in the total number of contacts in the assemblies. This observation can also be related to the size of the REV, which depends on the fine contents. While the condition is largely met for small and large fine contents with 100,000 grains, this is less the case for intermediate fine contents in which local microstructure rearrangements reflect more in the macroscopic response. This is observed in Fig. 16 through the evolution of the boundary term in the stress decomposition. In Appendix A, we showed that our sample size of 100,000 particles can be considered as a REV.

In the third regime (i.e. $f_c \ge 45\%$), relatively minimal ΔZ is observed; the strain softening within this regime transitions from moderate to minimal. ΔZ values in the last two regimes are higher for dense assemblies in comparison to the loose assemblies, for the same $f_c$ value; ΔZ<0 for the loose assemblies in regime 3 indicating an ultimate contact gain. While Skempton & Brogan (1994) and Thevanayagam et al., (2002) provided a conceptual delineation of the regimes based on void ratio, here, we provide a mechanistic delineation of the regimes with our analysis of contact evolution and strain softening. In Table 3, we present a summary of the unique trends observed for each of the three regimes identified based on our analysis of the macromechanical and the micromechanical characteristics of the studied gap-graded assemblies.
<figure><img src="image_8.png" /><figcaption>**Fig. 8** (a) Change in coordination number during shearing, ΔZ (filled circle for dense assemblies and open circles for loose assemblies). ΔZ>0 means contact lost and ΔZ<0 means contact gained during shearing. (b) excess friction angle (strain softening), $\phi_p-\phi_c$ for dense assemblies of binary mixture.</figcaption></figure>

<table><caption>**Table 3:** Summary of regime based on macromechanical and micromechanical characteristics</caption><thead><tr><th>Δ in parameters<br>as the f<sub>c</sub> ↑</th><th>Regime 1<br>(f<sub>c</sub> &lt; f<sub>c</sub><sup>th</sup>)</th><th>Regime 2<br>(f<sub>c</sub><sup>th</sup> &lt; f<sub>c</sub> &lt; 45%)</th><th>Regime 3<br>(f<sub>c</sub> &gt; 45%)</th></tr></thead><tbody><tr><td>e<sub>ini</sub>; e<sub>crit</sub></td><td>Decreases</td><td></td><td>Increases</td></tr><tr><td>(q/p)<sub>max</sub></td><td>Increases</td><td></td><td>Decreases</td></tr><tr><td>(-dε<sub>v</sub>/dε<sub>a</sub>)<sub>max</sub></td><td>Increases</td><td></td><td>Decreases</td></tr><tr><td>Z<sub>m</sub></td><td>Decreases</td><td>Increases</td><td>Increases marginally to a plateau</td></tr><tr><td>Z<sub>mc</sub>; Z<sub>mf</sub></td><td>No change</td><td>Increases</td><td>Increases marginally to a plateau</td></tr><tr><td>ΔZ</td><td>No change</td><td>Increases</td><td>Decreases marginally to a plateau</td></tr><tr><td>φ<sub>p</sub>-φ<sub>c</sub></td><td>Increases</td><td>Plateau</td><td>Decreases</td></tr></tbody></table>

Note: When the data in Fig. 8 are considered, the boundary between Regime 2 and Regime 3 range between $f_c \in [35\% ; 45\%]$ depending on density. This range will also depend on the particle size distribution considered.

## 4. Particle based stress-transmission

Using the Love-Weber particle stress tensor as described in Nicot et al., (2013), the mean stress per particle, $\sigma_{ij}^p$, is calculated as:

$$ \sigma_{ij}^p = \frac{1}{Vp} R_p \sum_{c=1}^{N_c^p} f_j^c \mathbf{n}_i \qquad (3) $$
<figure><img src="image_9.png" /><figcaption>**Fig. 8** (a) Change in coordination number during shearing, ΔZ (filled circle for dense assemblies and open circles for loose assemblies). ΔZ>0 means contact lost and ΔZ<0 means contact gained during shearing. (b) excess friction angle (strain softening), φ<sub>p</sub>-φ<sub>c</sub> for dense assemblies of binary mixture.</figcaption></figure>

<table><caption>**Table 3:** Summary of regime based on macromechanical and micromechanical characteristics</caption><thead><tr><th>Δ in parameters<br>as the f<sub>c</sub> ↑</th><th>Regime 1<br>(f<sub>c</sub> &lt; f<sub>c</sub><sup>th</sup>)</th><th>Regime 2<br>(f<sub>c</sub><sup>th</sup> &lt; f<sub>c</sub> &lt; 45%)</th><th>Regime 3<br>(f<sub>c</sub> &gt; 45%)</th></tr></thead><tbody><tr><td>e<sub>ini</sub>; e<sub>crit</sub></td><td>Decreases</td><td></td><td>Increases</td></tr><tr><td>(q/p)<sub>max</sub></td><td>Increases</td><td></td><td>Decreases</td></tr><tr><td>(-dε<sub>v</sub>/dε<sub>a</sub>)<sub>max</sub></td><td>Increases</td><td></td><td>Decreases</td></tr><tr><td>Z<sub>m</sub></td><td>Decreases</td><td>Increases</td><td>Increases marginally to a plateau</td></tr><tr><td>Z<sub>mc</sub>; Z<sub>mf</sub></td><td>No change</td><td>Increases</td><td>Increases marginally to a plateau</td></tr><tr><td>ΔZ</td><td>No change</td><td>Increases</td><td>Decreases marginally to a plateau</td></tr><tr><td>φ<sub>p</sub>-φ<sub>c</sub></td><td>Increases</td><td>Plateau</td><td>Decreases</td></tr></tbody></table>

Note: When the data in Fig. 8 are considered, the boundary between Regime 2 and Regime 3 range between f<sub>c</sub> ∈ [35%; 45%] depending on density. This range will also depend on the particle size distribution considered.

## 4. Particle based stress-transmission

Using the Love-Weber particle stress tensor as described in Nicot et al., (2013), the mean stress per particle, σ<sub>ij</sub><sup>p</sup>, is calculated as:

$$ \sigma_{ij}^{p} = \frac{1}{Vp} R_p \sum_{c=1}^{N_c^p} f_j^c \mathbf{n}_i \quad (3) $$
where the summation runs over all contacts c made by the particle p. $V^p$ is the particle volume, $R_p$ is the radius of the particle, $N_c^p$ is the number of contacts involving the particle, $f_j^c$ is the contact force vector and $n_i$ is the unit branch vector between two particles in contact. The stress tensor at the material point scale (i.e. for a large number of particles), $\sigma_{ij}$, is given as a weighted average of the particle stresses (by particle volume):

$$ \sigma_{ij} = \frac{1}{V} \sum_{p=1}^{N_p} V^p \sigma_{ij}^p \quad (4) $$

where V is the volume of the entire system and $N_p$ is the number of particles in the system. $\sigma_{ij}$ can be decomposed to the contribution of each particle category, i.e. fine and coarse, as:

$$ \sigma_{ij} = \frac{1}{V} \sum_{p=1}^{N_p^c} V^p \sigma_{ij}^p + \frac{1}{V} \sum_{p=1}^{N_p^f} V^p \sigma_{ij}^p \quad (5) $$

where $N_p^c$ and $N_p^f$ are the number of coarse and fine particles in the system, respectively<fnref n="9" />. The mean effective stress, p, for the entire system is calculated as:

$$ p = \frac{\sigma_{xx} + \sigma_{yy} + \sigma_{zz}}{3} \quad (6) $$

The proportion of the total mean stress transmitted by the finer fraction, $\alpha_f$, is calculated as $p^f/p$, and similarly, the proportion of the total mean stress transmitted by the coarser fraction, $\alpha_c = 1 - \alpha_f$, is calculated as $p^c/p$, where $p^f$ and $p^c$ are the mean stress transmitted by the finer and the coarser populations, respectively.

Fig. 9 shows the evolution of the $\alpha_f$ during shearing for the dense and loose assemblies of binary mixtures<fnref n="10" />. The $\alpha_f$ for the dense assemblies having $f_c \ge 30\%$ decreased significantly (with more pronounced decrease observed for $35\% \le f_c \le 50\%$ which corresponds with regime 2 in Table 3) from the start of shearing until $\epsilon_a \approx 0.08$ beyond which the change in the $\alpha_f$ became minimal (Fig. 9a). $\alpha_f$ for the loose assemblies generally exhibited minimal changes from the start to the end of shearing (Fig. 9b). In agreement with the data presented by Sufian et al., (2021) for sand-silt mixtures subjected to constant mean stress triaxial compression, the

<footnote n="9">⁹ Note that $\sigma_{ij}^c$ and $\sigma_{ij}^f$ correspond to the contribution of coarse and fine particles to the total stress and not to the mean stress in coarse and fine particles. To compute the stress in the coarse particles, we consider both the coarse-coarse and the coarse-fine contacts. Similarly, to compute the stress in the fine particles, we consider both the fine-fine and the coarse-fine contacts.</footnote>

<footnote n="10">¹⁰ In Appendix B and Fig. A2, we show that the evolution of the deviatoric stress transmitted by the fines is generally similar to the evolution of the mean stress presented in Fig. 8.</footnote>
where the summation runs over all contacts c made by the particle p. $V^p$ is the particle volume, $R_p$ is the radius of the particle, $N_c^p$ is the number of contacts involving the particle, $f_j^c$ is the contact force vector and $n_i$ is the unit branch vector between two particles in contact. The stress tensor at the material point scale (i.e. for a large number of particles), $\sigma_{ij}$, is given as a weighted average of the particle stresses (by particle volume):

$$ \sigma_{ij} = \frac{1}{V} \sum_{p=1}^{N_p} V^p \sigma_{ij}^p \quad (4) $$

where V is the volume of the entire system and $N_p$ is the number of particles in the system. $\sigma_{ij}$ can be decomposed to the contribution of each particle category, i.e. fine and coarse, as:

$$ \sigma_{ij} = \frac{1}{V} \sum_{p=1}^{N_p^c} V^p \sigma_{ij}^p + \frac{1}{V} \sum_{p=1}^{N_p^f} V^p \sigma_{ij}^p \quad (5) $$

where $N_p^c$ and $N_p^f$ are the number of coarse and fine particles in the system, respectively<fnref n="9" />. The mean effective stress, p, for the entire system is calculated as:

$$ p = \frac{\sigma_{xx} + \sigma_{yy} + \sigma_{zz}}{3} \quad (6) $$

The proportion of the total mean stress transmitted by the finer fraction, $\alpha_f$, is calculated as $p^f/p$, and similarly, the proportion of the total mean stress transmitted by the coarser fraction, $\alpha_c = 1 - \alpha_f$, is calculated as $p^c/p$, where $p^f$ and $p^c$ are the mean stress transmitted by the finer and the coarser populations, respectively.

Fig. 9 shows the evolution of the $\alpha_f$ during shearing for the dense and loose assemblies of binary mixtures<fnref n="10" />. The $\alpha_f$ for the dense assemblies having $f_c \ge 30\%$ decreased significantly (with more pronounced decrease observed for $35\% \le f_c \le 50\%$ which corresponds with regime 2 in Table 3) from the start of shearing until $\epsilon_a \approx 0.08$ beyond which the change in the $\alpha_f$ became minimal (Fig. 9a). $\alpha_f$ for the loose assemblies generally exhibited minimal changes from the start to the end of shearing (Fig. 9b). In agreement with the data presented by Sufian et al., (2021) for sand-silt mixtures subjected to constant mean stress triaxial compression, the

<footnote n="9">⁹ Note that $\sigma_{ij}^c$ and $\sigma_{ij}^f$ correspond to the contribution of coarse and fine particles to the total stress and not to the mean stress in coarse and fine particles. To compute the stress in the coarse particles, we consider both the coarse-coarse and the coarse-fine contacts. Similarly, to compute the stress in the fine particles, we consider both the fine-fine and the coarse-fine contacts.</footnote>

<footnote n="10">¹⁰ In Appendix B and Fig. A2, we show that the evolution of the deviatoric stress transmitted by the fines is generally similar to the evolution of the mean stress presented in Fig. 8.</footnote>
observation made here indicates the transmission of stresses from the fine particles to the coarse particles or a concentration of stress on the coarse particles, in the dense assemblies, at the early stage of shearing. This phenomenon occurred such that for each $f_c$, a unique $\alpha_f$ exists for both the initially dense and loose assemblies at large strains (i.e. at the critical state).

<figure><img src="image_9.png" /><figcaption>**Fig. 9** Evolution of the proportion of mean stress transmitted by the fines during shearing for (a) dense assemblies (b) loose assemblies<fnref n="11" /></figcaption></figure>

Fig. 10 shows the contribution of the finer fraction, $\alpha_f$, and the coarser fraction, $\alpha_c$, to the total mean stress for the studied assemblies. At all stages of shearing considered (Fig. 10a-d), the fine particles did not contribute to stress transmission (i.e. $\alpha_f \approx 0$) until the $f_c^{th}$ beyond which there was a progressive increase in the $\alpha_f$ and a consequential reduction in $\alpha_c$. The primary stress-transmitting skeleton at $f_c < f_c^{th}$ is therefore the coarser fraction. Based on a conceptual analysis proposed in the experimental studies on granular mixtures (Thevanayagam et al., 2002; Thevanayagam & Mohan, 2000), the mechanical behaviour of the mixtures should be primarily controlled by the fines at $f_c \ge f_c^{th}$. Vallejo (2001) in his experimental study on binary granular mixtures suggests that at $f_c>60\%$, the fines generally become the primary stress-transmitting matrix while the coarser fraction plays basically no role in stress transmission. The evidence from our micromechanical analysis shows that while the fines began to transmit stress at the $f_c^{th}$, they did not become the primary stress-transmitting fabric until $f_c$ becomes higher than an equivalent fines content $f_c^{eq}$. Here, $f_c^{eq}$ is the estimated $f_c$ at which the total mean

<footnote n="11">¹¹ By comparing the *q/p* data (Figs. 3 & 4) with the $\alpha_f$ data (Fig. 9), it is observed that the fluctuation in the *q/p* data for the assemblies having $f_c \in [30\% ; 40\%]$ is not reflected in the relatively smooth $\alpha_f$ data. This is because the fluctuations evident in the *q/p* data is masked behind the apparently smooth $\alpha_f$ data as a result of the normalisation of the stress transmitted by the fines with the total mean stress ($p^f/p$).</footnote>
observation made here indicates the transmission of stresses from the fine particles to the coarse particles or a concentration of stress on the coarse particles, in the dense assemblies, at the early stage of shearing. This phenomenon occurred such that for each $f_c$, a unique $\alpha_f$ exists for both the initially dense and loose assemblies at large strains (i.e. at the critical state).

<figure><img src="image_10.png" /><figcaption>**Fig. 9** Evolution of the proportion of mean stress transmitted by the fines during shearing for (a) dense assemblies (b) loose assemblies<fnref n="11" /></figcaption></figure>

Fig. 10 shows the contribution of the finer fraction, $\alpha_f$, and the coarser fraction, $\alpha_c$, to the total mean stress for the studied assemblies. At all stages of shearing considered (Fig. 10a-d), the fine particles did not contribute to stress transmission (i.e. $\alpha_f \approx 0$) until the $f_c^{th}$ beyond which there was a progressive increase in the $\alpha_f$ and a consequential reduction in $\alpha_c$. The primary stress-transmitting skeleton at $f_c < f_c^{th}$ is therefore the coarser fraction. Based on a conceptual analysis proposed in the experimental studies on granular mixtures (Thevanayagam et al., 2002; Thevanayagam & Mohan, 2000), the mechanical behaviour of the mixtures should be primarily controlled by the fines at $f_c \ge f_c^{th}$. Vallejo (2001) in his experimental study on binary granular mixtures suggests that at $f_c>60\%$, the fines generally become the primary stress-transmitting matrix while the coarser fraction plays basically no role in stress transmission. The evidence from our micromechanical analysis shows that while the fines began to transmit stress at the $f_c^{th}$, they did not become the primary stress-transmitting fabric until $f_c$ becomes higher than an equivalent fines content $f_c^{eq}$. Here, $f_c^{eq}$ is the estimated $f_c$ at which the total mean

<footnote n="11">¹¹ By comparing the *q/p* data (Figs. 3 & 4) with the $\alpha_f$ data (Fig. 9), it is observed that the fluctuation in the *q/p* data for the assemblies having $f_c \in [30\% ; 40\%]$ is not reflected in the relatively smooth $\alpha_f$ data. This is because the fluctuations evident in the *q/p* data is masked behind the apparently smooth $\alpha_f$ data as a result of the normalisation of the stress transmitted by the fines with the total mean stress ($p^f/p$).</footnote>
387

pressure is equally shared between the finer and the coarser fraction. For the studied binary mixtures, $f_c^{eq}$ is estimated to range from 48% to 53% for the dense assemblies at the initial state and the peak; and reaches 55% for the loose assemblies at the initial state and for all assemblies at the critical state. Contrary to Vallejo (2001), while the fines transmit the largest fraction of the total mean stress for $f_c > f_c^{eq}$, the coarse particles continue to make a secondary and yet significant contribution to the total mean stress for much larger fine contents (at least up to $f_c$=70% studied here).

388

389

390

391

392

393

394

Since particles stresses result from contact forces, we estimate in Fig. 11, the proportion of the total number contacts shared by the fines<fnref n="12" />, $C_n^f$, and the coarse particles, $C_n^c$. Comparing these geometric statistics to the stress contributions enables a better understanding of the trends observed in Fig. 10. Again, the trends in Fig. 11 can be classified into two regimes delineated by the threshold fines content, $f_c^{th}$. For $f_c < f_c^{th}$, $C_n^f$ generally increased with $f_c$ but remained lower than $C_n^c$, indicating that in this regime, the coarse particles are involved in a larger proportion of the total number of contacts than the fines. For $f_c \ge f_c^{th}$, $C_n^f$ increased monotonically with $f_c$ and was higher than $C_n^c$; in this case, the fines became dominant in terms of their proportion of the total number of contacts. By comparing Fig. 10 and 11, we observe that although the fines became dominant in terms of their proportion of the total number of contacts from $f_c = f_c^{th}$, they did not primarily contribute to the total mean stress until $f_c > f_c^{eq}$. This indicates that most of the contacts involving fine grains do not carry large forces for $f_c^{th} < f_c < f_c^{eq}$.

395

396

397

398

399

400

401

402

403

404

405

406

<footnote n="12">¹² The contacts made by the coarse particles include both the coarse-coarse and the coarse-fine contacts while the contacts made by the fine particles include both the fine-fine and the coarse-fine contacts. The proportion of the total number of contacts shared by the fines, $C_n^f$, is calculated as: $C_n^f = \frac{2N_{f-f}+N_{c-f}}{2N_{f-f}+2N_{c-f}+2N_{c-c}}$.</footnote>
387

pressure is equally shared between the finer and the coarser fraction. For the studied binary mixtures, $f_c^{eq}$ is estimated to range from 48% to 53% for the dense assemblies at the initial state and the peak; and reaches 55% for the loose assemblies at the initial state and for all assemblies at the critical state. Contrary to Vallejo (2001), while the fines transmit the largest fraction of the total mean stress for $f_c > f_c^{eq}$, the coarse particles continue to make a secondary and yet significant contribution to the total mean stress for much larger fine contents (at least up to $f_c$=70% studied here).

388

389

390

391

392

393

394

Since particles stresses result from contact forces, we estimate in Fig. 11, the proportion of the total number contacts shared by the fines<fnref n="12" />, $C_n^f$, and the coarse particles, $C_n^c$. Comparing these geometric statistics to the stress contributions enables a better understanding of the trends observed in Fig. 10. Again, the trends in Fig. 11 can be classified into two regimes delineated by the threshold fines content, $f_c^{th}$. For $f_c < f_c^{th}$, $C_n^f$ generally increased with $f_c$ but remained lower than $C_n^c$, indicating that in this regime, the coarse particles are involved in a larger proportion of the total number of contacts than the fines. For $f_c \ge f_c^{th}$, $C_n^f$ increased monotonically with $f_c$ and was higher than $C_n^c$; in this case, the fines became dominant in terms of their proportion of the total number of contacts. By comparing Fig. 10 and 11, we observe that although the fines became dominant in terms of their proportion of the total number of contacts from $f_c = f_c^{th}$, they did not primarily contribute to the total mean stress until $f_c > f_c^{eq}$. This indicates that most of the contacts involving fine grains do not carry large forces for $f_c^{th} < f_c < f_c^{eq}$.

395

396

397

398

399

400

401

402

403

404

405

406

<footnote n="12">¹² The contacts made by the coarse particles include both the coarse-coarse and the coarse-fine contacts while the contacts made by the fine particles include both the fine-fine and the coarse-fine contacts. The proportion of the total number of contacts shared by the fines, $C_n^f$, is calculated as: $C_n^f = \frac{2N_{f-f}+N_{c-f}}{2N_{f-f}+2N_{c-f}+2N_{c-c}}$.</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 10 Proportion of mean stress transmitted by the fines and coarse particles at (a) the initial state for dense assemblies (b) the initial state for the loose assemblies (c) the peak (d) the critical state (unique trends for both dense and loose assemblies)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 11 Proportion of the total contacts belonging to the finer and coarser fraction at (a) the initial state for dense assemblies (b) the initial state for the loose assemblies (c) the peak (d) the critical state (unique trends for both dense and loose assemblies)</figcaption></figure>

407

408

409

410

411

412

413

414

415

416
<figure><img src="image_11.png" /><figcaption>Fig. 10 Proportion of mean stress transmitted by the fines and coarse particles at (a) the initial state for dense assemblies (b) the initial state for the loose assemblies (c) the peak (d) the critical state (unique trends for both dense and loose assemblies)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11 Proportion of the total contacts belonging to the finer and coarser fraction at (a) the initial state for dense assemblies (b) the initial state for the loose assemblies (c) the peak (d) the critical state (unique trends for both dense and loose assemblies)</figcaption></figure>

407

408

409

410

411

412

413

414

415

416
Figs. 12 and 13 show the cumulative distributions of the individual stresses sustained by the particles normalised by the average particle stresses in the whole assembly ($P^p/\overline{P^p}$), for the coarser and finer fraction respectively, at the initial states, the peak and the critical state. This is aimed at understanding the range and state of particle stresses existing within the particle size fractions in each assembly at different shearing stages. A wider range of stresses exist within the finer fraction (Fig. 13) than the coarser fraction (Fig. 12). For the coarser fraction, the proportion of the particles with stress values below the average stress (i.e. $P^p/\overline{P^p} \le 1$) was similar for all $f_c$ values considered. Below the average stress, the cumulative distribution shifts to the right as $f_c$ increased; the reverse was observed above the average stress. The width of the distributions for the coarser fraction tend to become smaller as $f_c$ increased (Fig. 12a-d). This indicates a more even stress distribution within the coarser particles as $f_c$ increased. The proportion of the fines with stress values below the average (i.e. $P^p/\overline{P^p} \le 1$) decreased as $f_c$ increased (Fig. 13).

The proportion of the fines bearing no stress (i.e. $P^p/\overline{P^p} = 0$), referred to as the non-active fines, $R_f$, are illustrated by the vertical line at $P^p/\overline{P^p} = 0$ in Fig. 13 and are plotted in Fig. 14a. Fig. 14b shows the change in the percentage of the non-active fines present within the assemblies during shearing, $\Delta R_f$ (i.e. the difference in the $R_f$ at the initial states and the critical state). Within the assemblies having $f_c<f_c^{th}$, virtually all the fines (98.3% to 99.9%) are non-active. This confirms that at $f_c<f_c^{th}$, the finer fraction does not contribute to stress transfer. For $f_c^{th} < f_c < 45\%$, $R_f$ decreases more significantly than for $f_c > 45\%$ (Fig. 14a). The drop in the proportion of non-active fine (within $f_c^{th} < f_c < 45\%$) depends significantly on the relative density of the assemblies (dense or loose) and the stage of shearing (initial, peak or critical state). Three regimes are identified following the analysis here. For $f_c<f_c^{th}$, the fines are largely dormant whatever the sample state<fnref n="13" />. Also, $\Delta R_f \approx 0$ (Fig. 14b), indicating shearing had no effect on the activity of the fines. For $f_c^{th} \le f_c < 45\%$, the proportion of active fines is very sensitive to the packing density and the stage of shearing (Fig. 14a). Also, $\Delta R_f$ is highest (Fig. 14b), suggesting that more of the active fines initially present within the assemblies became non-active by the end of shearing (38%-85% and 11-12% for dense and loose assemblies, respectively). For $f_c \ge 45\%$, the proportion of non-active fines is rather low and exhibit a limited sensitivity to the packing density and the stage of shearing.

<footnote n="13">¹³ Note that, depending on the sample preparation procedure, we could imagine having fewer non-active fines at initial state. This is for instance the case when samples are prepared in the lab using the moist tamping method that gather fines around contacts between coarse grains because of capillary effects.</footnote>
Figs. 12 and 13 show the cumulative distributions of the individual stresses sustained by the particles normalised by the average particle stresses in the whole assembly ($P^p/\overline{P^p}$), for the coarser and finer fraction respectively, at the initial states, the peak and the critical state. This is aimed at understanding the range and state of particle stresses existing within the particle size fractions in each assembly at different shearing stages. A wider range of stresses exist within the finer fraction (Fig. 13) than the coarser fraction (Fig. 12). For the coarser fraction, the proportion of the particles with stress values below the average stress (i.e. $P^p/\overline{P^p} \le 1$) was similar for all $f_c$ values considered. Below the average stress, the cumulative distribution shifts to the right as $f_c$ increased; the reverse was observed above the average stress. The width of the distributions for the coarser fraction tend to become smaller as $f_c$ increased (Fig. 12a-d). This indicates a more even stress distribution within the coarser particles as $f_c$ increased. The proportion of the fines with stress values below the average (i.e. $P^p/\overline{P^p} \le 1$) decreased as $f_c$ increased (Fig. 13).

The proportion of the fines bearing no stress (i.e. $P^p/\overline{P^p} = 0$), referred to as the non-active fines, $R_f$, are illustrated by the vertical line at $P^p/\overline{P^p} = 0$ in Fig. 13 and are plotted in Fig. 14a. Fig. 14b shows the change in the percentage of the non-active fines present within the assemblies during shearing, $\Delta R_f$ (i.e. the difference in the $R_f$ at the initial states and the critical state). Within the assemblies having $f_c<f_c^{th}$, virtually all the fines (98.3% to 99.9%) are non-active. This confirms that at $f_c<f_c^{th}$, the finer fraction does not contribute to stress transfer. For $f_c^{th} < f_c < 45\%$, $R_f$ decreases more significantly than for $f_c > 45\%$ (Fig. 14a). The drop in the proportion of non-active fine (within $f_c^{th} < f_c < 45\%$) depends significantly on the relative density of the assemblies (dense or loose) and the stage of shearing (initial, peak or critical state). Three regimes are identified following the analysis here. For $f_c<f_c^{th}$, the fines are largely dormant whatever the sample state<fnref n="13" />. Also, $\Delta R_f \approx 0$ (Fig. 14b), indicating shearing had no effect on the activity of the fines. For $f_c^{th} \le f_c < 45\%$, the proportion of active fines is very sensitive to the packing density and the stage of shearing (Fig. 14a). Also, $\Delta R_f$ is highest (Fig. 14b), suggesting that more of the active fines initially present within the assemblies became non-active by the end of shearing (38%-85% and 11-12% for dense and loose assemblies, respectively). For $f_c \ge 45\%$, the proportion of non-active fines is rather low and exhibit a limited sensitivity to the packing density and the stage of shearing.

<footnote n="13">¹³ Note that, depending on the sample preparation procedure, we could imagine having fewer non-active fines at initial state. This is for instance the case when samples are prepared in the lab using the moist tamping method that gather fines around contacts between coarse grains because of capillary effects.</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 12 Cumulative distribution of normalised particle mean stress for coarser fraction (a) dense; initial state (b) loose; initial state (c) dense; peak state (d) critical state for both dense (solid line) and loose (dashed line) assemblies</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 13 Cumulative distribution of normalised particle mean stress for finer fraction (a) dense; initial state (b) loose; initial state (c) dense; peak state (d) critical state for both dense (solid line) and loose (dashed line) assemblies</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12 Cumulative distribution of normalised particle mean stress for coarser fraction (a) dense; initial state (b) loose; initial state (c) dense; peak state (d) critical state for both dense (solid line) and loose (dashed line) assemblies</figcaption></figure>

447

448

449

450

451

<figure><img src="image_14.png" /><figcaption>Fig. 13 Cumulative distribution of normalised particle mean stress for finer fraction (a) dense; initial state (b) loose; initial state (c) dense; peak state (d) critical state for both dense (solid line) and loose (dashed line) assemblies</figcaption></figure>

452

453

454

455

456
<figure><img src="image_14.png" /><figcaption>Fig. 14 (a) Percentage of non-active fines, $R_f$, at different $f_c$ (b) Change in the percentage of non-active fines during shearing (filled circle for dense assemblies and open circles for loose assemblies). $\Delta R_f>0$ and $\Delta R_f<0$ mean higher and lower percentages of inactive fines at the end of shearing in comparison to the initial state, respectively.</figcaption></figure>

5. Contact based stress-transmission

Using the Love-Weber contact based stress tensor including Bagi boundary term (Bagi, 1999), the contributions of all contact types to the total mean stress transmitted by an entire granular system is given as:

$$ \sigma_{ij}^{c} = \frac{1}{V} \sum_{c=1}^{N_c} f_{j}^{c} l_{i} \qquad (11) $$

where $N_c$ is the total number of contacts in the granular system, $V$ is the volume of the system, $l$ is the branch vector ($\|l\| = R_1 + R_2$ for a sphere 1 to sphere 2 contact; $\|l\| = R_p$ for a sphere $p$ to wall contact), and $f^c$ is the contact force vector. This stress tensor can be decomposed based on the contact type as:

$$ \sigma_{ij}^{k} = \frac{1}{V} \sum_{c=1}^{N_c^k} f_{j}^{c} l_{i} \qquad (12) $$

where $k \in \{c-c, c-f, f-f, s-wl\}$. c-c, c-f, f-f, and s-wl denote coarse to coarse, coarse to fine, fine to fine and sphere to wall contacts, respectively. The proportion of the total mean stress from each contact type is given by $\alpha_k = p^k/p$. Note that the boundary term is expected to vanish as soon as the sample domain is sufficiently large (which should be the case when the REV condition is met).
<figure><img src="image_15.png" /><figcaption>Fig. 14 (a) Percentage of non-active fines, $R_f$, at different $f_c$ (b) Change in the percentage of non-active fines during shearing (filled circle for dense assemblies and open circles for loose assemblies). $\Delta R_f>0$ and $\Delta R_f<0$ mean higher and lower percentages of inactive fines at the end of shearing in comparison to the initial state, respectively.</figcaption></figure>

5. Contact based stress-transmission

Using the Love-Weber contact based stress tensor including Bagi boundary term (Bagi, 1999), the contributions of all contact types to the total mean stress transmitted by an entire granular system is given as:

$$ \sigma_{ij}^{c} = \frac{1}{V} \sum_{c=1}^{N_c} f_{j}^{c} l_{i} \qquad (11) $$

where $N_c$ is the total number of contacts in the granular system, $V$ is the volume of the system, $l$ is the branch vector ($\|l\| = R_1 + R_2$ for a sphere 1 to sphere 2 contact; $\|l\| = R_p$ for a sphere $p$ to wall contact), and $f^c$ is the contact force vector. This stress tensor can be decomposed based on the contact type as:

$$ \sigma_{ij}^{k} = \frac{1}{V} \sum_{c=1}^{N_c^k} f_{j}^{c} l_{i} \qquad (12) $$

where $k \in \{c-c, c-f, f-f, s-wl\}$. c-c, c-f, f-f, and s-wl denote coarse to coarse, coarse to fine, fine to fine and sphere to wall contacts, respectively. The proportion of the total mean stress from each contact type is given by $\alpha_k = p^k/p$. Note that the boundary term is expected to vanish as soon as the sample domain is sufficiently large (which should be the case when the REV condition is met).
Fig. 15 shows the evolution of the proportion of the stress transmitted by each contact type during shearing, for the dense and loose assemblies. At large strains (i.e., $\epsilon_a > 0.3$), the contribution of the coarse to coarse contacts to the total mean stress, $\alpha_{c-c}$, decreases with $f_c$ such that $\alpha_{c-c} \approx 0$ for $f_c > 40\%$ (Fig. 15a&d). The contribution of the fine-fine contacts to the total mean stress, $\alpha_{f-f}$, increases with $f_c$ for $f_c > f_c^{th}$ while $\alpha_{f-f} \approx 0$ for $f_c < f_c^{th}$ (Fig. 15c&f). The contribution of the coarse to fine contacts, $\alpha_{c-f}$, initially increases with $f_c$ until $f_c = 40\%$ and thereafter decreases with further increase in $f_c$ (Fig. 15b&e). During the early stage of shearing (i.e., $\epsilon_a < 0.1$), we observed a decrease in $\alpha_{c-c}$ and $\alpha_{f-f}$ which is accompanied by a commensurate increase in the $\alpha_{c-f}$. This indicates a redistribution of stress from the c-c and f-f contacts to the c-f contacts. Although the coarse particles lost stress as indicated in the decrease in $\alpha_{c-c}$, they are also gained stress from the fines as indicated in both the increase in $\alpha_{c-f}$ and the decrease in $\alpha_{f-f}$. Ultimately, there was a redistribution of stress from the fines to the coarse particles as indicated in Fig. 9a.

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

<caption>**Fig. 15** Evolution of the proportion of mean stress transmitted by contact types during shearing for (a) c-c contacts; dense (b) c-f contacts; dense (c) f-f contacts; dense (d) c-c contacts; loose (e) c-f contacts; loose (f) f-f contacts; loose</caption>
Fig. 15 shows the evolution of the proportion of the stress transmitted by each contact type during shearing, for the dense and loose assemblies. At large strains (i.e., $\epsilon_a > 0.3$), the contribution of the coarse to coarse contacts to the total mean stress, $\alpha_{c-c}$, decreases with $f_c$ such that $\alpha_{c-c} \approx 0$ for $f_c > 40\%$ (Fig. 15a&d). The contribution of the fine-fine contacts to the total mean stress, $\alpha_{f-f}$, increases with $f_c$ for $f_c > f_c^{th}$ while $\alpha_{f-f} \approx 0$ for $f_c < f_c^{th}$ (Fig. 15c&f). The contribution of the coarse to fine contacts, $\alpha_{c-f}$, initially increases with $f_c$ until $f_c = 40\%$ and thereafter decreases with further increase in $f_c$ (Fig. 15b&e). During the early stage of shearing (i.e., $\epsilon_a < 0.1$), we observed a decrease in $\alpha_{c-c}$ and $\alpha_{f-f}$ which is accompanied by a commensurate increase in the $\alpha_{c-f}$. This indicates a redistribution of stress from the c-c and f-f contacts to the c-f contacts. Although the coarse particles lost stress as indicated in the decrease in $\alpha_{c-c}$, they are also gained stress from the fines as indicated in both the increase in $\alpha_{c-f}$ and the decrease in $\alpha_{f-f}$. Ultimately, there was a redistribution of stress from the fines to the coarse particles as indicated in Fig. 9a.

<figure><img src="image_16.png" /><figcaption>**Fig. 15** Evolution of the proportion of mean stress transmitted by contact types during shearing for (a) c-c contacts; dense (b) c-f contacts; dense (c) f-f contacts; dense (d) c-c contacts; loose (e) c-f contacts; loose (f) f-f contacts; loose</figcaption></figure>
500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

Fig. 16 shows the proportion of the mean stress transmitted by each contact types at the initial states, the peak and the critical state. The trends observed between the contact types<fnref n="14" /> can be classified into four regimes based on the $f_c$. In the first regime ($f_c<f_c^{th}$), the contribution of the fine to fine contacts to the total mean stress, $\alpha_{f-f} \approx 0$, indicating that the f-f contacts do not contribute to mean stress below the $f_c^{th}$. In this regime, the condition $\alpha_{c-c} > \alpha_{c-f} > \alpha_{f-f}$ generally holds, hence, the regime is referred to as coarse-dominated, following Vallejo (2001) and Sufian et al. (2021). Within the second regime referred to as transitional coarse-dominated ($f_c^{th} < f_c < f_c^{tr}$, where $f_c^{tr}$ is referred here as transitional fines content where $\alpha_{c-c} = \alpha_{f-f}$), the condition $\alpha_{c-f} > \alpha_{c-c} > \alpha_{f-f}$ generally holds. In the third regime where $f_c^{tr} < f_c < f_c^{eq}$ ($f_c^{eq}$ was defined in Section 4 as the $f_c$ at which the total mean stress is equally shared between the coarser and the finer fractions) and referred to as transitional fines-dominated, the condition $\alpha_{c-f} > \alpha_{f-f} > \alpha_{c-c}$ holds. The transitional fines content, $f_c^{tr}$, delineating the transitional regimes (i.e. the transitional-coarse dominated and the transitional-fines dominated regimes) was found around $f_c = 40\%$ for the assemblies and the shearing stages considered here (Fig. 16a-d). The transitional regime based on the micromechanical analysis ranged from $f_c^{th}$ to $f_c^{eq}$ and is therefore within $f_c = 30\%-48\%$ and $f_c = 35\%-55\%$ for the dense and loose assemblies, respectively. For $f_c > f_c^{eq}$ (i.e., in the fourth regime), the condition $\alpha_{f-f} > \alpha_{c-f} > \alpha_{c-c}$ holds, hence, the regime is fines-dominated. While the classification in the experimental study by Vallejo (2001) was based on void ratio, we present here a similar classification based on the micromechanical analysis conducted, with distinct delineations in agreement with the 3-D DEM study on binary mixtures having size ratio, $\lambda=4$, reported by Minh et al., (2014).

<footnote n="14">¹⁴ Excluding the contributions to the mean stress from the contacts between the spherical particles and the walls, $\alpha_{s-wl}$, which were generally less than 0.15. We observed that $\alpha_{s-wl}$ values for the underfilled assemblies ($f_c<f_c^{th}$) were generally higher than the overfilled assemblies. This suggests that a greater number of particles is required to ensure a representative element volume (REV) in underfilled assemblies in comparison to the overfilled assemblies.</footnote>
Fig. 16 shows the proportion of the mean stress transmitted by each contact types at the initial states, the peak and the critical state. The trends observed between the contact types<fnref n="14" /> can be classified into four regimes based on the $f_c$. In the first regime ($f_c<f_c^{th}$), the contribution of the fine to fine contacts to the total mean stress, $\alpha_{f-f} \approx 0$, indicating that the f-f contacts do not contribute to mean stress below the $f_c^{th}$. In this regime, the condition $\alpha_{c-c} > \alpha_{c-f} > \alpha_{f-f}$ generally holds, hence, the regime is referred to as coarse-dominated, following Vallejo (2001) and Sufian et al. (2021). Within the second regime referred to as transitional coarse-dominated ($f_c^{th} < f_c < f_c^{tr}$, where $f_c^{tr}$ is referred here as transitional fines content where $\alpha_{c-c} = \alpha_{f-f}$), the condition $\alpha_{c-f} > \alpha_{c-c} > \alpha_{f-f}$ generally holds. In the third regime where $f_c^{tr} < f_c < f_c^{eq}$ ($f_c^{eq}$ was defined in Section 4 as the $f_c$ at which the total mean stress is equally shared between the coarser and the finer fractions) and referred to as transitional fines-dominated, the condition $\alpha_{c-f} > \alpha_{f-f} > \alpha_{c-c}$ holds. The transitional fines content, $f_c^{tr}$, delineating the transitional regimes (i.e. the transitional-coarse dominated and the transitional-fines dominated regimes) was found around $f_c = 40\%$ for the assemblies and the shearing stages considered here (Fig. 16a-d). The transitional regime based on the micromechanical analysis ranged from $f_c^{th}$ to $f_c^{eq}$ and is therefore within $f_c = 30\%-48\%$ and $f_c = 35\%-55\%$ for the dense and loose assemblies, respectively. For $f_c>f_c^{eq}$ (i.e., in the fourth regime), the condition $\alpha_{f-f} > \alpha_{c-f} > \alpha_{c-c}$ holds, hence, the regime is fines-dominated. While the classification in the experimental study by Vallejo (2001) was based on void ratio, we present here a similar classification based on the micromechanical analysis conducted, with distinct delineations in agreement with the 3-D DEM study on binary mixtures having size ratio, $\lambda=4$, reported by Minh et al., (2014).

<footnote n="14">¹⁴ Excluding the contributions to the mean stress from the contacts between the spherical particles and the walls, $\alpha_{s-wl}$, which were generally less than 0.15. We observed that $\alpha_{s-wl}$ values for the underfilled assemblies ($f_c<f_c^{th}$) were generally higher than the overfilled assemblies. This suggests that a greater number of particles is required to ensure a representative element volume (REV) in underfilled assemblies in comparison to the overfilled assemblies.</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 16** Proportion of mean stress transmitted by contact type at (a) the initial state for dense assemblies b) the initial state for the loose assemblies (c) at the peak (d) the critical state (unique trends for both dense and loose assemblies). Void ratio data represented with dashed grey lines are plotted on the right side of each subplot to help identify the different regimes.</figcaption></figure>

In a similar fashion to Fig. 11, Fig. 17 shows the proportion of all contacts belonging to each contact types at the different stages of shearing. The trends between the contact types here could be generally grouped into two categories delineated by the threshold fines content<fnref n="15" />.

- For $f_c<f_c^{th}$, the proportions of the total contacts belonging to the f-f contacts, $C_{f-f}^t$, were lower than those belonging to the c-c contacts, $C_{c-c}^t$. The condition, $C_{c-c}^t > C_{c-f}^t > C_{f-f}^t$ generally holds, except for $f_c=30\%$ in the loose samples (Fig. 17b) and at the critical state (Fig. 17d) where the condition $C_{c-f}^t > C_{c-c}^t > C_{f-f}^t$ holds.

- For $f_c \ge f_c^{th}$, the f-f contacts dominate the assemblies, therefore the condition $C_{f-f}^t > C_{c-f}^t > C_{c-c}^t$ holds.

<footnote n="15">¹⁵ The proportions of all contact types belonging to the wall, $C_{wl}^t$, were less than 0.17 and 0.08 for $f_c<f_c^{th}$ and for $f_c>f_c^{th}$, respectively.</footnote>
<figure><img src="image_17.png" /><figcaption>**Fig. 16** Proportion of mean stress transmitted by contact type at (a) the initial state for dense assemblies b) the initial state for the loose assemblies (c) at the peak (d) the critical state (unique trends for both dense and loose assemblies). Void ratio data represented with dashed grey lines are plotted on the right side of each subplot to help identify the different regimes.</figcaption></figure>

In a similar fashion to Fig. 11, Fig. 17 shows the proportion of all contacts belonging to each contact types at the different stages of shearing. The trends between the contact types here could be generally grouped into two categories delineated by the threshold fines content<fnref n="15" />.

- For $f_c<f_c^{th}$, the proportions of the total contacts belonging to the f-f contacts, $C_{f-f}^t$, were lower than those belonging to the c-c contacts, $C_{c-c}^t$. The condition, $C_{c-c}^t > C_{c-f}^t > C_{f-f}^t$ generally holds, except for $f_c=30\%$ in the loose samples (Fig. 17b) and at the critical state (Fig. 17d) where the condition $C_{c-f}^t > C_{c-c}^t > C_{f-f}^t$ holds.

- For $f_c \ge f_c^{th}$, the f-f contacts dominate the assemblies, therefore the condition $C_{f-f}^t > C_{c-f}^t > C_{c-c}^t$ holds.

<footnote n="15">¹⁵ The proportions of all contact types belonging to the wall, $C_{wl}^t$, were less than 0.17 and 0.08 for $f_c<f_c^{th}$ and for $f_c>f_c^{th}$, respectively.</footnote>
From the threshold fines content, the $C_{c-c}^t$ values were significantly low (ranging from 2.1e-2 to 9e-5) at all states, and was zero for $f_c=70\%$ at the critical state (Fig. 17d). This shows that the fines disperse the coarse particles (i.e. $C_{c-c}^t$ becomes minimal) from the threshold fines content, $f_c^{th}$, which is in agreement with the 3D DEM study on binary mixtures by Minh et al., (2014), and not beyond $f_c^{th}$ as suggested in the experimental study on sandy-gravels by Skempton & Brogan (1994). The suggestion by Salgado et al., (2000) in their experimental study, that the fines controls the mechanical behaviour of sand-silt mixtures when the sand particles are completely floating ($C_{c-c}^t \approx 0$) in the silts is not supported by the data shown here. While the fines disperse the coarse particles from $f_c^{th}$, they did not become the primary stress-transmitting matrix until $f_c^{eq}$ is reached as explained in Section 4. In the transitional regime, there is the possibility that some fines get trapped between the coarse particles while being loosely connected to the rest of the fines in the assemblies.

Also, while for $f_c<f_c^{th}$, the $C_{c-f}^t \approx 0$ in the dense assemblies (Fig. 17a), in the loose assemblies (Fig. 17b), $C_{c-f}^t > 0$. This suggests that, in agreement with Thevanayagam et al., (2002) and Shire, (2014), the fines in the dense assemblies are confined within the voids between the coarse particles with little interaction with the coarse matrix; in the loose assemblies, the fines interact more with the coarse particles. The fact that the fines are interacting more with the coarse particles in the loose case may indicate the fines are trapped within the coarse particles thereby creating larger voids which may require more fines to fill in comparison to the dense assemblies. This observation is consistent with Salgado et al., (2000) and Lade & Yamamuro (1997) who found more fines on the surfaces of the coarse particles in their loose assemblies, in comparison to the dense assemblies, and suggested that this is responsible for the larger drop in the $e_{min}$ than in the $e_{max}$ for a given increase in the fines content. The drop in the $e_{max}$ is mitigated than in the $e_{min}$ because the fines interacting with the coarse particles helps to create larger voids hence the lower drop in the void ratio in the loose case.
From the threshold fines content, the $C_{c-c}^t$ values were significantly low (ranging from 2.1e-2 to 9e-5) at all states, and was zero for $f_c=70\%$ at the critical state (Fig. 17d). This shows that the fines disperse the coarse particles (i.e. $C_{c-c}^t$ becomes minimal) from the threshold fines content, $f_c^{th}$, which is in agreement with the 3D DEM study on binary mixtures by Minh et al., (2014), and not beyond $f_c^{th}$ as suggested in the experimental study on sandy-gravels by Skempton & Brogan (1994). The suggestion by Salgado et al., (2000) in their experimental study, that the fines controls the mechanical behaviour of sand-silt mixtures when the sand particles are completely floating ($C_{c-c}^t \approx 0$) in the silts is not supported by the data shown here. While the fines disperse the coarse particles from $f_c^{th}$, they did not become the primary stress-transmitting matrix until $f_c^{eq}$ is reached as explained in Section 4. In the transitional regime, there is the possibility that some fines get trapped between the coarse particles while being loosely connected to the rest of the fines in the assemblies.

Also, while for $f_c<f_c^{th}$, the $C_{c-f}^t \approx 0$ in the dense assemblies (Fig. 17a), in the loose assemblies (Fig. 17b), $C_{c-f}^t > 0$. This suggests that, in agreement with Thevanayagam et al., (2002) and Shire, (2014), the fines in the dense assemblies are confined within the voids between the coarse particles with little interaction with the coarse matrix; in the loose assemblies, the fines interact more with the coarse particles. The fact that the fines are interacting more with the coarse particles in the loose case may indicate the fines are trapped within the coarse particles thereby creating larger voids which may require more fines to fill in comparison to the dense assemblies. This observation is consistent with Salgado et al., (2000) and Lade & Yamamuro (1997) who found more fines on the surfaces of the coarse particles in their loose assemblies, in comparison to the dense assemblies, and suggested that this is responsible for the larger drop in the $e_{min}$ than in the $e_{max}$ for a given increase in the fines content. The drop in the $e_{max}$ is mitigated than in the $e_{min}$ because the fines interacting with the coarse particles helps to create larger voids hence the lower drop in the void ratio in the loose case.
<figure><img src="image_17.png" /></figure>

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

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

<figure><img src="image_20.png" /><figcaption>**Fig. 17** Proportion of the total contacts belonging to each category of contact type at (a) the initial state for dense assemblies (b) the initial state for the loose assemblies (c) the peak (d) the critical state (unique trends for both dense and loose assemblies). Void ratio data represented with dashed grey lines are plotted on the right side of each subplot to help identify the different regimes.</figcaption></figure>

**6. Stress-based skeleton void ratio**

The idea of sand skeleton void ratio was introduced in the experimental study by Kuerbis (1989) in order to understand the behaviour of silty sands with various fines content under undrained triaxial shearing. This proposition, which was based on a conceptual analysis of the fabric of soil mixtures, have been applied in other experimental and numerical studies (Benahmed et al., 2015; Chang & Yin, 2011; Lade & Yamamuro, 1997; Ni et al., 2004; Pitman et al., 1994; Rahman et al., 2008; Thevanayagam, 1998; Thevanayagam et al., 2002; Vaid, 1994). In order to understand the critical state behaviour of the gap-graded assemblies subjected to drained shearing in this study, we propose a skeleton void ratio based on the stress sustained by individual particles. The skeleton void ratio, *e*<sub>skel</sub>, is here defined as the void ratio of a granular system when the fines with particle stress, *P*<sup>*p*</sup>, lower than a threshold stress value are regarded as part of the void space and are therefore excluded. Note that we propose to exclude only the fines (and not the coarse particles) for our skeleton void ratio calculation. This is based on the fact that loosely stressed fine particles may be eroded through the constrictions
<figure><img src="image_18.png" /></figure>

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

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

<figure><img src="image_21.png" /><figcaption>**Fig. 17** Proportion of the total contacts belonging to each category of contact type at (a) the initial state for dense assemblies (b) the initial state for the loose assemblies (c) the peak (d) the critical state (unique trends for both dense and loose assemblies). Void ratio data represented with dashed grey lines are plotted on the right side of each subplot to help identify the different regimes.</figcaption></figure>

## 6. Stress-based skeleton void ratio

The idea of sand skeleton void ratio was introduced in the experimental study by Kuerbis (1989) in order to understand the behaviour of silty sands with various fines content under undrained triaxial shearing. This proposition, which was based on a conceptual analysis of the fabric of soil mixtures, have been applied in other experimental and numerical studies (Benahmed et al., 2015; Chang & Yin, 2011; Lade & Yamamuro, 1997; Ni et al., 2004; Pitman et al., 1994; Rahman et al., 2008; Thevanayagam, 1998; Thevanayagam et al., 2002; Vaid, 1994). In order to understand the critical state behaviour of the gap-graded assemblies subjected to drained shearing in this study, we propose a skeleton void ratio based on the stress sustained by individual particles. The skeleton void ratio, $e_{skel}$, is here defined as the void ratio of a granular system when the fines with particle stress, $P^p$, lower than a threshold stress value are regarded as part of the void space and are therefore excluded. Note that we propose to exclude only the fines (and not the coarse particles) for our skeleton void ratio calculation. This is based on the fact that loosely stressed fine particles may be eroded through the constrictions
of the soil matrix whereas coarse particles cannot as they are usually larger than the constriction size (Garner & Fannin, 2010; ICOLD, 2017; A. W. Skempton & Brogan, 1994; Thevanayagam & Mohan, 2000). In Figure 17 we show the variation of the proportion of the total mean stress transmitted by the “skeleton” finer fraction, $\alpha_f^{skel}$, for different stress thresholds, where the threshold stress is defined as:

$$ \text{Threshold stress} = \underset{\substack{\text{threshold} \\ \text{coefficient}}}{x} \times \underset{\substack{\text{Average} \\ \text{particle stress}}}{\overline{P^p}} \qquad (10) $$

In Fig. 18, the threshold coefficient, $x$, ranges from 0 to 2.

- Unsurprisingly, the contribution of the fine remains negligible whatever the threshold value for $f_c<f_c^{th}$ (the contribution for all the fines, $\alpha_f$, is already negligible).

- For $f_c^{th} \le f_c \le 45\%$, generally no change is observed in $\alpha_f^{skel}$ until $x = 0.1$, which means that the fine grains transmitting a pressure $P^p$ lower than 10% of the mean pressure $\overline{P^p}$, contribute only marginally to the stress transmitted by the finer fraction.

- For $f_c>45\%$, no change was observed in $\alpha_f^{skel}$ until $x > 0.25$.

From the above analysis, we can retain $x=0.1$ as a reasonable threshold to define the skeleton fines for all fines content.

Fig. 19 shows the proportion of the fines left in an assembly after each exclusion. These “skeleton” fines are included in the solid volume while computing the stress-based skeleton void ratio, $e_{skel}$ plotted in Fig. 20. For $f_c \ge f_c^{th}$, the proportion of the fines included as part of the assembly skeleton decreases monotonically as the threshold stress increases (Fig. 19a-d). The mechanical void ratio, $e_{mech}$ has been computed in the literature as an alternative void ratio index to the void ratio, $e$, of granular materials (Liu et al., 2022; Otsubo, 2016). $e_{mech}$ excludes the particles with $c \le 1$. In Fig. 20, we compare the evolutions of $e$, $e_{mech}$ and $e_{skel}$ (for $x = 0.1$), $e_{skel}$; $x = 0.1$, for the different fines contents. These evolutions are put in parallel with the peak and the critical state stress ratios.

As shown in Figure 19, $e_{skel}$; $x = 0.1$, shows a negative correlation with the critical state strength, $(q/p)_{crit}$ (Fig. 20b & d), and no correlation with the peak strength, $(q/p)_{max}$ (Fig. 20a). Instead, at the peak, $(q/p)_{max}$ correlates negatively with $e$, (Fig. 20a & c), in agreement with established relationship between the initial void ratio and the peak strength in the literature (Adesina et al., 2023; Holubec & D’Appolonia, 1973; Ng, 2004). In order to determine the best performing void ratio index among the indexes considered here, we estimated the goodness of the fit, $R^2$, for the relationship between the void ratio values determined using the indexes and the strength exhibited by the assemblies (Table A.1 in the Appendix C). It is obvious that
of the soil matrix whereas coarse particles cannot as they are usually larger than the constriction size (Garner & Fannin, 2010; ICOLD, 2017; A. W. Skempton & Brogan, 1994; Thevanayagam & Mohan, 2000). In Figure 17 we show the variation of the proportion of the total mean stress transmitted by the “skeleton” finer fraction, $\alpha_f^{skel}$, for different stress thresholds, where the threshold stress is defined as:

$$ \text{Threshold stress} = \underset{\substack{\text{threshold} \\ \text{coefficient}}}{x} \times \underset{\substack{\text{Average} \\ \text{particle stress}}}{\overline{P^p}} \qquad (10) $$

In Fig. 18, the threshold coefficient, $x$, ranges from 0 to 2.

- Unsurprisingly, the contribution of the fine remains negligible whatever the threshold value for $f_c<f_c^{th}$ (the contribution for all the fines, $\alpha_f$, is already negligible).

- For $f_c^{th} \le f_c \le 45\%$, generally no change is observed in $\alpha_f^{skel}$ until $x = 0.1$, which means that the fine grains transmitting a pressure $P^p$ lower than 10% of the mean pressure $\overline{P^p}$, contribute only marginally to the stress transmitted by the finer fraction.

- For $f_c>45\%$, no change was observed in $\alpha_f^{skel}$ until $x > 0.25$.

From the above analysis, we can retain $x=0.1$ as a reasonable threshold to define the skeleton fines for all fines content.

Fig. 19 shows the proportion of the fines left in an assembly after each exclusion. These “skeleton” fines are included in the solid volume while computing the stress-based skeleton void ratio, $e_{skel}$ plotted in Fig. 20. For $f_c \ge f_c^{th}$, the proportion of the fines included as part of the assembly skeleton decreases monotonically as the threshold stress increases (Fig. 19a-d). The mechanical void ratio, $e_{mech}$ has been computed in the literature as an alternative void ratio index to the void ratio, $e$, of granular materials (Liu et al., 2022; Otsubo, 2016). $e_{mech}$ excludes the particles with $c \le 1$. In Fig. 20, we compare the evolutions of $e$, $e_{mech}$ and $e_{skel}$ (for $x = 0.1$), $e_{skel}$; $x = 0.1$, for the different fines contents. These evolutions are put in parallel with the peak and the critical state stress ratios.

As shown in Figure 19, $e_{skel}$; $x = 0.1$, shows a negative correlation with the critical state strength, $(q/p)_{crit}$ (Fig. 20b & d), and no correlation with the peak strength, $(q/p)_{max}$ (Fig. 20a). Instead, at the peak, $(q/p)_{max}$ correlates negatively with $e$, (Fig. 20a & c), in agreement with established relationship between the initial void ratio and the peak strength in the literature (Adesina et al., 2023; Holubec & D’Appolonia, 1973; Ng, 2004). In order to determine the best performing void ratio index among the indexes considered here, we estimated the goodness of the fit, $R^2$, for the relationship between the void ratio values determined using the indexes and the strength exhibited by the assemblies (Table A.1 in the Appendix C). It is obvious that
$e_{skel}$; $x = 0.1$ provides the best fit ($R^2=0.84$) thereby yielding the best prediction of $(q/p)_{crit}$, while the best performing void ratio index for $(q/p)_{max}$ is $e (R^2=0.90$ for underfilled and $R^2=0.97$ for the overfilled assemblies). This finding is consistent with the fact that the microstructure of granular assemblies (including the force chain network) evolve rapidly under shearing at the critical state (Deng et al., 2022; Wautier et al., 2018). Consequently, the non-active and the marginally active fines do not have the chance to interact with the force chains and should therefore not be considered as part of the void when assessing the mechanical behaviour. On the contrary, at the peak, although the non-active and the marginally active fines do not contribute significantly to the total mean stress, they serve as a support for the force chain network and as a result constitute an important part of the microstructure of the assemblies (Tordesillas et al., 2010; Wautier et al., 2018; Zhu et al., 2016a; Zhu et al., 2016b). Therefore, they have to be taken into account in the void ratio calculation.

<figure><img src="image_21.png" /><figcaption>**Fig. 18** Estimation of the proportion of the mean stress transmitted by the finer fraction when particles with particle stress lower than a threshold are excluded (a) dense; initial state (b) loose; initial state (c) dense; peak state (d) critical state (unique trends for both dense and loose assemblies)</figcaption></figure>
$e_{skel}$; $x = 0.1$ provides the best fit ($R^2=0.84$) thereby yielding the best prediction of $(q/p)_{crit}$, while the best performing void ratio index for $(q/p)_{max}$ is $e (R^2=0.90$ for underfilled and $R^2=0.97$ for the overfilled assemblies). This finding is consistent with the fact that the microstructure of granular assemblies (including the force chain network) evolve rapidly under shearing at the critical state (Deng et al., 2022; Wautier et al., 2018). Consequently, the non-active and the marginally active fines do not have the chance to interact with the force chains and should therefore not be considered as part of the void when assessing the mechanical behaviour. On the contrary, at the peak, although the non-active and the marginally active fines do not contribute significantly to the total mean stress, they serve as a support for the force chain network and as a result constitute an important part of the microstructure of the assemblies (Tordesillas et al., 2010; Wautier et al., 2018; Zhu et al., 2016a; Zhu et al., 2016b). Therefore, they have to be taken into account in the void ratio calculation.

<figure><img src="image_22.png" /><figcaption>**Fig. 18** Estimation of the proportion of the mean stress transmitted by the finer fraction when particles with particle stress lower than a threshold are excluded (a) dense; initial state (b) loose; initial state (c) dense; peak state (d) critical state (unique trends for both dense and loose assemblies)</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 19 Proportion of skeleton fines as a function of the threshold coefficient (a) dense; initial state (b) loose; initial state (c) dense; peak state (d) critical state (unique trends for both dense and loose assemblies)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 20 (a) void ratio indexes at the peak and $(q/p)_{max}$ (b) void ratio indexes at the critical state and $(q/p)_{crit}$, for different $f_c$ values. Relationship between (c) $e$ and $(q/p)_{max}$ (d) $e_{skel}$; $x=0.1$ and $(q/p)_{crit}$.</figcaption></figure>

635

636

637

638

639

640

641

642

643

644

645
<figure><img src="image_23.png" /><figcaption>Fig. 19 Proportion of skeleton fines as a function of the threshold coefficient (a) dense; initial state (b) loose; initial state (c) dense; peak state (d) critical state (unique trends for both dense and loose assemblies)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 20 (a) void ratio indexes at the peak and $(q/p)_{max}$ (b) void ratio indexes at the critical state and $(q/p)_{crit}$, for different $f_c$ values. Relationship between (c) $e$ and $(q/p)_{max}$ (d) $e_{skel}$; $x=0.1$ and $(q/p)_{crit}$.</figcaption></figure>
646 **7. Summary and Conclusions**

647 This study presents a micromechanical evaluation of the regimes delineating the behaviour of
648 gap-graded granular assemblies of different fines content, using discrete element simulations.
649 The existence of the regimes delineated by milestone fines content was assessed using the
650 macromechanical and micromechanical characteristics of the assemblies, and the contributions
651 of the particle size fractions and contact types to the total mean stress. The key findings and
652 conclusions are presented below.

653 (i) Two regimes were identified based on the macroscopic characteristics of the
654 assemblies, three regimes based on particle scale analysis and four regimes based
655 on contact scale analysis. We provide distinct delineations of the regimes by fines
656 content, in an original manner. Fig. 21 presents a summary of all the regimes
657 identified in this study. Table 3 highlights the unique trends characterising the three
658 regimes identified based on the macromechanical and the micromechanical
659 analyses conducted. We showed that the boundaries delineating the identified
660 regimes depend on density and stress state. We acknowledge that the boundaries
661 can also vary with a change in particle size ratio (and more generally a change in
662 coarse and fine particle size distributions) which is beyond the scope of our study.

663 (ii) The micromechanical analyses adopted overcome the limitation of experimental
664 studies where it is difficult to ascertain the claims emanating from the conceptual
665 analysis presented in the studies, since individual particles stresses or particle
666 dispersion cannot be easily determined in experiments. The claims that the fines
667 become the primary stress-transmitting matrix at the threshold fines content or at a
668 limiting fines content; or that stress is equally sheared by the coarse and the fine
669 particles in the transitional zone (Vallejo, 2001) or at f<sub>c</sub>>35% (Shire et al., 2014)
670 are not supported by the evidence provided here. Instead, we found that the coarse
671 particles are dispersed by the fines (i.e. where coarse to coarse contacts are minimal
672 or non-existent) from the threshold fines content. In addition, the fines do not
673 contribute to the total mean stress below the threshold fines content, they play a
674 secondary role in the transitional zone (f<sub>c</sub> ∈ [30%; 55%] depending on density and
675 stress state), and only become the primary stress-transmitting matrix beyond the
676 transitional zone (i.e. when f<sub>c</sub> > f<sub>c</sub><sup>eq</sup>; f<sub>c</sub><sup>eq</sup> ∈ [48%; 55%] depending on density and
677 stress state).
646 **7. Summary and Conclusions**

647 This study presents a micromechanical evaluation of the regimes delineating the behaviour of
648 gap-graded granular assemblies of different fines content, using discrete element simulations.
649 The existence of the regimes delineated by milestone fines content was assessed using the
650 macromechanical and micromechanical characteristics of the assemblies, and the contributions
651 of the particle size fractions and contact types to the total mean stress. The key findings and
652 conclusions are presented below.

653 (i) Two regimes were identified based on the macroscopic characteristics of the
654 assemblies, three regimes based on particle scale analysis and four regimes based
655 on contact scale analysis. We provide distinct delineations of the regimes by fines
656 content, in an original manner. Fig. 21 presents a summary of all the regimes
657 identified in this study. Table 3 highlights the unique trends characterising the three
658 regimes identified based on the macromechanical and the micromechanical
659 analyses conducted. We showed that the boundaries delineating the identified
660 regimes depend on density and stress state. We acknowledge that the boundaries
661 can also vary with a change in particle size ratio (and more generally a change in
662 coarse and fine particle size distributions) which is beyond the scope of our study.

663 (ii) The micromechanical analyses adopted overcome the limitation of experimental
664 studies where it is difficult to ascertain the claims emanating from the conceptual
665 analysis presented in the studies, since individual particles stresses or particle
666 dispersion cannot be easily determined in experiments. The claims that the fines
667 become the primary stress-transmitting matrix at the threshold fines content or at a
668 limiting fines content; or that stress is equally sheared by the coarse and the fine
669 particles in the transitional zone (Vallejo, 2001) or at f<sub>c</sub>>35% (Shire et al., 2014)
670 are not supported by the evidence provided here. Instead, we found that the coarse
671 particles are dispersed by the fines (i.e. where coarse to coarse contacts are minimal
672 or non-existent) from the threshold fines content. In addition, the fines do not
673 contribute to the total mean stress below the threshold fines content, they play a
674 secondary role in the transitional zone (f<sub>c</sub> ∈ [30%; 55%] depending on density and
675 stress state), and only become the primary stress-transmitting matrix beyond the
676 transitional zone (i.e. when f<sub>c</sub> > f<sub>c</sub><sup>eq</sup>; f<sub>c</sub><sup>eq</sup> ∈ [48%; 55%] depending on density and
677 stress state).
(iii) We found that the threshold stress which determines whether a fine particle transmits a marginal or a significant stress with respect to the contribution of the fines to the total mean stress is a fraction of the average particle stress, $\overline{Pp}$, within an assembly; where the threshold stress values are $0.1 \times \overline{Pp}$ for $f_c^{th} \le f_c \le 45\%$, and $0.25 \times \overline{Pp}$ for $f_c > 45\%$. The threshold stress concept proposed in this study can be useful in determining the particles that constitute important mesostructures (i.e. important force-chain networks) in gap-graded granular assemblies. The skeleton void ratio determined based on the threshold stress correlates with the critical state strength of the assemblies, while no correlation was found between the critical state strength and alternative void ratio indexes such as the void ratio and the mechanical void ratio at the critical state. The marginal stress transmitting fines (marginally active fines) determined based on the threshold stress may constitute a part of the particles susceptible to internal erosion, in addition to the non-inactive fines in a gap-graded assembly (particles with zero or only one contact).

In addition to the standard underfilled and overfilled regimes, future DEM studies could focus on the two transitional regimes identified here, and seek to understand the role played by marginally active fines i) in the susceptibility of gap-graded assemblies to internal erosion and ii) in the triggering of mechanical instabilities. As an additional perspective, the distinct characteristics of the regimes shown in this study will prove useful in the development of micromechanical models for gap-graded materials in which the typical mesoscale grain arrangements can be tailored for each regime.

<table><caption>Fig. 21 Regime identification by macroscale, particle scale and contact scale characteristics</caption><thead><tr><th>Macroscale</th><th>Underfilled</th><th colspan="2">Overfilled</th></tr></thead><tbody><tr><td>Particle scale</td><td>α<sub>f</sub> ≈ 0<br>Coarse-dominated</td><td>α<sub>f</sub> &lt; α<sub>c</sub> (The coarse contribute primarily to total mean stress while the fines play a secondary role)</td><td>α<sub>f</sub> &gt; α<sub>c</sub> (The fines play a primary role)</td></tr><tr><td>Contact scale</td><td>α<sub>c-c</sub> &gt; α<sub>c-f</sub> &gt; α<sub>f-f</sub><br>Coarse-dominated</td><td>α<sub>c-f</sub> &gt; α<sub>c-c</sub> &gt; α<sub>f-f</sub><br>Transitional<br>coarse-dominated</td><td>α<sub>c-f</sub> &gt; α<sub>f-f</sub> &gt; α<sub>c-c</sub><br>Transitional<br>fines-dominated</td></tr></tbody></table>
(iii) We found that the threshold stress which determines whether a fine particle transmits a marginal or a significant stress with respect to the contribution of the fines to the total mean stress is a fraction of the average particle stress, $\overline{Pp}$, within an assembly; where the threshold stress values are $0.1 \times \overline{Pp}$ for $f_c^{th} \le f_c \le 45\%$, and $0.25 \times \overline{Pp}$ for $f_c > 45\%$. The threshold stress concept proposed in this study can be useful in determining the particles that constitute important mesostructures (i.e. important force-chain networks) in gap-graded granular assemblies. The skeleton void ratio determined based on the threshold stress correlates with the critical state strength of the assemblies, while no correlation was found between the critical state strength and alternative void ratio indexes such as the void ratio and the mechanical void ratio at the critical state. The marginal stress transmitting fines (marginally active fines) determined based on the threshold stress may constitute a part of the particles susceptible to internal erosion, in addition to the non-inactive fines in a gap-graded assembly (particles with zero or only one contact).

In addition to the standard underfilled and overfilled regimes, future DEM studies could focus on the two transitional regimes identified here, and seek to understand the role played by marginally active fines i) in the susceptibility of gap-graded assemblies to internal erosion and ii) in the triggering of mechanical instabilities. As an additional perspective, the distinct characteristics of the regimes shown in this study will prove useful in the development of micromechanical models for gap-graded materials in which the typical mesoscale grain arrangements can be tailored for each regime.

<table><caption>Fig. 21 Regime identification by macroscale, particle scale and contact scale characteristics</caption><thead><tr><th>Macroscale</th><th>Underfilled</th><th colspan="2">Overfilled</th></tr></thead><tbody><tr><td>Particle scale</td><td>α<sub>f</sub> ≈ 0<br>Coarse-dominated</td><td>α<sub>f</sub> &lt; α<sub>c</sub> (The coarse contribute primarily to total mean stress while the fines play a secondary role)</td><td>α<sub>f</sub> &gt; α<sub>c</sub> (The fines play a primary role)</td></tr><tr><td>Contact scale</td><td>α<sub>c-c</sub> &gt; α<sub>c-f</sub> &gt; α<sub>f-f</sub><br>Coarse-dominated</td><td>α<sub>c-f</sub> &gt; α<sub>c-c</sub> &gt; α<sub>f-f</sub><br>Transitional<br>coarse-dominated</td><td>α<sub>c-f</sub> &gt; α<sub>f-f</sub> &gt; α<sub>c-c</sub><br>Transitional<br>fines-dominated</td></tr></tbody></table>
712 Notations

713 C Total number of contacts in a granular system

714 $C_{c-c}^t$ The proportion of all contact types belonging to the c-c contacts

715 $C_{c-f}^t$ The proportion of all contact types belonging to the c-f contacts

716 $C_{f-f}^t$ The proportion of all contact types belonging to the f-f contacts

717 $C_{wl}^t$ The proportion of all contact types belonging to the wall

718 $C_n^f$ Proportion of the total number contacts belonging to the finer fraction

719 $C_n^c$ Proportion of the total number contacts belonging to the coarser fraction

720 $d_{50}$ Mean size of the fine grains

721 $D_{50}$ Mean size of the coarse grains

722 $e_{max}$ Maximum void ratio

723 $e_{min}$ Minimum void ratio

724 $e_{crit}$ Void ratio at the critical state

725 $e_{hc}$ Void ratio of the host sand

726 $e_{hf}$ Void ratio of the silt

727 $\varepsilon_a$ Axial strain (in y-direction)

728 $\varepsilon_q$ Deviatoric strain

729 $\varepsilon_{vol}$ Volumetric strain

730 $f_c$ Fines content

731 $f_{th}$ Threshold fines content

732 $f_c^{th}$ Threshold fines content

733 $f_c^L$ Limiting fines content

734 $f_j^c$ Contact force vector

735 $l$ Unit branch vector

736 $N_p$ Number of particles in the granular system

737 $N_c$ Number of contacts in a granular system

738 $N_p^c$ Number of coarse particles in the system

739 $N_p^f$ Number of fine particles in the system

740 $N_c^p$ Number of contacts involving a particle

741 $N_0$ Number of particles with no contact

742 $N_1$ Number of particles with one contact

743 $N_{c\ge1}$ Number of coarse particles with one or more contacts

744 $N_{f\ge1}$ Number of fine particles with one or more contacts

745 $NC_f$ Normalised total contacts for the coarser fraction

746 $NC_f$ Normalised total contacts for the finer fraction

747 $n_i$ Unit branch vector

748 $p$ Mean effective stress

749 $p^c$ Mean stress transmitted by the coarser fraction

750 $p^f$ Mean stress transmitted by the finer fraction

751 $p^k$ Mean stress transmitted by a contact type

752 $p^p$ Particle mean stress

753 $\overline{p^p}$ Average particle mean stress

754 $\alpha_c$ Proportion of the total mean stress transmitted by the coarser fraction

755 $\alpha_f$ Proportion of the total mean stress transmitted by the finer fraction

756 $\alpha_k$ Proportion of the total mean stress transmitted a contact type

757 $(q/p)_{max}$ Stress ratio at the peak

758 $(q/p)_{crit}$ Stress ratio at the critical state
712 Notations

713 $C$ Total number of contacts in a granular system

714 $C_{c-c}^t$ The proportion of all contact types belonging to the c-c contacts

715 $C_{c-f}^t$ The proportion of all contact types belonging to the c-f contacts

716 $C_{f-f}^t$ The proportion of all contact types belonging to the f-f contacts

717 $C_{wl}^t$ The proportion of all contact types belonging to the wall

718 $C_n^f$ Proportion of the total number contacts belonging to the finer fraction

719 $C_n^c$ Proportion of the total number contacts belonging to the coarser fraction

720 $d_{50}$ Mean size of the fine grains

721 $D_{50}$ Mean size of the coarse grains

722 $e_{max}$ Maximum void ratio

723 $e_{min}$ Minimum void ratio

724 $e_{crit}$ Void ratio at the critical state

725 $e_{hc}$ Void ratio of the host sand

726 $e_{hf}$ Void ratio of the silt

727 $\varepsilon_a$ Axial strain (in y-direction)

728 $\varepsilon_q$ Deviatoric strain

729 $\varepsilon_{vol}$ Volumetric strain

730 $f_c$ Fines content

731 $f_{th}$ Threshold fines content

732 $f_c^{th}$ Threshold fines content

733 $f_c^L$ Limiting fines content

734 $f_j^c$ Contact force vector

735 $l$ Unit branch vector

736 $N_p$ Number of particles in the granular system

737 $N_c$ Number of contacts in a granular system

738 $N_p^c$ Number of coarse particles in the system

739 $N_p^f$ Number of fine particles in the system

740 $N_c^p$ Number of contacts involving a particle

741 $N_0$ Number of particles with no contact

742 $N_1$ Number of particles with one contact

743 $N_{c\ge1}$ Number of coarse particles with one or more contacts

744 $N_{f\ge1}$ Number of fine particles with one or more contacts

745 $NC_f$ Normalised total contacts for the coarser fraction

746 $NC_f$ Normalised total contacts for the finer fraction

747 $n_i$ Unit branch vector

748 $p$ Mean effective stress

749 $p^c$ Mean stress transmitted by the coarser fraction

750 $p^f$ Mean stress transmitted by the finer fraction

751 $p^k$ Mean stress transmitted by a contact type

752 $p^p$ Particle mean stress

753 $\overline{p^p}$ Average particle mean stress

754 $\alpha_c$ Proportion of the total mean stress transmitted by the coarser fraction

755 $\alpha_f$ Proportion of the total mean stress transmitted by the finer fraction

756 $\alpha_k$ Proportion of the total mean stress transmitted a contact type

757 $(q/p)_{max}$ Stress ratio at the peak

758 $(q/p)_{crit}$ Stress ratio at the critical state
<table>
  <tr>
    <td>759</td>
    <td>q</td>
    <td>Deviatoric stress</td>
  </tr>
  <tr>
    <td>760</td>
    <td>q<sup>f</sup></td>
    <td>Deviatoric stress transmitted by the finer fraction</td>
  </tr>
  <tr>
    <td>761</td>
    <td>q<sup>c</sup></td>
    <td>Deviatoric stress transmitted by the coarser fraction</td>
  </tr>
  <tr>
    <td>762</td>
    <td>q<sub>f</sub></td>
    <td>Proportion of the deviatoric stress transmitted by the finer fraction</td>
  </tr>
  <tr>
    <td>763</td>
    <td>R-squared</td>
    <td>Goodness of fit</td>
  </tr>
  <tr>
    <td>764</td>
    <td>R<sub>d</sub></td>
    <td>Ratio of the mean size of the coarse grains to the mean size of the fine grains</td>
  </tr>
  <tr>
    <td>765</td>
    <td>R<sub>p</sub></td>
    <td>Radius of particle</td>
  </tr>
  <tr>
    <td>766</td>
    <td>R<sub>1</sub></td>
    <td>Radius of first particle in contact</td>
  </tr>
  <tr>
    <td>767</td>
    <td>R<sub>2</sub></td>
    <td>Radius of second particle in contact</td>
  </tr>
  <tr>
    <td>768</td>
    <td>R<sub>f</sub></td>
    <td>Proportion of non-active fines</td>
  </tr>
  <tr>
    <td>769</td>
    <td>e<sub>skel</sub></td>
    <td>Skeleton void ratio</td>
  </tr>
  <tr>
    <td>770</td>
    <td>e<sub>skel</sub>;x</td>
    <td>Skeleton void ratio determined using a threshold coefficient</td>
  </tr>
  <tr>
    <td>771</td>
    <td>Skl<sub>f</sub></td>
    <td>Proportion of fines included in skeleton</td>
  </tr>
  <tr>
    <td>772</td>
    <td>x</td>
    <td>Threshold coefficient</td>
  </tr>
  <tr>
    <td>773</td>
    <td>μ<sub>prep</sub></td>
    <td>Inter-particle friction coefficient during assembly preparation (isotropic compression)</td>
  </tr>
  <tr>
    <td>774</td>
    <td>μ</td>
    <td>Inter-particle friction coefficient</td>
  </tr>
  <tr>
    <td>775</td>
    <td>λ</td>
    <td>Size fraction</td>
  </tr>
  <tr>
    <td>776</td>
    <td>V</td>
    <td>Volume of system</td>
  </tr>
  <tr>
    <td>777</td>
    <td>V<sup>p</sup></td>
    <td>Volume of particle</td>
  </tr>
  <tr>
    <td>778</td>
    <td>σ<sub>zz</sub></td>
    <td>Normal stress in z-direction</td>
  </tr>
  <tr>
    <td>779</td>
    <td>σ<sub>yy</sub></td>
    <td>Normal stress in y-direction</td>
  </tr>
  <tr>
    <td>780</td>
    <td>σ<sub>xx</sub></td>
    <td>Normal stress in x-direction</td>
  </tr>
  <tr>
    <td>781</td>
    <td>σ<sub>ij</sub></td>
    <td>Second order stress tensor</td>
  </tr>
  <tr>
    <td>782</td>
    <td>σ<sub>ij</sub><sup>c</sup></td>
    <td>Second order stress tensor for all contacts</td>
  </tr>
  <tr>
    <td>783</td>
    <td>σ<sub>ij</sub><sup>p</sup></td>
    <td>Second order stress tensor for a particle</td>
  </tr>
  <tr>
    <td>784</td>
    <td>σ<sub>ij</sub><sup>c</sup></td>
    <td>Average mean stress transmitted by the coarser fraction</td>
  </tr>
  <tr>
    <td>785</td>
    <td>σ<sub>p</sub><sup>f</sup></td>
    <td>Average mean stress transmitted by the finer fraction</td>
  </tr>
  <tr>
    <td>786</td>
    <td>σ<sub>p</sub><sup>f</sup></td>
    <td>Average mean stress transmitted by the finer fraction</td>
  </tr>
  <tr>
    <td>787</td>
    <td>φ</td>
    <td>Friction angle</td>
  </tr>
  <tr>
    <td>788</td>
    <td>φ<sub>p</sub></td>
    <td>Peak friction angle</td>
  </tr>
  <tr>
    <td>789</td>
    <td>φ<sub>c</sub></td>
    <td>Critical state friction angle</td>
  </tr>
  <tr>
    <td>790</td>
    <td>ψ<sub>max</sub></td>
    <td>Maximum dilatancy angle</td>
  </tr>
  <tr>
    <td>791</td>
    <td>Z</td>
    <td>Mean coordination number</td>
  </tr>
  <tr>
    <td>792</td>
    <td>Z<sub>ini</sub></td>
    <td>Mean coordination number at the initial state (after isotropic compression)</td>
  </tr>
  <tr>
    <td>793</td>
    <td>Z<sub>crit</sub></td>
    <td>Mean coordination number at the critical state</td>
  </tr>
  <tr>
    <td>794</td>
    <td>Z<sub>g</sub></td>
    <td>Geometrical coordination number</td>
  </tr>
  <tr>
    <td>795</td>
    <td>Z<sub>gc</sub></td>
    <td>Geometrical coordination number for coarse particles</td>
  </tr>
  <tr>
    <td>796</td>
    <td>Z<sub>gf</sub></td>
    <td>Geometrical coordination number for fine particles</td>
  </tr>
  <tr>
    <td>797</td>
    <td>Z<sub>m</sub></td>
    <td>Mechanical coordination number</td>
  </tr>
  <tr>
    <td>798</td>
    <td>Z<sub>mc</sub></td>
    <td>Mechanical coordination number for coarse particles</td>
  </tr>
  <tr>
    <td>799</td>
    <td>Z<sub>mf</sub></td>
    <td>Mechanical coordination number for fine particles</td>
  </tr>
  <tr>
    <td>800</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>801</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>802</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>803</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>804</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>805</td>
    <td></td>
    <td></td>
  </tr>
</table>
<table>
  <tr>
    <td>759</td>
    <td>q</td>
    <td>Deviatoric stress</td>
  </tr>
  <tr>
    <td>760</td>
    <td>q<sup>f</sup></td>
    <td>Deviatoric stress transmitted by the finer fraction</td>
  </tr>
  <tr>
    <td>761</td>
    <td>q<sup>c</sup></td>
    <td>Deviatoric stress transmitted by the coarser fraction</td>
  </tr>
  <tr>
    <td>762</td>
    <td>q<sub>f</sub></td>
    <td>Proportion of the deviatoric stress transmitted by the finer fraction</td>
  </tr>
  <tr>
    <td>763</td>
    <td>R-squared</td>
    <td>Goodness of fit</td>
  </tr>
  <tr>
    <td>764</td>
    <td>R<sub>d</sub></td>
    <td>Ratio of the mean size of the coarse grains to the mean size of the fine grains</td>
  </tr>
  <tr>
    <td>765</td>
    <td>R<sub>p</sub></td>
    <td>Radius of particle</td>
  </tr>
  <tr>
    <td>766</td>
    <td>R<sub>1</sub></td>
    <td>Radius of first particle in contact</td>
  </tr>
  <tr>
    <td>767</td>
    <td>R<sub>2</sub></td>
    <td>Radius of second particle in contact</td>
  </tr>
  <tr>
    <td>768</td>
    <td>R<sub>f</sub></td>
    <td>Proportion of non-active fines</td>
  </tr>
  <tr>
    <td>769</td>
    <td>e<sub>skel</sub></td>
    <td>Skeleton void ratio</td>
  </tr>
  <tr>
    <td>770</td>
    <td>e<sub>skel</sub>;x</td>
    <td>Skeleton void ratio determined using a threshold coefficient</td>
  </tr>
  <tr>
    <td>771</td>
    <td>Skl<sub>f</sub></td>
    <td>Proportion of fines included in skeleton</td>
  </tr>
  <tr>
    <td>772</td>
    <td>x</td>
    <td>Threshold coefficient</td>
  </tr>
  <tr>
    <td>773</td>
    <td>μ<sub>prep</sub></td>
    <td>Inter-particle friction coefficient during assembly preparation (isotropic compression)</td>
  </tr>
  <tr>
    <td>774</td>
    <td>μ</td>
    <td>Inter-particle friction coefficient</td>
  </tr>
  <tr>
    <td>775</td>
    <td>λ</td>
    <td>Size fraction</td>
  </tr>
  <tr>
    <td>776</td>
    <td>V</td>
    <td>Volume of system</td>
  </tr>
  <tr>
    <td>777</td>
    <td>V<sup>p</sup></td>
    <td>Volume of particle</td>
  </tr>
  <tr>
    <td>778</td>
    <td>σ<sub>zz</sub></td>
    <td>Normal stress in z-direction</td>
  </tr>
  <tr>
    <td>779</td>
    <td>σ<sub>yy</sub></td>
    <td>Normal stress in y-direction</td>
  </tr>
  <tr>
    <td>780</td>
    <td>σ<sub>xx</sub></td>
    <td>Normal stress in x-direction</td>
  </tr>
  <tr>
    <td>781</td>
    <td>σ<sub>ij</sub></td>
    <td>Second order stress tensor</td>
  </tr>
  <tr>
    <td>782</td>
    <td>σ<sub>ij</sub><sup>c</sup></td>
    <td>Second order stress tensor for all contacts</td>
  </tr>
  <tr>
    <td>783</td>
    <td>σ<sub>ij</sub><sup>p</sup></td>
    <td>Second order stress tensor for a particle</td>
  </tr>
  <tr>
    <td>784</td>
    <td>σ<sub>ij</sub><sup>c</sup></td>
    <td>Average mean stress transmitted by the coarser fraction</td>
  </tr>
  <tr>
    <td>785</td>
    <td>σ<sub>p</sub><sup>f</sup></td>
    <td>Average mean stress transmitted by the finer fraction</td>
  </tr>
  <tr>
    <td>786</td>
    <td>σ<sub>p</sub><sup>f</sup></td>
    <td>Average mean stress transmitted by the finer fraction</td>
  </tr>
  <tr>
    <td>787</td>
    <td>φ</td>
    <td>Friction angle</td>
  </tr>
  <tr>
    <td>788</td>
    <td>φ<sub>p</sub></td>
    <td>Peak friction angle</td>
  </tr>
  <tr>
    <td>789</td>
    <td>φ<sub>c</sub></td>
    <td>Critical state friction angle</td>
  </tr>
  <tr>
    <td>790</td>
    <td>ψ<sub>max</sub></td>
    <td>Maximum dilatancy angle</td>
  </tr>
  <tr>
    <td>791</td>
    <td>Z</td>
    <td>Mean coordination number</td>
  </tr>
  <tr>
    <td>792</td>
    <td>Z<sub>ini</sub></td>
    <td>Mean coordination number at the initial state (after isotropic compression)</td>
  </tr>
  <tr>
    <td>793</td>
    <td>Z<sub>crit</sub></td>
    <td>Mean coordination number at the critical state</td>
  </tr>
  <tr>
    <td>794</td>
    <td>Z<sub>g</sub></td>
    <td>Geometrical coordination number</td>
  </tr>
  <tr>
    <td>795</td>
    <td>Z<sub>gc</sub></td>
    <td>Geometrical coordination number for coarse particles</td>
  </tr>
  <tr>
    <td>796</td>
    <td>Z<sub>gf</sub></td>
    <td>Geometrical coordination number for fine particles</td>
  </tr>
  <tr>
    <td>797</td>
    <td>Z<sub>m</sub></td>
    <td>Mechanical coordination number</td>
  </tr>
  <tr>
    <td>798</td>
    <td>Z<sub>mc</sub></td>
    <td>Mechanical coordination number for coarse particles</td>
  </tr>
  <tr>
    <td>799</td>
    <td>Z<sub>mf</sub></td>
    <td>Mechanical coordination number for fine particles</td>
  </tr>
  <tr>
    <td>800</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>801</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>802</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>803</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>804</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>805</td>
    <td></td>
    <td></td>
  </tr>
</table>
## Appendix A. REV scale analysis

In order to confirm that the sample size of 100,000 particles used for all $f_c$ in this study is a representative element volume (REV), in Fig. A.1, we show that increasing the assembly size to 200,000 particles for both the dense and loose assemblies having $f_c$ = 30% did not significantly influence the stress-strain responses. The slight sensitivity to assembly size observed in the volumetric strains was also reported in the DEM study on REV for granular materials by Adesina et al., (2022) and can be attributed to the small variation in the initial void ratio of the assemblies. The REV test is shown here for $f_c$ = 30% because it represents the threshold fines content and is the fines content with the highest coarse-fine particle interaction (Fig. 17) which indicates the largest REV size.

<figure><img src="image_24.png" /><figcaption>**Fig. A.1** Effect of assembly size on (i) the stress-strain responses for dense assemblies ($\mu_{prep}$=0.03) and (ii) loose assemblies ($\mu_{prep}$=0.5) (iii) volumetric strain responses for the dense assemblies and (iv) the loose assemblies having $f_c$ = 30%</figcaption></figure>
## Appendix A. REV scale analysis

In order to confirm that the sample size of 100,000 particles used for all $f_c$ in this study is a representative element volume (REV), in Fig. A.1, we show that increasing the assembly size to 200,000 particles for both the dense and loose assemblies having $f_c$ = 30% did not significantly influence the stress-strain responses. The slight sensitivity to assembly size observed in the volumetric strains was also reported in the DEM study on REV for granular materials by Adesina et al., (2022) and can be attributed to the small variation in the initial void ratio of the assemblies. The REV test is shown here for $f_c$ = 30% because it represents the threshold fines content and is the fines content with the highest coarse-fine particle interaction (Fig. 17) which indicates the largest REV size.

<figure><img src="image_25.png" /><figcaption>**Fig. A.1** Effect of assembly size on (i) the stress-strain responses for dense assemblies ($\mu_{prep}$=0.03) and (ii) loose assemblies ($\mu_{prep}$=0.5) (iii) volumetric strain responses for the dense assemblies and (iv) the loose assemblies having $f_c$ = 30%</figcaption></figure>

806

807

808

809

810

811

812

813

814

815

816

817

818

819

820

821

822

823

824

825

826
## Appendix B. Fine contribution to deviatoric stress

Fig. A.2 shows the evolution of the proportion of the deviatoric stress, $q$, transmitted by the fines during shearing, for dense and loose assemblies of different fines content, $f_c$. Provided that all stress tensors are diagonal in the frame (x, y, z), the deviatoric stress for the entire system is calculated as:

$$q = \sigma_{zz} - (\sigma_{xx} + \sigma_{yy})/2 \quad (A1)$$

Similar to the mean stress (see Section 4), the proportion of the total deviatoric stress transmitted by the finer fraction, $q_f$, is calculated as $q^f/q$, while the proportion of the total deviatoric stress transmitted by the coarser fraction, $q_f$, is calculated as $q^c/q$, where $q^f$ and $q^c$ are the deviatoric stress transmitted by the finer and the coarser populations, respectively.

As expected, the $q_f$ data is generally similar to the $\alpha_f$ data presented in Fig. 9. It is worthy of note that the $q_f$ data presented in Fig. A2 starts after the initial state, at $\epsilon_a = 0.02$, since $q = 0$ at the initial state. We observed that for $f_c \ge 30\%$, the magnitude of $q_f$ is lower than the magnitude of $\alpha_f$ during shearing. This indicates that the stress is more isotropic in the finer fraction than in the coarser fraction.

<figure><img src="image_25.png" /><figcaption>Fig. A.2 Evolution of the proportion of deviatoric stress transmitted by the fines during shearing for (a) dense assemblies (b) loose assemblies</figcaption></figure>
## Appendix B. Fine contribution to deviatoric stress

Fig. A.2 shows the evolution of the proportion of the deviatoric stress, $q$, transmitted by the fines during shearing, for dense and loose assemblies of different fines content, $f_c$. Provided that all stress tensors are diagonal in the frame $(x, y, z)$, the deviatoric stress for the entire system is calculated as:

$$q = \sigma_{zz} - (\sigma_{xx} + \sigma_{yy})/2 \quad (A1)$$

Similar to the mean stress (see Section 4), the proportion of the total deviatoric stress transmitted by the finer fraction, $q_f$, is calculated as $q^f/q$, while the proportion of the total deviatoric stress transmitted by the coarser fraction, $q_c$, is calculated as $q^c/q$, where $q^f$ and $q^c$ are the deviatoric stress transmitted by the finer and the coarser populations, respectively.

As expected, the $q_f$ data is generally similar to the $\alpha_f$ data presented in Fig. 9. It is worthy of note that the $q_f$ data presented in Fig. A2 starts after the initial state, at $\epsilon_a = 0.02$, since $q = 0$ at the initial state. We observed that for $f_c \ge 30\%$, the magnitude of $q_f$ is lower than the magnitude of $\alpha_f$ during shearing. This indicates that the stress is more isotropic in the finer fraction than in the coarser fraction.

<figure><img src="image_26.png" /><figcaption>Fig. A.2 Evolution of the proportion of deviatoric stress transmitted by the fines during shearing for (a) dense assemblies (b) loose assemblies</figcaption></figure>
## Appendix C. Sensitivity analysis to the threshold coefficient used in the definition of the skeleton void ratio.

In Table A.1 and Fig. A.3, we show that at the critical state, $e_{skel}$; $x = 0.1$ exhibited a negative correlation with $(q/p)_{crit}$ and yielded the highest R-squared value ($R^2=0.84$), while both the $e$ and $e_{mech}$ do not correlate with $(q/p)_{crit}$. $e_{skel}$; $x = 0.05$ and $e_{skel}$; $x = 0.25$ also yielded a negative correlation with $(q/p)_{crit}$, with $R^2=0.73$ and $R^2=0.76$, respectively. $e_{skel}$; $x = 1$ yielded a poor correlation ($R^2=0.08$) with $(q/p)_{crit}$, indicating that using a threshold coefficient, $x = 1$, for the skeleton void ratio calculation for the gap-graded assemblies, resulted in the exclusion of important stress-transmitting fines, especially from the overfilled assemblies. $x = 1$ is the standard for defining force-chains for narrowly-graded assemblies (Peters et al., 2005). In agreement with earlier studies (Adesina et al., 2023; Holubec & D’Appolonia, 1973; Ng, 2004), at the peak, $e$, correlates negatively with $(q/p)_{max}$, and yielded the highest R-squared value ($R^2=0.97$), followed by $e_{mech}$ ($R^2=0.66$). $e_{skel}$ correlated poorly with $(q/p)_{max}$.

<table><caption>**Table A.1** Assessment of void ratio indexes with the strength value, $(q/p)$</caption><thead><tr><th>Void ratio index</th><th>e</th><th>e<sub>mech</sub></th><th>e<sub>skel</sub>; x=0.05</th><th>e<sub>skel</sub>; x=0.10</th><th>e<sub>skel</sub>; x=0.25</th><th>e<sub>skel</sub>; x=1.0</th></tr></thead><tbody><tr><td>R<sup>2</sup><sub>peak</sub></td><td>0.965</td><td>0.656</td><td>0.119</td><td>0.074</td><td>-</td><td>9.55e-5</td></tr><tr><td>R<sup>2</sup><sub>crit</sub></td><td>-</td><td>0.208</td><td>0.725</td><td>0.844</td><td>0.764</td><td>0.08</td></tr></tbody></table>

Note: The $R^2$ values here are determined from inverse relationship between strength and void ratio. ‘-’ means strength increased with void ratio as against established relationship between strength and void ratio.

<figure><img src="image_26.png" /><figcaption>**Fig. A.3** Relationship between (a) void ratio indexes and peak $q/p$, $(q/p)_{max}$ (b) void ratio indexes and critical state $q/p$, $(q/p)_{crit}$</figcaption></figure>
## Appendix C. Sensitivity analysis to the threshold coefficient used in the definition of the skeleton void ratio.

In Table A.1 and Fig. A.3, we show that at the critical state, $e_{skel}$; $x = 0.1$ exhibited a negative correlation with $(q/p)_{crit}$ and yielded the highest R-squared value ($R^2=0.84$), while both the $e$ and $e_{mech}$ do not correlate with $(q/p)_{crit}$. $e_{skel}$; $x = 0.05$ and $e_{skel}$; $x = 0.25$ also yielded a negative correlation with $(q/p)_{crit}$, with $R^2=0.73$ and $R^2=0.76$, respectively. $e_{skel}$; $x = 1$ yielded a poor correlation ($R^2=0.08$) with $(q/p)_{crit}$, indicating that using a threshold coefficient, $x = 1$, for the skeleton void ratio calculation for the gap-graded assemblies, resulted in the exclusion of important stress-transmitting fines, especially from the overfilled assemblies. $x = 1$ is the standard for defining force-chains for narrowly-graded assemblies (Peters et al., 2005). In agreement with earlier studies (Adesina et al., 2023; Holubec & D’Appolonia, 1973; Ng, 2004), at the peak, $e$, correlates negatively with $(q/p)_{max}$, and yielded the highest R-squared value ($R^2=0.97$), followed by $e_{mech}$ ($R^2=0.66$). $e_{skel}$ correlated poorly with $(q/p)_{max}$.

<table><caption>**Table A.1** Assessment of void ratio indexes with the strength value, $(q/p)$</caption><thead><tr><th>Void ratio index</th><th>e</th><th>e<sub>mech</sub></th><th>e<sub>skel</sub>; x=0.05</th><th>e<sub>skel</sub>; x=0.10</th><th>e<sub>skel</sub>; x=0.25</th><th>e<sub>skel</sub>; x=1.0</th></tr></thead><tbody><tr><td>R<sup>2</sup><sub>peak</sub></td><td>0.965</td><td>0.656</td><td>0.119</td><td>0.074</td><td>-</td><td>9.55e-5</td></tr><tr><td>R<sup>2</sup><sub>crit</sub></td><td>-</td><td>0.208</td><td>0.725</td><td>0.844</td><td>0.764</td><td>0.08</td></tr></tbody></table>

Note: The $R^2$ values here are determined from inverse relationship between strength and void ratio. ‘-’ means strength increased with void ratio as against established relationship between strength and void ratio.

<figure><img src="image_27.png" /><figcaption>**Fig. A.3** Relationship between (a) void ratio indexes and peak $q/p$, $(q/p)_{max}$ (b) void ratio indexes and critical state $q/p$, $(q/p)_{crit}$</figcaption></figure>
875 Acknowledgement

876 This work was supported by the National Research Agency (ANR), France. It is a product of
877 the research project titled “Micro-scale driven Models for Discrete for Discrete Materials
878 (MiMoDiM)”, with the reference, ANR-22-CE46-0001. A CC-BY public copyright license has
879 been applied by the authors to the present document and will be applied to all subsequent
880 versions up to the Author Accepted Manuscript arising from this submission, in accordance
881 with the grant’s open access conditions.

882

883 Credit Authorship

884 Peter Adesina: Conceptualization, Methodology, Investigation, Formal analysis,
885 Visualization, Writing – original draft.

886 Antoine Wautier: Funding acquisition, Conceptualization, Methodology, Supervision,
887 Writing – review & editing.

888 Nadia Benahmed: Conceptualization, Methodology, Supervision, Writing – review & editing.

890

891 Compliance with ethical standards

892 **Conflict of interest** The authors declare that there is no conflict of interests regarding the
893 publication of this article. Publication has been approved by all authors. None of the material
894 presented in the paper is submitted or published elsewhere.

895

896

897 References

898

899 Adesina, P. A., O'Sullivan, C., & Wang, T. (2023). DEM study on the effect of particle shape
900 on the shear behaviour of granular materials. *Computational Particle Mechanics*,
901 *Accepted*.

902 Adesina, P., O'Sullivan, C., Morimoto, T., & Otsubo, M. (2022). Determining a
903 representative element volume for DEM simulations of samples with non-circular
904 particles. *Particuology*, **68**, 29–43. https://doi.org/10.1016/j.partic.2021.10.007

905 Adesina, P., O'Sullivan, C., & Wang, T. (2024). Understanding the interplay between particle
906 shape, grading and sample density on the behaviour of granular assemblies: A DEM
907 approach. *Granular Matter*, **26**(1). https://doi.org/10.1007/s10035-023-01383-2

908 Bagi, K. (1999). Microstructural stress tensor of granular assemblies with volume forces.
909 *Journal of Applied Mechanics, Transactions ASME*, **66**(4), 934–936.
910 https://doi.org/10.1115/1.2791800

911 Benahmed, N., Nguyen, T. K., Hicher, P. Y., & Nicolas, M. (2015). An experimental
912 investigation into the effects of low plastic fines content on the behaviour of sand/silt
913 mixtures. *European Journal of Environmental and Civil Engineering*, **19**(1), 109–128.
914 https://doi.org/10.1080/19648189.2014.939304

915 Bolton, M. D. (1986). The strength and dilatancy of sands. *Geotechnique*, **36**(1), 65–78.

916 Chang, C. S., & Yin, Z. Y. (2011). Micromechanical modeling for behavior of silty sand with
875 Acknowledgement

876 This work was supported by the National Research Agency (ANR), France. It is a product of
877 the research project titled “Micro-scale driven Models for Discrete for Discrete Materials
878 (MiMoDiM)”, with the reference, ANR-22-CE46-0001. A CC-BY public copyright license has
879 been applied by the authors to the present document and will be applied to all subsequent
880 versions up to the Author Accepted Manuscript arising from this submission, in accordance
881 with the grant’s open access conditions.

882

883 Credit Authorship

884 **Peter Adesina:** Conceptualization, Methodology, Investigation, Formal analysis,
885 Visualization, Writing – original draft.

886 **Antoine Wautier:** Funding acquisition, Conceptualization, Methodology, Supervision,
887 Writing – review & editing.

888 Nadia Benahmed: Conceptualization, Methodology, Supervision, Writing – review & editing.

890

891 Compliance with ethical standards

892 **Conflict of interest** The authors declare that there is no conflict of interests regarding the
893 publication of this article. Publication has been approved by all authors. None of the material
894 presented in the paper is submitted or published elsewhere.

895

896

897 References

898

899 Adesina, P. A., O'Sullivan, C., & Wang, T. (2023). DEM study on the effect of particle shape
900 on the shear behaviour of granular materials. *Computational Particle Mechanics*,
901 *Accepted*.

902 Adesina, P., O'Sullivan, C., Morimoto, T., & Otsubo, M. (2022). Determining a
903 representative element volume for DEM simulations of samples with non-circular
904 particles. *Particuology*, **68**, 29–43. https://doi.org/10.1016/j.partic.2021.10.007

905 Adesina, P., O'Sullivan, C., & Wang, T. (2024). Understanding the interplay between particle
906 shape, grading and sample density on the behaviour of granular assemblies: A DEM
907 approach. *Granular Matter*, **26**(1). https://doi.org/10.1007/s10035-023-01383-2

908 Bagi, K. (1999). Microstructural stress tensor of granular assemblies with volume forces.
909 *Journal of Applied Mechanics, Transactions ASME*, **66**(4), 934–936.
910 https://doi.org/10.1115/1.2791800

911 Benahmed, N., Nguyen, T. K., Hicher, P. Y., & Nicolas, M. (2015). An experimental
912 investigation into the effects of low plastic fines content on the behaviour of sand/silt
913 mixtures. *European Journal of Environmental and Civil Engineering*, **19**(1), 109–128.
914 https://doi.org/10.1080/19648189.2014.939304

915 Bolton, M. D. (1986). The strength and dilatancy of sands. *Geotechnique*, **36**(1), 65–78.

916 Chang, C. S., & Yin, Z. Y. (2011). Micromechanical modeling for behavior of silty sand with
917 influence of fine content. *International Journal of Solids and Structures*, 48(19), 2655-
918 2667. https://doi.org/10.1016/j.ijsolstr.2011.05.014

919 Cubrinovski, M., & Ishihara, K. (2002). Maximum and minimum void ratio characteristics of
920 sands. *Soils and Foundations*, 42(6), 65–78. https://doi.org/10.3208/sandf.42.6_65

921 Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular
922 assemblies. *Geotechnique*, 30(3), 331–336. https://doi.org/10.1680/geot.1980.30.3.331

923 da Cruz, F., Emam, S., Prochnow, M., Roux, J., & Chevoir, F. (2005). Rheophysics of dense
924 granular materials : Discrete simulation of plane shear flows. *PHYSICAL REVIEW*,
925 72(February), 1–17. https://doi.org/10.1103/PhysRevE.72.021309

926 Deng, N., Wautier, A., Tordesillas, A., Thiery, Y., Yin, Z. Y., Hicher, P. Y., & Nicot, F.
927 (2022). Lifespan dynamics of cluster conformations in stationary regimes in granular
928 materials. *Physical Review E*, 105(1), 1–15.
929 https://doi.org/10.1103/PhysRevE.105.014902

930 Evans, D. M., & Zhou, S. (1995). Liquefaction Behaviour of Sand-Gravel Composites.
931 *Journal of Geotechnical and Geoenvironmental Engineering*, 121(3), 287–298.

932 Garner, S. J., & Fannin, R. J. (2010). Understanding internal erosion: a decade of research
933 following a sinkhole event. *International Journal on Hydropower and Dams*, 17(3), 93–
934 98.

935 Holubec, I., & D'Appolonia, E. (1973). Effect of particle shape on the Engineering Properties
936 of Granular Soils. *Evaluation of Relative Density and Its Role in Geotechnical Projects*
937 Involving Cohesionless Soils*, ASTM STP 523, 304–318.
938 https://doi.org/10.1051/epjconf/201714006006

939 ICOLD. (2013). *Internal Erosion of Existing Dams, Levees and Dykes, and Their
940 Foundations*. In: Bridle, R. and Fell, R., Eds., *Bulletin 164, Volume 1: Internal Erosion
941 Processes and Engineering Assessment, International Commission on Large Dams*,
942 Paris.

943 ICOLD, I. C. on L. dams. (2017). *ICOLD Bulletin 164 on Internal erosion of existing dams*,
944 *leeves and dikes, and their foundations*.

945 Jiang, M. D., Yang, Z. X., Barreto, D., & Xie, Y. H. (2018). The influence of particle-size
946 distribution on critical state behavior of spherical and non-spherical particle assemblies.
947 *Granular Matter*, 20(4), 1–15. https://doi.org/10.1007/s10035-018-0850-x

948 Kuerbis, R. H. (1989). Effect of gradation and fines content on the undrained response of
949 sand. In *Department of Civil Engineering, The University of British Columbia*.
950 https://doi.org/10.1016/0148-9062(90)90054-6

951 Lade, P. V., Liggio, C. D., & Yamamuro, J. A. (1998). Effects of Non-Plastic Fines on
952 Minimum and Maximum Void Ratios of Sand. *Geotechnical Testing Journal*, 21(4),
953 336–347. https://doi.org/10.1520/gtj11373j

954 Lade, P. V., & Yamamuro, J. A. (1997). Effects of nonplastic fines on static liquefaction of
955 sands. *Canadian Geotechnical Journal*, 34(6), 918–928. https://doi.org/10.1139/t97-052

956 Li, W., Chu, Y., Deng, G., Cai, H., Xie, D., & Lee Lee, M. (2023). Study of shear induced
957 stress redistribution in gap-graded soils by discrete element method. *Computers and
958 Geotechnics*, 156(June 2022), 105248. https://doi.org/10.1016/j.compgeo.2023.105248

959 Li, Y., Otsubo, M., Ghaemi, A., Dutta, T. T., & Kuwano, R. (2022). Transition of gap-graded
960 soil fabric – shear wave measurements and dispersion relation. *Soils and Foundations*,
961 62(1), 101092. https://doi.org/10.1016/j.sandf.2021.101092

962 Liu, D., Morimoto, T., Carraro, J. A. H., & O'Sullivan, C. (2022). A semi-empirical re-
963 evaluation of the influence of state on elastic stiffness in granular materials. *Granular
964 Matter*, 24(2). https://doi.org/10.1007/s10035-022-01215-9

965 Minh, N. H., Cheng, Y. P., & Thornton, C. (2014). Strong force networks in granular
966 mixtures. *Granular Matter*, 16(1), 69–78. https://doi.org/10.1007/s10035-013-0455-3
917 influence of fine content. *International Journal of Solids and Structures*, 48(19), 2655-
918 2667. https://doi.org/10.1016/j.ijsolstr.2011.05.014

919 Cubrinovski, M., & Ishihara, K. (2002). Maximum and minimum void ratio characteristics of
920 sands. *Soils and Foundations*, 42(6), 65–78. https://doi.org/10.3208/sandf.42.6_65

921 Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular
922 assemblies. *Geotechnique*, 30(3), 331–336. https://doi.org/10.1680/geot.1980.30.3.331

923 da Cruz, F., Emam, S., Prochnow, M., Roux, J., & Chevoir, F. (2005). Rheophysics of dense
924 granular materials : Discrete simulation of plane shear flows. *PHYSICAL REVIEW*,
925 72(February), 1–17. https://doi.org/10.1103/PhysRevE.72.021309

926 Deng, N., Wautier, A., Tordesillas, A., Thiery, Y., Yin, Z. Y., Hicher, P. Y., & Nicot, F.
927 (2022). Lifespan dynamics of cluster conformations in stationary regimes in granular
928 materials. *Physical Review E*, 105(1), 1–15.
929 https://doi.org/10.1103/PhysRevE.105.014902

930 Evans, D. M., & Zhou, S. (1995). Liquefaction Behaviour of Sand-Gravel Composites.
931 *Journal of Geotechnical and Geoenvironmental Engineering*, 121(3), 287–298.

932 Garner, S. J., & Fannin, R. J. (2010). Understanding internal erosion: a decade of research
933 following a sinkhole event. *International Journal on Hydropower and Dams*, 17(3), 93–
934 98.

935 Holubec, I., & D'Appolonia, E. (1973). Effect of particle shape on the Engineering Properties
936 of Granular Soils. *Evaluation of Relative Density and Its Role in Geotechnical Projects*
937 Involving Cohesionless Soils*, ASTM STP 523, 304–318.
938 https://doi.org/10.1051/epjconf/201714006006

939 ICOLD. (2013). *Internal Erosion of Existing Dams, Levees and Dykes, and Their
940 Foundations*. In: Bridle, R. and Fell, R., Eds., *Bulletin 164, Volume 1: Internal Erosion
941 Processes and Engineering Assessment, International Commission on Large Dams*,
942 Paris.

943 ICOLD, I. C. on L. dams. (2017). *ICOLD Bulletin 164 on Internal erosion of existing dams*,
944 *leeves and dikes, and their foundations*.

945 Jiang, M. D., Yang, Z. X., Barreto, D., & Xie, Y. H. (2018). The influence of particle-size
946 distribution on critical state behavior of spherical and non-spherical particle assemblies.
947 *Granular Matter*, 20(4), 1–15. https://doi.org/10.1007/s10035-018-0850-x

948 Kuerbis, R. H. (1989). Effect of gradation and fines content on the undrained response of
949 sand. In *Department of Civil Engineering, The University of British Columbia*.
950 https://doi.org/10.1016/0148-9062(90)90054-6

951 Lade, P. V., Liggio, C. D., & Yamamuro, J. A. (1998). Effects of Non-Plastic Fines on
952 Minimum and Maximum Void Ratios of Sand. *Geotechnical Testing Journal*, 21(4),
953 336–347. https://doi.org/10.1520/gtj11373j

954 Lade, P. V., & Yamamuro, J. A. (1997). Effects of nonplastic fines on static liquefaction of
955 sands. *Canadian Geotechnical Journal*, 34(6), 918–928. https://doi.org/10.1139/t97-052

956 Li, W., Chu, Y., Deng, G., Cai, H., Xie, D., & Lee Lee, M. (2023). Study of shear induced
957 stress redistribution in gap-graded soils by discrete element method. *Computers and
958 Geotechnics*, 156(June 2022), 105248. https://doi.org/10.1016/j.compgeo.2023.105248

959 Li, Y., Otsubo, M., Ghaemi, A., Dutta, T. T., & Kuwano, R. (2022). Transition of gap-graded
960 soil fabric – shear wave measurements and dispersion relation. *Soils and Foundations*,
961 62(1), 101092. https://doi.org/10.1016/j.sandf.2021.101092

962 Liu, D., Morimoto, T., Carraro, J. A. H., & O'Sullivan, C. (2022). A semi-empirical re-
963 evaluation of the influence of state on elastic stiffness in granular materials. *Granular
964 Matter*, 24(2). https://doi.org/10.1007/s10035-022-01215-9

965 Minh, N. H., Cheng, Y. P., & Thornton, C. (2014). Strong force networks in granular
966 mixtures. *Granular Matter*, 16(1), 69–78. https://doi.org/10.1007/s10035-013-0455-3
Ng, T.-T. (2004). Triaxial Test Simulations with Discrete Element Method and Hydrostatic Boundaries. *Journal of Engineering Mechanics*, 130(10), 1188-1194.
https://doi.org/10.1061/(asce)0733-9399(2004)130:10(1188)

Ni, Q., Tan, T. S., Dasari, G. R., & Hight, D. W. (2004). Contribution of fines to the compressive strength of mixed soils. *Geotechnique*, 54(9), 561-569.
https://doi.org/10.1680/geot.2004.54.9.561

Nicot, F., Hadda, N., Guessasma, M., Fortin, J., & Millet, O. (2013). On the definition of the stress tensor in granular media. *International Journal of Solids and Structures*, 50(14-15), 2508-2517. https://doi.org/10.1016/j.ijsolstr.2013.04.001

Otsubo, M. (2016). *Particle Scale Analysis of Soil Stiffness and Elatic Wave Propagation*.

Peters, J. F., Muthuswamy, M., Wibowo, J., & Tordesillas, A. (2005). Characterization of force chains in granular material. *Physical Review E - Statistical, Nonlinear, and Soft Matter Physics*, 72(4), 1-8. https://doi.org/10.1103/PhysRevE.72.041307

Pitman, T. D., Robertson, P. K., & Sego, D. C. (1994). Influence of fines on the collapse of loose sands. *Canadian Geotechnical Journal*, 31(5), 728-739.
https://doi.org/10.1139/t94-084

Rahman, M. M., Lo, S. R., & Gnanendran, C. T. (2008). On equivalent granular void ratio and steady state behaviour of loose sand with fines. *Canadian Geotechnical Journal*, 45(10), 1439-1456. https://doi.org/10.1139/T08-064

Salgado, R., Bandini, P., & Karim, A. (2000). Shear Strength and Stifness of Silty Sand. *Journal of Geotechnical and Geoenvironmental Engineering*, 126(5), 451-462.

Sarkar, D., Goudarzy, M., König, D., & Wichtmann, T. (2020). Influence of particle shape and size on the threshold fines content and the limit index void ratios of sands containing non-plastic fines. *Soils and Foundations*, 60(3), 621-633.
https://doi.org/10.1016/j.sandf.2020.02.006

Shire, T. (2014). *Micro-scale Modelling of Granular Filters*. Thesis submitted to the Department of Civil and Engineering Civil and Environmental Engineering, Imperial College, London, United Kingdom.

Shire, T., O'Sullivan, C., & Hanley, K. J. (2016). The influence of fines content and size-ratio on the micro-scale properties of dense bimodal materials. *Granular Matter*, 18(3).
https://doi.org/10.1007/s10035-016-0654-9

Shire, T., O'Sullivan, C., Hanley, K. J., & Fannin, R. J. (2014). Fabric and Effective Stress Distribution in Internally Unstable Soils. *Journal of Geotechnical and Geoenvironmental Engineering*, 140(12), 04014072.
https://doi.org/10.1061/(ASCE)GT.1943-5606.0001184

Skempton, A., & Brogan, J. (1994). Experiments on piping in sandygravels. *Géotechnique*, 44(3), 449-460.

Skempton, A. W., & Brogan, J. M. (1994). Experiments on piping in sandy gravels. *Geotechnique*, 44(3), 449-460. https://doi.org/10.1680/geot.1995.45.3.565

Šmilauer et al, V. (2021). *Yade Documentation 3rd ed*. The Yade Project (p. DOI:10.5281/zenodo.5705394 (http://yade-dem.org/do).

Sufian, A., Artigaut, M., Shire, T., & O'Sullivan, C. (2021). Influence of Fabric on Stress Distribution in Gap-Graded Soil. *Journal of Geotechnical and Geoenvironmental Engineering*, 147(5), 1-14. https://doi.org/10.1061/(asce)gt.1943-5606.0002487

Thevanayagam, S. (1998). Effect of Fines and Confining Stress on Undrained Shear Strength of Silty Sands. *Journal of Geotechnical and Geoenvironmental Engineering*, 124(6), 479-491. https://doi.org/10.1061/(asce)1090-0241(1999)125:11(1024)

Thevanayagam, S., & Mohan, S. (2000). Intergranular state variables and stress-strain behaviour of silty sands. *Geotechnique*, 50(1), 1-23.
https://doi.org/10.1680/geot.2000.50.1.1
967 Ng, T.-T. (2004). Triaxial Test Simulations with Discrete Element Method and Hydrostatic

968 Boundaries. *Journal of Engineering Mechanics*, **130**(10), 1188-1194.

969 https://doi.org/10.1061/(asce)0733-9399(2004)130:10(1188)

970 Ni, Q., Tan, T. S., Dasari, G. R., & Hight, D. W. (2004). Contribution of fines to the
971 compressive strength of mixed soils. *Geotechnique*, **54**(9), 561–569.

972 https://doi.org/10.1680/geot.2004.54.9.561

973 Nicot, F., Hadda, N., Guessasma, M., Fortin, J., & Millet, O. (2013). On the definition of the
974 stress tensor in granular media. *International Journal of Solids and Structures*, **50**(14–
975 15), 2508–2517. https://doi.org/10.1016/j.ijsolstr.2013.04.001

976 Otsubo, M. (2016). *Particle Scale Analysis of Soil Stiffness and Elatic Wave Propagation*.

977 Peters, J. F., Muthuswamy, M., Wibowo, J., & Tordesillas, A. (2005). Characterization of
978 force chains in granular material. *Physical Review E - Statistical, Nonlinear, and Soft
979 Matter Physics*, **72**(4), 1–8. https://doi.org/10.1103/PhysRevE.72.041307

980 Pitman, T. D., Robertson, P. K., & Sego, D. C. (1994). Influence of fines on the collapse of
981 loose sands. *Canadian Geotechnical Journal*, **31**(5), 728–739.

982 https://doi.org/10.1139/t94-084

983 Rahman, M. M., Lo, S. R., & Gnanendran, C. T. (2008). On equivalent granular void ratio
984 and steady state behaviour of loose sand with fines. *Canadian Geotechnical Journal*,
985 45(10), 1439–1456. https://doi.org/10.1139/T08-064

986 Salgado, R., Bandini, P., & Karim, A. (2000). Shear Strength and Stifness of Silty Sand.
987 *Journal of Geotechnical and Geoenvironmental Engineering*, **126**(5), 451–462.

988 Sarkar, D., Goudarzy, M., König, D., & Wichtmann, T. (2020). Influence of particle shape
989 and size on the threshold fines content and the limit index void ratios of sands
990 containing non-plastic fines. *Soils and Foundations*, **60**(3), 621–633.

991 https://doi.org/10.1016/j.sandf.2020.02.006

992 Shire, T. (2014). *Micro-scale Modelling of Granular Filters*. Thesis submitted to the
993 Department of Civil and Engineering Civil and Environmental Engineering, Imperial
994 College, London, United Kingdom.

995 Shire, T., O'Sullivan, C., & Hanley, K. J. (2016). The influence of fines content and size-
996 ratio on the micro-scale properties of dense bimodal materials. *Granular Matter*, **18**(3).
997 https://doi.org/10.1007/s10035-016-0654-9

998 Shire, T., O'Sullivan, C., Hanley, K. J., & Fannin, R. J. (2014). Fabric and Effective Stress
999 Distribution in Internally Unstable Soils. *Journal of Geotechnical and
1000 Geoenvironmental Engineering*, **140**(12), 04014072.

1001 https://doi.org/10.1061/(ASCE)GT.1943-5606.0001184

1002 Skempton, A., & Brogan, J. (1994). Experiments on piping in sandygravels. *Géotechnique*, **44**(3), 449–460.

1003 Skempton, A. W., & Brogan, J. M. (1994). Experiments on piping in sandy gravels. *Geotechnique*, **44**(3), 449–460. https://doi.org/10.1680/geot.1995.45.3.565

1004 Šmilauer et al, V. (2021). *Yade Documentation 3rd ed. The Yade Project* (p.
1005 DOI:10.5281/zenodo.5705394 (http://yade-dem.org/do).

1006 Sufian, A., Artigaut, M., Shire, T., & O'Sullivan, C. (2021). Influence of Fabric on Stress
1009 Distribution in Gap-Graded Soil. *Journal of Geotechnical and Geoenvironmental
1010 Engineering*, **147**(5), 1–14. https://doi.org/10.1061/(asce)gt.1943-5606.0002487

1011 Thevanayagam, S. (1998). Effect of Fines and Confining Stress on Undrained Shear Strength
1012 of Silty Sands. *Journal of Geotechnical and Geoenvironmental Engineering*, **124**(6),
1013 479-491. https://doi.org/10.1061/(asce)1090-0241(1999)125:11(1024)

1014 Thevanayagam, S., & Mohan, S. (2000). Intergranular state variables and stress-strain
1015 behaviour of silty sands. *Geotechnique*, **50**(1), 1–23.

1016 https://doi.org/10.1680/geot.2000.50.1.1
1017 Thevanayagam, S., Shenthan, T., Mohan, S., & Liang, J. (2002). Undrained Fragility of Clean Sands, Silty Sands, and Sandy Silts. *Journal of Geotechnical and Geoenvironmental Engineering*, 128(10), 849-859. https://doi.org/10.1061/(asce)1090-0241(2002)128:10(849)

1021 Thornton, C. (2000). Numerical simulations of deviatoric shear deformation of granular media. *Geotechnique*, 50(1), 43-53. https://doi.org/10.1680/geot.2000.50.1.43

1023 Thornton, C., & Anthony, S. J. (1998). Quasi-static deformation of particulate media. *Phil. Trans. R. Soc. Lond. A*, 356, 2763-2782.

1025 Thornton, Colin. (2015). Granular Dynamics, Contact Mechanics and Particle System Simulations: A DEM study. *Granular Dynamics, Contact Mechanics and Particle System Simulations: A DEM Study*, 24, 1-195. https://doi.org/10.1007/978-3-319-18711-2

1029 Tordesillas, A., Walker, D. M., & Lin, Q. (2010). Force cycles and force chains. *Physical Review E - Statistical, Nonlinear, and Soft Matter Physics*, 81(1), 5-7.
https://doi.org/10.1103/PhysRevE.81.011302

1032 Vaid, Y. P. (1994). Liquefaction of silty soils. *In: Proceedings of Ground Failures Under Seismic Conditions. Geotech. Spec. Publ.* 44, ASCE, New York, 1-16.

1034 Vaid, Y. P., & Sasitharan, S. (1992). The strength and dilatancy of sand. *Canadian Geotechnical Journal*, 29, 522-526. https://doi.org/10.1680/geot.1987.37.4.517

1036 Vallejo, L. E. (2001). Interpretation of the limits in shear strength in binary granular mixtures. *Canadian Geotechnical Journal*, 38(5), 1097-1104.
https://doi.org/10.1139/cgj-38-5-1097

1039 Wang, T., Liu, S., Wautier, A., & Nicot, F. (2022). Updated skeleton void ratio for gravelly sand mixtures considering the effect of grain size distribution. *Canadian Geotechnical Journal*, 59(11), 12-23. https://doi.org/10.1139/cgj-2021-0443

1042 Wautier, A., Bonelli, S., & Nicot, F. (2018). Micro-inertia origin of instabilities in granular materials. *International Journal for Numerical and Analytical Methods in Geomechanics*, 42(9), 1037-1056. https://doi.org/10.1002/nag.2777

1045 Xiao, Y., Xiang, J., Liu, H., & Ma, Q. (2017). Strength-dilatancy relation of sand containing non-plastic fines. *Geotechnique Letters*, 7(2), 204-210.
https://doi.org/10.1680/jgele.16.00144

1048 Yang, S., Lacasse, S., & Sandven, R. (2005). Determination of the Transitional Fines Content of Mixtures of Sand and Non-plastic Fines. *Geotechnical Testing Journal*, 29(2), 102-107. https://doi.org/10.1520/GTJ14010

1051 Yin, Z. Y., Zhao, J., & Hicher, P. Y. (2014). A micromechanics-based model for sand-silt mixtures. *International Journal of Solids and Structures*, 51(6), 1350-1363.
https://doi.org/10.1016/j.ijsolstr.2013.12.027

1054 Zhu, H., Nguyen, H. N. G., Nicot, F., & Darve, F. (2016). On a common critical state in localized and diffuse failure modes. *Journal of the Mechanics and Physics of Solids*, 95, 112-131. https://doi.org/10.1016/j.jmps.2016.05.026

1057 Zhu, H., Nicot, F., & Darve, F. (2016). Meso-structure organization in two-dimensional granular materials along biaxial loading path. *International Journal of Solids and Structures*, 96, 25-37. https://doi.org/10.1016/j.ijsolstr.2016.06.025

1060 Zuo, L., & Baudet, B. A. (2015). Determination of the transitional fines content of sand-non plastic fines mixtures. *Soils and Foundations*, 55(1), 213-219.
https://doi.org/10.1016/j.sandf.2014.12.017

1063
1064
1065
1017 Thevanayagam, S., Shenthan, T., Mohan, S., & Liang, J. (2002). Undrained Fragility of Clean Sands, Silty Sands, and Sandy Silts. *Journal of Geotechnical and Geoenvironmental Engineering*, 128(10), 849-859. https://doi.org/10.1061/(asce)1090-0241(2002)128:10(849)
1018
1019
1020
1021 Thornton, C. (2000). Numerical simulations of deviatoric shear deformation of granular media. *Geotechnique*, 50(1), 43-53. https://doi.org/10.1680/geot.2000.50.1.43
1022
1023 Thornton, C., & Anthony, S. J. (1998). Quasi-static deformation of particulate media. *Phil. Trans. R. Soc. Lond. A*, 356, 2763-2782.
1024
1025 Thornton, Colin. (2015). Granular Dynamics, Contact Mechanics and Particle System Simulations: A DEM study. *Granular Dynamics, Contact Mechanics and Particle System Simulations: A DEM Study*, 24, 1-195. https://doi.org/10.1007/978-3-319-18711-2
1026
1027
1028
1029 Tordesillas, A., Walker, D. M., & Lin, Q. (2010). Force cycles and force chains. *Physical Review E - Statistical, Nonlinear, and Soft Matter Physics*, 81(1), 5-7.
1030
1031 https://doi.org/10.1103/PhysRevE.81.011302
1032 Vaid, Y. P. (1994). Liquefaction of silty soils. *In: Proceedings of Ground Failures Under Seismic Conditions. Geotech. Spec. Publ.* 44, ASCE, New York, 1-16.
1033
1034 Vaid, Y. P., & Sasitharan, S. (1992). The strength and dilatancy of sand. *Canadian Geotechnical Journal*, 29, 522-526. https://doi.org/10.1680/geot.1987.37.4.517
1035
1036 Vallejo, L. E. (2001). Interpretation of the limits in shear strength in binary granular mixtures. *Canadian Geotechnical Journal*, 38(5), 1097-1104.
1037
1038 https://doi.org/10.1139/cgj-38-5-1097
1039 Wang, T., Liu, S., Wautier, A., & Nicot, F. (2022). Updated skeleton void ratio for gravelly sand mixtures considering the effect of grain size distribution. *Canadian Geotechnical Journal*, 59(11), 12-23. https://doi.org/10.1139/cgj-2021-0443
1040
1041
1042 Wautier, A., Bonelli, S., & Nicot, F. (2018). Micro-inertia origin of instabilities in granular materials. *International Journal for Numerical and Analytical Methods in Geomechanics*, 42(9), 1037-1056. https://doi.org/10.1002/nag.2777
1043
1044
1045 Xiao, Y., Xiang, J., Liu, H., & Ma, Q. (2017). Strength-dilatancy relation of sand containing non-plastic fines. *Geotechnique Letters*, 7(2), 204-210.
1046
1047 https://doi.org/10.1680/jgele.16.00144
1048 Yang, S., Lacasse, S., & Sandven, R. (2005). Determination of the Transitional Fines Content of Mixtures of Sand and Non-plastic Fines. *Geotechnical Testing Journal*, 29(2), 102-107. https://doi.org/10.1520/GTJ14010
1049
1050
1051 Yin, Z. Y., Zhao, J., & Hicher, P. Y. (2014). A micromechanics-based model for sand-silt mixtures. *International Journal of Solids and Structures*, 51(6), 1350-1363.
1052
1053 https://doi.org/10.1016/j.ijsolstr.2013.12.027
1054 Zhu, H., Nguyen, H. N. G., Nicot, F., & Darve, F. (2016). On a common critical state in localized and diffuse failure modes. *Journal of the Mechanics and Physics of Solids*, 95, 112-131. https://doi.org/10.1016/j.jmps.2016.05.026
1055
1056
1057 Zhu, H., Nicot, F., & Darve, F. (2016). Meso-structure organization in two-dimensional granular materials along biaxial loading path. *International Journal of Solids and Structures*, 96, 25-37. https://doi.org/10.1016/j.ijsolstr.2016.06.025
1058
1059
1060 Zuo, L., & Baudet, B. A. (2015). Determination of the transitional fines content of sand-non plastic fines mixtures. *Soils and Foundations*, 55(1), 213-219.
1061
1062 https://doi.org/10.1016/j.sandf.2014.12.017
1063
1064
1065