A micromechanics-based model for sand-silt mixtures

Zhen-Yu Yin, Jidong Zhao, Pierre-Yves Hicher

► To cite this version:

Zhen-Yu Yin, Jidong Zhao, Pierre-Yves Hicher. A micromechanics-based model for sand-silt mixtures. International Journal of Solids and Structures, 2014, 51 (6), pp.1350-1363. 10.1016/j.ijsolstr.2013.12.027. hal-01004847

HAL Id: hal-01004847

https://hal.science/hal-01004847

Submitted on 8 Apr 2017

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.

Public Domain
# A micromechanics-based model for sand-silt mixtures

Zhen-Yu Yin<sup>a,b</sup>, Jidong Zhao<sup>c</sup>, Pierre-Yves Hicher<sup>a,d</sup>

<sup>a</sup>LUNAM University, Ecole Centrale de Nantes, UMR CNRS GeM, Nantes, France

<sup>b</sup>Department of Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

<sup>c</sup>Department of Civil and Environmental Engineering, Hong Kong University of Science and Technology, Hong Kong, China

<sup>d</sup>College of Civil Engineering, Tongji University, Shanghai, China

Experimental observations have shown the significant impact of fines or coarse grains on the behavior of sand-silt mixtures. To describe the behavior of sand-silt mixtures under both drained and undrained conditions, this paper presents a mathematical model based on a micromechanical approach. The novelty of this model is the introduction of the equivalent mean size and the evolution of the position of the critical state line with fines content for various sand-silt mixtures. The predictive capability of the model was evaluated by comparing the model simulations with experimental results on undrained triaxial tests of Foundry sand-silt mixtures with fines content, $f_c$ = 0–100% and Ottawa sand-silt mixtures with fines content $f_c$ = 0–50%, and on drained triaxial tests of Hong Kong Completely Decomposed Granite (HK-CDG) mixtures before and after erosion. The predicted local behavior in the contact planes has also been examined. It shows that all local contact planes are mobilized to different degrees in terms of local stress and strain and that a few active contact planes contribute dominantly to the deformation of the assembly, leading to an anisotropic global behavior when the soil is subjected to external loading.

**Keywords:**

Constitutive relations

Elasto-plasticity

Micromechanics

Sand

Silts

Critical state

## 1. Introduction

In recent years, the stress-strain behavior of sand-silt mixtures has been increasingly studied. Numerous experimental investigations (Pitman, Robertson, & Sego, 1994; Zlatovic & Ishihara, 1995; Lade & Yamamuro, 1997; Thevanayagam & Mohan, 2000; Salgado, Bandini, & Karim, 2000; Polito & Martin, 2001; Thevanayagam, Shenthan, Mohan, & Liang, 2002; Xenaki & Athanasopoulos, 2003; Naeini & Baziar, 2004; Yang, Sandven, & Grande, 2006; Murthy, Loukidis, Carraro, Prezzi, & Salgado, 2007; Papadopoulou & Tika, 2008; Belkhatir, Arab, Schanz, Missoum, & Della, 2011; etc.) have indicated that both fines and coarse grains have a significant effect on the physical properties of sand-silt mixtures such as index void ratios, relative density as well as on the mechanical behavior of the mixtures including instability, critical state, strength and stress-dilatancy. Effective theoretical approaches characterizing the behavior of sand-silt mixtures may provide useful information for many relevant problems such as the study of internal erosion in soils. To the knowledge of the authors, however, a unified approach has not yet been proposed.

From a micromechanical point of view, the different cases of fines content correspond indeed to microstructures with different

packing arrangements for a granular mix, which can be illustrated in Fig. 1 for cases ranging from coarse grains skeleton (sand) to coarse grains separated by fines and to fines skeleton (silt). Notably, for the case of coarse grains skeleton, the coarse grain contacts play a primary role, whereas with the increase of fines content, the fine grain contacts begin to play an increasingly important role as the coarse grains become dispersed in a system wherein they provide a secondary reinforcement effect. Furthermore, based on the references from Table 1, a transition zone can be seen before the soil mixture behavior is entirely governed by the fine grains (index void ratios, position of CSL around $f_c$ = 25–30% on various sand-silt mixtures). Approaches based on micromechanics appear to be a rational way of addressing the issues relevant to sand-silt mixtures. In this regard, the Discrete Element Method (DEM) has proved to be effective in modeling the mechanical behavior of a granular material and in capturing the microstructure information from the particle scale of the system (Guo & Zhao, 2013). However, at the current stage it remains difficult to apply DEM directly to the simulation of sand-silt mixtures, since adding fines may substantially increase the number of particles in a simulation, which may also lead to computational difficulties due to an excessive difference between the sand particles and fines. More recently, Chang and Yin (2011) have proposed a micromechanics-based model for silty sand in which the fines content is relatively low (less than 20%) and the soil behavior is primarily dominated by the coarse grain skeleton. However, the behavior of a general
<figure><img src="image_1.png" /><figcaption>Fig. 1. Sand-silt mixture packing with different fines content ($f_c$: the fines content).</figcaption></figure>

<table><caption>**Table 1**
Minimum void ratio and initial critical state void ratio with values of parameters for various sand-silt mixtures.</caption><thead><tr><th rowspan="2">Materials</th><th colspan="3">Coarse grains (sand)</th><th colspan="3">Fine grains (silt)</th><th colspan="3">For e<sub>min</sub> or e<sub>cr0</sub></th></tr><tr><th>D<sub>50</sub>/mm</th><th>e<sub>min</sub></th><th>e<sub>cr0</sub></th><th>d<sub>50</sub>/mm</th><th>e<sub>min</sub></th><th>e<sub>cr0</sub></th><th>a</th><th>m</th><th>f<sub>th</sub></th></tr></thead><tbody><tr><td>Foundry sand-silt<sup>a</sup></td><td>0.25</td><td>0.608</td><td>0.795</td><td>0.01</td><td>0.627</td><td>0.86</td><td>-0.4</td><td>0.5</td><td>0.25</td></tr><tr><td>Hokksund sand-silt<sup>b</sup></td><td>0.44</td><td>0.574</td><td>0.93</td><td>0.032</td><td>0.731</td><td>1.22</td><td>-0.2</td><td>0.45</td><td>0.25</td></tr><tr><td>Assyros sand-silt<sup>c</sup></td><td>0.3</td><td>0.594</td><td>0.805</td><td>0.02</td><td>0.664</td><td>0.77</td><td>-0.3</td><td>0.35</td><td>0.35</td></tr><tr><td>Monterey sand-silt<sup>d</sup></td><td>0.43</td><td>0.730</td><td>-</td><td>0.03</td><td>0.624</td><td>-</td><td>-0.6</td><td>0.55</td><td>0.23</td></tr><tr><td>Nevada 50/200 sand-silt<sup>e</sup></td><td>0.18</td><td>0.573</td><td>-</td><td>0.05</td><td>0.754</td><td>-</td><td>0.4</td><td>0.5</td><td>0.25</td></tr><tr><td>Ottawa 50/200 sand-silt<sup>e</sup></td><td>0.2</td><td>0.548</td><td>-</td><td>0.05</td><td>0.753</td><td>-</td><td>0</td><td>0.7</td><td>0.3</td></tr></tbody></table>

sand-silt mixture with higher fines content is important for different soils and applications. Currently, the model of Chang and Yin (2011) cannot be applied to predict the “transitional” mixture whose behavior is partially and entirely governed by the fine grains.

This paper therefore proposes a micromechanics-based model to predict the behavior of sand-silt mixtures following the previous work of Chang and Yin (2011), and is organized as follows. The micromechanical model for silty sand is first briefly described. The model is then extended to include the effects of both fines and coarse grains, through which the deformation of an assembly can be obtained by integrating the movement of the inter-particle contacts in all orientations. The model calibrations are based on experimental results available in the literature. Experimental results on Foundry sand-silt mixtures, Ottawa sand-silt mixtures and HK-CDG mixtures (Hong Kong Completely Decomposed Granite) before and after erosion are used for evaluating the model's capability for predicting the behavior of sand-silt mixtures under both drained and undrained conditions.

## 2. Micromechanics based model of silty sand

The model to be proposed in this paper is based on the micromechanics model of sand first proposed by Chang and Hicher (2005) and further developed by Chang and Yin (2010a), Chang and Yin (2010b, Chang and Yin (2011), Chang, Yin, and Hicher (2011),

Yin, Chang, and Hicher (2010) and Yin and Chang (2012) accounting for various features including anisotropy, instability, cyclic loading effects and fines effects. In this approach sand is accounted as an assembly of particles, and an inter-particle model defines the relationship between inter-particle forces and movements in a contact plane between two grains. The stress-strain relationship is then obtained by integrating the particle movements and inter-particle forces. Normally such a model needs to be established using a phenomenological approach, since the packing detail is not known a priori, and the particle sizes, shapes as well as their alignments can be randomly distributed. It is instructive to begin by a brief introduction of the features of this micromechanics approach. In this approach, the elastic stiffness for two particles is assumed to follow a revised Hertz-Mindlin formulation (Hicher & Chang, 2006),

$$k_n = k_{n0} \left( \frac{f_n}{f_{ref}} \right)^n; \quad k_r = k_{rR} k_{n0} \left( \frac{f_n}{f_{ref}} \right)^n \quad (1)$$

where $f_{ref}$ is the reference value, $f_{ref} = (3V/Nl)$, $p_{at}$, $p_{at}$ is atmosphere pressure equal to 101.3 kPa, $V/N$ is the total number of contacts per unit volume, $l$ is the branch length between the two particles, and $f_n$ is the contact force in normal direction. Three material constants $k_{n0}$, $k_{rR}$ and $n$ are needed as input. To describe the plastic inter-particle behavior, a Mohr-Coulomb type yield function is defined in the contact force space ($f_n$, $f_s$, $f_t$), where the subscripts $n$, $s$, and $t$ are the local coordinates on each contact as shown in Fig. 2).

<footnote>a Thevanayagam et al. (2002).</footnote>

<footnote>b Yang et al. (2006).</footnote>

<footnote>c Papadopoulou and Tika (2008).</footnote>

<footnote>d Polito and Martin (2001).</footnote>

<footnote>e Lade and Yamarama (1997).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 2. Local coordinate at an inter-particle contact in the micromechanics sand model.</figcaption></figure>

$$F(f_n, f_r, \kappa) = \frac{f_r}{f_n} - \kappa(\delta_r^p) \quad (2)$$

where $\kappa(\delta_r^p)$ is a hardening/softening function; the shear force $f_r = \sqrt{f_s^2 + f_t^2}$ and the plastic sliding $\delta_r^p = \sqrt{(\delta_s^p)^2 + (\delta_t^p)^2}$ are defined. The hardening function is defined by a hyperbolic curve in $\kappa-\delta_r^p$ plane:

$$\kappa = \frac{k_{p0} \tan \phi_p \delta_r^p}{f_n \tan \phi_p + k_{p0} \delta_r^p} \quad (3)$$

where the initial slope of the hyperbolic curve is assumed as $k_{p0} = k_{pR}k_n$, with $k_{pR}$ as input.

The plastic sliding at an inter-particle contact often occurs with an upward or downward movement, which is shear-induced dilation/contraction, expressed as follows,

$$\frac{d\delta_n^p}{d\delta_r^p} = D \left( \tan \phi_{pt} - \frac{f_r}{f_n} \right) \quad (4)$$

where D is the material constant for stress dilatancy and $\tan\phi_{pt}$ represents the obliquity at which the dilation is zero. This assertion of a coupling effect due to plastic sliding of two particles is not supported by direct observations on the microfabric of soils. It is rather an hypothesis made from the observed behavior of soil specimen.

To account for the interactions among neighboring particles, a density state, $e_c/e$ is introduced in reference to the void ratio at critical state, $e_c$ which is a function of the mean effective stress of packing $p'$:

$$e_c = e_{cr0} - \lambda \left( \frac{p'}{p_{at}} \right)^{\zeta} \quad (5)$$

where $e_{cr0}$ is the initial value of the critical state line; $\zeta$ and $\lambda$ are material constants (Li & Wang, 1998).

Based on experimental results on various silty sands, the location of the critical state line is significantly influenced by the fines content. Taking that into account, the initial value of the critical state line $e_{cr0}$ shown in Eq. (5) is given as a function of $f_c$ by:

$$e_{cr0} = e_{hc,cr0}(1 - f_c) + af_c \quad (6)$$

where $e_{hc,cr0}$ is the initial value of the critical state line for host sand. Thus, combining Eq. (6) with Eq. (5), the critical void ratio for silty sand with different fines contents can be obtained.

The density state is used to adjust the apparent friction angle in light of the experimental observation by Biarez and Hicher (1994):

$$\tan \phi_p = \frac{e_c}{e} \tan \phi_\mu \quad (7)$$

For dense packing ($e_c/e > 1$), the apparent inter-particle friction angle, $\phi_p$ is greater than the internal friction angle, $\phi_\mu$. When the packing structure dilates, the degree of interlocking and the

apparent frictional angle are reduced, which results in a strain-softening phenomenon. For loose packing ($e_c/e < 1$), the apparent frictional angle, $\phi_p$ is smaller than the internal friction angle, $\phi_\mu$.

The value of $\tan\phi_{pt}$ in Eq. (4) is also a function of density state in light of observations on the phase transformation line. It is given by:

$$\tan \phi_0 = \left(\frac{e_c}{e}\right)^{-1} \tan \phi_\mu \quad (8)$$

This relationship allows a dense packing to dilate at an earlier stage of shear loading (a similar approach using an exponential form can be found in Manzari and Dafalias (1997)) Li and Dafalias (2000).

The stress-strain relationship for an assembly can be determined from integrating the behaviour of inter-particle contacts in all orientations (Chang & Hicher 2005; Yin et al. 2010). Thus, a sufficiently large contact number $N$ for a representative packing volume $V$ is needed to sum up global stresses and strains from forces and displacements of micro planes. The total number of contacts per unit volume can be approximately related to the void ratio. Chang and Yin (2011) have treated silty sand as an inter-granular structure with an equivalent inter-grain void ratio $e_{eq0}$. The total number of contacts per unit volume is expressed as:

$$\frac{N}{V} = \frac{3C_n}{\pi d^3 (1 + e_{eq0})} \frac{(1 + e_{eq0})e_{eq0}}{(1 + e_{eq})e_{eq}} \quad (9)$$

where $d$ is the mean particle size of packing and $C_n$ is the average co-ordination number, which can be approximated for silty sand by (Chang & Misra 1990).

$$C_n = 13.28 - 8e_{eq0} \quad (10)$$

where $e_{eq}$ is expressed by the current void ratio and fines content, $f_c$ according to Chang and Yin (2011),

$$e_{eq} = \frac{e - af_c}{1 - f_c} \quad (11)$$

where $a$ is a material constant depending on the fabric structure of the soil mixture (e.g. grading, particle shape). The equation can be reduced to the skeleton void ratio defined by Vaid (1994) when $a = -1$, and implies that for $f_c = 0\%$, $e_{eq}$ converges to the void ratio of sand when micromechanical model extended to sand-silt mixtures

## 2.1. Equivalent inter-grain void ratio and contact density

For the case of a coarse grain skeleton, adding fines may render the mechanical behavior of such a mixture stronger than that of the host coarse grain soil. According to Chang and Yin (2011), the equivalent inter-granular contact index, $e_{eq}$ can be expressed by Eq. (11).
For the case of a fine grains skeleton, the reinforcement effect by the coarse grains was introduced by Thevanayagam et al. (2002) to obtain an equivalent inter-fine void ratio as the index of active contacts, given by:

$$e_{eq} = \frac{e}{f_c + \frac{1-f_c}{(R_d)^m}} \qquad (12)$$

where R_d is the ratio of the mean size of coarse grains, D_50 to the mean size of fine grains, d_50 and m is a coefficient depending on grain characteristics and fine grain packing (0 < m < 1). Eq. (12) implies that for f_c = 100%, e_{eq} converges to the void ratio of fine grains.

A transition zone can be seen before the soil mixture behavior is entirely governed by the fine grains (index void ratio in Fig. 3, position of CSL in Fig. 4, around f_c = 25-35%). By taking into account the transition zone from the coarse grains skeleton to the fine grains skeleton, a unified equivalent inter-grain void ratio can be defined by fines content by combining Eq. (11) with Eq. (12), and using the hyperbolic tangent function tanh x = (e^{2x} - 1)/(e^{2x} + 1):

$$e_{eq} = \frac{e - af_c}{1 - f_c} \frac{1 - \tanh[\xi(f_c - f_{th})]}{2} + \frac{e}{f_c + \frac{1 - f_c}{(R_d)^m}} \times \frac{1 + \tanh[\xi(f_c - f_{th})]}{2} \qquad (13)$$

where ξ is a material constant controlling the evolution rate of the transition zone (see details in the sequel); f_th is the threshold fines content from the coarse grains skeleton to the fine grains skeleton.

The total number of contacts per unit volume for sand-silt mixtures adopts a form similar to Eq. (9), but with a further modification of the mean particle size, d, by employing the following form of hyperbolic tangent function:

$$d = D_{50} \frac{1 - \tanh[\xi(f_c - f_{th})]}{2} + d_{50} \frac{1 + \tanh[\xi(f_c - f_{th})]}{2} \qquad (14)$$

Eq. (14) indicates that d varies with f_c. For f_c = 0%, d = D_50 (mean particle size of sand) while for f_c = 100%, d = d_50 (mean particle size of silt). Consequently, the equivalent inter-grain contact density, N/V can be obtained by combining Eqs. (9), (10), (13), and (14).

## 2.2. Unified critical void ratio

According to Chang and Yin (2011), the void ratio of silty sand, e can be expressed by the void ratio of the coarse grains, e_hc and the fines content, f_c, as:

$$e = e_{hc}(1 - f_c) + af_c \qquad (15)$$

For sandy silt, we adopt the form of Eq. (12) for e_{eq} and assume the void ratio of the soil mixture, e can be expressed by the void ratio of the fine grains, e_{hf} and the fines content, f_c:

$$e = e_{hf} \left( f_c + \frac{1 - f_c}{(R_d)^m} \right) \qquad (16)$$

If the transition zone from the coarse grains skeleton to the fine grains skeleton is taken into account, the void ratio of the soil mixture can be expressed by combining Eq. (15) and (16) for soil mixtures from sand to silt in an analogous manner to the current void ratio, e:

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

The above formulation can be further verified by comparing the minimum void ratios of the soil mixtures for various fines content, i.e., by replacing e by e_min, e_hc by e_hc,min and e_hf by e_hf,min. The minimum void ratio, for example, can be expressed as:

$$e_{min} = [e_{hc,min}(1 - f_c) + af_c] \frac{1 - \tanh[\xi(f_c - f_{th})]}{2} + e_{hf,min} \left( f_c + \frac{1 - f_c}{(R_d)^m} \right) \frac{1 + \tanh[\xi(f_c - f_{th})]}{2} \qquad (18)$$

Six sand-silt mixtures were selected to examine the minimum void ratio expressed by Eq. (18). The results are shown in Fig. 3 (circles representing data) and summarized in Table 1. Note that the material constants a, m and f_th can be directly measured from the results. Four different values of ξ were used for comparison. Based on curve fittings from Fig. 3, ξ = 20 appear to be generally acceptable for all cases of the sand-silt mixtures. Fig. 3 also shows

<figure><img src="image_3.png" /><figcaption>Fig. 3. Evolution of the minimum void ratio with fines content for different values of ξ.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 4. Critical state lines of three sand-silt mixtures: (a) Critical state lines in e-logp' plan, (b) Initial critical void ratio versus fines content, and (c) Normalized critical state lines.</figcaption></figure>

that the unified formulation of the void ratio is suitable for describing the index void ratio of the soil mixture with different fines content.

Similar to $e_{min}$, the fines content also influences significantly the location of the critical state line of the soil mixtures. Fig. 4(a) shows the critical state lines in the e-logp' plane measured from triaxial tests on Foundry sand-silt mixtures by Thevanayagam et al. (2002), on Hokksund sand-silt mixtures by Yang et al. (2006), and Assyros sand-silt mixtures by Papadopoulou and Tika (2008), respectively. Based on measurements of these CSLs, Eq. (5) is adopted for the location of the critical state line. All three parameters in Eq. (5) can be determined for soil mixtures with different fines content. The same values of $\zeta$ and $\lambda$ can be assumed for each soil mixtures with almost the same minerals for sand and silt. Note that the use of the same values of $\zeta$ and $\lambda$ is not always necessary but practically convenient. Consequently, only the initial critical state void ratio, $e_{cr0}$ changes with fines content. Fig. 4(b) shows the variation of the measured values of $e_{cr0}$ with $f_c$ for all sand-silt mixtures. It is shown that the evolution trend and the threshold fines content for $e_{cr0}$ are similar to those for $e_{min}$. We hereby adopt a similar form of Eq. (18) for $e_{cr0}$,

$$e_{cr0} = [e_{hc,cr0}(1 - f_c) + af_c] \frac{1 - \tanh[\zeta(f_c - f_{th})]}{2} + e_{hf,cr0} \left(f_c + \frac{1 - f_c}{(R_d)^m}\right) \frac{1 + \tanh[\zeta(f_c - f_{th})]}{2} \quad (19)$$

where $e_{hc,cr0}$ and $e_{hf,cr0}$ are initial critical state void ratios for coarse grains and fine grains, respectively.

It is worth noting that the values of $a$ and $m$ for the minimum void ratio can be adopted for $e_{cr0}$ as well, based on the selected soil mixtures. However, the values of $a$ and $m$ are not necessarily the same for these two particular aspects of the material behavior. This result can be regarded as broadening the assumption that $e_{max}$ and $e_{min}$ are specific points of the critical state line for specific values of $p'$ as proposed by Biarez and Hicher (1994).

The $e_{cr0}$ obtained by Eq. (19) is used to normalize the critical state line with different fines contents. All normalized curves in Fig. 4(c) demonstrate that Eq. (19), with the same set of values of $a$ and $m$ as the minimum void ratios, is appropriate.

## 2.3. Further remarks on the model and summary of model parameters

All above mentioned new void ratios are summarized in Table 2. Arithmetic examples were given to further explain these definitions based on Foundry sand-silt mixture ($a = -0.4$, $m = 0.5$, $f_{th} = 0.25$, $R_d = 25$, $\xi = 20$, in Table 1): (1) for a given current void ratio $e = 0.65$, according to Eq. (13) the equivalent inter-grain void ratio is $e_{eq} = 0.843$ for silty sand with $f_c = 0.1$ and $e_{eq} = 0.707$ for sandy silt with $f_c = 0.9$; (2) for a given $e = 0.65$, according to Eq. (17) the void ratio of the coarse grains is $e_{hc} = 0.792$ for silty sand with $f_c = 0.1$ (the influence of second term with $e_{hf}$ is negligible)
<table><caption>**Table 2**
Definitions of new void ratios.</caption><thead><tr><td>Void ratio</td><td>Definition</td><td>Equation number</td></tr></thead><tbody><tr><td>e<sub>eq</sub></td><td>Equivalent inter-grain void ratio</td><td>(13)</td></tr><tr><td>e<sub>hc</sub></td><td>Current void ratio for pure coarse grains</td><td>(17)</td></tr><tr><td>e<sub>hf</sub></td><td>Current void ratio for pure fine grains</td><td>(17)</td></tr><tr><td>e<sub>hc,min</sub></td><td>Minimum void ratio for pure coarse grains</td><td>(18)</td></tr><tr><td>e<sub>hf,min</sub></td><td>Minimum void ratio for pure fine grains</td><td>(18)</td></tr><tr><td>e<sub>hc,cro</sub></td><td>Initial critical state void ratio for pure coarse grains</td><td>(19)</td></tr><tr><td>e<sub>hf,cro</sub></td><td>Initial critical state void ratio for pure fine grains</td><td>(19)</td></tr></tbody></table>

and void ratio of the fine grains is $e_{hf} = 0.707$ for sandy silt with $f_c = 0.9$ (the influence of first term with $e_{hc}$ is negligible); (3) for the coarse grains with minimum void ratio $e_{hc,min} = 0.608$ and the fine grains with $e_{hf,min} = 0.627$, their mixture with different $f_c$ can have different $e_{min}$ shown in Fig. 4(1-b) according to Eq. (18); (4) similar to (3), for the coarse grains with initial critical state void ratio $e_{hc,cro} = 0.795$ and the fine grains with $e_{hf,cro} = 0.86$, their mixture with different $f_c$ can have different $e_{cro}$ shown in Fig. 4(1-b) according to Eq. (19).

Since pure sand and pure silt can have different natures (mineral, shape, roughness, etc.), the fines content can also influence other properties (friction angle, small strain stiffness, slope of CSL, etc.) of the soil mixture. For instance, for friction angle $\phi_{\mu}$, the value of $\phi_{\mu}$ for silty sand is influenced by the fines content (Salgado et al. 2000). If the friction angles for pure sand and pure silt are identical or almost identical, a unique value of $\phi_{\mu}$ can be retained for all the mixtures. If this is not the case, we could try to suggest a continuous evolution of the friction angle from sand to silt. A linear or nonlinear evolution of $\phi_{\mu}$ from $\phi_{\mu,sand}$ to $\phi_{\mu,silt}$ with the increasing amount of fines, as shown in previous sections, could be established based on experimental measurements of the studied mix. This could be an increasing or a decreasing function according to the relative values of these two friction angles. At the present time there is a lack of experimental data able to back this idea, so this remains hypothetical. Note that in this paper we consider a more general case that pure sand and pure silt have almost identical natures.

In summary, the proposed model has two CSL related parameters ($e_{hc,cro}$ and $e_{hf,cro}$) instead of one parameter ($e_{cro}$), and three additional parameters $a$, $m$ and $f_{th}$. All these additional parameters can be determined in a straightforward way (Figs. 3 and 4 and Table 1). Elastic and plastic parameters of the model can be easily determined from either drained or undrained isotropic consolidated triaxial tests on sand or silt specimens (Yin et al., 2010; Yin & Chang, 2012) (see Table 3).

### 3. Model simulations and discussion

Triaxial tests on three different sand-silt mixtures were selected to examine the model performance on different aspects: (1) general undrained behavior through undrained triaxial tests on Foundry sand-silt mixtures with various void ratios and fines content varying from 0 to 100%; (2) static liquefaction through undrained triaxial tests on Ottawa sand-silt mixtures with fines content varying from 0 to 50%; (3) general drained behavior through drained triaxial tests on HK-CDG mixtures before and after internal erosion. The grain size distribution for all selected materials is shown in Fig. 5.

#### 3.1. Undrained behavior of sand-silt mixtures

Undrained triaxial tests on sand-silt mixtures with different fines content (Foundry sand $f_c = 0\%$, non plastic crushed silica fines

$f_c = 100\%$, and soil mixtures $f_c = 7, 15, 25, 40, 60\%$) by Thevanayagam et al. (2002) were selected to validate the model performance. The selected tests include three series: Series A with an increasing inter-granular void ratio ($e_g = 0.598-0.947$) and a decreasing inter-fine void ratio ($e_f = 1.84-0.77$), Series B with the same inter-granular void ratios ($e_g = 0.8$) and Series C with the same inter-fine void ratio ($e_f = 0.9$). CSL related parameters ($e_{hc,cro}$, $e_{hf,cro}$, $\zeta$ and $\lambda$, see Fig. 4-1) and fines content related parameters ($a, m$ and $f_{th}$, see Fig. 3a) were calibrated by the experimental results. Other parameters were determined from the triaxial tests on pure sand. In the absence of isotropic compression tests and measurements pertaining to small strain stiffness, $k_{n0} = 70$ was determined from the initial slope of the stress-strain curve, with assumed values $k_{fR} = 0.5$ and $n = 0.5$; the internal friction angle $\phi_{\mu} = 30^{\circ}$ was obtained from triaxial tests; $k_p = 0.15$ and $D = 1$ were obtained by fitting the stress-strain and the stress path curves. All determined values of parameters (Table 2) were used to predict other undrained triaxial tests on Foundry sand-silt mixtures with various fines contents.

#### 3.1.1. Test series A

Figs. 6(a and b) show a comparison between the measured and predicted results for silty sand samples ($f_c = 0-25\%$). In general, the sample of host sand ($f_c = 0\%$) behaves as a typical dense material and the samples of silty sand with 15% and 25% fines behave as loose materials, even though the samples with 0% and 15% fines have the same void ratio and the sample with 25% fines has an even smaller void ratio. This contrast of behavior for the two samples indicates that the void ratio is not the appropriate density index for predicting the behavior for silty sand. It is however found that the inter-granular void ratios of the four samples increase ($e_g = 0.60, 0.71, 0.87$ and $0.947$) as the fines content increases ($f_c = 0, 7\%, 15\%$ and $25\%$). This correlation indicates that the inter-granular void ratio would reflect better the fabric structure, as suggested by Thevanayagam et al. (2002).

Figs. 6(c and d) compare the measured results with predictions for sandy silt samples ($f_c = 40-100\%$) with the sample of $f_c = 25\%$ as a reference. At different initial void ratios ($e_0 = 0.425-0.77$), the samples exhibit a decreasing inter-fine void ratio with the increase of fines content ($e_f = 1.84-0.77$). The pure silt specimen shows a dilative behavior, whereas the sandy silt specimens with 40% and 60% fines content exhibit a contractive behavior, even though the void ratio of the silt sample is bigger than that of the sandy silt samples. This observation further confirms that the inter-fine void ratio may serve as a better index than the void ratio for predicting the behavior of sandy silt. Notably, however, the predicted loading paths for high $e_0$ values deviate considerably from the experimental data.

#### 3.1.2. Test series B

Fig. 7 shows the model predictions in comparison with the experimental results for silty sand samples at nearly the same $e_g = 0.8$, corresponding to void ratio $e_0 = 0.8, 0.672$ and $0.53$ for $f_c = 0, 7\%,$ and $15\%$, respectively. The loose clean sand specimen
<table><caption>**Table 3**
Values of model parameters for selected sand-silt mixtures.</caption><thead><tr><th>Group</th><th>Parameters</th><th>Foundry sand-silt</th><th>Ottawa sand-silt</th><th>HK-CDG</th></tr></thead><tbody><tr><td rowspan="3">Elastic parameters</td><td>k<sub>n0</sub></td><td>70</td><td>70</td><td>30</td></tr><tr><td>k<sub>fR</sub></td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><td>n</td><td>0.5</td><td>0.5</td><td>1</td></tr><tr><td rowspan="3">Plastic parameters</td><td>k<sub>p</sub></td><td>0.15</td><td>0.12</td><td>0.18</td></tr><tr><td>D</td><td>1</td><td>1</td><td>1</td></tr><tr><td>φ<sub>μ</sub></td><td>30</td><td>32</td><td>40</td></tr><tr><td rowspan="4">CSL related parameters</td><td>e<sub>hc,cro</sub></td><td>0.795</td><td>0.805</td><td>0.65</td></tr><tr><td>e<sub>hf,cro</sub></td><td>0.86</td><td>1.03</td><td>0.8</td></tr><tr><td>ζ</td><td>0.68</td><td>0.196</td><td>0.45</td></tr><tr><td>λ</td><td>0.03</td><td>0.081</td><td>0.03</td></tr><tr><td rowspan="3">Fines content related parameters</td><td>a</td><td>-0.4</td><td>0</td><td>-0.5</td></tr><tr><td>m</td><td>0.55</td><td>0.7</td><td>0.23</td></tr><tr><td>f<sub>th</sub></td><td>0.28</td><td>0.3</td><td>0.27</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Grain size distribution of selected soil mixtures.</figcaption></figure>

(FC: 0%) reaches static liquefaction after undrained shearing and loses almost all its strength at around 20% axial strain. On the other hand, the two silty sand specimens (FC: 7%, FC: 15%) are dilative. In particular, the $f_c$ = 7% case shows an initial post-peak softening and a strength increase to hardening behavior, which clearly demonstrates the impact of the fines content on the overall behavior of the mixture. At the same inter-granular void ratio, an increase in the fines content tends to render the soil to be more dilative and to increase its undrained strength. According to this series of tests, the inter-granular void ratio alone may not be sufficient as an index for predicting the mechanical behavior of sand-silt mixtures.

### 3.1.3. Test series C

Fig. 8 shows comparisons between the experimental and numerical results for sandy silt samples at almost the same inter-fine void ratio $e_f$ = 0.9, corresponding to the void ratio $e_0$ = 0.373, 0.54 and 0.879 for $f_c$ = 40%, 60%, and 100%, respectively. As seen in the Figure, the model can reproduce reasonably well the secondary reinforcement effect by the coarse grains: (1) the specimen with $f_c$ = 40% is significantly stronger than the other two specimens (with $f_c$ = 60% and 100%) due to the significant reinforcement by the coarse grains; (2) the sandy silt specimen with $f_c$ = 60% is slightly stronger than the silt specimen with $f_c$ = 100% since the coarse grains are relatively lower and are dispersed in the sample to provide adequate secondary reinforcement, and the behavior of the mixture is primarily affected by the inter-fine

contacts. According to the tests of Series C, the inter-fine void ratio alone is not sufficient as an index for predicting the behavior of the sand-silt mixture.

### 3.2. Static liquefaction of sand-silt mixtures

Lade and Yamamuro (1997) performed undrained triaxial tests on Ottawa sand-silt mixtures with different fines content varying from 0% to 50% to evaluate the effect of fines on static liquefaction. Unlike tests on Foundry sand-silt mixtures, the selected test samples show a general increasing trend of relative density with the increase of the fines content (see Fig. 12, $D_r$ from 15% to 41%). The model parameters were calibrated as follows: the fines content related parameters ($a$, $m$ and $f_{th}$) were measured from the minimum void ratios, shown in Fig. 3(f); the CSL related parameters ($\zeta$, $\lambda$) by Murthy et al. (2007) for Ottawa clean sand were adopted, and others ($e_{hc,cro}$, $e_{hf,cro}$) were determined from all selected undrained tests; the friction angle $\phi_\mu$ was determined from stress paths of triaxial tests; the same values of elastic parameters of Foundry sand were adopted due to the lack of data; the plastic parameters ($G_p$, $D$) were determined by curve fitting on the triaxial test of pure sand ($f_c$ = 0%). All determined parameters summarized in Table 2 were used to predict other tests.

A comparison between the experimental and numerical results presented in Fig. 9 clearly shows that the presence of fines can greatly increase the potential of static liquefaction of soil mixtures with density increasing from $D_r$ = 15% to 41%, in which the inter-granular and inter-fine void ratios are not kept constant. Lade and Yamamuro (1997) have indicated that the fines appear to create a “particle structure” in the soil which makes the soil more vulnerable to liquefaction, even through the density of the soil increases from 15% up to 41%. This observation is inconsistent with the soil behavior commonly observed in sand, i.e., the soil should exhibit more dilatant behavior as density increases. The proposed micromechanical model accounts for the interactions among neighboring particles by considering the fines content-dependent CSL with its related density state effect. Thus it can accurately capture the effect of this liquefiable “particle structure” mentioned by Lade and Yamamuro (1997) and reproduce their experiments reasonably well. It is hereby clear that neither the relative density nor the void ratio may provide satisfactory indications of the liquefaction potential in sand-silt mixtures.

All these results (undrained behaviors and static liquefaction) also show that neither the inter-granular void ratio nor the inter-fine void ratio can be a relevant density index if one wants to analyze the behavior of sand-silt mixtures with $f_c$ varying from 0 to 100%. They demonstrate the interest of using the equivalent void ratio as proposed in this paper to unify all the mixtures from pure
<figure><img src="image_6.png" /></figure>

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

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

<figure><img src="image_9.png" /><figcaption>Fig. 6. Comparison between experimental results and simulations for Foundry Sand-silt mixtures from sand to silt: (a, c) effective stress path, and (b, d) deviatoric stress versus axial strain.</figcaption></figure>

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

<figure><img src="image_11.png" /><figcaption>Fig. 7. Comparison between experimental results and simulations for Foundry silty sand with an inter-granular void ratio of 0.8: (a) effective stress path, and (b) deviatoric stress versus axial strain.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Comparison between experimental results and simulations for Foundry sandy silt with an inter-fine void ratio of 0.9: (a) effective stress path, and (b) deviatoric stress versus axial strain.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Comparison between experimental results and simulations for Ottawa sand-silt mixtures: (a) effective stress path, and (b) deviatoric stress versus axial strain.</figcaption></figure>

sand to pure silt with a unique expression of the density index introduced in the proposed model.

### 3.3. Drained behavior of eroded sand-silt mixtures

Chang and Zhang (2011) developed a stress-controlled erosion apparatus to investigate the initiation and development of suffusion (internal erosion) under complex stress states and to study the effect of suffusion on the soil's stress-strain behavior. The experimental set-up consists of a triaxial cell, a pressurized water supply system, a soil collection system, and a water collection system. A hollow base pedestal is designed for the eroded soil and seepage water to pass through. Three series of erosion tests were carried out under three stress conditions: isotropic stress, triaxial compression (i.e., the major principal stress is along the seepage direction), and triaxial extension. After the proposed stress state had been applied to the specimen, the pressurized water supply and the water collection system were connected to the triaxial cell. During the internal erosion testing, the applied hydraulic gradient, vertical displacement, lateral deformation, outflow rate, and eroded soil mass were measured continuously. Conventional drained triaxial compression or extension test was then carried

out to study the stress-strain behavior of the soil subjected to internal erosion.

A series of drained triaxial tests on samples before and after erosion on a gap-graded HK-CDG mixtures (see Chang 2012) were selected for simulations: BM-C-50, BM-C-100, and BM-C-200 are conventional drained triaxial tests under confining stresses of 50, 100, and 200 kPa, respectively on identical samples before erosion test; GS-C-50, GS-C-100, and GS-C-200 are conventional drained triaxial tests under confining stresses of 50, 100 and 200 kPa, respectively, on isotropically eroded samples; GS-C-4 and GS-C-5 are drained triaxial tests under a confining stress of 50 kPa on anisotropically eroded samples (under $q = 50$ and 100 kPa, respectively).

The model parameters were calibrated as follows (see Fig. 10): (a) the elastic parameters ($k_{n0}$ and $n$) were determined from an isotropic compression test with an assumed value of $k_{fR}$ = 0.5; (b) the friction angle $\phi_{\mu}$ was determined from stress paths of triaxial tests; (c) the CSL related parameters ($\zeta$, $\lambda$) were measured from drained triaxial tests. It was difficult to determine these two parameters because strain localizations appeared systematically at the end of the tests (see Chang 2012); due to a lack of experimental data the assumption of parallel critical state lines has been made, in
<figure><img src="image_14.png" /><figcaption>Fig. 10. Comparison between experimental results and simulations for HK decomposed granites before and after internal erosion: (a, c) deviatoric stress versus axial strain, and (b, d) void ratio versus axial strain.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 11. Comparison between experimental results and simulations for HK decomposed granites before and after internal erosion: (a, c) deviatoric stress versus axial strain, and (b, d) void ratio versus axial strain.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 12. Local stress–strain behavior of silty sand with an inter-granular void ratio of 0.8.</figcaption></figure>

agreement with the results presented in Fig. 4 for various sand-silt mixtures. (d) due to the lack of data concerning Chang's tests, the values of the other parameters ($e_{hc,cro}$, $e_{hf,cro}$, $a$, $m$ and $f_{th}$) had to be assumed based on published data of other materials in Figs. 3 and 4, only $e_{cro}$ was measured (see Fig. 10b); (e) finally, the plastic parameters ($G_p$, $D$) were determined by curve fitting on the triaxial test of HK-CDG mixtures before erosion. All determined parameters summarized in Table 2 were used to predict the other tests.

Presented in Fig. 11 is a comparison between the experimental and the numerical results of drained triaxial tests on samples before and after erosion. Both experimental and numerical results show that the loss of fine particles during internal erosion affects significantly the mechanical behavior of the soil. When the soil loses a significant amount of fine particles, its stress–strain behavior changes from a dilutive to a contractive behavior. When a part of the solid fraction is removed from the material, the void ratio increases and, as a consequence, the sliding resistance decreases, which can create additional inter-particle displacements, and a greater deformation of the granular assembly. This has been observed experimentally during the erosion of the specimen: at constant stress, the specimen is undergoing straining. Under an isotropic stress state, the specimen volume decreases during erosion, whereas under an anisotropic stress state, the axial strain is also increasing during erosion and this increase is more pronounced for higher deviatoric stresses.

The proposed model was able to reproduce the general trend of the soil behavior until the peak strength which corresponds in the experiments to the development of strain localization within the specimen. In some cases, discrepancies between experimental and numerical results could be found, possibly due to the difficulty of measuring accurately certain material properties, such as the

elastic stiffness or the friction angle, as a function of the fines content. Another aspect which is not taken into account in our simulations is that the state of the sample after the erosion process is not known. Some heterogeneity within the specimen could be produced by the combined effects of transport and filtration of the fines particles. No detail of this particular aspect was given by Chang. A future paper will be devoted to a detailed application of our model to the problem of internal erosion.

### 3.4. Micromechanical analysis for undrained behavior of sand-silt mixtures

In this section, we investigate the local stress–strain behavior for selected contact planes with different orientations $\theta = 0°, 18°, 28°, 45°, 55°, 72°$ and $90°$ (see Fig. 2, $x$ representing vertical direction). The undrained triaxial tests on silty sand samples at $e_g = 0.8$ (Fig. 7) and on sandy silt samples at $e_f = 0.9$ (Fig. 7) of Foundry sand-silt mixtures were selected as examples for illustrating the local behavior in contact planes (see Figs. 12 and 13), wherein the stress–strain behavior at the local scale is presented to examine the impact of fines to silty sand at $e_g = 0.8$ and the impact of coarse grains to sandy silt at $e_f = 0.9$ on the stress–dilatancy and undrained strength. General observations can be made from these Figs.:

* Figs. 12(a) and 13(a) present the local stress paths for the selected contact orientations (from a-1 to a-5, see also Fig. 2) for all cases ($f_c = 0%, 7\%$ and $15\%$ for silty sand, and $f_c = 40%, 60\%$ and $100\%$ for sandy silt). For the $0°$ and $90°$ contact planes, the shear stress is null under triaxial condition. For the others, the stress state with the maximum stress ratio is located in the $55°$ contact plane. Generally
<figure><img src="image_17.png" /><figcaption>Fig. 13. Local stress-strain behavior of sandy silt with an inter-fine void ratio of 0.9.</figcaption></figure>

for all contact planes, both the fines for silty sand and the coarse grains for sandy silt contribute reinforcement to the shear stress, which agrees well with the global behavior observed.

◇ The local shear stress-strain curves shown in Figs. 12(b) and 13(b) clearly indicate that all contact planes are mobilized with varied degrees. The planes with the largest movements are located around the orientation of 55° (close to $\pi/4 + \phi_{\mu}2 = 60°$). These active contact planes contribute largely to the overall deformation of the specimen. Among all the selected contact planes, only the 55° contact plane has a local shear stress-strain curve similar to the global one. For all contact planes, adding fines in silty sand or adding coarse grains in sandy silt generally results in a higher level of shear stress, which is consistent with the global behavior.

◇ Figs. 12(c) and 13(c) present the local normal stress-strain curves, in which normal strains are mostly composed by elastic strains due to the change of normal stresses and plastic strains induced by the dilatancy law during shearing. The results show that, although the global volumetric strain is constrained to be zero during undrained shearing, the local normal strains develop in every contact plane to different degrees. In general, the higher fines content in silty sand or more coarse grains in sandy silt leads to a higher normal stress level which agrees with the global behavior (see Fig. 7).

Evidently, for samples with a different fines content, the local contact planes demonstrate quite different responses, and the local

stress-strain behavior is significantly affected by the fines content as well. These observations are in good agreement with the global behavior. The local stresses and strains at the contact planes are highly non-uniform and evolve constantly throughout the triaxial test, which may reasonably account for the development of anisotropy induced by external loads since some properties are dependent of local stresses and strains. We could not find in the literature experimental analyses at the microscale which could support this analysis. In this context, DEM analyses could be useful, but up to now the mixture of fine and coarse grains cannot be simulated due to high computational cost. This could be the topic of future researches. What we wanted to illustrate with these simulations was how the model works at the microscale and that the behavior at grain scale could be significantly affected by the fines content.

## 4. Conclusions

The micromechanical stress-strain model of Chang and Hicher (2005) has been adapted based on the approach proposed by Chang and Yin (2011) in order to reproduce the behavior of silt-sand mixtures. The model takes into account the influence of fines in silty sand and the influence of coarse grains in sandy silt on soil dilatancy and strength by formulating the evolution of the critical state line of sand-silt mixtures from sand to silt based on experimental results. All model parameters can be determined from conventional triaxial tests and index void ratio tests.

The proposed model has been employed to simulate the behaviors of Foundry sand-silt mixtures from sand ($f_c = 0\%$) to silt
(f_c = 100%) and Ottawa sand-silt mixtures with different fines contents (f_c = 0-50%) in undrained triaxial tests, as well as HK-CDG mixtures before and after erosion in drained triaxial tests. A good comparison has been observed between the model predictions with experimental results in all three cases. The model captures well the influence of the fines content from sand to silt on the overall behavior of soil mixtures. Thus, it may prove to be useful for simulating internal erosion in soil.

The predicted local behavior in the contact planes has also been examined in the case of undrained triaxial tests with different fines contents to clarify the fines effect in silty sand and the coarse grains effect in sandy silt. The local stress-strain responses in the contact planes agree well with the global behavior of soil mixtures from sand to silt. It is shown that all local contact planes are mobilized to different degrees in terms of local stress and strain. A few active contact planes contribute dominantly to the deformation of the assembly, while most contact planes only play a marginal role in the overall specimen strain. While the local stresses and strains are evidently highly non-uniform from the model simulation, an anisotropic global behavior may be reasonably expected when a soil is subjected to external shear. Hence, the micromechanical model provides an effective means to probe the cross-scale behaviors in soils.

## Acknowledgements

This study was financially supported by the National Natural Science Foundation of China (Grant No. 41372285, 511611 30523), the Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20110073120012), the Research Grants Council of Hong Kong (RGC/GRF 623609), the ANR-NSFC grant for the project RISMOGEO and the project C2D2 LEVEES financed by the French Ministry of Ecology. Part of the research was conducted during a stay of the first author at Hong Kong University of Science and Technology as visiting professor.

## Appendix A1. Elasto-plastic relationship

The incremental expression of the local behavior at inter-particle contact is given by:

$$ \dot{f}_{i}^{\alpha} = k_{ij}^{\alpha p} \dot{\delta}_{j}^{\alpha} \quad (A1) $$

## Appendix A2. Macro micro relationship

The stress-strain relationship for an assembly can be determined from integrating the behavior of inter-particle contacts in all orientations. During the integration process, a relationship is required to link the macro and micro variables. Using the static hypotheses, we obtain the relation between the macro strain and inter-particle displacement (here, we do not consider the finite strain condition)

$$ \dot{u}_{j,i} = A_{ik}^{-1} \sum_{\alpha=1}^{N} \dot{\delta}_{j}^{\alpha} l_{k}^{\alpha} \quad (A2) $$

where $\dot{\delta}_j$ is the relative displacement between two contact particles and the branch vector $l_k$ is the vector joining the centers of two contact particles. It is noted that contact particles include both direct contact and indirect contact of neighboring particles associate with a Voronoi polyhedron. For convenience, we let N be the total number of contact orientations. The variables $\dot{\delta}_j^\alpha$ and $l_k^\alpha$ are defined, respectively as the averaged values of $\dot{\delta}_j$ and $l_k$ for all contacts belonging to the $\alpha^{th}$ orientation. The fabric tensor in Eq. (A2) is defined as:

$$ A_{ik} = \sum_{\alpha=1}^{N} l_{i}^{\alpha} l_{k}^{\alpha} \quad (A3) $$

Using the principle of energy balance, the mean force on the contact plane of each orientation are:

$$ \dot{f}_{j}^{\alpha} = \dot{\sigma}_{ij} A_{ik}^{-1} l_{k}^{\alpha} V \quad (A4) $$

The stress increment $\dot{\sigma}_{ij}$ can be obtained by the contact forces and branch vectors for contacts in all orientations:

$$ \dot{\sigma}_{ij} = \frac{1}{V} \sum_{\alpha=1}^{N} \dot{f}_{j}^{\alpha} l_{i}^{\alpha} \quad (A5) $$

## Appendix A3. Computation scheme

The problem is defined by the following:

Initially, we know the global variables ($\sigma_{ij}$ and $\epsilon_{ij}$) for the assembly and the local variables ($f_j^\alpha$ and $\delta_j^\alpha$) for each contact orientation. For a given loading increment, which can be stress control, strain control or mixed mode, out of the 12 variables ($\Delta\sigma_{ij}$ and $\Delta\epsilon_{ij}$), 6 of them are unknown. The objective is to determine all global variables ($\sigma_{ij}$ and $\epsilon_{ij}$) and local variables ($f_j^\alpha$ and $\delta_j^\alpha$) at the end of load increment. For a system with N inter-particle orientations, the number of unknown is 3 N for $f_j^\alpha$ and 3 N for $\delta_j^\alpha$. The total number of unknown is 3 N + 3 N + 6.

The following constraints must be satisfied:

(1) The local constitutive equation. Since there are three equations for each contact plane orientation, the total number of equations is 3N, N being the total number of inter-particle orientations.

(2) The static hypothesis between global stress and local forces: the number of equations is 3N.

(3) The strain definition between global strain and local displacement: the number of equations is 6.

The total number of unknowns is the same as the total number of equations. Therefore, the solution can be determined.

Using the above equations, the following relationship between stress and strain can be obtained:

$$ \dot{u}_{i,j} = C_{ijmp} \dot{\sigma}_{mp}; \quad \text{where} \quad C_{ijmp} = A_{ik}^{-1} A_{mn}^{-1} V \sum_{\alpha=1}^{N} (k_{jp}^{\alpha p})^{-1} l_{k}^{\alpha} l_{n}^{\alpha} \quad (A6) $$

The summation can be replaced by an integral over orientations. The integral can lead to a closed-form solution for the elastic modulus of randomly packed equal-size particles. However, in the elastic plastic behavior, due to the non-linearity nature of the local constitutive equation, a numerical calculation with iterative process is necessary to carry out the summation. In order to facilitate the numerical calculation, the orientations are selected to coincide with the locations of Gauss integration points in a spherical coordinate. Summation over these orientations with the Gauss weighting factor for each orientation is equivalent to determining the integral over orientations. We found that the results were more accurate by using a set of fully symmetric integration points. From a study of the performance of using different numbers of orientations, we found N ≥ 74 to be adequate.

For a strain-controlled test, the previous equation is not useful especially at the after-peak range of strain-softening. In this case, a method of "elastic predictor-plastic corrector" was adopted to obtain the solution. For a mix-mode loading condition, additional process of distributing the unbalance stresses is needed.
## References

Belkhatir, M., Arab, A., Schanz, T., Missoum, H., Della, N., 2011. Laboratory study on the liquefaction resistance of sand-silt mixtures: effect of grading characteristics. *Granular Matter* 13 (5), 599–609.

Biarez, J., Hicher, P.Y., 1994. *Elementary Mechanics of Soil Behaviour*. Balkema, Rotterdam, The Netherlands.

Chang, D.S., Zhang, L.M., 2011. A stress-controlled erosion apparatus for studying internal erosion in soils. *Geotech. Test*. J. 34 (6), 1945–7545.

**Chang, D. S., 2012.** Internal erosion and overtopping erosion of earth dams and landslide dams. PhD thesis of The Hong Kong University of Science and Technology.

Chang, C.S., Hicher, P.Y., 2005. An elastic-plastic model for granular materials with microstructural consideration. *Int. J. Solids Struct*. 42 (14), 4258–4277.

Chang, C.S., Misra, A., 1990. Application of uniform strain theory to heterogeneous granular solids. *J. Eng. Mech. ASCE* 116 (10), 2310–2328.

Chang, C.S., Yin, Z.-Y., 2010a. Modeling stress-dilatancy for sand under compression and extension loading conditions. *J. Eng. Mech. ASCE* 136 (6), 777–786.

Chang, C.S., Yin, Z.-Y., 2010b. Micromechanical modelling for inherent anisotropy in granular materials. *J. Eng. Mech. ASCE* 136 (7), 830–839.

Chang, C.S., Yin, Z.-Y., 2011. Micromechanical modeling for behavior of silty sand with influence of fine content. *Int. J. Solids Struct*. 48 (19), 2655–2667.

Chang, C.S., Yin, Z.-Y., Hicher, P.Y., 2011. Micromechanical analysis for inter-particle and assembly instability of sand. *J. Eng. Mech. ASCE* 137 (3), 155–168.

Guo, N., Zhao, J.D., 2013. The Signature of shear-induced anisotropy in granular media. *Comput. Geotech*. 47, 1–15.

Hicher, P.Y., Chang, C.S., 2006. Anisotropic nonlinear elastic model for particulate materials. *J. Geotech. Geoenviron. Eng. ASCE* 132 (8), 1052–1061.

Lade, P.V., Yamamuro, J.A., 1997. Effects of nonplastic fines on static liquefaction of sands. *Can. Geotech*. J. 34 (6), 918–928.

Li, X.S., Dafalias, Y.F., 2000. Dilatancy for Cohesionless Soils. *Geotechnique* 50 (4), 449–460.

Li, X.S., Wang, Y., 1998. Linear representation of steady-state line for sand. *J. Geotech. Geoenviron. Eng. ASCE* 124 (12), 1215–1217.

Manzari, M.T., Dafalias, Y.F., 1997. A critical state two-surface plasticity model for sands. *Geotechnique* 47 (2), 255–272.

Murthy, T.G., Loukidis, D., Carraro, J.A.H., Prezzi, M., Salgado, R., 2007. Undrained monotonic response of clean and silty sands. *Géotechnique* 57 (3), 273–288.

Naeini, S.A., Baziar, M.H., 2004. Effect of fines content on steady-state strength of mixed and layered samples of sand. *Soil Dyn. Earthquake Eng*. 24 (3), 181–187.

Papadopoulou, A., Tika, T., 2008. The effect of fines on critical state and liquefaction resistance characteristics of non-plastic silty sands. *Soils Found*. 48 (5), 713–725.

Pitman, T.D., Robertson, P.K., Sego, D.C., 1994. Influence of fines on the collapse of loose sands. *Can. Geotech*. J. 31 (5), 728–739.

Polito, C.P., Martin, J.R., 2001. Effects of nonplastic fines on the liquefaction resistance of solids. *J. Geotech. Geoenviron. Eng. ASCE* 127 (7), 408–415.

Salgado, R., Bandini, P., Karim, A., 2000. Shear strength and stiffness of silty sand. *J. Geotech. Geoenviron. Eng. ASCE* 126 (5), 451–462.

Thevanayagam, S., Mohan, S., 2000. Intergranular state variables and stress-strain behaviour of silty sands. *Géotechnique* 50 (1), 1–23.

Thevanayagam, S., Shenthan, T., Mohan, S., Liang, J., 2002. Undrained fragility of clean sands, silty sands and sandy silts. *J. Geotech. Geoenviron. Eng. ASCE* 128 (10), 849–859.

**Vaid, Y.P., 1994.** Liquefaction of silty soils. Proc., Ground failures under seismic conditions, *Geotech. Spec. Publ.* 44, ASCE, New York, 1–16.

Xenaki, V.C., Athanasopoulos, G.A., 2003. Liquefaction resistance of sand-silt mixtures: an experimental investigation of the effect of fines. *Soil Dyn. Earthquake Eng*. 23 (3), 183–194.

Yang, S.L., Sandven, R., Grande, L., 2006. Instability of sand-silt mixtures. *Soil Dyn. Earthquake Eng*. 26 (2–4), 183–190.

Yin, Z.-Y., Chang, C.S., Hicher, P.Y., 2010. Micromechanical modelling for effect of inherent anisotropy on cyclic behaviour of sand. *Int. J. Solids Struct*. 47 (14–15), 1933–1951.

Yin, Z.-Y., Chang, C.S., 2012. Stress-dilatancy for sand under loading and unloading conditions. *Int. J. Numer. Anal. Meth. Geomech*. 37 (8), 855–870.

**Zlatovic, S., Ishihara, K., 1995.** On the influence of nonplastic fines on residual strength. In *Proceedings of IS-TOKY095, The First International Conference on Earthquake Geotechnical Engineering, Tokyo, 14–16 November 1995*. Edited by K. Ishihara. A.A. Balkema, Rotterdam, the Netherlands. vol. 1, pp. 239–244.