Modelling the poroplastic damageable behaviour of earthen materials

This paper presents a new constitutive model on the poroplastic behaviour of earthen materials accounting for stiffness degradation, using the approach of continuum damage mechanics. The poroplastic behaviour is modelled based on the bounding surface plasticity (BSP) theory and the concept of effective stress while isotropic damage is modelled using a scalar variable. Plastic flow and damage evolution occur simultaneously in a coupled process which take into account the impact of suction. The model was successfully validated against results of triaxial compression tests performed at different relative humidities and confining pressures. Despite the relatively small number of material parameters, this model can reproduce the essential features of earthen materials behaviour observed experimentally: suction-induced hardening and stiffening, post-peak softening, as well as the progressive transition from contractant to dilatant volumetric behaviour. Use of the BSP theory allows to reproduce a smooth stress–strain relation as experimentally observed, instead of an abrupt change upon plastic yielding predicted by classic elastoplastic models. Furthermore, the present model also furnishes a quantitative description on the degradation of elastic properties hitherto not accounted for, thanks to the additional scalar damage variable.


Introduction
Unfired crude (untreated) earth is a hygroscopic porous material that contains generally a small quantity of active clay minerals. Currently, considering the context of sustainable construction, reduction of energy consumption and greenhouse gas emission, it is regarded as a promising non-industrial construction material. Several construction technics using earthen materials have been invented, leading to monolithic walls, like rammed earth and cob, or to masonry units, like adobe or compressed earth blocks [1]. After constructions, earthen buildings being always subject to atmospheric conditions possess relatively low water contents and high suction which impact significantly on their mechanical behaviour.
In order to study the behaviour of unsaturated earthen materials, a series of large-scale laboratory tests were carried out by several researchers [2][3][4][5]. The main conclusion drawn from all these experiments is that the mechanical properties of earthen materials, whatever their origin and their implementation technic, are strongly influenced by their water content: even a small increase of the latter leads to a significant reduction of both the failure strength and Young's modulus E, while the irreversible plastic strain is increased. Apart from plastic deformations and suction effects, a progressive but non-negligible degradation of Young's modulus E with increasing stress level were observed in recent study [6,7]. This stiffness reduction is conjectured to be a consequence of damage induced by initiation and development of micro-cracks. The above experimental observations motivate the development of a new constitutive model which can simulate the poroplastic behaviour of earthen materials while accounting for damage effects.
To begin with, the main drawback of classic elastoplastic models is their erroneous prediction of a sharp transition at yield from elastic to elasticplastic behaviour (jump of apparent stiffness), as well as their difficulty to reproduce realistic volumetric behavior (typically a smooth transition from contractancy to dilatancy). Models developed using the Bounding Surface Plasticity (BSP) theory can overcome these limitations. This theory, firstly proposed by Mroz [8] and Dafalias et al. [9,10] to simulate metal behaviors, was successively applied by Bardet [11] on saturated sands and later on followed by Gajo and Muirwood [12] as well as by Yu and Khong [13]. BSP theory was later extended to partially saturated materials in order to describe suction effects [14][15][16][17][18][19]. In addition, BSP theory is also applied on geomaterials other than soils due to its versatility, such as rockfills and clay-fouled ballast [20,21].
In parallel, there have been extensive studies on the development of coupled elastoplastic damage models [22][23][24][25][26]. However, most of them are devoted to modeling dry materials without considering coupled hydromechanical interactions, not to mention the capillarity effects due to partial saturation. In recent years, several researchers have attempted to couple damage with hydromechanical interactions. For instance, a coupled poro-elastoplastic damage model is proposed by Shao et al. [27] and extended to partially saturated conditions in order to study coupled hydromechanical damageable behaviours in dryingwetting processes. Zhang et al. [28] developed a unified elastic-viscoplastic damage model to describe long-term hydromechanical behavior of argillites under unsaturated condition. Bui et al. presented a new constitutive model accounting simultaneously for the impact of damage on hydraulic and mechanical properties of unsaturated poroplastic geomaterials, by means of the thermodynamic framework for partially saturated media [29]. The present work follows a similar path of this last reference.
After a summary of experimental investigations showing the main trends on the mechanical behaviour of earthen materials, a new constitutive model is presented for unsaturated earthen materials. It is based on the BSP theory and the effective stress concept. It adopts a non-associative flow rule and a simple radial mapping rule, and takes into account an isotropic damage variable. Plastic straining and damage evolution occur as two coupled processes. The influence of suction on plastic straining and damage evolution is taken into account. Finally, the performance of this model is investigated by comparing the numerical simulation with a series of triaxial tests results on compacted earth samples at different hydraulic conditions and confining pressures.

Summary of experimental investigations
The raw earth materials used in this study came from Limonest, in the south-eastern part of France and will here-after be refered to as Lim. It possesses a natural clay content of 35%. To study the impact of clay content on the hydromechanical behaviour, sand is added to obtain two other materials, refered to as Mix1 and Mix2 with clay contents respectively of 17 and 26%. The optimum condition of the Lim earth from proctor test was used to prepare specimens. Besides, in order to make the comparison of different soils meaningful and easily interpretable, the optimum condition (12.5%, 1.95 g/cm 3 ) is also applied to Mix1 and Mix2. Cylindrical samples (3,5 cm diameter and 7 cm high) were fabricated using a home-made mold which allows a double compaction (compression applied simultaneously at both ends of the sample) of the raw earth at controlled displacement rate. After fabrication, samples were cured at three different relative humidity conditions (23, 75, and 97%) at 23°C until reaching equilibrium. The total suctions, denoted by s t , corresponding to the three values of relative humidity (RH) can be calculated by Kelvin's law: where q w and M w are the density and the molar mass of water, R the perfect gas constant and T the absolute temperature. Since distilled water was used in the tests, the solute concentration is expected to be very low; hence, the osmotic suction should be negligible. The total suction in the above relation is therefore reduced to matric suction s. In the rest of this paper, the term 'suction' will therefore refer to 'matric suction'. The matric suctions corresponding to the three values of RH (23, 75, and 97%) imposed at 23°C are 201, 39 and 4 MPa respectively. A series of drained triaxial compression tests were performed at three confinement pressures (0, 100 and 600 kPa). Non-contact sensors were installed to measure the axial and radial strains. Two types of mechanical loadings were applied at constant relative humidity conditions. More details on the experimental device can be found in the paper [7]. The first loading type is classic in triaxial tests: axial load is applied at a constant displacement rate until failure. After measuring the maximum deviator stress at failure (noted f c ) for one particular specimen using this type-1 loading, the second loading type consisting of loading-unloading cycles at prescribed stress levels, respectively 20, 40, 60 and finally 80% of f c was conducted on two other specimens prepared under identical conditions. The second loading-unloading cycle was used to estimate the evolution of Young's modulus E with the level of axial stress applied, assuming unloading is reversible. Young's modulus E is therefore computed from the following equation: where Dr cycle xx ; De cycle xx are the differences respectively of the axial stress and axial strain between the maximal and minimal of a load cycle. The loadings were applied by imposing a constant axial displacement rate of 0.002 mm/s (corresponding to 0.003%/s) in loading and unloading. Figure 1 illustrates the typical stress-strain behaviour (Mix1 at 100 kPa confining pressure) extracted from the results of a series of triaxial compression tests at 3 different relative humidities. Firstly, all tests exhibit a smooth stress-strain behaviour from the beginning to the end, with non linearity appearing practically from the beginning. Secondly, large residual strains are systematically observed after unloading for all type-II tests. These observations indicate the appearance of plastic (irreversible) strains even at low stress levels. In addition, a post-peak softening behaviour is captured for the case of low relative humidity (e.g. RH = 23%). Last but not the least, the typical stress-strain behaviour of earthen materials is sensitive to relative humidity conditions, including the amplitudes of residual strain and f c . All the remarks mentioned above are also observed for Lim and Mix2. Figure 2a shows values of Young's modulus E measured during the unloading-reloading cycle at 20% f c and the evolution with relative humidity. Young's modulus decreases drastically with the rise of relative humidity. Results presented in Fig. 2b on the evolution of Young's modulus E versus stress level indicate that there is a progressive and significant degradation of the material stiffness when the stress increases. This degradation of stiffness is faster at low stress and slower at high stress. Consistently with the stiffness degradation, macro-cracks oriented in the direction of major principal stress appeared during the compression test.
To summarize, results of our experimental investigation indicate that the basic mechanical behaviour of unsaturated earthen materials can be characterized by irreversible plastic straining and stiffness degradation. The first behaviour is commonly modelled using Axial strain ε xx (mm/mm) classic elastoplasticity theory, despite its unrealistically prediction of a stiffness-discontinuity on the stress-strain curve at yield. To improve this, it is proposed to adopt the approach of Bounding Surface Plasticity (BSP). The stiffness degradation due to microcracks has not been considered in the literaure on earthen materials constructions. To address this issue, it is proposed to follow the approach of continuum damage mechanics.

A poro-elastoplastoplastic damageable constitutive model
As mentionned in the previous section, the construction of the new constitutive model requires to address two mechanisms: a poroplastic mechanism capable of simulating the suction-dependent plastic straining and a damage mechanism to capture the stiffness degradation. The poroplastic mechanism in this new model follows the same approach as the CASMNS model developed by Lai et al. [19] for unsaturated soils. The formulation is based on the concepts of effective stress and bounding surface plasticity theory. To accurately simulate volume changes (e.g. transition from compressibility to dilatancy), a state-dependent non-associative flow rule is used, by adopting the plastic potential of Yu [30]. The starting point is the definition of an effective stress to account for effects of partial saturation.

Partial saturation and effective stress
There are two main classical modelling approaches on partially saturated soils. All of them need two independent stress variables. The classic BBM model is based on net stress and suction whereas a few others use an effective stress and suction [14,15,31]. Some recently appeared models proposed to use the degree of saturation instead of suction as a generalised stress variable [32][33][34]. This new approach seems to be promising, but the feedback is still limited. In this paper, we follow the approach based on an effective stress and suction. At partially saturated states, the effective stress tensor r 0 ij can be expressed as [35]: where r ij is the total stress tensor, u a is the pore air pressure, s is the difference between pore air and pore water pressure (u a À u w ), referred as suction, d ij is the second order identity tensor, and v is defined as the effective stress parameter and generally is state dependent. In all the following, a prime (') above a stress variable will denote the effective stress counterpart. In many practical cases, air pressure can be identified as the atmospheric pressure, taken as the reference datum for stresses and fluid pressures, hence, Eq. (3) can be simplified to: In this paper, only the configuration of conventional triaxial tests is of concern, exhibiting cylindrical symmetry, hence axial and radial directions are always the directions of principal stresses and strains. Let r 1 denotes the total axial stress, r 2 ¼ r 3 the total radial stress, e 1 the axial strain and e 2 ¼ e 3 the radial strain. This cylindrical symmetry allows to work with bivectors instead of 3-vectors; however, the various transformation and stiffness matrices lose their symmetry. The classical notation p; q ð Þ is used to denote the mean and equivalent deviatoric stresses. It is evident that: hence, we only need to refer to q in future. Several definitions for v have been proposed in the existing literature: some researchers identified the effective stress parameter with degree of saturation S r [5,36]; while others argued that v should be a function of suction [37]. In this study, the functional form (9) below firstly proposed by Khalili and Khabbaz [38] has been adopted: where a is a material constant and s e represents the airentry suction which can be determined from water retention measurements. It is a classic result that when the failure stress states are plotted in the stress space p 0 ; q ð Þ, with an adequately chosen function v to be applied in the translation (4), they lie approximately on a unique line. Khalili and Khabbaz showed that this result is obtained for a large class of unsaturated soils for the functional form Eq. (6) with a ¼ 0:55 [38]. Based on our experimental results on compacted earth samples, the particular value a ¼ 0:85 appears to lead to the same result for the three different earths studied, as illustrated in Fig. 3 for the case of Lim. This value of a was adopted in the subsequent numerical studies.
For future manipulations, it is useful to introduce the following concise vectorial notations to denote axial/radial (i.e. principal) components and mean/ deviatoric components: where e v is the volumetric strain and e q is the deviatoric strain.
The mean/deviatoric effective stresses are linked to their axial/radial counterparts by the following relations (which evidently also holds for total stresses): The following equation holds for the conjugate strain components: Within the framework of elastoplasticity, strains are supposed to be decomposable as the sum of an elastic (superscripts e) and a plastic component (superscripts p): For the damage mechanism, a damage variable D is introduced as an internal variable to simulate the effects of microcracks. For the sake of simplicity, only isotropic damage is considered here and D 2 0; 1 ð Þis a scalar variable.

General concept of bounding surface plasticity
The concept of bounding surface plasticity (BSP), also known as two-surface plasticity, was first originated by Mroz [8] and Dafalias et al. [9,10], and applied by Bardet [11] and Yu et al. [13] to soils. Russell and Khalili [39] appear to be the first who applied the BSP theory to model unsaturated soils. Similar to the bounding surface models, sub-loading surface models are also widely used for modelling complicated plasticity behaviour for soils. Based on Hashiguchi's initial framework [40], Yao et al. proposed a simple but very useful unified hardening equation [41] and corresponding models to consider complicated soil behaviour like hardening behaviour related to thermal effects [42] and hardening behaviour related to unsaturation [43].
The major advantage of bounding surface plasticity (BSP) theory is the ability to reproduce a progressive appearance of irreversible strains hence a smooth stress-strain behaviour at increasing loads, thereby avoiding an abrupt change predicted by classical elastoplastic models when yielding occurs. For this purpose, two surfaces are introduced generally: a Loading Surface (LS) f ¼ 0 containing the current stress point r 0 and a bounding surface (BS)F ¼ 0. As loading increases so that the current stress point advances in the direction of the BS, inducing additional plastic strains, the BS itself evolves through a hardening mechanism that depends on the accumulated plastic strains. The central idea of BSP is that the plastic strain rate, which governs the rate of hardening hence the movement of the BS, depends on the distance separating the current and image stress point and the BS. The plastic strain rate is designed to be small at large distances from the BS and increases when approaching the latter. In order to define this ''distance'', it is necessary to define an image stress point r 0 on the BS. The ''distance'' in question will then be the one separating r 0 and r 0 . Among those mapping rules (giving r 0 as a function of r 0 ), one frequently used is called ''radial mapping'' introduced by Dafalias due to its simplicity [9]. It consists of extrapolating the position vector linking the origin to the current stress until it intersects the image stress on the Bounding Surface as illustrated in Fig. 4.

Elastic mechanism accounting for damage
It is assumed that the loading surface always lies inside the bounding surface and that there is no ''elastic domain'' such that strain increments contain generally an elastic and a plastic component. To simplify, an isotropic elastic mechanism is adopted, commonly used in other constitutive models [14][15][16]: C e is the elastic compliance matrix which is usually expressed in terms of the bulk and shear moduli K and G. To account for the effects of damage on the stressdependent elastic properties, the following expressions are adopted based on the previous results of [23,27,29]: where D is the isotropic damage variable, m the Poisson ratio, j a material constant, and e 0 the initial void ratio.

Bounding surface, loading surface and hardening mechanism
In this model, a whale-head shaped bounding surface (BS) is adopted. It writes: In the above equation, M represents the slope of critical state line (CSL) in p 0 À q plane, while n and r are two other model constants. The first constant n specifies the shape of the bounding surface and the second constant r is a space ratio used to control the intersection point of the CSL with the BS. p p is a variable analogue to the preconsolidation pressure in the classic Cam Clay model that determines the current size and position of the bounding surface. For simplicity, the loading surface (LS) is defined simply by homotheticity relative to the BS, as illustrated in Fig. 4: After some simple manipulations, we obtain: where b is the scale factor between the BS and the LS, b 2 0; 1 ð Þ, and g is the stress ratio. An isotropic hardening mechanism based on the Cam Clay model is applied here. In the presence of suction and damage effects, the preconsolidation pressure, noted p p , is assumed to take the following form, inspired from previous work [27,31]: In the above relation, 1 À D ð Þrepresents a softening effect due to damage, whereas the other factor 1 þ k 1 Â l s ð Þ ð Þ accounts for suction-hardening, k 1 being a model constant and the function l s ð Þ is defined in Eq. (19) here-below. At positive suction, the preconsolidation pressure at full saturation p 0 is supposed to vary with the plastic volumetric strain according to the following relation: The suction-dependent coefficient of virgin (i.e. elastoplastic) compressibility k s ð Þ is assumed to take the following form: here k 0 is the limiting value of k s ð Þ at full saturation, and j is the elastic coefficient of compressibility (applicable for unloading), assumed suction-independent. The function l s ð Þ plays an important role here and must be chosen properly, ensuring continuity at full saturation and the capacity to simulate wetting collapse. The following form for l s ð Þ is adopted in this work: where the suction-dependent effective stress parameter v is defined in Eq. (6).
3.5 Plastic potential, non-associative flow rule and plastic modulus The dilatancy, defined as the ratio between the plastic volumetric strain increment and the equivalent plastic deviatoric strain increment, has an important bearing on the soil behaviour. This state-dependent parameter can be incorporated in the definition of the plastic potential. In the present model, the following stressdependent relation proposed by Yu is adopted [30]: In the above equation, m is another model constant which governs the plastic potential, whereas constants M and n as well as the stress ratio g have already been defined. This equation implies the dilatancy depends only on the stress ratio and becomes zero at the critical stress ratio g ¼ M. A plastic flow rule which satisfies the above dilatance ratio can be written as: where m is a unit vector with components obtained from Eq. (20): and dk is the plastic multiplier to be determined from the consistency condition. In this study, plastic flow is coupled with damage evolution under a general loading condition. Therefore, the plastic strain rate and damage evolution rate have to be determined simultaneously. The consistency condition on the bounding surface requires: The reference stress r 0 ¼ 1 MPa has been added to make the dimension of different terms homogeneous. The plastic volumetric strain rate have already been defined in Eq. (21) and the damage evolution rate dD will be defined in Eq. (34). The above equation can be recast into a slightly more compact form by normalising the pseudo-gradient of the BS. The following notation is introduced: wherer is a pseudo-gradient operator and jj:jj represents the classic Euclidean norm of a vector. Similarly, the pseudo-gradient and the associated unit vector for the LS can be defined: The classic assumption that stress increments at the current stress point and at the image stress point give the same plastic multiplier allows to write [9]: dk ¼ n p dp 0 þ n q dq þ n D r 0 dD þ n s dS H where H is plastic modulus on the loading surface, and H b the plastic modulus on the bounding surface. When the current stress point is on the BS, hence coinciding with its image stress point, the combination of the Eqs. (21), (23), (24) and (26) leads to the following expression of the plastic modulus: At this point, we will introduce the simplifying assumption that suction will remain constant, ds ¼ 0, which corresponds to the experimental conditions. Combining Eqs. (13), (16), (17) and (22) with (27), the hardening modulus at the image point can be further developed to: In the more general case when current stress point lies strictly inside the bounding surface, the plastic modulus H is the sum of H b and another term, denoted by H d . This additive term depends in some way on the distance d between r 0 and r 0 . In this present model, H d is postulated to be a function of the scale factor b defined in Eq. (15), inspired from the previous work of Bardet and Morvan et al. [11,15]. The final form of the plastic modulus H is written as follows: in which h is a material constant. Note that the above expression of the plastic modulus satisfies the following two conditions: corresponding respectively to the two limiting cases when the current stress lies on the BS and when it is far from the BS.

Damage mechanism
Similar to the case of plasticity, damage kinetics is controlled by a pseudo-potential of dissipation. It is often described by a scalar-valued damage criterion function which depends on the damage variable D and/ or the conjugate damage force Y D . Although the latter can be derived directly from the framework of thermodynamics, this usually leads to overly complex expressions. An alternative simpler family of damage evolution laws, which has been successfully applied to concrete and rocks under compressive deviatoric stresses is based on the tensile part of principal strains.
Following the approach of Mazars [23], the tensile part of the principal strain, noted e þ i , is defined by: Note that the above definition is consistent with the sign convention of positive compressive strains. Moreover, in conventional triaxial compression tests, the axial compressive strain e 1 is positive while the radial tensile strain e 3 is negative.
The damage force Y D is assumed to depend on the tensile parts of principal strains defined above and a simple damage criterion is proposed in the following form based on [26,28,29]: The function r D; s ð Þ represents the damage threshold which depends on the current values of damage and suction. More specifically, r 0 defines the initial damage threshold and r 1 s ð ÞD accounts for the drift due to damage evolution and suction, with the parameters k 3 controlling the suction-dependent sensitivity. The damage evolution rate is determined through the classic ''flow'' rule: in which the damage multiplier dk D should be determined from the consistency condition of damage criterion, fulfilling the Kuhn-Tucker condition: Specifically, in the case of conventional triaxial tests at constant suction ds ¼ 0, the damage evolution rate can be deduced as the following expression: 3.7 Elastoplastic compliance matrix accounting for damage For the particular case of constant suction (the only case considered in this paper) ds ¼ 0, the plastic strain rate in Eq. (21), on account of Eq. (26), can be rewritten in matrix form: where C p is the plastic compliance matrix. Summation of Eqs. (11) and (37) leads us to: Taking into account the damage evolution rate dD deduced from Eq. (36) and combining with Eqs. (8) and (9), the following relationship between principle strain and stress increments is obtained after some simple manipulations: with: where C epD represents the elastoplastic compliance matrix accounting for damage. The above equations define completely the relation between strain and stress increments.
4 Model validation against experimental data and discussion

Determination of model parameters
The proposed bounding surface plastic damage model requires 14 parameters for its complete definition. They characterize particular aspects of material behaviour: (1) m and j to describe the elastic behaviour.
(2) Seven constants on the plastic behaviour: k 0 and C define the position of critical state line (CSL) in the e À ln p 0 plane while M defines the slope of CSL in p 0 À q plane; n and r specify the shape of bounding surface; h controls the hardening modulus; and m controls the stress-dilatancy. (3) Another four model parameters are needed to account for suction effects: s e and a define the effective stress used in this model, k 1 and k 2 account for suction effects on hardening and volumetric deformability respectively. (4) Finally, the constant k 3 controls damage evolution. The above model parameters can be determined from a few classic tests in the laboratory. m, j and k 0 can be obtained from classic oedometric and triaxial consolidation tests at full saturation under drained conditions: j is based on the slope of elastic unloading-reloading line on e À ln p 0 plane, m is the Poisson's ratio (can be determined from axial vs radial strain rates during elastic unloading) and k 0 is detemined from the loading part of the isotropic triaxial consolidation tests expressed in e À ln p 0 plane. C and M associated with the definition of critical state line can be deduced using triaxial tests at different controlled suctions. The parameter k 2 defining the reduction of volumetric deformability k s ð Þ with suction can be assessed by carrying out virgin oedometric compression tests at different controlled suctions. s e is identified as air-entry suction determined from the classical water retention curve and a unique value a ¼ 0:85 was adopted for the three different soils as mentioned in the Sect. 3.1. The remaining parameters: h; m; n; r; k 1 and k 3 have to be determined by trial and error using experimental results.

Model implementation
To implement the constitutive model proposed in this paper, a simple computer program was developed under Matlab environment. The performance of this model was investigated by comparing the numerical simulation with a series of triaxial tests results for earthen materials (Lim and Mix1), at different hydraulic conditions (23,75, and 97% relative humidity) and two confining pressure conditions (100 and 600 kPa). The parameters used for the study are presented in Table 1.  experiments) and transitions from pre-peak hardening to post-peak softening were also well simulated. As shown in Fig. 5b, the volumetric strain versus axial strain e v -e 1 as predicted by this model is also in good agreement with the experimental data. This confirms that the model can indeed reproduce the complex volumetric behaviour: transition from contraction to dilation state. The model also captures the suppressed dilation induced by the increasing relative humidity (suction decreasing).

Stress-strain behaviour and volumetric evolution
At lower hydraulic state RH = 23%, some deviations are observed both in stress strain curve and volumetric evolution, especially for the post-peaking stage. One possible explanations is that some form of localization may appear, which cannot be simulated by a purely continuous modelling approach. Note that this kind of deviation is also faced with in previous modelling work [27,29]. Figure 6a, b shows the results of triaxial compression tests and model simulations of Lim at 600 kPa confining pressure. In q-e 1 plane, the results predicted by the model are quitredicted by this model is alson e ve 1 plot, it is seen the dilatancy reproduced by the model  is overestimated compared with the test results especially for the RH = 23% case. To improve the model's performance, a suction-dependency plastic potential is considered to use in the future work.
By making a comparison between Figs. 5 and 6, it is worth noting that the model also captures the increased deviator stress and overall reduced dilation with an increase in confining pressure. As confining pressure increases, particle crushing begins to dominate plastic deformation, resulting in a reducing tendency for a peak in the shear resistance and for a volumetric expansion.
Triaxial compression tests on Mix1 at three different hydraulic states (RH = 23%, RH = 75% and RH = 97%) and at 100 kPa confining pressure were conducted in the laboratory. A good agreement is obtained between numerical simulation and experimental data as illustrated in Fig. 7. Figure 8 presents evolutions of Young's modulus versus radial strain during triaxial tests on Mix1 at confining pressure of 100 kPa and two relative humidity conditions (RH = 75% and RH = 97%). There is an acceptable agreement between numerical predictions and experimental data, indicating the developed model qualitatively describes the general trend observed in test data: elastic modulus decreases with the damage induced by the radial strain. Considering there is only an additive parameter k 3 to describe the damage evolution rate, the simulation results is acceptable.

Degradation of elastic property
In addition, it is worth noticing the degradation rate of Young's modulus at RH = 75% (higher suction, lower water content) is higher than that at RH = 97% (lower suction, higher water content), which is conform to the physical intuition that higher water content should result in higher ductility.
Finally, in order to illustrate the effects of damage on the stress-strain behaviour, numerical simulations with different values of k 3 are presented for Mix1 at 75% RH and 100 kPa confining pressure (Fig. 9). With decreasing of k 3 , the damage evolution rate becomes higher. From the simulation results, it could be observed that the response of compacted earthen materials becomes more and more brittle when the effect of damage on plastic flow is stronger.

Conclusions
A new constitutive model for unsaturated earthen materials is proposed accounting simultaneously for the impact of suction and confining pressure on mechanical properties. This model is formulated using the concept of effective stress and bounding surface plasticity (BSP) theory under a critical state framework. It adopts an whale-head shaped loading and bounding surfaces, a simple radial mapping rule, a non-associative flow rule which generally gives a better description on volumetric behaviour, and suction-dependent hardening controlled by plastic volumetric strains. At the present stage of development, only isotropic damage is considered. Damage Axial strain ε 1 Volumetric strain ε v RH=23% RH=75% RH=97% RH=23%(Experiment data) RH=75%(Experiment data) RH=97%(Experiment data) evolution rate is assumed to be driven by tensile strain and restrained by suction. Plastic flow and damage evolution occur as coupled processes.
The model has been used to simulate the triaxial compression tests subjected to different relative humidities and confining pressures. Good agreement is generally obtained between model predictions and experimental data. Of fundamental importance in practice, this model only requires 14 independent soil parameters for its definition and makes it ideally suited for engineering applications. Despite the minimal number of parameters, the model is able to reproduce the essential trends in the behaviour of partially saturated earthen materials: suction-induced hardening and stiffening, smooth stress-strain behaviour, post-peak softening ae well as contractancydilatancy transition. Furthermore, the present model qualitatively describes a degradation of elastic properties observed in experimental data. The model in its current form is intended for monotonic loading, an extension of this model to cyclic loading will be the subject of future work.