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

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 modeling 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.


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], stressdilatancy [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 intergranular friction between rock grains. While at a high sc beyond a threshold value sc ð Þ th , the mechanical behavior of SRM is primarily governed by friction characteristics of soil grains.
Although some pioneering works have been done on the constitutive modeling 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 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 variables. 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 stressbearing 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: 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: where C ur ¼ ðD 60 Þ rock =ðD 10 Þ rock is the coefficient of non-uniformity for the rock fraction,C us ¼ ðd 60 Þ soil =ðd 10 Þ soil is the coefficient of non-uniformity for the soil fraction, and R d ¼ ð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:  [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) is equal to the skeleton void ratio of the SRM. This important finding is essential to unify descriptions of SRMs with varying soil contents.

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 Þ n plane: where e cr is critical state void ratio, e C is the theoretical critical void ratio at the atmospheric pressure, p a is the atmospheric pressure, k is the magnitude of the slope, and n 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 (r 3 = 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 (C ur ¼ 2:36; C us ¼ 7:01 and R d =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 moves 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 C and k are shown in Fig. 2. e C and k 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.
Critical state data of SRMs with different sc are plotted in e sk -p plane in Fig. 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 (Fig. 3a) and soil-dominated structure (Fig. 3b). It can be seen in Fig. 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 Fig. 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.

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 stressstrain 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.

Elastic behavior
According to standard elasto-plasticity, a total strain increment de is additively split into an elastic strain increment de e and a plastic strain increment de p : The volumetric and deviatoric elastic strain increments are given, respectively, by: where K is the elastic bulk modulus, G is the elastic shear modulus, p0 ¼ r 1 þ 2r 3 =3 is the mean effective stress, and q ¼ r 1 À r 3 is the deviatoric stress (r 1 , r 2 , and r 3 are the principal stress values). G can be estimated following [19]: where G 0 is an elastic material constant, p a is the atmospheric pressure, and e is the void of the considered material (soil or rock). The elastic bulk moduli K is given by the following relation: where m is the Poison's ratio. The estimation of G 0 and m can be found in Appendix.

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: where g 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 g c ¼ g. 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 point [31] which are not considered in the present study.
According to the theory of plasticity, a loading index dL can be defined as: where K p is the plastic hardening modulus and is expressed by Li and Dafalias (2000) [9] as: In the above equation, n is a positive model parameter, h ¼ h 1 À h 2 e 0 with h 1 and h 2 are model parameters and e 0 the initial void ratio, M is the stress ratio at a critical state and w ¼ e À e cr is the state parameter defined by Been and Jefferies (1985) [1] with e the current void ratio and e cr the critical state void ratio for the current p'. When using skeleton void ratio instead of the void ratio, a skeleton state parameter w sk is introduced as w sk ¼ e sk À ðe cr Þ sk , as seen in Fig. 4. Note that when softening occurs (K p \0) the plastic multiplier dL is positive (i.e. plasticity is active) when the stress ratio decreases.
Since soil-rock mixtures are typical non-associated materials, the non-associated flow rule is adopted to define the plastic strain increments, as follows: where de p d and de p v are the plastic deviatoric strain increment and the plastic volumetric strain increment, respectively. The ratio D ¼ de p v = de p d evaluates the amplitude of the dilatancy and can be chosen as follows: where d 0 and m are model parameters. For associated flow rule, Therefore, for dL [ 0, the following expressions can be derived for both incremental deviatoric and volumetric strains: Equations (18) and (19) set up the relationship between stain and stress increments. Eventually, an elasto-plastic constitutive relation is derived in a matrix form: 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 Eq. (20) is given in the plane dp 0 ; dq À Á and the dual plane de v ; de q À Á . It can be generalized to any kind of incremental stress and strain tensors by assuming coaxiality between incremental stress and strain.

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 m; (2) four critical state parameters, i.e. M; e C ; k and n; (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 is adopted to replace global void ratio e in the above-mentioned equations. Accordingly, skeleton state parameter w sk is used in the place of w. 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 scspecific model parameters.

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], Pure soil (sc = 1) 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 rockdominated 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. SRMs exhibit post peak strain softening and volumetric expansion at a low confining pressure i.e., r 3 = 150 kPa. While at higher confining pressures, i.e., r 3 = 300 kPa and 600 kPa, SRMs present strain hardening and the volumetric contraction. As shown in Fig. 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.

Stress-strain-volume behavior
Some discrepancies between simulated and experimental curves may be observed when sc is close to (sc) th , i.e., when sc = 0.3 and 0.5. One possible reason for this phenomenon is that when sc is close to (sc) th , 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.

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 u m is defined as: where r 1 and r 3 are major and minor principal stresses, respectively. In axisymmetric conditions, sinu m ¼ 3g=6 þ g. Figure 7 shows the comparisons between test results and model predictions in terms of the relationship between mobilized friction angle u m and void ratio e. It can be seen that the model captures well the variations of u m with e for SRMs under different confining pressures, although the model slightly overestimates u 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).  6 Effect of soil content on essential behavior of SRMs

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 a and flow rule direction b can be introduced as illustrated in Fig. 8 in the p-q plane and dual de v -de d plane for a non-associated flow rule with on a yield surface corresponding to Mohr-Coulomb criterion.
Knowing the model yield function from Eq. (12), the normal direction to the yield surface a and flow rule direction b are shown in Fig. 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 g, a of SRMs with different sc are the same, indicating a nondependence of the normal yield surface direction upon sc. On the other hand, b increases with sc (positive for dilation and negative for contraction), indicating that SRMs exhibit more dilantancy with the increase in 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. The mesoscale origin of plastic deformation in SRMs is illustrated in Fig. 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 Fig. 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.
The conceptual model of Fig. 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 rockdominated SRMs, the non-associated character is less pronounced with the increase in sc, as shown in Fig. 11. It should be noted that the normal directions of yield and plastic potential surfaces shown in Fig. 11 are computed directly from DEM results, while the normal to yield surface in Fig. 9 depends on the expression of yield surface selected (e.g., Eq. (12)).
Similarly, in the proposed model, soil content is found to affect the flow rule direction, but the mechanical state and the yield surface of rock-dominated SRMs are not affected. Indeed, this is an intrinsic property of the elasto-plastic model of Li and Dafalias (2000) [9] used in this study. The model assumes that the yield surface follows the stress state at any time (see Eq. 12 where the current stress ratio is used). Consequently, the yield surface is independent of the soil content. The consistency between model predictions and DEM simulations proves the ability of this constitutive model to properly reflect the non-associated properties of SRMs with sc.

Effect of sc on peak and critical state friction angle
As can be seen in Fig. 7, u m increases to a peak value u p , i.e., peak friction angle, and then decreases to the critical state friction angle u cr at relatively low confining pressures. Values of u p and u cr of SRMs with varying sc are displayed in Fig. 12. It can be seen that u cr almost keeps constant for both rock-dominated structure (at low sc) and soil-dominated structure (at high sc), while u 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, u cr is primarily governed by intergranular friction between rock grains. Therefore, u cr is constant for all rock-dominated SRMs. Likewise, for soil-dominated structures, a soil skeleton exists, and u cr is primarily governed by friction characteristics of soil grains. Therefore, u cr is also constant for all soil-dominated  SRMs. The reason why u 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, u 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 u p decreases with sc for rock-dominated SRMs and increases with sc for soil-dominated SRMs.

Closure remarks
The main contribution of this work is the proposal of a simple method to predict the stress-strain responses of SRMs with varying soil contents. This method incorporates the concept of an updated skeleton void ratio. Using this concept, SRMs with different sc should no longer be regarded as different materials with their own sc-dependent model parameters. Only model parameters of pure rock and of pure soil are required to describe the stress-strain response SRMs with varying sc. Our proposal to adopt the skeleton void ratio is generic and can be applied to other constitutive frameworks, ranging from phenomenological models to micro-mechanically based models provided that skeleton parameters (for instance, e sk and w sk ) are adopted. In this manuscript, the constitutive framework proposed by Li and Dafalias (2000) [9] is adopted as an example to demonstrate the effectiveness of this method. Extending this investigation toward other constitutive frameworks will be considered in future work. The proposed method demonstrates a satisfying ability for predicting stress-strain responses of SRMs with different soil contents and confining pressures. In addition, it successfully reflects the non-associativity behavior of rockdominated SRMs. Some discrepancies between simulated and experimental curves are observed when sc is close to (sc) th . One possible reason for this phenomenon is that a dual skeleton structure composed of both rock and soil grains is formed in SRM when sc is close to (sc) th . This could be improved in future work by characterizing the complex dual skeleton structure in DEM simulations or X-ray tomography images for instance. As all the model parameters in the brackets are known, h is determined based on dq À de q curves along drained triaxial loading paths. Then parameters h 1 and h 2 can be obtained by equation h ¼ h 1 À h 2 e 0 .