Zhan Yu

Jian-Fu Shao
email: jian-fu.shao@polytech-lille.fr

Minh-Ngoc Vu

Gilles Armand

Numerical study of thermo-hydro-mechanical responses of in situ heating test with phase-field model

Keywords: THM coupling, damage, plastic deformation, phase-field method, COx claystone, Heating test

This study is part of numerical simulations performed on an in situ heating test conducted by the French National Radioactive Waste Management Agency (Andra) at the Meuse/Haute-Marne Underground Research Laboratory (URL) to study the thermal-hydromechanical behavior of the host Callovo-Oxfordian COx claystone in quasi real conditions, through the international research project DECOVALEX. The emphasis is put on the pore pressure increase generated by the heat released from high level radioactive waste and the consequence on the damage evolution of the host rock. A phase-field method is proposed to describe the evolution of damaged zones around the heating borehole. Both tensile and shear cracks are taken into account. The evolution of damage is coupled with temperature variation, rock deformation, and fluid pressure change. Moreover, the structural anisotropy of elastic properties, permeability, and heat conductivity of the host rock is taken into account. Numerical results are compared with in-situ measurements.

Introduction

In the context of underground disposal of radioactive waste, rocks act as the geological barrier for the confinement of radioactive agents. The geological barrier is subjected to hydromechanical loading during and after construction of storage facilities and to thermal loading due to heat emitted by radioactive canisters. It is then crucial to investigate both the short and long thermal and hydromechanical (THM) responses of the geological barrier. In France, the Callovo-Oxfordian (COx) claystone formation has been selected as the potential geological barrier for both intermediate-level long-lived waste (IL-LLW) and high-level waste (HLW) [START_REF] Armand | Main outcomes from in situ thermo-hydro-mechanical experiments programme to demonstrate feasibility of radioactive high-level waste disposal in the callovo-oxfordian claystone[END_REF]. Extensive laboratory tests at sample scale have first been performed for the characterization of basic THM behaviors of the COx claystone, for instance [START_REF] Gasc-Barbier | Creep behavior of bure clayey rock[END_REF][START_REF] Mohajerani | A laboratory investigation of thermally induced pore pressure in the callovo-oxfordian claystone[END_REF][START_REF] Menaceur | The thermo-mechanical behaviour of the callovooxfordian claystone[END_REF][START_REF] Liu | Effects of deviatoric stress and structural anisotropy on compressive creep behavior of a clayey rock[END_REF][START_REF] Liu | Moisture effects on damage and failure of claystone under compression[END_REF][START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of callovo-oxfordian claystone: from experimental studies to model calibration and validation[END_REF][START_REF] Desbois | Deformation in cemented mudrock (callovooxfordian clay) by microcracking, granular flow and phyllosilicate plasticity: insights from triaxial deformation, broad ion beam polishing and scanning electron microscopy[END_REF][START_REF] Liu | Effects of relative humidity and mineral compositions on creep deformation and failure of a claystone under compression[END_REF][START_REF] Liu | Multi-step triaxial compressive creep behaviour and induced gas permeability change of clay rich rock[END_REF][START_REF] Liu | Mechanical behavior of claystone in lateral decompression test and thermal effects[END_REF]. On the other hand, a series of in situ experiments have been and are being realized under the coordination of the Andra in the Meuse/Haute-Marne Underground Research Laboratory (URL) [START_REF] Armand | Main outcomes from in situ thermo-hydro-mechanical experiments programme to demonstrate feasibility of radioactive high-level waste disposal in the callovo-oxfordian claystone[END_REF]. The goals of those experiments are the characterization of thermohydromechanical responses of the host rock during the different phases of geological disposal, and then to demonstrate the feasibility of constructing and operating a disposal site in the COx claystone formation. In particular, heat emitted by HLW packages causes temperature rise and pore pressure increase in the surrounding saturated porous rock essentially due to the difference of the thermal expansion coefficient between pore water (2.3 10 -4 K -1 at 20 • C and 7.2 10 -4 K -1 at 90 • C) and that of solid skeleton (1.5 10 -5 K -1 for the COx claystone) [START_REF] Mohajerani | A laboratory investigation of thermally induced pore pressure in the callovo-oxfordian claystone[END_REF][START_REF] Braun | Theoretical analysis of pore pressure diffusion in some basic rock mechanics experiments[END_REF][START_REF] Vu | Thermal pressurization coefficient of anisotropic elastic porous media[END_REF]. The induced pressure built-up is followed by a slow water flow phase owing to the low permeability of the host rock (1.0 10 -21 m 2 to 2.0 10 -20 m 2 ) [START_REF] Seyedi | On the thermo-hydro-mechanical pressurization in callovo-oxfordian claystone under thermal loading[END_REF][START_REF] Vu | Thermo-poro-mechanical coupled processes during thermal pressurization around nuclear waste repository[END_REF]. Several heating experiments from small-scale to full-scale have been conducted in the URL to study the THM behavior of the COx claystone in quasi real conditions, as well as to assess different aspects related to the conception and design optimization of disposal facilities [START_REF] Armand | Main outcomes from in situ thermo-hydro-mechanical experiments programme to demonstrate feasibility of radioactive high-level waste disposal in the callovo-oxfordian claystone[END_REF][START_REF] Conil | In situ investigation of the thm behavior of the callovo-oxfordian claystone[END_REF][START_REF] Bumbieler | Feasibility of constructing a full scale radioactive high-level waste disposal cell and characterization of its thermo-hydro-mechanical behavior[END_REF]. Among them, the full scale experiment, called ALC1604, has been conducted in the URL since 2009 in order to assess the constructability and the behavior of a full scale HLW cell subjected to a heat power similar to that of real waste packages. Thermo-mechanical responses of buried steel casing and THM responses(temperature, pore pressure, deformation) of the embedding claystone have been monitored [START_REF] Bumbieler | Feasibility of constructing a full scale radioactive high-level waste disposal cell and characterization of its thermo-hydro-mechanical behavior[END_REF]. Numerical modeling of those responses is one of the tasks defined for the DECOVALEX2019 projet, especially the Task E [START_REF] Seyedi | Upscaling thm modelling from small-scale to full-scale in-situ experiment in the callovooxfordian claystone[END_REF].

During the last decades, different kinds of numerical methods and computer codes have been developed for modeling THM coupling problems related to various engineering applications, for instance [START_REF] Taron | Thermal-hydrologic-mechanical-chemical processes in the evolution of engineered geothermal reservoirs[END_REF][START_REF] Tsang | Coupled thermo-hydro-mechanical processes in the near field of a high-level radioactive waste repository in clay formations[END_REF][START_REF] Yasuhara | Coupled thermo-hydro-mechanicalchemical modeling by incorporating pressure solution for estimating the evolution of rock permeability[END_REF][START_REF] Ogata | Modeling of coupled thermal-hydraulicmechanical-chemical processes for predicting the evolution in permeability and reactive transport behavior within single rock fractures[END_REF][START_REF] Zhou | Water saturation effects on dynamic fracture behavior of sandstone[END_REF] just to mention a few. To address the specific THM coupling issues related to the geological disposal of radioactive waste, an international research project DECOVALEX has been created and is working for a few decades. A number of advances have been achieved and a summary insight can be found in [START_REF] Tsang | Decovalex project: from 1992 to 2007[END_REF][START_REF] Birkholzer | Decovalex-2015: an international collaboration for advancing the understanding and modeling of coupled thermo-hydro-mechanical-chemical (thmc) processes in geological systems[END_REF][START_REF] Birkholzer | 25 years of decovalexscientific advances and lessons learned from an international research collaboration in coupled subsurface processes[END_REF]. More specifically, both continuous and discontinuous approaches have been proposed for solving THM coupling problems in saturated and partially porous media with elastic and inelastic models [START_REF] Garitte | Modelling the mont terri he-d experiment for the thermal-hydraulicmechanical response of a bedded argillaceous formation to heating[END_REF][START_REF] Garitte | Evaluation of the predictive capability of coupled thermo-hydro-mechanical models for a heated bentonite/clay system (he-e) in the mont terri rock laboratory[END_REF][START_REF] Graupner | Comparative modelling of the coupled thermal-hydraulic-mechanical (thm) processes in a heated bentonite pellet column with hydration[END_REF]. Those methods have been compared through different benchmarks. Other methods have been developed for modeling THM coupling in fractured rocks and with growth of induced cracks [START_REF] Rutqvist | A numerical study of thm effects on the near-field safety of a hypothetical nuclear waste repository-bmt1 of the decovalex iii project. part 3: Effects of thm coupling in sparsely fractured rocks[END_REF][START_REF] Mcdermott | Application of hybrid numerical and analytical solutions for the simulation of coupled thermal, hydraulic, mechanical and chemical processes during fluid flow through a fractured rock[END_REF][START_REF] Hokr | Real-case benchmark for flow and tracer transport in the fractured rock[END_REF][START_REF] Lu | Calibration of water-granite interaction with pressure solution in a flow-through fracture under confining pressure[END_REF]. Different kinds of methods have been used for dealing with the activation, nucleation and growth of cracks and fractures. For example, enriched finite element methods have been first developed to account for displacement discontinuities at the elementary level with the help of enriched shape functions [START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations, part 1: fundamentales[END_REF]. On the other hand, extended finite element method (XFEM) have been developed to deal with displacement discontinuities by using nodal enrichment techniques [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF][START_REF] Zeng | Numerical study of hydraulic fracture propagation accounting for rock anisotropy[END_REF][START_REF] Zeng | Study of hydraulic fracturing in an anisotropic poroelastic medium via a hybrid edfm-xfem approach[END_REF]. However, macroscopic fracturing of most rocks materials is generally due to the localization or coalescence of diffuse micro-cracks (damage). The proper description of the transition from diffuse damage to localized cracks remains an open challenge. In order to address this issue, the phase field method has been developed and attracted more and more attention. It is initially based on the variational principle for linear fracture mechanics [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and the optimal approximation methods of functionals with jumps [START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via t-convergence[END_REF][START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]. It has firstly been applied to modeling of tensile crack propagation in elastic materials [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. A number of extensions have been proposed [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF]. Due to the relative simplicity of numerical implementation in standard finite element computer codes, the phase-field method has successfully been applied to various engineering problems including multi-physics coupling [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. part ii. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elasticplastic solids[END_REF], finite deformation [START_REF] Borden | A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects[END_REF], coupling between damage and plasticity [START_REF] Fang | Phase field fracture in elasto-plastic solids: variational formulation for multi-surface plasticity and effects of plastic yield surfaces and hardening[END_REF]. However, the phase-field models based on the variational formulation are suitable for brittle cracking in elastic materials. For ductile cracking in plastic materials, in particular with a non-associative flow rule, a new formulation has been developed for coupled phase-field and plastic modeling using a rigorous thermodynamics framework [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF]. This new formulation is able to conveniently deal with ductile fracturing in plastic geological materials with a non-associative plastic flow rule. However, in most previous studies, only tensile cracks have been taken into account and generally related to tensile strain energy. But for rock-like materials subjected to compressive stresses, both tensile and shear cracks can be induced. Therefore, the objective of the present study is to propose a new phase-field method in order to be able to describe such complex cracking scenarios. On the other hand, in the framework of the project DECO-VALEX2019, and in particular of the Task E, different kinds of numerical methods have been adopted by different participating teams for the prediction of THM responses of Callovo-Oxfordian (COx) claystone during the in-situ heating test ALC1604. In the present study, the phase-field method is adopted and implemented in the finite element framework for the description of initiation and propagation of cracks. The description of cracking is further coupled with thermal and hydromechanical solutions. By using this coupled method, the THM responses of COx claystone during the excavation and heating phases of the ALC1604 heating experiment are estimated. In particular, the evolutions of damaged zones around the excavated borehole are described and analyzed. Numerical results are compared with in-situ measurements.

Summary of phase-field method

Regularized crack density function

We consider here rock-like materials subjected to complex loading paths including tensile and compressive stresses. According to previous experimental studies related to damage and cracking process in those materials [START_REF] Wong | Micromechanics of faulting in westerly granite[END_REF][START_REF] Wong | Analysis of crack coalescence in rock-like materials containing three flawsfpart i: experimental approach[END_REF][START_REF] Paterson | Experimental rock deformation-the brittle field[END_REF], both tensile and shear cracks can initiate and propagate depending on stress path and materials heterogeneities. In general, under tensile stresses or in triaxial extension conditions, tensile cracks represent the main process. Under triaxial compression conditions, in particular when the confining stress is high, closed shear cracks are induced and coupled with frictional sliding along cracks. Each set of cracks is denoted by Γ α (α = t, s respectively for the tensile and shear cracks). The total area of crack surfaces for each set is given by:

A Γα = Γα dA (1)
In practice, it is quasi impossible to follow individually the evolution of multiple cracks. Therefore, the cracks are approximated by a crack density functional Γ dα so that the total area of crack surfaces is calculated by [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF]:

A Γα ∼ = A Γ dα = Ω Γ dα (d α , ∇d α )dΩ (2)
The auxiliary variable d α ∈ [0, 1] is called the phase-field variable to represent the current damage state such that d α = 0 denotes the undamaged state and d α = 1 the entirely damaged state. Physically, the phase-field variable is equivalent to the damage variable introduced in classical continuum damage models. For instance, in order to describe the initiation and propagation of tensile and shear micro-cracks, two independent damage variables have been introduced in some continuum damage models [START_REF] Mazars | Continuum damage theory-application to concrete[END_REF][START_REF] Chen | Modelling of plastic deformation and damage in cement-based material subjected to desiccation[END_REF]. It is noted that the approximate total area of crack surfaces is calculated by a volume integral over the fixed initial geometrical domain Ω. Further, the crack density functional Γ dα is defined in a non-local manner through the gradient of phase-field variable ∇d α . In this study, the crack density functional proposed in [START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via t-convergence[END_REF] is adopted:

Γ dα (d α , ∇d α ) = d 2 α 2l d + l d 2 ∇d α .∇d α (3)
l d > 0 is a length scale parameter introduced for the regularization of localized phase-field problem.

It is here considered that crack propagation is an irreversible process. The evolution rate of total crack surface area in an arbitrary volume V ⊂ Ω is non-negative and it is assumed to be finite:

ȦΓ dα = d dt V Γ dα (d α , ∇d α )dV = V Γdα (d α , ∇d α )dV ≥ 0 (4)
The fulfilment of the non-negative condition for any arbitrary volume V implies that the evolution rate of crack density functional is non-negative: [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF]:

Γdα = d dt Γ dα (d α , ∇d α ) = d α l d ḋα + l d ∇d α .∇ ḋα ≥ 0 (5)

Governing equation of phase-field variable

The phase-field method applied to brittle crack propagation in elastic solids has generally been based on the variational approach to fracture mechanics [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. The phase-field variable is seen as the solution of the Euler equation of a variational problem [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. For modeling coupling between phase-field evolution and plastic deformation, a new approach has been proposed by [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF]. That approach is based on the balance equations of forces and energy as well as the irreversible thermodynamics. By applying that approach respectively to the tensile and shear cracks, the governing equations of two phase-field variables can be expressed as follows [START_REF] Yu | A double phase-field method for tensile and shear cracking in rock-like materials with plastic deformation[END_REF]:

   -∂ψ ∂d t ( e , V p , d t , d s ) -g t c [ d t l d -l d div(∇d t )] = 0, -∂ψ ∂d s ( e , V p , d t , d s ) -g s c [ d s l d -l d div(∇d s )] = 0 (6)
g t c and g s c denote the critical fracture energy for the tensile and shear cracks respectively.V p represents the set of internal plastic hardening variables. The scalar-valued functional ψ( e , V p , d t , d s ) is the total stored energy per unit volume. It is assumed that the total stored energy can be expressed as the sum of an elastic part and plastic part. Each part of the stored energy are dependent on two phase-field variables. Thus, one gets:

ψ( e , V p , d t , d s ) = W e ( e , d t , d s ) + W p (V p , d t , d s ) (7)
According to [START_REF] Liu | Moisture effects on damage and failure of claystone under compression[END_REF], the driving forces controlling the evolutions of phase-field variables are defined by the derivatives of the stored energy functional. Based on experimental data showing degradation of elastic properties by induced damage, it appears convenient to express the elastic strain energy of the damaged material as the degraded form of that of the undamaged material by including effects of phase-field variables. On the other hand, in order to consider the initial (structural) anisotropy in elastic behavior of rocks such as the COx claystone, the elastic strain energy of the undamaged material is here decomposed into a spherical part and a deviatoric part:

W 0 e ( e ) = W 0 sph ( e ) + W 0 dev ( e ) = 1 2 e : C 0 sph : e + 1 2 e : C 0 dev : e (8)
The spherical and deviatoric parts of the elastic stiffness tensor, C 0 sph and C 0 dev , are defined as follows:

   C 0 sph = k 0 reuss 1 1; C 0 dev = C 0 -C 0 sph (9)
C 0 denotes the total initial elastic stiffness tensor of the undamaged material. The coefficient k 0 reuss is called the equivalent Reuss bulk modulus which is a function of elastic coefficients.

As an examples, we consider here an orthotropic material characterized by 3 elastic moduli (E 1 , E 2 , E 3 ), 3 Poisson's ratios (ν 12 , ν 23 , ν 31 ) and three shear moduli (G 12 , G 23 , G 31 ). In that case, the Reuss bulk modulus is given by:

1 k 0 reuss = ( 1 E 1 + 1 E 2 + 1 E 3 ) -2( ν 12 E 1 + ν 23 E 2 + ν 31 E 3 ) (10)
On the other hand, the effects of open or closed micro-cracks on the elastic behavior of the damaged material are generally different. In order to consider that difference, the following crack status coefficient r ± e is introduced:

   if tr( e ) ≥ 0 : r + e = 1, r - e = 0 if tr( e ) < 0 : r + e = 0, r - e = 1 (11)
It is assumed that the spherical part of elastic strain energy (related to bulk modulus in an isotropic material) is affected by the tensile crack only while the deviatoric part (related to the shear modulus in an isotropic material) is affected by both the tensile and shear cracks. The elastic strain energy function of the damaged material is thus expressed by:

W e ( e , d t , d s ) = {r + e h 1 (d t ) + r - e }W 0 sph + {r + e h 2 (d t ) + r - e h 3 (d s )}W 0 dev ( 12
)
Three functions h 1,2,3 are introduced to define the degradation of the elastic strain energy.

Based on previous studies [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field fe implementations[END_REF][START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF]. The following common form is adopted:

           h 1 (d t ) = (1 -d t ) 2 h 2 (d t ) = (1 -d t ) 2 h 3 (d s ) = (1 -d s ) 2 (13)
In plastic materials with damage evolution, the total plastic strain energy σ : ˙ p can be decomposed into two parts: one part is related to the plastic dissipation D p and another part is stored for plastic hardening A p Vp . It is assumed that the stored plastic energy can contribute to the evolution of phase-field variables. However, it is not easy to exactly calculate the quantity of stored plastic energy involved in the damage evolution. Based on previous studies [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF] and for the sake of simplicity, the stored plastic energy contributing to the damage evolution is calculated as a portion of the total plastic strain energy. Therefore, the stored plastic energy for the undamaged material is calculated by:

Ẇ 0 p ( ˙ p ) = ησ : ˙ p (14)
0 η 1 is called the Taylor-Quinny coefficient. It is further assumed that the effects of the plastic stored energy on the evolution of open or closed micro-cracks can be different. Thus, like for the elastic strain energy, the following crack status coefficient is introduced for the plastic stored energy:

   if tr( p ) ≥ 0 : r + p = 1, r - p = 0 if tr( p ) < 0 : r + p = 0, r - p = 1 (15)
Moreover, it is generally observed that the plastic part of stored energy contributes to the evolution of phase-field variables only in the case of dilatant plastic volumetric strain. Therefore, the plastic part of the stored energy of damaged material is given by:

Ẇp ( p , d t , d s ) = {r + e h 2 (d t ) + r - e h 3 (d s )}r + p Ẇ 0 p + r - p Ẇ 0 p (16)
The accumulated value of W p is calculated by the summation of each incremental value during the whole loading history. By substituting ( 12) and ( 16) for [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of callovo-oxfordian claystone: from experimental studies to model calibration and validation[END_REF] and making the derivatives of ψ( e , V p , d t , d s ) with respect to d t and d s , the governing equations ( 6) for the two phase-field variables become:

   -r + e h 1 W 0 sph -r + e h 2 (W 0 dev + r + p W 0 p ) -g t c ( 1 l d d t -l d div(∇d t )) = 0, -r - e h 3 [W 0 dev + r + p W 0 p ] -g s c ( 1 l d d s -l d div(∇d s )) = 0 (17)
On the other hand, according to [START_REF] Yu | A double phase-field method for tensile and shear cracking in rock-like materials with plastic deformation[END_REF], the non-negative fracture dissipation rates for the tensile and shear cracks are expressed as follows:

D f -t = g t c Γt d ≥ 0, D f -s = g s c Γs d ≥ 0 (18) with Γt d = d t l d ḋt + l d ∇d t .∇ ḋt ≥ 0, Γs d = d s l d ḋs + l d ∇d s .∇ ḋs ≥ 0 ( 19
)
As mentioned in [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF] for a single phase-field variable, the fracture dissipation rates for the two phase-field variables both contain a non-local term due to the fact that Γt d (and Γs d ) depends on ∇d t (or ∇d s ). This non-local nature of fracture dissipation makes that the solution to the phase-field problem is not mesh sensitive. Further, as g t c > 0 and g s c > 0 , the fact that Γt d ≥ 0 (and Γs d ≥ 0) implies D f -t ≥ 0 (and D f -s ≥ 0), and ensures the fulfilment of the second thermodynamics principle. This physically means that the creation of crack surfaces is an irreversible process and the crack surface density function is not decreasing. Considering the definition of Γt d (and Γs d ) given in [START_REF] Taron | Thermal-hydrologic-mechanical-chemical processes in the evolution of engineered geothermal reservoirs[END_REF], the above non-local thermodynamics consistency enforces that the local phase-field variables are also not decreasing ( ḋt ≥ 0 and ḋs ≥ 0). According to [START_REF] Bumbieler | Feasibility of constructing a full scale radioactive high-level waste disposal cell and characterization of its thermo-hydro-mechanical behavior[END_REF], the crack irreversibility implies that the phase-field variables increase only if the values of energy driving forces increase. Inspired by the concept of energy history functional proposed in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field fe implementations[END_REF], the following history functionals are defined to represent the maximum values of energy driving forces reached during the loading history:

           r + e W 0 sph -H sph ≤ 0, Ḣsph ≥ 0, Ḣsph (r + e W 0 sph ) = 0 r + e (W 0 dev + r + p W 0 p ) -H t dev ≤ 0, Ḣt dev ≥ 0, Ḣt dev (r + e (W 0 dev + r + p W 0 p )) = 0 r - e (W 0 dev + r + p W 0 p ) -H s dev ≤ 0, Ḣs dev ≥ 0, Ḣs dev (r - e (W 0 dev + r + p W 0 p )) = 0 (20)
Accordingly, the governing equations of the phase-field variables are rewritten as follows:

-h 1 H sph -h 2 H t dev -g t c ( 1 l d d t -l d div(∇d t )) = 0, -h 3 H s dev -g s c ( 1 l d d s -l d div(∇d s )) = 0 (21)
Finally, the standard derivative of the elastic strain function leads to the stress tensor in the damaged material:

σ = ∂ψ ∂ e = C b (d t , d s ) : e (22)
The elastic stiffness tensor of the damaged material (in drained condition for porous materilas) is given by:

C b (d t , d s ) = [r + e h 1 (d t ) + r - e ]C 0 sph + [r + e h 2 (d t ) + r - e h 3 (d s )]C 0 dev (23)

Coupling of THM and damage fields

In this work, the porous medium is considered an open thermodynamic system that exchanges fluid mass and heat with the outside. The state variables used are therefore the elastic strain tensor e , the plastic strain tensor p , the fluid mass variation per unit of initial volume of the porous medium m, the temperature variation θ, as well as the two phase-field variables d t and d s . By postulating the existence of a thermodynamic potential for this open system in terms of these state variables [START_REF] Coussy | Mechanics of Porous Continua[END_REF][START_REF] Coussy | [END_REF], the thermo-poro-elastic constitutive relations are given by:

(σ -σ 0 ) = C b (d t , d s ) : e -B(p -p 0 ) -A b θ (24) (p -p 0 ) = M (-B : e + m e ρ f 0 ) + 3α m M θ (25) (s -s 0 ) = s 0 m m e + α b : (σ -σ 0 ) -(3α m -B : α b )(p -p 0 ) + C b σ T 0 θ (26)
σ 0 , p 0 , s 0 and T 0 denote the values of stress, fluid pressure, entropy and temperature at the reference state. The variation of temperature is defined as θ = T -T 0 . ρ f 0 is the volumetric density of fluid at the reference state. m e denotes the elastic variation of fluid mass per unit volume of the initial configuration. The scalar coefficient M is the Biot modulus and B is the second order tensor of Biot coefficients. α b is the tensor of thermal dilatation coefficients in drained conditions while A b denotes the tensor of thermoelastic coupling coefficients. C b σ is the volumetric heat for constant stress under drained conditions. The coefficient α m denotes the differential thermal dilation of saturated porous medium. In this study, though the initial anisotropy of elastic properties of COx claystone is taken into account, due to the difficulty to get feasible experimental data and for the sake of simplicity, the tensors of Biot's coefficients and thermal dilation coefficients are taken as isotropic ones. Therefore, one has B = bI and α b = α b I. Accordingly, the tensor of thermoelastic coupling coefficients becomes A b = 3K b α b I, with K b being the equivalent drained bulk modulus which is dependent on the damage variable as described below. α m is the differential thermal dilation coefficient whose expression is detailed below.

The total stress field verifies the following equilibrium equations at all points of the porous medium:

div(σ) + f v = 0 (27)
The fluid flow in the porous medium is described by the Darcy's conduction law and the mass balance equation. Together with the constitutive relations, the diffusivity equation for fluid flow is given by:

div( k µ .∇p) = 1 M ∂p ∂t + b ∂ kk ∂t -3α m ∂θ ∂t (28)
k is the tensor of permeability coefficients while µ the dynamic viscosity of fluid. On the other hand, by using the Fourier's diffusion law and the energy balance equation, the heat diffusivity equation can be written as follows:

λdiv(∇θ) = C b σ ∂θ ∂t -3α m T 0 ∂p ∂t + 3α b K b T 0 ∂ kk ∂t -k∇p(C p ∇θ) (29)
In this Equation, K b is the equivalent drained bulk modulus of the damaged material and defined by

K b = k reuss = k 0 reuss (1 -d t ) 2 for r + e = 1.
The induced damage affects the thermal and hydromechanical properties of materials. The degradation of elastic properties is already described by the constitutive relations given above. In this study, it is further assumed that the permeability tensor is affected by the tensile damage which represents the open cracks. Without making a rigorous homogenization procedure, the following empirical relation is adopted:

k = k initial exp(βd t ) (30)
The coefficient β controls the evolution rate of permeability. Further, the porosity and the Biot's coefficient are also considered as functions of the tensile damage:

φ(d t ) = φ initial + (1 -φ initial )d t (31) b(d t ) = b initial + (1 -b initial )d t (32)
Accordingly, the Biot's modulus varies with the tensile damage and is calculated by:

1 M (d t ) = (b -φ) K s + φ K f ( 33
)
with K s and K f being the bulk moduli of the solid phase and fluid, both taken as constant.

Finally, the differential thermal dilation coefficient is given by:

α m (d t ) = (b -φ)α b + φα f ( 34
)
α f is the thermal dilation coefficient of the fluid.

Numerical implementation in finite element method

The THM and damage fields are determined by using the standard finite element method. Due to the transitional regime of fluid flow and heat conduction and also the mechanical nonlinearity related to damage evolution and plastic deformation, the boundary values problems related to each field are solved by an incremental method. The total loading history is divided into a finite number of time steps. For the current step, the time increment is denoted as ∆t.

By making the geometrical meshing and using the standard shape functions related to the physical fields to determine and for the element type selected. One gets the following system of linear equations to solve for the displacement field [START_REF] Shao | Study of massive water injection by thermoporomechanical coupling model[END_REF]:

                           R uu ∆U n+1 + C up ∆P n+1 + C uT ∆T n+1 = ∆F n+1 e R uu = Ω ( t B u )C b (d t , d s )B u dV C up = Ω ( t B u )(-b(d t ))N p dV C uT = t Ω B u (-3α b K b (d t ))N p dV ∆F e = s t ( t N u )∆tdS + Ω ( t N u )∆f v dV (35)
U, F e , ∆P and ∆T are respectively the column matrix of increment nodal displacements, nodal mechanical forces, fluid pressure and temperature. In the above relations, N u denotes the matrix of shape functions for approximation of the displacement field, while B u is the matrix of derivatives of shape functions to calculate the corresponding strain field. N p is the matrix of shape function for approximation of the fluid pressure field.

Similarly, the system of discrete equations for the pore pressure field for the current loading step is given by:

                                   (∆t n+1 R pp + M pp )∆P n+1 + C pu ∆U n+1 + C pT ∆T n+1 = ∆t n+1 (-R pp P n + F n+1 ω ) R pp = Ω ( t B p )(k(d t )/µ)B p dV M pp = Ω ( t N p )( 1 M (d t ) )N p dV C pu = -t C up C pT = Ω ( t N p )(-3α m (d t ))N p dV F ω = -sω ( t N p )w n+1 .ndS (36)
B p is the matrix of derivatives of shape functions N p to calculate the corresponding gradient of fluid pressure. Finally, the system of equations for the temperature field for the current loading step is given by:

                                     (∆t n+1 R T T + M T T )∆T n+1 + C T u ∆U n+1 + C T p ∆P n+1 = ∆t n+1 (-R T T T n + F n+1 q ) R T T = Ω ( t B p )( λ T 0 )B p dV C T u = -t C uT M T T = Ω ( t N p )( C b σ T 0 )N p dV C T p = t C pT F q = -sq ( t N p )(q n+1
.n/T 0 )dS (37) By using the same finite element mesh as for the THM fields, the boundary values problems (see Appendix for the mathematical formulations) for the determination of two damage variables are also transformed into the systems of discrete equations. Namely, for the tensile phase field, one gets:

           K d t d t = F d t K d t = Ω {(g t c l d + 2H sph + 2H t dev )N T d N d + g t c l d B T d B d }dΩ F d t = Ω 2(H sph + H t sph )N T d dΩ (38)
And for the shear phase field, one has:

           K d s d s = F d s K d s = Ω {(g t c l d + 2H s dev )N T d N d + g t c l d B T d B d }dΩ F d s = Ω 2H s dev N T d dΩ (39)
N d and B d are the shape functions and their derivatives for the two phase fields.

As mentioned above, some THM parameters are functions of two damage variables. Therefore, THM problems are inherently coupled with phase-field problems. However, based on preliminary numerical studies, it is found that an explicit decoupled algorithm appears very efficient to solve these coupled problems. That means that the values of all THM parameters affected by the damage processes at the current loading step (n+1) are calculated by using the values of two phase fields at the previous steps n.

Behavior of COx claystone in conventional triaxial compression test

Before the application of phase-field method to THM modeling of the ALC1604 experiment, the mechanical behavior of the Callovo-Oxfordian claystone in conventional triaxial compression test is studied first. The objective here is to capture the transition from diffuse damage to localized cracks in the tested sample by using the proposed phase-field method. The COx claystone generally exhibits an elastic-plastic behavior [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of callovo-oxfordian claystone: from experimental studies to model calibration and validation[END_REF]. However, by putting the emphasis on the damage and cracking modeling, the plastic behavior is here described by a simple model including the Drucker-Prager linear yield function, an isotropic hardening law and a non-associated plastic potential. They are summarized as follows:

           Yield function: f p = q + α p (p -C) ≤ 0 Hardening law: α p = α 0 + (α m -α 0 ) γp B+γp Plastic potential: G = q + η c (C -p)ln( C-p P ) (40)
p = tr(σ)/3 is the mean stress. q = 3(s : s)/2 is usually called as the deviatoric stress in rock mechanics (while called as the equivalent tensile or von Mises stress in continuum mechanics), with s being the tensor of deviatoric stress. C is the yield threshold in hydrostatic tension. α p denotes the current frictional coefficient varying from the initial value α 0 to the asymptotic one α m . The parameter B controls the plastic hardening rate. And η c is a dilatancy coefficient. γ p is the equivalent plastic shear strain used as the plastic hardening variable.

On the other hand, for the description of cracking process by using the phase-filed method, it is needed to identify the critical fracture energy parameters g t c and g s c . In general, their values cannot be directly measured from conventional laboratory tests. As an alternative method, the value of g t c can be calibrated from the peak stress in a uniaxial tension test or the peak force in a three or four points bending test. The value of g s c is generally fitted from the peak stress of uniaxial or conventional triaxial compression tests. However, according to experimental studies [START_REF] Evans | The brittle-ductile transition in rocks: Recent experimental and theoretical progress, The Brittle-Ductile Transition in Rocks[END_REF][START_REF] Wong | The brittle-ductile transition in porous rock: A review[END_REF], rock-like materials exhibit a transition from brittle to ductile failure with the increase of confining stress. Moreover, according to micro-mechanical studies [START_REF] Hu | A micromechanics-based elastoplastic damage model for rocks with a brittle-ductile transition in mechanical response[END_REF][START_REF] Zhao | A micromechanics-based plastic damage model for quasi-brittle materials under a large range of compressive stress[END_REF], the brittle-ductile transition can be related to the increase of the local frictional coefficient of closed micro-cracks surfaces with confining stress. In this study, the brittle-ductile transition of the COx claystone is described by the increase of g s c with confining stress. Inspired by the previous study on granite [START_REF] Wang | Modeling of damage and cracking in heterogeneous rock-like materials by phase-field method[END_REF], the following empirical relation is adopted: g s c = g s c0 e bc(σc/σr) , with g s c0 = 4.25 × 10 -4 kN/mm, b c = 0.0977; σ c denotes the confining stress while σ r = 1M P a is a fixed normalizing coefficient. The model's parameters used in the simulation of conventional triaxial compression tests are given in Table 1 Each tested sample is considered as a small structure subjected to the boundary conditions as shown in Figure 1(a). The vertical (also called axial) displacement is blocked at the bottom surface and prescribed on the top surface. The axial displacement increment adopted here is ∆u = 2 × 10 -3 (mm). Moreover, in order to facilitate the initiation of macroscopic crack, a weak region is placed inside the sample, as shown in Figure 1(b). The size of weak region is equal to that of one element. In the weak region, the values of critical fracture energies g t and g s are reduced to 1% of their initial values. All simulations are performed under plane strain conditions. The reference mesh is composed of 5000 quadrilateral elements and 5151 nodes (called mesh 1).

(a) (b) In Figure 2, average axial and radial strains versus deviatoric stress curves are presented for three conventional triaxial compression tests with different confining pressures. A good agreement is obtained between the numerical results and experimental data found in [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of callovo-oxfordian claystone: from experimental studies to model calibration and validation[END_REF] for both the pre-and post-peak regimes. The consequence of induced cracks on the macroscopic behavior is well captured. The distributions of tensile and shear damage variables are presented in Figure 3 for the test with a confining pressure of P c = 12 MPa. For the sake of comparison, the two damage components are presented with the same scale bar. It is found that under conventional triaxial compression conditions, the tensile damage remains very small and difficult to distinguish in the Figure . The shear cracking is the dominating process controlling the progressive failure of the sample.

In order to verify the mesh sensitivity of numerical results, a comparative calculation is performed by using 1250 elements (called mesh 2) for the test with 6 MPa confining stress. The obtained averaged stress-strain curves are presented in Figure 4. One can see that the macroscopic mechanical responses obtained with the two mesh sizes nearly coincide. It seems that the numerical results obtained with the proposed phase-field model are almost insensitive to mesh size. It is worth noticing that the simulation of conventional triaxial compression tests presented above shows the efficiency of the proposed phase-field model for the description of mechanical behavior of the COx claystone. It would be desirable to verify the efficiency of the proposed model for the description of hydromechanical behavior through suitable laboratory tests, for instance, conventional triaxial compression tests with fluid injection. Unfortunately, such tests are not available on the COx claystone. The efficiency of the proposed hydromechanical coupling model is verified through the in situ ALC1604 experiment presented in the next section. Further, in the previous study [START_REF] Stanchits | Fracturing of porous rock induced by fluid injection[END_REF], the authors studied cracking patterns by monitoring acoustic emission (AE) activity during water injection into initially dry stressed cylindrical specimens of Flechtingen sandstone. Narrow AE activity concentration bands were identified. Those bands were in agreement with the macroscopic fractures observed on the tested samples. The formation of such localized cracked bands is well captured by the proposed phase-field model. However, the orientation and thickness of cracked bands should depend on the initial heterogeneities (inclusions and pores for instance) of rocks. This feature will be investigated in our future works.

THM analysis of ALC1604 test

The ALC1604 experiment is an ongoing in situ heating test conducted by Andra in the Meuse/Haute-Marne underground research laboratory. The experiment is a full scale representation of a single high-level waste cell in the Callovo-Oxfordian claystone. As shown in Figure 5, the ALC1604 borehole was drilled from the GAN gallery and is oriented parallel to the major horizontal principal stress σ H . More details of the ALC1604 heating test can be found in [START_REF] Bumbieler | Feasibility of constructing a full scale radioactive high-level waste disposal cell and characterization of its thermo-hydro-mechanical behavior[END_REF] published in this same special issue. Only some key input data are given here to facilitate the discussion of numerical results obtained by using the phase-field model.

The main heating phase started on 18 April 2013, at a constant nominal power of 220 W/m for the 15 m occupied by the heater elements, at a depth of between 10m and 25m in the cell from the wall of GAN drift. This value has been designed to reach a maximum temperature of 90 o C at the casing wall after 2 years of heating.

Input data of modeling

In this study, the ALC1604 experiment is simulated by using a two-dimensional (2D) plane strain configuration considering a representative cross section as given in Figure 6(a). Due to the double symmetry with respect to the horizontal and vertical axes, only a quarter of the full section is considered. Further, considering the symmetry conditions between two neighboring alveolus, the adiabatic and impermeable conditions are prescribed on the right external boundary as well as on the horizontal and vertical axes while the constant values of pore pressure and temperature are prescribed at the top boundary. About the mechanical boundary conditions, the constant initial stresses are prescribed on the right and top boundary surfaces. The reference finite element mesh (called mesh 1) contains 1560 elements and 1650 nodes and is shown in Figure 6(b). More precisely, 24 elements are used to mesh the 1/4 part of the ALC1604 wall. The size of the smallest elements used in the reference mesh is 0.0271×0.0147m. The scale length for the phase-field model is taken as l d =0.0147m, which is equal to the smallest size of element. The values of two critical fracture energies can exhibit some size effects and also be influenced by water saturation degree and temperature. Unfortunately, no in situ measurements of these parameters are available. Considering that the elastic modulus and mechanical strength of in situ rock masses are usually smaller than those of rock samples used in laboratory tests due to the presence of initial defects, the values of g t c and g s c used in the in situ heating experiment are taken as smaller than those used in conventional triaxial compression tests and equal to: g t c = g s c = 4.0 -4 kN/mm. Further, these values were verified by preliminary calculations so that the extent of damaged zones around the borehole is close to that observed in situ experiments [START_REF] Armand | Main outcomes from in situ thermo-hydro-mechanical experiments programme to demonstrate feasibility of radioactive high-level waste disposal in the callovo-oxfordian claystone[END_REF][START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of callovo-oxfordian claystone: from experimental studies to model calibration and validation[END_REF][START_REF] Conil | In situ investigation of the thm behavior of the callovo-oxfordian claystone[END_REF].

(a) The COx claystone is considered as a transversely isotropic material with different values of elastic modulus, permeability, and thermal conductivity in the parallel (horizontal) and perpendicular (vertical) directions to bedding planes. In Table 2, the initial values of the input parameters are presented. And we use the same set of plastic parameters as that used in the simulation of triaxial compression tests.

ALC cell r = 0.375m 25m 25m x y . 0  w n . 0  q n . 0  w n . 0  q n . 0  q n . 0  w n 0 y u  0 x u  12.7 y MPa    4.7 p MPa  20 T C   12.4 x MPa    (b)

Elastic parameters

Young's modulus

E =5GPa E ⊥ =3GPa Poisson's ratio ν =0.2 ν ⊥ =0.35 Hydraulic parameters Permeability k 0 =2 × 10 -20 m 2 k 0⊥ =1 × 10 -20 m 2 Biot coefficient b 0 =0.6 Porosity φ 0 =0.16
Thermal parameters Thermal conductivity About the strategy for considering THM and damage coupling, a sequenced explicit algorithm is adopted. It is assumed that the evolutions of fluid pressure and rock deformation do not significantly affect the heat diffusion problem. Therefore, at each time step, the thermal problem is first solved. The obtained values of temperature are used input for the hydraulic and mechanical problems. These two problems are then solved simultaneously by an iterative updating method. Finally, the two damage problems are solved to update the values of phase-field variables. Their values are fuzzed and used as input data for updating the THM parameters which are used in the resolution of the THM problems at the next time step.

λ =2W m -1 : K -1 λ ⊥ =1.33W m -1 : K -1
Based on the in situ drilling speed and on the preliminary analysis of excavation front effect, the excavation of the ALC1604 borehole is assumed to be mechanically and hydraulically completed during 24 hours on the first day. During this period, the stress traction and fluid pressure are linearly decreased from their initial values to zero (atmospheric pressure). The scenarios of the stress traction and fluid pressure reduction on the borehole wall are presented in Figure 7(a)(b). After the excavation, there is a waiting period before heating. During the waiting period and heating phases, the stress traction and fluid pressure are maintained to the atmospheric pressure. The heating starts on the 176th day after the excavation. The real heating power in the three-dimensional configuration is converted to the equivalent heating power for the 2D model. In practice, a 3D thermal simulation is first performed using the real power of the in situ heating experiment to determine the distribution of temperature in the whole domain. The distribution on the cross section selected for 2D thermo-hydromechanical (THM) simulations is also extracted. Then, a 2D thermal simulation is conducted on this section using an equivalent heating power, which is adjusted in order to obtain very similar distributions of temperature on the cross section between the 2D and 3D calculations. For the present study, the retained equivalent heating power is 220/3 W/m. Based on preliminary calculations on influences of time increment size, in particular on evolutions of fluid pressure, the following time increments are used in the THM simulations: 0-1 day: 1 hour (excavation) 1-176 days: 1 day (waiting period) 176-186 days: 1 hour (heating starts at day 176; small time increment is used as the temperature increases quickly during this period) 186-2160 days: 1 day (larger time increment is adopted as the temperature increase rate becomes smaller) (a)

Main numerical results

In order to compare numerical results with in-situ measurements, the THM responses at selected measuring points are presented and discussed. At first, the sensor set of the ALC4003 borehole is considered. The evolutions of temperature at three selected points are presented and compared with the measured values in Figure 8. It is found that the temperature increases quickly during the first period of heating and progressively evolves towards a stationary value. At the sensor 03, there is a good agreement between numerical prediction and in situ measurement. However, at the sensors 01 and 02, relatively large differences are observed. It seems that the numerical results under-estimate during the first heating period and over-estimate during the second heating period the temperature increase. These differences can be due to multiple reasons. For instance, as the real 3D heating problem is here simplified to the 2D one, the heat diffusion along the third axis (along the heating cell) is neglected. On the other hand, according to some previous studies, an improved thermal model should be used to better describe thermo-mechanical coupling by considering local thermal non-equilibrium [START_REF] Gelet | A thermo-hydro-mechanical coupled model in local thermal nonequilibrium for fractured hdr reservoir with double porosity[END_REF]. In Figures 9 and 10, one shows the numerical results of pore pressure and the comparisons with the experimental data measured by the sensors ALC1616-02 and ALC1617-02. These two sensors are respectively placed on the horizontal and vertical axes of the ALC borehole. It is seen that there are two distinct phases for the fluid pressure evolution. The fluid pressure evolution due to the excavation process between the instances A and B, and that during the heating phase after the instance B. In general, the excavation can induce an increase (or decrease) of pore pressure due to a compressive (or dilatant) volumetric deformation of the claystone. Then this pressure perturbation disappears progressively by the fluid diffusion in the rock. When the heating starts, there is a sharp increase of pore pressure mainly due to the differential thermal dilation coefficients between the claystone and water. The pore pressure progressively decreases with time due to the hydraulic flow process even if the heating continues. These basic trends are correctly reproduced by the numerical modeling. More precisely, the distribution of pore pressure at the end of excavation is shown in Figure 11. It is first found that in spite of the symmetrical mechanical loading conditions (the two principal stresses in the plane are quasi identical), the pore pressure exhibits a clearly dissymmetrical distribution. This is due to the multiple factors including the anisotropic property of elastic modulus and permeability and the dissymmetry of excavation induced damage zone. An overpressure zone appears in the horizontal direction at some distance from the alveole while the pore pressure along the vertical axis is systematically reduced. The ALC1616-02 sensor is inside the overpressure zone. The value of pore pressure at the position of this sensor increases during the excavation. The numerical result captures this pressure increase tendency. But the amplitude of the increase obtained is much lower than that measured by the sensor. Several reasons can be mentioned. For instance, the volumetric compaction due to both elastic and plastic strains induced by excavation is underestimated. The local Biot's modulus of COx claystone is underestimated. Further investigations are needed to understand the reasons for this scatter and to improve the numerical prediction. On the other hand, at the position of the ALC1617-02 sensor, the pore pressure is decreased by the excavation. Although the numerical simulation reproduces this trend but clearly under-estimates the amplitude of decrease. Again, as for the point ALC1616-02, it is needed to perform further investigations to improve the numerical predictions. For instance, fully three-dimensional simulations can be realised to identify possible geometrical effects on the pore pressure evolutions. As an advantage of the proposed double phase-field method, the distributions of both tensile and shear damages can be obtained and presented in Figure 12 at the end of excavation. Again, due to the structural anisotropy of elastic properties and permeability of the COx claystone, though the initial stresses are quasi isotropic, one gets anisotropic damaged zones between the horizontal and vertical directions. It is found that the damaged zone is clearly localized along the horizontal axis. Further, the excavation induced damaged zone is dominated by the tensile cracks. This result on the damaged zone distribution seems in good concordance with the in-situ observation of cracked zones around the alveole ALC3005, as shown in Figure 13. One can see that tensile cracks are concentrated in the horizontal direction. However, in this study, only concentrated damaged zones are obtained. Localized cracks as those observed in-situ are not reproduced by the numerical modeling. Some important improvements should be introduced. One idea is to consider the spatial variability of mineralogical compositions, and then of the mechanical properties, of the COx claystone in the numerical modeling. In this way, the initiation and localization of damage in some weak zones can be enhanced to capture the onset and propagation of localized cracks. On the other hand, improved return-mapping algorithms of the plastic and/or damage function and time stepping procedures can be envisaged to enhance the numerical localization of plastic strain and damage, as discussed in some previous works as explained in [START_REF] Heinze | Modeling porousrockfracturinginducedby fluid injection[END_REF][START_REF] Heinze | Systematic study of the effects of mass and time scaling techniques applied in numerical rock mechanics simulations[END_REF]. [68]

As an example of the evolutions of pore pressure and induced damages during the heating phase, their distributions on 1190 days after heating are presented in Figure 14. One can see that due to the adiabatic and impermeable conditions used on the external boundaries, the heat conduction generates an increase of temperature which in its turn induces a rise of pore pressure, even in the zone far from the heating cell, for instance at 4m to 6m. Further, due to the anisotropy of heat conductivity and permeability, the pressure distribution remains anisotropic. The value of overpressure is higher in the vertical direction than in the horizontal one (note that the horizontal permeability is higher than the vertical one). On the other hand, the damage process is continuously dominated by the tensile cracks, as presented in Figure 15. It is seen that the extents of strong damaged zones (say the value of damage variable is higher than 0.8) remain nearly unchanged with respect to those at the end of excavation. Only the extents of weak damaged zones exhibit a small extension of about 20% due to the heating process. Further, even if the amplitude of the shear damage is smaller than that of the tensile one, its extent is slightly larger than that of the tensile one. However, we should specify that the free displacement condition is assumed at the cell wall, which may lead to an overestimation of damaged zones. In the case when the steel casing is included in the modeling, the heating induced additional evolution of damages zones may not occur. As for the simulation of conventional triaxial compression tests, the mesh sensitivity of numerical results is also verified for the numerical study of the ALC1604 heating experiment. For this purpose, a second mesh (called mesh 2) is used and composed of 4560 elements and 4704 nodes. The comparisons of numerical results between the two meshes are presented in Figures 8 and9, respectively for the evolutions of temperature and pore pressure. In Figures 17 and18, the comparisons of tensile and shear damages along the horizontal and vertical axes are presented. One can observe that the differences between the two meshes are very small. It seems that the numerical results obtained with the phase-field model show a small mesh sensitivity.

Conclusion

In this study, we have developed a new phase-field method for modeling both the tensile and shear cracks in rock-like materials. The proposed model also takes into account the coupling between crack evolution and plastic deformation. The proposed phase-field method is further coupled with THM problems. All coupled THM and phase-field problems are solved by using the standard finite element method framework. The proposed method is first applied to modeling the mechanical behavior of the COx claystone in triaxial compression tests. It is found that the phase-field method is able to describe the transition from diffuse damage to localized cracking as well as the macroscopic mechanical responses of the claystone in both the pre-and post-peak regimes. The proposed phase-field method is then coupled with THM modeling and applied to study the THM responses during an in-situ heating experiment. In a general way, the numerical results well reproduce the main trends of in-situ measurements in terms of the evolutions of temperature, pore pressure and induced damaged zones. However, some quantitative differences exist on the evolutions of temperature and pore pressure. It is found that the initial anisotropy of the claystone in terms of elastic properties, permeability and heat conductivity plays a key role in the THM responses. The damaged zones due to excavation and induced by heating process are mainly controlled by tensile cracks. The amplitude of damage and the extent of damaged zones are moderately increased by the heating process. Several important improvements can be envisaged in future works. For instance, improved return-mapping algorithms of the plastic and/or damage function and time stepping method should be considered to better capture the localization of plastic strain and damage. Natural spatial variability of mineral compositions and then of mechanical and physical properties of the COx claystone should also be taken into account. Three-dimensional simulations are also desirable in order to avoid errors induced by twodimensional simplifications.

Acknowledgement:

DECOVALEX is an international research project comprising participants from industry, government and academia, focusing on development of understanding, models and codes in complex coupled problems in sub-surface geological and engineering applications; DECOVALEX-2019 is the current phase of the project. The authors appreciate and thank the DECOVALEX-2019 Funding Organizations Andra, BGR/UFZ, CNSC, US DOE, ENSI, JAEA, IRSN, KAERI, NWMO, RWM, S ÚRAO, SSM and Taipower for their financial and technical support of the work described in this report. The statements made in the report are, however, solely those of the authors and do not necessarily reflect those of the Funding Organizations.

Appendix:

The boundary value problem for the tensile phase field is defined as follows:

       h 1 H sph + h 2 H t dev + g t c ( 1 l d d t -l 2 d ∆d t ) = 0 in Ω ḋt ≥ 0 in Ω d t = 1 on Γ t ∇d t • n = 0 on ∂Ω (41)
The boundary value problem for the shear phase field is defined as follows:

      

Figure 1 :

 1 Figure 1: Conventional triaxial compression tests: (a) Geometry and boundary conditions of samples; (b) Location of weak element

Figure 2 :

 2 Figure 2: Axial and lateral strains versus differential stress in conventional triaxial compression tests with different confining pressures (2, 6, 12 MPa): comparisons between numerical results (continuous line) with consideration of plastic deformation and the experimental data (dotted lines with triangles) reported in[START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of callovo-oxfordian claystone: from experimental studies to model calibration and validation[END_REF]

Figure 3 :- 1 Figure 4 :

 314 Figure 3: Distributions of tensile and shear crack fields simulated with consideration of plastic deformation in the conventional triaxial compression test with Pc=12 MPa, at four levels of axial strain: (a) ε a = 2.7%; (b) ε a = 2.8%; (c) ε a = 2.85%; (d) ε a = 2.9%

Figure 5 :

 5 Figure 5: General view of the ALC experiment

Figure 6 :

 6 Figure 6: (a) Geometrical domain and boundary conditions for 2D plane strain modeling; (b) Illustration of reference mesh (called mesh1)

Figure 7 :

 7 Figure 7: Reduction of the (a) stress traction coefficient and (b) the pore pressure at the cell wall as a function of time.

Figure 8 :

 8 Figure 8: (a) Locations of the sensors 01-03 around the ALC4003; Temperature evolutions at (b) sensor 01: d=1.06m; (c) sensor 02: d=1.98m; (d) sensor 03: d=2.98m.

Figure 9 :Figure 10 :

 910 Figure 9: Pore pressure evolution at the point 1616-02 in ALC1616 (horizontal)

Figure 11 :

 11 Figure 11: Distribution of pore pressure at t=1 day (at the end of excavation).

Figure 12 :

 12 Figure 12: Distributions of (a) tensile damage and (b) shear damage at t=1 day (end of excavation).

Figure 13 :

 13 Figure 13: Picture of cracked zones around the alveole ALC3005.

Figure 14 :

 14 Figure 14: Distribution of pore pressure after t=1190 day of heating

Figure 15 :Figure 16 :

 1516 Figure 15: Distribution of (a) tensile damage and (b) shear damage at t=1190 days after heating

of mesh 1 d c of mesh 1 d t of mesh 2 d c of mesh 2

 2

Figure 17 :of mesh 1 d c of mesh 1 d t of mesh 2 d c of mesh 2

 172 Figure 17: Comparisons of tensile and shear damages along the horizontal axis between two meshes at t=1190 days after heating

Figure 18 :

 18 Figure 18: Comparisons of tensile and shear damages along the vertical axis between two meshes at t=1190 days after heating

h 3 H s dev + g s c ( 1 l d d s -l 2 d

 2 ∆d s ) = 0 in Ω ḋs ≥ 0 in Ω d s = 1 on Γ s ∇d s • n = 0 on ∂Ω (42)References:

Table 1 :

 1 . Typical values of mechanical parameters for COx claystone and numerical simulation input data

	Elastic parameters:
	E(GPa)	4
	ν	0.2
	Plastic parameters:
	C(MPa)	8.8
	B	1.2×10 -4
	α 0	0.25
	α m	1.1
	η c	0.5
	Crack fields parameters:
	g t c (kN/mm)	5.1 × 10 -4
	g s c (kN/mm)	g s c = g s c0 e bc(σc/σr) ; g s c0 = 4.25 × 10 -4 kN/mm, b c = 0.0977
	l d (mm)	0.37
	Axial displacement increment:
	∆u(mm)	2 × 10 -3

Table 2 :

 2 Representative THM parameters of COx claystone used in ALC modeling.