# The mechanical behaviour of gap-graded soils under triaxial compression and extension loading paths

Peter Adesina, Antoine Wautier, Nadia Benahmed

► To cite this version:

Peter Adesina, Antoine Wautier, Nadia Benahmed. The mechanical behaviour of gap-graded soils under triaxial compression and extension loading paths. Acta Geotechnica, 2025, 20 (8), pp.3811-3834. 10.1007/s11440-025-02600-3. hal-05228338

HAL Id: hal-05228338

https://hal.inrae.fr/hal-05228338v1

Submitted on 7 Nov 2025

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
THE MECHANICAL BEHAVIOUR OF GAP-GRADED SOILS UNDER TRIAXIAL
COMPRESSION AND EXTENSION LOADING PATHS

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

¹,²INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182

Aix-en-Provence, France

³Department of Civil and Mineral Engineering, University of Toronto, 35 St. George St., Toronto, ON, Canada, M5S 1A4.

*Corresponding authors. E-mail addresses: peter.adesina@inrae.fr, peter.adesina@utoronto.ca
(P. Adesina); antoine.wautier@inrae.fr (Antoine Wautier); nadia.benahmed@inrae.fr (Nadia Benahmed)

**Abstract**

The influence of mechanical loading paths on the characteristics of gap-graded granular assemblies was investigated using the discrete element method (DEM). Dense and loose gap-graded assemblies with finer fraction content, $f_c$, ranging from 0-100% were prepared and subjected to drained triaxial compression and extension loading paths. After examining key macroscale quantities, micromechanical analyses were conducted to elicit the particle-scale characteristics including the evolution of the fabric of the assemblies under the different loading paths. The results of the DEM analysis confirm the validity of the Mohr-Columb failure criteria at the critical state. While the mobilised friction angle at the peak is higher under extension than in compression, no significant difference was obtained in the critical state friction angle for both loading paths. Despite the higher mean stress transmitted by the gap-graded assemblies under compression in comparison to extension, the contribution of the finer particles to the total mean stress is not significantly influenced by the loading paths. Our data show that the variation in the fabric of granular assemblies under different loading paths does not always stem from an initial inherent anisotropy. Fabric anisotropy is marginally higher under extension than in compression despite having an initial isotropic fabric.

**Keywords: Gap-graded assemblies, triaxial compression/extension, macros micromechanical behaviour, DEM**

<footnote n="1">¹Post doctoral researcher, INRAE, Aix-Marseille University; UMR RECOVER, 3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence, France (Now Post doctoral Fellow, Department of Civil and Mineral Engineering, University of Toronto, 35 St. George St., Toronto, ON, Canada M5S 1A4). Email: peter.adesina@inrae.fr; peter.adesina@utoronto.ca; 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; nadia.benahmed@inrae.fr</footnote>
32

# 1. Introduction

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

Mechanical loading path is one of the important extrinsic factors that dictate the engineering behaviour of soils. In the past, several attempts have been made to elicit the mechanical behaviour of soils under different loading paths including drained and undrained monotonic loading, cyclic loading, and varying intermediate stress conditions. In their experimental study, Salvatore et al. [1] identified differences in the behaviour of uniformly graded Houston sand when subjected to triaxial compression and extension loading paths and concluded that the classical critical state framework which was developed from triaxial compression tests is not valid for the constitutive modelling of behaviour under other loading paths. Yamamuro et al., [2] attributed the lower shear strength obtained for uniformly graded Cambria sand under triaxial extension to strain localisation and geometrical instability which progresses more readily during extension tests than in compression. When subjected to undrained triaxial extension, Toyoura sand exhibited a more contractive and softer behaviour than under triaxial compression [3]. Differences in the friction angle of Cambria sand subjected to drained triaxial compression and extension tests with high confining pressure were reported by Yamamuro and Lade [4]. These identified influence of loading paths show that predicting the behaviour of sand using the characteristics developed under a single loading path can be misleading.

49

50

51

52

53

54

55

56

57

58

Most of the earlier investigations on the influence of triaxial compression and extension loading paths involve experimental studies where narrowly graded or uniformly graded sands were considered [1–5]. Meanwhile, only a few studies have considered the influence of loading path on the mechanical behaviour of gap-graded soils, which are prevalent in nature and are used in several civil engineering applications. The mechanical behaviour of gap-graded soils which contain a mixture of finer and coarser fractions of particles differ from uniformly graded soils [6], and depends on the finer fraction content, $f_c$. In the experimental work by Kuwano et al., [7] and Sanjei et al., [8] where the influence of triaxial compression and torsional shear on the mechanical behaviour of eroded and un-eroded gap-graded soils were investigated, only underfilled assemblies with $f_c \le 20\%$ were considered.

59

60

61

62

63

In the present study, we investigate the mechanical behaviour of gap-graded assemblies with finer fraction content, $f_c \in [0\%; 100\%]$, under triaxial compression and extension loading paths, using the numerical discrete element method (DEM). Our objectives are to: (i) identify the key characteristics of gap-graded soils under compression and extension loading paths (ii) elicit the contribution of the finer fraction to the overall behaviour of gap-graded soils when
64 subjected to these loading paths (iii) determine the influence of the loading paths on the fabric
65 of gap-graded assemblies.

66 Earlier experimental work and constitutive modelling of uniformly graded sands suggest that
67 the critical state strength under triaxial compression differ from the strength under extension
68 because of the inherent anisotropy in sands [1,3,9]. In this study, we take advantage of the
69 unique capability of DEM to investigate the influence of the fabric on the shear strength of
70 gap-graded assemblies when subjected to triaxial compression and extension loading paths.
71 The details of the numerical simulation approach adopted are presented in Section 2. In Section
72 3, we show the influence of the finer fraction on the strength and dilatancy and the microscopic
73 characteristics of gap-graded assemblies subjected to the different loading paths. In Sections 4
74 and 5, we investigate the influence of the loading paths on the contributions of the finer and
75 coarser particle fractions and contact types to the total mean stress transmitted by the
76 assemblies. Finally, the contact normal fabric anisotropy and its link to the macroscale
77 characteristics of the assemblies is explored in Section 6.

78

## 2. Numerical Simulation approach

80 DEM simulations were conducted using the open source code YADE [10] with interparticle
81 interactions following a simple linear elasto-plastic law proposed by Cundal and Strack [11]
82 and illustrated in Fig. 1. In this law, a contact is created when two particles overlap. Two linear
83 springs and a slider provide the link between the two particles in contact. A normal force, $F_n$,
84 is defined as proportional to the overlap, $u_n$, between the particles, where the normal stiffness
85 constant, $k_n$, is given as $k_n = E \frac{2R_1R_2}{R_1+R_2}$, with E a material modulus, and $R_1, R_2$ the radii of the
86 particles. Also, a tangential force, $F_t$, which is proportional to the relative tangential
87 displacement, $u_t$, between the two particles is introduced, where the tangential stiffness
88 constant, $k_t$, is a fraction of the normal stiffness constant ($k_t = \alpha k_n$). While the normal
89 displacement is defined on the current geometry, the tangential displacement is defined
90 incrementally, with $\Delta u_t$ being the relative tangential displacement (translation and rotation) of
91 the contact point between two geometrical configurations. The linear behaviour of $F_t$ is
92 governed by the Columb friction condition based on an internal friction angle, $\phi$, which
93 determines the largest possible ratio, $F_n/F_t$ (see Fig. 1). After the interparticle contact forces
94 are computed, the induced particle displacements are determined by integrating Newton's
95
96

second law of motion using the DEM code YADE. The parameters for the simulations are
detailed in Table 1.

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

97
98
99

<caption>Fig. 1 Elasto-frictional contact law used in discrete element method simulations</caption>

100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118

The assemblies with finer fraction contents, $f_c = 0\%$ and $f_c = 100\%$ were generated to have a monodisperse grading which contained coarser particles with diameter $D_c=3.075$ mm, and finer particles with diameter, $D_f=0.375$ mm, respectively. According to the soil classification system in [12], $D_f=0.375$ mm is classified as sand while $D_c=3.075$ mm is classified as gravel. Therefore, in this paper, the assemblies having $f_c = 0\%$ and $f_c = 100\%$ will be referred to as clean gravel and clean sand, respectively. The sand-gravel mixtures with $f_c \in [10\%; 90\%]$ were generated to have a bimodal grading with a particle size ratio, $\lambda = D_c/D_f = 8.2$. This size ratio is sufficiently large for the finer particles to migrate through the constriction between the coarser particles if subjected to a hydraulic gradient [13–15]. We considered a bimodal grading to ensure that the results obtained emanate solely from the interactions between the coarser and finer fractions. The number of the finer and coarser particles for each $f_c$ value is presented in Table 2. All assemblies contained 100,000 spherical particles (Fig. 2). This sample size is considered a representative element volume (REV) for all finer fraction content following the parametric study presented in Adesina et al., [16] (see Appendix B for particle-scale REV analysis). Shi et al [17] specifically studied the REV for gap-graded assemblies within 3-D rigid boundaries. They found that REV size for gap-graded assemblies decreases with an increase in the finer fraction content, $f_c$, and recommended that the sample size ratio, SSR = $D/d_{max} = 8-3$, for $f_c = 10\% - 85\%; \lambda = 8$, where D is the initial sample width and $d_{max}$ is the maximum particle diameter. The SSR values for the assemblies in the present study
119

($\lambda$ = 8.2) (see Table 2), are higher (by 38%-46%) than the values recommended by Shi et al [17]. This ensured that the assemblies were large enough to minimise the influence of any potential strain and stress non-uniformity close to the rigid walls, on the simulation results.

120

121

122

The different samples considered in this study were first subjected to isotropic compression at 100 kPa within six frictionless walls. During this step, the interparticle friction coefficient, $\mu$, was controlled to prepare assemblies of varying initial densities. Specifically, $\mu$=0.03 was used to generate dense assemblies, while $\mu$=0.5 was used for loose assemblies. The void ratios for the dense and the loose assemblies are denoted as $e_{min}$ and $e_{max}$ respectively. It is important to note that these void ratio values cannot be directly mapped to the values obtained using standard experimental procedures; rather, they follow the standard practice adopted in DEM to generate assemblies of varying initial densities. Isotropic compression was terminated upon reaching the target stress and upon attaining the unbalanced force ratio<fnref n="2" /> limit of 0.01 imposed to ensure static equilibrium.

123

124

125

126

127

128

129

130

131

132

Following Jiang et al. [18], density scaling was employed in our simulations to increase the critical timestep, $\Delta t_{cr} \propto r\sqrt{\rho/E}$, where r is the radius of the smallest particle in the granular system, $\rho$ is the particle density and E is the Young's modulus. While particle density was scaled by 1000 times the original value for the assemblies having $f_c \in [10\% ; 100\%]$, for the assemblies having $f_c = 0\%$ which contained only the coarser particles, particle density was scaled by 10 times the original value. This was such that $\Delta t_{cr}$ values for all assemblies were within the same order of magnitude ($\Delta t_{cr} = 3.27 \times 10^{-5}$ s for $f_c \in [10\% ; 100\%]$; $\Delta t_{cr} = 2.68 \times 10^{-5}$ s for $f_c = 0\%$). A fraction of this critical timestep values ($0.9\Delta t_{cr}$) was applied to ensure numerical stability. It has been shown in earlier studies that applying density scaling enhances computational efficiency without compromising the quasi-static conditions of DEM simulations, provided the inertia number of the simulations remains sufficiently small [19,20].

133

134

135

136

137

138

139

140

141

142

143

Following isotropic compression, each assembly was subjected to drained triaxial compression in simulation series 1, and drained triaxial extension in simulation series 2. For both series of simulations, a constant confining pressure of $\sigma_{xx}=\sigma_{yy}=100$ kPa was applied in both the x-direction and y-direction. A strain rate of 0.01 s⁻¹ and -0.01 s⁻¹ was applied in the z-direction for the triaxial compression and extension loading paths, respectively<fnref n="3" />. A numerical damping

144

145

146

147

<footnote n="2">² The maximum resultant force on the different grains divided by the mean contact force.</footnote>

<footnote n="3">³ Throughout this paper, soil mechanics conventions are used with compression being positive and extension being negative.</footnote>
coefficient of 0.05 is used for both simulation series. The simulations conditions imposed in this study were carefully chosen to ensure that the simulations were quasi-static, with inertia number, $I = \dot{\epsilon}d\sqrt{\frac{\rho}{p}} < 10^{-3}$, where $\dot{\epsilon}$ is the shear rate, $d$ is the particle diameter, $\rho$ is the particle density and $p$ is the confining pressure [21,22]. For our simulations, $I_{dmax} = 1.60 \times 10^{-4}$ and $I_{dmin} = 1.95 \times 10^{-5}$ for the coarser and the finer grains, respectively. During shearing, $\mu=0.5$ was used for all simulations, regardless of the $\mu$ value used during the sample preparation stage$^4$. A total of 60 shearing simulations were conducted for this study. Fifteen fines contents, ranging from 0% to 100%, were considered for two relative densities (dense and loose) and two loading paths (compression and extension).

<table><caption>**Table 1** Simulation parameters</caption><thead><tr><td>Parameter</td><td>Value</td></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>, for f<sub>c</sub> = 0%</td><td>2700 × 10<sup>1</sup></td></tr><tr><td>Particle density, ρ, kg/m<sup>3</sup> for f<sub>c</sub> ∈ [10%; 100%]</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>

<footnote n="4">⁴ Given that the assemblies in this study are made of spherical particles and no rolling resistance was applied in our simulations, $\mu = 0.5$ was used during shearing to achieve macroscopic friction angles similar to the values obtained for sands. However, we note that the $\mu$ value used in this study is higher than the range ($\mu = 0.12 - 0.35$) determined experimentally for quartz sand particles [72]. Also, it has been shown in DEM studies (e.g. Huang et al [73]) that particle interaction may be dominated by particle rolling for $\mu \ge 0.5$, as opposed to sliding for $\mu < 0.5$.</footnote>
<table><caption>171 Table 2 Number of particles for each assembly</caption><thead><tr><th>Percentage of fine particles by weight</th><th>No. of coarse particles</th><th>No. of fine particles</th><th>Total</th><th>Sample size ratio</th><th>Percentage of fine particles by weight</th><th>No. of coarse particles</th><th>No. of fine particles</th><th>Total</th><th>Sample size ratio</th></tr></thead><tbody><tr><td>0%</td><td>100,000</td><td>-</td><td></td><td>44.0</td><td>45%</td><td>222</td><td>99,778</td><td></td><td>6.6</td></tr><tr><td>10%</td><td>1,612</td><td>98,388</td><td></td><td>11.3</td><td>50%</td><td>182</td><td>99,818</td><td></td><td>6.4</td></tr><tr><td>15%</td><td>1,021</td><td>98,979</td><td></td><td>9.7</td><td>60%</td><td>121</td><td>99,879</td><td></td><td>6.1</td></tr><tr><td>20%</td><td>723</td><td>99,277</td><td>100,000</td><td>8.7</td><td>70%</td><td>78</td><td>99,922</td><td>100,000</td><td>5.8</td></tr><tr><td>25%</td><td>543</td><td>99,457</td><td></td><td>7.9</td><td>80%</td><td>45</td><td>99,955</td><td></td><td>5.6</td></tr><tr><td>30%</td><td>423</td><td>99,577</td><td></td><td>7.4</td><td>90%</td><td>20</td><td>99,980</td><td></td><td>5.5</td></tr><tr><td>35%</td><td>337</td><td>99,663</td><td></td><td>7.1</td><td>100%</td><td>-</td><td>100,000</td><td></td><td>5.4</td></tr><tr><td>40%</td><td>272</td><td>99,728</td><td></td><td>6.8</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

172

173

<figure><img src="image_3.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(f)</figcaption></figure>

174

175

176

177

178

179

180

<caption>Fig. 2 Dense specimens with different percentages of fines content: (a) $f_c$ = 0 % (b) $f_c$ = 10 % (c) $f_c$ = 30 % (d) $f_c$ = 50 % (e) $f_c$ = 70 % (f) $f_c$ = 100 %. (The particle sizes in each assembly are not drawn to the same scale. Sample size decreases as $f_c$ increases so that all assemblies have the same total number of particles).</caption>
### 3. Macroscale and micromechanical characteristics of the binary mixtures with varying fines content under drained triaxial compression and extension.

#### 3.1 Macroscale characteristics

Fig. 3 shows representative loading paths for dense and loose assemblies having $f_c$ = 30%, 50% and 70%, under drained triaxial compression and extension. As expected, the slopes of both loading paths in the $q-p$ plane, $\Delta q/\Delta p$ = 3 since the shearing is triaxial, where the deviatoric stress $q = \sigma_{zz} - (\sigma_{yy} + \sigma_{xx})/2$ under compression, $q = (\sigma_{yy} + \sigma_{xx})/2 - \sigma_{zz}$ under extension, and the mean stress, $p = (\sigma_{zz} + \sigma_{yy} + \sigma_{xx})/3$ (axisymmetric conditions). The loading paths are similar for all $f_c$ values. The maximum point on the loading paths for the dense assemblies marks the peak shear strength, $(q/p)_{peak}$, while the maximum point for the loose assemblies represents the critical state shear strength, $(q/p)_{crit}$. Although the evolution of the loading paths is similar for all $f_c$ values, $(q/p)_{peak}$ and $(q/p)_{crit}$ values vary with $f_c$ (see Fig. 5).

<figure><img src="image_9.png" /><figcaption>**Fig. 3** Stress-path in $q-p$ plane (a) $f_c$=30% (b) $f_c$=50% (c) $f_c$=70%. Solid gray lines mark the peak lines while dashed gray lines mark the critical state lines in both compression and extension</figcaption></figure>

Fig. 4a shows the initial void ratios for the dense and loose assemblies of varying fines contents. As expected, similar initial void ratio values were obtained for the assemblies having $f_c$ = 0 and $f_c$ = 1.0, at both the dense and loose states. In agreement with earlier experimental and numerical studies on granular mixtures [23–31], initial increase in $f_c$ decreases both $e_{min}$ and
$e_{max}$, as the finer particles fill the voids between the coarser particles. This initial reduction in void ratio continues until the threshold fines content, $f_c^{th}$, where the voids between the coarser particles have been fully filled by the finer fraction. Beyond the $f_c^{th}$, further increase in $f_c$ results in an increase in void ratio, as the finer particles separate the coarser fraction. The $f_c^{th}$ for dense assemblies ($f_c = 0.30$) is lower than for the loose assemblies ($f_c = 0.35$), which is consistent with the presence of larger pores in loose assemblies. As shown in earlier studies, $f_c^{th}$ also depends on the assembly grading and particle shape [23, 27, 32–35] which are beyond the scope of this study.

Fig. 4b shows the influence of $f_c$ on the void ratios at the peak, $e_{peak}$ for the dense assemblies, and the critical state void ratios, $e_{crit}$, for both dense and loose assemblies, under triaxial compression and extension loading paths. As expected, $e_{peak}$ is lower than $e_{crit}$ since dense assemblies continue to dilate during shearing (i.e. exhibit an increase in void ratio and volume), until reaching the critical state where both the void ratio and volume remains constant (no visible shear bands were observed in these dense assemblies). The influence of $f_c$ on both $e_{peak}$ and $e_{crit}$ follows the same trend as the initial void ratio $e_{min}$ and $e_{max}$ (Fig. 4a), although the threshold fines content shifted from $f_c^{th}=30\%$ at the peak to $f_c^{th}=35\%$ at the critical state. No significant influence of the loading paths was observed as similar $e_{peak}$ and $e_{crit}$ values were obtained at each $f_c$ for both triaxial compression and extension tests.

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

<figure><img src="image_11.png" /></figure>
**Fig. 4** Effect of fines content, $f_c$, on (a) the initial void ratios, $e$ (b) the void ratios at the peak, $e_{peak}$, and the critical state $e_{crit}$. Fitting curves correspond to the conceptual model proposed by Yin et al., [36] for sand-silt mixtures

The stress-strain and volumetric strain responses obtained under both triaxial compression and extension (see Figs. 16-19 in Appendix A) agree with established trends reported for sands in the geomechanics literature. The stress ratio $q/p$ for the initially dense assemblies increases to a peak and then softens until a plateau which corresponds to the critical state is reached. As expected, these initially dense assemblies ultimately exhibited a dilative volumetric response ($\epsilon_{vol} < 0$) under the triaxial compression and extension loading paths (Figs. 16 & 17). However, the initial contraction that precedes the dilation exhibited by the dense assemblies under triaxial compression was not observed under extension (see insets of Figs. 16d-f & Figs. 17d-f in Appendix A). The stress ratio $q/p$ for the initially loose assemblies increased monotonically from the start of shearing to the critical state without exhibiting significant softening. In general, these initially loose assemblies exhibited a contractive volumetric response ($\epsilon_{vol} > 0$) under both triaxial compression and extension (Figs 20 & 21).

Fig. 5 shows the stress ratio values at the peak, $(q/p)_{peak}$, for the dense assemblies and the critical state stress ratio, $(q/p)_{crit}$, for both dense and loose assemblies under triaxial compression and extension loading paths, for all $f_c$ values. The corresponding plots for the peak friction angle, $\phi_{peak}$, and the critical state friction angle, $\phi_{crit}$, are presented in Fig. 5c and Fig. 5d, respectively, where $\phi = \sin^{-1}[(\sigma_{zz} - \sigma_{yy})/(\sigma_{zz} + \sigma_{yy})]$ under triaxial compression and $\phi = \sin^{-1}[(\sigma_{yy} - \sigma_{zz})/(\sigma_{zz} + \sigma_{yy})]$ under extension. Given the fluctuations observed at the critical state, especially for $f_c \in [10\%; 40\%]$, $(q/p)_{crit}$, $\phi_{crit}$ and $e_{crit}$ are average values over $\epsilon_a > 0.35$. For both compression and extension loading paths, $(q/p)_{peak}$ and $\phi_{peak}$ initially increase monotonically as $f_c$ increases until $f_c^{th}$, and then decrease with further increase in $f_c$. This trend, which reflects the void ratio of the assemblies, is well established in previous experimental and numerical studies on soil mixtures [16,28]. The $(q/p)_{peak}$ and $\phi_{peak}$ values for the clean gravel and clean sand assemblies are similar. For all $f_c$ values, $(q/p)_{peak}$ is higher in compression than in extension (Fig. 5a) in agreement with the experimental study by Salvatore et al., [1] for uniformly graded sands. However, because the relation between friction and stress ratio is different in compression and extension, the reverse is true for $\phi_{peak}$ (Fig. 5c), as also shown in earlier experimental data on glass beads and sands [5,37–42]. The fact that the peak friction in compression is different than in extension indicates
255

that the Mohr-Coulomb failure criterion does not hold perfectly at the peak as friction depends on the intermediate principal stress component.

256

257

Generally, $(q/p)_{crit}$ and $\phi_{crit}$ decrease monotonically until a minimum value is reached at $f_c = 40\%$ (Figs. 5b&d). Further increase in $f_c$ resulted in a monotonic increase such that $(q/p)_{crit}$ and $\phi_{crit}$ values for assemblies having $f_c \ge 60\%$ are similar. The assemblies having $f_c \ge 60\%$ fall within the finer fraction-dominated regime according to the DEM study conducted by Adesina et al., [16] on gap-graded assemblies subjected to drained triaxial compression. The data presented here shows that an increase in $f_c$ within this overfilled, fine dominated regime does not have a significant influence on $(q/p)_{crit}$ and $\phi_{crit}$; all gravel-sand mixtures within the regime have a critical state strength similar to that of the clean sand. The marginal increase in $(q/p)_{crit}$ and $\phi_{crit}$ observed within $f_c = 10\%-20\%$, especially under triaxial compression agrees with earlier experimental studies on Ottawa sand and Carmague silty sand subjected to drained triaxial compression [43–45].

258

259

260

261

262

263

264

265

266

267

268

The influence of $f_c$ on the critical state shear strength appears counterintuitive when compared with $e_{crit}$ in that the initial general reduction in $(q/p)_{crit}$ and $\phi_{crit}$ as $f_c$ increases is accompanied by a reduction in $e_{crit}$; ditto for the subsequent increase in $(q/p)_{crit}$ and $\phi_{crit}$ (Fig. 5b). However, in the soil mechanics literature, a reduction in shear strength is generally associated with an increase in void ratio. It has been shown in the DEM study by Adesina et al., [16] that the influence of $f_c$ on $(q/p)_{crit}$ can be better understood using the stress-based skeleton (SBS) void ratio proposed in the study. Adesina et al., [16] established a linear correlation between the SBS void ratio at the critical state and $(q/p)_{crit}$, and showed that no correlation exist between $e_{crit}$ and $(q/p)_{crit}$.

269

270

271

272

273

274

275

276

277

For all $f_c$ values considered, $(q/p)_{crit}$ values under triaxial compression are higher than the values obtained under triaxial extension. However, no significant difference was observed in the $\phi_{crit}$ values obtained under both loading paths, especially for $f_c \in [25; 50\%]$ (i.e. within the stress-based transitional regime identified in [16]). These data show that the dependence of the friction angle on the loading path diminishes as the material approaches critical state, such that a unique friction angle independent of the loading path is obtained. While our DEM data suggest that the Mohr-Coulomb criterion is not valid at the peak, the validity of the criterion at the critical state is not violated.

278

279

280

281

282

283

284

285

The influence of the finer fraction on the friction angle of the gap-graded assemblies considered is quantified as $\phi/\phi^{f_c=0\%}$ and presented in Fig. 6. As a result of the finer fraction, the peak friction angle increases by about 60% for $f_c=30\%$, while the critical state friction angle

286

287
288

decreases by 60% for $f_c$=40%, under both triaxial compression and extension loading paths.
289
The similar $\phi/\phi^{f_c=0\%}$ values obtained for the triaxial compression and extension loading paths
290
indicate that the influence of the finer fraction on the friction angle of the assemblies is
291
independent of the loading path. This is true for the peak and the critical state alike. This
292
observation suggests an isotropic contribution of the finer fraction to the stress under both
293
loading paths.

295

296

297

298

<figure><img src="image_12.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d)</figcaption></figure>

294

**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}$ (c) friction angle at peak, $\phi_{peak}$ (d) critical state friction angle, $\phi_{crit}$.

295

296

297

298
<figure><img src="image_16.png" /><figcaption>Fig. 6 The influence of finer fraction on the friction angle at the peak and the critical state</figcaption></figure>

Dense granular materials dilate (i.e. increase in volume) when subjected to shearing. This dilatancy of dense granular materials is characterised by breaking of interlocks and loss of interparticle contacts, in addition to a simultaneous strain softening [46,47]. The higher maximum dilatancy rate, $(-d\epsilon_v/d\epsilon_a)_{max}$, obtained under compression in comparison to extension (Fig. 7b) indicates a tendency for higher strain softening under triaxial compression than in extension, where strain softening $=(q/p)_{peak} - (q/p)_{crit}$ (Fig. 7a). Furthermore, the higher dilatancy exhibited under compression also suggests that higher amount of energy is required to break particle interlocks within the assemblies subjected to compression in comparison to extension. This relatively more pronounced dilatancy observed under compression indicates more kinematic constraints in compression than extension. This tendency can be attributed to the increase in pre-peak mean pressure in compression, whereas the pre-peak mean pressure decreases in extension. Dilatancy increases with packing density. This is true even for assemblies of binary mixtures such that the maximum dilatancy angle, $\psi_{max}$, decreases linearly with an increase in the initial void ratio, where $\psi_{max} = \sin^{-1}[((-d\epsilon_v/d\epsilon_a)_{max})/(2 + (-d\epsilon_v/d\epsilon_a)_{max})]$. However, the rate of the decline in dilatancy with an increase in void ratio is faster under compression than in extension (Fig. 7c). In agreement with earlier experimental and numerical studies on clean sands and sand-silt mixtures [47–50], we confirm a linear relationship between the peak friction angle and the maximum dilatancy angle for the assemblies of binary mixtures studied (Fig. 8).
The comparison between the friction angle and the dilatancy angle is central in understanding the non-associated character of granular plasticity. The difference $\phi_{peak} - \psi_{max}$ can be used to quantify the extent to which the flow rule is non-associated. Based on the second-order work criterion for mechanical stability [51], the greater the difference between these angles, the more likely is a material to fail. Therefore, the slope of $\phi_{peak}-\psi_{max}$ with respect to the fine content can indicate whether the finer fraction have a stabilisation effect. Slope > 1 indicates an increasing deviation from the associated condition as the finer fraction increases and can be related to a decrease in material stability. We observed that under compression, the slope of $\phi_{peak}-\psi_{max}$ is 0.80 for the underfilled assemblies (i.e. $f_c<f_c^{th}$). This implies an increase in material stability as $f_c$ increases, in agreement with Wang et al. [51]. Conversely, increasing $f_c$ within the overfilled assemblies ($f_c>f_c^{th}$) may constitute a decrease in material stability as the slope = 1.36. When considering the triaxial extension case, the slope $\phi_{peak}-\psi_{max}$ is equal to 1.81 and 3.47 for the underfilled assemblies and overfilled assemblies, respectively. This suggests that an increase in $f_c$ may result in a less stable material when subjected to extension in comparison to compression.

<figure><img src="image_17.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(c)</figcaption></figure>
339

**Fig. 7** Effect of fines content on the: (a) strain softening, $\phi_{peak} - \phi_{crit}$ (b) maximum dilatancy rate, $(-d\epsilon_v/d\epsilon_a)_{max}$ for dense assemblies subjected to triaxial compression and extension. Relationship between (c) maximum dilatancy angle, $\psi_{max}$ and initial void ratio, $e_{min}$, for dense assemblies of binary mixture subjected to triaxial compression and extension

340

341

342

343

344

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

345

346

347

348

349

350

**Fig. 8** Relationship between peak friction angle, $\phi_{peak}$, and maximum dilatancy angle, $\psi_{max}$ for dense assemblies of binary mixture subjected to triaxial compression and extension. The data points using pentagram and hexagram as markers represents the data for the threshold finer fraction content, $f_c^{th}$

351

## 3.2 Micromechanical characteristics

352

353

354

355

356

357

According to Thornton [52] and Edwards [53], the limiting condition for interlocking equilibrium in granular assemblies of frictional spheres is set at a mechanical coordination number, $Z_m = 4$. $Z_m$ is the average number of contacts made by individual particles in a granular assembly when the rattlers, i.e. the particles with less than two contacts, which do not contribute to the contact force network are excluded. $Z_m < 4$ indicate a packing too weakly connected to be at an interlocking equilibrium. In this case, internal kinematic degrees of
freedom exists so that the particles with less than four contacts may move easily. Following
Thornton (2000), $Z_m$ is defined as<fnref n="5" />:

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

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. Fig. 9a shows the value of $Z_m$ for the whole assembly considering all particle size fractions at the initial, the peak, and the critical states. At all shearing stages and loading paths considered, $Z_m$ decreases with an increase in $f_c$ until reaching the $f_c^{th}$, indicating a decline in interlocking. Beyond $f_c^{th}$, $Z_m$ increases monotonically with an increase in $f_c$ such that a plateau is reached for $f_c > \approx 60\%$. A similar trend is shown in the DEM study by Liu et al., [54] for sand-silt mixtures having particle size ratio, $\lambda = 8.4$ and 18.1, and at the initial states. In general, $Z_m \ge 4$ for all assemblies at the initial states and the peak, indicating these assemblies were in interlocking equilibrium. However, at the critical state, $Z_m < 4$ for the assemblies having $f_c = 35\%$ and $f_c = 40\%$, indicating a transition from an interlocked to a less interlocked state during shearing. These assemblies exhibited relatively large fluctuation in their $q/p$ and $\epsilon_{vol}$ data at the critical state, as shown in Figs. 16-19 in Appendix A.

In Fig. 9b & c, we quantify the average mechanical contact<fnref n="6" /> gained or lost in each assembly during shearing by determining $\Delta Z_m^{peak} = Z_m^{peak} - Z_m^{ini}$, and $\Delta Z_m^{crit} = Z_m^{crit} - Z_m^{ini}$ where $Z_m^{ini}$, $Z_m^{peak}$ and $Z_m^{crit}$ are the mechanical coordination numbers at the initial state, the peak and the critical state, respectively. $\Delta Z_m > 0$ means contact gained and $\Delta Z_m < 0$ means contact lost during shearing. We observed that at the peak (Fig. 9b), the dense assemblies under triaxial extension had lost more mechanical contacts than under compression. However, for all $f_c$ considered, $\Delta Z_m$ values under both compression and extension loading paths were similar at the critical state (Fig. 9c). The initial delay in contact loss (with respect to the axial strain) during compression can be attributed to the initial contraction exhibited by the dense assemblies during triaxial compression which was not observed during extension (see insets of Figs. 16d-f & Figs. 17d-f). This indicates that an initially more connected interparticle network with a stronger particle interlock is created under compression than in extension. A stronger particle interlock requires higher amount of energy to break during dilation, hence, the higher dilatancy and lower amount of contact loss observed during compression than in extension (Fig. 7b &

<footnote n="5">⁵ 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="6">⁶ Mechanical contacts are the contacts made by particles with more than two contacts</footnote>
9b). Also, if less amount of energy is required to break particle interlock as indicated by the lower dilatancy exhibited under extension, it is reasonable to expect a more rapid contact loss under extension than in compression, at the peak.

<figure><img src="image_21.png" /><figcaption>**Fig. 9** Effect of fines content on the (a) mechanical coordination number, $Z_m$, for all particles. Change in mechanical coordination number during shearing, $\Delta Z_m$, (b) dense assemblies at the peak (c) dense and loose assemblies at the critical state</figcaption></figure>

**4. Contribution of particle fractions and contact types to stress transmission under triaxial compression and extension loading paths**

The proportion of the stress transmitted by the coarser and finer fractions have been used in the literature on gap-graded soils to understand the role played by each particle size fraction on the mechanical behaviour of the soils [26,27,55]. In this study, we aim to determine the influence of the loading path (i.e. triaxial compression or extension here) on the proportion of the total mean stress transmitted by the coarser and the finer fractions. Following the Love-Weber particle stress tensor as described in Nicot et al., (2013), the stress on an individual particle, $\sigma_{ij}^p$, is calculated as:
407

$$ \sigma_{ij}^{p} = \frac{1}{V^p} R_p \sum_{c=1}^{N_c^p} f_j^c \mathbf{n}_i \qquad (2) $$

408

409

410

411

412

where the summation is performed over all contacts c made by the particle p. $V^p$ is the volume of the particle, $R_p$ is the particle radius, $N_c^p$ is the number of contacts made by the particle, $f_j^c$ is the contact force vector and $\mathbf{n}_i$ is the unit branch vector between two particles in contact. The stress tensor at the assembly scale, $\sigma_{ij}$, is given as a weighted average of the stresses sustained by the particles:

413

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

414

415

where $V$ is the volume of the entire assembly and $N_p$ is the number of particles in the assembly. We can decompose $\sigma_{ij}$ to the contribution of the finer and the coarser fraction as:

416

$$ \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 \qquad (4) $$

417

418

where $N_p^c$ and $N_p^f$ are the number of coarser and finer particles in the granular assembly, respectively<fnref n="7" />.

419

420

Supposing the stress of interest is the mean effective stress, the total mean effective, $p$, transmitted by the entire assembly of particles is determined as:

421

$$ p = \frac{\sigma_{xx} + \sigma_{yy} + \sigma_{zz}}{3} \qquad (5) $$

422

423

424

425

Therefore, the proportion of the total mean stress transmitted by the finer fraction, $\alpha_f$, is determined as $p^f/p$. Similarly, the proportion of the total mean stress transmitted by the coarser fraction, $\alpha_c = 1 - \alpha_f$, is determined as $p^c/p$, where $p^f$ and $p^c$ are the mean stress transmitted by the finer and coarser fractions, respectively.

426

427

428

429

In Fig. 10, we show the contributions to the total mean stress, of the finer fraction, $\alpha_f$, and the coarser fraction, $\alpha_c$, under triaxial compression and extension loading paths, at the initial state, the peak, and the critical state. As detailed in Adesina et al., [16], the contributions of the particle size fractions to the total mean stress can be classified into the coarse-dominated

<footnote n="7">⁷ 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>
regime where $\alpha_f \approx 0$, the transitional regime where $0 < \alpha_f < \alpha_c$ and the fines-dominated regime where $\alpha_f > \alpha_c$. For all $f_c$ value, $\alpha_f$ and $\alpha_c$ are similar under both triaxial compression and extension (with $\alpha_f$ being only slightly lower in extension than in compression within $f_c^{th} < f_c < 80\%$) (Fig. 10c and d). This indicate that there is no significant influence of the loading paths on the contribution of the particle fractions to the total mean stress. This observation is supported by the data in Fig. 11 which shows no significant influence of the loading paths on the total number of contacts made by the coarser and finer fraction. The lower total number of contacts made by the finer fraction within $f_c^{th} < f_c < 80\%$, under triaxial extension in comparison to the compression loading path (Fig. 11c) is consistent with the data in Fig. 10c where the proportion of the mean stress transmitted under triaxial extension is lower than in compression for similar range of $f_c$.

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

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

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

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

- coarse; compression

- coarse; extension

- fine; compression

- fine; extension
443
444
445
446

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

447
448
449
450
451
452
453
454
455
456

Similar to the analysis on the contribution of the particle size fractions to the total mean stress, we estimated the contribution of the contact types to the total mean stress to understand the influence of loading paths on contact interactions within the assemblies. The Love-Weber stress tensor including Bagi boundary term [57], give the contributions of all contact types to the total mean stress transmitted by an entire granular assembly as:

457

$$ \sigma_{ij}^{c} = \frac{1}{V} \sum_{c=1}^{N_{c}} f_{j}^{c} l_{i} \qquad (6) $$
where $N_c$ is the total number of contacts in the granular assembly, $V$ is the assembly volume, $l$ is the branch vector between elements in contact ($\|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. A decomposition of this contact-based stress tensor to the contributions of each contact type can be performed using:

$$ \sigma_{ij}^k = \frac{1}{V} \sum_{c=1}^{N_c^k} f_j^c l_i \quad (7) $$

where $k \in \{\text{c-c}, \text{c-f}, \text{f-f}, \text{s-wl}\}$ with c-c, c-f, f-f, and s-wl denoting coarse to coarse, coarse to fine, fine to fine and sphere to wall contacts, respectively. The contribution made by each contact type to the total mean stress is therefore given as $\alpha_k = p^k/p$. Note that the boundary term is expected to be negligible when the REV condition is met.

Fig. 12 shows the proportion of the total mean stress transmitted by each contact types, at the initial states, the peak and the critical state for both triaxial compression and extension loading paths. In agreement with earlier experimental and numerical studies on gap-graded assemblies [27–29], the contributions made by each contact type to the total mean stress can be classified into four categories detailed in Adesina et al., [16] and summarised here as:

- The coarse-dominated zone where $\alpha_{c-c} > \alpha_{c-f} > \alpha_{f-f}$

- The transitional coarse-dominated zone where $\alpha_{c-f} > \alpha_{c-c} > \alpha_{f-f}$

- The transitional fine-dominated zone where $\alpha_{c-f} > \alpha_{f-f} > \alpha_{c-c}$

- The fines-dominated zone where $\alpha_{f-f} > \alpha_{c-f} > \alpha_{c-c}$

Similar to the observation made in Fig. 10, the compression and extension loading paths did not significantly influence the contributions of the contact types to the total mean stress. Within $f_c^{th} < f_c < 80\%$, the contribution of the fine-fine contacts to the total mean stress was only slightly lower under triaxial extension in comparison to the compression loading path (Fig. 12c & d). The contribution of the coarse-fine contacts to the total mean stress was also slightly lower in extension than in compression within $f_c^{th} < f_c < 50\%$, especially at the peak.

As noted earlier, provided that the REV condition is met, the contribution to the total mean stress from the wall-sphere interactions, $\alpha_{s-wl}$, should be negligible. At the critical state, $\alpha_{s-wl} < 0.03$ for the clean gravel ($f_c = 0\%$) and the clean sand ($f_c = 100\%$) assemblies indicating that these assemblies are representative. $\alpha_{s-wl}$ is relatively higher within $10\% < f_c \le 40\%$ where particle to wall interaction is observed to increase to a peak before declining with further
488

increase in $f_c$. This observation agrees with the data in the DEM study by Shi et. al [17] where it has been shown that the RVE sizes for gap-graded granular assemblies increases as the coarser fraction content increases, i.e. as $f_c$ decreases. The very similar $\alpha_{s-wl}$ values obtained under the compression and extension loading paths imply that the REV condition is not significantly influenced by the loading path.

489

490

491

492

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

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

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

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

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

493

494

495

496

497

498

499

**Fig. 12** 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

Fig. 13 shows the proportion of the total contacts belonging to each contact type at the initial states, the peak and the critical state and for both triaxial compression and extension loading paths. The trends observed here can be grouped into three categories:

502

- The coarse-coarse contact dominated zone where $C_{c-c}^t > C_{c-f}^t > C_{f-f}^t$

503

- The transitional coarse-fine contact dominated zone where $C_{c-f}^t > C_{c-c}^t > C_{f-f}^t$
- The fine-fine contact dominated zone where $C_{f-f}^t > C_{c-f}^t > C_{c-c}^t$

The REV analysis presented in Adesina et al., [16] confirms that the 100,000 particles used satisfy the REV condition for all $f_c$ considered in this study. The higher particle to wall interaction obtained for the assemblies having $10\%<f_c\ll 40\%$ (Fig. 13c & d) in comparison to other assemblies confirms that a larger sample size may be required to satisfy the REV for lower finer fraction content.

<figure><img src="image_32.png" /><figcaption>**Fig. 13** 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>
## 5. Influence of triaxial compression and extension loading paths on the fabric of gap-graded assemblies

A number of previous studies have emphasised the influence of inherent fabric on the critical state characteristics of sands when subjected to different loading paths [9,58–60]. Meanwhile, the critical state theory first postulated by Casagrande [61] was based on the assumption that the large deformations required to reach the critical state are such that any influence of initial fabric and loading path would have vanished. This means that a unique behaviour is expected at the critical state regardless of loading path or initial inherent anisotropy. In order to investigate the influence of loading path on fabric in this study, we consider the evolution of the contact normal anisotropy of gap-graded assemblies prepared to an isotropic stress state and subjected to drained triaxial compression and extension shearing. An initial isotropic fabric exempts the influence of initial inherent anisotropy and ensures that the results obtained emanate solely from the influence of the loading paths. A contact normal fabric tensor ($\Phi_{cij}$) defined by Satake [62] is calculated as:

$$ \Phi_{cij} = \frac{1}{N_c} \sum_{c=1}^{N_c} n_{ci}^c n_{cj}^c \quad (8) $$

where $N_c$ is the total number of contacts in the granular assembly and the vector $n_c$ gives the contact normal orientation.

Following Barreto et al. [63], the contact normal fabric anisotropy, $a_{cn}$, for 3-D assemblies is given as:

$$ a_{cn} = \frac{1}{\sqrt{2}} \sqrt{(\Phi_{c1} - \Phi_{c3})^2 + (\Phi_{c3} - \Phi_{c2})^2 + (\Phi_{c1} - \Phi_{c2})^2} \quad (9) $$

where $\Phi_{c1}$, $\Phi_{c2}$ and $\Phi_{c3}$ are the major, intermediate and the minor eigenvalues of the tensor. The fabric anisotropy varies from zero when the assemblies are perfectly isotropic to unity when the major axes of all contacts align along the major principal stress.

Fig. 14 shows the evolutions of $a_{cn}$ during shearing for the dense gap-graded assemblies subjected to triaxial compression (Fig. 14a-c) and triaxial extension (Fig. 14d-f) loading paths. All interparticle contacts are considered and grouped into the coarse-dominated, the transitional and the fines-dominated regimes following Adesina et al., [16]. In general, we observed that under both triaxial compression and extension loading paths, $a_{cn}$ evolves in a similar manner as the stress ratio, $q/p$ (see Figs. 16-19 in Appendix A). This is in agreement with the data presented in previous DEM studies where narrowly-graded granular assemblies were
547

considered [63–66]. The evolution of $a_{cn}$ can be divided into an initial linear phase (indicated in the insets of Fig. 14 which show $a_{cn}$ at small strain, i.e. $\epsilon_a \le 0.002$), a non-linear phase during which $a_{cn}$ for dense assemblies increases to a peak and thereafter decreases, and a critical state where $a_{cn}$ reaches a plateau. We observed that the fluctuation in the $q/p$ data for the underfilled and transitional assemblies (see Figs. 16 - 19 in Appendix A) is reflected in the $a_{cn}$. This observation signifies a rapidly evolving and less stable fabric for the underfilled and transitional assemblies. This is also reflected in their coordination number at the critical state (see Fig. 9a and discussion in Section 3.2).

548

549

550

551

552

553

554

555

Fig. 15 shows the maximum $a_{cn}$, $a_{cn}^{max}$, and the critical state $a_{cn}$, $a_{cn}^{crit}$, obtained during shearing for all assemblies. $a_{cn}^{max}$ initially increases with an increase in $f_c$ within the underfilled assemblies (i.e. $f_c<f_c^{th}$). Further increase in $f_c$ resulted in a drastic decrease in $a_{cn}^{max}$ followed by a monotonic increase (Fig. 15a). Both $a_{cn}^{max}$ and $a_{cn}^{crit}$ reflect the void ratio of the assemblies. In general, the lower the void ratio, the lower the anisotropy exhibited by the assemblies especially at the critical state (Fig. 15b). This is in agreement with the finding in the DEM study by Adesina et al [67] where assemblies of varying $f_c$ and $\lambda$ values were subjected to drained triaxial compression within periodic boundaries. The anisotropic tendency of a granular packing subjected to shearing is influenced by the available space for particle rearrangement within the packing. A denser packing where particles are closely packed would have less room for particle rearrangement than a looser packing with more space which allows for greater particles rearrangement. The general reduction in $a_{cn}$ as $f_c$ initially increases can be attributed to a reduction in the room for particle arrangement consequent upon the reduction in void ratio. The subsequent increase in $a_{cn}$ follows from the increase in the room for particle rearrangement as the void ratio increases. The evolution of anisotropy as $f_c$ increases is influenced by the stress contribution of the most abundant contact type within the granular assemblies [67]. At the critical state, the contribution to the mean stress of the coarse-coarse contacts, i.e. the most abundant within the underfilled assemblies (Fig. 13d), decreases with an increase in $f_c$ (Fig. 12d), thereby resulting in a more isotropic contact network and a decrease in the overall critical state $a_{cn}$. However, the contribution of the fine-fine contacts, the most abundant contact type within the overfilled assemblies, increases monotonically with an increase in $f_c$, thereby yielding a more anisotropic contact network and an increase in the critical state $a_{cn}$.

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

We observed that $a_{cn}$ is, in general, marginally higher under triaxial extension than under compression. This observation may indicate that during triaxial extension, less contacts are

578
opened perpendicular to the principal stress directions in comparison to compression where more contacts are opened perpendicular to the major principal stress axis. Our data show that critical state fabric is influenced by the loading path. Despite having the same initial fabric, the contact normal anisotropy values obtained during triaxial compression differs from the values obtained during extension. Therefore, the notion of a unique critical state fabric that is independent of loading path is not supported by the DEM data presented here. While inherent fabric can influence the critical state characteristics of granular materials as shown in the constitutive analysis by Li et al., [9], the sensitivity of the critical state fabric to the loading path obtained in this study does not stem from an initial inherent fabric. Rather, it is the variation in the loading path that induces a change in the fabric and the critical state characteristics of the assemblies.

<figure><img src="image_33.png" /><figcaption>**Fig. 14** Effect of fines content, $f_c$, on the evolution of contact normal fabric anisotropy for dense assemblies: (a) $f_c$ = 0%-25%; triaxial compression (b) $f_c$ = 30%-40%; triaxial compression (c) $f_c$ = 45%-100%; triaxial compression (d) $f_c$ = 0%-25%; triaxial extension (b) $f_c$ = 30%-40%; triaxial extension (c) $f_c$ = 45%-100%; triaxial extension</figcaption></figure>

590

591

592

593

594

595

596

597

598
<figure><img src="image_34.png" /><figcaption>**Fig. 15** Effect of fines content, $f_c$, on (a) maximum contact normal fabric anisotropy during shearing, $a_{cn}^{max}$, (b) critical state contact normal anisotropy, $a_{cn}^{crit}$</figcaption></figure>

**6. Conclusions**

This study is aimed at identifying the key characteristics that differentiate gap-graded granular assemblies when subjected to triaxial compression and extension loading paths. Dense and loose gravel-sand mixtures of fines contents, $f_c \in [10\%; 90\%]$, were prepared along with clean gravel ($f_c = 0\%$) and clean sand ($f_c = 100\%$) assemblies. Each assembly was then subjected to drained triaxial compression and triaxial extension test series. The key characteristics that differentiate the behaviour of the assemblies under triaxial compression and extension are itemised below:

1. Our DEM simulations show that the evolutions of the shear stress under triaxial compression and extension loading paths are similar for the assemblies of binary mixtures studied. However, higher stress ratio and dilatancy were obtained under triaxial compression than in extension, in agreement with earlier experimental studies on narrowly graded and monodisperse sands and glass beads [1,41]. In contrast, while the mobilised friction at the peak is higher under extension than in compression in agreement with an earlier experimental studies [5,37,39–42], the mobilised friction values for both loading paths were similar at the critical state. Our DEM data confirms the validity of the Mohr-Coulomb criterion at the critical state.

2. While dense assemblies under both triaxial compression and extension loading paths dilated under shearing as expected, the initial contraction that precedes the dilation during triaxial compression was not observed under extension. A denser packing is therefore obtained before dilatancy under compression such that dilatancy is higher for the triaxial compression loading path in comparison to extension. This observed
624

625

626

627

628

629

3. Despite the observed differences in the overall stress transmitted by the gap-graded assemblies under triaxial compression and extension, the contributions of the finer and coarser fractions to the total mean stress remain very similar under both loading paths. Only a negligible reduction in the contribution of the finer fraction to the total mean stress was obtained under triaxial extension in comparison to triaxial compression, especially for the assemblies within $f_c^{th} < f_c < 80\%$ and at the peak state.

630

631

632

633

634

635

4. The identified differences in the characteristics of gap-graded assemblies subjected to drained triaxial compression and extension loading paths have implications for the development of generalised framework for predicting the behaviour of granular mixtures. Our DEM data reinforce the need to adequately capture the influence of loading paths (including varying intermediate principal stresses, q constant, simple shear, cyclic loading, etc) to develop more robust predictive models in the future.

636

637

638

639

640

641

642

Recommendations for future studies

643

The influence of the grading characteristics (e.g. $f_c$ and $\lambda$) on the critical state line (CSL) may be less apparent in the q-p' plane (see Fig. 3), in comparison with the e-p' plane. While the CSL for gap-graded assemblies in e-p' plane is beyond the scope of this study, it is important to consider this in future studies, particularly the combined influence of the finer grain content, the particle size ratio and the loading path on the critical state line in the e-p' plane. This will shed light on whether the loading path modifies the critical state characteristics of soil mixtures with different gradings. It will also enable the assessment of the appropriate void ratio index (skeleton, equivalent or mechanical ratios proposed in [14-16,32,36,44,69,70]) for constructing the CSL, and to understand the influence of loading path on the CSL in e-p' plane for different void ratio indexes.

644

645

646

647

648

649

650

651

652

653

While the influence of intermediate stress ratio, $b$, has been considered in the literature for uniformly graded assemblies [3,37-42], studies on the interplay between grading and intermediate stress are lacking for gap-graded assemblies. Focusing solely on triaxial

654

655

656
657 compression and extension loading paths in the present study enabled the simulation of a
658 wide range of finer fraction content. The influence of intermediate stress ratio and other
659 stress paths including plane strain condition which are beyond the scope of the present
660 study should be considered in future studies.

661

662 Notations

<table><tr><td>663</td><td>C</td><td>Total number of contacts in a granular system</td></tr><tr><td>664</td><td>C<sub>p</sub><sup>c</sup></td><td>Total number of contacts made by the coarser particles</td></tr><tr><td>665</td><td>C<sub>p</sub><sup>f</sup></td><td>Total number of contacts made by the finer particles</td></tr><tr><td>666</td><td>C<sub>c-c</sub><sup>t</sup></td><td>The proportion of all contact types belonging to the c-c contacts</td></tr><tr><td>667</td><td>C<sub>c-f</sub><sup>t</sup></td><td>The proportion of all contact types belonging to the c-f contacts</td></tr><tr><td>668</td><td>C<sub>f-f</sub><sup>t</sup></td><td>The proportion of all contact types belonging to the f-f contacts</td></tr><tr><td>669</td><td>C<sub>wl</sub><sup>t</sup></td><td>The proportion of all contact types belonging to the wall</td></tr><tr><td>670</td><td>D<sub>c</sub></td><td>Diameter of the coarser grains</td></tr><tr><td>671</td><td>D<sub>f</sub></td><td>Diameter of the finer grains</td></tr><tr><td>672</td><td>e<sub>max</sub></td><td>Maximum void ratio</td></tr><tr><td>673</td><td>e<sub>min</sub></td><td>Minimum void ratio</td></tr><tr><td>674</td><td>e<sub>peak</sub></td><td>Void ratio at the peak state</td></tr><tr><td>675</td><td>e<sub>crit</sub></td><td>Void ratio at the critical state</td></tr><tr><td>676</td><td>ε<sub>a</sub></td><td>Axial strain (in z-direction)</td></tr><tr><td>677</td><td>ε<sub>vol</sub></td><td>Volumetric strain</td></tr><tr><td>678</td><td>f<sub>c</sub></td><td>Fines content</td></tr><tr><td>679</td><td>f<sub>c</sub><sup>th</sup></td><td>Threshold fines content</td></tr><tr><td>680</td><td>f<sub>j</sub><sup>c</sup></td><td>Contact force vector</td></tr><tr><td>681</td><td>l</td><td>Unit branch vector</td></tr><tr><td>682</td><td>N<sub>p</sub></td><td>Number of particles in the granular system</td></tr><tr><td>683</td><td>N<sub>c</sub></td><td>Number of contacts in a granular system</td></tr><tr><td>684</td><td>N<sub>c</sub><sup>k</sup></td><td>Number of kth contact type in a granular system</td></tr><tr><td>685</td><td>N<sub>p</sub><sup>c</sup></td><td>Number of coarser particles in the system</td></tr><tr><td>686</td><td>N<sub>p</sub><sup>f</sup></td><td>Number of finer particles in the system</td></tr><tr><td>687</td><td>N<sub>c</sub><sup>p</sup></td><td>Number of contacts involving a particle</td></tr><tr><td>688</td><td>N<sub>0</sub></td><td>Number of particles with no contact</td></tr><tr><td>689</td><td>N<sub>1</sub></td><td>Number of particles with one contact</td></tr><tr><td>690</td><td>n<sub>i</sub></td><td>Unit branch vector</td></tr><tr><td>691</td><td>p</td><td>Mean effective stress</td></tr><tr><td>692</td><td>p<sup>c</sup></td><td>Mean stress transmitted by the coarser fraction</td></tr><tr><td>693</td><td>p<sup>f</sup></td><td>Mean stress transmitted by the finer fraction</td></tr><tr><td>694</td><td>p<sup>k</sup></td><td>Mean stress transmitted by a contact type</td></tr><tr><td>695</td><td>α<sub>c</sub></td><td>Proportion of the total mean stress transmitted by the coarser fraction</td></tr><tr><td>696</td><td>α<sub>f</sub></td><td>Proportion of the total mean stress transmitted by the finer fraction</td></tr><tr><td>697</td><td>α<sub>k</sub></td><td>Proportion of the total mean stress transmitted a contact type</td></tr><tr><td>698</td><td>(q/p)<sub>peak</sub></td><td>Stress ratio at the peak</td></tr><tr><td>699</td><td>(q/p)<sub>crit</sub></td><td>Stress ratio at the critical state</td></tr><tr><td>700</td><td>q</td><td>Deviatoric stress</td></tr></table>
<table>
  <tr>
    <td>701</td>
    <td>R<sub>p</sub></td>
    <td>Radius of particle</td>
  </tr>
  <tr>
    <td>702</td>
    <td>R<sub>1</sub></td>
    <td>Radius of first particle in contact</td>
  </tr>
  <tr>
    <td>703</td>
    <td>R<sub>2</sub></td>
    <td>Radius of second particle in contact</td>
  </tr>
  <tr>
    <td>704</td>
    <td>μ</td>
    <td>Inter-particle friction coefficient</td>
  </tr>
  <tr>
    <td>705</td>
    <td>λ</td>
    <td>Size fraction</td>
  </tr>
  <tr>
    <td>706</td>
    <td>V</td>
    <td>Volume of system</td>
  </tr>
  <tr>
    <td>707</td>
    <td>V<sup>p</sup></td>
    <td>Volume of particle</td>
  </tr>
  <tr>
    <td>708</td>
    <td>σ<sub>zz</sub></td>
    <td>Normal stress in z-direction</td>
  </tr>
  <tr>
    <td>709</td>
    <td>σ<sub>yy</sub></td>
    <td>Normal stress in y-direction</td>
  </tr>
  <tr>
    <td>710</td>
    <td>σ<sub>xx</sub></td>
    <td>Normal stress in x-direction</td>
  </tr>
  <tr>
    <td>711</td>
    <td>σ<sub>ij</sub></td>
    <td>Second order stress tensor</td>
  </tr>
  <tr>
    <td>712</td>
    <td>σ<sub>ij</sub><sup>c</sup></td>
    <td>Second order stress tensor for the coarser grains</td>
  </tr>
  <tr>
    <td>713</td>
    <td>σ<sub>ij</sub><sup>f</sup></td>
    <td>Second order stress tensor for the finer grains</td>
  </tr>
  <tr>
    <td>714</td>
    <td>σ<sub>ij</sub><sup>p</sup></td>
    <td>Second order stress tensor for a particle</td>
  </tr>
  <tr>
    <td>715</td>
    <td>σ<sub>ij</sub><sup>c</sup></td>
    <td>Second order stress tensor for all contact types</td>
  </tr>
  <tr>
    <td>716</td>
    <td>φ</td>
    <td>Friction angle</td>
  </tr>
  <tr>
    <td>717</td>
    <td>φ<sub>peak</sub></td>
    <td>Peak friction angle</td>
  </tr>
  <tr>
    <td>718</td>
    <td>φ<sub>crit</sub></td>
    <td>Critical state friction angle</td>
  </tr>
  <tr>
    <td>719</td>
    <td>Φ<sub>cij</sub></td>
    <td>Second order contact normal fabric tensor</td>
  </tr>
  <tr>
    <td>720</td>
    <td>Φ<sub>c1</sub></td>
    <td>Major eigenvalue of the contact normal fabric tensor</td>
  </tr>
  <tr>
    <td>721</td>
    <td>Φ<sub>c2</sub></td>
    <td>Intermediate eigenvalue of the contact normal fabric tensor</td>
  </tr>
  <tr>
    <td>722</td>
    <td>Φ<sub>c3</sub></td>
    <td>Minor eigenvalue of the contact normal fabric tensor</td>
  </tr>
  <tr>
    <td>723</td>
    <td>n<sub>c</sub></td>
    <td>Unit vector representing contact normal orientation</td>
  </tr>
  <tr>
    <td>724</td>
    <td>a<sub>cn</sub></td>
    <td>Contact normal anisotropy</td>
  </tr>
  <tr>
    <td>725</td>
    <td>a<sub>cn</sub><sup>max</sup></td>
    <td>Maximum contact normal anisotropy</td>
  </tr>
  <tr>
    <td>726</td>
    <td>a<sub>cn</sub><sup>crit</sup></td>
    <td>Contact normal anisotropy at the critical state</td>
  </tr>
  <tr>
    <td>727</td>
    <td>ψ<sub>max</sub></td>
    <td>Maximum dilatancy angle</td>
  </tr>
  <tr>
    <td>728</td>
    <td>Z<sub>m</sub></td>
    <td>Mechanical coordination number</td>
  </tr>
  <tr>
    <td>729</td>
    <td>Z<sub>m</sub><sup>ini</sup></td>
    <td>Mechanical coordination number at the initial state</td>
  </tr>
  <tr>
    <td>730</td>
    <td>Z<sub>m</sub><sup>peak</sup></td>
    <td>Mechanical coordination number at the peak</td>
  </tr>
  <tr>
    <td>731</td>
    <td>Z<sub>m</sub><sup>crit</sup></td>
    <td>Mechanical coordination number at the critical state</td>
  </tr>
</table>

733

734

735

736

737

738

739

740

741

742
743

744

745

746

747

## Appendix A. Shearing responses under triaxial compression and extension loading paths

748

749

750

751

752

753

754

755

756

757

758

759

760

761

762

763

764

765

Figs. 16 and 17 show the evolution of the stress ratio $q/p$, and the volumetric strain ($\epsilon_{vol}$), against the major principal axial strain<fnref n="8" /> ($\epsilon_a$), for dense assemblies of varying fines content subjected to drained triaxial compression and extension, respectively. Figs. 18 and 19 show the corresponding plots for the loose assemblies. The assemblies were sheared until $\epsilon_a \approx 0.5$ to analyse their characteristics at the critical state. The shearing responses in Figs. 16-19 were divided into three classes and subplots based on the findings in Adesina et al., [16] where the assemblies with $f_c = 0\%-25\%$, $f_c = 30\%-40\%$ and $f_c = 45\%-100\%$ fall within the coarse-dominated, the transitional and the fines-dominated regimes, respectively. The inset in each subplot shows the influence of $f_c$ on the shearing responses at small strain (i.e. $\epsilon_a \le 0.02$).

Excluding the clean gravel assembly, ($f_c = 0$), higher fluctuations were observed in the large strain shearing responses of the coarse-dominated and the transitional assemblies which contained lower number of contacts in comparison to the fines-dominated assemblies [16,71].

In general, the evolutions of the shearing responses under both triaxial compression and extension loading paths are generally similar. The main difference between the loading paths was identified in the volumetric response in that the initial contraction exhibited by the dense assemblies under compression was not found under extension. The implications of this behaviour on the dilatancy of the assemblies under both loading paths is discussed in Section 3.2.

<footnote n="8">⁸ The axial strain in the loading direction, $\epsilon_a$ (i.e. $\epsilon_{zz}$) is negative under compression and positive under extension. In Figs. 16-19, $\epsilon_a$ is made positive for both compression and extension to facilitate the comparison of the stress-strain and volumetric strain evolution under both loading paths.</footnote>
<figure><img src="image_35.png" /><figcaption>766 Fig. 16 Effect of fines content, $f_c$, on the stress-strain and the volumetric strain responses for
768 dense assemblies under drained triaxial compression (a,d) $f_c$ = 0%-25% (b,e) $f_c$ = 30%-40%
769 (c,f) $f_c$ = 45%-100%.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>771 Fig. 17 Effect of fines content, $f_c$, on the stress-strain and the volumetric strain responses for
772 dense assemblies under drained triaxial extension (a,d) $f_c$ = 0%-25% (b,e) $f_c$ = 30%-40% (c,f)
773 $f_c$ = 45%-100%.</figcaption></figure>
<figure><img src="image_37.png" /><figcaption>776 Fig. 18 Effect of fines content, f_c, on the stress-strain and the volumetric strain responses for
777 loose assemblies under drained triaxial compression (a,d) f_c = 0%-25% (b,e) f_c = 30%-40%
778 (c,f) f_c = 45%-100%.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>779 Fig. 19 Effect of fines content, f_c, on the stress-strain and the volumetric strain responses for
780 loose assemblies under drained triaxial extension (a,d) f_c = 0%-25% (b,e) f_c = 30%-40% (c,f)
781 f_c = 45%-100%.</figcaption></figure>

783

784

785
## Appendix B. Particle-scale REV analysis

In the earlier DEM study by Adesina et al [16] where only triaxial compression was considered, it was shown that the number of particles required to satisfy REV follows a bell-shaped curve, with the maximum number of particles required around the threshold finer fraction content, $f_c^{th}$. Given that $f_c^{th} = 30\%$ for the dense assemblies in this study, Fig. 20 shows the cumulative distributions of the individual particle stresses normalised by the average particle stress in an assembly ($P^p/\bar{P}^p$), for dense assemblies having $f_c = 30\%$, and subjected to triaxial compression. As shown in Fig. 20, increasing the number of particles from 100,000 to 200,000 particles did not significantly influence the distribution of stresses on both the coarser (Fig. 20a) and finer particles (Fig. 20b) in the assemblies. This further confirms that the number of coarse grains in the assemblies is sufficient to capture the statistical distribution of particle stresses. Assembly-scale REV analysis showing the stress-strain responses for dense and loose assemblies with $f_c = 30\%$ and having 100,000 and 200,000 particles are presented in Adesina et al [16]. While the analysis in Fig. 20 is for the assemblies under compression, we expect a similar result under extension given that loading path did not significantly influence the proportion of mean stress transmitted by the coarser and finer particle fractions as shown in Fig. 12.

<figure><img src="image_39.png" /><figcaption>**Fig. 20** Effect of sample size on the cumulative distribution of normalised particle mean stress for: (a) coarser particles (b) finer particles within dense assemblies with $f_c = 30\%$ under triaxial compression.</figcaption></figure>
811 Acknowledgement

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

818

819 Credit Authorship

820 Peter Adesina: Conceptualization, Methodology, Investigation, Formal analysis,
821 Visualization, Writing – original draft.

822 Antoine Wautier: Funding acquisition, Conceptualization, Methodology, Supervision,
823 Writing – review & editing.

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

826

827 Compliance with ethical standards

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

831

832 References

833

834 1. Salvatore E, Modoni G, Andò E, Albano M, Viggiani G. Determination of the critical
835 state of granular materials with triaxial tests. Soils Found. 2017;57(5):733-44.

836 2. Yamamuro JA, Lade P V. Strain Localization in Extension Tests on Granular Materials. J Eng Mech. 1995;121(7):828-36.

838 3. Yoshimine M, Ishihara K, Vargas W. Effects of principal stress direction and
839 intermediate principal stress on undrained shear behavior of sand. Soils Found.
840 1998;38(3):179-88.

841 4. Yamamuro JA, Lade P V. Drained Sand Behavior in Axisymmetric Tests at High Pressures. J Geotech Eng. 1996;122(2):109-19.

843 5. Corfdir A, Sulem J. Comparison of extension and compression triaxial tests for dense sand and sandstone. Acta Geotech. 2008;3(3):241-6.

845 6. Strahler A, Stuedlein AW, Arduino PW. Stress-Strain Response and Dilatancy of Sandy Gravel in Triaxial Compression and Plane Strain. J Geotech Geoenvironmental Eng. 2016;142(4):1-11.

848 7. Kuwano R, Santa Spitia LF, Bedja M, Otsubo M. Change in mechanical behaviour of gap-graded soil subjected to internal erosion observed in triaxial compression and torsional shear. Geomech Energy Environ [Internet]. 2021;27:100197. Available from:
849
850
851 https://doi.org/10.1016/j.gete.2020.100197

852 8. Chitravel S, Otsubo M, Kuno M, Kuwano R. Post-erosion mechanical responses of internally unstable gap-graded soil under drained torsional simple shear and triaxial compression. Soils Found [Internet]. 2022;62(6):101224. Available from: https://doi.org/10.1016/j.sandf.2022.101224

856 9. Li XS, Dafalias YF. Constitutive Modeling of Inherently Anisotropic Sand Behavior. J Geotech Geoenvironmental Eng. 2002;128(10):868-80.

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

860 11. Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Geotechnique. 1979;30(3):331-6.

862 12. Craig RF. Craig's Soil Mechanics. Spon Press, Taylor & Francis Group, LONDON AND NEW YORK. 2004. 102-105 p.

864 13. Lade P V., Liggio CD, Yamamuro JA. Effects of Non-Plastic Fines on Minimum and Maximum Void Ratios of Sand. Geotech Test J. 1998;21(4):336-47.

866 14. Thevanayagam S, Shenthan T, Mohan S, Liang J. Undrained Fragility of Clean Sands, Silty Sands, and Sandy Silts. J Geotech Geoenvironmental Eng. 2002;128(10):849-59.

868 15. Rahman MM, Lo SR, Gnanendran CT. On equivalent granular void ratio and steady state behaviour of loose sand with fines. Can Geotech J. 2008;45(10):1439-56.

870 16. Adesina P, Wautier A, Benahmed N. A micromechanics-based classification of the regimes delineating the behaviour of gap-graded soils. Comput Geotech [Internet]. 2024;168. Available from: https://doi.org/10.1016/j.compgeo.2024.106165

873 17. Shi X, He Z, Zhao J, Liu J. Determination of the size of representative volume element for gap-graded granular materials. Powder Technol [Internet]. 2024; Available from: https://doi.org/10.1016/j.eneco.2024.107302

876 18. Jiang MD, Yang ZX, Barreto D, Xie YH. The influence of particle-size distribution on critical state behavior of spherical and non-spherical particle assemblies. Granul Matter [Internet]. 2018;20(4):1-15. Available from: https://doi.org/10.1007/s10035-018-0850-x

880 19. Thornton C. Numerical simulations of deviatoric shear deformation of granular media. Geotechnique. 2000;50(1):43-53.

882 20. Thornton C, Anthony SJ. Quasi-static deformation of particulate media. Phil Trans R Soc Lond A. 1998;356:2763-82.

884 21. da Cruz F, Emam S, Prochnow M, Roux J noël, Chevoir F. Rheophysics of dense granular materials : Discrete simulation of plane shear flows. Phys Rev. 2005;72(February):1-17.

887 22. Shire T, Hanley KJ, Stratford K. DEM simulations of polydisperse media: efficient contact detection applied to investigate the quasi-static limit. Comput Part Mech [Internet]. 2020;8(4). Available from: https://doi.org/10.1007/s40571-020-00361-2

890 23. Zuo L, Baudet BA. Determination of the transitional fines content of sand-non plastic fines mixtures. Soils Found [Internet]. 2015;55(1):213-9. Available from: http://dx.doi.org/10.1016/j.sandf.2014.12.017

893 24. Kuerbis RH. Effect of gradation and fines content on the undrained response of sand. Department of Civil Engineering, The University of British Columbia. 1989.

895 25. Lade P V., Yamamuro JA. Effects of nonplastic fines on static liquefaction of sands. Can Geotech J. 1997;34(6):918-28.

897 26. Li W, Chu Y, Deng G, Cai H, Xie D, Lee Lee M. Study of shear induced stress redistribution in gap-graded soils by discrete element method. Comput Geotech [Internet]. 2023;156(June 2022):105248. Available from: https://doi.org/10.1016/j.compgeo.2023.105248
901 27. Sufian A, Artigaut M, Shire T, O'Sullivan C. Influence of Fabric on Stress Distribution in Gap-Graded Soil. J Geotech Geoenvironmental Eng. 2021;147(5):1-14.

902 28. Vallejo LE. Interpretation of the limits in shear strength in binary granular mixtures. Can Geotech J. 2001;38(5):1097-104.

905 29. Minh NH, Cheng YP, Thornton C. Strong force networks in granular mixtures. Granul Matter. 2014;16(1):69-78.

907 30. Li Y, Otsubo M, Ghaemi A, Dutta TT, Kuwano R. Transition of gap-graded soil fabric - shear wave measurements and dispersion relation. Soils Found [Internet]. 2022;62(1):101092. Available from: https://doi.org/10.1016/j.sandf.2021.101092

910 31. Cubrinovski M, Ishihara K. Maximum and minimum void ratio characteristics of sands. Soils Found. 2002;42(6):65-78.

912 32. Wang T, Liu S, Wautier A, Nicot F. Updated skeleton void ratio for gravelly sand mixtures considering the effect of grain size distribution. Can Geotech J. 2022;59(11):12-23.

915 33. Shire T, O'Sullivan C, Hanley KJ, Fannin RJ. Fabric and Effective Stress Distribution in Internally Unstable Soils. J Geotech Geoenvironmental Eng [Internet]. 2014;140(12):04014072. Available from: http://ascelibrary.org/doi/10.1061/%28ASCE%29GT.1943-5606.0001184

919 34. Evans DM, Zhou S. Liquefaction Behaviour of Sand-Gravel Composites. J Geotech Geoenvironmental Eng. 1995;121(3):287-98.

921 35. Sarkar D, Goudarzy M, König D, Wichtmann T. Influence of particle shape and size on the threshold fines content and the limit index void ratios of sands containing non-plastic fines. Soils Found. 2020;60(3):621-33.

924 36. Yin ZY, Zhao J, Hicher PY. A micromechanics-based model for sand-silt mixtures. Int J Solids Struct [Internet]. 2014;51(6):1350-63. Available from: http://dx.doi.org/10.1016/j.ijsolstr.2013.12.027

927 37. Li X, Ma J, Lv L. Glass Beads Test with True Triaxial Stress Path Achieved by Conventional Triaxial Apparatus. Fractal Fract. 2024;8(4).

929 38. Lade P V., Wang Q. Analysis of Shear Banding in True Triaxial Tests on Sand. J Eng Mech. 2001;127(8):762-8.

931 39. Lade P V., Duncan JM. Cubical Triaxial Tests on Cohesionless Soil. ASCE J Soil Mech Found Div. 1973;99(SM10):793-781.

933 40. Wang Q, Lade P V. Shear Banding in True Triaxial Tests and Its Effect on Failure in Sand. J Eng Mech. 2001;127(8):754-61.

935 41. Li K, Li X, Chen Q, Nimbalkar S. Laboratory Analyses of Noncoaxiality and Anisotropy of Spherical Granular Media under True Triaxial State. Int J Geomech. 2023;23(9):1-13.

938 42. Alshibli K, Williams H. A True Triaxial Apparatus for Soil Testing with Mixed Boundary Conditions. Geotech Test J. 2005;28(6):534-43.

940 43. Salgado R, Bandini P, Karim A. Shear Strength and Stifness of Silty Sand. J Geotech Geoenvironmental Eng. 2000;126(5):451-62.

942 44. Chang CS, Yin ZY. Micromechanical modeling for behavior of silty sand with influence of fine content. Int J Solids Struct [Internet]. 2011;48(19):2655-67. Available from: http://dx.doi.org/10.1016/j.ijsolstr.2011.05.014

945 45. Benahmed N, Nguyen TK, Hicher PY, Nicolas M. An experimental investigation into the effects of low plastic fines content on the behaviour of sand/silt mixtures. Eur J Environ Civ Eng. 2015;19(1):109-28.

948 46. Santamarina JC, Cho GC. Soil behaviour: The role of particle shape. In: Advances in Geotechnical Engineering: The Skempton Conference - Proceedings of a Three Day Conference on Advances in Geotechnical Engineering, organised by the Institution of
951 Civil Engineers. 2004. p. 604–17.

952 47. Adesina P, O'Sullivan C, Wang T. Understanding the interplay between particle shape, grading and sample density on the behaviour of granular assemblies: A DEM approach. Granul Matter [Internet]. 2024;26(1). Available from: https://doi.org/10.1007/s10035-023-01383-2

955 48. Bolton MD. The strength and dilatancy of sands. Geotechnique. 1986;36(1):65–78.

957 49. Vaid YP, Sasitharan S. The strength and dilatancy of sand. Can Geotech J. 1992;29:522–6.

959 50. Xiao Y, Xiang J, Liu H, Ma Q. Strength–dilatancy relation of sand containing non-plastic fines. Geotech Lett. 2017;7(2):204–10.

961 51. Wang T, Wautier A, Liu S, Nicot F. How fines content affects granular plasticity of under-filled binary mixtures. Acta Geotech [Internet]. 2022;17(6):2449–63. Available from: https://doi.org/10.1007/s11440-021-01430-3

964 52. Thornton C. Granular Dynamics, Contact Mechanics and Particle System Simulations: A DEM study. Granul Dyn Contact Mech Part Syst Simulations A DEM study. 2015;24:1–195.

967 53. Edwards SF. The equations of stress in a granular material. Phys A Stat Mech its Appl. 1998;249(1–4):226–31.

969 54. Liu D, Morimoto T, Carraro JAH, O'Sullivan C. A semi-empirical re-evaluation of the influence of state on elastic stiffness in granular materials. Granul Matter [Internet]. 2022;24(2). Available from: https://doi.org/10.1007/s10035-022-01215-9

972 55. Nie Z, Zhu Y, Zou J, Gong J, Liu S. DEM study of the microscopic characteristics and internal stability of binary mixtures. Powder Technol [Internet]. 2019;352:314–24. Available from: https://doi.org/10.1016/j.powtec.2019.04.077

975 56. Nicot F, Hadda N, Guessasma M, Fortin J, Millet O. On the definition of the stress tensor in granular media. Int J Solids Struct. 2013;50(14–15):2508–17.

977 57. Bagi K. Microstructural stress tensor of granular assemblies with volume forces. J Appl Mech Trans ASME. 1999;66(4):934–6.

979 58. Dafalias YF, Manzari MT. Simple Plasticity Sand Model Accounting for Fabric Change Effects. J Eng Mech. 2004;130(6):622–34.

981 59. Papadimitriou GA, Dafalias YF, Yoshimine M. Plasticity Modelling of the Effect of Sample Preparation Method on Sand Response. Soils Found. 2005;45(2):109–23.

983 60. Li XS, Dafalias YF. Anisotropic Critical State Theory: Role of Fabric. J Eng Mech. 2012;138(3):263–75.

985 61. Casagrande A. Characteristics of Cohesionless Soil Affecting the Stability of Slopes and Earth Fills. J. Boston Soc. Civ. Eng. 1925–1940, 257–276, reprinted in contribution to soil mechanics; 1936.

988 62. Satake M. Fabric tensor in granular materials. In: In: PA Vermeer & HJ Luger (eds) IUTAM Symposium on Deformations and Failure of Granular Materials 1982 Rotterdam:Balkema. 1982. p. 63–8.

991 63. Barreto D, O'Sullivan C, Zdravkovic L. Quantifying the evolution of soil fabric under different stress paths. AIP Conf Proc. 2009;1145:181–4.

993 64. Adesina PA, O'Sullivan C, Wang T. DEM study on the effect of particle shape on the shear behaviour of granular materials. Comput Part Mech [Internet]. 2023; Available from: https://doi.org/10.1007/s40571-023-00632-8

996 65. Maeda K, Sakai H, Kondo A, Yamaguchi T, Fukuma M, Nukudani E. Stress-chain based micromechanics of sand with grain shape effect. Granul Matter. 2010;12(5):499–505.

999 66. Rothenburg L, Bathurst RJ. Micromechanical features of granular assemblies with planar elliptical particles. Geotechnique. 1992;42(1):79–95.
1001 67. Adesina P, O'Sullivan C, Wautier A. The influence of particle size ratio on the
1002 mechanical behaviour of gap-graded granular assemblies under drained triaxial
1003 compression: A DEM study. Comput Geotech [Internet]. 2025;179(August
1004 2024):106987. Available from: https://doi.org/10.1016/j.compgeo.2024.106987

1005 68. Wu W, Kolymbas D. On some issues in triaxial extension tests. Geotech Test J.
1006 1991;14(3):276-87.

1007 69. Thevanayagam S. Intergrain contact density indices for granular mixes - I: Framework.
1008 Earthq Eng Eng Vib. 2007;6(2):123-34.

1009 70. Ni Q, Tan TS, Dasari GR, Hight DW. Contribution of fines to the compressive
1010 strength of mixed soils. Geotechnique. 2004;54(9):561-9.

1011 71. Adesina P, O'Sullivan C, Morimoto T, Otsubo M. Determining a representative
1012 element volume for DEM simulations of samples with non-circular particles.
1013 Particuology [Internet]. 2022;68:29-43. Available from:
1014 https://doi.org/10.1016/j.partic.2021.10.007

1015 72. Senetakis K, Coop MR, Todisco MC. Tangential load-deflection behaviour at the
1016 contacts of soil particles. Geotech Lett. 2013;3(2):59-66.

1017 73. Huang X, Hanley K, O'Sullivan C, Kwok CY. Exploring the influence of interparticle
1018 friction on critical state behaviour using DEM. Int J Numer Anal Methods Geomech
1019 [Internet]. 2014;38:1276-97. Available from:
1020 http://onlinelibrary.wiley.com/doi/10.1002/nag.527/abstract%5Cnhttps://web.natur.cun
1021 i.cz/uhigug/masin/download/MTVC_IJNAMG06-pp.pdf

1022