Solenn Le Pense
email: lepense@hw.ac.uk

Chloé Arson

Ahmad Pouya

A fully coupled damage-plasticity model for unsaturated geomaterials accounting for the ductile-brittle transition in drying clayey soils

Keywords: Constitutive laws, Damage, Plasticity, unsaturated geomaterials, Effective stress

This paper presents a hydro-mechanical constitutive model for clayey soils accounting for damageplasticity couplings. Specific features of unsaturated clays such as confining pressure and suction effects on elastic domain and plastic strains are accounted for. A double effective stress incorporating both the effect of suction and damage is defined based on thermodynamical considerations, which results in a unique stress variable being thermodynamically conjugated to elastic strain.

Coupling between damage and plasticity phenomena is achieved by following the principle of strain equivalence and incorporating the double effective stress into plasticity equations. Two distinct criteria are defined for damage and plasticity, which can be activated either independently or simultaneously. Their formulation in terms of effective stress and suction allows them to evolve in the total stress space with suction and damage changes. This leads to a direct coupling between damage and plasticity and allows the model to capture the ductile/brittle behaviour transition occurring when clays are drying. Model predictions are compared with experimental data on Boom Clay, and the flexibility of the model is illustrated by presenting results of simulations in which either damage or plasticity dominates the coupled behaviour.

Introduction

Pressing needs for sustainable structures and safe geological repositories require the development of reliable models to predict the behaviour of natural geomaterials (e.g., soils, rocks) and engineered materials (e.g. compacted backfill materials, cement-based materials, ceramics...).

One of the recurring modelling challenges is the prediction of deformation, stiffness and strength of porous media with a solid matrix containing clay minerals.

Experimental evidence show that clayey soils can exhibit either a brittle or a ductile behaviour [START_REF] Dehandschutter | Brittle fractures and ductile shear bands in argillaceous sediments: inferences from Oligocene Boom Clay (Belgium)[END_REF]. The transition between both behaviours depends on multiple factors including moisture content [START_REF] Al-Shayea | The combined effect of clay and moisture content on the behavior of remolded unsaturated soils[END_REF]. Under deviatoric loading, clayey soils can undergo large permanent strains. Their properties, such as stiffness , strength, or permeability, are also known to be subject to changes after being submitted to hydric or mechanical solicitations. In clayey soils, these changes are related to several physical phenomena, such as the deterioration of cemented bonds, hence the destructuration of the material, or the change in water content, resulting in a decrease of suction-induced bonding. Sophisticated plasticity models proposed for clayey soils allow capturing suction hardening and wetting collapse. However these models are not suited for stiffer or bonded materials, which can undergo both plastic deformation and stiffness degradation. A phenomenological variable d, called "damage", can be defined at the continuum scale to quantify the energy dissipated by stiffness degradation. Note that d represents the effects of multiple microscopic dissipative processes that lead to the loss of adhesion between material surfaces, such as micro-crack propagation or debonding due to an increase of water saturation. Coupling damage and plasticity in a thermodynamically consistent framework raises many issues when one wants to ensure thermodynamical consistency while keeping the model simple enough to allow for easy calibration and incorporation into a numerical code. Models coupling damage and plasticity are often material and loading path specific, and difficult to generalise to a broader category of problems related to the coupled effects of mechanical stress and suction in unsaturated clay-bearing porous media. One of the fun-damental issues that needs to be addressed is the choice of thermodynamic variables, in particular the stress variable involved in the yield and damage criteria.

State-of-the art models are often designed to fit experimental data for specific materials subjected to specific stress-paths. By contrast, the modelling approach presented in this paper is aimed to predict the transition between brittle and ductile deformation regimes for different fabrics, clay contents and hydro-mechanical stress paths. Model calibration and numerical implementation are facilitated by the low number of constitutive parameters employed in the formulation (14 parameters in total, 8 for the mechanical part of the model and 6 for the hydraulic part). The proposed framework is flexible so that each component can be refined if one wants to adapt it to a specific material. The described model is designed in order to be adaptable to a wide range of geomaterials, ranging from stiff clayey soils to mixtures of clay and sand. To illustrate the versatility of the framework, the model was calibrated against experimental data obtained for a variety of geomaterials including Boom clay and mixtures of clay soil and sand.

The work presented in this paper provides a general method to couple damage and plasticity in porous materials that have a clay-bearing solid matrix. Clay minerals are expected to play a critical role in the deformation and retention properties of the damaged medium. Section 2 reviews the main modelling strategies available to date to model hydro-mechanical plasticity and damage in unsaturated porous media, and introduces the main concepts and states variables used to account for the effect of suction and damage. Then, the concept of double effective stress is introduced in section 3, and its coupling with damage and plasticity is developed. In Section 4, the behaviour of the mechanical model is analysed, as well as its limits and its sensitivity to the main parameters. Section 5 presents the simulation of different laboratory tests performed on unsaturated clay-bearing geomaterials. The comparison between models predictions and experimental data reported in the literature is used as a basis to assess the performance of the model.

The sign convention used is the one of soil mechanics. Compressive stresses and strains take therefore positive values.

Damage in Unsaturated Clay-Bearing Porous Media

Pore-Scale Hydro-Mechanical Couplings

In this study, we are interested in modelling multiphasic media made of a solid skeleton containing pores filled with a mixture of liquid and gas. The difference between gas and liquid pore pressures, s = u a -u w , is called suction. In the case in which air remains equal to the atmospheric pressure, water pressure is negative and suction takes a positive value. The air-water interface (called meniscus) starts to curve when suction increases. The radius of the meniscus decreases when suction increases, and once it becomes as small as the pore throats, air can invade the porous structure, which becomes unsaturated. The combination of the water surface tension and the negative pore water pressure results in a force that tends to pull the soil grains towards one another. The resulting force on the solid skeleton is similar to a compressive stress [START_REF] Santamarina | Soil Behavior at the Microscale: Particle Forces, in: Soil Behavior and Soft Ground Construction[END_REF]. An increase in suction will therefore lead to a decrease of the total volume (shrinkage), and wetting soils (i.e. decreasing suction) will usually make them swell. Suction also contributes to stiffen the soil against external loading thanks to grain bonding induced by water menisci in tension. The additional component of normal force at the contact will also prevent slippage between grains and thus increases the external force needed to cause plastic strains [START_REF] Ridley | The mechanical behaviour of an unsaturated compacted silty clay[END_REF]. However, when wetting a soil under constant mechanical loading, the resaturation destroys the bonds formed by water menisci and may induce an irrecoverable volumetric compression (called collapse) [START_REF] Muñoz Castelblanco | Some aspects of the compression and collapse behaviour of an unsaturated natural loess[END_REF]. These main characteristics of unsaturated soils mechanical behaviour are represented in figure 1. Changes in suction may also induce irreversible processes (plasticity or damage) during a drying process [START_REF] Wang | Irreversible deformation and damage in argillaceous rocks induced by wetting/drying[END_REF][START_REF] Alonso | Environmental degradation of claystones[END_REF].

State variables for unsaturated porous media

Unsaturated soil models are usually extensions of saturated soil ones. The most widely used of them is the Cam-clay model, first developed by [START_REF] Roscoe | On The Yielding of Soils[END_REF] and later modified by [START_REF] Roscoe | On the generalized stress-strain behaviour of wet clay[END_REF]. Extension to unsaturated states requires the definition of specific state variables. A comprehensive review of the existing stress frameworks can be found in the paper of [START_REF] Nuth | Effective stress concept in unsaturated soils: Clarification and validation of a unified framework[END_REF]. from [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF]. p ′ is the net mean stress. [START_REF] Houlsby | The work input to an unsaturated granular material[END_REF] demonstrated that, assuming the incompressibility of the solid matrix and the water phase, the work input to an unsaturated soil can be written as:

ẇ = [σ -(S r u w + (1 -S r )u a )I] : ε -(u a -u w )φ Ṡ r (1)
where σ is the total stress tensor, ε the total strain rate tensor, u a and u w the air and water pore pressures, φ the porosity, S r the degree of saturation, and I the identity matrix.

This formulation leads to the introduction of two state variables respectively conjugated to the strain rate, ε, and to the degree of saturation rate Ṡ r .

The stress quantity related to the strain increment is,

σ * = σ -(S r u w + (1 -S r )u a )I = σ -u a I + (u a -u w )S r I = σ net + sS r I (2)
which is a particular form of Bishop's effective stress [START_REF] Bishop | The effective stress principle[END_REF] in which the χ factor is taken equal to 1. This stress has been used by many authors and has been attributed different names, such as the average skeleton stress tensor [START_REF] Jommi | Remarks on the constitutive modelling of unsaturated soils[END_REF], the constitutive stress [START_REF] Sheng | Finite element formulation and algorithms for unsaturated soils. Part I: Theory[END_REF] or the generalised effective stress [START_REF] Laloui | On the use of the generalised effective stress in the constitutive modelling of unsaturated soils[END_REF]. In the following, the term constitutive stress will be used.

Other expressions have been proposed for this constitutive stress, accounting for the energy of the air-water interface [START_REF] Pereira | Adaptation of existing behaviour models to unsaturated states: application to CJS model[END_REF][START_REF] Nikooee | Effective Stress in Unsaturated Soils: A Thermodynamic Approach Based on the Interfacial Energy and Hydromechanical Coupling[END_REF], the different levels of porosity [START_REF] Alonso | A microstructurally based effective stress for unsaturated soils[END_REF], or the compressibility of the solid matrix through the Biot's coefficient [START_REF] Chateau | Micromechanics of saturated and unsaturated porous media[END_REF][START_REF] Jia | Elastoplastic damage modelling of argillite in partially saturated condition and application[END_REF]. However, for the sake of simplicity, we will keep the simple expression of equation 2, although the framework could easily accommodate a different expression for the constitutive stress.

According to equation 1, a second suction-related state variable, work-conjugated to the increment of degree of saturation is required. It will be called modified suction in the following and is written as:

s * = φs (3)
Constitutive equations are derived from an energy potential,

σ * = ∂ψ ∂ε e (4)
s * = - ∂ψ ∂S r (5)
in which ψ is the Helmholtz free energy, and ε e is the elastic strain tensor.

Due to the presence of the term sS r in the constitutive stress expression, the relationship between suction and degree of saturation has a great impact on the soil unsaturated mechanical behaviour.

For the sake of simplicity, hysteresis effects will be neglected in the present study, and the degree of saturation is expressed as a bijective function of the modified suction:

S r = f (s * ) (6)
The expression of [START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils[END_REF], in which the modified suction is used in place of suction, is used for the water retention properties:

S r = 1 1 + (α vg s * ) n vg m vg (7)
where α vg , n vg and m vg are parameters calibrated to fit experimental data.

It is worth noting that this stress framework choice provides many advantages in terms of numerical implementation. Indeed, for a degree of saturation S r = 1, equation 2 becomes σ * = σ -u w I, and Terzaghi's saturated effective stress is recovered. This smooth transition between unsaturated and saturated states makes it easy to simulate the behaviour of a soil submitted to negative as well as positive water pressures with a single model.

Dissipative mechanisms: Hydro-mechanical plasticity and damage models

Under deviatoric loading, clayey soils can undergo large permanent strains. Their properties, such as stiffness , strength, or permeability, are also known to be subject to changes after being submitted to hydric or mechanical solicitations. In clayey soils, these changes can be related to the destructuration of the material or to the deterioration of water bonds induced by tension in the menisci. Several approaches have been used to model this degradation. Some models have been developed which assume elastic moduli to be functions of the amount of plastic strains [START_REF] Hueckel | Coupling of elastic and plastic deformations of bulk solids[END_REF][START_REF] Sulem | Elastoplastic modelling of Red Wildmoor sandstone[END_REF][START_REF] Gajo | A model for stress and plastic strain induced nonlinear, hyperelastic anisotropy in soils[END_REF]. However, these models usually (except for [START_REF] Sulem | Elastoplastic modelling of Red Wildmoor sandstone[END_REF]) do not incorporate a strength reduction with plastic straining. Other models, developed for so-called structured, bonded or sensitive clays, focus on the increase in the size of the yield surface due to structure, which decreases at large strains to recover the yield surface of the reconstituted material [START_REF] Rouainia | A kinematic hardening constitutive model for natural clays with loss of structure[END_REF][START_REF] Kavvadas | A constitutive model for structured soils[END_REF][START_REF] Liu | A structured Cam Clay model[END_REF][START_REF] Nova | A constitutive model for bonded geomaterials subject to mechanical and/or chemical degradation[END_REF][START_REF] Baudet | A constitutive model for structured clays[END_REF][START_REF] Karstunen | Effect of Anisotropy and Destructuration on the Behavior of Murro Test Embankment[END_REF]. A parameter is then introduced to account for the degradation of structure, which evolves with plastic strains in the aforementioned models. These models do not account for the concomitant degradation of elastic stiffness, although it has been shown to decrease during loading for stiff clays such as claystones [START_REF] Chiarelli | Modeling of elastoplastic damage behavior of a claystone[END_REF]. Boom Clay samples extracted around a gallery have also shown a reduced small-strain shear modulus in the excavation damage zone [START_REF] Dao | Impact of excavation damage on the thermo-hydro-mechanical properties of natural Boom Clay[END_REF] Another approach, the one that is used in this paper, is to use the framework of Continuum Damage Mechanics, first developed for metals and later extended to concrete and rocks. This approach assumes that the degradation of material properties is due to the initiation and propagation of microcracks in rocks. This approach has been used for concrete behaviour modelling [START_REF] Grassl | Damage-plastic model for concrete failure[END_REF] as well as for stiff cemented clays [START_REF] Einav | Coupled damage and plasticity models derived from energy and dissipation potentials[END_REF]. Several approaches were proposed to model the evolution of stiffness and the accumulation of irreversible deformation induced by anisotropic damage [START_REF] Arson | Generalized stress variables in Continuum Damage Mechanics[END_REF]. Up to now, few attempts have been made to model damage in unsaturated geomaterials. Some models have been developed which consider damage in unsaturated geomaterials [START_REF] Arson | A mixed damage model for unsaturated porous media[END_REF], damage-plasticity couplings in saturated geomaterials [START_REF] Chiarelli | Modeling of elastoplastic damage behavior of a claystone[END_REF][START_REF] Conil | Poroplastic damage model for claystones[END_REF][START_REF] Yu | A Transversely Isotropic Damage Model for Boom Clay[END_REF], damage and viscoplasticity in unsaturated geomaterials [START_REF] Dufour | A thermodynamically consistent framework for saturated viscoplastic rock-materials subject to damage[END_REF], and even damage-plasticity couplings in unsaturated geomaterials [START_REF] Hoxha | Saturated and unsaturated behaviour modelling of Meuse-Haute/Marne argillite[END_REF][START_REF] Jia | Elastoplastic damage modelling of argillite in partially saturated condition and application[END_REF]. However, these models, initially formulated for rocks, ignore some specific important features of clayey soil behaviour, such as the dependence of elastic moduli to pressure. Moreover, damage-plasticity models proposed for rocks so far fail at predicting the transition between ductile and brittle behaviour associated with suction increase. [START_REF] Vaunat | Bond degradation and irreversible strains in soft argillaceous rock[END_REF] developed a model for bonded granular soils, based on microstructural considerations, which is able to reproduce both strength and stiffness degradation coupled with elastoplasticity. It has been later extended to other materials and loading scenarios [START_REF] Pinyol | A constitutive model for soft clayey rocks that includes weathering effects[END_REF][START_REF] Yang | A constitutive model for unsaturated cemented soils under cyclic loading[END_REF][START_REF] Cardoso | A constitutive model for compacted expansive and bonded marls[END_REF]. This model, however, has been developed for a very specific class of geomaterials. By contrast, our modelling framework can be applied to a broader range of materials, from quasi-brittle stiff clays for which damage occurs in the elastic strain domain, to soft clays for which plasticity is the predominant dissipative mechanism.

Principle of effective stress in Continuum Damage Mechanics

Introduced by [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF], the effective stress in the sense of damage mechanics is based on the fact that the resisting section decreases when micro-cracks develop. This approach can be extended to a broad range of unsaturated clay-bearing geomaterials, in which damage also represents the loss of bonding due to water menisci. Note for instance that in some thermodynamic frameworks [START_REF] Coussy | Revisiting the thermodynamics of hardening plasticity for unsaturated soils[END_REF], air-water interfaces are part of the apparent solid skeleton.

A scalar damage variable, d, is defined as an average of the proportion of damaged surfaces in the material. d ranges from d = 0 for an intact material to d = 1 for a totally damaged material with no residual resistance. Assuming that damage is isotropic and affects similarly all components of the stress tensor, the effective stress tensor then becomes

σ = σ 1 -d (8)
More complex expressions could be used in place of equation 8, in order to accommodate more sophisticated behaviours, such as anisotropic damage.

In this study, damage is a phenomenological variable that describes the combined effects of multiple mechanisms on elastic stiffness degradation.

Hydro-mechanical damage-plasticity model based on the concept of double effective stress

Introduction of a double effective stress accounting for suction and damage effects

The two previous sections allowed us to introduce two quantities describing the stress applied on the solid matrix. On the one hand, the constitutive stress, in unsaturated soils, takes into account the effect of water menisci in tension, acting like a compressive stress on the solid matrix. On the other hand, the effective stress, in the sense of damage mechanics, enables us to account for the decreasing material surface sustaining mechanical loads, resulting from the creation of microcracks and the modification of water bonds. There is a need to define a new quantity, representing the stress applied on the solid matrix when the material is affected by both suction and damage simultaneously. This quantity will be called the double effective stress and is assumed to control the porous material mechanical behaviour.

Two simple combinations of the previous effective stresses can be imagined to incorporate both damage and suction into this double effective stress, σ * :

σ * 1 = σ -u a I + sS r I 1 -d = σ * 1 -d (9) σ * 2 = σ 1 -d -u a I + sS r I = σ -u a I + sS r I (10)
To choose between these two expressions, we assume that a damaged sample submitted to a change in suction should behave differently compared to the intact sample. This hypothesis has been considered by other authors, such as [START_REF] Carmeliet | Poromechanical modelling of shrinkage and damage processes in unsaturated porous media[END_REF], who consider that damaged materials experience more swelling when wetted that intact ones.

Assuming that the total applied stress, the gas pressure, as well as damage are kept constant, the change in the double effective stress due to a suction increment would be:

σ * 1 = (s Ṡ r + S r ṡ)I 1 -d (11) σ * 2 = (s Ṡ r + S r ṡ)I (12)
Thanks to the use of a hyperelastic formulation (section 3.3), elastic strains are directly related to the double effective stress. The strain change due to suction change would therefore be the same for an intact and a damaged sample for the second expression.

We will thus choose the first expression for the double effective stress:

σ * = σ -u a I + sS r I 1 -d (13)
We define the following quantities:

• Mean stress: p = 1 3 tr(σ)

• Deviatoric stress tensor:

σ d = σ -pI • Deviatoric stress: q = 3 2 σ d : σ d
Then the double effective triaxial variables are:

p * = p -u a + sS r 1 -d (14) q * = q = q 1 -d (15)
It can be noted that with this definition of the double effective stress, suction effects are isotropic, and thus don't have any impact on the deviatoric stress.

The existence of a double effective stress, in which suction and damage effects on mechanical behaviour are included, is a key assumption in the following modelling developments. In the following sections, we will study how this double effective stress allows for damage and suction effects on elastic and dissipative behaviours to be reproduced.

Expression of Helmholtz free energy

We assume that the material state is described by the values of the following state variables:

The elastic strain, ε e , the degree of saturation, S r , damage, d, and a hardening variable, χ. We assume that elastic, plastic and hydraulic potential energy functions are decoupled and that processes are isothermal. We propose the following form for Helmholtz free energy :

ψ = ψ(ε e , S r , d, χ) = ψ e (ε e , d) + ψ l (S r ) + ψ p (ε p , χ) (16)
In order to build a damage constitutive model, an extra assumption has to be added to the concept of effective stress. Concerning the damage-elastic part of Helmholtz free energy, we choose to use the form proposed by [START_REF] Ju | On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects[END_REF],

ψ e (ε e , d) = ψ e 0 (ε e )(1 -d) (17)
which, after derivation gives the following expression of the constitutive stress:

σ * = ∂ψ e ∂ε e = (1 -d) ∂ψ e 0 ∂ε e (18)
The double effective stress is therefore related to elastic strains through the following constitutive relationship:

σ * = σ * 1 -d = ∂ψ e 0 ∂ε e (19)
The relationship given in equation 19 implies that, in a damaged material, the double effective stress will be linked to elastic strains with the same relationships that the constitutive stress in an intact material. This is the principle of strain equivalence defined by [START_REF] Lemaitre | Aspect phénoménologique de la rupture par endommagement[END_REF] which states that the strain associated with a damaged state under the applied stress is equivalent to the strain associated with its undamaged state under the effective stress. The principle of strain equivalence is illustrated in figure 2.

This approach has the advantage of being easily extensible to plasticity, by replacing stresses by double effective stresses in classical equations.

Elasticity

For the sake of simplicity, the elasticity is assumed to be linear in the following developments.

However, experimental evidence show that bulk and shear moduli of geomaterials increase with confining pressure, which may have an important effect on the material behaviour, especially if a large confining pressure range is considered. The present framework can be adapted to non linear elasticity. In order to ensure the conservation of the elastic deformation energy, it is necessary to formulate the model within the framework of hyper-elasticity [START_REF] Zytynski | On modelling the unloading-reloading behaviour of soils[END_REF]. Challenges related to the degradation of pressure dependent elastic moduli in porous material were discussed in Le Pense (2014).

Incorporating linear elasticity into equation 19 gives:

         p * q          =           K 0 0 3G                    ε e v ε e s          (20)
in which K is the bulk modulus, G the shear modulus, ε e v the volumetric elastic strain, and ε e s the deviatoric elastic strain. This gives the following apparent stiffness matrix when expressed in terms of constitutive stresses:

         p * q          =           K(1 -d) 0 0 3G(1 -d)                    ε e v ε e s          (21)
Coupling of the principle of strain equivalence with a damaged effective stress therefore leads to a degradation of apparent elastic moduli with damage, without the need to explicitly express them as functions of the damage parameter.

Damage onset and evolution

Since suction has no effect on deviatoric stress, we will not consider the anisotropy induced by damage in this paper. We adopt Drucker-Prager damage criterion, which is expressed in terms of double effective stresses, so as to follow the principle of strain equivalence: in which C 1 is a hardening parameter. The lower C 1 , the faster d will increase with deviatoric stress. The C 2 coefficient allows the dependence on confining pressure to be accounted for. Indeed, geomaterials are known to be more brittle at low confining pressure and more plastic at high confining pressures. C 0 enables the modification of the damage threshold.

f d = q -C 2 p * -C 0 -C 1 d = 0 (22) a. p ~* q ~--C 0 +C 1 d C 2 d ➚
The shape of the damage criterion in the double effective stress space is given in figure 3a for different values of damage. It can be seen that, when damage increases, the material is hardening with respect to effective stresses.

Expressed in total stresses, equation 22 becomes:

q 1 -d -C 2 p + sS r 1 -d -C 0 -C 1 d = 0 (23) q -C 2 (p + sS r ) -(1 -d)(C 0 + C 1 d) = 0 (24)
The shape of the corresponding damage criterion in the total stress space is given in figure 3.

Figure 3b shows the evolution of the damage criterion with suction. Although suction does not have an effect on the damage criteria in the double effective stress space, suction increases the stress value for which damage is initiated when considering total stresses.

Figure 3c shows the evolution of the damage criterion with damage. It can been seen that, although the intact fraction of the material is hardening, an apparent softening behaviour appears after a certain value of damage is reached, when considering total stresses.

Deriving equation 22 gives the consistency condition,

0 = ḟd = ∂ f d ∂ p * ṗ * + ∂ f d ∂ q q + ∂ f d ∂d ḋ ( 25
)
from which the damage evolution law can be deduced:

ḋ = A d ( σ) : σ * ( 26
)
where

A d = 1 C 1 - C 2 3 I + σd q .
This expression of the damage evolution rate as well as the damage criterion implies that damage initiation and evolution are solely related to elastic strains. Note that some authors assumed that damage is initiated by an accumulation of plastic strains. By contrast, we decoupled damage and plasticity, which allows modelling a wide range of materials, subject to damage propagation only, or plastic dissipation only, or both damage and plastic dissipation.

Coupled damage and plasticity model: suction hardening and damage softening

According to [START_REF] Jommi | Remarks on the constitutive modelling of unsaturated soils[END_REF], extending a poromechanical model from saturated to unsaturated materials requires the two following steps:

• the substitution of the average skeleton stress for effective stress

• introduction in the basic saturated elastoplastic model of the modifications necessary to take into account the effects of the interfaces on the overall mechanical behaviour According to [START_REF] Ju | On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects[END_REF], plasticity occurs only in the undamaged counterpart of the bulk, and the expression of plastic flow for the damaged material can be obtained by using effective properties and effective stress in the expression of plastic flow for undamaged materials. Therefore, the characterisation of the plastic response should be formulated in the damaged effective stress space and the stress tensor should be replaced by the damaged stress tensor, σ, into the equations of plasticity. This follows the principle of strain equivalence.

Similarly to the Barcelona Basic Model (BBM) [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF], the most widely used model for unsaturated soils, we use the modified Cam-Clay model (Burland, 1965) as a basis to predict plasticity in saturated geomaterials. Based on Jommi's and Ju's recommendations, we formulate the yield criterion in terms of double effective stress to extend the model to damaged and unsaturated geomaterials. In addition, we introduce a dependence of the yield criterion to suction:

The yield surface is therefore taken of the following form:

f p = q2 -M 2 p * (p * c (p 0 , s) -p * ) ( 27
)
in which p * c is the preconsolidation pressure, which is a function of suction and the saturated preconsolidation pressure, p 0 (equation 34).

Cam-clay models have been developed in the framework of Critical State Soil Mechanics [START_REF] Roscoe | On The Yielding of Soils[END_REF]. The critical state concept states that soils and other granular materials, if continuously distorted until they flow as a frictional fluid, will come into a well-defined critical state. At the onset of the critical state, shear distortions occur without any further changes in mean stress, deviatoric stress or void ratio. The critical state is described in the ( p * , q) plane by the line of equation:

q = M p * (28)
In some recent models (such as BBM), non-associate flow rules are adopted to predict plastic volumetric strains under a variety of stress paths. For the sake of simplicity, we considered an associate flow rule (like in the original Cam-Clay model). The plastic potential is defined as:

g p = f p = q2 -M 2 p * (p * c -p * ) (29)
The plastic flow rule is:

εp = Λp ∂g p ∂ σ * = Λp ∂g p ∂ p * I 3 + ∂g p ∂q 3 σd 2 q (30) in which ∂g p ∂ σ * = M 2 (2 p * -pc * ) 3 I + 3 σd (31)
The hardening law is defined as: In order to reproduce the extension of the elastic domain with suction, the preconsolidation pressure is sought in the form of a function of suction and saturated preconsolidation pressure (p 0 ):

ṗ0 = p 0 λ -κ εp v ( 32
p * c = p * c (p 0 , s) (33)
When drawn in the (p, s) plane, this curve is called the Loading-Collapse (LC) curve. Many different expressions were proposed for the equation of the LC curve [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF][START_REF] Jommi | Remarks on the constitutive modelling of unsaturated soils[END_REF][START_REF] Buisson | Coupling of hydraulic hysteresis and stressstrain behaviour in unsaturated soils[END_REF][START_REF] Sheng | A constitutive model for unsaturated soils : thermomechanical and computational aspects[END_REF][START_REF] Sun | Elastoplastic prediction of hydro-mechanical behaviour of unsaturated soils under undrained conditions[END_REF]. Experimental characterisation of clay-bearing unsaturated materials is often driven by the determination of BBM parameters. In order to facilitate the calibration of our model, we chose the LC equation proposed by [START_REF] Sheng | A constitutive model for unsaturated soils : thermomechanical and computational aspects[END_REF], because Sheng et al.'s approach is the closest we found to the BBM model. The equation of the LC curve is:

p * c = p r p 0 p r λ-κ λs-κ + sS r (34) λ s = λ[(1 -r) exp(-βs) + r] ( 35
)
The shape of the yield criterion in the double effective stress space is given in figure 4a. As expected, the yield surface in the double effective stress space does depend on suction, but not on damage.

Expressed in total stresses, the equation of the yield surface (equation 27) becomes

q 2 -M 2 (p + sS r ) (1 -d)p * c -(p + sS r ) = 0 (36)
and the equation of the critical state line becomes

q = M(p + sS r ) (37)
Figure 4b shows the evolution of the yield surface with suction. With respect to total stresses, the elastic domain increases with suction. Suction also induces an apparent cohesion.

Figure 4c shows the evolution of the yield surface with damage. Damage has a softening effect on the plastic behaviour. Although plastic and damage dissipative potentials were assumed to be decoupled, the assumption of a double effective stress, associated with the principle of strain equivalence, allows for a direct damage-plasticity coupling. Indeed, although damage and plasticity criteria are expressed in terms of the double effective stress, and consequently do not depend explicitly on damage and suction, they evolve with damage and suction in the total stress space.

. The following section will illustrate how the proposed model behaves for its mechanical part, based on specific sets of parameters.

Illustration of the mechanical damage-plastic behaviour

The model has been designed in a flexible way, which enables the independent refinement of its basic components (effective stresses, elasticity, damage and plasticity equations) to fit specific materials behaviours. Analysis of the model behaviour, as well as its validation, will focus on clayey geomaterials, such as Boom clay, since these materials exhibit simultaneously a strong plastic behaviour, as well as damage.

Summary of Boom Clay data from the literature

Boom Clay has been selected as a possible host rock for deep radioactive waste disposal in Belgium. It is considered as an overconsolidated plastic clay.

Most experimental studies published on Boom Clay focus on the characterisation of physical properties (such as retention and permeability), or on the hydro-mechanical response of the soil in unsaturated conditions. Very few measures were done to study the degradation of stiffness with stress and suction, despite the proven existence of an Excavation Damaged Zone (EDZ) around underground openings. Excavation induced fractures were observed around galleries [START_REF] Bastiaens | SELFRAC: Experiments and conclusions on fracturing, selfhealing and self-sealing processes in clays[END_REF][START_REF] Bernier | Twenty-five years' geotechnical observation and testing in the Tertiary Boom Clay formation[END_REF][START_REF] Van Marcke | Excavation induced fractures in a plastic clay formation: Observations at the HADES URF[END_REF]. Damage was also studied by means of seismic [START_REF] Bastiaens | SELFRAC: Experiments and conclusions on fracturing, selfhealing and self-sealing processes in clays[END_REF] and acoustic [START_REF] Lavrov | Acoustic emission in host-rock material for radioactive waste disposal: comparison between clay and rock salt[END_REF] measurements. Boom clay can exhibit both ductile and brittle behaviours [START_REF] Dehandschutter | Brittle fractures and ductile shear bands in argillaceous sediments: inferences from Oligocene Boom Clay (Belgium)[END_REF]. More recent studies using advanced imaging technique provided evidence of cracks in Boom Clay samples [START_REF] Bésuelle | A Laboratory Experimental Study of the Hydromechanical Behavior of Boom Clay[END_REF] The transition between the failure modes depends strongly on the confining pressure and is also influenced by the water content [START_REF] Al-Shayea | The combined effect of clay and moisture content on the behavior of remolded unsaturated soils[END_REF] and by the overconsolidation ratio.

Boom clay has been extensively studied either from experiments on undisturbed natural samples, or on samples prepared by compaction from Boom clay powder. Many experimental data on saturated natural Boom clay are available in the literature [START_REF] Baldi | Developments in modelling of thermo-hydro-geomechanical behaviour of Boom clay and clay-based buffer materials[END_REF][START_REF] Coll | Endommagement des roches argileuses et perméabilité induite au voisinage douvrages souterrains[END_REF][START_REF] Sultan | Yielding and plastic behaviour of Boom clay[END_REF]. However, concerning the unsaturated behaviour, most of the studies have been made on compacted [START_REF] Bernier | Suction-controlled experiments on Boom clay[END_REF][START_REF] Romero | Characterisation and thermo-hydro-mechanical behaviour of unsaturated Boom clay: an experimental study[END_REF] or remoulded [START_REF] Al-Mukhtar | The fabric of a clay soil under controlled mechanical and hydraulic stress states[END_REF] samples, and only a few on undisturbed samples [START_REF] Cui | Suction effects in deep Boom Clay block samples[END_REF][START_REF] Vecchia | Some remarks on the hydro-mechanical constitutive modelling of natural and compacted Boom clay[END_REF].

Moreover, mechanical tests at different suctions are limited to oedometer and isotropic compression tests. By comparing experiments on natural and compacted samples, [START_REF] Vecchia | Some remarks on the hydro-mechanical constitutive modelling of natural and compacted Boom clay[END_REF] concluded that the same constitutive framework seems to be applicable to natural Boom clay and to the material compacted from the clay powder. However, mechanical parameters have to be adapted for different microstructures.

Boom clay is a more complex material than other clay stones such as Callovo-Oxfordian argillites, which exhibit a less plastic behaviour. The following simulations will demonstrate that the modelling approach that we proposed above is suitable to predict the brittle/ductile transition in unsaturated geomaterials in which the behaviour is strongly influenced by plastic deformation, confining pressure, and water content.

Data available to calibrate our model involves tests performed on cores of different origins, taken at different depths. The mineral composition of the samples varied greatly from one experiment to the other, and therefore, a high variability was noted in the mechanical and physical properties. In the following numerical study, we used material parameters that fell in the range of values reported in the literature, and we adapted the set of parameters to the different soils tested, in order to match experimental test results. Some values found in the literature for elasticity, plasticity, and retention parameters are given in table 1. No similar damage model has been found which would allow to determine the range of values for our damage parameters.

Table 1: List of the model parameters and their range of values as found in the literature Franc ¸ois et al. The set of data chosen to study the sensitivity of the damage-plastic model to the different parameters is given in table 2.

Elasticity

Plasticity Damage Initial

state

E ν M λ -κ p 0 C 0 C 1 C 2 p MPa MPa MPa
MPa MPa 300 0.4 1 0.05 5.5 0 10 0.5 3

Table 2: Parameters chosen as a basis for the parametric study

Damage model

Although the damage part of the model has been chosen to be formulated with the minimum number of parameters, and to be based on Drucker-Prager, no identical model, expressed in terms of the damaged effective stress has be found in the literature. The behaviour of this model, and the range of parameters for which it gives sensible results is studied in this section.

It can be seen in figure 5a that the model can exhibit a radial contraction under triaxial loading for certain sets of parameters.

This feature appears for sets of parameters which do not respect equation 38 (details of the calculation are given in Appendix A):

3(1 -2ν) 1 + ν C 2 3 (38)
in which ν is the Poisson's ratio, and C 2 the slope of the damage criterion.

Figure 5 shows the effect of the damage parameters C 2 (slope of the damage criterion) and C 1 (hardening parameter) on the stress-strain curves as well as the evolution of damage with axial strain. As expected from the theoretical developments presented earlier, the damage model can reproduce a hardening behaviour followed by a softening behaviour. It can be seen that variations of C 2 result mainly in a modification of the damage threshold, whereas changing C 1 modifies the damage evolution rate. High values of C 1 therefore result in higher peak stress values.

Damage-plasticity coupling

The plasticity part of the model is similar to BBM, expressed in terms of the unsaturated constitutive stress. It has been widely studied in the literature and sensitivity to its parameters will not be detailed here. 20 when only damage or plasticity is considered. Two cases are presented, one for which plasticity is the dominant dissipative phenomena, using the parameters of table 2 (figure 6a,c,e), and one for which damage is dominant (figure 6b,d,f).

0 1 2 3 4 5 6 7 8 -4 -2 0 2 4 6 8 10 q (MPa) Effect of C 2 Effect of C 1 a. b. c. d. ε a (%) ε r (%) C2=0.3 C2=0.5 C2=0.7 -4 -2 0 2 4 6 8 10 Effect of C 2 Effect of C 1 a. b. c. d. ε a (%) ε r (%) ε r (%) ε a (%) C1=5 MPa C1=10 MPa C1=20 MPa 0 25 50 75 100 -4 -2 0 2 4 6 8 10 damage (%) ε a (%) Effect of C 2 Effect of C 1 a. b. c. d. ε a (%) ε r (%) ε r (%) ε a (%) -4 -2 0 2 4 6 8 10 ε a (%) Effect of C 2 Effect of C 1 a. b. c. d. ε a (%) ε r (%) ε r (%) ε a (%)
For the plasticity dominated case, it can be seen on figure 6a that the softening behaviour after deviatoric stress peak, characteristic of the damage model, is absent for the coupled model. The coupled model stress-strain behaviour also follows the same pattern as the plasticity model. This shows that, when using model parameters suitable for Boom clay, plasticity dominates damage effects. Figure 6c shows, however, that damage is triggered, and develops up to 20%. When looking only at the stress-strain curve, one could mistake the non-linear behaviour as being the results of plasticity effects only. It should therefore be noted that the analysis of the strain-stress measurement only could hide the appearance of damage. The main effect of damage in that case, is to decrease the apparent yield stress of the material, as seen in figure 6a, which has a negligible effect on the final amount of plastic strain (figure 6e).

In order to demonstrate the capabilities of the model, figure 6b shows the behaviour of the model when damage dominates plasticity effects. The parameters chosen for this example are mostly the one of table 2, but with C 1 = 5 MPa and λκ = 0.01, values chosen to increase the influence of damage, and decrease the influence of plasticity. These results show that the presented model can also reproduce a damage dominated behaviour, with a stress-strain behaviour of the coupled model similar to the one of the damage model. Plasticity effects result in greater strains, but have no effect on the peak deviatoric stress.

These illustrative simulations show that the current model is highly versatile and, depending on the set of parameters chosen, can reproduce damage-plasticity couplings, dominated either by plasticity or damage behaviours.

Simulation of hydro-mechanical experiments on Boom Clay

This section aims at comparing simulation results, using the model developed in the previous sections, with hydro-mechanical experiments results on clayey soils from the literature. As mentioned previously, parameters are adjusted to fit specific tests, but are chosen to lie within the range of values reported in the literature.

Elastic swelling

A significant advantage of the constitutive stress approach is the ability to capture suction induced strains without the need of extra parameters in addition to mechanical and retention parameters.

To illustrate this feature, an oedometer swelling test was simulated on a clay material (experimental data from [START_REF] Volckaert | Thermal-hydraulic-mechanical and geochemical behaviour of the clay barrier in radioactive waste repositories (model development and validation[END_REF]). The vertical stress was kept constant (σ v = 0.1 MPa) while suction was decreased from 230 to 0 MPa. The mechanical and retention properties chosen in the simulation are given in table [START_REF] Carmeliet | Poromechanical modelling of shrinkage and damage processes in unsaturated porous media[END_REF]. Indeed, the test is suction-controlled, therefore the volume of the sample can change as water tends to fill the pores during the wetting phase. Damaged samples are more compliant than undamaged samples: the resistance of the solid skeleton to pore filling and expansion is less in damaged materials, which tend to swell more than undamaged samples during wetting. This behaviour has been observed in oedometric swelling experiments on callovo-oxfordian argillite samples [START_REF] Mohajerani | Oedometric compression and swelling behaviour of the Callovo-Oxfordian argillite[END_REF], where it is seen that swelling capacity increases with damage.

Triaxial tests on saturated samples at different confining pressures

Triaxial drained compression tests with unloading-reloading cycles are simulated for two confining pressures (3 MPa and 4 MPa). The experimental data (from [START_REF] Baldi | Developments in modelling of thermo-hydro-geomechanical behaviour of Boom clay and clay-based buffer materials[END_REF]) show the influence of the confining pressure on the deviatoric response. A degradation of the elastic modulus can also be seen from the unloading-reloading curves.

Elastic, plastic and damage parameters chosen in this study are summarised in table 4. The preconsolidation pressure is taken equal to 6 MPa, which is in the range of values observed on samples from underground laboratories. The other mechanical parameters are chosen to fit the experimental results reported in [START_REF] Baldi | Developments in modelling of thermo-hydro-geomechanical behaviour of Boom clay and clay-based buffer materials[END_REF].

Figure 8 shows the comparison between experimental and numerical stress/strain curves. The main trends observed in the laboratory are captured by the model. For instance, it is noted that the stiffness measured during the unloading paths is less than the stiffness measured during the first loading paths. As expected, the loading stress supported by the sample before damage propagation is higher at higher confining pressure. However, the model does not capture well the smooth transition between elastic and plastic behaviour. This limitation of the model can be explained by the use of Cam-clay model, in which elasticity is assumed for all states of stress inside the yield surface. This behaviour could be improved by using more advanced versions of the Cam-Clay model, such as bounding surface plasticity [START_REF] Dafalias | Bounding Surface Plasticity. I: Mathematical Foundation and Hypoplasticity[END_REF] or continuous hyperplasticity [START_REF] Puzrin | Fundamentals of kinematic hardening hyperplasticity[END_REF]. The volumetric behaviour could also be improved by using non-associated flow rules.

Figure 8g shows the corresponding stress paths in the double effective stress space. It can be seen that for low confining pressure, the stress path attains the damage criterion earlier, which allows for more damage to be developed before the critical state is reached. The activation of the two competitive dissipation phenomena, damage and plasticity, depends on the confining pressure.

Simulation of the ductile/brittle transition with suction increase

Although no experimental data have been found in the literature about Boom clay, [START_REF] Al-Shayea | The combined effect of clay and moisture content on the behavior of remolded unsaturated soils[END_REF] showed that materials with high clay content exhibit a ductile/brittle behaviour transition when their water content decreases (see figure 9a). Ductile behaviour is characterised by the ability to sustain large plastic strains during plastic hardening. Brittle behaviour is characterised by abrupt failure at a well-defined peak strength with strong softening. Figure 9a also shows higher shear strength for low water contents.

Although the experimental data from Al-Shayea ( 2001) are difficult to interpret and the exact experimental procedure can not be reproduced in simulation due to the lack of data (retention properties, preconsolidation pressure, unloading-reloading curves), we will show that our model can reproduce a similar transition between a ductile and brittle behaviour when suction increases.

Triaxial compression tests under constant suction (0 MPa, 0.5 MPa, 1 MPa) are simulated. The confining pressure is taken equal to 200 kPa, and since samples are compacted in the experiments taken as reference, the preconsolidation pressure is taken equal to 500 kPa. The complete list of parameters chosen for this study are given in table 5. The stress-strain curves obtained for different suctions are given in figure 9b. The corresponding effective stress paths can be seen on figure 9e.

The evolution of damage and plastic strains with axial strain is given in figure 9c-d. Figure 9b shows that our model can adequately reproduce the transition from a ductile behaviour for low suction, to a brittle behaviour for higher suctions.

At low suction, the plastic yield stress is low, and the plastic criterion is reached before the damage criterion. This leads to the development of large plastic strains, and damage remains low because of the lack of increase of the deviatoric stress . At higher suctions, the elastic domain is enlarged. The damage criterion is therefore reached before the plastic criterion. The deviatoric stress, and therefore damage, reaches higher values before the triggering of plasticity.

Discussion and conclusions

A constitutive modelling framework allowing for damage-plasticity couplings in unsaturated porous media has been proposed. This framework is based on the assumption of a double ef- The presented modelling framework, based on the combined assumptions of the existence of a double effective stress and the principle of strain equivalence, presents many advantages. The numerical implementation is straightforward and is able to accommodate different plasticity and damage models without the need of heavy code modifications.

Once implemented into a finite element code, this modelling framework will enable the mod-

Figure 1 :

 1 Figure 1: Influence of suction on volumetric compression and volume changes due to wetting and drying, adapted

Figure 2 :

 2 Figure 2: Principle of strain equivalence

Figure 3 :

 3 Figure 3: Shape of damage criterion. (a) Double effective stress space, (b and c) Total stress space (b. effect of suction, c. effect of damage)

Figure 4 :

 4 Figure 4: Shape of yield surface. (a) Double effective stress space, (b and c) Total stress space (b. effect of suction, c. effect of damage)

a

  Results for natural Boom Clay, parameters for drying and wetting curves b Calculated from K = k i (1 + e)p (non linear elasticity) for an initial state e = 0.59 and p = 4.4 MPa c equivalent to λκ 1 + e in the cited references

Figure 5 :

 5 Figure 5: Effect of C 2 (a,c) and C 1 (b,d) on the stress-strain curves (a,b) and damage evolution as a function of axial strain (c,d)

Figure 6 :

 6 Figure 6: Comparison of damage, plastic, and coupled model behaviour, for a plasticity dominated (a,c,e) and a damage dominated (b,d,f) coupling. (a,b) Stress-strain curve, (c,d) Damage evolution, (e,f) Plastic strain evolution

Figure 7 :

 7 Figure 7: Volumetric swelling strains in oedometric conditions under a vertical load σ v = 0.1 MPa (compared with experimental data from Volckaert et al. (1996))

Figure 8 :

 8 Figure 8: (a-d) Triaxial test data from Baldi et al. for confinements (a,c) 3 MPa and (b,d) 4 MPa, compared with simulation results (a-f). (g) Double effective stress paths

Figure 9 :

 9 Figure 9: Stress-strain behaviour of a clay-sand mixture (40 % clay) at different moisture contents: (a) Experimental data from Al-Shayea (2001). (b) Present model. (c) Damage evolution. (d) Plastic strains. (e) Double effective stress path

Table 3 :

 3 3. Elasticity and retention parameters for the swelling testSwelling strains (volumetric strains) computed for an intact as well as a damaged material are represented on figure 7. Damage is assumed to remain constant during the test, and the suction state to be homogeneous within the sample.Knowledge of the water retention properties in addition to mechanical stiffness parameters allows us to reproduce adequately the elastic swelling behaviour observed during wetting. Moreover, a different swelling behaviour is observed for intact and damaged samples, which is in accordance with other works

	Elasticity		Retention
	K	ν	S r	α vg	n vg	m vg
	MPa			kPa -1
	200	0.3	0	0.28.10 -3 2.3	0.21

Table 4 :

 4 Boom Clay mechanical parameters

Acknowledgements

The authors would like to thank ENS Cachan for funding the PhD thesis of the first author.

Appendix A. Damage Model, limits for which the sample experiences radial contraction during a triaxial test

For a triaxial stress state, the elasticity law reads:

in which E is the Young's modulus, and ν the Poisson's ratio. Subscript a is related to the axial direction, and r to the radial direction.

Damage increment for a triaxial stress state:

We want ε r 0. Since ḋ 0:

We can see that this expression is always true when there is no confinement, i.e. for σ 0 = 0 (true if C 2 3).

Otherwise:

(

If one wants it to be true for complete damage, i.e. d = 1, this gives the following relationship between the slope of the damage criterion, C 2 , and the Poisson's ratio, ν:

For a given set of parameters, it is also possible to determine from which value of damage radial contraction will start: .12)


Cristian Ruiz

Joseph Emeras

Emmanuel Jeanvoine

Lucas Nussbaum

Distem: Evaluation of Fault Tolerance and Load Balancing Strategies in Real HPC Runtimes through Emulation

Keywords: Experimentation, HPC runtimes, Fault tolerance, Load balancing, Emulation

The era of Exascale computing raises new challenges for HPC. Intrinsic characteristics of those extreme scale platforms bring energy and reliability issues. To cope with those constraints, applications will have to be more flexible in order to deal with platform geometry evolutions and unavoidable failures. Thus, to prepare for this upcoming era, a strong effort must be made on improving the HPC software stack. This work focuses on improving the study of a central part of the software stack, the HPC runtimes. To this end we propose a set of extensions to the Distem emulator that enable the evaluation of fault tolerance and load balancing mechanisms in such runtimes. Extensive experimentation showing the benefits of our approach has been performed with three HPC runtimes: Charm++, MPICH, and OpenMPI.

I. INTRODUCTION

During the last decade, numerical simulations have become a key component of most advances in today's society, science and engineering. Those simulations, computing interactions of many complex models, leverage the computing power of modern HPC infrastructures. The increasing scale and volume of those simulations is a driver for infrastructure evolution, and before 2020, HPC infrastructures will likely reach Exascale. Several projects around the world have studied the challenges of computing at Exascale and numerous improvements need to be achieved. In particular, The International Exascale Software Project [START_REF] Dongarra | The International Exascale Software Project Roadmap[END_REF] explored those challenges from the software stack point of view. One identified direction is the improvement of the HPC runtimes to ease the programming of Exascale applications. Among other directions, the report recommends to focus on the support of platform heterogeneity, load balancing, and resilience. Experimentation is fundamental in the HPC context. Evaluating current HPC runtimes under complex realistic conditions is critical for application developers and researchers.

This work contributes to enabling experimental evaluation of HPC runtimes. Starting from Distem, which is a versatile emulator for studying distributed systems, we designed an emulator suitable for the evaluation of HPC runtimes, enabling specifically: (1) introduction of heterogeneity and dynamic imbalance among the computing resources; (2) introduction of failures. Those features provide runtime designers with the ability to experiment their prototypes under a large range of conditions. This paper is organized as follows. First, Section II describes how the Distem emulator was extended to address the specific challenges of the evaluation of HPC runtimes. Then, Section III presents experiments that demonstrate the benefit of such a tool for HPC runtime evaluation. Finally, Section IV compares our approach to related works, and Section V concludes this paper and discusses our future work.

II. DESIGN OF AN EMULATOR TO STUDY HPC RUNTIMES

This section first highlights the important points when studying HPC runtimes and presents the difficulties when evaluating them. Then, it presents the design of Distem and how it enables testing HPC runtimes in heterogeneous environments, possibly under evolving conditions and in presence of failures.

A. Study of HPC runtimes

To design interesting features for studying HPC runtimes, we first have to delve what is important to achieve such studies and why it is difficult to achieve.

Compared to typical computing environments, the main characteristic of HPC environments is the high number of computing resources and its critical density. Furthermore, typical HPC workloads imply a large amount of computations that are likely to be executed for several days or months. Because of the high number of resources, the probability that failures hit one or more resources while computing is very high. High density of the resources also increases failures since it is complicated to extract the generated heat from computers and since heat is known to damage electronic components. Thus, it is quite obvious that HPC runtimes have to deal with failures.

To cope with the heat, hardware manufacturers and OS developers offer several solutions that usually aim at reducing the computing power dynamically (for instance by reducing the frequency of some CPU cores, or by disabling some cores). For classical applications this has almost no impact but parallel codes usually used in HPC applications would be disrupted. Indeed, a typical parallel application runs the same workload on all the computing units using network communications to synchronize and exchange required data among the computing units. To avoid inactivity of some computing CPU while other are performing a computation with a reduced computing power, HPC runtimes also have to implement load-balancing strategies. In this paper, we will focus on the way to study fault-tolerance and load-balancing features in HPC runtimes.

When implementing such features, HPC runtime designers and developers tackle the difficulty to ensure that runtimes are working as expected and to tune them according to the execution platform.

Regarding fault-tolerance, fault injectors can be used to validate the recovery mechanisms. Those injectors can be embedded in the runtime but this often suffer from a lack of realism, like the one implemented in Charm++ [START_REF] Kale | CHARM++: A Portable Concurrent Object Oriented System Based on C++[END_REF]. Indeed, the runtime can stop some running processes but it cannot simulate an hardware or a network failure since it does not control the kernel. Some frameworks exist with various level of realism. However, realistic ones imply intrusive kernel modifications [START_REF] Coti | Os-level failure injection with systemtap[END_REF]. Regarding load-balancing, ad-hoc mechanisms are usually used, leveraging CPU scaling to reduce the frequency of some CPU. This approach lacks accuracy as it is hard to perform such operations on a lot of nodes at the required time to show precisely the benefit of a given load-balancing strategy. Furthermore, CPU scaling cannot set any frequency on any processor: only few frequencies are usually available, which can be an issue to subtly study a load-balancing strategy.

B. A framework for the evaluation of HPC applications 1) General overview of Distem: The features, general architecture and evaluation of Distem were presented in [START_REF] Sarzyniec | Design and Evaluation of a Virtual Experimental Environment for Distributed Systems[END_REF], but are shortly summarized here. Distem is able to emulate a heterogeneous virtual platform on top of a homogeneous cluster in order to meet the experimenter's requirements. A virtual platform is composed of virtual nodes and virtual networks.

On the networking side, Distem can emulate a complete and complex network topology. Virtual nodes can belong to different virtual networks and can act as gateways between several networks. Network links can be configured in order to provide specific performance (latency, available bandwidth). On the computing side, Distem can use two methods to emulate slower nodes in the virtual platform: a) CPU scaling and b) a specific method designed to work independently from the native supported frequencies of a CPU.

Distem also proposes two operating modes. In the realscale mode, the total number of emulated resources will be the same as the number of physical resources on the platform. In the node-folding mode, Distem uses lightweight virtualization to launch several virtual nodes (vnodes in Distem's terminology) on the same physical resource, that will be shared among them. Lightweight virtualization is based on Linux containers which is a strong advantage since applications are executed on the same hardware as classical execution. This is not really the case when running applications on the top of classical virtualization solutions.

2) Emulate heterogeneous and evolving experimental conditions: We have previously shown [START_REF] Sarzyniec | Design and Evaluation of a Virtual Experimental Environment for Distributed Systems[END_REF] that Distem can be used to create virtual heterogeneous platforms from a homogeneous cluster, with vnodes with different frequencies and number of cores, different IO capabilities and network links with different bandwidth and latency. In the context of this work, we extended Distem so that those characteristics of the vnodes could be updated dynamically without redeploying the virtual platform. This is useful to achieve complex experiments where the platform is modified, like it could happen in the reality. Thus, evaluating CPU-or communication-aware loadbalancers of HPC runtimes is possible and easy.

3) Support for fault-tolerance: Handling failures is crucial for HPC runtimes, thus being able to generate realistic failures is required to perform accurate experiments. Unexpected loss of nodes in a computing infrastructure is normally due to a hardware failure. Nodes can be lost in three different ways:

• Graceful: the node is shut down cleanly, using an operating system command. All applications are notified and they have some graceful period to perform some cleaning tasks (usually some seconds). • Soft: The node is forced to shut down. Applications running on the node as in the previous case are notified, however the node is shut down immediately. The application do not have the possibility to react. • Hard: The node failed abruptly. The operating system is not able to terminate cleanly. The peer application are not notified (network communications are still seen as active until the failure is detected by the peer, using a timeout).

Distem is able to emulate these three behaviors using Linux containers and Netfilter. More specifically, the failures are implemented as follows:

• Graceful: classical stop of LXC containers (lxc-stop), this sends a TERM signal to all the processes belonging to the container. If applications are able to handle specifically this signal, they can perform some cleaning operations.

• Soft: forced stop a LXC containers (lxc-stop -k),
this sends a KILL signal to all the processes belonging to the container. The effect is immediate. • Hard: we drop all the packets received and sent by a vnode using iptables.

Speaking in the fault-tolerance terminology, Graceful and Soft failures are fail-stop failures, Hard failures are a kind of omission failures. Performance failures can be emulated as well by setting a very high latency on some virtual network interface or by radically decreasing the CPU performance of some nodes, implying very slow responses to requests. Distem is not able yet to generate byzantine failures but it could be added, if you consider the network point of view, leveraging Netfilter. Indeed, libnetfilter_queue offers an API to alter packets received by the kernel and allows to re-inject them.

4) Event injection framework:

Introducing a platform modification manually is convenient to test a given feature but does not allow complex and realistic experimentation. Thus, being able to modify the virtual platform automatically and in a deterministic way during an experiment is required. For that purpose, we added an event injection framework.

The following automatic modifications are supported for a given set of vnodes: (1) modification of the CPU frequency; (2) modification of the network interfaces capabilities (latency and bandwidth); (3) start and stop (in the three different presented ways). The events can be specified in two ways. First, it is possible to use an event trace that specifies which modification occurs at which date (relatively to the start of the experiment). Second, it is possible to define automatically the date of event arrival according to various probability distributions. Currently, uniform, exponential and Weibull distributions are supported. The virtual platform description (CPU frequency modification, network link modification, and virtual node arrival/disappearance) can be updated dynamically and upon request during an experiment.

Thanks to this framework, the orchestration of complex and realistic platform modifications can be performed in order to achieve advanced experiments in a great variety of contexts.

III. EXPERIMENTAL EVALUATION

This section validates the benefits of Distem in the context of HPC runtimes evaluation. For that purpose we performed four kinds of experiments. First, we show how three HPC runtimes detect failures and behave when a process of the application is hit by a failure (Section III-A). Second, we validate the fault injection mechanisms of Distem depending on the degree of over-subscription (Section III-B). Third, we evaluate the fault-tolerance strategies implemented in Charm++ using a predetermined failure trace and probabilistic failures following several MTB values (Section III-C). Finally, we evaluate the behavior of some Charm++ load-balancers using Distem to create heterogeneous platforms and to dynamically change performance of nodes (Section III-D).

Experimental evaluation often suffers from poor repeatability [START_REF] Gómez | Replications Types in Experimental Disciplines[END_REF], which is a major problem to trust the results. Using Distem alleviates this issue as emulated platforms can be reproduced identically between two experiments, even when performed on different physical platforms. For the sake of reproducibility we also provide our experiments' output logs and scripts to generate the figures. These are available at 1 .

Experimental Setup:

The experiments have been performed on three clusters of the Grid'5000 testbed2 :

• Graphene, in Nancy, 144 nodes ( We used the following software stack: Debian Jessie, Linux kernel version 3.16, Charm++ version 6.6.0, OpenMPI 1.8.5 and MPICH 3.1. We used two applications for Charm++: a) Jacobi3D a 7-point 3-dimensional stencil that computes the transmission of heat on a three dimensional space. It implements in-memory checkpoint-restart [START_REF] Zheng | FTC-Charm++: An In-Memory Checkpoint-Based Fault Tolerant Runtime for Charm++ and MPI[END_REF] and b) Stencil3D a 6-point 3-dimensional stencil that iteratively average values in a 3D grid. This application has been chosen since it is relevant to observe the benefit of Charm's LB functionality. For the MPI based runtimes, we used the set of NAS benchmarks. In this section, we use Distem to assess the failure detection capacity of HPC runtimes. The experiment consist in injecting a failure with Distem and observe if the failure is detected and what is the action taken by the runtime. We use Jacobi3D application for Charm++ and the NAS benchmarks for the MPI based runtimes. The parameters and specific application are not mentioned as the results remained unchanged for a large set of applications and parameters.

A. Evaluating failure detection of HPC runtimes

The results are summarized in Table I. Charm++ detects Graceful and Soft failures. Charm++ reacts in the same way to both kind of failures: it rolls back the execution to the previous checkpoint. In the first case the runtime catches the signal but it ignores it. When Hard failures are injected Charm++ is unable to detect them and the application freezes 3 . OpenMPI detects the three types of failures. However, the runtime does not abort the execution. For the first two types of failures the runtime clean up all the processes running on the different nodes but the main application blocks. In the presence of a Hard failure the runtimes does not clean the processes running on the other machines 4 . MPICH detects the three types of failures. It aborts the execution when a failure is detected which is the expected behavior.

Testing HPC runtimes under the presence of failures is very important and it has to be integrated into the software development cycle. Distem simplified the uncovering of problems in the failure handling for widely used HPC runtimes.

B. Validity of fault injection mechanism

In this experiment, we validate the Distem fault injection mechanism. To do so, we compare three mechanisms of fault injection: a) the mechanism built into Charm++, b) real faults using a classical OS reboot (/sbin/reboot), c) Distem fault injection mechanism when running 1,2, or 4 containers per physical machine.

We run 1000 iterations of Jacobi3D using 64 physical machines in the Graphene cluster. The checkpoint frequency was set up to each 100 iterations. A failure trace was defined, that injects faults after the first, second and third checkpoint. Faults here, refer to Soft shutdown. We inject a total of 10 failures that were distributed in a way that the probability of making replica nodes crash is reduced. The same trace of failures is used for the three mechanisms. In Table II we observe that Charm++ fault injection provides a very optimistic and unrealistic scenario where the application always finishes its execution successfully. We include in the table the walltime of the application in order to show the overhead introduced by Distem which goes from 6.17% to 9.25% due to the use of containers and the degree of over-subscription. Additionally, it is shown that the walltime remains in an acceptable range to be meaningful and to be comparable with real injection. In this experiment, we show that fault injection mechanisms at the application level have to be complemented with a more realistic approach.

We demonstrate that Distem offers realistic experimental conditions and that the level of realism is maintained in the presence of over-subscription.

C. Study of fault tolerance strategies in Charm++

In this experiment, we use Distem to evaluate the fault tolerance capacity of the Charm++ runtime [START_REF] Kale | CHARM++: A Portable Concurrent Object Oriented System Based on C++[END_REF]. The experiment has been divided in two parts. The first part uses a static and predetermined trace of failures like described in the previous section whereas the second part uses events that are generated following a probabilistic law. We use the event generator described in Section II-B2 which is able to manage these two aforementioned scenarios.

1) Trace of predetermined failures:

For this experiment, we use the Graphene cluster. We run 1000 iterations of Jacobi3D on 128 cores (32 nodes, 4 vnodes per node, 1 core per vnode) with a problem size of 1024×1024×1024 and a decomposition of 128 3 and configured with different frequencies for the checkpoint: {10, 25, 50, 100, 120, 130, 140} iterations. This experiment evaluates the cost and trade-off of different checkpoint frequencies. This is summarized in the resulting application walltime (Figure 1). The most important overheads are at both extremes of the range: 10 and 140 which represent the overhead of performing a checkpoint very often and the overhead of loosing a large number of iterations. The overhead due to the loss of iterations is around 20% larger than the loss due to checkpointing. This experiment illustrates how Distem can be used to fully automate the tunning of fault tolerant applications by choosing the best checkpointing frequency. This is necessary given that the cost of checkpointing will vary according to the application.

2) Trace of failures that follow MTBF: For this experiment, we use the Griffon cluster. We run 1000 iterations of Jacobi3D on 256 cores (32 nodes, 8 vnodes per node, 1 core per vnode) with a problem size of 2048 × 2048 × 1024 and a decomposition of 256 3 . We set the frequency of checkpointing to 200 iterations.

In this experiment instead of using a deterministic trace we use a probabilistic distribution: the Weibull distribution. This distribution is suitable for approximating MTBF behaviors [START_REF] Robert | Models for fault-tolerance at very large scale[END_REF], [START_REF] Koren | Fault-tolerant systems[END_REF]. We use decreasing Weibull's scale parameter from 1 year (one of the classical MTBF values for one processor) to 1 day (to simulate faulty machines) and we use this event management feature on every vnode to ensure that they can fail independently and following the selected MTBF. We also choose a small shape parameter k = 0.5 (typical values being 0.5 and 0.7) to ensure that failures will arrive rapidly. With this, we increase the probability of encountering a failure in our experiments without changing the scale.

Figure 2 presents the results of the study of Charm's fault tolerance (FT) feature with different MTBF values. One can observe that the overhead of replication is quite low compared to the gain in case of failure (as already shown in the previous experiment). The replication only induces a slowdown of a few seconds on a 20 minutes run. It is interesting to observe that a single failure (in the case of a 365 days MTBF) leads to a noticeable loss of time. When comparing with Table III, which presents the number of iterations lost at a given failure 5 , we see that for a 180 days MTBF, the total number of iterations lost is larger than for a 90 days MTBF but with a total run time being smaller. This is due to the fact that failure detection and restart preparation took between 10 and 20 seconds instead of the 3 seconds measured in the other cases.

MTBF

In this experiment, thanks to Distem's failure generation framework, we were able to easily simulate node failures with different MTBF values and thus to observe Charm's FT efficiency. Using a probabilistic law to generate events is suitable to achieve realistic evaluations. However, depending on the time and the scale of an experiment, results might not be fully reproducible because of the probabilistic aspect.

D. Evaluating load balancing strategies in Charm++

In this experiment we use Distem to generate heterogeneous conditions from a homogeneous platform and to generate dynamic perturbation on the nodes. With this setup, we illustrate Charm's load balancing (LB) behavior in heterogeneous conditions. We used 8 physical machines of the Paravance cluster. For each machine a vnode is configured per core, making a total of 128 vnodes. We run 100 iterations of Stencil3D application with a 1024 3 problem size and a 128 3 decomposition. We create two different scenarios:

• heterogeneous: A platform where half of the vnodes have a CPU clock reduced by Distem from 2.4GHz to 1.2GHz. This case is an interesting study of what would be the result of running Stencil3D on 2 cluster partitions with different CPU clock speeds. • dynamic: A homogeneous platform where the available CPU power of a sub-part of the vnodes is dynamic. Each two minutes, 1/8 of vnodes are downclocked by Distem for two minutes, then set back to their original frequency. This scenario can mimic a cooling failure in a data center or CPU performance variation due to oversubscription [START_REF] Gupta | Improving HPC Application Performance in Cloud through Dynamic Load Balancing[END_REF].

For both scenarios we evaluate two load balancers (Re-fineLB and Hybrid) provided by Charm++. Therefore we test the following variations: LBOff: load balancing is deactivated. Charm++ objects will not be migrated and no processor imbalance check will be done by Charm++ runtime. RefineLB: centralized load balancer that moves objects away from the most overloaded processors to balance the workload. It also limits the number of objects being migrated. Hybrid: distributed load balancer that uses a hierarchical strategy.

Figure 3 shows the CPU utilization for each of the 128 vnodes in the two scenarios and the different variations. Figure 3a shows the CPU utilization for an heterogeneous hardware with not load balancer activated (LBOff ). It is observed that half of the vnodes have a lower CPU utilization which is the result of running a tightly coupled application in an heterogeneous platform. Figures 3b,3c depict a smooth curve for the CPU utilization as a result of the activation of the load balancers. Although Hybrid shows a better CPU utilization, the final runtime of the application is bigger than the normal execution time (LBOff ). This is due to CPU utilization generated by the load balancer itself. RefineLB shows a CPU utilization of 5.67% and a reduction of 6.40% in the walltime. Figure 3d shows the CPU utilization for a dynamic platform. Again the curve of utilization is smoothed by both load balancers. We observe here, as in the previous case that RefineLB behaves better than Hybrid. In this experiment, the vnode 0 is downclocked periodically. For the test without LB, it is visible that vnode 0 is globally slower than others and, it delays vnodes 1, 2 and 3 that show a lower CPU utilization than the others. This experiment demonstrated that Distem enables experimenters to easily simulate perturbations and heterogeneity of nodes in order to evaluate load balancing strategies.

IV. RELATED WORKS

There exists a large number of works that aim at building experimental environments providing controlled performance of CPU or network. These are more extensively detailed in our previous works [START_REF] Buchert | Methods for Emulation of Multi-core CPU Performance[END_REF] and [START_REF] Sarzyniec | Design and Evaluation of a Virtual Experimental Environment for Distributed Systems[END_REF]. However, none of them provide the complete integration of the chain of tools for the full control of all the environment characteristics as Distem does.

In the context of HPC runtime evaluation, no standard method is used to run experiments. As the Charm++ team is seen by most as doing state-of-the-art research in the field of HPC runtimes, we use them as an example and survey some of their recently published papers to see how their experiments have been carried out. [START_REF] Menon | A Distributed Dynamic Load Balancer for Iterative Applications[END_REF] proposes to evaluate a distributed load balancer for Charm++ using an ad-hoc simulator and a real application trace from IBM BG/Q. Building one's own simulator and extracting trace from real application is time consuming. Simulation is valuable to deal with extreme-scale platform, but it might not represent all the details of real applications. [START_REF]A 'Cool' Way of Improving the Reliability of HPC Machines[END_REF] describes a temperature-aware load balancer for Charm++. Their experiments leverage the DVFS capabilities of processors to study the load balancer behavior. DVFS is not always available or is available only for few frequencies; the same study could have been carried out with Distem without specific hardware support. Indeed, Distem is able to emulate DVFS at core level. A similar experiment using DVFS is presented in [START_REF] Menon | Thermal aware automated load balancing for HPC applications[END_REF]. However, it could not be performed at the core level since DVFS was only available at the socket level. Such a limitation does not exist on a Distem platform. [START_REF] Ni | ACR: Automatic Checkpoint/Restart for Soft and Hard Error Protection[END_REF] and [START_REF] Ni | Hiding Checkpoint Overhead in HPC Applications with a Semi-Blocking Algorithm[END_REF] present fault tolerance mechanisms on top of Charm++. Studying such properties of the HPC runtimes is probably one of the most difficult things. Generating realistic failures is not easy: one possibility is to mimic a failure in the application; another is to kill the application's processes, or to shutdown some nodes. Distem removes the burden of managing this manually by providing users with an event injection framework. This framework allows the experimenter to define realistic failure conditions, like MTBF, frequently used as a failure model in the literature. [START_REF] Totoni | Toward Runtime Power Management of Exascale Networks by On/Off Control of Links[END_REF] proposes to study network link utilization of different applications according to the underlying network topology. To create new network topologies, this work required advanced network management capabilities to switch off some links. Such experiment could be performed on any platform with Distem since virtual network topologies can be redefined.

V. CONCLUSIONS AND FUTURE WORKS

In this paper, we presented how Distem was extended to support the evaluation of HPC runtimes and applications. After detailing some relevant design parts of Distem with respect to the HPC context (heterogeneous and evolving experimental conditions, realistic fault-injections), we carried out several experiments, studying two key properties of HPC runtimes that are included on the targeted challenges to reach before 2020: (1) fault-tolerance and (2) load-balancing properties of HPC runtimes. A lot of other experiments could have been carried out but the goal was to prove the usefulness of Distem with respect to enabling those experiments. In each experiment, we were able to draw conclusions that might be useful for the designers of runtimes. For instance, new load balancing strategies could be implemented in Charm++ in order to fit the next-generation HPC platform characteristics. Being able to execute experiments on a large set of platform configurations in a repeatable way, as made possible by Distem, is a sound basis to design and improve the HPC runtimes of the future.

Future works will follow three directions. First, we will carry out experiments related to I/O evaluation. As seen in the paper, it would be an interesting feature to evaluate mechanisms using on-disk checkpointing. Distem already has support for I/O throttling that leverages a Linux feature called blkio but extensive evaluations still have to be achieved and maybe some design refactoring will have to be performed since in some cases the throttling might not work as expected (case of buffered I/O for instance). Second, we would like to emulate more node characteristics, such as power consumption, heat dissipation in order to perform realistic experiment that focus on future HPC and Cloud challenges like power-aware scheduling. We could also imagine introducing an energy cost model to be able to evaluate fancy scheduling/loadbalancing algorithms like follow-the-sun algorithms. Finally, we look forward to working with HPC runtime and application developers in order to contribute to designing a better next generation of HPC software.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Overhead of different checkpoint frequencies under the presence of failures

Fig. 3 :

 3 Fig. 3: Figures (a), (b) and (c) represent the execution in a virtual infrastructure where half of the nodes CPU clock is reduced to half of the frequency. Figures(d), (e) and (f) represent the execution in a virtual infrastructure where the CPU clock of some machines is varying.

TABLE I :

 I Failure detection. C refers to the roll-back of the application to the previous checkpoint, H refers to the fact that processes hang, E refers to the termination of MPI processes

TABLE II :

 II Percentage of successful application executions

TABLE III :

 III Number of iterations lost in failures.

https://gforge.inria.fr/frs/?group id=3519

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).

Charm++ developer team was notified about this behavior

This problem has been fixed in the latest OpenMPI version 1.10.1

A loss of 0 iteration is due to 2 contiguous failures, in this case, the 2nd failure does not impact the number of iterations lost.