Open access peer-reviewed chapter

A State-Dependent Constitutive Model for Unsaturated Rockfill Materials

Written By

Liujiang Wang and Zhongzhi Fu

Submitted: September 26th, 2019 Reviewed: May 21st, 2020 Published: September 18th, 2020

DOI: 10.5772/intechopen.92902

Chapter metrics overview

342 Chapter Downloads

View Full Metrics


This chapter presents a state-dependent elastoplastic constitutive model for both saturated and unsaturated rockfill materials. The model, which is developed within an extended critical-state framework, uses two independent stress state variables: total stress and total suction. The loading-collapse (LC) curve proposed by Oldecop and Alonso for unsaturated rockfills is used herein. A unified hardening parameter, which could consider the effects of stress level, internal state (density) and relative humidity, is introduced to describe the state-dependent dilatancy of saturated and unsaturated rockfill materials. The details of the model formulation and parameters determination are described and reported. Numerical simulations on the triaxial tests, such as the drained shear tests on the saturated specimens with different initial dry densities, shear tests on the specimens with different relative humidity and wetting deformation tests under constant vertical strain rate, have been carried out using the proposed model. The numerical results show that the stress-strain relationships at both loose and dense, saturated and unsaturated states can be properly modelled with a single set of parameters. Additionally, the proposed model can also capture some other key features such as the strain-softening behaviour at the dense state and low confining stress, the sudden stress relaxing subjected to the flooding under a constant vertical strain.


  • constitutive relations
  • unsaturated rockfills
  • collapse deformation
  • initial density

1. Introduction

It is well recognised from the experimental studies and engineering practice that the influence of water on the mechanical behaviour of rockfill materials is significant. A lot of laboratory experiments carried out in the past showed that the rockfill materials may undergo significant amounts of strain upon flooding [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. This phenomenon is usually called ‘collapse deformation’, referring to a strain increment not related to the load increment. On the other hand, a remarkable influence of the compacted density and confining pressure on the strength and deformation behaviours of rockfill materials was observed in some large-scale triaxial tests [11, 12, 13, 14, 15]. The experimental results on rockfill materials at a dense state showed the remarkable strain-softening and dilatancy behaviours, while the rockfill materials at a loose state exhibited the strain hardening and volumetric contraction behaviours. This state-dependent behaviour was also observed in the triaxial tests on rockfill materials used in Guanyinyan dam that will be introduced in the next section. Furthermore, compared with the fully saturated rockfill materials, a more obvious state-dependent behaviour was observed for unsaturated rockfills in the relative humidity-controlled triaxial tests [7, 16]. The rockfills subjected to the larger suction exhibited more remarkable strain-softening and dilatancy at a low confining pressure or with a small initial void ratio. Thus, presenting a constitutive model which could consider the collapse deformation and state-dependent behaviour for unsaturated rockfill materials is necessary.

Rockfill materials have been widely used in the constructions of rockfill dams, and they should qualitatively be defined very well before beginning to the design stage [17, 18, 19, 20]. The collapse settlement of upstream shell for clay core rockfill dam occurs during the impounding, which has attracted extensive attentions from both engineers and scientists since 1970s [3], for it may cause cracks on the dam crest and result in damage to impermeable systems. Thereby, proposing a satisfactory technique for the prediction of collapse settlement is clearly important. Nobari and Duncan [3] first introduced collapse effects into finite element (FE) analysis to predict collapse settlement in rockfill dams. This calculation starts by performing the FE analysis using a set of material parameters corresponding to the dry conditions; then the collapse is numerically simulated in two stages. First, the stress change caused by the saturation at the constant strain is determined using the experimental data. In the second stage, the nodal forces that restore the equilibrium are applied, and the computed strains reproduce the collapse strains. Naylor et al. [21] generalised this method for arbitrary constitutive models, which required the knowledge of two sets of constitutive parameters for dry and saturated conditions. However, numerous experiments have invalidated the implicit assumption that the sequence of loading and wetting does not affect the final state of materials and led to several attempts to establish other approaches [22, 23, 24, 25]. Shen and Wang [22] expressed the magnitudes of the wetting-induced volumetric strain and deviatoric strain using empirical functions, including both confining pressure and deviatoric stress level during the wetting process. By virtue of the assumption of coaxiality of the strain rate tensor and the stress tensor, the strain components are estimated and translated into nodal forces, which are further applied to the finite elements to obtain the wetting deformation of the dam. Recently, this so-called single curve method is widely used in dam engineering in China due to its practicality. Nevertheless, the point to underline here, in the above two procedures, is essentially a computational device that is not necessarily associated with a physical reality [26, 27]. The mechanical properties of rockfill materials are closely related to the breakage properties of rock particles [4, 28] and these breakage properties are sometimes remarkably affected by the degree of saturation [1, 10, 29]. Hence, the physical mechanism underlying the collapse deformation seems to be the difference of breakage rates between saturated and unsaturated rock particles. Therefore, introducing the effect of relative humidity on the rockfill performance is by no means an academic one to model the collapse deformation. Oldecop and Alonso [1] first proposed an elastoplastic model that was valid to express compression properties of unsaturated rockfills considering the influence of relative humidity. Thereafter, within the framework of the Barcelona Basic Model [30] for unsaturated soils, they extended this compression model to the three-dimensional stress state and applied it in the analyses of Beliche and Lechago rockfill dams [31, 32]. Furthermore, Chávez and Alonso [7] developed a work-hardening model to describe the effect of particle breakage at different relative humidities. Besides, Kohgo et al. [33] proposed a three-dimensional elastoplastic model for unsaturated rockfills according to a modification of their generalised elastoplastic model for unsaturated soils. Bauer [34] extended his hypoplastic model by introducing a moisture-dependent solid hardness and taking the stress relaxation caused by the change of solid hardness into consideration. The state-dependent behaviour was not considered in these models.

On the other hand, due to the obvious state-dependent behaviour of sands, many state-dependent models have been proposed [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45]. Among them, the state-dependent model proposed by Li and Dafalias [41], with only a set of model parameters, can precisely represent the complex stress-strain relationship of Toyoura sand over a wide range of densities and pressures. Within the framework of this model, Xiao et al. [46] developed a generalised elastoplastic model to describe the state-dependent dilatancy for rockfill materials. Sun et al. [47] introduced the fractional order method to describe the state-dependent behaviour of granular materials without using the plastic potential function. However, to the best of our knowledge, these models cannot describe the combination effect of the initial dense state and prevailing relative humidity on the stress-strain relationship of rockfill materials.

The purpose of this chapter is to present an elastoplastic model for unsaturated rockfill materials by introducing a unified hardening parameter which could well capture the essential physical phenomena behind the observed collapse deformation and state-dependent dilatancy for saturated and unsaturated rockfill materials. The conceptual bases for the model and the basic formulation are given in Section 2. Thereafter, the proposed model is incorporated into the coupled flow-deformation analysis FE code [48, 49, 50, 51] using the explicit stress integration algorithm for the elastoplastic model of unsaturated soils [52, 53]. Finally, the comparisons are conducted using the experimental data from a series of triaxial tests on saturated and unsaturated rockfill materials, and the capability of the proposed model to capture the state-dependent behaviour and collapse deformation is validated.


2. Constitutive framework and formulations

2.1 Basic constitutive variables

Due to the large permeability for rockfills, water never fills the large rockfill voids unless structures are submerged. Environmental changes and rainfall can only modify the relative humidity. The extreme case of rockfill flooding is a concern in the rockfill structures which become inundated, such as the upstream rockfill shells of the zoned earth and rockfill dams. Hence, a suitable stress space to describe the isotropic compression states of unsaturated rockfill materials is (p, RH), where pis the mean total stress and RH is the relative humidity. In fact, the relative humidity is the ratio of the vapour pressure present in the air to the vapour pressure when the air is saturated with water vapour. The relative humidity in the gas phase and the matric suction, s, in the rockfill pore water are related by the psychrometric relationship [52]:


where Ris the gas constant, Tis the absolute temperature of the reference system, Mwis the molecular mass of water and ρwis the density of water at the reference temperature. πis called the osmotic suction, which is due to the presence of solutes in the rockfill pore water. The sum ψ = s + πis called the total suction. Total and matric suction would be equal with each other in the case when the rockfill pores contained only pure water with no solutes, for instance, rockfill materials applied in the dams. The total stress and total suction are chosen as the basic constitutive variables herein.

2.2 Normal compression model

Considering an isotropic compression test on unsaturated rockfills, in which a rockfill sample at a given relative humidity (subsequently maintained) is subjected to p-load increments along virgin states, which is in full correspondence with the compression behaviour of saturated rockfills, the specific volume will be given by


where νis the specific volume, eis the void ratio and Nis the intercept of the normal compression lines with the ν-axis when p = 0. The compression index λ(ψ) is assumed to be a function of the total suction and can be interpolated from the compression indices at the fully saturated state (ψ = 0) and very dry state. For rockfill materials, the very dry state is usually defined as the water content close to zero, generally corresponding to the water content less than 0.45% in practice. A schematic representation of Eq. (1) can be discretized into a series of normal compression lines, as shown in Figure 1.

Figure 1.

Normal compression lines for saturated and unsaturated rockfill materials.

As mentioned by McDowell and Bolton [53], the plastic deformation of the granular materials can be attributed to the two deformation mechanisms. Under low stress, the plastic deformation is due to the particle rearrangement. The second mechanism is called clastic yielding, and it is attained when the applied stress causes the onset of particle breakage. According to this concept and corresponding experimental results, an elastoplastic compressibility model for rockfill was developed by Oldecop and Alonso [1]. Under isotropic compression, a mean threshold total stress, py, which marks the beginning of particle breakage, was introduced in the model. Below this threshold stress, the influences of suction on the compression index were not measured in the compression tests. On the other hand, the effect of water action on the compression behaviour was strong when beyond this threshold stress, which is the source of the collapse deformation. Note that, the difference of compression index below and beyond this threshold stress is vanished when the material is in the very dry state. Therefore, the compression index for unsaturated rockfills can be expressed as


where λi is the slope of the normal compression line when particle rearrangement is active only, and (λi + λd) is the slope of the normal compression line when both plastic deformation mechanisms (particle rearrangement and particle breakage) are active. According to the previous literature [1], the compressibility index accounting for the clastic yielding is linearly related to the logarithm of total suction through a material parameter:


where λ0dis the maximum clastic compressibility index in the saturated condition, patm is the atmospheric pressure and αψ is a model parameter. The elastic volume change for rockfill materials within unloading-reloading (URL) paths are given by


where the slope κis assumed to be independent of the water action. In addition, there was a moderate swelling behaviour observed in the experiments on the unsaturated rockfills, which was caused by the increase in water content. This swelling is assumed to be a reversible elastic deformation, and a linear relationship is assumed between the swelling volume change and the logarithm of total suction:


where κψ is the total suction-based expansion/compression index. For simplicity, κψ is assumed to be independent of the stress level.

2.3 Loading-collapse yield curve

For the isotropic compression condition, the yield stress, p0, was defined by the following expression such as [1]:


Eq. (7), which is referred to as loading-collapse (LC) yield surface, describes the relationship between the equivalent yield stress (p0defines the position of the yield curve, and it was identified as the yield stress for a very dry rockfill) and the yield stress (p0) in unsaturated states, which is controlled by the total suction. By introducing the compressibility index for the clastic yielding expressed in Eq. (4) into Eq. (7), the LC yield surface is then obtained. A plot of the LC yield curves for different values of p0for Pancrudo slate that were tested by Oldecop and Alonso [1] is given in Figure 2. There are two regions in which the yield stress, p0, does not depend on the total suction, which correspond to the stress levels below the clastic yield stress py (p0< 0.29 MPa) and to the very dry states (ψ > 67 MPa), respectively.

Figure 2.

Load-collapse curves for different values ofp0.

Thereafter, a simple volumetric hardening is used here to follow the evolution of p0:


where His the unified hardening parameter related to the plastic volumetric strain and state-dependent behaviour, which will be introduced in Section 2.5.

2.4 Yield surface

Yao et al. [54] developed a simple but robust constitutive model (the original UH model) to reproduce the mechanical behaviours in the isotropic and triaxial states for both normally consolidated (NC) and overconsolidated (OC) soils, such as non-elastic deformation in reloading, peak strength, shear-dilation and strain-softening. In his model, a unified hardening parameter (H) was proposed, and an associated flow rule was adopted for the robust and convenience of numerical implementation. Hence, in accordance with the aim of simplicity, the original UH model is adopted as the base for developing a new state-dependent constitutive model of unsaturated rockfill materials. Thus, the yield and plastic potential surfaces in the stress space (p, qand ψ) are defined as follows:


where qis the deviator stress and Mis the critical-state slope.

With respect to the rockfill materials, the critical-state slope, M, is no longer a constant. Experimental results on rockfill materials show that the critical state friction angle is dependent on the mean stress [19, 55, 56, 57] that decreases with an increase in p. Moreover, according to the RH-controlled triaxial tests on rockfill materials in Chávez and Alonso [7], Mwas also found to be dependent on the RH, and the relationship between Mand logσ3for samples with different RH seemed to be a series of parallel lines. From the experimental data [7, 16], the value of Mincreases with the increase of total suction (decrease RH). Besides, considering the critical stress ratio is impossible to be negative with the increase of mean stress, the following function is adopted for the sake of simplicity:


where Mo and Mres are the initial and residual stress ratios for rockfills in saturated condition under low (0.1 MPa) and high confining pressures, respectively; αM and βψ are material constants. Figure 3 shows that Eq. (10) can well describe the change of the critical stress ratio with mean pressure and relative humidity for unsaturated rockfill materials.

Figure 3.

Variations of critical stress ratio with mean stress and relative humidity: (a) crushed cambric slate and (b) Pancrudo slate.

Using above equations, a three-dimensional view of the yield surfaces in the (p, qand ψ) space is given in Figure 4, where ψi is the total suction corresponding to the very dry state.

Figure 4.

Three-dimensional yield surface in stress space (p,qandψ).

2.5 Unified hardening parameter

The unified hardening parameter was first proposed to provide a unified description of the mechanical behaviours for both clays and sands by Yao et al. [58]. It was then modified and employed to model the hardening process for both normally consolidated and overconsolidated clays [54, 59]. The unified hardening parameter, which is adopted to describe the hardening of the yield surface of rockfill materials can be written as


where ηis the stress ratio, η = q/p; εvpis the plastic volumetric strain; and Mf is the potential failure stress ratio, which represents the potential peak strength that changes with the compaction density and relative humidity. Since dHis always larger than or equal to zero, the following conclusions can be drawn from Eq. (11):

  1. 0 ≤ η < M(negative dilatancy condition): dεvp>0.

  2. η = M(characteristic state condition): dεvp=0.

  3. M ≤ η < Mf(positive dilatancy condition): dεvp<0.

As indicated above, the dilatancy of rockfill materials can be reasonably described by this hardening parameter.

The next is introducing the determination of potential failure stress ratio. For rockfill materials, the peak friction angle, ϕp, on a contact plane is dependent on the degree of interlocking by neighbouring particles, which can be related to the state of the packing void ratio [60]:


where ϕμ is the internal friction angle corresponding to the critical state, which has the relationship with the critical stress ratio Mas shown in Eq. (17); and mis a material constant when granular materials are in the fully saturated condition, but it varies with the confining pressure and suction in unsaturated conditions according to the experimental observations, which can be corrected using the following function:


where σ3 is the confining pressure; m1, m2 and n are material constants. From Eq. (13), we can find that parameter mis kept as m1 in the saturated condition, but decreases with the increasing confining pressure when in the unsaturated condition. Thus, the nonlinear variation of the peak strength with the changes of suction and confining pressure can be well described using this function.

Furthermore, ec in Eq. (12) is the void ratio corresponding to the critical state, which is a function of the mean stress. However, for the rockfills in the unsaturated condition, the influence of the total suction on the critical void ratio should be considered as well. From the experimental data in [7, 16], the relation between critical stress ratio and mean stress for saturated and unsaturated rockfills seems to be two parallel lines. Thus, the critical void ratio can be expressed based as [61].


where ζis the material constant; (eref, pref) is a reference point on the critical state line eref is a function of total suction, which nonlinearly increases with the increase of total suction. For the sake of simplicity, the following function is used to express the variation of reference critical void ratio with total suction:


where eref0means the reference critical void ratio in the saturated condition; and αe is the material constant, which describes the increase rate of reference critical void ratio with total suction.

In the triaxial compression condition, the stress ratios, Mfand M, can be expressed by using peak friction angle (ϕp) and internal friction angle (ϕμ) as


Combining Eqs. (12), (16) and (17), the relation between the potential stress ratio (Mf) and critical state stress ratio (M) can then be obtained.

Substituting Eq. (16) into Eq. (11), the hardening parameter is then obtained, which can reflect the influence of the compacted density and relative humidity on the mechanical behaviour of unsaturated rockfill materials. Considering the peak frictional angle, ϕp, is generally greater than ϕμ, especially in the dense and very dry states, the dilatancy behaviour for the rockfill materials can thus be well captured using the modified hardening parameter. In addition, once an obvious dilatancy occurs for the rockfill in the low confining pressure, the degree of interlocking and the peak frictional angle are reduced, and the strain-softening phenomenon can also be well described by using this hardening parameter.

2.6 Elastic moduli

Considering the nonlinear changes of the elastic moduli with the changes of degree of compaction and stress level, the elastic shear modulus, G, is calculated by using the following empirical Equation [46]:


where G0 is a material constant. In this equation, the current void ratio, e, is used instead of the initial void ratio. The elastic bulk modulus, K, is equal to


where ν is the Possion’s ratio.

2.7 Determination of model parameters

The application of the model requires the information on the following stress states and parameters.

  1. Initial state: initial stresses (pi, qi and si), initial void ratio, e0, and initial mean yield stress, p0.

  2. Parameters directly associated with the isotropic compression behaviour: py, threshold yield mean stress for the onset of clastic phenomena; λi, normal compressibility index when instantaneous deformation mechanism is active only; λ0d, maximum clastic compressibility index for fully saturated conditions; κ, compressibility coefficient along elastic (unloading-reloading) stress paths and αs, compressibility parameter, which controls the rate of stiffness increase with total suction.

  3. Parameters directly associated with the critical state: M0, maximum critical stress ratio for the mean stress approaching to 0.1 MPa; Mres, residual critical stress ratio corresponding to the large mean stress; αM, controls the rate of critical stress ratio decrease with mean stress; βψ, controls the rate of critical stress ratio increase with total suction; (eref0, pref) is a reference point on the critical state line in the saturated conditions; χ, controls the rate of critical void ratio decrease with mean stress and αe, controls the rate of critical void ratio increase with total suction.

  4. Parameters directly associated with elastic moduli: G0, material constant associated with shear modulus within the elastic domain; ν, Possion’s ratio, which is assumed to be independent of suction and taken as 0.2 ∼ 0.3 for rockfill materials and κψ, compressibility coefficient for changes in suction within the elastic region.

  5. Other default parameters: m1, controls the rate of peak strength change with void ratio in the saturated conditions; m2, minimum value of parameter, m, corresponding to the high confining stress at very dry state and n, controls the vary of mwith suction and confining pressure.

In general, the determination of the model parameters will require relative humidity-controlled testing methods, and suggested stress paths for different sets of parameters are the following:

  1. Tests that involve isotropic compression (loading and unloading) for very dry specimen and saturated ones; they provide data to find λi, p0, λ0dand κ.

  2. Tests that involve isotropic compression for very dry specimen subsequent with flooding under low applied stress and constant stress p0; they provide data to find κψ and py. As mentioned in the previous study [1], py can be obtained using the following function:


where εexpansion is the measured total expansion strain due to flooding under low applied stress and εcollapse is the measured total collapse strain due to flooding under applied stress p0.

  1. Tests that involve isotropic compression under different relative humidity from what a linear relationship between collapse strains and the logarithm of rockfill water content can be observed; Oldecop and Alonso [1] provided a method to find αψ using the following function


where χψ is the experimental coefficient relating collapse strain with total suction at constant stress p0 (p0 > py).

  1. Triaxial tests (loading, unloading and reloading) at different relative humidities and saturated condition, and they provide data to find M0, Mres, αM, βψ, eref0, pref, χ, αe and G0.

Note that the default material parameters, m1, m2 and n, should be determined according to the back analysis using the experimental data by fixing the previous determined parameters. Besides, the above-mentioned test is the minimum experimental programme, so the more tests are required to determine the more reliably model parameters.


3. Model validations

In order to verify the capability and feasibility of the proposed constitutive model, it was incorporated into the hydro-mechanical coupling FEM code [48] that was developed in the Institute of Hydraulic Structures of Hohai University [49, 50]. The explicit stress-integration algorithm for unsaturated soils proposed by Sheng et al. [51, 62] was used for the numerical implementation of the proposed constitutive model. The triaxial tests on the saturated specimens with different initial void ratio, triaxial tests with different relative humidity and triaxial wetting tests under constant vertical strain rate were simulated by means of a hexahedral element. The confining pressure was applied against the upper and lateral boundaries, a constant vertical displacement rate was then applied to the upper boundary for shearing and the axial stresses and volumetric strains were computed finally. For the tests at dry state, a constant high suction, which represents the laboratory conditions, was maintained in the specimen. Flooding was simulated by converting negative pore-water pressures (total suctions) of element nodes into zero.

3.1 Behaviour of rockfill compacted to various initial void ratios

At first, the model was used to predict drained triaxial tests on saturated rockfills with different initial void ratios. The mixture of weak and strong weathered limestone, which was selected as the rockfill material for the Guanyinyan dam (Panzhihua, China), was used in the triaxial tests. The gradation of the tested material was scaled down by making the specimen’s cumulative grain size distributions parallel to the gradation curve of the in situ rockfill material combining with the equivalent substitution method. Figure 5 shows the grain size distribution. The initial void ratio was controlled as 0.24, 0.28 and 0.37 in the sample preparation, respectively. Before shearing, the samples were saturated and isotropically compressed to the confining pressures ranging from 0.2 to 1.2 MPa. The model parameters, which are listed in Table 1, were determined using the procedure introduced in the above section based on the experimental stress-strain relationships. Note that, rockfill specimens for laboratory tests were prepared by compacting the material in thin layers to obtain uniform samples [63]. Thus, a significant pre-consolidation effective mean stress (p0= 0.3 MPa), associated with compaction, was used in the model to simulate the triaxial tests.

Figure 5.

Grain size distribution curves.

Elastic parametersIsotropic compression parametersCritical state parametersDefault parameters
G0 = 125λi-κ = 0.001 MPa−1M0 = 1.95m1 = 0.36
ν = 0.28λ0d= 0.02 MPa−1Mres = 1.55m2 = 0.1
κψ = 0αψ = 0.001 MPa−1αM = 1.2 MPa−1n = 3.5
py = 0.01βψ = 0.01
p0= 0.3 MPaeref0= 0.51
e0 = 0.24, 0.28 and 0.37pref = 0.02 MPa
ζ = 0.105
αe = 0.03

Table 1.

Model parameters for tested limestone.

The comparisons between the numerical results and experimental measurements at different confining pressures are shown in Figure 6. It can be found that proposed model is satisfied to capture the key features of the loose and dense saturated specimens with one unique set of soil parameters. For instance, the denser of the specimen, the larger the deformation modulus and shear strength. Besides, it clearly shows that the rockfill with a smaller initial void ratio (denser state) dilates more remarkably due to the difference between the peak stress ratio and critical stress ratio is much greater, especially at a low confining pressure. With the increase of initial void ratio, the difference between the peak stress ratio and critical stress ratio decreases, which leads to the decrease of volumetric dilatancy degree. It is worth to note that there is an obvious strain-softening phenomenon observed for the dense specimen under the low confining pressures. To our satisfaction, the proposed model can also predict this strain-softening to some extent, although there is a little difference which may be caused by the difficult to determine the parameters corresponding to the critical state.

Figure 6.

Measured and predicted stress-strain relationship (q vs. εa and εv vs. εa) for rockfill with different initial void ratios under confining pressures: (a)σ3 = 0.2 MPa, (b)σ3 = 0.4 MPa, (c)σ3 = 0.8 MPa and (d)σ3 = 1.2 Mpa.

Figure 7 shows the predicted relationship between the potential failure stress ratio, Mf, and stress ratio, η, for three different initial void ratios when σ3 = 0.2 MPa, which can be used to explain the strain-softening phenomenon. From the initial state to the final state, the state η = Mf (i.e., the peak strength state) appears, shown as the intersections of the 450 line and the Mf (η) lines for different initial void ratios. Before the peak strength state is reached, η < Mf, the rockfill material undergoes a hardening process and remains in a stable state. After the peak strength state, ηbecomes slightly larger than Mf, and then the rockfill material undergoes a softening process that is unstable. It can be found that the part of η > Mf becomes more remarkable with the decrease of the initial void ratio, which explains the appearance of the most significant strain-softening process for the densest specimen (e0 = 0.24). Additionally, Figure 8 presents the influence of the confining pressure on the changes of potential failure stress ratio with the changes of stress ratio when e0 = 0.24. It can be found from this figure that Mf (η) line for the larger confining pressure is always in the left side of the 450 line (η < Mf) from the initial state to the final state. Thus, there is no softening phenomenon appearing when there is a large confining pressure applied on the rockfills, and the deviatoric stress only increases up to the strength corresponding to the critical state.

Figure 7.

Relationships between potential failure stress ratio, mf, and stress ratio, η, for three different initial void ratios whenσ3 = 0.2 MPa.

Figure 8.

Relationships between potential failure stress ratio, mf, and stress ratio, η, for different confining pressures whene0 = 0.24.

Noting that, the difference of stress-strain curve between the measured and predicted is relatively larger with the increase of confining pressure, which is shown in Figure 6(c) and (d). Compared with that experimental stress-strain response, the predicted deviatoric stress increases more slowly with the increase of axial strain. From our analysis, this may be caused by the significant particle breakage under the high confining pressures for the rockfills used in the tests, which led to the considerable increase of compressibility modulus before shearing. However, the influence of the particle breakage on the compressibility modulus of rockfill materials has not been considered in the proposed model. In the construction of model, a linear relationship in the e-p plane is assumed, which means the compression index is kept unchanged throughout the isotropic and triaxial compression process. Even so, the quantitative results are in general in good agreement with the experimental measurements. Finally, from above comparisons, it is indicated that the behaviours of the state-dependent dilatancy and strain-softening for rockfill materials can be well captured by the proposed model.

3.2 Behaviour of rockfill with various RH

In order to illustrate the capability of proposed model to predict the stress-strain relationships of unsaturated rockfills, the results of RH-controlled triaxial tests on the crushed cambric slate were adopted as comparison. These tests were conducted by Chavez and Alonso [7], which involved two series of tests. The first series was performed in the fully saturated condition. After compaction, the samples were confined under pressures of 0.1, 0.3, 0.5 and 0.8 MPa and then sheared under the constant confining pressure. In a second series, samples were compacted in a similar manner as saturated ones and maintained at a constant low relative humidity (RH = 36%, equivalent to the total suction, s = 142 MPa). After that, similar stress paths were applied to these dry samples. The model parameters are given in Table 2. Most of parameters were taken from the Chavez and Alonso [7], while the rest were determined according to the experimental stress-strain curves.

Elastic parametersIsotropic compression parametersCritical state parametersDefault parameters
G0 = 80λi-κ = 0.001 MPa−1M0 = 2.1m1 = 0.36
ν = 0.29λ0d= 7.48 × 10−2 MPa−1Mres = 1.45m2 = 0.1
κψ = 0.4 × 10−3αs = 7.86 × 10−3 MPa−1αM = 1.8 MPa−1n = 3.5
py = 0.01 MPaβψ = 0.01
p0= 0.3 MPaeref0= 1.0
e0 = 0.6pref = 0.01 MPa
ζ = 0.1
αe = 0.03

Table 2.

Model parameters for crushed cambric slate.

Figure 9 shows the predicted and measured stress-strain curves in the conditions of fully saturated and dry (RH = 36%). For the saturated sample under a low confining pressure (0.1 MPa), the deviatoric stress increases with the axial strain up to a constant value for axial strains reaching around 10% and then remains constant, which shows the hardening process. However, for the dry sample with the same confining pressure, a peak strength is observed for axial strain reaching the value between 4 and 6%, and a little strain-softening appears after that. On the other hand, when the sample is sheared under the high confining pressures, the increase of deviatoric stress with the axial strain continues at the end of triaxial tests for both saturated and dry samples. Thus, it is difficult to determine the critical state condition in this case. It can be found a good agreement of deviatoric curves (qvs. εa) between modelling and experimental data at both saturated and dry states. Besides, the unloading-reloading process can be well predicted as well with the proposed model. However, some difficulties are observed in reproducing the volumetric response by this model. Computed dilation rates, especially at the low confining pressures, tend to underestimate measured values. Experimental measured dilation rates are significantly larger than that predicted, especially for the unsaturated rockfills under the low confining pressures. There are two reasons accounting for this difference. The first is possibly due to the assumption of critical state conditions in the determination of model parameters. In fact, the attainment of the critical state conditions for a gravel material is probably outside the possibilities of a triaxial test. Very larger strains and a fully broken gravel structure are probably required, but strain localization is also the failure mode in most of tests, which will complicate the interpretation. Thus, the proper observation of critical state conditions would require the larger deformations which were not achieved in the triaxial tests. The second relevant reason is the shape of particles in a gravel. The tested rockfill material has elongated particles with sharp edges, so interlocking and dilation seem very likely to happen even after some significant breakage and deformation, and this also increases the difficulty to modelling the volumetric strains.

Figure 9.

Comparison between model predications and experimental measurements: Triaxial tests on rockfill in the fully saturated and dry states (RH = 36%). (a) σ3 = 0.1 MPa. (b) σ3 = 0.3 MPa. (c) σ3 = 0.5 MPa. (d) σ3 = 0.8 MPa.

3.3 Behaviour of rockfill during shearing and wetting

Finally, two sets of triaxial wetting tests on rockfills considering the suction reduction [32] were adopted as comparison to verify the capability to predict the collapse deformation for the proposed model. The tested materials were schists and greywacke, respectively, which were used as the rockfill materials for the inner and outer shell of Beliche Dam. Tests were performed on a dry material first (presumably equilibrated with the relative humidity prevailing in the laboratory, for total suction equal to 20 MPa).

At some stage during the triaxial tests, specimens were then flooded. The applied vertical strain rate was maintained during the sample inundation and subsequent straining. The stress-suction path actually experienced by the specimens in the stress space (p, q and ψ) is illustrated in Figure 10. The compacted specimen is initially loaded isotropically under very dry condition (point A). The specimen is then loaded vertically at constant suction until the limiting condition (point B) is reached. Thereafter, it was flooded with the constant of vertical strain rate maintained. The volumetric collapse induces a sudden loss of vertical stress for the reduction of suction to zero, and Point C is reached. As the total suction reaching zero, the yield envelope experiences a reduction in size. As the vertical strain rate is maintained, the specimen starts resisting again along the loading path and eventually reaches a new yield condition at point D, and then a final limiting state at point E. The model parameters for schist and greywacke are listed in Table 3, for most of the parameters are referred from Alonso et al. [32], and the rest are determined based on the corresponding experimental results combined with back-analysis method.

Figure 10.

Sketch showing stress-suction path applied in triaxial tests initially dry and later flooded when limiting conditions are reached.

Elastic parametersIsotropic compression parametersCritical state parametersDefault parameters
Inner schists/outer greywacke
G0 = 80/100λi-κ = 0.025/0.01 MPa−1M0 = 1.75/1.8m1 = 0.36/0.4
ν = 0.3/0.3λ0d= 0.028/0.01 MPa−1Mres = 1.2/1.6m2 = 0.16/0.2
κψ = 0/0αs = 7.0/2.0 (×10−3 MPa−1)αM = 0.1/0.5 MPa−1n = 3.5/3.5
py = 0.01/0.01 MPaβψ = 0.03/0.01
p0= 0.3/0.3 MPaeref0= 0.72/0.87
e0 = 0.538/0.538pref = 0.01/0.01 MPa
ζ = 0.22/0.15
αe = 0.042/0.01

Table 3.

Model parameters for inner schists and outer greywacke.

The predicted and measured responses of triaxial wetting tests on inner schists and outer greywacke are represented in Figures 11 and 12. It can be found that an overall good agreement is achieved between the experimental observations and the model predictions. The measured stress-strain curves in the dry condition, the transient loss of deviatoric stress associated with flooding and the final recovery of strength in the saturated condition can be well captured with the proposed model. However, the dilatancy is underestimated, especially for samples under the low confining pressures. It is worth to note that the shear dilatancy response in the shearing process after flooding for rockfills of outer greywacke can be reproduced by the proposed model, although the predicted compression volumetric strain is a little larger than that measured. In general, the proposed model can reasonably predict the collapse deformation of the unsaturated rockfill materials subjected to the flooding, although it also has a little shortage to predict the volumetric strain.

Figure 11.

Triaxial tests on compacted schists from inner shell: (a)qvs.εaand (b)εvvs.εa.

Figure 12.

Triaxial tests on compacted greywacke from outer shell: (a)qvs.εaand (b)εvvs.εa.


4. Conclusions

For modelling the strength and deformation behaviour of rockfill materials over various ranges of densities and relative humidity, a state-dependent elastoplastic constitutive model was developed in this study. This model is formulated within an extended critical-state framework by using two independent stress state variables: total stress and total suction. The compressibility model for rockfill materials proposed by Oldecop and Alonso [1] was adopted in this model. The exponential function was used to describe the decrease of the critical stress ratio with the mean stress, while the logarithmic increases of the critical stress ratio and critical void ratio with the total suction were assumed. By substituting these functions into the relation function of peak friction angle with critical friction angle for granular materials proposed by Biarez and Hicher [63], a unified hardening parameter, which could reflect the state-dependent dilatancy, strain-softening and unsaturated effects, was presented for modelling stress-strain relationships of rockfill materials. In addition, the relevant testing procedures have been demonstrated to calibrate and obtain the required model parameters.

The capabilities of the proposed model were illustrated by modelling three series of triaxial tests performed on different rockfill materials (drained triaxial test on saturated samples compacted to various l void ratio, triaxial test with relative humidity controlled and triaxial wetting test). An agreement was obtained between the experimental and numerical results. For the drained triaxial tests on the saturated samples with different initial void ratios, both stress-strain relationship and volumetric strain can be well predicted, which indicates the capability of the proposed model to consider state-dependent behaviour. Through the comparison with the results of relative humidity-controlled triaxial test, the proposed model can be found to well consider the influence of the relative humidity on the stress-strain relationship and volumetric strain of rockfill materials. In addition, the proposed model can well simulate the feature of stress-relaxing according to the comparison with the results of triaxial wetting test. However, the dilatancy is underestimated with the proposed model, especially for the rockfills at the dry state and under the low confining pressures. In general, the model presents an advantage to predict the state-dependent behaviour and collapse deformation of unsaturated rockfills, such as the strain-softening and dilatancy at the dense state. Despite the progress achieved, there is still a room for improvements to be made, for instance, to consider the influence of particle breakage on the changes of compressibility and critical state at different relative humidities.


  1. 1. Oldecop LA, Alonso EE. A model for rockfill compressibility. Géotechnique. 2001;51(2):127-139
  2. 2. Sowers GF, Williams RC, Wallace TS. Compressibility of broken rock and settlement of rockfills. In: Proceedings of the 6th ICSMFE; Montreal. Vol. 2. 1965. pp. 561-565
  3. 3. Nobari ES, Duncan JM. Effect of reservoir filling on stresses and movements in earth and rockfill dams. Report No. TE-72-1. Department of Civil Engineering, University of California; 1972
  4. 4. Marsal RJ. Mechanical properties of rockfill. In: Hirschfeld RC, Poulos SJ, editors. Embankment Dam Engineering. Casagrande Volume. New York: John Wiley & Sons; 1973. pp. 109-200
  5. 5. Oldecop LA, Alonso E(AP d A). Suction effects on rockfill compressibility. Géotechnique. 2003;53(2):213-234
  6. 6. Oldecop LA, Alonso EE. Testing rockfill under relative humidity control. Geotechnical Testing Journal. 2004;27(3):269-278
  7. 7. Chávez C, Alonso EE. A constitutive model for crushed granular aggregates which includes suction effects. Soils and Foundations. 2003;43(4):215-227
  8. 8. Chávez C, Romero E, Alonso EE. A rockfill triaxial cell with suction control. Geotechnical Testing Journal. 2009;32(3):219-231
  9. 9. Bauer E, Fu ZZ, Liu SH. Hypoplastic constitutive modelling of wetting deformation of weathered rockfill materials. Frontiers of Architecture and Civil Engineering in China. 2010;4(1):78-91
  10. 10. Kohgo Y, Asano I, Hayashida Y. Mechanical properties of unsaturated low quality rockfills. Soils and Foundations. 2007a;47(5):947-959
  11. 11. Marsal RJ. Large scale testing of rockfill materials. Journal of the Soil Mechanics and Foundations Division. 1967;93(2):27-43
  12. 12. Charles JA, Watts KS. The influence of confining pressure on the shear strength of compacted rockfill. Géotechnique. 1980;30(4):353-367
  13. 13. Dine BSE, Dupla JC, Frank R, et al. Mechanical characterization of matrix coarse-grained soils with a large-sized triaxial device. Canadian Geotechnical Journal. 2010;47(4):425-438
  14. 14. Weng MC, Chu BL, Ho YL. Elastoplastic deformation characteristics of gravelly soils. Journal of Geotechnical and Geoenvironmental Engineering. 2013;139(6):947-955
  15. 15. Xiao Y, Liu H, Chen Y, et al. Bounding surface model for rockfill materials dependent on density and pressure under triaxial stress conditions. Journal of Engineering Mechanics. 2014a;140(4):04014002
  16. 16. Alonso EE, Romero EE, Ortega E. Yielding of rockfill in relative humidity-controlled triaxial experiments. Acta Geotechnica. 2016;11(3):455-477
  17. 17. Varadarajan A, Sharma KG, Abbas SM, et al. Constitutive model for rockfill materials and determination of material constants. International Journal of Geomechanics. 2006;6(4):226-237
  18. 18. Honkanadavar NP, Sharma KG. Testing and modeling the behavior of riverbed and blasted quarried rockfill materials. International Journal of Geomechanics. 2014;14(6):04014028
  19. 19. Fu Z, Chen S, Peng C. Modeling cyclic behavior of rockfill materials in a framework of generalized plasticity. International Journal of Geomechanics. 2014;14(2):191-204. DOI: 10.1061/(ASCE)GM.1943-5622.0000302
  20. 20. Tosun H. Discussion of “construction of concrete-faced rockfill dams with weak rocks” by Hao-Feng Xing, Xiao-Nan Gong, Xiao-Guang Zhou, and Hai-Feng Fu. Journal of Geotechnical and Geoenvironmental Engineering. 2006;132(6):778-785. DOI: 10.1061/ASCE1090-02412006132:6778
  21. 21. Naylor DJ, Tong SL, Shahkarami AA. Numerical modelling of saturation shrinkage. In: Proceedings of 3rd International Symposium on Numerical Models in Geomechanics; Niagara Falls. 1989. pp. 636-648
  22. 22. Shen ZJ, Wang JP. Numerical simulation of construction behaviour of clay core dam and its movement due to reservoir impounding. Hydro-Science and Engineering. 1988;4(1):47-64
  23. 23. Li GY, Wang LS, Mi ZK. Research on stress-strain behaviour of soil core rockfill dam. Chinese Journal of Rock Mechanics and Engineering. 2004;23(8):1363-1369
  24. 24. Fang XS. Test study and numerical simulation on wetting deformation of gravel sand [dissertation of masteral degree]. Nanjing: Hohai University; 2005
  25. 25. Wei S. Study on wetting deformation behaviour and numerical model of coarse-grained materials [dissertation of doctoral degree]. Nanjing: Hohai University; 2006
  26. 26. Naylor DJ, Maranha JR, Neves EMD, et al. A back-analysis of Beliche dam. Géotechnique. 2001;51(4):377-381
  27. 27. Fu ZZ, Chen SS, Liu SH. Hypoplastic constitutive modelling of the wetting induced creep of rockfill materials. Science China Technological Sciences. 2012;55(7):2066-2082
  28. 28. Marachi ND, Chan CK, Seed HB, et al. Strength and deformation characteristics of rockfill materials. Report No. TE-69-5. Department of Civil Engineering, University of California; 1969
  29. 29. Terzaghi K. Discussion on settlement of salt springs and Lower Bear River concrete face dams. Transactions of the American Society of Civil Engineers. 1960;123(2):139-148
  30. 30. Alonso EE. A constitutive model for partially saturated soils. Géotechnique. 1990;40(3):405-430
  31. 31. Alonso EE, Olivella S, Pinyol NM. A review of Beliche Dam. Géotechnique. 2005;55(4):267-285
  32. 32. Alonso EE, Olivella S, Pinyol NM, et al. Modelling the response of Lechago earth and rockfill dam. Géotechnique. 2011;61(5):387-407
  33. 33. Kohgo Y, Asano I, Hayashida Y. An elastoplastic model for unsaturated rockfills and its simulations of laboratory tests. Soils and Foundations. 2007;47(5):919-929
  34. 34. Bauer E. Hypoplastic modelling of moisture-sensitive weathered rockfill materials. Acta Geotechnica. 2009;4(4):261-272
  35. 35. Jefferies MG. Nor-sand: A simple critical state model for sand. Géotechnique. 1993;43(1):91-103
  36. 36. Bauer E. Calibration of a comprehensive hypoplastic model for granular materials. Soils and Foundations. 1996;36(1):13-26
  37. 37. Gudehus G. A comprehensive constitutive equation for granular materials. Journal of the Japanese Geotechnical Society. 1996;36(1):1-12
  38. 38. Wu W, Bauer E, Kolymbas D. Hypoplastic constitutive model with critical state for granular materials. Mechanics of Materials. 1996;23(1):45-69
  39. 39. Manzari MT, Dafalias YF. A critical state two-surface plasticity model for sands. Géotechnique. 1997;47(2):255-272
  40. 40. Gajo A, Wood DM. Severn-Trent sand: A kinematic-hardening constitutive model: The q-p formulation. Géotechnique. 1999;49(5):595-614
  41. 41. Li XS, Dafalias YF. Dilatancy for cohesionless soils. Géotechnique. 2000;50(4):449-460
  42. 42. Wang ZL, Dafalias YF, Li XS, et al. State pressure index for modeling sand behavior. Journal of Geotechnical & Geoenvironmental Engineering. 2002;128(6):511-519
  43. 43. Yao YP, Sun DA, Luo T. A critical state model for sands dependent on stress and density. 2004;28(4):323-337
  44. 44. Loukidis D, Salgado R. Modeling sand response using two-surface plasticity. Computers and Geotechnics. 2009;36(1–2):166-186
  45. 45. Chakraborty T, Salgado R, Loukidis D. A two-surface plasticity model for clay. Computers and Geotechnics. 2013;49(APR.):170-190
  46. 46. Xiao Y, Liu H, Chen Y, et al. Testing and modeling of the state-dependent behaviors of rockfill material. Computers and Geotechnics. 2014;61(1):153-165
  47. 47. Sun Y, Gao Y, Zhu Q. Fractional order plasticity modelling of state-dependent behaviour of granular soils without using plastic potential. International Journal of Plasticity. 2017;(3):53-69
  48. 48. Liu SH, Wang LJ, Wang ZJ, Bauer E. Numerical stress-deformation analysis of cut-off wall in clay-core rockfill dam on thick overburden. Water Science and Engineering. 2016;9(3):219-226
  49. 49. Wang LJ, Liu SH, Fu ZZ, Li Z. Coupled hydro-mechanical analysis of slope under rainfall using modified elasto-plastic model for unsaturated soils. Journal of Central South University. 2015;22(5):1892-1900
  50. 50. Wang L, Liu S, Liao J, et al. Field load tests and modelling of soft foundation reinforced by soilbags. Geosynthetics International. 2019:1-36
  51. 51. Sheng D, Sloan SW, Gens A, et al. Finite element formulation and algorithms for unsaturated soils. Part 1: Theory. International Journal for Numerical and Analytical Methods in Geomechanics. 2003a;27(9):745-765
  52. 52. Coussy O. Mechanics of Porous Continua. Chichester, UK: Wiley; 1995
  53. 53. McDowell GR, Bolton MD. On the micromechanics of crushable aggregates. Géotechnique. 1998;48(5):667-679
  54. 54. Yao YP, Hou W, Zhou AN. UH model: Three-dimensional unified hardening model for overconsolidated clays. Géotechnique. 2009;59(5):451-469
  55. 55. Indraratna B. Large-scale triaxial testing of greywacke rockfill. Géotechnique. 1993;43(1):37-51
  56. 56. Frossard E, Dano C, Hu W, et al. Rockfill shear strength evaluation: A rational method based on size effects. Géotechnique. 2012;62(5):415-427
  57. 57. Xiao Y, Liu HL, Zhu JG, et al. Modeling and behaviours of rockfill materials in three-dimensional stress space. Science China Technological Sciences. 2012;55(10):211-226
  58. 58. Yao YP, Sun DA, Matsuoka H. A unified constitutive model for both clay and sand with hardening parameter independent on stress path. Computers and Geotechnics. 2008;35(2):210-222
  59. 59. Yao YP, Zhou AN. Non-isothermal unified hardening model: A thermo-elasto-plastic model for clays. Géotechnique. 2013;63(15):1328-1345
  60. 60. Biarez J, Hicher P-Y. Elementary Mechanics of Soil Behaviour. Rotterdam, Netherlands: A.A. Balkema; 1994
  61. 61. Hicher PY, Chang CS. A microstructural elastoplastic model for unsaturated granular materials. International Journal of Solids and Structures. 2007;44(7–8):2304-2323
  62. 62. Sheng D, Smith DW, Sloan SW, et al. Finite element formulation and algorithms for unsaturated soils. Part II: Verification and application. International Journal for Numerical and Analytical Methods in Geomechanics. 2003b;27(9):767-790
  63. 63. Naylor DJ, Maranha DNE, Veiga Pinto AA, et al. Prediction of construction performance of Beliche dam. Géotechnique. 1986;36(3):359-376

Written By

Liujiang Wang and Zhongzhi Fu

Submitted: September 26th, 2019 Reviewed: May 21st, 2020 Published: September 18th, 2020