Zhan Yu

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

Yue Sun

Meng Wang

Minh-Ngoc Vu

Carlos Plua

Numerical analysis of hydro-thermal fracturing in saturated rocks by considering material anisotropy and micro-structural heterogeneity

Keywords: Thermal cracking, material anisotropy, heterogeneous rocks, claystone, phase-field method, radioactive waste

This study is devoted to numerical modeling of cracking process induced by temperature change in saturated porous rocks in the context of geological disposal of radioactive waste. Effects of material anisotropy and heterogeneity are taken into account. The macroscopic elastic properties are determined from two steps of homogenization by considering pores and mineral inclusions at two different scales. An extended phase-field model is proposed to describe the initiation and propagation of localized cracks. Two damage variables are introduced to conveniently represent both tensile and shear cracks. New damage evolution criteria are defined by incorporating the pore pressure effect. Three application examples are presented. By assuming a random distribution of pores and inclusions, the efficiency of the proposed model for capturing the progressive cracking process is first verified in a triaxial compression test. The thermal cracking process in an anisotropic and heterogeneous sample is then investigated. The respective influences of elastic anisotropy and spatial variability of pores and inclusions are outlined. Finally, the proposed model is applied to a series of real laboratory thermal cracking tests. Both hydromechanical responses and cracking evolution patterns are investigated. Numerical results are compared with experimental measurements.

The main mechanisms involved in the thermal cracking process are highlighted.

Introduction

Geological disposal is considered as a possible solution of radioactive waste management in many countries. In this context, it is primordial to investigate short and long term thermo-hydromechanical response of geological and engineered barriers. 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]. During several decades now, the Andra has coordinated a series of research programs for the investigation of thermo-hydromechanical properties of the COx claystone, including laboratory and in situ experiments, theoretical and numerical modeling [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].

Among various important issues to be investigated, the understanding and characterization of cracking processes represents a particular interest. Indeed, initial and induced cracks can significantly affect the transport properties of host rock and then the confining performance of geological barrier. Cracks can be generated by different factors, disturbances induced by gallery or borehole excavation, and pore over-pressure due to temperature rise.

The thermal cracking is particularly concerned by the radioactive waste which generates heat power during a long period. Indeed, the heat emitted by HLW packages causes temperature rise and pore pressure increase in the surrounding saturated porous rock as the thermal expansion coefficient of pore water (2.3 ×10 -4 K -1 at 20 • C and 7.2 ×10 -4 K -1 at 90 • C) is higher than that of solid skeleton (1.5 ×10 -5 K -1 for the COx claystone). The excessive pore pressure increase can lead to rock cracking [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]. At the same time, the basic mechanical properties of clayey rocks such as the COx claystone can also be influenced by temperature change [START_REF] Menaceur | The thermo-mechanical behaviour of the Callovo-Oxfordian claystone[END_REF][START_REF] Liu | Mechanical behavior of claystone in lateral decompression test and thermal effects[END_REF]. Moreover, in order to study the thermo-hydromechanical (THM) behavior of the COx claystone in quasi real conditions and to assess different aspects related to the conception and design optimization of disposal facilities, several heating experiments (among others) from small-scale to full-scale have been conducted in the Underground Research Laboratory (URL) of Andra [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]. In particular, the monitored THM responses obtained in the full scale experiment, called ALC1604 [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], have been the main object of numerical modeling in the Task E of the DECOVALEX2019 projet [START_REF] Seyedi | Upscaling THM modelling from small-scale to full-scale in-situ experiment in the Callovo-Oxfordian claystone[END_REF]. The participating teams have used different numerical methods and constitutive models. For more details, the readers can refer to the special issue devoted to this topic [START_REF] Birkholzer | DECOVALEX-2019: An international collaboration for advancing the understanding and modeling of coupled thermo-hydro-mechanical-chemical (THMC) processes in geological systems[END_REF]. The phase-field method [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture Computer Methods in[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] has also been used in that project [START_REF] Yu | Numerical study of thermo-hydro-mechanical responses of in situ heating test with phase-field model[END_REF]. By solving a boundary value problem for crack field, this method is able to describe continuously the transition from diffuse damage to localized cracks in three-dimensional configurations. It has been applied to different kinds of materials and problems, such as 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], and plastic damage coupling [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][START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow Computer[END_REF]. However, in most previous studies, homogeneous rock materials have been considered.

The COx claystone is characterized by a multi-scale heterogeneity [START_REF] Robinet | Effects of mineral distribution at mesoscopic scale on solute diffusion in a clay-rich rock: Example of the Callovo-Oxfordian mudstone (Bure, France)[END_REF][START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: from experimental studies to model calibration and validation Computers and Geotechnics[END_REF]. For the sake of simplicity, two main families of materials heterogeneity should be taken into account: interparticle pores inside the sub-microscopic clay aggregate, and mineral inclusions, mainly quartz and calcite in the microscopic porous clay matrix. It is known that the mechanical properties of the COx claystone can be affected by the porosity and mineral contents of quartz and calcite [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: from experimental studies to model calibration and validation Computers and Geotechnics[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]. The cracking process, in particular, the onset and localization of cracks, can also depend on the micro-structural parameters. At the same time, the COx claystone also exhibits a transversely isotropic behavior. The material anisotropy affects also the cracking process [START_REF] Liu | Effects of deviatoric stress and structural anisotropy on compressive creep behavior of a clayey rock Applied Clay[END_REF][START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: from experimental studies to model calibration and validation Computers and Geotechnics[END_REF][START_REF] Zhao | A constitutive model for anisotropic clay-rich rocks considering micro-structural composition[END_REF].

The objective of this study is to develop an extended and improved phase-field model to describe hydro-thermal damage and cracking in saturated COx claystone, by considering material anisotropy and heterogeneity. For this purpose, a new theoretical framework is set up for the evolution description of tensile and shear crack fields in coupled thermohydromechanical conditions. The effects of pore pressure and temperature changes are properly taken into account. Further, the material heterogeneity is considered by using an anisotropic elastic model issued from two steps of homogenization by successively involving pores and mineral inclusions at two different scales. The efficiency of the proposed model is assessed by numerical examples and laboratory mechanical and thermal cracking tests.

Phase-field method for thermo-hydromechanical problems

We consider here a saturated porous medium that occupies the volume Ω with the external boundary ∂Ω. It is subjected to the body force f b in Ω, the surface force t N on the part of external boundary ∂Ω f and the prescribed displacement u on the complementary part external boundary ∂Ω u . Meanwhile, the porous medium is also subjected to the fluid flux ω on the external boundary ∂Ω ω and to the prescribed fluid pressure p on the external boundary ∂Ω p . At last, the temperature is assumed to be the same for the fluid and for the rock skeleton. It has the heat flux q on the external boundary ∂Ω q and the prescribed temperature change θ on the external boundary ∂Ω θ . The problem to be solved is to determine the displacement (strain and stress), pore fluid pressure and temperature fields, as well as the process of initiation and propagation of cracks within Ω during the entire loading history.

Regularized crack topology

The phase-field method for quasi-static brittle fracture finds its theoretical background from the variational principle proposed in [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. The basic idea of phase-field method is to approximate the sharp crack topology by the regularized smeared crack topology [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] by introducing a scalar-valued auxiliary variable d(x), which constitutes the so-called phase (or damage ) field. This variable takes the unit value on the crack surface and progressively vanishes away from it. The phase-field can be determined by solving its own boundary value problem. Further, the differential equations governing the boundary value problem are the Euler equations of the variational problem consisting to minimize the total crack density (surface) functional in the cracked body. The phase-field d(x) can also be seen as equivalent to the scalar damage variable used in the continuum damage mechanics.

Complex cracking modes are generally observed in rock-like materials under compressiondominating stresses [START_REF] Wong | Micromechanics of faulting in westerly granite[END_REF][START_REF] Wong | Analysis of crack coalescence in rock-like materials containing three flaws-part I: experimental approach[END_REF]. In order to easily deal with such mixed-mode cracks, in this work, two phase fields are introduced to represent the two basic families of cracks, i.e. tensile and shear cracks (d t for tensile crack and d s for shear crack). Each phase field is described by the following boundary problem [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]:

         d α -l 2 d △d α = 0 in Ω d α (x) = 1 on Γ ∇d α (x) • n = 0 on ∂Ω (1)
with α = t, s. l d is a length scale parameter controlling the width of smeared cracks.

According to [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], the above boundary value problem is equivalent to the following variational one:

d α (x) = Arg{ inf d∈S d Γ l d (d α )}; α = t, s (2)
with S d = {d α |d α (x) = 1 at x ∈ Γ}. Γ l d (d) denotes the total crack surface area, defined by:

Γ l d = Γ t l d (d t ) + Γ s l d (d s ) = Ω {γ t (d t , ∇d t ) + γ s (d s , ∇d s )}dV (3)
Γ α l d (α = t, s) represent the contributions of tensile and shear cracks to the total crack area, respectively. The tensile and shear crack density (surface area per unit volume) γ α (d α , ∇d α ) are defined as:

γ α (d α , ∇d α ) = 1 2 { 1 l d (d α ) 2 + l d ∇d α .∇d α } ; α = t, s (4)
In this way, the energy density per unit volume requested to create the crack is composted by tensile crack energy and shear crack energy:

w c (d t , d s , ∇d t , ∇d s ) = g t c γ t (d t , ∇d t ) + g s c γ s (d s , ∇d s ) (5)
with the material toughness for the tensile and shear cracks g t c and g s c . Further details of this double-phase field method can be found in [START_REF] Yu | Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling Computers and Geotechnics[END_REF].

Variational framework

According to the variational principle proposed in [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], the crack propagation is governed by the minimization of an energy functional, which is composed of the stored (or available) energy and that used to create the crack. This principle is here extended to saturated porous materials subjected to thermal-hydromechanical loads. The total energy functional is expressed in the following general form:

E(ε e , m, θ, d t , d s ) = Ω w e (ε e , m, θ, d t , d s )dΩ + Ω w c (d t , d s , ∇d t , ∇d s )dΩ (6)
w e denotes the elastic strain energy of cracked material, which is a function of elastic strain tensor ε e , fluid mass change per unit initial volume m and variation of temperature θ, as well as two damage variables.

Elastic free energy

For an undamaged saturated porous medium, the constitutive relations due to the thermo-poroelastic theory [START_REF] Olivier | Poromechanics[END_REF][START_REF] Cheng | [END_REF] can be expressed as: By using these constitutive relations, the elastic free energy of an undamaged saturated porous medium can be expressed in the following form:

σ -σ 0 = C 0 b : ε e -B(p -p 0 ) -A b θ (7)
p -p 0 = M (-B : ε e + m ρ 0 f ) + 3α m M θ (8)
s -s 0 = s 0 m m + α b : (σ -σ 0 ) -(3α m -B : α b )(p -p 0 ) + C b σ T 0 θ (9)
w 0 e (ε e , p, θ) = 1 2 σ b : ε e + 1 2 (p -p 0 ) 2 M - 1 2 C b σ T 0 θ 2 (10)
σ b denotes the Biot (elastic) effective stress tensor, which is defined by:

σ b = (σ -σ 0 ) + B(p -p 0 ) (11)
In order to better describe two cracking processes, some physical mechanisms are here considered. In most rocks, the evolution of tensile crack is generally induced by tensile stress while that of shear crack is driven by shear or deviatoric stresses. For this reason, the Biot effective stress tensor is decomposed into a positive (tensile) part and a negative (compression) part such as σ b = σ b+ +σ b-. Accordingly, the elastic strain energy is rewritten as follows:

w 0 e (ε e , p, θ) = w 0+ e (ε e ) + w 0- e (ε e ) + 1 2 (p -p 0 ) 2 M - 1 2 C b σ T 0 θ 2 (12)
with

   w 0+ e = 1 2 σ b+ : ε e w 0- e = 1 2 σ b-: ε e (13)
The decomposition of Biot effective stress tensor is due to the operators P ± σ : (a detailed description is given in [START_REF] Yu | Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling Computers and Geotechnics[END_REF])

   σ b+ = P + σ : σ b σ b-= P - σ : σ b (14)
As for the previous study [START_REF] Yu | Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling Computers and Geotechnics[END_REF], it is assumed that the tensile damage affects the positive stresses-related part while the shear damage influences the negative stresses-related part of the elastic strain energy. The effects of induced damages on the rest of poroelastic properties are neglected, such as Biot coefficient b and modulus M as well as the specific heat C b σ . In this way, the equation ( 12) can be rewritten as:

w e (ε e , p, θ, d t , d s ) = h t (d t )w 0+ e + h s (d s )w 0- e + 1 2 (p -p 0 ) 2 M - 1 2 C b σ T 0 θ 2 (15)
by using a common form of degradation function:

h α (d α ) = (1 -k)(1 -d α ) 2 + k ; α = t, s (16)
where k is a small positive value to avoid numerical instability after the material is fully broken (d α = 1).

Evolution of crack fields

In this subsection, we present the main equations governing the evolution of two crack fields. A detailed description can be found in the previous study [START_REF] Yu | Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling Computers and Geotechnics[END_REF]. By solving the minimization problem [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF][START_REF] Wu | Chapiter one -Phase field modeling of fracture Advances in Applied Mechancis[END_REF] of E given in [START_REF] Liu | Mechanical behavior of claystone in lateral decompression test and thermal effects[END_REF], the governing equations for two crack fields (α = t, s) can be written as:

                       - ∂w e ∂d α -g α c δ d α γ α = 0 , ḋα > 0 , in Ω - ∂w e ∂d α -g α c δ d α γ α ≤ 0 , ḋα = 0 , in Ω ∂γ ∂∇d α .n = 0 , on δΩ (17)
The evolution of two crack fields are described by the following equations:

         -h ′ t (d t )w 0 e+ -g t c 1 l d d t -l d div(∇d t ) = 0 , ḋt > 0 -h ′ s (d s )w 0 e--g s c 1 l d d s -l d div(∇d s ) = 0 , ḋs > 0 (18)
According to these equations, the tensile crack evolution is driven by the elastic strain energy w 0 e+ which is a function of tensile (positive) Biot effective stresses σ b+ (see Equations ( 13)). At the same time, the evolution of shear crack should be driven by the elastic strain energy w 0 e-which is attributed to compressive (negative) Biot effective stresses σ b-. However, the result due to the variational approach cannot correctly reflect the shear cracking mechanism in rocks. Indeed, it is known that the shear cracking is mainly driven by the maximum shear stress and compressive mean stress. Based on this physical mechanism and by extending the classical Mohr-Coulomb criterion [START_REF] Yu | Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling Computers and Geotechnics[END_REF], an alternative driving force w s -is proposed and given as follows:

w s -= 1 2G ⟨ ⟨σ t 1 ⟩ --⟨σ t 3 ⟩ - 2cosφ + ⟨σ t 1 ⟩ -+ ⟨σ t 3 ⟩ - 2 tanφ -c⟩ 2 + ( 19
)
with the bracket ⟨.⟩ ± such as:

   ⟨a⟩ + = a, ⟨a⟩ -= 0, a ≥ 0 ⟨a⟩ + = 0, ⟨a⟩ -= a, a < 0 ( 20
)
σ t 1 and σ t 3 are the major and minor Terzaghi effective principal stress [START_REF] Debuhan | On the validity of the effective stress concept for assessing the strength of saturated porous materials: a homogenization approach[END_REF][START_REF] Lydzba | Stress equivalence principle for saturated porous media Comptes Rendus Mecanique[END_REF]. c and φ denote the cohesion and friction angle of closed cracks. By substituting w s -for Equations ( 18), the evolution criteria for two crack fields are now expressed as:

             -h ′ t (d t )w 0 e+ -g t c 1 l d d t -l d div(∇d t ) = 0 , ḋt > 0 -h ′ s (d s )w s --g s c 1 l d d s -l d div(∇d s ) = 0 , ḋs > 0 (21)
Based on the concept of [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations[END_REF] for the irreversible process of crack evolution, two energy history functionals are defined as follows:

   H t (t) = max[w 0 e+ (τ )] τ ∈[0,t] H s -(t) = max[w s -(τ )] τ ∈[0,t] (22)
Using these energy history functionals, Equations ( 21) can be rewritten as:

             -h ′ t (d t )H t -g t c 1 l d d t -l d div(∇d t ) = 0 -h ′ s (d s )H s --g s c 1 l d d s -l d div(∇d s ) = 0 ( 23
)

Numerical implementation in finite element method

In this thermo-hydromechanical problem considering tensile and shear crack, five coupled physical unknown fields should be determined by using the framework of finite element method. The weak forms of static equilibrium, hydraulic diffusion and heat diffusion are expressed as:

Ω δε : C b (d t , d s ) : εdV - Ω δε : (δp)IdV - Ω δε : (3α b K b (d t , d s )T )IdV = Ω f t.δudS (24) Ω k(d t ) µ ∇p.∇(δp)dV = Ωω k(d t ) µ δp.∇pndS - Ω 1 M ∂p ∂t δpdV - Ω b ∂ε kk ∂t δpdV + Ω 3α m ∂θ ∂t δpdV (25)
Ω λ∇T.∇(δθ)dV = Ωq λδT ∇T.ndS - Ω C b ε ∂θ ∂t δθdV - Ω (3α b K b (d t , d s )T 0 ) ∂ε kk ∂t δθdV + Ω (3α m T 0 ) ∂p ∂t δθdV + Ω k(d t ) µ ∇p.(C p ∇θ)δθdV (26)
Following Equations ( 23), the weak forms of governing equations of the tensile and shear crack fields can be written as:

Ω {(g t c /l d + 2H t )d t δd t + g t c ∇d t • ∇(∂d t )}dV = Ω 2H t ∂d t dV (27) Ω {(g s c /l d + 2H s -)d s δd s + g s c ∇d s • ∇(∂d s )}dV = Ω 2H s -∂d s dV (28)
From these weak forms, the THM fields are obviously affected by the variation of

C b (d t , d s ), k(d t ) and K b (d t , d s ):
The drained elastic stiffness tensor of damaged material C b is given by:

C b (d t , d s ) = h t (d t )P + σ + h s (d s )P - σ : C 0 b ( 29
)
with C 0 b being the elastic stiffness tensor of undamaged material. The permeability tensor of cracked porous medium k(d t ) is defined as:

k(d t ) = k 0 exp(β k d t ) ( 30
)
with the initial permeability tensor of intact porous medium k 0 and the parameter β k controlling its evolution. The drained bulk modulus K b (d t , d s ) can be calculated according to

C b (d t , d s ).
For the sake of simplicity, the rest of the material parameters are not affected by the damage variables such as: the dynamic viscosity µ and the tensor of thermal conductivity λ.

Based on elementary approximation with suitable shape functions, the following systems of equations to be solved are obtained:

         R uu ∆U + C up ∆P + C uT ∆T = ∆F e C pu ∆U + (∆tR pp + M pp )∆P + C pT ∆T = ∆t(-R pp P + ∆F ω ) C T u ∆U + C T p ∆P + (∆tR T T + M T T )∆T = ∆t(-R T T T + ∆F q ) (31)
The increment nodal values of displacement△U, pore pressure△P and temperature△T can be determined by a fully coupled method at each time step. The detailed expression of all matrices can be found in the previous study [START_REF] Yu | Numerical study of thermo-hydro-mechanical responses of in situ heating test with phase-field model[END_REF].

On the other hand, the discrete system of equations for the tensile crack can be written as:

         K d t d t = F d t K d t = Ω {(g t c /l d + 2H t )N T p N p + g t c l d B T p B p }dV F d t = Ω 2H t N T p dV (32)
and for the shear crack:

         K d s d s = F d s K d s = Ω {(g s c /l d + 2H s -)N T p N p + g s c l d B T p B p }dV F d s = Ω 2H s -N T p dV (33)
As mentioned above, the damage effects are taken into account in this THM problem. The five fields are inherently coupled.

In this work, an explicit decoupled algorithm is chosen to solve the THM coupled problems with relatively small time increment. Moreover, the so-called Alternate Minimization (AM) solver proposed in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF] is used to couple the displacement field and phase fields because of its good numerical robustness.

Consideration of material heterogeneity

In most rock-like materials, the initiation or nucleation of new cracks is strongly motivated by the material heterogeneity. For instance, cracks can initiate at interfaces between stiff inclusions and weak matrix phase, or from some weak zones with high porosity. In a number of previous studies, a numerical trick is generally used. Weak elements are arbitrarily placed inside samples or structures in order to facilitate the nucleation and guide the propagation of cracks. In order to more physically reflect the effect of material heterogeneity on macroscopic mechanical properties, a micro-mechanics based approach is adopted in this work. More precisely, by using an analytical homogenization method, macroscopic elastic properties are determined as explicit functions of porosity and inclusion volume fraction at the microscopic scale.

As an example, we consider here the Callovo-Oxfordian (COx) claystone which has been widely investigated in France in the context of geological disposal of radioactive waste.

According to the previous studies [START_REF] Robinet | Effects of mineral distribution at mesoscopic scale on solute diffusion in a clay-rich rock: Example of the Callovo-Oxfordian mudstone (Bure, France)[END_REF], two representative material scales can be considered.

Mineral particles, mainly quartz, calcite and pyrite, are randomly distributed in a quasi continuous clay matrix at the mesoscopic scale (hundreds of micrometers). The majority of pores are found inside the clay matrix at the microscopic scale (less than micrometer).

Therefor, two micro-structural parameters are adopted here, the porosity of the clay matrix f p and the volume fraction of mineral inclusions f i which are defined as:

f p = Ω p Ω p + Ω m ( 34
)
f i = Ω in Ω i + Ω p + Ω m (35)
Ω is the volume of representative volume element, while Ω p , Ω i and Ω m are respectively the volumes occupied by pores, inclusions and solid clay matrix.

Effective elastic properties

The macroscopic elastic properties are determined by performing two steps of linear homogenization. Considering the inclusion-matrix morphology, the standard Mori-Tanaka scheme is used for each step of homogenization [START_REF] Mori | Averages stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF]. At the first step, the effective elastic tensor of porous matrix C mp is determined by taking into account of porosity effect. One gets:

C mp = (1 -f p )C m : [(1 -f p )I + f p (I -P p : C m ) -1 ] -1 (36)
P p (C m ) is the fourth order Hill tensor for ellipsoidal pores [START_REF] Giraud | Application of results on Eshelby tensor to the determination of effective poroelastic properties of anisotropic rocks-like composites[END_REF][START_REF] Zhao | A constitutive model for anisotropic clay-rich rocks considering micro-structural composition[END_REF], C m is the elastic tensor of solid matrix, and I is the fourth order unit tensor. At the mesoscopic step of homogenization, the macroscopic elastic tensor of heterogeneous rocks is determined after including the influence of mineral inclusions:

C hom = C mp + [f i (C in -C mp ) : D in ] : [I + f i (C in -I)] -1 (37)
with

D in = [I + P i : (C in -C mp )] -1 (38)
P i (C mp ) is the Hill tensor for spherical inclusions [START_REF] Giraud | Application of results on Eshelby tensor to the determination of effective poroelastic properties of anisotropic rocks-like composites[END_REF][START_REF] Zhao | A constitutive model for anisotropic clay-rich rocks considering micro-structural composition[END_REF]. For the sake of simplicity, different families of mineral particles are merged into an equivalent inclusion phase, and its elastic stiffness tensor is denoted as C in .

Description of material heterogeneity

With the homogenization-based models, the macroscopic elastic properties are now functions of porosity and inclusion volume fraction at two different scales. Therefore, the spatial variability of macroscopic elastic properties of COx claystone are attributed to the nonuniform distribution of porosity and inclusion concentration. In some previous works, for instance [START_REF] Chu | Stochastic multiscale modeling with random fields of material properties defined on nonconvex domains[END_REF], random micro-structures were generated by using Monte-Carlo realizations and the elastic properties fields at the mesoscopic scale were determined by using a movingwindow homogenization method. A non-Gaussian distribution was proposed for the variability of elastic tensor. Similarly, we assume here a random distribution of porosity f p and volumetric fraction of mineral inclusions f i by using the standard Weibull distribution function:

φ α = m α β α ( f α β α ) mα-1 exp[-( f α β α ) m α ] , α = p, i (39)
where f α is the volume fraction of pores (α = p) and mineral inclusions (α = i). β α is the scale parameter of the distribution which represents the mean value of random variable. m α is the homogeneity index of the material.

Mechanical behavior of COx claystone

In this section, the proposed elastic phase-field model is applied to studying the mechanical behavior of the Callovo-Oxfordian claystone. The emphasis is put on the full mechanical response during a triaxial compression test by considering the nucleation, propagation and localization of cracks. For this purpose, the tested sample is seen as a small structure composed of heterogeneous material. The sample geometry and boundary conditions are presented in Figure 1. Two-dimensional plane strain configuration is adopted. The sample is divided into 20000 quadrilateral elements with 20301 nodes. The elements with 0.37 mm size are placed by a structured way. In this study, the samples are drilled in the perpendicular direction to the bedding planes.

The input parameters are mainly based on previous studies [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: from experimental studies to model calibration and validation Computers and Geotechnics[END_REF][START_REF] Zhao | A constitutive model for anisotropic clay-rich rocks considering micro-structural composition[END_REF]. A transversely isotropic elastic behavior is considered for both the solid clay matrix and COx claystone.

The reference values of elastic properties for the solid clay matrix are: E s ∥ = 3.6 GPa, [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: from experimental studies to model calibration and validation Computers and Geotechnics[END_REF]. For the sake of simplicity, the same value of homogeneity index m = 1.5 is adopted for both pores and inclusions. In this way, the Weibull distribution of inclusion volume fraction and porosity are calculated and illustrated in Figure 1. According to 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] Borden | A phase-field description of dynamic brittle fracture Computer Methods in[END_REF][START_REF] Wu | Chapiter one -Phase field modeling of fracture Advances in Applied Mechancis[END_REF], the scale length l d is conveniently taken as being 1 to 3 times of the smallest element size. In our study, it is equal to twice of the smallest element size: l d = 7.4 × 10 -4 m. In addition, the tensile material toughness g t c can be physically related to the uniaxial tensile strength of tensile crack σ t . The following relation is widely used:

E s ⊥ = 2.4 GPa, ν ∥ = 0.3, ν ∥⊥ = 0.
g t c = 256σ 2 t l d 27E ⊥ (40)
It is worth noticing that this relation was initially proposed for isotropic materials. It is here extended to transversely isotropic ones. For this purpose, the perpendicular elastic modulus E ⊥ is used. For the COx claystone, the uniaxial tensile strength is evaluated to be σ t = 3

MPa. This leads to g t c = 15 N/m. The shear damage evolution is controlled by the material toughness g s c , friction angle φ and cohesion c. Their values can be calibrated from the peak deviatoric stresses obtained in triaxial compression tests with different confining pressures.

Some details can be found in [START_REF] Yu | Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling Computers and Geotechnics[END_REF]. In this work, the following values are retained: g s c = 950 N/m, φ = 15 o and c = 0.1 MPa.

Parameters Values

Elastic parameters In Figure 3, we show the distributions of tensile and shear damage at four different loading stages during the triaxial compression test with a confining pressure of 2 MPa. It is observed that the shear damage evolves more quickly and localizes earlier than the tensile one. Therefore, shear cracking is the dominant process under the triaxial compression condition. But due to the coupling effect, the tensile damage also initiates and evolves. The cracking patterns (localized damage bands) obtained by the phase-field model are consistent with common experimental observations. It is worth noticing that, in the majority of previous studies, a weak element was inserted in the mesh to enhance the onset of damage localization. But in the present work, the process of localized damage is naturally due to the material heterogeneity, without the need for any numerical tips.

E s ∥ = 3.6 GPa, E s ⊥ = 2.4 GPa, ν ∥ = 0.3, ν ∥⊥ = 0.3, G ∥⊥ = 1.

Study of hydro-thermal fracturing

In this section, we shall present numerical modeling of cracking process induced by temperature change in saturated porous materials. A numerical test is first considered to show the coupling mechanisms between thermal-hydro-mechanical responses and cracking process. Then, a laboratory heating test is investigated.

A numerical test

We consider here a cylindrical sample of saturated porous medium with a radius of 10mm, as shown in Figure 4. The external boundary is first subjected to a confining stress of σ 0 rr = -12 MPa. An initial pore pressure of p 0 = 5 MPa is also prescribed. Then, a constant heat power of Q=0.44 W is applied at the cylinder center point for t ≥ 0. The external boundary surface is considered as impermeable adiabatic during the heating period. Due to the double symmetry of this plane strain problem, only a quarter of the cross section is considered, and it is meshed by 50632 quadrilateral elements with 50990 nodes. The elements with a average size of 0.4 mm are placed by a structured (regular) way. In Figure 5, we show the distributions of temperature, pore pressure and tensile damage at four different heating times respectively of 600, 880, 1080 and 1200 minutes. It is worth noticing that in this specific case, it is expected that the tensile damage should be the main cracking process due to the thermal induced pore pressure raise. Thus, only the tensile damage distribution is presented in this figure. One can see that due to the small size of the sample, the heat and liquid diffusion takes place quite rapidly in the whole sample.

The distributions of temperature and pore pressure become nearly uniform quite quickly.

However, it is very interesting to observe that the tensile damage field is clearly not uniform. This is the consequence of material heterogeneity. Moreover, at the first two heating steps, namely at 600 and 880 minutes, the thermal-induced pore pressure is clearly smaller than the prescribed confining stress of 12 MPa. It means that the average effective stress still remains in compression. But one gets already important distributions of tensile damage.

At some particular points, the damage variable value is even close to 1. This shows that due to material heterogeneity, there exist a number of weak zones in the sample, facilitating the appearance of cracks. Therefore, the material heterogeneity plays a primordial role in the cracking process. Finally, at the two last heating steps, multiple localized fractures are formed in the sample. The proposed phase-field model captures this process well.

Figure 6 shows the results at 1200 minutes of heating by considering an isotropic material.

As said before, the distributions of temperature and pore pressure are nearly uniform, and they are not clearly affected by the material anisotropy. But more interestingly, even the difference on the tensile damage field is also very small. This seems to indicate that the tensile damage evolution in this case is mainly influenced by the material heterogeneity rather than anisotropy.

In the present micro-mechanics based elastic model, the macroscopic elastic properties are affected by the porosity and inclusions at two different scales. In order to investigate the respective effects of these two kinds of micro-structural parameters, two comparative calculations are performed by respectively considering the porosity or inclusion fraction heterogeneity alone. The obtained results of tensile damage field at 1200 minutes of heating are presented in Figure 7. Significant differences are obtained. In the case of heterogeneous porosity but homogeneous inclusion fraction, one obtains a small number of but large localized fractures. Moreover, there are more cracks along the horizontal direction than the vertical one. It seems that in this situation, the material anisotropy effect is significant.

Inversely, for the case with homogeneous porosity but heterogeneous inclusion fraction, one gets a high number of but narrow fractures. There is no clear effect of material anisotropy.

Temperature

Pore pressure Tensile damage

Analysis of laboratory tests

In the previous studies [START_REF] Braun | Theoretical Analysis of Pore Pressure Diffusion in Some Basic Rock Mechanics Experiments[END_REF][START_REF] Braun | Thermo-hydro-mechanical behavior of the Callovo-Oxfordian claystone: Effects of stress paths and temperature changes[END_REF], the authors have performed a series of laboratory tests aiming at the characterization of thermal induced cracking in saturated COx claystone.

Those tests were conducted on cylindrical samples under undrained conditions. The geometry, initial and boundary conditions are presented in Figure 8. Three tests, denoted as EXT1, EXT2 and EXT3, are here considered. The initial values of radial, axial stresses (σ r 0 , σ a 0 ), pore pressure (p 0 ) and temperature (T 0 ) of these tests are given in Table 2. The elastic and damage parameters, the average porosity and inclusion volume fraction as well as the homogeneity index are the same as those given in the previous section on the triaxial compression test. The other input parameters as given in Table 3.

However, the volumetric thermal expansion coefficient of water α w is taken as a function of temperature by using the following empirical relation [START_REF] Zhutovsky | Evaluation of the Thermal Expansion Coefficient Using Non-Destructive Testing Conference paper[END_REF]: Permeability (perpendicular to bedding)

α w = (-
m 2 k v 1.33 × 10 -20
Thermal conductivity (parallel to bedding)

W • m -1 • K -1 λ h 1.95
Thermal conductivity (perpendicular to bedding)

W • m -1 • K -1 λ v 1.28
Heat capacity of rock

J • kg -1 • K -1 C p 800
Thermal expansion of rock

K -1 α b 1.5 × 10 -5
Bulk modulus of water Pa K w 2.2 × 10 9

Dynamic viscosity of water

Pa • s µ w 1 × 10 -3 Density of water kg • m -3 ρ w 1 × 10 3
Heat capacity of water

J • kg -1 • K -1 C w 4180
After setting the initial conditions given in Table 2, the samples were subjected to a heating process at the external boundary surfaces. The variations of temperature in three tests are presented in Figure 9. The evolution of axial stress on the upper boundary surface was also monitored during the heating, and the obtained results are provided in the same figure. In the present numerical study, the experimental values of temperature and axial stress are used as input data of boundary conditions.

The thermal cracking problem is here solved under the axi-symmetrical conditions. The half-cross section of sample is meshed by 2375 (95×25) quadrilateral elements and 2496 nodes, as shown in Figure 10. The elements with a size of 0.4 mm are created by a structured way. At the same time, by using the material heterogeneity coefficient given above, the distributions of porosity and inclusion volume fraction inside the meshed sample are calculated and presented also in Figure 10. The scale length l d is equal to the 1.5 time of the element size: l d = 0.6 mm. The material toughness parameters are calculated with Equation (40). In order to quantify the critical stresses driving the macroscopic thermal failure, the variation of Terzaghi effective stresses were given in the experimental data respectively on the upper and lateral surfaces of samples [START_REF] Braun | Theoretical Analysis of Pore Pressure Diffusion in Some Basic Rock Mechanics Experiments[END_REF][START_REF] Braun | Thermo-hydro-mechanical behavior of the Callovo-Oxfordian claystone: Effects of stress paths and temperature changes[END_REF]. The same values are here calculated from numerical results and compared with the experimental data in Figure 12. There is also a good qualitative agreement between numerical and experimental results despite some quantitative differences, in particular for the radial effective stress. In general, the radial effective stresses remain or increase in the compression zone. On the contrary, due to the significant rise of pore pressure, the initial compressive axial effective stresses progressively decrease and become tensile ones. It is interesting to observe that the experimental values of peak tensile effective axial stresses in three tests are quite close to each other and around 3 MPa. But the numerical results are generally higher than the experimental ones, and reach 4 MPa the tests in EXT1 and EXT3. This feature will be discussed later in the sensitivity study. In addition to the hydromechanical response, the proposed phase-field method is able to describe the onset and propagation of cracks in terms of damage localization. Only the tensile damage d t appears in the tests. For instance, the distributions of tensile damage in the test EXT3 are presented in Figure 13 for four different heating instances. The value d t = 1 denotes the totally broken areas, marked by the red color. It is found that due to the material heterogeneity, the damage localization starts from some weak zones. With the heating process, a main crack is generated and propagates in a sub-horizontal direction.

Some secondly cracks are also created. The numerical cracking patterns seem to be consistent with that observed in the tested sample as shown in Figure 13. (MPa) at t=9058 s for test EXT3

Sensitivity study

In the preceding section, the numerical results for a selected case are compared with the experimental data, demonstrating a high degree of consistency of the proposed phase field model. This section aims to further evaluate the sensitivity of the model by examining three aspects: i) the influence of randomly generated heterogeneous distributions of material parameters, ii) the impact of mesh size, and iii) the effect of tensile damage driving force.

In order to investigate the sensitivity of material heterogeneity, three additional simula- Furthermore, the mesh dependence is also studied by considering two additional finer meshes: 9500 (190×50) and 38000 (380×100) elements, respectively. For each mesh, three random distributions of material properties are used. The variations of the axial Terzaghi effective stress are shown for each of these two meshes in Figure 17. For the sake of readability, one of the three randoms generations is selected to illustrate the crack pattern. It is worth noticing that in the calculations using the two different meshes, the scale length l d is defined as 1.5 times the element size and the corresponding toughness parameter g t c (l d ) is then decreased according to Equation [START_REF] Chu | Stochastic multiscale modeling with random fields of material properties defined on nonconvex domains[END_REF]. This leads to a decrease in the thickness of obtained cracks. Given that the main cracks always appear quasi horizontally, they can be seen as tensile cracks. Further, the curves of effective stress versus temperature are almost identical for the different cases. It can be concluded that the random distribution of material heterogeneity does not significantly impact the macroscopic responses of the studied sample.

tions
Figure 18 presents a comparison of all nine calculations previously mentioned (3 generations x 3 mesh finesse levels). The maximum difference in the peak axial stress among these calculations is approximately 7%. This suggests that the mesh size and random distribution have little impact on the macroscopic response of the tested sample, and the numerical results can be considered as representative. As a summary of above comparisons, it is shown that the tensile damage is the principal mechanism in three thermal cracking tests. But the tensile damage is driven by the positive strain energy w 0+ e = 1 2 σ b+ : ε e defined in equations [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF], where σ b+ is the positive part of Biot effective stress. In order to investigate the influence of the tensile damage driving force on the hydromechanical responses and cracking process, four different kinds of forms are here considered:

1. case 1, w 0+ e = 1 2 σ b+ : ε e , σ b = (σ -σ 0 ) + B(p -p 0 ) 2. case 2, w 0+ e = 1 2 σ t+ : ε e , σ t = (σ -σ 0 ) + I(p -p 0 ) 3. case 3, w 0+ e = 1 2 σ bθ : ε e , σ bθ = (σ -σ 0 ) + B(p -p 0 ) + A b θ 4. case 4, w 0+ e = 1 2 σ tθ : ε e , σ tθ = (σ -σ 0 ) + I(p -p 0 ) + A b θ
The case 1 is the reference one presented above. Three additional series of calculations are performed by using the other forms of driving energy for the tensile damage. As the most relevant comparisons between the four cases, the variations of axial Terzaghi effective stress with temperature rise are calculated and compared with experimental data in Figure 19. It is found that one obtains almost the same results for all cases during the first stage of heating when the crack density is small. However, important differences appear when approaching the peak value and in the post-peak zone. From a quantitative point of view, it seems that the peak value obtained in case 2 is the closest one to the experimental data. In this case, the positive part of Terzaghi effective stress is used as the tensile damage driving force. It is also seen that the contribution of thermal term A b θ involved in the cases 3 and 4 leads to a decrease of the peak axial stress probably due to the fact the temperature change rather enhances the compressive part of stress. From these results, it seems the tensile cracking is mainly driven by the positive part of Terzaghi effective stress. However, there are still some quantitative scatters between the numerical and experimental results.

As mentioned above, the numerical results represent the volumetric averages of non-uniform stress and pore pressure fields while the experimental data are taken on the sample boundary.

Moreover, due to the technical limitations, the post-peak responses were not measured in three tests. The maximum values given correspond to the macroscopic failure of samples. It is highlighted that at the scale of tested samples, the onset of localized cracks is strongly correlated to the heterogeneous distributions of pores and inclusions. By considering such material heterogeneity, it is not needed to introduce artificial weak element to enhance the onset of localized cracks. The obtained cracking patterns are quite consistent with experimental observations.

There is a clear interaction between the material anisotropy and heterogeneity in the cracking process. It is found that the spatial variation of inclusion distribution attenuates the influence of bedding planes on the orientation of localized cracks. In the case with non-uniform distribution of porosity only, the localized cracks are clearly guided by the orientation of bedding planes.

Under thermo-hydromechanical loading, the pore pressure rise due to temperature change plays an essential role in cracking process, but it is affected by the permeability variation due to open cracks. There is a strong coupling mechanism between the thermo-hydromechanical responses and cracking process.

The proposed model is applied to the analysis of laboratory heating tests in which macroscopic fractures were observed. Again, due to material heterogeneity, the local fields of stresses, strains and pore pressure are clearly not uniform inside the samples. This leads to the onset of localized cracks. The cracking patterns obtained in numerical simulations are overall consistent with that observed in the tested samples. With the temperature increase, the tensile cracking is the principal mechanism, which is mainly driven by tensile Terzaghi effective stresses. The present work can be improved by considering plastic deformation in clayey rocks. For the sake of simplicity, the detailed expressions of analytical solutions are not given here and can be found in [START_REF] Booker | Consolidation around a spherical heat source International Journal of Solids and Structures[END_REF]. In Figures 21 and22, we show the comparisons between numerical results and analytical solutions for temperature, pore pressure, displacement and stress. A very good agreement is obtained.

σ 0 ,

 0 p 0 , s 0 , T 0 and ρ 0 f denote the values of stress, fluid pressure, entropy, temperature and volumetric density of fluid at the initial reference configuration. The variation of temperature is defined as θ = T -T 0 . 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 drained thermo-elastic coupling coefficients. C b σ is the volumetric specific 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, for the sake of simplicity, the tensors of Biot coefficients, thermal dilation coefficients and drained thermo-elastic coupling coefficients are taken as isotropic ones. Therefore, one has B = bI, α b = α b I and A b = 3K b α b I, with K b being an equivalent bulk modulus and I the second order unit tensor. In order to deal with anisotropic elastic materials, the so-called Reuss equivalent bulk modulus is here used: K b = k Reuss (d t , d s ).The Biot modulus is accordingly given by 1/M = (b -ϕ)/K m + ϕ/K f , with K m being the bulk modulus of solid matrix, K f that of fluid and ϕ porosity. The differential thermal dilation coefficient is given by α m = (b -ϕ)α b + ϕα f with the thermal dilation coefficient of the fluid α f .

  3 and G ∥⊥ = 1.02 GPa. The elastic behavior of mineral inclusions is assumed to be isotropic and characterized by E = 98 GPa and ν = 0.15. The mean values of inclusion volume fraction and of porosity are respectively: β p = 0.16 and β i = 0.4. With these values, the homogenized elastic properties can be obtained: 28, ν hom ∥⊥ = 0.33 and G hom ∥⊥ = 2.7 GPa. These macroscopic elastic properties are close to those reported for COx claystone in

Figure 1 :

 1 Figure 1: Geometry of sample and boundary conditions (left), spatial distributions of porosity (middle) and inclusion volume fraction (right) for triaxial compression test

Figure 2

 2 Figure 2 presents the numerical and experimental results for the triaxial compression test with a confining pressure of 2 MPa. It is found that both the pre-and post peak responses are correctly reproduced by the phase-field model. However, from a quantitative point of view, the strains in the pre-peak part are underestimated by the model. This is due to the fact that only the elastic damage behavior is taken into account in the present model as the emphasis is put on hydro-thermal fracturing. It is also worth noticing that the numerical results in Figure 2 represent the averaged strains of the tested samples while the experimental data are the local strains measured by gauges.

Figure 2 :

 2 Figure 2: Axial and lateral strains versus differential stress in triaxial compression test with a confining pressure of 2 MPa: comparisons between numerical results (continuous line) and experimental data (dotted lines) reported in[START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: from experimental studies to model calibration and validation Computers and Geotechnics[END_REF]

Figure 3 :

 3 Figure 3: Distributions of tensile and shear damage in the triaxial compression test with 2 MPa confining pressure, at four levels of axial strain: (a) ε a = 0.95%; (b) ε a = 1.00%; (c) ε a = 1.05%; (d) ε a = 1.10%

Figure 4 :

 4 Figure 4: Geometry and boundary conditions of numerical heating test

Figure 5 :

 5 Figure 5: Distributions of temperature ( o C), pore pressure (MPa) and tensile damage at four heating steps of 600, 880, 1080 and 1200 minutes

Figure 6 :Figure 7 :

 67 Figure 6: Distributions of temperature ( o C), pore pressure (MPa) and tensile damage at 1200 minutes of heating by considering an isotropic porous medium

6 ×

 6 10 -6 T 4 + 1.66 × 10 -3 T 3 -0.197796T 2 + 16.862446T -64.319951) × 10 -6 (41)

Figure 8 :

 8 Figure 8: Geometry, initial and boundary conditions of thermal cracking tests

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Experimental variations of temperature (left) and of axial stress (right) during heating phase, used as input boundary conditions in numerical simulations

2 Figure 12 :

 212 Figure 12: Variations of axial and radial effective stress in three tests

Figure 13 :

 13 Figure 13: Distributions of tensile damage d t for test EXT3 at five heating times; comparison with experimental observation

Figure 14 :

 14 Figure 14: Distribution of temperature ( o C), pore pressure (MPa), radial, axial strains (%) and stresses

  are carried out for the test of EXT3, using different random distributions of porosity and mineral inclusion, as depicted in 15. The average porosity and volumetric fraction of inclusion remain unchanged from the previous simulation. For simplicity, the focus is limited to the variation of axial Terzaghi effective stress and distribution of tensile damage, see Figure16. Compared with the previous simulation which is considered as a reference, the curves are quite similar. The difference in peak value between the three calculations is less than 5%. This suggests that the numerical results are meaningful of the macroscopic responses of the sample. Although the location of cracks may vary due to the different distributions, the main cracks are primarily horizontal, indicating that the tensile damage is dominating in this test. It is worth noting that no vertical cracking patterns are observed in any of the calculations using different random distributions of material parameters.

(a) Random field n o 1 (b) Random field n o 2 (c) Random field n o 3 Figure 15 :n o 2 and n o 3 Figure 16 :

 12315316 Figure 15: Random distribution for porosity (left) and volumetric fraction of inclusion (right): cases n o 1, n o 2 and n o 3

Figure 17 :Figure 18 :

 1718 Figure 17: Comparisons of axial Terzaghi effective stress versus temperature between different random fields of material parameters in test EXT3, by using mesh of 9500 (left) and 38000 elements (right); and related distribution of tensile damage d t at the peak stress stat

Figure 19 :

 19 Figure 19: Comparisons of variations of axial Terzaghi effective stress between 4 cases in three tests EXT1, EXT2 and EXT3

Figure 20 : 4 . 5 × 1 . 7 -1 α s 4 . 2 × 10 - 5 Dynamic viscosity of water Pa • s µ w 1 × 10 - 3 Density of water kg • m -3 ρ w 1 × 10 3

 20451742105110313 Figure 20: Geometry and boundary condition of 3D heating benchmark

Figure 21 :Figure 22 :

 2122 Figure 21: Evolution of temperature, pore pressure, displacement (along x-axis) and normal stress (along x-axis) with time: comparisons between numerical and analytical results at two different locations

Table 1 :

 1 Mechanical parameters of COx claystone

Table 2 :

 2

	EXT1	-11.8	-12.6	4.9	25
	EXT2	-5.6	-4.8	2.1	24
	EXT3	-12.4	-7.2	4.0	35

Initial condition in three thermal cracking tests

[START_REF] Braun | Thermo-hydro-mechanical behavior of the Callovo-Oxfordian claystone: Effects of stress paths and temperature changes[END_REF]

Sample σ r 0 (MPa) σ a 0 (MPa) p 0 (MPa) T 0 ( o C)

Table 3 :

 3 Input parameters for thermal cracking tests[START_REF] Braun | Thermo-hydro-mechanical behavior of the Callovo-Oxfordian claystone: Effects of stress paths and temperature changes[END_REF]

	Density of solid grains	kg • m -3	ρ s	2770
	Biot's coefficient	-	b	0.8
	Permeability (parallel to bedding)	m 2	k h	4.0 × 10 -20

Acknowledgment:

This work is jointly supported by the French National Agency for radioactive waste management (ANDRA) and the DECOVALEX project. 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. The authors appreciate and thank the DECOVALEX-2023 Funding Organizations for their financial and technical support of the work described in this paper. The statements made in the paper are, however, solely those of the authors and do not necessarily reflect those of the Funding Organizations.

Appendix: Verification with analytical solution

In order to verify the accuracy of computing code for basic thermo-hydromechanical coupling problem, we consider here a three-dimensional heating benchmark for which an analytical solution is found in [START_REF] Booker | Consolidation around a spherical heat source International Journal of Solids and Structures[END_REF]. As shown in Figure 20, we consider a cubic domain with a length of 30 m. The cube is heated at its center with a power of 700 W. Due to the symmetry conditions, only a 1/8 of this cube is meshed. The material of cube is an isotropic, linear and homogeneous porous medium. The materials parameters are those of the COx claystone as shown in Table 4. Additionally, the water is assumed to be incompressible and the Biot coefficient is equal to 1.