domain: spi phys

text:
Guilhem Rival
email: guilhem.rival@onera.fr

Eric Dantras
email: eric.dantras@univ-tlse3.fr

Thierry Paulmier
email: thierry.paulmier@onera.fr

Electronic irradiation effects on the electrical behaviour of a highperformance space-used thermoplastic polymer

Introduction

In space environment, satellite materials are subjected to several environmental constraints like high vacuum, temperature cycling or irradiation by space particles (electrons, protons, …). In case of polymer materials, irradiation by electrons presents two major issues. Since they are electrical insulators, the accumulation of electrons on their surface will lead to electrostatic discharge phenomena that can be responsible of electrical failures. On the other hand, high-energy electrons go through the polymers and transfer part of their energy to the material, which induces different phenomena like Radiation Induced Conductivity or physicochemical ageing. The initial physical properties of the polymer can then be modified. In this context, numerous studies have investigated these phenomena for many space-used polymers in order to anticipate their behaviour once in orbit [START_REF] Hanna | Radiation induced conductivity in space dielectric materials[END_REF], [START_REF] Roggero | Electrical conductivity of a silicone network upon electron irradiation: influence of formulation[END_REF]. This study focuses on PolyEtherEtherKetone (PEEK), a high-performance thermoplastic polymer which has recently but increasingly been used for space applications. Thus, it is necessary to analyse the effects of electronic irradiation on its physico-chemical structure and on its physical properties.

Results and Discussion

100 µm thick PEEK samples were irradiated, under high vacuum and at room temperature, in the SI-RENE irradiation facility installed at ONERA, Toulouse. Two irradiations campaigns have been carried out to inject ionizing doses of 12 MGy and 34 MGy in samples. Pristine and irradiated PEEK samples have then been analysed by Broadband Dielectric Spectroscopy. The DC conductivity data extracted for each sample at low frequency and high temperature, according to the Jonscher's law [START_REF] Jonscher | The "universal" dielectric response[END_REF], have been reported in an Arrhenius plot (see Figure ). With an increasing ionizing dose, DC conductivity of samples decreases indicating of a limited charge transport after irradiation. By using other analysis technics like Differential Scanning Calorimetry (DSC), this decrease has been associated with a crosslinking phenomenon observed through the increase of glass transition temperatures. Thus, the crosslinking nodes act as new traps for electric charges. Moreover, DSC also showed a degradation of the crystalline phase and therefore, it can be supposed that defects created at the amorphous / crystal interfaces decrease conductivity as well.

Conclusions

Irradiated PEEK samples showed a decrease of their DC conductivity associated with new charge traps created during irradiation. Thanks to different analysis, these new traps have been identified as crosslinking nodes. Thus, charge leakage being limited, this will lead to higher electrical potential on the polymer and therefore, to an increase of the risk of electrostatic discharges initiation.

1 )Figure :

 1: Figure: Arrhenius plot of DC conductivity of pristine and irradiated PEEK samples

Acknowledgement

The authors would like to thank the Région Occitanie for their financial support in this project.


Bo Sundman
email: bo.sundman@gmail.com

Christophe Sigli

Matthias Stratmann

Romain Le Tellier

Ursula R Kattner

Mauro Palumbo

Suzana G Fries

Keywords: Computational Thermodynamics, Phase Transformations, Simulations, Parallel Computing

Thermodynamic data are needed for all kind of simulations of materials processes. For kinetic models, for example, thermodynamic quantities, such as chemical potentials, driving forces for precipitation, thermodynamic factors for converting mobilities to diffusion coefficients, solubilities in different phases etc., are needed. It is very important that these data are consistent and accurate, also for predicting metastable states during phase transformations. Equally important is that the computation time must be reasonably short.

Frequently thermodynamic data are pre-evaluated into "lookup tables" to speed up calculations. This creates additional uncertainties in the data as they must be interpolated or extrapolated and conditions may differ from those originally assumed for the creating the lookup table. It would be much better if a full local equilibrium calculation could be performed when

Introduction

The background for the development of a free thermodynamic software for multicomponent calculations using the CALPHAD method (CALculation of PHAse Diagrams) is described by Sundman et al. There is extensive literature describing the use of thermodynamic data for materials and process development. Although somewhat outdated, an excellent review can be found in the book by Miodownik and Saunders [3].

Especially the development of phase field techniques [4,5,[START_REF] Warnken | Phase-field modeling of as-cast microstructure evolution in nickel-based superalloys[END_REF] and other techniques, such as fluid dynamics [START_REF] Schneider | Modeling of micro-and macrosegregation and freckle formation in single-crystal nickel-base superalloy directional solidification[END_REF] and finite element methods [START_REF] Banerjee | Coupling of phase diagram calculations for muliticomponent alloys with solidification micromodels is casting simulation software[END_REF] for simulating microstructure evolution has increased the demand for fast and reliable calculations of chemical potentials, driving forces, mobilities etc. Within the framework of Integrated Computational Materials Engineering (ICME) [START_REF]Integrated computational materials engineering: a transformational discipline for improved competitiveness and national security[END_REF] the necessary tools are coupled to accomplish the desired simulations and it is the intention that the OC software becomes one of these tools.

The OC software has been extensively tested on Windows and Linux systems and although there are limitations and occasional problems the overall performance is comparable with commercial thermodynamics software for materials science, such as FactSage [START_REF] Bale | [END_REF], Pandat [11] or Thermo-Calc [12].

In addition, the OC software can perform many equilibrium calculations in parallel, including calls to OCASI. The OCASI software is still under development but has been used in some applications to demonstrate how it can be integrated in simulations to utilize accurate and consistent ther- A development version with some new features and bugfixes is available at [http://www.github.com].

Computational Thermodynamics

Computational thermodynamics consists of software and databases for calculating phase diagrams and other kinds of thermodynamic data. The databases are a collection of parameterized model descriptions of the phases as function of constitution and state variables. The model descriptions are obtained from the thermodynamic assessment (or optimization) of individual systems. The basic techniques are described in the book by Lukas et al. [START_REF] Lukas | Computational Thermodynamics, the Calphad method[END_REF].

An important advantage from using an assessed thermodynamic database is that all data calculated with such a database are consistent, including phase solubilities, melting temperatures, heat capacities, chemical potentials etc.

State variables

A thermodynamic system is described by a set of state variables like T, P, V for temperature, pressure and volume. The amounts of the components can be given in moles, N, mole fractions, x, or mass fractions, w and the equilibrium state is at a minimum of the appropriate state function like internal energy, U, Helmholtz energy, A or the Gibbs energy, G.

Many properties can be derived from the state functions and some can also be used to control the system, such as chemical potentials, µ, enthalpies, H etc. All state variables and properties available in OC are listed in Appendix 2.

Thermodynamic models

In CALPHAD-type thermodynamic databases the Gibbs energy of each phase is described by a model reflecting its structure. Several different kinds of models are described by Lukas et al. [START_REF] Lukas | Computational Thermodynamics, the Calphad method[END_REF]. The simplest model is that for an ideal gas:

G gas M = i y i • G i + RT y i ln(y i ) + RT ln(P ) (1)
where • G i is the Gibbs energy of formation of the gas molecule i from the reference state of the elements, y i is the constituent fraction of the molecule in the gas and P the pressure. R is the gas constant, T the absolute temperature and y i ln(y i ) is the random configurational entropy.

The Gibbs energy model is per mole formula unit of the gas as the number of atoms will vary with the species. The mole fractions of the elements, x j , are given by

x j = i b ji y i k i b ki y j (2)
where j and k denote elements and i the molecules. b ji is the stoichiometric factor of element j in molecule i.

More complex phase models with sublattices for interstitials, chemical ordering or defects, including ionic constituents are implemented in OC. The interactions between the constituents on a sublattices can be described by excess parameters and a model to describe the contribution of magnetic ordering depending on the Curie or Néel temperature and the Bohr magneton number to the heat capacity is available. All these models are special cases of Compound Energy Formalism (CEF) [START_REF] Hillert | [END_REF] and are implemented in OC as well as the partially ionic 2-sublattice liquid model [17] for short range ordering (SRO) in the liquid phase. More models will be implemented with time and since the source code of OC is available skilled users can implement their own models.

Equilibrium calculations

The equilibrium algorithm used in OC was proposed by Hillert [18] and implemented by Jansson [19] in the POLY module of Thermo-Calc software and in the PMLFKT software by Lukas [20]. The implementation of this algorithm in OC is described in [21]. In this context it is sufficient to mention that the data structure, described in [21], is designed to allow calculation of several equilibria in parallel and it has an additional grid minimizer to ensure that the calculated equilibrium is a global and not a local one. It can also detect miscibility gaps.

A special feature of Hillerts algorithm [18] is the possibility to prescribe that a phase should be stable. This feature is essential for the ability to calculate phase diagrams and the direct calculation of a transition, such as the melting temperature of an alloy.

Parallel calculations

Calculation of multiple equilibria in parallel is not only useful for simulations but also for assessments of model parameters. Such assessments usually include a large number of experimental data, each representing a value determined at equilibrium. The assessment procedure calculates these equilibria many times varying the selected model parameters in order to find the best fit. Traditionally the equilibria are calculated sequentially but with the OC software it is possible to calculate them in parallel.

When tested on a PC with 4 CPU, each with 2 threads, the gain in speed for a calculation with 400 experimental equilibria in parallel was almost a factor of 4, i.e. no significant overhead.

The occurrence of a miscibility gap in a phase requires special treatment for parallel processing. A miscibility gap means that the phase is stable with two different compositions which requires that an extra composition set must be created in OC, see section 3.7. This must be done at the same time in all equilibria otherwise the data structure would become inconsistent. Thus the threads must be stopped while the extra composition set is created.

Normally this is not a problem during an assessments as miscibility gaps should be known from the experimental data and the user can manually create a sufficient number of composition sets before starting the parallel session. However, miscibility gaps or order/disorder transformations can create problems when encountered during simulations.

Additional data depending on T, P and phase composition

Thermodynamic data are important but not sufficient to simulate a phase transformation. For example, to consider the diffusion of elements their mobilities together with the thermodynamic factor need to be described to obtain the diffusion coefficient. In addition, there are also interface energies, interface mobilities, elastic constants and other phase based properties that may be of interest [22].

In the OC software such property descriptions can be stored together with the thermodynamic data for each phase as functions of T, P and constitution of the phase. The values of these properties are calculated together with the thermodynamic properties during the equilibrium calculation and can be accessed by the application software.

For example, the natural logarithm of the mobility of Fe in the face centered cubic (FCC) phase of the binary Fe-Ni system, in OC denoted MQ&Fe(FCC), can be described by a linear function:

MQ&Fe(FCC) = x FCC Fe MQ&Fe(FCC, Fe) + x FCC Ni MQ&Fe(FCC, Ni) (3)
where x FCC

MQ&Fe(FCC, Ni) = Q Fe RT ln(ν Fe ) (4)
where Q Fe is the activation energy and ν Fe the pre-exponential factor. As mobilities can vary several orders of magnitude it is a better approximation to describe the composition dependence by a linear relation of its logarithm.

If sufficient experimental or theoretical data are available, non-linear terms can be added to eq. 3.

Together with the thermodynamic factor the mobility can be transformed to a diffusion coefficient needed for the simulation of microstructure evolution in diffusion controlled processes.

The OCASI application software interface

Software interfaces for coupling thermodynamic calculations with application software following a proposed standard [START_REF] Eriksson | 2nd Colloquium on Process Simulations[END_REF] implemented in slightly different ways in ChemApp [START_REF] Petersen | [END_REF] and Thermo-Calc TQ [START_REF]Thermo-Calc Software AB[END_REF]. The idea behind this proposed standard is to overcome the differences between the different thermodynamic software packages which can make the implementation into application software difficult. The methodology with many separate subroutines and functions performing specific tasks makes the implementation straightforward and used also for OCASI.

A full documentation of OCASI is under preparation but a preliminary version is available at the websites mentioned above. Again, as the source code is available any user can directly access the source code and also data which are not primarily available. However, any such access should be done with caution as the data structure and source code may change in a later update whereas the OCASI interface should remain stable. A summary of the available subroutines is given in Appendix 1.

Fortran and C++ data structure compatibility

Using the ISO C standard it is possible to access data structures defined in the OC software written in Fortran from C++ and other languages. For example, the data for a composition set of a phase in OC are described by:

Interfacing with Java

Based on Fortran and C(++) compatibility, it is also possible to interface OC functionalities in the Java language thanks to the Java Native Interface (JNI) [START_REF]Java Native Interface[END_REF]. For the PROCOR software [START_REF] Tellier | Phenomenological analyses of corium propagation in LWRs: the PROCOR software platform[END_REF], such an interfacing has been developed for C. Cardon Ph.D. thesis regarding the development of a Cahn-Hilliard based multicomponent multiphase diffusion model [START_REF] Cardon | [END_REF]. This application is briefly described in section 4.4. In addition to the Cahn-Hilliard model requirements, equilibrium calculation functionalities have been interfaced for the coupling of other physical models to evaluate some physical parameters using the phase compositions obtained from OC equilibrium calculations.

Reading the database

The thermodynamic data for a system will normally be read from a database. The OC software supports the TDB format [START_REF]Thermo-Calc Software AB[END_REF] proposed originally by SGTE [START_REF] Ansara | The Scientific Group Thermodata Europe (SGTE), Computer Handling and Dissemination of Data[END_REF] and is used by several other software, most notably by Thermo-Calc and DICTRA [12]. The application software allows to select the elements for which data that should be extracted from the database and then to suspend phases that are not of interest for the application.

The application software must initiate the OC memory structure by a special subroutine call. This returns a pointer to an equilibrium record in the dynamic memory structure which contains information about the constitution and calculated results for all the phases as described in section 3.1.

Only after this initialization other subroutines, such as reading the database, can be called. The database has data like the models of phases and their thermodynamic model parameters and these data are stored in the static part of the OC memory structure since they are independent of the external conditions. The external conditions, the constitution of the phases and the calculated results are in the dynamic data memory. A schematic representation of both static and dynamic memory is shown if Fig. 1.

Setting conditions

The conditions are set individually by a call to a subroutine. In the simplest case these represent the temperature, T , pressure, P , and amount of the components, N i ,. The number of conditions needed to calculate the equilibrium is given by Gibbs phase rule where f is the degrees of freedom, n the number of components and p is the number of stable phases. Normally the number of stable phases is the result of the equilibrium calculation but OC allows that one or more phases are prescribed to be stable. If there are no prescribed stable phases, a system with a single component, n = 1, must have 3 conditions, for example T, P

f = n + 2 -p (5)
and N where N is the size of the system. For a binary system 4 conditions must be set and so on.

The amount of a component i can be given in moles, N i , or mass and if there is a major component they are often given in mole or mass fractions or per cent and the total size of the system is given as N. Chemical potentials or activities of the components can also be used instead of amounts.. It is also possible to specify the volume or the enthalpy of the system, or as already mentioned, specify that a phase has to be stable. It is important that all conditions are independent of one another. However, this does not mean that it is not possible to set irrational conditions, i.e. those that do not represent a possible equilibrium state. A full list of state variables in OC is given in Appendix 2.

Creating additional equilibria

Additional equilibrium records can be created to represent local equilibria in a simulation or assessment. These may be calculated sequentially or in parallel as they have independent sets of conditions and related data. Each equilibrium record has a pointer to its own dynamic data structure and this pointer must be supplied in the calls to most of the OCASI routines to indicate which equilibrium should be used.

Calculations

When the degree of freedom according to eq. 5 is zero the equilibrium can be calculated by calling a subroutine. The calculation procedure in OC is described in detail in [21]. In the call it is possible to specify whether the global minimizer should be used or the calculation should start from the last calculated set of stable phases and constitutions.

Composition sets and phase tuples

The number of phases in the system when reading from the database is equal to the number of phases with unique names and structures in the database. However, as phases can have miscibility gaps or several types of ordering transformations some of the phases may become stable with two or more unique compositions. In some software the user must specify this in advance but in OC the global grid-minimizer will detect miscibility gaps automatically. A phase may thus have several independent composition sets and each of these is called a phase tuple. The composition set index of a phase is specified after the phase name by the character "#" followed by a digit, for example LIQUID#2.

In some cases, for example phases with order/disorder transitions, it is important to know the state of ordering of a particular ohase tuple. For such cases composition sets can be created explicitly and they can be given a unique name and default composition that can be used to identify the ordered state.

The first composition set created when the phase is read from the database has no explicit set indication. A composition set created automatically by the grid-minimizer will have a suffix " AUTO" added to its name (before the # sign).

In OCASI there is an array of phase tuples where each index represents the combination of a phase and composition set. Thus the number of phase tuples may change whenever the grid-minimizer is called.

Extracting results

The normal way to obtain calculated results is to call a subroutine using state variable symbols from the list in Appendix 2. The values of properties that are modeled separately, like the Curie temperature, mobilities etc. and that have been read from the database can be obtained in the same way.

Using the pointer to the equilibrium data structure provided when reading the database it possible to extract individual calculated results of an equilibrium. For example, if the equilibrium record pointer is called "ceq" and as the array of records for each phase is called "phase varres", the value of the Gibbs energy per mole formula unit of the phase after an equilibrium calculation can be obtained as

gfu = ceq%phase_varres(lokcs)%gval(1,1)
where the % sign is used in the new Fortran standard to indicate items in a record. The lokcs variable is the index of the composition set of the phase and this is stored in the phase tuple record and thus provides an easy link between the phase tuple index and the phase data. For example, if the calculated mobility of component 1 of the same phase is stored in the location "ipy", see section 3.1, its value is obtained by mob(1) = ceq%phase_varres(lokcs)%gval(1,ipy)

Great care must be taken that the relevant data is extracted. The application program must not change the OC data structure.

Examples of application software

Several simple examples of how to use the OC software interface are provided with the source code. Some of the examples presented here may require the use of proprietary software components.

The beginnings of a reactor module

A process can sometimes be simulated by dividing it into sections, each of which represent a local equilibrium, and a transport of matter between these by diffusion or other means. The equilibrium calculation can be done by OC in a function "localequil" with the following arguments: Here the todo argument specifies VERBOSE to generate extra output, then to open the database file cosi.tdb and read the data for the system with the elements in the argument selel and calculate the equilibrium at the specified values of T, P and amounts of the elements. After a successful calculation the amounts of all elements in the gas, n(gas,*), their chemical potentials, mu(*), the total enthalpy, h, and the enthalpy of the gas, h(gas), will be provided in the outval array argument. The state variables that can be used are listed in Appendix 2. If there is an error the function will return a non-zero number.

The localequil function will remember whether it has already read a database and if no database is specified in a subsequent call the routine will use the static data that it has already stored. In the application software the amounts of elements in the gas and the other phases may be redistributed between different local equilibria to simulate a process. The range of options in the localequil function can easily be extended. One of the intentions of this function is to provide an application programmer a simple access of the OCASI interface.

Simulation of homogenization of an alloy AA7449

The OpenCalphad Software has been recently coupled to the Constellium thermodynamic package, ProPhase, using OCASI in order to use the faster equilibrium routine included in OC as well as the new parallelization recently implemented. ProPhase includes a long range diffusion subroutine which is used when simulating a solidification path with back diffusion or a homogenizing treatment of an as cast structure. The long-range diffusion within the FCC aluminum solid solution is simulated using an explicit finite difference scheme, assuming local equilibrium in each volume element, calculating the solute concentration in solid solution, and using Fick's diffusion equation in the solid solution to compute solute exchange between adjacent volume elements. For each time step, this routine performs many independent equilibrium calculations, one for each volume element, and is very well suited for parallelization. The simulation of a laboratory heat treatment at 472 • C of a multicomponent aerospace 7449 alloy (an Al-Mg-Zn-Cu-Fe-Si alloy, see Table 1) after casting has been used to measure the gain in speed provided by OpenCalphad in comparison with the previous software. The cast was performed in the laboratory using the Aluminum Association TP-1 grain refiner test crucible. The corresponding cast structure, as shown in Fig. 2 was simulated by a modified Scheil calculation with back diffusion. For this particular case, a two dimensional cylindrical geometry was chosen with a 61 µm secondary dendrite arm spacing. Care must be taken to mesh the dendritic cylinder with sufficiently small enough radius increments especially towards the end of solidification where the concentration gradients may be very steep. The evolution of composition gradients has been studied for a heat treatment consisting in a linear heating of 40 • C/h up to 472 • C followed by holding at 472 • C for 10 h. The corresponding predicted and measured incipient melting temperatures and integrated enthalpies for the DSC first melting peak are compared in Fig. 3. Overall, good agreement is obtained indicating that the proposed physical description of the phenomena occurring during solidification and homogenizing is acceptable. With our previous software such a calculation took 3.5 days on one processor. By using the parallel OpenCalphad routine the total computation time was reduced to 3.5 hours on 12 processors.

Phase-field models and OCASI

Many modern thermodynamic consistent phase-field models [START_REF] Tiaden | The multiphase-field model with an integrated concept for modelling solute diffusion[END_REF][START_REF] Kim | Phase-field model for binary alloys[END_REF][START_REF] Steinbach | Phase-field model with finite interface dissipation[END_REF][START_REF] Zhang | Phase-field model with finite interface dissipation: Extension to multi-component multi-phase alloys[END_REF] rely on thermodynamic and kinetic properties that can be obtained or ap- To be able to compute the equations governing phase-transformation and diffusion in the phase-field model with finite interface dissipation [START_REF] Zhang | Phase-field model with finite interface dissipation: Extension to multi-component multi-phase alloys[END_REF], Gibbs energies of the phases, G α m , their derivatives ∂G α m /∂y α,(s) i

and the atomic mobilities M α i can be calculated with an input consisting of the temperature T , pressure P and compositions y α,(s) i through the OCASI interface. This does not only speed up the simulations and increases the accuracy of the results, it also allows for a quick simulation setup of multi-component multi-phase simulations, as no analytical implementation of the Gibbs energies or parallelization of phase equilibria calculations, has to be carried out. This is demonstrated here for a small microstructure simulation example with two austenitic steel samples with different alloy composition in contact.

The two samples containing iron, manganese, silicon and carbon are chosen to represent the experiments on uphill diffusion by Darken [START_REF] Darken | Diffusion of Carbon in Austenite with a Discontinuity in Composition[END_REF]. Apart from the chosen thermodynamic and kinetic databases [START_REF]MatCalc thermodynamic and kinetic databases under Open Database Licenses[END_REF] under the Open Database License, only the temperature, the time step and the grid spacing had to be specified as in Tab. 2

The results of the quaternary simulation are shown in Fig. 4 for an experiment with a duration of 10 days at constant 1323 K and are compared with the results of a DICTRA sharp-interface simulation as well as measurements Thermodynamic database MatCalc mc fe v2.057.tdb [START_REF]MatCalc thermodynamic and kinetic databases under Open Database Licenses[END_REF] Mobility database MatCalc mc fe v2.008.ddb [START_REF]MatCalc thermodynamic and kinetic databases under Open Database Licenses[END_REF] from Darke, [START_REF] Darken | Diffusion of Carbon in Austenite with a Discontinuity in Composition[END_REF]. The results obtained with only openly available tools and data show excellent agreement with the results from the commercial DIC-TRA software.

Beyond that, the coupling of OpenCalphad and phase-field software like OpenPhase has been proven to be an efficient, accurate and easy to use tool for meso-scopic microstructure simulations.

Mesoscopic modelling of a liquid miscibility gap

For the study of severe accidents in light water nuclear reactors (LWR), the OpenCalphad software has been used in the development of a Cahn-Hilliard based diffusion model for treating the liquid-liquid miscibility gap associated with the corium (oxidic and metallic liquid materials present during a reactor core meltdown) in a reactor vessel lower head.

In-vessel corium is a complex thermodynamic system because of the phases that are present. In addition to partially oxidized cladding and fuel materials (UO 2 , ZrO 2 and Zr), the corium contains stainless steel elements (Fe, Ni, Cr) from the reactor's internal structures. In particular, the U-O-Zr-Fe system exhibits a liquid-liquid miscibility gap responsible for a liquid phases stratification of tremendous importance on the corium propagation

= a ∂g ∂x i x j =i - j κ i,j ∇ 2 x j (6)
where a and κ i,j are parameters of the model related to the interface thickness (a "numerical" parameter in this mesoscopic modeling) and interface tension while g is the free energy density (in Jm -3 ) of the homogeneous system.

In the binary U-O case, the dg dx O term was calculated applying the chain rule to the liquid Gibbs energy G liq m and its derivatives ∂G liq m /∂y liq i obtained from the OCASI interface [START_REF] Cardon | [END_REF].

For the ternary and quaternary systems, the additional hypothesis that local equilibrium of the redox chemical reactions is instantaneously reached is made. Therefore, ∂g ∂x i x j =i can be related to the chemical potentials obtained from "local" equilibrium calculations performed using the OCASI interface where the phase separation associated with the miscibility gap is not taken into account (i.e. the grid minimizer is turned off, cf. Section 2.3).

In order to illustrate this modeling, a 1D calculation for the U-O binary case is presented here; the domain length is L = 10.0 cm and the temperature and pressure conditions are T = 3200 K and p = 1 bar. Fig. 5 schematically depicts the initial configuration and the steady-state of the system. The The initial configuration corresponds to an oxidic phase (uranium dioxide UO 2 ) above a metallic phase (in this case metallic uranium U). The diffusion process leads to a steady state configuration with two suboxidized phases.

Fig. 6 shows the profile of oxygen molar fraction x O in the initial state and at steady state where it is verified that, far from the interface, the oxygen molar fraction is equal to the values given by an equilibrium calculation.

  [1]. A very limited first version of the OpenCalphad (OC) software was released in 2013. In February 2015 a second version capable also of multicomponent single equilibrium and some phase diagram calculations was released. During 2015 extensive new facilities were added and the software has become much more stable and both version 2 and 3 are available on [http://www.opencalphad.org]. Version 3 includes parallel equilibrium calculations using the OpenMP [http://www.openmp.org] library, a rudimentary assessment module and an Application Software Interface, OCASI, callable from C++ as well as Fortran and other programming languages. OC is written in the new Fortran standard.

  modynamic data. Examples and full documentation of the source code as well as macro files showing the interactive use of the software and the use of OCASI are available on the OC web site [http://www.opencalphad.org].

Fe and x

 x FCC Ni are the mole fractions of Fe and Ni in FCC, respectively, MQ&Fe(FCC, Fe) is the mobility of Fe in pure FCC Fe and MQ&Fe(FCC, Ni) is the mobility of a single Fe atom in pure FCC Ni, both of which may depend on T as:

Figure 1 :

 1 Figure 1: The static (a) and dynamic (b) parts of the OC memory. Each equilibrium record in (b) is independent and can be calculated in parallel.

Figure 2 :Figure 3 :

 23 Figure 2: Simulated concentration gradients after solidification of aerospace alloy AA7449

Figure 4 :

 4 Figure 4: Simulated composition profiles after 10 days of tempering at 1050 • C, with experimental measurements by Darken [36] as circles.

[

  38]. This stratification phenomenon results from the combination of multicomponent multiphase species diffusion and gravitational material movement due to buoyancy effects initiated by Rayleigh-Taylor instabilities. The natural convection in a corium pool can be simulated by Computational Fluid Dynamics (CFD) and a Cahn-Hilliard based mesoscopic modeling of the species diffusion appears to be a promising way to account for the stratification phenomenon in detailed simulations. This R&D effort regarding a detailed modeling of the thermohydraulic/thermochemical behavior of invessel corium has been undertaken at CEA Cadarache in support of integral models development in the context of the PROCOR software platform for sensitivity/uncertainty studies regarding corium propagation.A Cahn-Hilliard based diffusion model was first developed for the U-O binary case[START_REF] Cardon | [END_REF] and is under development for U-O-Zr and U-O-Zr-Fe. Under the assumption of constant molar volume, the model for a n-component system is based on Cahn-Hilliard evolution equations for molar fractions of n-1 components with diffusion potentials, μi , expressed as: μi

  interface thickness has been set to 1 cm.

Figure 5 :

 5 Figure 5: U-O system initial configuration and steady state

Figure 6 :

 6 Figure 6: Profile of molar oxygen fraction x O versus position z

  The todo argument passes a number of instructions to localequil, nsel gives the number of elements specified in argument selel, tp provides the values of T and P and nmz the element amounts. The argument outspec specifies which state variable values will be returned in outval. A call of this function to calculate a local equilibrium for the system Si-O-C at 1800 K and 1 bar with a specified amount of moles of the elements can be:

	! todo character variable	with "commands" to be executed
	! nsel integer	number of components
	! selel character*2 array	element names
	! tp	double precison array(2)	with values of T and P
	! nmz double precision array(nsel) with molar content for elements
	! outspec character variable	with information which data to extract
	! outval double precision array(*) with extracted values
	! if there is no error the function value returned is zero
	integer function localequil(todo,nsel,selel,tp,nmz,outspec,outval)
	!
	! This rouine calculates an equilibrium and return values

nsel=3 selel(1)='C '; selel(2)='O '; selel(3)='Si' tp(1)=1.80D2; tp(2)=1.0D5; nmz(si)=1.0D0; nmz(2)=2.0D0; nmz(3)=1.8D0

todo='verbose open cosi.tdb ' outspec='n(gas,*) mu(*) h h(gas) ' ierr=localequil(todo,nsel,selel,tp,nmz,outspec,outval)

Table 1 :

 1 Composition range of AA7449 according to the Aluminum Association

	Si	Fe	Cu Mn Mg Zn Ti+Zr rest
	min 0	0	1.4 0	1.8 7.5 0	Al
	max 0.12 0.15 2.1 0.2 2.7 8.7 0.25	Al

Table 2 :

 2 Simulation parameters

	Box size	50x10x10
	Grid spacing	0.001 m
	Time step	1000 s
	Temperature T	1323 K

Appendix 2

The available state variables and their notation in OC. All of them can be used to extract calculated values and many also to set conditions.

In this notation N is the total amount of moles, N(O) is the total amount of oxygen in moles, N(GAS,O) is the amount of moles of O in the gas phase.

NM(O) or X(O) is the overall mole fraction of O. The total number of moles of components in the gas phase is NP (GAS).

The heat capacity is not included in this list since it can be calculated as a "dot derivative, H.T" , as explained in [21].


A Fabbri

H K K Wong

B Lai

T A Bui
email: t.bui@plaxis.com

D Branque
email: denis.branque@entpe.fr

Constitutive elasto-plastic model for fine soils in the unsaturated to saturated saturation zone

Keywords: Quasi-saturated soils, poro-elasto-plasticity, entrapped air in porous media

Constitutive elasto-plastic model for fine soils in the unsaturated to saturated saturation zone

Introduction

For the design and the construction of geostructures, the knowledge of the response of soils subjected simultaneously to mechanical and hydraulic loadings is essential. In general, soils can be considered as multiphase porous media composed of the solid phase and porous space. Two or even more fluid phases can be simultaneously present in the porous space. For most cases in geotechnics, the two fluid phases are a liquid phase (water with other species) and a gas phase (air, water vapor). The hydromechanical behavior of soils depends not only on the constitutive behaviour of the solid skeleton, but also on the interaction between different phases.

Many experimental results have shown the important role of the gaseous phase in the hydromechanical behaviour of soils [START_REF] Gallage | Laboratory measurement of hydraulic conductivity functions of two unsaturated sandy soils during drying and wetting processes[END_REF][START_REF] Francois | Compression tests on a sandy silt at different suction and temperature levels[END_REF][START_REF] Ho | Volume change indices during loading and unloading of an unsaturated soil[END_REF][START_REF] Fleureau | Behavior of clayey soils on drying?wetting paths[END_REF]. And this effect was already taken into account in numerous constitutive elasto-plastic models [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF][START_REF] Dg | The 1999 r.m. hardy lecture: the implementation of unsaturated soil mechanics into geotechnical engineering[END_REF][START_REF] Nuth | Effective stress concept in unsaturated soils: clarification and validation of unified framework[END_REF][START_REF] Sheng | Review of fundamental principles in modelling unsaturated soil behaviour[END_REF][START_REF] Zhou | An advanced hydro-mechanical constitutive model for unsaturated soils with different initial densities[END_REF]. Actually, at degrees of saturation lower than 85%, the gaseous phase is continuous. The coexistence of liquid water (a wetting liquid) and gaseous air (non-wetting gas) leads to the surface tension phenomenon, and a difference between their pressures, commonly called "suction". It results in an apparent attraction between grains, thus increasing the intergranular normal contact force. This has the effect of increasing the shear strength and the stiffness of the soil skeleton. If a fair number of theoretical and experimental studies have explored the behavior of unsaturated soils, it is not the case for the transition between partially saturated and fully saturated states, called the "quasi-saturated state", which is almost never considered. The lack of researches in this domain is not in line with its importance. Indeed, earthwork materials are usually compacted to their optimum Proctor density. Under this condition, the saturation degree is close to 80% and the gaseous phase, which is no longer continuous, takes the form of air pockets and bubbles entrapped within the liquid phase [START_REF] Christiansen | Effect of entrapped air upon the permeability of soils[END_REF]. The presence of this entrapped air seems to significantly influence the behavior of the soil. Several studies indicate that the entrapped air affects the soil's hydraulic and mechanical properties [START_REF] Stonestrom | Water content dependence of trapped air in two soils[END_REF][START_REF] Faybishenko | Hydraulic behavior of quasi-saturated soils in the presence of entrapped air: laboratory experiments[END_REF][START_REF] Boutonnier | Comportement hydromécanique des sols fins proches de saturation cas des ouvrages en terre: coefficient B, déformations instantanées et différées[END_REF]. However, at the degree of saturation mentioned above, it appears that neither unsaturated nor saturated formalisms can accurately reproduce the behavior of fine compacted soils. For that purpose, a first coupled hydro elastic model that takes into account the physical-mechanical interactions between different phases as well as the kinematics of each constituent (liquid water, dissolved air, gaseous air and solid grains) has been recently developed [START_REF] Lai | Poroelastic behaviour of fine compacted soils in the unsaturated to saturated transition zone[END_REF]. In particular, this model accounts for the interfacial tension, migration of gaseous and liquid phases, which have important impacts on the mechanical behaviour. The development leads to a system of highly non-linear partial differential equations which can be solved numerically using the finite element method. This model developed allows simulating the quasi-saturated soils behaviors and the continuous transition between the saturation domains. In order to complete the model presented in [START_REF] Lai | Poroelastic behaviour of fine compacted soils in the unsaturated to saturated transition zone[END_REF], the main objective of this paper is to present an extension of the theoretical developments to the elastoplastic behaviors of soils. At first, physical descriptions of air entrapment process in soil and air content evolution are presented. A simplification of the mass conservation equation, based on the analysis of characteristic time of the diffusion process of dissolved air is proposed. The second part of the paper is dedicated to the development of elasto-plastic law of behavior. For that purpose, a definition of the generalized effective stress for all saturation domains is introduced, as well as some assumptions on the evolution of loading surface with suction in the quasi-saturated domain. At last, the new model has been implemented in a numerical code Hydromech written in C++, developed originally by [START_REF] Pereira | Etude des couplages hydromécaniques et effets de non saturation dans les géomatériaux. Application aux ouvrages souterrains[END_REF], that has been used to simulate oedometer tests with different hydromechanical loading paths.

Description of entrapped air

Air bubbles and air pockets

The air phase is defined as "entrapped" when it is no longer in contact with the external atmosphere. As already discussed by Vaughan [START_REF] Vaughan | Observations on the behaviour of clay fill containing occluded air bubbles[END_REF], for a given entrapped air volume, all the air molecules can either be formulated as numerous spherical bubbles of radius r bubble embedded within the liquid phase or take a more complex geometry that closely follows the shape of the porous space to maximize its volume-to-surface ratio. A schematic diagram of the distinction between air bubbles and air pockets is reported in Figure 1.

The migration of air bubbles through the porous network can be due to their advection by the liquid phase or by their transport within the liquid phase, which can be regarded as a diffusion process. At constant pressure and temperature, they do not modify the free energy of the solid matrix. Their sole mechanical impact is to reduce the apparent compressibility of the liquid phase. The absence of interaction with the surface of the pores forces these bubbles to have a radius lower than the characteristic percolation size of the porous network, denoted by r perc . If not, their unobstructed movement across the porous space network becomes unlikely and the spherical geometry of the entrapped air volumes may be modified at their interactions with the surface of the pores. In this case, the entrapped air is assumed to form air pockets whose spatial migration can only take place by dissolution in liquid water followed by advective transport and Fickian diffusion of the dissolved air molecules (see, for example, [START_REF] Schuurman | The compressibility of an air/water mixture and theoretical relation between the air and water pressures[END_REF]). Due to evident geometrical constraints, the radius of the liquid-gas interface, r pocket , is lower than the radius of the pore, r p where the pocket is located. This description gives the following relation between the radii of the porous network and of the liquid-gas interfaces:

r bubble < r perc ≤ r pocket ≤ r p .
(1)

Equilibrium between gaseous and dissolved air molecules

The thermodynamic equilibrium between dissolved air and gaseous air is assumed to be reached anywhere at any time. Consequently, at constant pressure, the air pocket dissolution kinetics are driven by the velocity of dissolved air-molecule migration within the liquid phase, while the dissolution/bubbling processes induced by gas pressure variation are assumed to be instantaneous. Under the assumption of dilute solutions and perfect gas, the amount of air dissolved in liquid water is directly linked to the gaseous air pressure through Henry's law (cf. [START_REF] Lai | Poroelastic behaviour of fine compacted soils in the unsaturated to saturated transition zone[END_REF] for more details):

P G = K H M H2O M G meq , (2)
where P G is the gas pressure, meq is the mass of dissolved air per unit of liquid water mass at equilibrium with the entrapped air pockets, M H2O and M G are the molar mass of water and air, respectively, and K H is Henry's constant. At 20 o C, and assuming that air is composed of 20% O 2 and 80% N 2 , K H ≈7300 MPa. The use of the Henry's law requires the existence of stable air pockets within the pore network of the REV. To account for the full saturation case, we introduce m, which is the mass of dissolved air per unit of liquid water mass:

m = m aL m L , (3)
where m aL is the mass of dissolved air per unit of total initial volume (that is volume of solid matrix and pore space) and m L is the mass of liquid water per unit of total initial volume.

When the air in the gaseous state still exists within the REV, S L < 1, we have m = meq . When the fully saturated state is reached S L = 1 (that is when all the pore air in the REV has been dissolved), eq. ( 2) is no longer applicable. In this case, m ≤ meq . In this case the following conditions are considered:

meq ≥ m ≥ 0 ; 1 ≥ S L ≥ 0 ; ( m -meq ) (1 -S L ) = 0. ( 4
)
The relation [START_REF] Francois | Compression tests on a sandy silt at different suction and temperature levels[END_REF] indicates that the dissolution process is directly related to the gas pressure within the air pockets, which is in turn linked to the shape of the liquid-air interface through Young-Laplace's law:

P G -P L = 2γ L,G r , (5)
where P L is the liquid pressure, γ L,G is the surface tension of the liquid-gas interface, r is its radius, which is assumed to be spherical, r = r bubble for the air bubbles and r = r pocket for the air pockets. Since r bubble < r pocket , the combined use of eqs. ( 2) and ( 5) results in a higher concentration of dissolved air in the vicinity of the air bubbles than in the vicinity of the air pockets. As a consequence, dissolved air migrates from the bubbles toward the pockets, which ultimately makes the bubbles disappear in favor of the pockets. Because the characteristic diffusion time at the pore scale is much shorter than other characteristic times at the structural scale, the coalescence of small air pockets is assumed to be instantaneous.

In conclusion to this analysis, trapped air is assumed to take exclusively the form of air pockets (spherical or non-spherical) whose kinematics are similar to the solid skeleton's. To assess the geometry of these pockets, and especially of their interface with the surrounding liquid phase, it is necessary to make assumptions on how these pockets are created and their evolution.

From partial saturation to saturation: air entrapment mechanism

The physical origin of air entrapment can be explained by the diagram in Figure 2A: during wetting of an initially partially saturated soil, the smaller pores of the porous network are invaded by the liquid phase, which forces the evacuation of the gaseous phase out of the medium through the connected network of larger pores, which remain dried [START_REF] Faybishenko | Influence of entrapped on soils permeability: theory and experiment[END_REF]. When the saturation ratio increases above a certain threshold, denoted by S E e , which is close to 90% (cf. [START_REF] Faybishenko | Influence of entrapped on soils permeability: theory and experiment[END_REF], for example), all the pores with a radius lower than or equal to the percolation radius become saturated by water. At that stage, the gas's relative permeability reaches zero, but some air remains in the larger pores. This air can be considered as "entrapped" since all the capillaries that connect these pores to the rest of the porous network are saturated with water. The total volume of entrapped air within an REV with an initial volume of dΩ 0 just after connectivity of the gas phase is lost is equal to φ(1 -S E e )dΩ 0 , where φ is the Lagrangian porosity. It should be remembered that the capillary pressure, or suction, is defined as the difference between the nonwetting fluid (G) and the wetting fluid (L): P c = P G -P L . The value of P c needed to reach saturation S E e is denoted by P c,e and, assuming a spherical liquid-gas interface, a negligible thickness of the water adsorbed layer, and a wetting angle of 0 o , is equal to:

P c,e = 2γ L,G r perc . (6)
Rigorously, P c,e should be called "percolation suction" since it marks the transition between a continuous and a discontinuous gas phase. For the sake of simplicity, we assume here that P c,e can be confounded with entry air suction, which means that the hysteresis between drainage and wetting is ignored. According to this description, the interface radius of all the air pockets just after their creation during a wetting process is equal to r perc . After this stage, the geometry of the pockets may change due to dissolution/bubbling processes or to their mechanical deformation. Whatever these transformations, they take on a shape that maximizes their interface radius with the surrounding liquid phase. It is quite difficult to predict the evolution of the optimal interface radius with air content since the geometry of the air pocket may be quite complex. To overcome this problem, any decrease of air content is assumed to induce a reduction of the number of pockets rather than the decrease of their size; and vice-versa, the air content increase leads to the creation of new air pockets of optimized size, which may be notably induced by the so-called snap-off effect (see Figure 2B). This assumption seems to be in agreement with the microtomography observations made by [START_REF] Higo | Local waterretention behaviour of sand during drying and wetting process observed by micro x-ray tomography with trinarisation[END_REF]. It may be explained by the fact that the volume and the shape of air pockets are driven by the geometry of the pore where they are trapped. At equilibrium, the trapped air phase consists of numerous pockets with the same liquid-gas interface radii, which are as high as possible. When the air content decreases, some pockets may be more impacted by others. These pockets will then dissolve and feed the other (more stable) pockets so that they keep the same geometry.

In conclusion, all the stable pockets are assumed to have the same interface radius with the liquid phase, and this radius is assumed to remain constant and equal to r perc . A consequence of this assumption is that, even if the entrapped air is in the form of isolated pockets, its pressure is homogeneous within a RVE and the poromechanical mixture theory's approach, notably presented in [START_REF] Coussy | Mechanics and Physics of Porous Solids[END_REF], can be applied.

From saturation to partial saturation: air entrance mechanism

Creation of embedded air pockets and bubbles within the liquid phase during a drainage stage from the fully saturated state can result from two major mechanisms. The first one is the heterogeneous nucleation of dissolved air. However, a rough estimation using a formalism similar to the one presented in [START_REF] Osselin | Poromechanics of salt nucleation within a unsaturated reservoir rock[END_REF] leads to the occurrence of nucleation for a supersaturation of dissolved air within the range of 1.1-1.2, which is quite unlikely in geotechnical applications.

The second mechanism is the snap-off effect. However, this effect remains located near the external surfaces of the material and should occur for suction levels within the range of entry air suction. In conclusion, a direct transition, without the occurrence of trapped air, is assumed from the saturated to the unsaturated domains.

Conclusions on the definition of in-pore phases and the limits of saturation domains

For the sake of simplicity, the gravity term is not considered in the following development.

The above physical description leads to the consideration of two in-pore phases (liquid, L, and gas, G), each composed of two components: the liquid phase is made of liquid water (w) and dissolved air (aL), while the gas phase is made of dry air (aG) and water vapor (v). In addition, two types of distinct behavior have been identified for the gas phase: it can be either free if P atm -P L ≥ P c,e or trapped if P atm -P L ≥ P c,e . To facilitate the integration of this discontinuity into the equations, without loss of generality, the gaseous phase will conceptually be divided into two sets that can freely exchange within each other (see Figure 3). The first is the part that is invaded by water during the wetting process before the loss of connectivity of the gaseous phase. The gaseous phase within this volume is denoted with the superscript (m). The remaining part therefore corresponds to locations where the air is entrapped when the saturation attains the threshold S e . The air within this volume is denoted by the superscript (t). It is important to underline that this distinction between these two geometrical domains of the air phase is only made to simplify the derivations and that it is not necessary for the development of the model. Using these notations, the saturation ratio can be expressed as:

S L = 1 -S (m) aG -S (t) aG , (7)
where S L is the saturation ratio, S

aG is the fraction of the volume presently occupied by the gaseous phase that is invaded by the liquid phase when S L = S e (or when P atm -P L = P c,e ), while S (t) aG represents the complementary part that remains to be occupied by the entrapped air when S L = S e (or when P atm -P L = P c,e ).

The porous medium is considered "unsaturated" as long as the gaseous phase remains continuous and connected. This happens when the difference between the atmospheric pressure and the in-pore liquid pressure is higher than the entry air suction. Within this domain, the relationship between suction and saturation degree is depicted by the retention curve. A fairly large number of closed-form equations exist to reproduce this curve accurately for different materials. In this study we use the prominent equation provided by Brooks and Corey [START_REF] Brooks | Hydraulic properties of porous media and their relation to drainage design[END_REF] for saturation ratios ranging between 0 and S e : The combination of ( 7) and ( 8) gives:

S L = S e P atm -P L P c,e -a for P atm -P L ≥ P c,e . ( 8
) (t) (t) (m)
When the gaseous phase becomes discontinuous (no "free" air is left), S (m) aG becomes null and the gas pressure can exceed P atm . According to the description of air pockets presented above, and assuming that p v p aG , the link between the gas pressure and the liquid pressure is:

P G = P L + P c,e ≈ p aG for P atm -P L ≤ P c,e , (10)
where p aG is the partial pressure of dry air.

The quantity of entrapped air can change due to dissolution or bubbling. When all of the entrapped air is dissolved, the material becomes fully saturated. The mathematical conditions that define the fully saturated states are S aG = 0 and P atm -P L ≤ P c,e . This description leads to the shape of the retention curve shown in Figure 3.

The impact of the compaction on the transition between the unsaturated and the quasi-saturated domains is taken into account through the variation of P c,e with the void ratio. The general tendency, which is observed on alluvials soils by [START_REF] Vanapalli | The influence of soil structure and stress history on the soil-water characteristics of a compacted till[END_REF], and on finer soils by [START_REF] Cabarkapa | Automated triaxial apparatus for testing unsaturated soils[END_REF], is a decrease of P c,e when the void ratio increase. But this relation is not unique and must be experimentally determined for each soil.

3. Behavior of the in-pore phases

Mass conservation equations

To predict the behavior of both liquid and gas phases, it is necessary to consider two distinct mass conservations equations. In this paper, similarly to what is done in [START_REF] Lai | Poroelastic behaviour of fine compacted soils in the unsaturated to saturated transition zone[END_REF], the mass conservations of water (both liquid and vapour) and air (both gaseous and dissolved) are chosen. It leads to the following system of equation (For more detail, see [START_REF] Lai | Poroelastic behaviour of fine compacted soils in the unsaturated to saturated transition zone[END_REF]):

ρ 0 L e 1 + e 0 A e ∂e ∂t + A p ∂P L ∂t - ∂S (t) aG ∂t = ∇ • (D L ∇P L ) , (11)
ρ 0 L e 1 + e 0 S L ∂ m ∂t + B p ∂P L ∂t + B e ∂e ∂t + β 1 ρ G ρ 0 L -m ∂S (t) aG ∂t = ∇ • ( mD L ∇P L + ρ L D e ∇ m) , (12)
with

A p = S L K L - ∂S (m) aG ∂P L A e = S L e - ∂S (m) aG ∂e B p = m S L K L - ∂S (m) aG ∂P L + β 1 ρ G ρ 0 L S (t) aG P G B e = m S L e - ∂S (m) aG ∂e + β 1 S (t) aG ρ G ρ 0 L 1 e + 1 P G ∂P c,e ∂e ,
where, ρ L is the density of liquid water, K L its bulk modulus, ρ G is the density of the gas, e is the void ratio and e 0 is the initial void ratio. ∂S (m) aG /∂e and ∂S (m) aG /∂P L can be calculated in accordance with the relation [START_REF] Zhou | An advanced hydro-mechanical constitutive model for unsaturated soils with different initial densities[END_REF]. D e is the effective diffusion coefficient of dissolved air within the liquid water and D L is the global permeability coefficient of water, which is equal to:

D L = K L + c v ρ L δ p , (13)
with c v the mass concentration of vapor per unit of gas volume, K L the transport coefficient of liquid water, and δ p the transport coefficient of vapor. They satisfy:

K L = ρ 0 L κ 0 κ L r η L ; δ p = D v e M H2O RT , (14)
where κ 0 is the material's intrinsic permeability, in m 2 ; κ L r is liquid water's relative permeability coefficient, which is a nondimensional quantity that varies between 0 (at a low saturation ratio) and 1 (when the material is fully saturated); η L is the dynamic viscosity of liquid water; D v e is the effective diffusion coefficient of water vapor within the air phase, R is the perfect gas constant, T the absolut temperature and M H2O the molar mass of water. For the variation of K L with saturation, the relation proposed by Martin et al. [START_REF] Martin | Bm 3.2: Febex mock-up test[END_REF] is used, due to its simplicity and its quite good consistency with experimental data on fine soils:

K L = K 0 L (S L ) 3 , (15)
where

K 0 L = ρ 0 L κ 0 /η L is the value of K L at full saturation.
For the variation of δ p , Millington's expression [START_REF] Millington | Gas diffusion in porous media[END_REF] is arbitrarily chosen:

δ p = δ 0 p (1 -S L ) 7/3 , (16)
where δ 0 p is the value of δ p for a dried material.

Finally, β 1 is defined by:

β 1 = Y(P c,e -P atm + P L ), (17)
with Y(x) = 0 for x < 0 and Y(x) = 1 otherwise (step function). In other words, β 1 = 0 for the unsaturated domain and β 1 = 1 for the quasi-saturated and fully saturated domains.

The final step to close the system is to consider the equilibrium between gaseous and dissolved air, as is depicted by the relation [START_REF] Fleureau | Behavior of clayey soils on drying?wetting paths[END_REF]. In other words, m should be equal to meq = (P G M G )/(M H2O K H ) in the unsaturated and quasisaturated domains, while S (t) aG = 0 in the fully saturated domain. A necessary and sufficient condition to fulfill these two conditions is to consider the following form for the variation of S (t) aG :

∂S (t) aG ∂t = β 2 C p ∂P L ∂t + C e ∂e ∂t -∇ • m K L + c v ρ L δ p ∇P L + ρ L D e ∇ m , (18) where
C p = ρ 0 L e 1 + e 0 meq S L K L - ∂S (m) aG ∂P L + β 1 P G meq S L + ρ G ρ 0 L S (t) aG C e = ρ 0 L e 1 + e 0 meq S L + β 1 S (t) aG ρ G ρ 0 L 1 e -β 1 1 P G ∂P c,e ∂e -meq ∂S (m) aG
∂e and

β 2 = (1 + e 0 ) sign(S aG ) eρ 0 L meq -β 1 ρG ρL , (19)
with sign(S aG )=0 when S aG = 0 (when the material is fully saturated) and sign(S aG )=1 when S aG > 0 (the material is either unsaturated or quasi-saturated).

Simplification of the system of equation

The analysis of the relations ( 12) and ( 18) leads to the definition of τ D , which is the characteristic time of the dissolved air diffusion:

τ D ≈ e 1 + e 0 L 2 D v e , ( 20
)
where L is the characteristic length.

If τ D is strongly higher than the characteristic time of the solicitation and/or observation, the diffusion processes can be ignored, and the system of equation ( 11), [START_REF] Faybishenko | Hydraulic behavior of quasi-saturated soils in the presence of entrapped air: laboratory experiments[END_REF], and ( 18) is simplified to:

ρ 0 L e 1 + e 0 A e ∂e ∂t + A p ∂P L ∂t - ∂S (t) aG ∂t = ∇ • (D L ∇P L ) , (21a)
∂S (t) aG ∂t = β 2 C p ∂P L ∂t + C e ∂e ∂t -∇ • ( meq D L ∇P L ) . (21b)

Numerical application of the mass conservation equations

To illustrate the system formed by equations ( 11), [START_REF] Faybishenko | Hydraulic behavior of quasi-saturated soils in the presence of entrapped air: laboratory experiments[END_REF], and (18), a first calculation was performed considering an undeformable porous material (i.e., for ∂e/∂t = 0). This calculation considered the 1D configuration of a cylindrical 5-cm-high sample totally isolated on its lateral surfaces. This was done using the data reported in Table 1 and the calculation was performed with a home-made finite volume code computed in C.

The sample was initially in the unsaturated domain, with a saturation ratio equal to S 0 L = 0.0012 and therefore:

P L (x, t = 0) = P 0 = P atm -P c,e S e S 0 L (1/a) = 10kPa (22a) m(x, t = 0) = m0 = M G P atm K H M H2O = 2.2 × 10 -5 (22b) S (t) aG (x, t = 0) = 1 -S e (22c)
At t = 0 the material was put in contact with liquid water at one end of the sample. Its pressure linearly increased from P 0 = 10 [kPa] to P atm in 864 s and remained constant until the end of the simulation. The dissolved air content the material is under the unsaturated domain. The evolution of the saturation ratio is driven by the relation [START_REF] Sheng | Review of fundamental principles in modelling unsaturated soil behaviour[END_REF], and m remains constant and equal to m0 . When the material becomes quasi-saturated, the shape of the saturation ratio evolution changes and it is driven by the compressibility and dissolution of the air phase. The pressure of the gas phase is no longer equal to P atm , which then increases the dissolved air content. As a result, a gradient of dissolved air content is created between the core of the material and its boundary, which will progressively dissolve all the entrapped air. The propagation of this dissolution front, which is driven by the diffusion process of dissolved air, is very slow (less than 2 mm in 150 h). This numerical example thus clearly shows that even though the diffusion of dissolved air must be taken into account to properly model the long term behavior of quasi-saturated soils, the simplified system of equation ( 12) and ( 18) may be acceptable when short term behavior is assessed.

Mechanical equilibrium and rheological law

Sign convention and definition of the mechanical variables

Most of the elasto-plastic laws for unsaturated soils were developed with the sign convention of soil mechanics. Consequently, in this paper compressive stresses are taken as positives and positive strains express contraction.To derive the rheological law, the classical scalar variables of the Cam Clay model, namely the volumetric strain v , the mean stress σ v , the deviatoric strain d and the deviatoric stress σ d , are used. They are linked to the stress and strain tensors through the relations:

v = tr( ) ; d = 2 3 d : d , (23a)
σ v = 1 3 tr(σ) ; σ d = 3 2 σ d : σ d , (23b)
where and σ are the strain and stress tensors while d and σ d are their deviatoric parts, which are defined by:

d = - v 3 δ ; σ d = σ -σ v δ, ( 23c
)
where δ is the second order identity tensor.

For an axisymmetric problem, the relations (23) simplified in the form:

v = 1 + 2 2 ; d = 2 3 | 1 -2 | , (24a)
σ v = 1 3 (σ 1 + 2σ 2 ) ; σ d = |σ 1 -σ 2 | , (24b)
where 1 and σ 1 are the principal strain and the principal stress along the axis of symmetry, while 2 and σ 2 are the principal stress and principal strain along the radial and orthoradial directions.

Consistently with the classical strain decomposition assumption, volumetric and deviatoric strains can be split in elastic (superscript e) and plastic (superscript p) parts:

v = e v + p v ; d = e d + p d . (25)
Under the soil mechanics sign convention the relation between the effective stress and the total stress is:

σ = σ -Pφ δ ; σ v = 1 3 (σ 1 + 2σ 2 ) -Pφ , ( 26
)
where σ is the effective stress tensor, σ v is the mean effective stress, Pφ is the equivalent pore pressure. In the saturated regime, Pφ = P L , and σ is equal to Terzaghi's effective stress. In the quasi-saturated and unsaturated regime, the effect of the interfacial energy is ignored, Pφ = (1 -χ(S L ))P G + χ(S L )P L , and σ is equal to Bishop's effective stress. Finally, since this study focuses on the high saturation ratio, the assumption χ(S L ) = S L is used. On the basis of these restrictions, the extensive use of the relation ( 17)that define the coefficient β 1 as well as the ones that specify the values of

S (t) aG , S (m)
aG and P G for the different domain of saturation (equations 7-10) and bellow), allows to write Pφ and its temporal derivation in the forms:

Pφ = (β 1 + (1 -β 1 )S L )P L + (1 -β 1 )(1 -S L )P atm + β 1 S (t) aG P c,e , (27)
∂ Pφ ∂t = [β 1 + (1 -β 1 )(1 -a)S L ] ∂P L ∂t + P c,e ∂S (t) aG ∂t , (28)
where we recall that a is the coefficient of the Brooks and Correy retention curve.

Elastoplastic rheological law

To describe the elastic behavior, the isotropic logarithmic elastic behavior of the Cam Clay model is adopted:

d v = dσ v K ; K = 1 + e 0 k σ v , (29a)
d d = dσ d 3G ; G = 3(1 -2ν) 2(1 + ν) K, (29b)
where e 0 denotes the initial void ratio, K and G are the bulk and the shear moduli, while ν, which is assumed to be constant, is the Poisson ratio. This elastic law leads to a linear variation of the volumetric strain (or void ratio) with the logarithm of the mean effective stress during an elastic isotropic compression, and the slope of this variation is equal to k (cf. Figure 5).

To describe the plastic behavior, the general framework of isotropic strainhardening poroelastoplasticity (cf. [START_REF] Coussy | Revisiting the thermodynamics of hardening plasticity for unsaturated soils[END_REF]) is used with the following expression for the loading surface:

f (σ v , σ d , σ c ) = σ 2 d -M 2 σ v (σ c -σ v ) . (30)
In plane (σ v , σ d ), it leads to the elliptic criterion drawn in the Figure (6A). The size of the loading surface is driven by a constant M , which is the slope of the critical state line, and the hardening dual variable σ c . From a physical point of view, σ c is the preconsolidation stress. That is, during an isotropic compression load, no plastic strain occurs as long as σ v remains strictly lower than σ c . As a consequence, in the plane (σ v , P c ), the variation of σ c , which is described by the loading collapse curve (cf. Figure 6B), directly limits the elastic zone. Within the wide range of formulations available, a slight modification of Alonso et al.'s formulation [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF] is chosen to describe the variation of σ c with P c :

σ c σ ref = σ 0 σ ref l sat -k l(Pc )-k . ( 31
)
In this expression, σ 0 is the value of σ c at full saturation. l sat is the stiffness parameter of saturated soil that corresponds to the slope of the void ratio variation with the logarithm of the effective mean stress during an isotropic virgin compression (i.e., for σ v which is higher than the initial preconsolidation stress denoted by σ 0 c in the Figure 5, and hence σ c = σ v ). The l(P c ) function in (31) accounts for the increase of resistance and stiffness with suction. The choice of this function's formulation is particularly important because it determines the shape of the loading surface, which then affects the capacity of the model to reproduce wetting collapse. When the material is quasi-saturated, on the basis of experimental studies on lateritic and saprolitic soils [START_REF] Cruz | Peculiarities of tropical lateritic and saprolitic soils used as construction materials, selection, control and acceptance criteria-dams. Progress report of committee on tropical and saprolitic soils[END_REF] as well as Gangetic silt and Canyon Dam clay [START_REF] Shahu | Effective stress behavior of quasi-saturated compacted cohesive soils[END_REF], l(P c ) is assumed to remain equal to l sat . According to the relation [START_REF] Magnan | Etude théorique et expérimentale de la compressibilité du fluide interstitiel dans un sol argileux presque saturé[END_REF], this assumption leads to σ c = σ 0 when P c ≤ P c,e . Despite the high number of experimental studies conducted on this subject, no clear tendencies have been observed when the material is unsaturated. In some studies, l increases P c , while in others it decreases [START_REF] Sheng | Review of fundamental principles in modelling unsaturated soil behaviour[END_REF]. In this paper, to be consistent with the relation [START_REF] Magnan | Etude théorique et expérimentale de la compressibilité du fluide interstitiel dans un sol argileux presque saturé[END_REF], a form similar to the one proposed by Alonso et al. [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF] is adopted. It follows that the formulation of l(P c ) for the three domains of saturation can be described by:

l(P c ) = (1 -β 1 )l sat (1 -r) e -b(Pc-Pc,e) + r + β 1 l sat , ( 32
)
where b and r are two material parameters and β 1 is defined by the relation [START_REF] Schuurman | The compressibility of an air/water mixture and theoretical relation between the air and water pressures[END_REF].

A direct consequence of this form of l(P c ) is that σ c is equal to σ 0 in the quasi-saturated regime (cf. Figure 6). The physical meaning of this result is that no capillary hardening (or softening) is assumed as long as the air phase is entrapped, thus, from the point of view of the solid skeleton, the soil is saturated. Even if this result seems quite natural, it still needs to be experimentally verified.

Finally, the evolution of the preconsolidation pressure with the hardening variable, namely p v , is assumed to follow the Cam Clay incremental law:

dσ 0 = 1 + e 0 l sat -k σ 0 d p v ⇒ ∂σ c ∂ p v = 1 + e 0 l(P c ) -k σ c . ( 33
)
Assuming the associated flow rules, the combination of the consistency condtion (i.e. df = 0 during yield) with the relations (30), ( 31), [START_REF] Lai | A new constitutive model of unsaturated soils using bounding surface plasticity (bsp) and a nonassociative flow rule[END_REF], and [START_REF] Monnet | Calibration of an unsaturated air-water-soil model[END_REF] provides the following relations for the plastic strains, in which ξ denotes the ratio between the deviatoric and mean effective stresses (ξ = σ d /σ v ):

∂ p v ∂t = o Λ ∂f ∂σ v , (34a)
∂ p d ∂t = o Λ ∂f ∂σ d , ( 34b
)
where the plastic multiplier is equal to: with

o Λ= 1 H ∂f ∂σ v ∂σ v ∂t + ∂f ∂σ d ∂σ d ∂t + ∂f ∂P c ∂P c ∂t , (34c)
H = (1 + e 0 )(σ v ) 3 M 4 -ξ 4 l(P c ) -k (34d) and ∂f ∂σ v = σ v M 2 -ξ 2 , ( 34e
)
∂f ∂σ d = 2σ v ξ, ( 34f
) ∂f ∂P c = (β 1 -1) (l sat -k) l sat (1 -r) b e -b(Pc-Pc,e) ln σ0 σ ref (l(P c ) -k) 2 (σ v ) 2 M 2 -ξ 2 .
(34g) The combination of ( 25), [START_REF] Shahu | Effective stress behavior of quasi-saturated compacted cohesive soils[END_REF], and (34) eventually provides:

d v d d =    1 K + L H ∂f ∂σ v 2 L H ∂f ∂σ v ∂f ∂σ d L H ∂f ∂σ v ∂f ∂Pc , L H ∂f ∂σ v 1 3G + L H ∂f ∂σ d 2 L H ∂f ∂σ d ∂f ∂Pc    •   dσ v dσ d dP c   , (35)
where L = 1 if f = 0 (elastoplastic evolution) and L = 0 if f < 0 (elastic evolution).

Mechanical equilibrium

If the body forces are neglected, the local mechanical equilibrium leads to the relation:

∇ • ∂σ ∂t + ∂ Pφ ∂t δ = 0. ( 36
)
The injection of the elastic incremental laws ( 29) in [START_REF] Fabbri | Effect of temperature on the sorption curves of earthen materials[END_REF] leads to:

∇ • E T : ∂( -p ) ∂t + ∂ Pφ ∂t δ = ∇ • C T : ∂ ∂t + ∂ Pφ ∂t δ = 0, ( 37
)
where E T is the tangent fourth-order tensor of elasticity for which E T : d e = K e v δ + 2G e d , and is related to the displacement field, denoted by u through the relation:

= -(1/2)( t ∇u + ∇u). Finally, C T is the apparent tangent fourth-order elasto-plastic tensor of stiffness.

Numerical implementation

Final set of equations

Within the limits of short-term behaviors for which the diffusion of dissolved air is negligible, the final set of the equation is composed of the two partial differential equations (12-21b) and [START_REF] Noiret | Hydromechanical behavior of tournemire argillites: Measurement of the poroelastic parameters and estimation of the intrinsic permeability by oedometric tests[END_REF], where the two unknowns are the liquid pressure P L and the displacement field u.

The resulting model uses 12 parameters that can be classified into the following three groups:

• The first group consists of the eight parameters of the elastoplastic rheological law. The elastic parameter k, the slope of the NCL at saturation (l sat ), and the preconsolidation pressure (σ 0 ) can be directly determined by conventional isotropic compression tests on saturated soils. They can also be indirectly determined from a saturated oedometer test. The parameters b and r, which express the variation of l with capillary pressure, can be estimated from two additional isotropic compressions (or oedometer) tests, but under unsaturated conditions. The slope of the CLS (M ) and the Poisson coefficient (ν) can be determined from at least one triaxial test (several triaxial tests at several confinement pressures provides a better assessment). Finally, the initial void ratio e 0 can be estimated from the mass variation between fully saturated and dried states.

• The second group contains the three parameters of the retention curve: P c,e , a, and S e .

• The last group contains the two transport parameters, namely K L and δ p , which can be determined from classical standardized tests. If the longterm behavior is considered, the knowledge of an additional parameter, namely the effective diffusion of dissolved air (D v e ), is required.

Discretisation

The following developments will be limited to 1D problems under oedometric conditions along the x axis. In that case, all quantities only depend on the spatial coordiante x and time t while u = ue x and thus xx = -∂u/∂x is the sole non null componant of . For the stress field, only σ xx intervenes in the unique equilibrium condition ∂σ xx /∂x = 0. The external surface force T = -σ • n, with n the outgoing normal vector, is the external surface force, also reduces to a scalar T (the x-component) at each of the two extremities. The multiplication of the two equations respectively by a virtual displacement field, denoted by u * = u * e x , and a virtual pressure field P * followed by an integration leads, after simplification and while neglecting the variation S m aG with e, to the following expression for the variational formulation:

Ω (u * ) : C T : ( u) - ∂u * ∂x ∂ Pφ ∂t dΩ = Γ u * ∂T ∂t dΓ, (38)
Ω φρ 0 L P * A p ∂P L ∂t - ∂S (t) aG ∂t + S L φ ∂ u ∂x + ∂P * ∂x D L ∂P L ∂x dΩ = Γ P * q L dΓ, (39)
where Ω is the total volume of the domain and Γ is its surface. The notations φ = e/(1 + e 0 ) and u = ∂u/∂t are used to simplfy the reading of expressions, while q L = D L ∇P L • n is the water influx at the external surface.

Using the classic finite element discretization, the integral over the domain Ω is the sum of integrals over each element Ω i . The displacements and pressure fields inside an element at any time t is expressed as a combination of the elementary nodal displacements via a set of shape functions:

u(x, t) = k N u i,j U i,j ; u * (x, t) = k N u i,j U * i,j , (40)
P L (x, t) = k N p i,k P i,j ; P * (x, t) = k N p i,k P * i,k . (41)
The choice of the type of element for the spacial discretization is essential. As noted by [START_REF] Zienkiewicz | The Finite Element Method[END_REF], the interpolation degree for the displacement fields must be strictly higher than whose of pressure fields. Here, the same approach as the one developed in [START_REF] Pereira | Etude des couplages hydromécaniques et effets de non saturation dans les géomatériaux. Application aux ouvrages souterrains[END_REF] is used, which consists in considering linear elements for pressure fields (k = 2) and quadratic elements for displacement fields (j = 3). It leads to the spacial discretization represented in the Figure 7.

Applying the above discretization to the variational equations, and after assembling the elementary contributions, while accounting for the arbitrariness of the virtual field, we get:

D ∂U ∂t + Q ∂P ∂t = F u S + F u T , (42) Local indexing
R ∂U ∂t + S ∂P ∂t + T P = F p S + F p q (43)
where U and P are the global matrices of nodal displacements and pressure at the structural level, with D, Q, , R, S, T , F u S , F u T , F p S and F p q coming from the assemblage of elementary matrices

D i , Q i , , R i , S i , T i , F u S,i , F u T ,i , F p S,i and F p
q,i of the form:

D i = Ωi (B u i ) t C T B u i dΩ ; Q i = - Ωi (B u i ) t ((1 -β 1 )(1 -a)S L + β 1 ) N p i dΩ ; F u S,i = Ωi (B u i ) t β 1 P c,e

∂S (t) aG

∂t dΩ ;

F u T ,i = Ωi (N u i ) t ∂T ∂t dΩ ; R i = Ωi (N p i ) t ρ 0 L S L B u i dΩ ; S i = Ωi (N p i ) t φρ 0 L A p N p i dΩ ; T i = Ωi (B p i ) t D L B p i dΩ ; F p S,i = Ωi (N p i ) t φρ 0 L ∂S (t) aG ∂t dΩ ; F p q,i = Ωi (N p i ) t q L dΩ ;
where

N u i = [N u i,1 N u i,2 N u i,3 ] ; N p i = [N p i,1 N p i,2 ] ; B u i = ∂N u i ∂x ; B p i = ∂N p i ∂x .

Resolution scheme

This poro-elastoplastic system of equations was implemented in the code Hydromec1D initially developed by [START_REF] Pereira | Etude des couplages hydromécaniques et effets de non saturation dans les géomatériaux. Application aux ouvrages souterrains[END_REF]. The resolution uses a classical iterative method to compute the stress correction due to the plastic strain, based on the elastic stiffness. For this 1D problem, where computation time is not an issue, a very small time step is used to ensure convergence.

Results and discussion

The consistency of the results provided by the model was tested on a comparison with the experimental data provided by Magnan and Dang [START_REF] Magnan | Etude théorique et expérimentale de la compressibilité du fluide interstitiel dans un sol argileux presque saturé[END_REF], which is illustrated in Figure 8. This experiment is described in detail in [START_REF] Lai | A new constitutive model of unsaturated soils using bounding surface plasticity (bsp) and a nonassociative flow rule[END_REF]. Roughly it consists in an undrained oedometer test, with an axial stress ranging from 28 kPa to 800 kPa, on a quasi-saturated 5-cm-high sample 6.5 cm in diameter. To simulate this test, the system of equation (12-21b) and ( 37) is solved using an adaptation of the finite element code Hydromech developed by [START_REF] Pereira | Etude des couplages hydromécaniques et effets de non saturation dans les géomatériaux. Application aux ouvrages souterrains[END_REF]. In this experiment, the soil sample was taken from a site in Cubzacles-Ponts, in southeastern France, whose parameters are summarized in Table 2 (data provided by [START_REF] Monnet | Calibration of an unsaturated air-water-soil model[END_REF]).

Pressure sensor

Two simulations were performed: with and without air dissolution. The results are reported in Figure 9. A good match is observed between the calculated excess pore pressures and those from the experimental measurements (cf. Figure 9A). Moreover, it can be observed that the consideration of air dissolution significantly affects the model's prediction. For this example, and with the set of material parameters adopted, the simulation that takes into account air dissolution leads to a better estimation of the void ratio evolution, whereas the simulation ignoring air dissolution underestimates the deformability (see Figure 9B). The difference between the two simulations is significant because the soil under consideration is highly deformable. These results highlight the importance to take into account the air dissolution phenomenon for a precise estimation of the deformability of quasi-saturated soils. The same kind of results were obtained using an elastic approach in [START_REF] Lai | Poroelastic behaviour of fine compacted soils in the unsaturated to saturated transition zone[END_REF]. However, contrary to what it was made in [START_REF] Lai | Poroelastic behaviour of fine compacted soils in the unsaturated to saturated transition zone[END_REF], where it was necessary to adjust the value of the elastic parameters, the elasto-plastic parameters used to obtain the results of the Figure 9 were direclty the ones provided by [START_REF] Monnet | Calibration of an unsaturated air-water-soil model[END_REF]. In other words, no calibration of model parameter was realized. This result gives thus some confidence on the accuracy of this model when the soil remains in the quasi-saturated domain. It is however not sufficient to test the model's ability to correctly reproduce the transitions (both spacial and temporal) among the several domains of saturation and to analyze the importance of taking into account the elasto-plastic processes.

For this purpose, a theoretical simulation was undertaken considering an oedometric imbibition (null radial displacement and flows) at constant axial stress. An initially unsaturated homogeneous 5-cm-thick sample was subjected to a monotonic increase of liquid pressure at the top boundary surface (x = 0), from -250 [kPa] to 900 [kPa] in 150 [s]. At the bottom boundary (x = 5cm), a noflux condition was considered. The total axial stress was kept constant equal to its initial value σ xx (x, t) = 1MPa. Finally, the mass concentration of dissolved air in the incoming liquid water was, at any time, equal to the concentration of dissolved air within the pores of the sample, which were directly in contact with the external surfaces. The results of this simulation are shown in Figure 10.

Until t = 37s the whole sample is unsaturated (US). The liquid saturation profile is directly linked to the liquid pressure through the retention curve [START_REF] Sheng | Review of fundamental principles in modelling unsaturated soil behaviour[END_REF] because the gas pressure is kept equal to P atm . At x = 0m, the transition between the unsaturated (US) and quasi-saturated (QS) domains occurs for t = 37s. This is illustrated by the variation in shape of both S L (t) (see Figure 10A) and P G (t) curves (see Figure 10B). The former stems from the modification of the link between S L and P L , which is no longer driven by [START_REF] Sheng | Review of fundamental principles in modelling unsaturated soil behaviour[END_REF] in the quasi-saturated domain, but by air dissolution and compressibility, as depicted by the relation (21b). The latter is the consequence of the assumptions that the gas pressure is equal to the atmospheric pressure in the unsaturated domain and that the difference between the gas pressure and the liquid pressure becomes constant in the quasi-saturated domain (see eq. ( 10)). In line with [START_REF] Millington | Gas diffusion in porous media[END_REF], this augmentation of the gas pressure induces a reduction of the effective stress. It explains the increase of slope observed in Figure 10D, which denotes an acceleration of the void ratio increase when the material becomes quasi-saturated. In addition, as it is shown in Figures 10C andD, the decrease in effective stress induces an elastoplastic transition, which occurs while the material is under the quasi-saturated regime of saturation. In consequence, an accurate modeling of the quasi-saturated domain appears to be important for an accurate prediction of the soil deformation. At x = 0m and t = 145s, the value of the saturation ratio reaches 1, and the material becomes fully saturated. From then on, the saturation remains constantly equal to 1. This analysis demonstrates the ability of the model to reproduce temporal transitions of saturation domains, irrespective of the mechanical behavior of the solid skeleton (either elastic or elastoplastic). Finally, the ability of the model to reproduce the spacial transitions of saturation domains is illustrated in Figure 10A, in which a coexistence of the three saturation domains can be observed at t = 150s.

These two numerical examples demonstrate some of the abilities of the theoretical framework developed. However, the extent to which this degree of complexity is necessary to properly design earthworks near saturation can be questioned. Concerning the entrapped air dissolution process, the first numerical example shows the strong effects on both soil deformability and the evolution of entrapped air content with in-pore pressure. Ignoring this phenomenon would underestimate the deformability of the soil and to overestimate the air content within the soil. It would thus be prejudicial in the assessment of rutting, for example.

A great asset of this model is that it allows continuous and smooth spatial and temporal transitions between the saturation domains. The importance of the transition between the unsaturated and the quasi-saturated states is quite obvious. Although the top layer of an earthwork (such as an earth dam or a road or railway embankment) is quite often in an unsaturated state, the deeper layers remain quasi-saturated. Thus, interfaces between the saturation domains will be present, and a simulation that considers only one saturation domain will not accurately estimate the global deformation of the earthwork. In addition, using physical parameters of real soil, the model predicts a transition between the quasi-saturated and the fully saturated domains for an excess pore water pressure of 300 kPa, which corresponds to a 30-m column of water. In all cases in which such excess pore pressures may be encountered, the transitions between the saturation domains should be taken into account.

Other models allowing one to take into account the entrapped air phase and to consider transitions of the saturation domain have already been developed in the past. Compared to these models, the main addition of this new model is to provide a better physical description of the behavior of entrapped air by considering that it is attached to the solid skeleton rather than to the fluid phase and that its pressure is driven by the liquid-gas interface curvature, which tends to be as small as possible. However, without additional experimental data, no definitive conclusion can be drawn on the most accurate strategy to model quasi-saturated soils.

In addition, this model is still subjected to some major limitations. At first, only isothermal configurations are considered. As it is underlined by [START_REF] Francois | ACMEG-TS: A constitutive model for unsaturated soils under non-isothermal conditions[END_REF][START_REF] Yao | Non-isothermal unified hardening model: a thermo-elasto-plastic model for clays[END_REF], temperature variations would affect the elasto-plastic behavior of the soil, its retention curve (similarly to what is presented in [START_REF] Fabbri | Effect of temperature on the sorption curves of earthen materials[END_REF]) and dissolution/bubbling processes of entrapped air. Thus, spatial and temporal variations of the temperature field may significantly change the impact of the entrapped air on the global behavior of the soil. Secondly, this model assumes that the material is isotropic, which is unlikely due to the directional compaction method commonly used. A more realistic model would therefore assume a transverse isotropic behavior such as the one depicted by [START_REF] Noiret | Hydromechanical behavior of tournemire argillites: Measurement of the poroelastic parameters and estimation of the intrinsic permeability by oedometric tests[END_REF] for Tournemire argillites. Among the same line, if clay soils are considered, it might not be entirely justified to assume that the compressibility of the solid matrix is negligible. It would be necessary to introduce the Biot's coefficient in the expression of the effective stress and to consider the variation of porosity at constant macroscopic strain caused by the deformation of solid particles [START_REF] Coussy | Mechanics and Physics of Porous Solids[END_REF]. Finally, for a global hydromechanical assessment, the impact of the pore space deformations (plastic, thermal, cracks, etc.) on the transfer properties, which may be far from negligible [START_REF] Monfared | Temperature and damage impact on the permeability of opalinus clay[END_REF], as well as the hysteresis of the retention curve in the US domain [START_REF] Zhou | A contact angle-dependent hysteresis model for soil-water retention behaviour[END_REF] should be taken into account.

Conclusion

The behaviour of quasi-saturated materials is an important factor to be considered when designing cuttings and embankments in which earthwork materials are compacted to the optimum proctor density. Typically, soil compaction is performed at the optimum Proctor or on the wet side of the optimum, which means that the soil is in a highly saturated state. In such case, the gas phase may become discontinuous and take the form of entrapped air pockets. In this context, the construction of a theoretical model for this type of soils was presented in this paper. It has required the consideration of various physical-mechanical phenomena and their couplings occurring within the tri-phasic medium consisting of the solid grains, liquid water containing dissolved air and the entrapped air pockets. In particular, the conditions leading to the presence of these pockets within the porous network, their evolution and their migration were discussed. Although they are trapped in the pores, their dissolution can take place. Dissolved air can migrate through the pore space, either by following the flow of the fluid or by diffusion. The development leads to a system of highly non-linear partial differential equations which was solved numerically using the finite element method. A simplification of the system for short term behaviors was then proposed. On the basis of this simplified system, a new elasto-plastic hydromechanical model has been developed that takes into account the physical-mechanical interactions between these different in-pore phases and the solid matrix. For that purpose, a bishop effective stress with ξ(S L ) = S L was considered, and no variation of preconsolidation effective stress and of the slope of the NCL were assumed when the soil is quasi-saturated. This new model has been implemented in a numerical code Hydromech written in C++, developed originally by [START_REF] Pereira | Etude des couplages hydromécaniques et effets de non saturation dans les géomatériaux. Application aux ouvrages souterrains[END_REF] that has been used to simulate oedometer tests with different hydromechanical loading paths. In particular, this code allows to simulate consistently the transition across different regimes of saturation, both with respect to space (progressive translation of a boundary between two neighbouring regimes) and to time (transition of one regime to another at a fixed material point); which constituted a difficult modelling problem at the start. Numerical studies carried out show that this model gives consistent results providing a clear demonstration of its ability to simulate with precision the hydro-mechanical behaviour of quasi-saturated soils containing entrapped air.

Figure 1 :

 1 Figure 1: Schematic representation of the entrapped air pockets and air bubbles within the porous network

Figure 2 :

 2 Figure 2: A: Schematic representation of the air entrapment process. B: Diagram of the snap-off effect.

Figure 3 :

 3 Figure 3: Shape of the retention curve and diagram of the unsaturated and quasi-saturated domains

S (m) aG = S e 1 -= 1 -

 11 P atm -P L P c,e -a for P atm -P L ≥ P c,e , S e for P atm -P L ≥ P c,e .

Figure 4 :

 4 Figure 4: Liquid pressure (A), dissolved air content (B), and saturation ratio (C) profiles at three times obtained from the numerical application of the mass conservation equations for an undeformable porous medium.

Figure 5 :

 5 Figure 5: Variation of the void ratio with the logarithm of the mean effective stress during isotropic compressions at several suction values. σ 0 0 and σ 0 c stand for the initial saturated and unsaturated preconsolidation stresses, respectively.

Figure 6 :

 6 Figure 6: Evolution of the loading surface with suction in the (σ v , σ d ) and (σ v , Pc) planes.

Figure 7 :

 7 Figure 7: Diagram of an element i with 3 interpolation node for the displacement (nodes 1, 2 and 3), two nodes for the pressure (nodes 1 and 2) and two Gauss nodes.

Figure 8 :

 8 Figure 8: Diagram of the consolidated-undrained oedometer test on clay samples from Cubzacles-Ponts conducted by [31].

Figure 9 :

 9 Figure 9: Comparison between the model results and the data from Magnan and Dang's experiment

Figure 10 :

 10 Figure 10: Results of the numerical simulation of oedometric imbibition. US, QS, and FS stand for the unsaturated, quasi-saturated, and fully saturated domains, respectively. A: Saturation ratio at three locations over time. B: Evolution at x = 0 of liquid and gas pressures. C: Effective stress at the x = 0 surface over time. C: Evolution of the void ratio as a function of the effective stress at x = 0. The dotted line represents the transition between the elastic (EL) and the elastoplastic (EP) behaviors. D: Void ratio at three locations over time.

Table 1 :

 1 Numerical values used for the simulation.

Table 2 :

 2 Material parameter of the soil from Cubzac-les-Ponts from[START_REF] Monnet | Calibration of an unsaturated air-water-soil model[END_REF].