A semi-empirical re-evaluation of the influence of state on elastic stiffness in granular materials

This study uses data acquired from three-dimensional discrete element method simulations to reconsider what measure of state can be used to predict stiffness in granular materials. A range of specimens with linear and gap-graded particle size distributions are considered and stiffness is measured using small amplitude strain probes. Analysis of the data firstly confirms that the void ratio, which is typically used as a measure of state in experimental soil mechanics, does not correlate well with shear stiffness. However, the empirical expressions developed by Hardin and his colleagues can capture variations in stiffness, provided an appropriate state variable is used. The study then highlights that the contribution of individual contacts to the overall stiffness is highly variable, depending on both the contact force transmitted and the particle size. Analyses explore how the stress transmission both within and between the different size fractions affects the overall stiffness. This heterogeneity in stiffness relates to the heterogeneity in the stress transmission amongst the different fractions. By considering the heterogeneity of stress distribution amongst different particle size fractions, a new semi-empirical stress-based state variable is proposed that provides insight into the factors that influence stiffness.


Introduction
The particle size distribution ( PSD ), is a basic characteristic of a granular material and is known to strongly affect its overall mechanical behaviour. The small-strain or elastic stiffness, G 0 , is important in geotechnical engineering design and in the interpretation of geophysical seismic test data. Based on studies of clean, natural quartz sands, Hardin & Richart [6] proposed the following empirical expression to describe the relationship between G 0 , the global void ratio, e , and the mean effective stress, p ′ : where the function F(e) quantifies how the density state of the material influences stiffness independently of the mean effective stress, and A and n are fitting parameters. For smooth elastic spheres, n can be taken as 0.33 following Hertzian theory [25]; the larger n values observed in experiments on sands may be attributed to finite surface roughness and non-spherical particle geometry [4,25]. Referring to Mitchell & Soga [21], a number of expressions for F(e) have been proposed in the literature. For both uniform and well-graded sands, experimental data support use of the form proposed by Hardin & Richart [6]: where c is an empirical constant which depends upon the particle characteristics. Based upon their experimental data, Hardin & Richart [6] suggested c = 2.97 for angular sands and c = 2.17 for rounded sands. In Eq. (2), e is taken as a measure of the material state. A number of researchers have highlighted a link between the PSD and the empirical (fitting) parameters adopted in Eqs. (1) and (2). Payan et al. [26], Menq [19], Wichtmann & Triantafyllidis [40] and Senetakis et al. [31] have all proposed relationships between A and the uniformity coefficient ( C u ).
It has long been accepted that in gap-graded soils finer particles may exist in the void space formed by the coarse fraction without transmitting stress and several variables have been proposed to describe the state of gap-graded soils.
(1) G 0 = AF(e)(p�) n (2) F(e) = (c − e) 2 1 + e DEM-generated data presented in Liu et al. [14] show that there can also be a finite number of inactive particles in materials that have continuous gradings with a high C u and trimodal materials comprising a mixture of three distinct size ranges. In these cases, it is difficult to support the use of e as a measure of state that can determine the mechanical behaviour. Mitchell (1976) assumed that the finer fraction does not contribute at all to stress transmission and proposed an intergranular void ratio, e g : where G f and G c are the specific gravities of finer and coarse grains, respectively. Given that G f and G c are equal, F finer is the proportion of the overall mass (or volume) that is made up of finer grains. More recently, researchers including Thevanayagam [36], Thevanayagam et al. [35] and Ni et al. [22] have argued that a valid measure of state should account for the proportion of the finer fraction that is active in stress transmission. Thevanayagam [36] therefore proposed an equivalent void ratio: where b denotes the proportion of the finer grains that actively engage in contacts. Some researchers (e.g. [5,30] proposed that b may depend upon F finer and the size ratio between the coarse particle diameter ( D c, 10 ) and finer particle diameter ( d f ,50 ). However, contrasting opinions are expressed in other studies (e.g. [42] around whether e * is an appropriate state variable because the contribution of the finer fraction to overall stress transmission cannot be measured experimentally. Recent experimental research has highlighted the challenges of directly applying Eqs. (1) and (2) to gap-graded soils. Wichtmann et al. [39] proposed correlations amongst c , n , A and F finer . Yang & Liu [41] highlighted the limitations of using Eq. (1) with gap-graded materials and proposed a relationship between A and F finer , i.e. A = 95.39e −F finer . Payan et al. [27] also observed a systematic variation in A with F finer (noting that they use a different functional form for F(e) ). Liu & Yang [16] indicated that curve fitting should be used to find A for a given F finer , and that at a given F finer , A depends on the difference in size between the coarser and finer fractions.
Working with discrete element method ( DEM ) data, Otsubo [23] proposed a mechanical void ratio, e m , which only considers active grains: where V v is the volume of voids, V ina is the volume of inactive particles (i.e. particles with no or 1 contact which do not transmit stress), and V a is the volume of the remaining "active" particles. The definition of e m is similar to that of e g (Eq. (3) above) in that both state variables e m and e g consider inactive particles as part of the void space. Considering a number of gap-graded DEM specimens, Otsubo et al. [24] found a correlation between e m and G 0 that holds over a range of F finer values, suggesting that e m may be a state variable that better correlates with G 0 , compared to e. Theoretical approaches have applied to derive general expressions for G 0 . Chang & Liao [2] proposed an Effective Medium Theory ( EMT ), which developed a relationship between G 0 and both the granular material properties and the stress level by assuming a uniform strain field (kinematic assumption). Also adopting the kinematic assumption proposed by Chang & Liao [2], Otsubo [23] proposed that the stiffness of a granular material is strongly associated with Z∕(1 + e) so that: where v p and G p are the particle Poisson's ratio and particle shear modulus, respectively. The coordination number Z , is the average number of contacts per particle. Otsubo et al. [24] found that for bimodal soil mixtures, the shear wave velocity, V s correlates more strongly with the ratio Z m ∕(1 + e m ) rather than Z∕(1 + e) . Equation (6) is based upon a mean response of the material and does not consider the details of the distribution of force amongst the different contacts. Following Hertzian contact mechanics, the stiffness of individual contacts depends upon the contact force.
In Eqs. (5) and (6), for gap-graded soils both active coarse and finer particles are considered equivalent in the calculation of Z m and e m . However recent studies [14,33,34]) have shown that active coarse and finer particles do not make equivalent contributions to stress transmission. Besides, other communities have studied the granular mixtures. Voivret et al. [38] have studied highly polydisperse packings and stated that the strong force chains tend to pass through the larger particles, with the smaller particles acting to support these strong force chains. Jongchansitto et al. [11] highlighted the significant effect of particle size and particle number ratios on stress transmission for granular composites. Shaebani et al. [32] also showed a strong relationship between microscale quantities and their macroscale counterparts for granular assemblies. Mean packing properties such as average coordination number may be strongly associated with its stiffness tensor.
Reflecting on the lack of consensus amongst prior research studies, and drawing on recent research into stress transmission in granular mixtures, this contribution aims to reassess what state variable can be used to determine G 0 . In contrast to previous studies, a range of gradations including both continuous and gap-graded specimens are considered. The paper firstly describes the DEM simulation approach adopted, prior to describing how the data were analysed and demonstrating how regression analyses were used to explore which state variable determines G 0 . Then the fundamental mechanisms that underlie the link between the PSD shape and G 0 are explored by quantifying the contribution of individual contacts to G 0 .

Numerical Simulation Approach
The DEM simulations were conducted on cubic specimens using a modified version of the open-source molecular dynamics code Granular LAMMPS [28]. Periodic boundary conditions were employed to reduce boundary effects (e.g. [33,37]. and so no gravitational body force was applied. A simplified Hertz-Mindlin contact model was used with a particle shear modulus of G p = 29.17 GPa and a particle Poisson's ratio of p = 0.2 following Huang et al. [8]. These parameters are similar to the elastic properties of quartz and have been used in a number of prior DEM studies (e.g. [8]. The grains were initially placed randomly in non-contacting positions using an in-house placement code, and then subjected to periodic isotropic compression to a target mean effective stress (p′) following [3]. Similar to the data presented in [14], specimen sizes of up to 636,871 particles were used to ensure representative element volumes were attained. Three p′ values (i.e. 100 kPa, 500 kPa, 1000 kPa) were applied. Once the target p′ value was attained, a relaxation stage was performed in order to reach a stable, quasi-static condition. It was confirmed that the coordination number, Z (the average value of contacts per grain), remained stable (i.e. any variation was less than 0.001) for approximately 500,000 simulation cycles prior to terminating the compression stage of the simulations. For all simulations presented, the unbalanced force ratio was lower than 0.001 at this stage [15].
Following Shire et al. [33], to investigate density effects, three inter-particle friction coefficients ( ) were used during isotropic compression: (1) "dense" specimens were generated using = 0.001; (2) "medium-dense specimens" were generated using = 0.1; (3) "loose" specimens were generated using = 0.3. This approach to specimen generation allows for a consistent generation of different packing configurations. However, the loose and dense specimens generated are not directly equivalent to physical laboratory specimens with relative densities of 0% and 100% respectively. Following isotropic compression, was set to be 0.3 and additional equilibration cycles were applied. Triaxial compression simulations with a constant radial stress and small increments in the axial strain were then carried out to explore the elastic behaviour of granular materials. While the simulation approach considers a highly ideal scenario, it allows for observation of the fundamental link between G 0 and the PSD as particle are spherical and all specimens are subjected to an isotropic stress condition. Figure 1 summarizes the gradings considered. Some of these gradations were also considered in [14]. Five "linear" specimens whose PSDs plot as a straight line (on the semilog axes) with C u values ranging between 1.2 and 4.8 are presented in Fig. 1a. The minimum particle diameter, d min , used in each linear specimen was 0.076 mm(i.e. a fine sand), the maximum particle diameter increased with increasing C u . For the bimodal gap-graded specimens (Fig. 1b), SR (i.e. D c,50 ∕d f ,50 ) values of 3.7, 8.4, 14.5 and 18.1 were considered, where D c,50 and d f ,50 are the median diameters of the coarse and finer fractions, respectively. The d min values were also 0.076 mm so that the maximum diameters, D max , increase with increasing SR. F finer values of 5%, 10%, 15%, 20%, 25%, 30%, 35% and 50% were considered. The SR and F finer are used to identify each specimen, for example, SR14.5F finer 50 indicates a bimodal specimen with SR = 14.5 and F finer = 50%. Trimodal specimens were considered in the study; these specimens are characterised by a F finer (by mass), a F coarser and an additional F int (Fig. 1(c)). Each trimodal specimen had a uniform finer particle diameter of 0.076 mm , an intermediate particle diameter of 0.425 mm and a coarse particle diameter of 0.985 mm . These specimens are named TriA_B , where A represents the specimen F finer and B represents the specimen F int so that Tri10_60 indicates a trimodal specimen with F finer = 10% and a F int = 60%. Overall, the study considered 5 linear PSDs , 32 bimodal gap-graded PSDs , and 22 trimodal gap-graded PSDs.
The e , e m , Z and Z m data for the linear specimens are presented in Fig. 2. As illustrated in Fig. 2a, similar to the observation by Youd [44], who considered experimentally measured maximum and minimum void ratios, the void ratios of the loose and dense specimens decrease with increasing C u , and a similar trend can be observed for e m . The offsets between the e and e m data show that even for specimens with a continuous grading, a finite proportion of the particles may not be active in stress transmission; the offset is greater for the loose specimens and increases with increasing C u . Figure 2b indicates that Z decreases with increasing C u , while, for a given density level, Z m does not vary much with C u . Figure 3 presents the e , e m , Z and Z m data for the bimodal gap-graded specimens. The observed trends in the variation of e and e m with F finer differ for the case of SR = 3.7 in comparison to the data for SR of 8.4 and 18.1; this can be explained by the fact that only when the SR > 6.5 can the finer grains in a purely bimodal material fit inside the voids formed by the coarse grains. Referring to Fig. 3a, for SR = 3.7 both the loose and dense e values initially decrease with increasing F finer and achieve a global minimum when F finer ≈ 20%. When F finer ≥ 20% , both the loose and dense e values exhibit a slight increase with increasing F finer . The  Fig. 3b, c show that the loose and dense e values initially reduce with increasing F finer to attain a global minimum and then subsequently increase with increasing F finer . These data agree with prior experimental observations (e.g. [45]. In contrast, the e m values initially increase with increasing F finer before exhibiting a sharp drop indicating a marked change of soil fabric; subsequently e m increases with increasing F finer . The variation in Z with F finer for SR = 3.7 is distinct from the trends observed for the other SR values. For SR = 3.7, Z decreases to attain a local minimum when F finer ≈ 15% , and subsequently increases with further increases in F finer (Fig. 3d). However, for a given packing density level, there is little variation in Z m with F finer . Figure 3e, f show that when SR ≥ 8.4 initially, the Z values reduce to a very low value, so that Z ≈ 0.05, indicating that a large number of inactive particles. Contrasting the Z values in Fig. 3e and f, with the e values in Fig. 3b and c, the change in Z is much steeper than the change in e once the finer grains start to engage in stress transmission. For all three SR values considered, the Z m values show little variation with F finer in either the dense or loose specimens.

Stiffness measured in small strain probes
The shear stiffness, G 0 , was measured from triaxial compression tests with a constant effective radial stress, ′ r , where a small increment in axial strain ( a ≤ 5 × 10 −6 ) was applied so that the specimen response could be considered isotropic and elastic. G 0 is calculated as: where q is the increment in deviatoric stress, and e q , the incremental elastic deviatoric strain is: e z = a is the applied strain in the z direction, e x and e y = r are the incremental strains in the x− and y− directions, respectively.
To illustrate how G 0 was determined, data for a representative specimen are presented in Fig. 4; the linear specimen considered has C u = 1.2, � 3 = 1000 kPa and is in the medium-dense density condition. Figure 4a illustrates the linear stress-strain behaviour observed for the applied small increment in axial strain ( a ≈ 5 × 10 −6 ). Figure 4b then shows the relationship between the deviatoric strain q and G 0 , confirming a constant G 0 for this increment in strain. Figure 5a shows the relationship between G 0 and e for the linear specimens at � 3 = 500 kPa . A state-dependent behaviour is evident; for a given C u , G 0 decreases systematically with increasing e . In line with the experimental research highlighted above, the data indicate that the relationship between G 0 and e is dependent upon C u ; the datasets for each C u are distinct and the data points shift to the left as C u increases. There is not a unique overall relationship between G 0 and e . In contrast, Fig. 5b shows a relatively good linear correlation between e m and G 0 , unifying the data for all 5 C u values considered. Figure 5c shows that, in line with EMT , a (weak) correlation exists when G 0 is plotted against Z∕(1 + e) . However, Fig. 5d shows that Z m ∕(1 + e m ) can better capture the variations in G 0 amongst the linear specimens. The same conclusion was drawn for the data obtained when ′ 3 = 100 kPa and 1000 kPa.   Figure 6 presents the variation in G 0 with e for three bimodal SRs (i.e. SRs of 3.7, 8.4, 18.1) when p′ is 500 kPa . As before the trends observed when SR = 3.7 differ from the cases where SR = 8.4 and 18.1. WhenSR ≥ 8.4 , the G 0 values vary with both e and F finer . As F finer increases, three distinct responses in G 0 can be observed. When 5% ≤ F finer ≤ 20% (data plotted in blue), the G 0 values exhibit similar sensitivity to changes in e and the data shift towards the left with increasingF finer . When20% ≤ F finer ≤ 35% , (data plotted in orange), the G 0 values are much more sensitive to changes ine . When F finer is approximately 50%, there is a noticeable reduction in the sensitivity of G 0 toe . These three categories identified ( 5 ≤ F finer ≤ 20%,20 ≤ F finer ≤ 35% and F finer > 35%) agree with the fabric classification boundaries proposed in Shire et al. [33] who considered the distribution of stress amongst the finer and coarser fractions. The agreement indicates that the distribution of stress amongst the different size fractions may have a significant influence on G 0 in the case of gap-graded soils. As before, similar patterns are observed for ′ 3 = 100 kPa and 1000 kPa. Figure 7 shows the performance of the four state variables considered in Fig. 5 for the bimodal specimens with SR of 8.4 at ′ 3 = 500 kPa. Figure 7a indicates no clear correlation between e and G 0 . In general, the G 0 − e relationship observed for both linear and bimodal gap-graded soils agrees with the prior DEM and experimental studies for granular materials. Figure 7b indicates some correlation between e m and G 0 , however the R 2 for this linear fit is approximately 0.80, this may be attributed to the inability of e m to differentiate contributions by particles or contacts transmitting different amounts of stress. Again, the correlation is poor where the original EMT ratio, i.e. Z∕(1 + e) is used (Fig. 7c). In contrast, the ratio Z m ∕(1 + e m ) gives a reasonable linear correlation, with a R 2 of 0.87. Figure 8 explores the applicability of Eq. (1) to bimodal gap-graded soils of SR of 8.4. In Fig. 8a the parameters A and c are identified by regression analysis using e , and taking n = 0.33; the correlation is very poor ( R 2 ≈ 0.46). However, Fig. 8b shows that if the ratio of Z m ∕(1 + e m ) is used in lieu of e in Eq. (2), (still assuming n = 0.33), a good fit to the data is attained, with R 2 ≈ 0.92. This result indicates that the ratio Z m ∕(1 + e m ) may be a better state variable to predict G 0 than e.

Bimodal specimens
Data for the trimodal specimens are presented in Fig. 9. The data on Fig. 9a confirm that for soil comprising a mixture of distinct particle sizes, e does not correlate well with G 0 . In contrast, Fig. 9b shows a linear correlation between Z m ∕(1 + e m ) and G 0 with R 2 ≈ 0.84. However, comparing the data on Fig. 5d with the data on Figs. 7d and 9b, the ratio Z m ∕(1 + e m ) gives a measurably stronger correlation    for the linear specimens than for the bimodal and trimodal specimens. Figure 10 combines data from all of the linear, bimodal and trimodal specimens considered. The poor predictive capacity of Eqs. (1) and (2) is observed when e is taken as the state variable (Fig. 10a). In contrast, the reasponable fit to the data can be attained by subsituting the ratio Z m ∕(1 + e m ) for e in Eq. (2). These data support the use of the functional forms presented by Hardin & Richart [6], however the controlled, ideal dataset generated here confirm that extrapolating the correlation of G 0 with e beyond the continuous PSDs that they considered is not appropriate for a unified study.

Heterogeneity of Contributions to Stiffness
If we accept the state variable that determines the stiffness of a soil is Z m ∕(1 + e m ) , we are not considering the fact, apart from the case of a lattice packing of uniform spheres, each contact has a different stiffness, and each particle transmits a different amount of stress. For the relatively simple smooth frictional spheres with uniform elastic properties that are considered here, the contact stiffness depends upon the force transmitted and the radii of the two transmitting particles [7,10]). A measure of state such is as Z m ∕(1 + e m ) assumes either (i) each active contact and each active particle has equal weight or (ii) that the average properties are representative. However, these assumptions may be appropriate for linear specimen and not for gap-graded soils. This could partially explain the reason why the ratio Z m ∕(1 + e m ) gives a measurably stronger correlation for the linear specimens than for the bimodal and trimodal specimens. Here we challenge the general applicability of any measure that considers Z m and e m by exploring how the particle size distribution influences the distribution of stress amongst the active particles and the distribution of stiffness amongst the force-transmitting contacts, especially for gap-graded soils.

Distribution of Contact Stiffnesses
A Hertzian contact model is used here and so the normal (incremental) stiffness ( k N ) depends on the diameters of the contacting particles and the contact force, as follows: where is the contact overlap and R * = (1∕R 1 + 1∕R 2 ) −1 , R 1 and R 2 are the radii of the contacting particles. The equivalent Young's modulus of the two contacting particles,E * p is: where E p1 and E p2 are the Young's moduli of the contacting particles, and v p1 and v p2 are the Poisson's ratios of the contacting particles. The tangential or shear contact stiffness is: where G * p , the equivalent shear modulus of the two contacting particles is: and G p1 and G p2 are the shear moduli of the contacting particles and the radius of the circular contact are a = √ R * . Figure 11a is a scatter plot of the normal stiffness, k N , values versus the normal force value for each contact point for the dense linear specimens. In all cases, the contact normal stiffness increases with increasing normal force, however there is a wide range of stiffness values and there is no clear difference amongst the linear specimens considered here. Figure 11b shows the cumulative distribution of the normal stiffnesses, while Fig. 11c shows the data normalised by the mean normal stiffness respectively. While no significant differences can be identified amongst these specimens, the range of k N values increases with increasing C u . Figure 12a-c are scatter plots of the k N values versus the normal force for each contact point for representative bimodal gap-graded specimens with SR = 8.4, and with F finer of 10%, 25% and 50% respectively. There is a clear difference between the three different contact types in these specimens. For all three F finer values considered, the maximum contact force between two coarse particles ( C − C ) is significantly larger than the maximum contact force between a coarser and finer particle ( C − F ) which in turn is larger than the maximum contact force between two finer particles ( F − F ). Each contact type exhibits a different relationship between normal force and normal stiffness, so that, at a given normal force, the C − C contacts have the highest normal stiffness and the F − F contacts have the lowest contact stiffness. This heterogeneity reflects a complex interdependency of force and stiffness; k n depends upon on the effective radius. However, assuming that there is no generation and no loss of contacts when a small increment in uniform strain is applied, the increase in force at the stiffer contacts will be larger. Figure 12d presents the cumulative distribution, by number of contacts, of the normal stiffness for the dense bimodal specimens with SR of 8.4; a distinct behaviour can be observed between the underfilled specimens (i.e. F finer of 5%, 10%, 15% and 20%) and overfilled conditions ( F finer of 25%, 30%, 35% and 50%). The significantly higher contact stiffness values when F finer < 20%, reflect the large stiffness values associated with the C − C contacts. At F finer = 20%, the cumulative distribution indicates a mix of small contact stiffnesses and large contact siffnesses. Then for F finer ≥ 20%, there is a more continuous distribution showing normal stiffness, while these finer related contacts show relatively small normal stiffness values in the overfilled case. Figure 12e presents the results for the loose condition with similar trends being identified. These data highlight the fact that contacts between particles of differenct sizes have measurable differences in stiffness; this has not always been considered in attempts to develop frameworks to predict G 0 .

Stiffness matrix approach to consider contributions of contacts to overall stiffness
While looking at the distribution of stiffness values amongst the contacts gives useful insight into the system, these values cannot be directly related to the overall shear stiffness and the orientation of the contacts relative to the overall deformation field is not considered. Contacts in a physical granular material form a relatively complex network and the properties of the network edges depend upon the contact orientations, the effective radii, and the transmitted contact forces. These characteristics cannot be directly captured considering only the contact stiffness values. To further advance our understanding we used a stiffness-matrix based approach to consider the contribution of each contact to the global stress field when a uniform strain field is assumed to describe the motion of all the particles. In this approach, following the procedure outlined in Itasca Consulting Group [9], a local stiffness matrix k ab is defined for the contact between particles a and b . For clarity, a representative 12 × 12 element stiffness matrix in the case where the contact normal is orientated along the x− axis is given here as: The increment in the contact force vector induced by an increment in displacement of the two grains is given by: where the incremental displacements of the two contacting particles are given by , and the vector f ab gives both the increments in forces and moments at the interaction between particles a and b.
Referring to Thornton [37], in a periodic cell, if the strain rate tensor applied to the periodic cell is given by ̇i j , then in a timestep Δt the incremental displacement of particle a due to the global strain field, Δx a i is Δx a i =̇i j x a i Δt , where x a i is the vector describing the location of particle a . In a similar manner, if an increment of strain Δ ij is applied to the system, the incremental displacement is given by: For simplicity, supposing that a deformation is small and assuming that the increment of branch vector Δl j ≈ 0 , ΔV ≈ 0 , and there are no newly generated contacts. If each contact experiences an incremental force, Δf c i , then the resulting increment in stress is given by: In engineering or Voigt notation the elasticity tensor D ij is given by: This is a symmetric tensor and only the upper right-hand side of the tensor is given here. Here the stiffness matrices for each contact were determined from the data generated in the DEM simulations. The particle displacements due a specified Δ ij were determined from Eq. (15) and the increment in force for each contact was determined using Eq. (14) before applying Eq. (17) to determine the increment in stress and hence elements of the stiffness matrix.
Based on the uniform strain approach assumption, the relative contributions of the different contact types to the stiffness can be estimated. Firstly, the overall stiffness D all ij is estimated, then the contribution of each of the networks can be isolated: CC contribution: For simplicity, supposing that a deformation is small and assuming that Δl j ≈ 0 , ΔV ≈ 0 , and no newly generated contacts: 13 Comparison between probe and matrix methods for linear specimens with different C u values: a G 0 values for these two methods when p′ is 500 kPa ; b Normalized stiffness of G 500kPa 0 ∕G 100kPa 0 For example, a strain field Δε ij > 0, Δε other = 0 is imposed to the system.
In this case, based on elastic theory, G 0 can be calculated as 0.25 × (D 13 + D 31 ) , where D 13 and D 31 are defined in Eq. (17). This stiffness matrix method therefore estimates the G 0 over the whole specimen. And the stiffness contribution of each fraction (i.e. coarse and finer fraction) can also be estimated. Figure 13 compares the stiffness values from the probe and the matrix method. For the dense specimens, as shown in Fig. 13a, the probe and the matrix methods give similar trends when the variation in G 0 with C u is considered. However, the stiffness values obtained from the matrix method consistently exceed those obtained from the probe method. This may be attributed to the uniform strain assumption underlying the stiffness matrix approach. Yimsiri & Soga [43] and Magnanimo et al. [17] amongst others have highlighted the inability of EMT to accurately capture the variation in stiffness amongst soil specimens. The limitations of this kinematic assumption [12,13,18] likely explain the overestimation of the stiffness values. However, upon normalization, the matrix approach predicts the same state-dependency stiffness as observed in the data acquired using the probe method. Referring to Fig. 13b, if the G 0 data obtained at 500 kPa are normalized by the data at 100 kPa , the responses are equivalent for both probe and matrix methods. Figure 14a shows the variation in the G 0 values from both probe and matrix methods for the bimodal specimens with SR = 8.4. As in the case of the linear specimens G matrix 0 values are consistently larger than the G 0 data obtained using the probe method. However, G 0 and G matrix 0 values exhibit similar trends. This is confirmed in Fig. 14b. Figure 14b shows that at a givenF finer , the G 0 values at 500 kPa normalized by the G 0 values at 100kPa , i.e. G 500kPa 0 ∕G 100kPa 0 , obtained using both methods are very similar. Both approaches predict very similar trends. Overall, while this matrix method overesti-matesG 0 , it can capture the variation in G 0 with state and the key factors that contribute to this overall stiffness.

Heterogeneity of stress amongst particles
To link the stiffness distribution to the particle-scale stress distribution, two approaches were used to quantify the distribution of stress amongst the different particle sizes in these specimens. In the first, particle-based approach, the overall stress tensor can be determined from the stress tensors for the individual particles as: where ′ ij is the effective overall stress tensor for the specimen, N t is the total number of stress-transmitting particles, a ij is the average stress tensor within particle a , and V a is the volume of particle a (e.g. [29]. The contribution of an individual particle, a , to the overall stress, ̂ a ij is then given by: In the second, contact-based approach, the contact force data were used to calculate the stresses in the virtual specimen so that the 3D stress tensor is given by: where N c,V is the total number of contacts in volume V , f c i is the force vector for contact c and l c j is the branch vector, i.e. the vector joining the centroids of the particles which contact at c (e.g. [1]. The contribution of contact c to the overall stress,̂ c ij is: Figure 15a shows the PSDs for all the linear specimens considered, where the particle diameter is normalized by the maximum particle size, i.e. d p /d p M . By adopting Eq. (22), the cumulative distribution of each particle to the macro-scale mean effective stress ( ′ ij ) is presented in Fig. 15b, defined as PSD stress . By adopting Eq. (23), the cumulative distribution of the contribution of each contact to macro-scale mean effective stress ( ′ ij ), CSD stress is presented in Fig. 15c. The x-axis is the branch vector length normalised by the can also be isolated. Figure 17a and b show that for the SR of 3.7, the contribution of the F − F contacts to the overall stiffness is negligble when F finer ≤ 15% , and subsequently increases with increasing F finer , while the contribution of the C − C contacts to the overall stiffness decreases with increasing F finer . In contrast, the contribution from the C − F contacts increases significantly and reaches approxmately 65% when F finer is 50%. For the SR of 8.4 (Fig. 17c, d), a different behaviour is observed when F finer ≤ 20% , where both the F − F and C − F contacts have a negligible contribution to the overall stiffness at all packing densities. Figure 17c and d also indicates that the C − F have a contribution of 80% to the overall stiffness when F finer is approximately 50%. These stiffness distribution data clearly indicate that different classes of contact make a distinct contribution to the stiffness, which again supports the idea that particle scale stress distribution may significantly affect the macro scale stiffness of gap-graded soils. These data put into doubt whether Z m and e m are relative measures of state to predict the behaviour of mixtures, since these two parameters consider all active contacts and grains as to be equivalent.  Figure 18a and c shows the distribution in the proportion of stress transmitted by these different contact classes for sample Tri25_15 and Tri25_45 , respectively. While Fig. 18b and d illustrates the distribution of the contribution to the overall stiffness amongst the different contact types for these two samples. Differences between CSD stress and CSD stiffness can be observed. For instance, for Tri25_15 in the dense condition, the C − F contacts transfer approximately 34% stress while the contribution of the C − F contacts to the overall stiffness is approximately 66%. A significant density effect is also observed for both the Tri25_15 and Tri25_45 specimens in considering both the stress and stiffness data. This strong density effect is similar to that observed for bimodal gapgraded soils, which may be one of the key characteristics for gap-graded soils.

Measure of state that accounts for heterogeneity in stress transmission
As discussed above, active finer and coarse particles are given equal weighting in the calculation of e m , while referring to Fig. 12, the C − C contacts are much stiffer. Accepting that the general framework proposed by Hardin & Richart [6] is valid, we empirically explored a measure of state in which the contribution of different contacts or particles is weighted. Identifying a contact-based measure of state is nontrivial, contacts are not associated with a volume.
Consequently, we focussed on developing a rational measure of state that weights the contribution of the particles by the stress they transmit by adopting the stress reduction factor proposed in Shire et al. [33]. Following Shire et al. [33], this factor quantifies the stress contribution amongst the different particle size fraction and it is used here to weight the contribution of the finer or coarse fraction in a modified void ratio, we term e . When ≤ 1 the relative contribution of finer fraction is lower than its volume fraction, the stress reduction factor of finer fraction, f can be then estimated as below: 19 Variation in stress reduction factors with F finer for bimodal specimens: a SR of 3.7; b SR of 8.4 The p ′ finer is the effective stress transferred by the finer fraction, n p is specimen porosity, N p,finer is the number of active finer particles. The new measure of state is then given as: where the V t is the total volume of specimen, the V c is the volume of active coarse grains, the V f is the volume of active finer grains. While when f ≥ 1 , as illustrated in Shire et al. [33], the material is overfilled, the stress contribution of finer fraction is larger than its volume fraction. The stress contribution of coarse fraction is smaller than its volume fraction The c values are also presented in Fig. 19; for instance, in the dense condition ( SR = 8.4) , when F finer is 35% and 50%, the stress transferred by the coarse fraction is lower than its volume fraction, the c rather than f is then presented in Fig. 19b. Figure 20a illustrates the variation in G 0 with e for all the specimens with SR = 8.4 at the ′ 3 of 500 kPa , the correlation can be approximated relatively well by a straight line with R 2 ≈ 0.91. In contrast, the correlations were much weaker when e ( R 2 ≈ 0.27 ) and e m ( R 2 ≈ 0.80 ) were considered (subsets of the data on Fig. 7). Figure 21a shows that the relationship between e and Z m ∕(1 + e m ) is linear (for R 2 = 0.92) considering all the bimodal and trimodal specimens. For the trimodal specimens, the active coarse and intermediate fractions are combined to estimate e . Figure 21b also show that a reasonable fit data is obtained using e as the state variable in Eq. (2) and taking n = 0.33 . These results clearly indicate that the distribution of stress amongst the different size fractions has a significant influence on G 0 in the case of gap-graded soils.

Perspective
The key measures of state identified here, i.e. e m , Z m , e cannot be directly quantified in experiments. However, the data presented here have implications for experimental studies. Firstly, the current results clearly demonstrate use of the global void ratio, e , in Eq. (2) will not give accurate estimates and cannot be used to predict G 0 because the proportion of inactive particles may be large, particularly in the case of soil mixtures. The efficacy of e also decreases when packing density shifts towards looser conditions. To complete our analysis, we explore the potential to fit the parameters of Eq. (2) to see if it can be used for a particular PSD . Table 1 summarizes the results of our regression analyses taking n = 0.33 and c = 1.28 (following Otsubo & O'Sullivan [25]. In agreement with Yang et al. [42] and Liu & Yang [16], the fitted A values decrease with increasing F finer for F finer < 25% and A increases when F finer ≥ 25% . In contrast, when 25% ≤ F finer ≤ 35% , the curve fitting is relatively poor (e.g. R 2 ≈ 0.27 when F finer = 30%) . The fabric of soils could change from being underfilled to overfilled when 25% ≤ F finer ≤ 35% . Therefore, the data indicate that for F finer outside of this transitional region, e may be adopted in Eqs. (1) and (2) provided A is determined for the particular PSD under consideration.

Conclusions
This study has presented a DEM investigation into the appropriate state variable to predict the small-strain stiffness of granular materials. A total of 218 specimens were considered including specimens with linear gradings and bimodal and trimodal gap-graded specimens, comprising two and three size fractions respectively. For each gradation, three packing densities were considered. The DEM simulations exploited high performance computing capabilities to achieve representative element volumes and used periodic boundaries to ensure homogeneous specimens.
The distribution of contact stiffness amongst the contacts in the specimens was considered from a particle-scale perspective. A theoretical stiffness-matrix method is used to analyse the contribution of each contact to the macro-scale stiffness. Even though the study is highly ideal (gravity free environment, spherical particles and isotropic fabric and stress), the following main points can be made based upon the data presented: For specimens with a linear grading, there is a clear link between the particle size distribution and its stiffness distribution. The state variables e m and Z m ∕(1 + e m ) are capable of describing the variation in stiffness amongst these linear specimens with greater accuracy than simply using e. For bimodal and trimodal gap-graded specimens, the variations in stiffness cannot be captured by correlating the data with e or e m . However, the ratio Z m ∕(1 + e m ) , which relies upon particle-scale data, can effectively describe the variations in G 0 amongst the specimens considered here.
A new stiffness-matrix approach was proposed to isolate the contributions of specified groups of particles to the overall stiffness. Based on this method and consideration of the distribution of contact stiffnesses, we showed that the C − C , C − F and F − F contacts each make measurably different contributions to the overall material stiffness.
A new semi-empirical state variable e was proposed to account for the non-uniformity in stress and contact forces within granular materials. The distribution of stress amongst the different size fractions is expressed by means of stress reduction factors and in turn has a significant influence on G 0 in the case of gap-graded soils.
It is not possible to measure e m , Z m or e experimentally. Our data support an empirical approach that uses regression analyses to determine the parameters in Eq. (1) for the particular material under consideration.