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

HAL
open science

# Constitutive model for soil-rock mixtures in the light of an updated skeleton void ratio concept

Tao Wang, Sihong Liu, Antoine Wautier, François Nicot

► To cite this version:

Tao Wang, Sihong Liu, Antoine Wautier, François Nicot. Constitutive model for soil-rock mixtures in the light of an updated skeleton void ratio concept. Acta Geotechnica, 2023, 18, pp.2991-3003. 10.1007/s11440-022-01756-6. hal-03852498

HAL Id: hal-03852498

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

Submitted on 15 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Constitutive model for soil-rock mixtures in the light of an

updated skeleton void ratio concept

Tao WANGª, Sihong LIUª,*, Antoine WAUTIERᵇ,ᶜ, François NICOTᵇ,ᵈ

ªHohai University, No.1 Xikang Road, Nanjing 210098, China.

ᵇUniversité Grenoble Alpes, INRAE, UR ETNA, 2 rue de la Papeterie-BP 76, F-38402

St-Martin-d'Hères, France.

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

Aix-en-Provence, Cedex 5, France.

ᵈUniversité Savoie Mont Blanc, EDYTEM (USMB/CNRS), Le Bourget du Lac, France

10

11 Email address:

wt@hhu.edu.cn (Tao WANG)

sihongliu@hhu.edu.cn (Sihong LIU, Corresponding author)

antoine.wautier@inrae.fr (Antoine WAUTIER)

francois.nicot@univ-smb.fr (François NICOT)

16

17

18

19

20
**Abstract:** As a type of special geological body, soil-rock mixtures (SRMs) are widely found in nature and used in civil engineering. Many structures, such as rockfill dams, highways and tunnels have used SRMs as building materials. Proper modelling of SRMs is of great importance to capture the complex behavior of this heterogeneous material. In this manuscript, a simple constitutive model incorporating the skeleton void ratio concept is proposed for SRMs with varying soil contents (*sc*). A prominent feature of the model is a unified description of the behavior of SRMs with varying *sc* such that only model parameters of pure rock and of pure soil are required. After calibration, the model shows a good capacity to predict the stress-strain response of SRMs under a wide range of *sc*, void ratios, and confining pressures. In particular, it captures well the non-associated behavior of rock-dominated SRMs with different *sc*. Furthermore, the *sc*-value is shown to modify the plastic flow direction of the material without influencing its yield surface.

**Keywords:** Soil-rock mixture; Skeleton void ratio; Soil content; Constitutive modelling; Critical state line; Flow rule direction; Non-associated behavior

# 1. Introduction

Soil-rock mixtures (SRMs) are heterogeneous materials composed of high-strength rocks, fine-grained soils and pores [37, 38, 42]. SRMs are widely encountered in geotechnical engineering, such as natural slopes [6], waste rocks and tailings from mining [8], clay-aggregate mixtures in rockfill dams [12] and tunnels [14]. The soil content (*sc*) is one of the most important factors governing the mechanical behavior of SRMs [3, 26, 29, 33, 36]. In recent years, many field tests, laboratory tests and numerical simulations have revealed that *sc* greatly affects the shear
strength [4, 16, 23, 37], failure modes [5, 15, 17], stress-dilatancy [4, 33], and critical state
parameters of SRMs [25, 33]. Experimental and numerical results both show that at a low *sc*, the
mechanical behavior of SRMs is primarily governed by inter-granular friction between rock grains.
While at a high *sc* beyond a threshold value (*sc*)<sub>th</sub>, the mechanical behavior of SRM is primarily
governed by friction characteristics of soil grains.

Although some pioneering works has been done on the constitutive modelling of SRMs [7,
17, 26], properly modeling the behavior of this kind of complex heterogeneous materials remains
an open challenge [35, 39, 40]. Predicting the mechanical properties SRMs for varying sc remains
a widely open issue.

1) In most constitutive models, SRMs with different *sc* are indeed treated as different
materials with their own *sc*-specific parameters [21, 24, 41]. For example, the shape of the critical
state line (CSL) of a SRM and its location depend on *sc*. Therefore, new laboratory tests are
required each time *sc* is updated in engineering projects, which results in a waste of time and
money.

2) Some models tried to build empirical equations to link model parameters with sc [23, 35].
For instance, some empirical equations have been established to fit CSLs of SRMs with sc based
on experimental data. This method usually introduces new parameters into the constitutive model,
i.e., parameters in empirical equations. But, these empirical equations lack of physical meaning
and their application to SRMs with different lithology, grain shape and gradation of rock and soils
remains questionable.

3) Although many numerical studies have analyzed the microstructure of SRMs, few
attempts have been made to incorporate micro- or meso-mechanisms into constitutive models. For
example, Wang et al (2021) [29] found that for rock-dominated SRMs, *sc* does not affect the normal direction of yield surface but it changes the flow rule direction. These findings could be reflected in constitutive models.

In soil constitutive models, the global void ratio *e* has been chosen as one of the main state variable. However, it was found to be an imperfect index to characterize the mechanical behavior of mixed soils like SRMs [29, 33]. This is because such a global index is not able to account for the non-active participation of the soil grains in the force transmission structure within a SRM. Alternatively, skeleton void ratio turns out to be a more appropriate index to reflect the density of SRMs. Skeleton void ratio corresponds to the void ratio of grains constituting the stress-bearing skeleton. In recent studies, this index has shown a strong potential to give unified descriptions of the behavior of SRMs with varying *sc* [29, 35].

The objective of this manuscript is to propose a simple method to predict the stress-strain response of SRMs with varying *sc* incorporating the updated skeleton void ratio concept. It is organized as follows. First, the updated skeleton void ratio for SRMs proposed by Wang et al. (2022) [29] is briefly reviewed. The advantage of using this skeleton void ratio index to characterize critical state lines of SRMs with varying *sc* is then shown. Next, a constitutive model incorporating the skeleton void ratio concept is proposed and the model is validated against experimental results. Eventually, effect of *sc* on the main properties of SRMs is studied and the model capabilities are discussed.

## 2. Review of an updated skeleton void ratio for SRMs

Wang et al. (2022) [29] proposed an updated skeleton void ratio index for SRMs. The main advantage of this updated skeleton void ratio is that it can consider the effect of gradations of both
soils and rocks. The skeleton void ratio proposed by Thevanayagam (2007) [27] is a special case of the updated skeleton void ratio in which mutual interaction between rock and soil grains during packing is neglected.

According to Wang et al. (2022) [29], the threshold soil content ($sc$)$_{th}$ that separates the rock-dominated structure and soil-dominated structure is:

$$ (sc)_{th} = \frac{e_r - b}{1 + e_r + e_s - a - b} \quad (1) $$

where $e_s$ and $e_r$ are the minimum void ratios of pure soil grains and pure rock grains, respectively. $a$ and $b$ are gradation-related parameters:

$$ a = (1 + e_s)\exp\left(\frac{-R_d^{0.5}}{C_{ur}C_{us}}\right) \quad (2) $$

$$ b = e_r\exp\left(\frac{-R_d}{C_{ur}C_{us}^{0.5}}\right) \quad (3) $$

where $C_{ur} = \frac{(D_{60})_{rock}}{(D_{10})_{rock}}$ is the coefficient of non-uniformity for the rock fraction, $C_{us} = \frac{(D_{60})_{soil}}{(D_{10})_{soil}}$ is the coefficient of non-uniformity for the soil fraction and $R_d = \frac{(D_{50})_{rock}}{(D_{50})_{soil}}$ is grain size disparity ratio.

For a rock-dominated structure, i.e., $sc < (sc)_{th}$, its skeleton void ratio $e_{sk}$ is given by:

$$ e_{sk} = \frac{e + sc}{1 - sc} - \frac{sc(1 + e_s)}{1 - sc} \exp\left(\frac{-R_d^{0.5}}{C_{ur}C_{us}}\right) \quad (4) $$

For a soil-dominated structure, i.e., $sc > (sc)_{th}$, its skeleton void ratio $e_{sk}$ expresses as:

$$ e_{sk} = \frac{e}{sc} - \frac{e_r(1 - sc)}{sc} \exp\left(\frac{-R_d}{C_{ur}C_{us}^{0.5}}\right) \quad (5) $$

All the parameters introduced in Equations (1)-(5) can be obtained from simple laboratory sieving and compaction tests.

Wang et al. (2022) [29] found that a rock-dominated (or soil-dominated) SRM has similar stress-strain responses with a pure rock specimen (or a pure soil specimen) if global void ratio of
the pure rock specimen (or pure soil specimen) equals to the skeleton void ratio of the SRM. This important finding is essential to unify descriptions of SRMs with varying soil contents.

## 3. Characterization of critical state lines of SRMs with updated skeleton void ratio index

The critical state is defined as the state at which the soil continues to deform at constant shear stress and constant volume. It has increasingly been used as a fundamental concept to characterize the strength and deformation properties of soils [2, 20, 22]. Li et al. (1998) [10] found that the critical state lines (CSLs) for cohesionless soils are straight lines in the e-(p/p_a)⁵ plane:

$$e_{cr} = e_{\Gamma} - \lambda(p/p_{a})^{\xi} \qquad (6)$$

Where ecr is critical state void ratio, eΓ is the theoretical critical void ratio at the atmospheric pressure, pa is the atmospheric pressure, λ is the magnitude of the slope, and ξ is the pressure exponent (with a typical value around 0.7).

Experimental data obtained from conventional drained triaxial tests on SRMs with different soil contents (sc=0, 0.1, 0.3, 0.5, 0.7 and 1) under different confining pressures (σ₃=150 kPa, 300 kPa and 600 kPa) conducted by Wang et al. (2022) [29] are adopted here to characterize the CSLs of SRMs with varying sc (Cur = 2.36, Cus = 7.01 and Rd=5.46). Detailed test procedures can be found in the quoted reference [29].

Figure 1 gives the CSLs of SRMs with sc ranging from 0 to 1 in e-p plane. It can be found that SRMs with different sc have different CSLs. The CSL move downward from the pure rock specimen until the soil content reaches (sc)th (around 0.54), and then the CSL would move upward to the position of pure soil specimen. According to Wang et al. (2022), all the SRMs were prepared with the same global void ratio. For a rock-dominated specimen, the skeleton void ratio
increases with sc and a more contractant behavior is observed. Therefore, the increase in sc widens the gap between current void ratio and critical void ratio and leads to the downward shift of the CSL. On the other hand, for a soil-dominated specimen, the increase in sc results in a larger skeleton void ratio and a more dilatant behavior, which consequently narrows the gap. Therefore, CSL will move upwards. Critical state line parameters e<sub>r</sub> and λ are shown in Figure 2. e<sub>r</sub> and λ are found to be functions of sc, which indicates that SRMs with different sc should be treated as different materials. Consequently, each SRM has its own CSL in e-p plane.

<figure><img src="image_3.png" /><figcaption>Fig. 1 Critical state lines of SRMs with different sc drawn in e-(p/p<sub>a</sub>)<sup>ξ</sup> plane</figcaption></figure>

<figure><img src="image_4.png" /></figure>
(a) $e_r$

(b) $\lambda$

<caption>Fig. 2 Evolution of critical state parameters with sc: (a) $e_r$; (b) $\lambda$</caption>

Critical state data of SRMs with different sc are plotted in $e_{sk}-p$ plane in Figure 3. Pure rock and pure soil materials have different CSLs because they have different grain shapes and gradations. Therefore, for a clearer view, we plot separately the critical state data of rock-dominated structure (Figure 3a) and soil-dominated structure (Figure 3b). It can be seen in Figure 3a, that the CSL of pure rock grains (sc=0) is shared with SRMs of different sc (sc=0, 0.1, 0.3 and $0.5<(sc)_{th} = 0.54$) for rock-dominated SRMs. Similarly, as shown in Figure 3b, for the soil-dominated specimens (sc=0.7, $1>(sc)_{th} = 0.54$), the CSL of the pure soil (sc=1) is shared with SRM with sc=0.7. Consequently, in the $e_{sk}-p$ plane, the critical state of all RSM materials is described either by the CSL of the rock or the SCL of the soil depending whether sc is below or above (sc)_{th}. This important finding is essential for establishing the constitutive model, as developed in the next section.

<figure><img src="image_5.png" /><figcaption>(a) Rock-dominated SRMs</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>(b) Soil-dominated SRMs
Fig. 1 Critical state lines of SRMs with different sc drawn in $e_{sk} - (p/p_a)^\xi$ plane</figcaption></figure>

## 4. A simple constitutive model for SRM using the skeleton void ratio index

The updated skeleton void ratio concept is incorporated into the state-dependent bounding surface model proposed by Li and Dafalias (2000) [9] to capture the stress-strain behaviors of SRMs with different sc. The fundamental of this model is that the SRMs and its host rock (or soil) with the same skeleton void ratio should exhibit the same stress-strain behavior, which has been reported by many researchers [28, 29, 35]. The fact that rock-dominated (or soil-dominated) SRMs share the same CSL with pure rock (or pure soil) also suggests to build a unified model for SRMs with varying sc from only model parameters of pure rock and pure soil. The constitutive framework is detailed in the following sections.

### 4.1 Elastic behavior

According to standard elasto-plasticity, a total strain increment $d\boldsymbol{\epsilon}$ is additively split into an elastic strain increment $d\boldsymbol{\epsilon}^e$ and a plastic strain increment $d\boldsymbol{\epsilon}^p$:

$$d\boldsymbol{\epsilon} = d\boldsymbol{\epsilon}^e + d\boldsymbol{\epsilon}^p \qquad (7)$$
The volumetric and deviatoric elastic strain increments are given respectively by:

$$d\epsilon_v^e = \frac{dp'}{K} \qquad (8)$$

$$d\epsilon_d^e = \frac{dq}{3G} \qquad (9)$$

where K is the elastic bulk modulus, G is the elastic shear modulus, $p' = \frac{\sigma_1+2\sigma_3}{3}$ is the mean effective stress and $q = \sigma_1 - \sigma_3$ is the deviatoric stress ($\sigma_1$, $\sigma_2$, and $\sigma_3$ are the principal stress values).

G can be estimated following [19]:

$$G = G_0 \frac{(2.97 - e)^2}{(1 + e)} \sqrt{p'p_a} \qquad (10)$$

where $G_0$ is an elastic material constant, $p_a$ is the atmospheric pressure, e is the void of the considered material (soil or rock).

The elastic bulk moduli K is given by the following relation:

$$K = G \frac{2(1+\nu)}{3(1-2\nu)} \qquad (11)$$

where $\nu$ is the Poisson's ratio.

The estimation of $G_0$ and $\nu$ can be found in Appendix.

## 4.2 Plastic behavior

For the sake of simplicity, the elasto-plastic model proposed by Li and Dafalias (2000) [9] in the triaxial compression stress space is adopted. In this model, the yield surface is given by:

$$f(p', q, \eta) = q - p'\eta_c = 0 \qquad (12)$$

where $\eta_c$ is the stress ratio when plasticity activates. In the model of Li and Dafalias (2000) [9] the yield surface is assumed to follow the stress state so that $\eta_c = \eta$. Note that this hypothesis holds as long as plasticity is activated. In case elastic unloading is considered, the model requires some additional mechanism such as a back-stress [13], or a memory of the reversal stress ratio
175 point [31] which are not considered in the present study.

176 According to the theory of plasticity, a loading index $dL$ can be defined as:

$$dL = \frac{1}{K_p} \left( \frac{\partial f}{\partial p'} dp' + \frac{\partial f}{\partial q} dq \right) = \frac{dq - \eta_c dp'}{K_p} = \frac{p'd\eta}{K_p} \quad \text{since } \eta_c = \eta \quad (13)$$

177 where $K_p$ is the plastic hardening modulus and is expressed by Li and Dafalias (2000) [9] as:

$$K_p = hG \left( \frac{M}{\eta} - \exp(n\psi) \right) = \frac{hG\exp(n\psi)}{\eta} (M\exp(-n\psi) - \eta) \quad (14)$$

178 In the above equation, $n$ is a positive model parameter, $h = h_1 - h_2 e_0$ with $h_1$ and $h_2$ are 179 model parameters and $e_0$ the initial void ratio, $M$ is the stress ratio at a critical state and 180 $\psi = e - e_{cr}$ is the state parameter defined by Been and Jefferies (1985) [1] with $e$ the current 181 void ratio and $e_{cr}$ the critical state void ratio for the current $p'$. When using skeleton void ratio 182 instead of the void ratio, a skeleton state parameter $\psi_{sk}$ is introduced as $\psi_{sk} = e_{sk} - (e_{cr})_{sk}$, as 183 seen in Figure 4. Note that when softening occurs ($K_p < 0$) the plastic multiplier $dL$ is positive 184 (i.e. plasticity is active) when the stress ratio decreases.

<figure><img src="image_7.png" /><figcaption>185 Fig. 4 Definition of skeleton state parameter $\psi_{sk}$ 186</figcaption></figure>

187 Since soil-rock mixtures are typical non-associated materials, the non-associated flow rule is 188 adopted to define the plastic strain increments, as follows:

$$d\varepsilon_d^p = dL = \frac{p'd\eta}{K_p} \quad (15)$$

$$d\varepsilon_v^p = DdL = \frac{Dp'd\eta}{K_p} \quad (16)$$
where $d\epsilon_d^p$ and $d\epsilon_v^p$ are the plastic deviatoric strain increment and the plastic volumetric strain increment, respectively. The ratio $D = \frac{d\epsilon_v^p}{|d\epsilon_d^p|}$ evaluates the amplitude of the dilatancy and can be chosen as follows:

$$D = \frac{d_0}{M} [M \exp(m\psi) - \eta] \quad (17)$$

where $d_0$ and $m$ are model parameters. For associated flow rule, $D = \frac{d\epsilon_v^p}{d\epsilon_d^p} = \frac{\frac{\partial f}{\partial p}}{\frac{\partial f}{\partial q}} = -\eta_c$.

Therefore, for $dL>0$, the following expressions can be derived for both incremental deviatoric and volumetric strains:

$$d\epsilon_d = d\epsilon_d^e + d\epsilon_d^p = \frac{dq}{3G} + \frac{P'd\eta}{K_p} = \left(\frac{1}{3G} + \frac{1}{K_p}\right) dq - \frac{\eta}{K_p} dp' \quad (18)$$

$$d\epsilon_v = d\epsilon_v^e + d\epsilon_v^p = \frac{dp'}{K} + D d\epsilon_q^p = \frac{D}{K_p} dq + \left(\frac{1}{K} - \frac{D\eta}{K_p}\right) dp' \quad (19)$$

Equations (18) and (19) set up the relationship between stain and stress increments.

Eventually, an elastoplastic constitutive relation is derived in a matrix form:

$$\left\{ \begin{matrix} dq \\ dp' \end{matrix} \right\} = \left[ \begin{matrix} 3G & 0 \\ 0 & K \end{matrix} \right] - \frac{h(dL)}{K_p + 3G - K\eta D} \left( \begin{matrix} 9G^2 & -3KG\eta \\ 3KGD & -K^2\eta D \end{matrix} \right) \left\{ \begin{matrix} d\epsilon_q \\ d\epsilon_v \end{matrix} \right\} \quad (20)$$

where $h(dL)$ is a Heaviside function with $h(L)=1$ for $dL>0$ and $h(dL)=0$, otherwise (plastic stain exists only when $dL$ is positive). Note that equation (20) is given in the plane $(dp', dq)$ and the dual plane $(d\epsilon_v, d\epsilon_q)$. It can be generalized to any kind of incremental stress and strain tensors by assuming coaxiality between incremental stress and strain.

## 4.3 Model parameters

To summarize, the proposed model includes eleven model parameters, all of them being calibrated from drained triaxial tests under different confining pressures: (1) Two elastic parameters, i.e. $G_0$ and $\nu$; (2) Four critical state parameters, i.e. $M$, $e_\Gamma$, $\lambda$ and $\xi$; (3) Two dilatancy parameters, i.e. $d_0$, $m$; (4) Three hardening parameters, i.e. $h_1$, $h_2$ and $n$.

It should be noted that no new parameters have been introduced into this approach. Skeleton
void ratio $e_{sk}$ of SRMs with different $sc$ are adopted to replace global void ratio $e$ in the above-mentioned equations. Accordingly, skeleton state parameter $\psi_{sk}$ is used in the place of $\psi$. By doing so, the behavior of SRMs with varying $sc$ can be predicted from only the model parameters of pure rock and pure soil. Thus, SRMs with different $sc$ should no longer be treated as different materials and only require own $sc$-specific model parameters.

## 5. Model performance

In this section, the experimental data obtained by Wang et al. (2022) [29] are adopted to test the performance of the constitutive model introduced above. In reference [29], conventional drained triaxial tests on soil-rock mixtures with $sc=0, 0.1, 0.3, 0.5, 0.7$ and 1 were conducted under three different confining pressures ($\sigma_3=150$ kPa, 300 kPa and 600 kPa).

The threshold soil content ($sc$)$_{th}$ of the SRM used in [26] is 0.54. The threshold soil content is larger than that reported by other researchers. The reason may be the use of angular gravel in their tests, which cause large voids among rock grains. Therefore, SRMs with $sc=0.1, 0.3$ and 0.5 are rock-dominated SRMs, whose behavior can be predicted from model parameters of pure rock, while stress-strain responses of soil-dominated SRM (i.e., $sc=0.7$) can be predicted from the model parameters of pure soil. The calibrated model parameters of the pure rock and pure soil used in [26] are listed in Table 1. The description of the calibration process of model parameters is reported in Appendix.

<table><caption>Table 1. Model parameters calibrated for pure rock ($sc=0$) and pure soil ($sc=1$).</caption><thead><tr><th>Elastic parameters</th><th>Critical state parameters</th><th>Dilatancy parameters</th><th>Hardening parameters</th></tr></thead></table>
<table><tbody><tr><td rowspan="4">Pure rock (sc=0)</td><td>G<sub>0</sub> = 200</td><td>M = 1.79</td><td rowspan="4">d<sub>0</sub> = 1.5<br/>m = -1.0</td><td rowspan="4">h<sub>1</sub> = 0.2<br/>h<sub>2</sub> = 0.3<br/>n = 3.9</td></tr><tr><td rowspan="3">v = 0.32</td><td>e<sub>Γ</sub> = 0.689</td></tr><tr><td>λ = 0.048</td></tr><tr><td>ξ = 0.70</td></tr><tr><td rowspan="4">Pure soil (sc=1)</td><td>G<sub>0</sub> = 50</td><td>M = 1.70</td><td rowspan="4">d<sub>0</sub> = 1.8<br/>m = 0.2</td><td rowspan="4">h<sub>1</sub> = 0.6<br/>h<sub>2</sub> = 2.0<br/>n = 7.0</td></tr><tr><td rowspan="3">v = 0.25</td><td>e<sub>Γ</sub> = 0.543</td></tr><tr><td>λ = 0.032</td></tr><tr><td>ξ = 0.70</td></tr></tbody></table>

## 5.1 Stress-strain-volume behavior

Figures 5 and 6 display comparisons between predicted results and experimental results in terms of stress ratios, volumetric strains and evolutions of global void ratios for SRMs with varying *sc*. It can be seen that the stress ratio and the volumetric strain can be fairly described by the model for drained shear responses of SRMs under a range of confining stresses and soil contents.

SRMs exhibit post peak strain softening and volumetric expansion at a low confining pressure i.e., $\sigma_3=150$ kPa. While at higher confining pressures, i.e., $\sigma_3=300$ kPa and 600 kPa, SRMs present strain hardening and the volumetric contraction. As shown in Figure 5, the model can capture both the stress and the volumetric strain behaviors of SRMs with varying *sc*, e.g., the strain hardening, the volumetric contraction, the strain softening, and the volumetric expansion. In addition, the predicted peak stress ratio, critical stress ratio, phase transformation point from contraction to dilatancy, ultimate volumetric strain, and evolution of void ratio agree well with experimental results.

Some discrepancies between simulated and experimental curves may be observed when *sc* is close to (*sc*)<sub>th</sub>, i.e., when *sc*=0.3 and 0.5. One possible reason for this phenomenon is that when *sc* is close to (*sc*)<sub>th</sub> a dual skeleton structure is formed in SRM which is composed of both rock and
soil grains. Therefore, some discrepancies may occur if we still regard the SRMs as pure

rock-skeleton structure or soil-skeleton structure.

<figure><img src="image_8.png" /><figcaption>Fig. 5 Comparison between predicted results and experimental results in terms of stress ratio versus axial strain and volumetric strain versus axial strain.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>(a) sc=0.1</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) sc=0.3</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) sc=0.5</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) sc=0.7</figcaption></figure>

<caption>Fig. 6 Comparison between predicted results and experimental results in terms of void ratio versus mean pressure</caption>

## 5.2 Mobilized friction angle

In this section, we assess whether our model can properly reflect the mobilized friction angle of SRMs with different sc and confining pressures. The mobilized friction angle $\varphi_m$ is defined as:

$$ \sin\varphi_m = \frac{\sigma_1 - \sigma_3}{\sigma_1 + \sigma_3} \qquad (21) $$

where $\sigma_1$ and $\sigma_3$ are major and minor principal stresses, respectively. In axisymmetric conditions, $\sin\varphi_m = \frac{3\eta}{6+\eta}$.

Figure 7 shows the comparisons between test results and model predictions in terms of the
relationship between mobilized friction angle $\varphi_m$ and void ratio $e$. It can be seen that the model

captures well the variations of $\varphi_m$ with $e$ for SRMs under different confining pressures, although

the model slightly overestimates $\varphi_m$ at high confining pressure of 600 kPa when $sc=0.3$ and 0.5.

This might be due to the formation of a dual skeleton structure in the specimen with both rock and

soil grains when $sc=0.3$ and 0.5 (close ($sc_{th}$ = 0.54).

<figure><img src="image_13.png" /><figcaption>Fig. 7 Comparison between predicted results and experimental results in terms of void ratio versus mobilized friction angle $\varphi_m$.</figcaption></figure>
## 6. Effect of soil content on essential behavior of SRMs

### 6.1 Effect of sc on the non-associate behavior

In order to investigate the effect of sc on the non-associated behavior of rock-dominated SRMs, the yield surface normal direction $\alpha$ and flow rule direction $\beta$ can be introduced as illustrated in Figure 8 in the p-q plane and dual $d\epsilon_v-d\epsilon_d$ plane for a non-associated flow rule with on a yield surface corresponding to Mohr-Coulomb criterion.

$$ \tan\alpha = \frac{q}{p} \qquad (22) $$

$$ \tan\beta = -\frac{d\epsilon_v}{d\epsilon_d} \qquad (23) $$

<figure><img src="image_14.png" /><figcaption>Fig. 8 Schematic diagram of normal to yield surface and flow rule direction. Both directions can be formulated with respect to angles $\alpha$ and $\beta$. The sign convention of soil mechanics is applied (positive stress in compression and positive strain in contraction).</figcaption></figure>

Knowing the model yield function from equation (12), the normal direction to the yield surface $\alpha$ and flow rule direction $\beta$ are shown in Figure 9 for rock-dominated SRMs ($sc=0, 0.1$
and 0.3) under confining pressures of 150 kPa, 300 kPa and 600 kPa. It can be seen that for each confining pressure, when compared at the same stress ratio η, α of SRMs with different sc are the same, indicating a non-dependence of the normal yield surface direction upon sc. On the other hand, β increases with sc (positive for dilation and negative for contraction), indicating that SRMs exhibit more dilantancy with the increase of sc. Therefore, the angle between normal direction of yield and flow rule direction decreases with sc, meaning the non-associate character is less pronounced with sc.

<figure><img src="image_15.png" /><figcaption>(a) $\sigma_3=150$ kPa</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) $\sigma_3=300$ kPa</figcaption></figure>

<figure><img src="image_17.png" /></figure>
(c) $\sigma_3=600$ kPa

**Fig. 9 Normal to yield surface and flow rule direction of rock-dominated SRMs with sc=0, 0.1 and 0.3 under confining pressure of 150 kPa, 300 kPa and 600 kPa.**

The mesoscale origin of plastic deformation in SRMs is illustrated in Figure 10. The macroscopic activation of the plastic behavior corresponds to substantial grain rearrangements resulting from the collapse of preexisting force chains oriented in the principal stress direction (the vertical direction in Figure 10). Once they collapse, the specimen shrinks in this direction together with smaller expansion in the lateral direction. When small grains fill the pore space, the vertical contraction decreases, whereas the lateral expansion is mostly unaffected. Therefore, the increase in soil content is expected to enhance dilatancy in SRMs specimens.

<figure><img src="image_18.png" /><figcaption>**Fig. 10 Illustration of the plastic deformation with and without soil grains [29]**</figcaption></figure>

The conceptual model of Figure 10 was proved with DEM results by Wang et al. (2021) [29] and Wautier et al. (2019) [32]. Wang et al. (2021) [29] found that for rock-dominated SRMs, *sc* does not affect the normal direction of yield surface but it changes the flow rule direction. For rock-dominated SRMs, the non-associated character is less pronounced with the increase of *sc*, as
290 shown in Figure 11. It should be noted that the normal directions of yield and plastic potential
291 surfaces shown in Fig. 11 are computed directly from DEM results, while the normal to yield
292 surface in Figure 9 depends on the expression of yield surface selected (e.g., Equation (12)).

293 Similarly, in the proposed model, soil content is found to affect the flow rule direction but the
294 mechanical state and the yield surface of rock-dominated SRMs are not affected. Indeed, this is an
295 intrinsic property of the elasto-plastic model of Li and Dafalias (2000) [9] used in this study. The
296 model assumes that the yield surface follows the stress state at any time (see Equation 12 where
297 the current stress ratio is used). Consequently, the yield surface is independent of the soil content.
298 The consistency between model predictions and DEM simulations proves the ability of this
299 constitutive model to properly reflect the non-associated properties of SRMs with sc.

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

300
301 **Fig. 11** Schematic diagram of yield surface *f* and plastic potential surface *g* (on the left) and their
302 normal directions **m** and **n** (on the right) for DEM specimens with different *sc* at the same stress
303 ratio η=0.43 (Wang et al., 2021 [29]).

304
6.2 Effect of sc on peak and critical state friction angle

305
As can be seen in Figure 7, φₘ increases to a peak value φₚ, i.e., peak friction angle, and
then decreases to the critical state friction angle $\phi_{cr}$ at relatively low confining pressures. Values of $\phi_p$ and $\phi_{cr}$ of SRMs with varying $sc$ are displayed in Figure 12. It can be seen that $\phi_{cr}$ almost keeps constant for both rock-dominated structure (at low $sc$) and soil-dominated structure (at high $sc$), while $\phi_p$ decreases with $sc$ for rock-dominated structure and increases with $sc$ for soil-dominated structure. For rock-dominated structures, as a rock skeleton exists, $\phi_{cr}$ is primarily governed by intergranular friction between rock grains. Therefore, $\phi_{cr}$ is constant for all rock-dominated SRMs. Likewise, for soil-dominated structures, a soil skeleton exists, and $\phi_{cr}$ is primarily governed by friction characteristics of soil grains. Therefore, $\phi_{cr}$ is also constant for all soil-dominated SRMs. The reason why $\phi_{cr}$ in rock-dominated structures is higher than that in soil-dominated structures is that the rock grains used in triaxial tests are more angular than soil grains. However, $\phi_p$ is related to the density (in particular $e_{sk}$) of the initial SRM samples. It is reported in [29] that $e_{sk}$ decreases with $sc$ for rock-dominated SRMs and increases with $sc$ for soil-dominated SRMs. This explains why $\phi_p$ decreases with $sc$ for rock-dominated SRMs and increases with $sc$ for soil-dominated SRMs.

<figure><img src="image_20.png" /><figcaption>Fig. 12 Changes in peak friction angle $\phi_p$ and critical friction angle $\phi_{cr}$ with $sc$</figcaption></figure>

($\sigma_3=150\text{kPa}$)
323 **7. Closure remarks**

324 The main contribution of this work is the proposal of a simple method to predict the
325 stress-strain responses of SRMs with varying soil contents. This method incorporates the concept
326 of an updated skeleton void ratio. Using this concept, SRMs with different sc should no longer be
327 regarded as different materials with their own sc-dependent model parameters. Only model
328 parameters of pure rock and of pure soil are required to describe the stress-strain response SRMs
329 with varying sc. Our proposal to adopt the skeleton void ratio is generic and can be applied to
330 other constitutive frameworks, ranging from phenomenological models to micro-mechanically
331 based models provided that skeleton parameters (for instance, esk and ψsk) are adopted. In this
332 manuscript, the constitutive framework proposed by Li and Dafalias (2000) [9] is adopted as an
333 example to demonstrate the effectiveness of this method. Extending this investigation toward
334 other constitutive frameworks will be considered in future work.

335 The proposed method demonstrates a satisfying ability for predicting stress-strain responses
336 of SRMs with different soil contents and confining pressures. In addition, it successfully reflects
337 the non-associativity behavior of rock-dominated SRMs. Some discrepancies between simulated
338 and experimental curves are observed when sc is close to (sc)th. One possible reason for this
339 phenomenon is that a dual skeleton structure composed of both rock and soil grains is formed in
340 SRM when sc is close to (sc)th. This could be improved in future work by characterizing the
341 complex dual skeleton structure in DEM simulations or X-ray tomography images for instance.

342 **Acknowledgements**

343 This work was supported by National Natural Science Foundation of China (Grant No.
344 U1765205 and No. 51979091). The authors also express their sincere thanks to the International
345 Research Network GeoMech (IRN GeoMech, CNRS) for promoting stimulating and convivial

346 interactions among researchers.

347 References

348 1 Been K, Jefferies MG (1985) A state parameter for soils. *Geotechnique* 35: 99-112

349 2 Casagrande A (1936) Characteristics of cohesionless soils affecting the stability of slopes and earth fills. J. Boston Society of Civil Engineers 23:13-32

351 3 Cen D, Huang D, Ren F (2017) Shear deformation and strength of the interphase between the soil-rock mixture and the benched bedrock slope surface. Acta Geotechnica 12: 391-413

353 4 Chen GX, Wu Q, Sun T, et al (2021) Cyclic behaviors of saturated sand-gravel mixtures under undrained cyclic triaxial loading. *Journal of Earthquake Engineering* 25:756-789

355 5 Dong H, Peng BC, Gao QF, et al (2021) Study of hidden factors affecting the mechanical behavior of soil-rock mixtures based on abstraction idea. *Acta Geotechnica* 16:595-611

357 6 Gao WW, Gao W, Hu RL, et al (2018) Microtremor survey and stability analysis of a soil-rock mixture landslide: a case study in Baidian town, China. *Landslides* 15:1951-1961

359 7 Gobbi, S., Santisi d'Avila, M. P., Lenti, L., Semblat, J. F., & Reiffsteck, P. (2022). Effect of active plastic fine fraction on undrained behavior of binary granular mixtures. *International Journal of Geomechanics*, 22(1), 06021035.

362 8 Khalili A, Wijewickreme D, Wilson GW (2010) Mechanical response of highly gap-graded mixtures of waste rock and tailings. Part I: Monotonic shear response. *Canadian Geotechnical Journal* 47:52-565

365 9 Li XS, Dafalias YF (2000) Dilatancy for cohesionless soils. *Geotechnique* 50: 449-460

366 10 Li XS, Wang Y (1998) Linear representation of steady-state line for soil. *Journal of geotechnical and geoenvironmental engineering* 124:1215-121

368 11 Liu HB, Zou DG (2013) Associated generalized plasticity framework for modeling gravelly soils considering particle breakage. *Journal of Engineering Mechanics* 139:606-615

370 12 Lu Y, Liu SH, Zhang YG, et al (2021) Hydraulic conductivity of gravelly soils with various coarse particle contents subjected to freeze-thaw cycles. *Journal of Hydrology* 598:126302

372 13 Manzari MT, Dafalias YF (1997) A critical state two-surface plasticity model for sands. *Géotechnique* 47:255-272

374 14 Napoli ML, Barbero M, Scavia C (2021) Tunneling in heterogeneous rock masses with a block-in-matrix fabric. *International Journal of Rock Mechanics and Mining Sciences* 138:104655.

377 15 Ng TT, Zhou W, Chang XL (2017) Effect of particle shape and fine content on the behavior of
378 binary mixture. *Journal of Engineering Mechanics* 143:C4016008

379 16 Phan QT, Bui HH, Nguyen GD, & Bouazza A (2021) Effect of particle rolling resistance on drained and undrained behaviour of silty sand. *Acta Geotechnica* 1-26

381 17 Porcino DD, Diano V, Triantafyllidis T, & Wichtmann T (2020) Predicting undrained static response of sand with non-plastic fines in terms of equivalent granular state parameter. *Acta Geotechnica* 15:867-882

384 18 Rahman, M. M., & Dafalias, Y. F. (2022). Modelling undrained behaviour of sand with fines and fabric anisotropy. *Acta Geotechnica*, 17(6), 2305-2324.

386 19 Richard FE (1970) Vibrations of soils and foundations. Prentice Hall. Inc., Englewood Cliffs, New Jersey

388 20 Roscoe KH, Schofield A, Wroth AP (1958) On the yielding of soils. *Geotechnique* 8:22-53

389 21 Saberi M, Annan CD, Konrad JM (2017) Constitutive modeling of gravelly soil-structure interface considering particle breakage. *Journal of Engineering Mechanics* 143:04017044

391 22 Schofield AN, Wroth P (1968) Critical state soil mechanics. London: McGraw-hill Vol. 310

392 23 Shi XS, Liu K, Yin J (2021) Analysis of mobilized stress ratio of gap-graded granular materials in direct shear state considering coarse fraction effect. *Acta Geotechnica* 6:1801-1814

395 24 Shi XS, Zhao JD, Gao YF (2021) A homogenization-based state-dependent model for gap-graded granular materials with fine-dominated structure. *International Journal for Numerical and Analytical Methods in Geomechanics* 45:1007-1028

398 25 Shi XS, Zhao JD, Yin JH, et al (2019) An elastoplastic model for gap-graded soils based on homogenization theory. *International Journal of Solids and Structures* 163:1-14

400 26 Sun, Z., Chu, J., & Xiao, Y. (2021). Formulation and implementation of an elastoplastic constitutive model for sand-fines mixtures. *International Journal for Numerical and Analytical Methods in Geomechanics*, 45(18), 2682-2708.

403 27 Thevanayagam S (2007) Intergrain contact density indices for granular mixes—I: Framework. *Earthquake engineering and engineering vibration* 6:123

405 28 Thevanayagam S, Mohan S (2000) Intergranular state variables and stress-strain behaviour of silty soils. *Geotechnique* 50:1-23

407 29 Wang T, Liu SH, Wautier A, Nicot F (2022) An updated skeleton void ratio for gravelly sand mixtures considering the effect of grain size distribution. *Canadian Geotechnical Journal*, 59(1): 12-23.

410 30 Wang T, Wautier A, Liu SH, Nicot F (2021) How fines content affects granular plasticity of under-filled binary mixtures, *Acta Geotechnica*, 2021: 1-15.

412 31 Wang ZL, Dafalias YF, Shen CK (1990) Bounding surface hypoplasticity model for sand. *Journal of engineering mechanics* 116: 983-1001
414 32 Wautier A, Bonelli S, Nicot F (2019) Rattlers' contribution to granular plasticity and
415 mechanical stability. *International Journal of Plasticity* 112:172-193

416 33 Vallejo LE (2021) Interpretation of the limits in shear strength in binary granular mixtures.
417 *Canadian Geotechnical Journal* 38:1097-1104

418 34 Xiao Y, Sun YF, Hanif KF (2015) A particle-breakage critical state model for rockfill material.
419 *Science China Technological Sciences* 58:1125-1136

420 35 Xiao Y, Sun YF, Liu HL, et al (2017) Model predictions for behaviors of sand-nonplastic-fines
421 mixtures using equivalent-skeleton void-ratio state index. *Science China Technological
422 Sciences* 60: 878-892

423 36 Xiong H, Yin Z Y, Zhao J, et al (2021) Investigating the effect of flow direction on suffusion
424 and its impacts on gap-graded granular soils. *Acta Geotechnica* 16: 399-419

425 37 Xu WJ, Qiang X, Hu RL (2011) Study on the shear strength of soil-rock mixture by large
426 scale direct shear test. *International Journal of Rock Mechanics and Mining Sciences*
427 48:1235-1247

428 38 Yang YT, Sun GH, Zheng H, et al (2019) Investigation of the sequential excavation of a
429 soil-rock-mixture slope using the numerical manifold method[J]. *Engineering Geology* 256:
430 93-109

431 39 Yin ZY, Huang HW, Hicher PY (2016) Elastoplastic modeling of sand-silt mixtures. *Soils and
432 Foundations* 56:520-532

433 40 Yin ZY, Zhao JD, Hicher PY (2014). A micromechanics-based model for sand-silt mixtures.
434 *International journal of solids and structures* 51:1350-1363

435 41 Zhang G, Zhang JM, Yu Y (2007) Modeling of gravelly soil with multiple lithologic
436 components and its application. *Soils and Foundations* 47:799-810

437 42 Zhao LH, Huang DL, Zhang SH, et al (2021) A new method for constructing finite difference
438 model of soil-rock mixture slope and its stability analysis. *International Journal of Rock
439 Mechanics and Mining Sciences* 138: 104605

440

441 **Appendix: Calibration of model parameters**

442 **(1) Elastic parameters**

443 The initial elastic shear modulus, *G*, can be obtained from the experimental data of deviatoric
444 stress, *q*, versus the deviatoric strain, *ε*<sub>s</sub> when the axial strain is lower than 0.2%. Rearrangement

445 of Equation (10) gives:
$$G_0 = G \frac{(1+e)}{(2.97-e)^2 \sqrt{P'P_a}} \qquad (A1)$$

The values of the elastic constant, $G_0$, at various confining pressures, can be determined from Equation (A1). The average value of $G_0$ under different confining pressures is adopted.

Based on Equations (7)-(11), the Poisson's ratio, $\nu$ can be obtained by:

$$\nu = \frac{9d\epsilon_s^e - 2d\epsilon_v^e}{18d\epsilon_s^e + 2d\epsilon_v^e} \approx \frac{9\epsilon_s - 2\epsilon_v}{18\epsilon_s + 2\epsilon_v} \qquad (A2)$$

## (2) Critical state parameters

$e_Γ$, $\lambda$ and $\xi$ can be determined by directly fitting the experimental data for the critical state line. The critical state stress ratio $M$ can be obtained by fitting critical state test data in $p'-q$ plane with a function of $q=Mp'$.

## (3) Dilatancy parameters

The parameter $m$ is determined from Equation (17) at a phase transformation state, at which $D=0$, and thus,

$$m = \frac{1}{\psi^d} \ln \frac{M^d}{M} \qquad (A3)$$

where $\psi^d$ and $M^d$ are the values of $\psi$ and $\eta$ at the phase transformation state.

Ignoring the small elastic strain, we have:

$$\frac{d\epsilon_v}{d\epsilon_q} \approx \frac{d\epsilon_v^p}{d\epsilon_q^p} = D = d_0 \left( \exp(m\psi) - \frac{\eta}{M} \right) \qquad (A4)$$

The parameter $d_0$ is determined based on the $d\epsilon_v-d\epsilon_q$ curve.

## (4) Hardening parameters

The parameter $n$ is determined by Equation (14) at a peak stress state, at which $K_p=0$:
$$n = \frac{1}{\psi^b} \ln \frac{M}{M^b} \qquad (A5)$$

461 where $\psi^b$ and $M^b$ are the values of $\psi$ and $\eta$ at the peak stress state.

462 Combining Equations (10) (13) and (14) for conventional drained tests ($dp'=dq/3$) yields:

$$\frac{dq}{d\epsilon_q} \approx \frac{dq}{d\epsilon_q^p} = \frac{K_p}{1 - \eta/3} = h \left\{ \frac{G_0(2.97 - e)^2 \sqrt{p'p_a} \left( \frac{M}{\eta} - \exp(n\psi) \right)}{(1 + e)(1 - \eta/3)} \right\} \qquad (A6)$$

463 As all the model parameters in the brackets are known, $h$ is determined based on $dq - d\epsilon_q$

464 curves along drained triaxial loading paths. Then parameters $h_1$ and $h_2$ can be obtained by

465 equation $h = h_1 - h_2 e_0$.