Zhan Yu

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

Gilles Duveau

Meng Wang

Minh-Ngoc Vu

Carlos Plua

Numerical modeling of gas injection induced cracking with unsaturated hydromechanical process in the context of radioactive waste disposal

Keywords: Radioactive waste disposal, COx claystone, damage, hydromechanical coupling, biphasic flow, phase-field method

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 hydromechanical behavior of the host Callovo-Oxfordian COx claystone in quasi real conditions through the international research project DECOVALEX. In this study, we present a numerical analysis of damage and cracking behavior due to gas effects on hydromechanical processes in unsaturated claystone, while considering the distribution of material heterogeneity. The proposed model is implemented in a finite element code designed to solve hydromechanical coupling problems under unsaturated conditions. The nucleation and propagation of cracks are described using an extended phase-field method, which takes into account the effects of gas and liquid pressure on the evolution of the phase-field. In particular, a macroscopic elastic model is determined using two steps of homogenization, which considers the effects of porosity and mineral inclusions. The spatial variability of these factors is modeled using the Weibull distribution function. Thus, the nucleation of cracks is directly influenced by the spatial distribution of material heterogeneity. The proposed model is applied to 3D benchmarks of gas injection in the context of radioactive waste disposal. The process, in which the overpressure-induced damage zone affects the behavior of the two-phase flow during gas injection, is well reproduced.

Introduction

The emergence of nuclear energy in many countries has brought forth the concern of effectively managing and disposing of radioactive waste within deep geological repositories. In France, the Callovo-Oxfordian (COx) claystone formation has been identified as a prospective geological barrier for the containment of both intermediate-level long-lived waste (IL-LLW) and high-level waste (HLW) [1,[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 fundamental thermo-hydro-mechanical (THM) behaviors of the COx claystone have been assessed and characterized based on extensive laboratory tests conducted at a sample scale [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] Braun | Thermo-hydro-mechanical behavior of the callovo-oxfordian claystone: Effects of stress paths and temperature changes[END_REF]. Furthermore, a series of in situ experiments have been conducted and are currently underway under the coordination of the Andra within the Meuse/Haute-Marne Underground Research Laboratory (MHM URL) [1]. In order to characterize the thermo-hydromechanical responses of the host rock, the heat emitted by high-level waste (HLW) packages induces temperature elevation and pore pressure increase within the surrounding saturated porous rock are studied in the previous numerical works of DECOVALEX projet [START_REF] Seyedi | Upscaling thm modeling from small-scale to fullscale in-situ experiments in the callovo-oxfordian claystone[END_REF][START_REF] Plúa | A reliable numerical analysis for large-scale modelling of a high-level radioactive waste repository in the callovo-oxfordian claystone[END_REF], with additional information on individual modeling efforts [START_REF] Plúa | Effects of inherent spatial variability of rock properties on the thermo-hydro-mechanical responses of a high-level radioactive waste repository[END_REF][START_REF] Guo | Influence of the interval length of hydraulic packer systems on thermally-induced pore pressure measurements in rock[END_REF][START_REF] Wang | Analysis of coupled thermal-hydromechanical processes in callovo-oxfordian clay rock: From full-scale experiments to the repository scale[END_REF][START_REF] Thatcher | Assessing the hydraulic and mechanical impacts of heat generating radioactive waste at the whole repository scale[END_REF][START_REF] Xu | A study of thermal pressurization and potential for hydro-fracturing associated with nuclear waste disposal in argillaceous claystone[END_REF][START_REF] Yu | Numerical study of thermo-hydro-mechanical responses of in situ heating test with phase-field model[END_REF].

Alternatively, in the context of radioactive waste disposal, under anoxic conditions, the corrosion of metallic materials results in the generation of hydrogen. The radioactive decay of the waste and water radiolysis play a role in the overall production of gas. When the gas production rate exceeds the diffusion rate, gas accumulation continues until the pressure reaches a point where it enables the establishment of a biphasic flow between water and gas within the pores of the surrounding rock. [START_REF] Ortiz | Gas generation and migration in boom clay, a potential host rock formation for nuclear waste storage[END_REF][START_REF] Talandier | La production et le transfert de gaz dans le stockage et dans la couche du callovo-oxfordien lien avec le transitoire hydraulique site de meuse/haute marne, Rapport Andra[END_REF][START_REF] Amann | Gas breakthrough experiments on fine-grained sedimentary rocks[END_REF]. In order to better represent the processes governing advective transport of gas, specifically focusing on low-permeability argillaceous repository host rocks and clay-based engineered barriers, several numerical models are constructed by the work [START_REF] Tamayo-Mas | Modelling advective gas flow in compact bentonite: Lessons learnt from different numerical approaches[END_REF][START_REF] Chittenden | Development of an upscalable hm model for representing advective gas migration through saturated bentonite[END_REF][START_REF] Damians | Modelling gas flow in clay materials incorporating material heterogeneity and embedded fractures[END_REF][START_REF] Dagher | Investigating models to represent gas transport in a swelling geomaterial[END_REF][START_REF] Kim | Discrete dilatant pathway modeling of gas migration through compacted bentonite clay[END_REF][START_REF] Harrington | Bentonite homogenisation during the closure of void spaces[END_REF].

Indeed, the fracturing process plays an important role for the coupling behavior of rock. The numerical simulation of fracture is a scientifically challenging issue, and it has significant advancements in recent decades. In framework of continuum mechanics, the discrete method and distributed approach attracted the most attention [START_REF] Borst | Discrete vs smeared crack models for concrete fracture : bridging the gap[END_REF]. In the discrete method considering the crack as discontinuities, the widely used ones are the enriched finite element method (EFEM) [START_REF] Oliyer | Continuum modelling of strong discontinuities in solid mechanics using damage models[END_REF][START_REF] Yue | Strong discontinuity fe analysis for heterogeneous materials: The role of crack closure mechanism[END_REF] and the extended finite element method (XFEM) [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF][START_REF] Fries | The extended/generalized finite element method: An overview of the method and its applications[END_REF][START_REF] Zeng | An extended finite element solution for hydraulic fracturing with thermohydro-elastic-plastic coupling[END_REF]. They employ enrichment techniques to represent the displacement discontinuities for element part (EFEM) and node part (XFEM) respectively. On the other hand, the distributed approach introduces smeared finite band with concomitant high local strains to describe the discontinuity of crack, e.g., the high order strain gradient models [START_REF] De Borst | Gradient-dependent plasticity: formulation and algorithmic aspects[END_REF][START_REF] Chambon | Plastic continuum with microstructure, local second gradient theories for geomaterials: localization studies[END_REF], non-local damage models [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Jirásek | Consistent tangent stiffness for nonlocal damage models[END_REF] and the phase field method in most recently [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF].

The phase-field method is an approach based on the concept of the variational approach to brittle fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], and it was numerically implemented by [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF]. This model was compared and discussed with gradient-enhanced continuum damage model [START_REF] De Borst | Gradient damage vs phase-field approaches for fracture: Similarities and differences[END_REF] and peridynamics model [START_REF] Wu | Chapiter onephase field modeling of fracture[END_REF][START_REF] Roy | Chapter 15 -phase field-based peridynamics damage model: Applications to delamination of composite structures and inelastic response of ceramics[END_REF] because of their close relations. The primary advantage of the phase field method is its capability to capture the continuous transition from diffuse damage to localized cracks. This eliminates the necessity for introducing a specific criterion to detect the initiation of new cracks. Furthermore, the method can be readily extended to address multiple cracks in three-dimensional problems. Most recently, a phase field model considering irreversible thermodynamics was proposed to describe the ductile cracking in plastic material in [START_REF] Choo | Coupled phase-field and plasticity modeling of geological materials: From brittle fracture to ductile flow[END_REF]. The micro-structure based elasto-plastic phase field model is developed by [START_REF] You | Brittle-ductile failure transition in geomaterials modeled by a modified phase-field method with a varying damage-driving energy coefficient[END_REF][START_REF] Meng | A novel phase-field model for mixed cracks in elastic-plastic materials incorporating unilateral effect and friction sliding[END_REF]. The developed phase field models were proposed to consider the cohesive fracture in the works of [START_REF] Paggi | Revisiting the problem of a crack impinging on an interface:a modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model[END_REF][START_REF] Chen | Phase-field modelling of cohesive fracture[END_REF][START_REF] Feng | Phase-field cohesive fracture theory: A unified framework for dissipative systems based on variational inequality of virtual works[END_REF].

Given that many rock-like materials experience predominantly compressive stresses, the initiation and propagation of cracks under compressive loads in such materials exhibit intricate and mixed patterns [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]. In order to capture the rock failure behavior under cyclic loading, a so-called hybrid phase field model was firstly proposed in work of [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]. And this model was widely used in the following works of mixed crack considering tensile and shear crack [START_REF] Zhou | Phase field modeling of brittle compressive-shear fractures in rock-like materials: A new driving force and a hybrid formulation[END_REF][START_REF] Yajun | A new hybrid phase-field model for modeling mixedmode cracking process in anisotropic plastic rock-like materials[END_REF]. Therefore, to conveniently represent tensile and shear cracks, [START_REF] Fei | Double-phase-field formulation for mixed-mode fracture in rocks[END_REF][START_REF] Yu | Modeling of mixed cracks in rock-like brittle materials under compressive stresses by a double-phase-field method[END_REF] have introduced two independent damage variables. On the other hand, given that the phase field model has its flexibility to deal with multi-physical coupling problems, [START_REF] Zhou | Phase field modeling of hydraulic fracture propagation in transversely isotropic poroelastic media[END_REF][START_REF] Chen | Hydraulic fracturing analysis in fluid-saturated porous medium[END_REF] have been performed on the hydraulic fracturing for saturated pore media. [START_REF] Yu | Numerical study of thermo-hydro-mechanical responses of in situ heating test with phase-field model[END_REF] have introduced constitutive balance equations of THM behavior coupled with phase field model in the context of radioactive waste disposal. [START_REF] Wang | Phase-field modeling of cracking process in partially saturated porous media and application to rainfall-induced landslides[END_REF] developed the phase field model to analysis the progressive failure in saturated and unsaturated porous rocks in rainfall-induced landslides problem.

The primary objective of this study is to elucidate the progressive failure phenomenon resulting from gas injection conducted within an excavated underground borehole in the context of radioactive waste disposal. In this study, we present a numerical analysis of the damage and cracking behavior induced by gas effects on hydromechanical processes in unsaturated claystone, taking into account the distribution of material heterogeneity. In DECOVALEX-2023 Task A2, the so-called PGZ in-situ gas injection experiment was performed by Andra at the MHM URL. As a work of this project, we propose a hydromechanical model which takes into account the effects of gas and liquid pressure on the evolution of the phase-field. The proposed model is applied to analysis the PGZ in-situ test, and reproduce the rock failure due to coupling behavior in a three-dimension simulation. As a novel contribution, this study presents an enhanced phase field model to simulate the cracking process in both saturated and unsaturated porous media, specifically considering scenarios involving high gas injection rates. To better capture the behavior of tensile, shear, and mixed cracks during the excavation and gas injection phases, we employ two independent damage variables. We propose a criterion for each cracking process that takes into account the influence of water and gas pore pressure, as well as the degree of saturation. Furthermore, we develop a macroscopic elastic model through a two-step homogenization process that incorporates the effects of porosity and mineral inclusions. The spatial variability of these factors is modeled using the Weibull distribution function. The initiation of cracks is observed to be influenced by the spatial distribution of material heterogeneity.

Phase-field model for biphasic flow hydromechanical problems

Basic assumptions

In this study, we assume that a porous medium is saturated by a liquid phase (indexed as lq) and a gas phase (indexed as gs). The gas mixture phase consists of liquid vapor (indexed as vp) and nitrogen (indexed as N 2), while dry air is neglected, as nitrogen is the gas used in the studied gas injection experiment. Moreover, solute nitrogen (indexed as N 2 s ) is taken into account to diffuse inside the liquid according to Henry's Law. Thus, the basic assumptions for the relationship between the liquid and gas phases need to be clarified. Since the pore is saturated by both liquid and gas, the relationship between porosity can be expressed as:

ϕ = ϕ lq + ϕ gs (1)
The degree of saturation is defined as:

     S lq = ϕ lq ϕ S gs = ϕgs ϕ (2)
The degrees of saturation of liquid and gas are denoted by S lq and S gs respectively, and they are related as S lq + S gs = 1. As for the gas phase, it is related to pressure as follows:

   p gs = p vp + p N 2 p c = p gs -p lq (3)
Here, p i represents the pressure of the liquid and gas components, where i can take on the values lq, gs, vp, and N 2. Additionally, p c is the capillary pressure. By approximating Dalton's law, we make the further assumption that:

   ρ gs = ρ N 2 + ρ vp ∇ρ gs = ∇ρ N 2 + ∇ρ vp (4)
where ρ i is the volumetric mass with i=gs, N 2, vp and ∇ stands for its gradient.Given that the nitrogen, liquid vapor and gas mixture are assumed to be perfect gases, one gets:

1 ρ i = RT M ol i p i (5)
where M ol i , R and T are respectively the gas molar mass, perfect gas constant and absolute temperature. Therefore, the liquid is assumed compressible which has:

dp lq ρ lq = dp lq c lq ( 6
)
with the liquid compressibility c lq . The constitutive relations for the liquid and solute nitrogen are written as:

∇ρ N 2 s = H N 2 ∇ρ N 2 (7)
where H N 2 represents the Henry's constant of nitrogen. Since only water is considered as the liquid phase in this study, water and its vapor are assumed to be in thermodynamic equilibrium. For the sake of simplicity, the capillary effect at the contact of water and vapor is neglected. Therefore, the equilibrium between water and its vapor can be expressed using Kelvin's law:

ρ vp ρ 0 vp = exp( M ol vp ρ lq RT (p lq -p gs )) (8)

Variational framework 2.2.1. Phase field model

Following the spirit of Griffith's theory, Francfort and Marigo proposed a framework for a phase-field model of quasi-static brittle fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. Considering the aforementioned porous medium domain Ω with a crack set Γ, the variational principle is the minimization of energy, which is composed of stored energy and crack energy:

E(ε e , m lq , m gs , Γ) = ∫ Ω/Γ w e (ε e , m lq , m gs , Γ)dΩ + ∫ Γ G c dΓ ( 9
)
where w e represents the density of elastic energy, and G c is the energy required to create a unit of crack. ε e represents the elastic strain, m lq and m gs represent the mass of liquid and gas per unit of initial volume. In our previous work [START_REF] Yu | Modeling of mixed cracks in rock-like brittle materials under compressive stresses by a double-phase-field method[END_REF], we proposed a variational model considering mixed crack by following the regularized form of Bourdin [START_REF] Bourdin | The variational approach to fracture[END_REF]:

E(ε e , m lq , m gs , d t , d s ) = ∫ Ω w e (ε e , m lq , m gs , d t , d s )dΩ+g t c ∫ Ω γ t (d t , ∇d t )dΩ+g s c ∫ Ω γ s (d s , ∇d s )dΩ (10)
The discrete crack set is approximated and smeared by the so-called phase field of tensile crack d t and shear crack d s , which vary between 0 (undamaged) and 1 (fully damaged), respectively. g t c and g s c represent the material toughness of tensile and shear damage. The tensile and shear crack densities per unit volume are given by:

γ α (d α , ∇d α ) = 1 2 { 1 l d (d α ) 2 + l d ∇d α .∇d α } ; α = t, s (11)
with the length scale parameter l d , which can control the width of smeared cracks. In this way, the failure process induced by tensile and shear behaviors can be considered. This mixed crack model is significant, especially for the studied rock-like materials.

Elastic free energy

Mechanical properties of most unsaturated soils and rocks are influenced by water saturation degree. Comprehensive experimental data exploring this influence is thoroughly presented in the work of [START_REF] Delwyn | Soil Mechanics for Unsaturated Soils[END_REF]. Building upon this, the study adopts an approach similar to [START_REF] Nuth | Effective stress concept in unsaturated soils: Clarification and validation of a unified framework[END_REF], which extends Biot's theory from saturated to unsaturated porous media. This extension involves calculating an equivalent interstitial pressure for multiphase media based on the pore pressure and saturation degree of each phase. Such a calculation allows for the formulation of an effective stress tensor, essential for poroelastic modeling in both saturated and unsaturated media, while incorporating the capillary effect in the effective stress definition. This methodology aligns with the principles outlined in the work of [60] on the thermodynamics of open systems, providing a robust theoretical foundation for the approach. In this way, for an intact porous medium which is assumed to saturated by liquid and gas phase, the constitutive relations can be expressed as:

σ -σ 0 = C 0 b : ε e -B i (p i -p 0 i ) (12)
p i -p 0 i = M ij (-B j : ε v + m j ρ j ) (13)
with i and j denote the studied three events: lq, vp and N 2. The stress tensor is denoted by σ, and p i represents the partial fluid pressure. The superscript . 0 presents the initial stat of the variables. The drained stiffness matrix of the intact porous medium is denoted by C 0 b , while m j represents the mass and ρ j represents the density related to fluid phase j. To simplify the analysis, the Biot coefficient tensor related to fluid phase B j is assumed to be isotropic, such that B j = b j I. This results in b lq = S lq b and b vp = b N 2 = S gs b. The Biot modulus tensor is represented by M ij . In the framework presented in Equation ( 9), according to Equation ( 12) and ( 13), the density of elastic free energy, not taking into account any damage, can be expressed as:

w 0 e (ε e , p i ) = 1 2 σ b : ε e + 1 2 (p i -p 0 i ) 2 M ij with i, j = lq, vp, N 2, (14)
To accommodate the impact of fluid pressure on the damage driving force in this poroelastic model, we employ the Bishop effective stress tensor, represented as:

σ b = σ + b i p i .
Furthermore, following our recent work of saturated hydromechanical coupling [START_REF] Yu | Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling[END_REF], to differentiate the contribution of stored energy to tensile and shear failure behavior, it is assumed that the extension behavior (denoted by . + ) initiates tensile crack, while the compression behavior (denoted by . -) creates shear crack. Thus, we have:

w 0 e (ε e , p i , d t , d s ) = h t (d t )w 0+ e + h s (d s )w 0- e + 1 2 (p i -p 0 i ) 2 M ij with i, j = lq, vp, N 2, (15) and 
  w 0+ e = 1 2 σ + b : ε e w 0- e = 1 2 σ - b : ε e (16)
where the effective stress is decomposed by positive and negative part:

σ b = σ + b + σ - b . One gets: { σ + = P + σ : σ σ -= P - σ : σ ( 17
)
The decomposition is carried out using an operator Pσ ± that depends on stress. The details of the operator can be found in [START_REF] Lubarda | Damage model for brittle elastic solids with unequal tensile and compressive strengths[END_REF][START_REF] Murakami | Continuum damage mechanics: a continuum mechanics approach to the analysis of damage and fracture[END_REF][START_REF] Yu | Modeling of mixed cracks in rock-like brittle materials under compressive stresses by a double-phase-field method[END_REF]. Given that the rock studied in this work is transversely isotropic, we apply stress decomposition rather than classical strain decomposition proposed by [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], which relies on the elastic parameters of isotropic materials. Moreover, this approach to stress decomposition has been verified to accurately reproduce the crack patterns observed in laboratory-scale tests, as detailed in our recent work [START_REF] Yu | Modeling of mixed cracks in rock-like brittle materials under compressive stresses by a double-phase-field method[END_REF]. In this way, the failure behavior is described using a common degradation function

h α (d α ) = (1 -k)(1 -d α ) 2 + k,
where α = t, s. k is a small positive number introduced to prevent computational errors once the material has completely failed, indicated by d α =1.

Evolution of crack fields

Since the double phase field of mixed crack has been introduced multiple times in our previous work, we have chosen to present only a few important equations to ensure readability. To calculate the damage variable d α (α = t, s), the minimization problem in Equation 10 must be solved:

                   - ∂w e ∂d α -g α c δ d α γ α = 0 , ḋα > 0 , in Ω - ∂w e ∂d α -g α c δ d α γ α ≤ 0 , ḋα = 0 , in Ω ∂γ ∂∇d α .n = 0 , on δΩ (18)
Furthermore, the governing equations can be determined 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 0- s -g s c { 1 l d d s -l d div(∇d s ) } = 0 , ḋs > 0 (19)
where w 0+ e is the extension energy as shown in Equation 16. w 0- s is the compressive-shear energy, which is assumed as a part of the compression energy w 0- e , one gets:

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

{ ⟨a⟩ + = a, ⟨a⟩ -= 0 when a ≥ 0 ⟨a⟩ + = 0, ⟨a⟩ -= a when a < 0 (21)
The cohesion c and friction angle φ are introduced based on the Mohr-Coulomb theory.

It is important to note that the modified Mohr-Coulomb form presented here was initially proposed in the study referenced as [START_REF] Li | A damage model for hard rock under stress-induced failure mode[END_REF][START_REF] Zhou | Phase field modeling of brittle compressive-shear fractures in rock-like materials: A new driving force and a hybrid formulation[END_REF], and then further developed in our previous study [START_REF] Yu | Modeling of mixed cracks in rock-like brittle materials under compressive stresses by a double-phase-field method[END_REF]. To indicate compressive-shear behavior, we consider the negative part of the major and minor Bishop effective stresses, ⟨σ b 1 ⟩ -and ⟨σ b 3 ⟩ -. Finally, we incorporate the concept proposed by [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field fe implementations[END_REF] to account for irreversible failure processes, yielding:

           -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 (22)
with the history energy function defined as:

{ H t (t) = max[w 0+ e (τ )] τ ∈[0,t] H s -(t) = max[w 0- s (τ )] τ ∈[0,t] (23)

Biphasic flow model considering damage

In order to describe the hydromechanical behavior of porous medium considering the gas diffusion, several basic constitutive equations should be introduced. As this work takes into account solute gas, the mass balance between liquid (water) and vapor, nitrogen, and solute nitrogen can be expressed as:

   ṁlq + ṁvp = -div( ⃗ w lq + ⃗ w vp ) ṁN2 + ṁN2d = -div( ⃗ w N 2 + ⃗ w N 2d ) (24)
where ⃗ w is the flow rate. Darcy's law is considered for both of fluid and gas flow:

⃗ w i ρ i = - k i µ i k r i ∇p i , i = lq or gs ( 25
)
where µ is the dynamic viscosity, k is the second-order tensor of permeability, and k is the relative permeability. We apply Fick's law to model gas diffusion as well as the diffusion of solute nitrogen in liquid:

     ⃗ w i ρi = ⃗ wgs ρgs + ρgs ρ i D vp N 2 τ ϕ gs ∇ ρ i ρgs , i = vp or N 2 ⃗ w N 2 s ρ N 2 = H N 2 ⃗ w lq ρ lq - ρ lq ρ N 2 D lq N 2 τ ϕ lq ∇ ρ N 2 s ρ lq ( 26
)
The diffusion coefficients for nitrogen in vapor and liquid are denoted as D vp N 2 and D lq N 2, respectively, while τ represents the tortuosity. By applying Equations ( 13), [START_REF] Borst | Discrete vs smeared crack models for concrete fracture : bridging the gap[END_REF], and (26), we can rewrite the mass balance equation (Equation ( 24)) as:

   C lqlq ∂p lq ∂t + C lqgs ∂pgs ∂t + C lqε ∂ε ∂t = div(R lqlq ∇p lq ) + div(R lqgs ∇p gs ) C gslq ∂p lq ∂t + C gsgs ∂pgs ∂t + C gsε ∂ε ∂t = div(R gslq ∇p lq ) + div(R gsgs ∇p gs ) (27)
The complete expressions for C ij (i, j = lq, gs, ϵ) and R ij (i, j = lq, gs) are presented in appendix. In this study, the hydraulic parameters of biphasic flow are based on the van Genuchten's law [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils1[END_REF], one gets:

S e = [1 + ( p c p r
) n ] -m , with S e = S lq -S res S max -S res [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF] where S e is the effective liquid saturation, S res and S max are the residual and maximum degrees of saturation, respectively. The Van Genuchten parameters n and m are related by the expression m = 1 -1/n. The Mualem-van Genuchten model is used to determine the liquid relative permeability, and the resulting equation is as follows:

k r lq = S 1/2 e [1 -(1 -S 1/m e ) m ] 2 (29)
And the gas relative permeability is a cubic function: k r gs = S 3 gs . Furthermore, to incorporate the damage effect in the constitutive equations for the mechanical, liquid, and gas fields, we have followed the approach outlined in our previous work [START_REF] Yu | Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling[END_REF]. This involves the degradation for the elastic stiffness matrix as indicated in Equation ( 12), and an increase in the permeability tensors in Equations ( 27), both of which are linked to the two damage variables we proposed. Thus, we have:

C b (d t , d s ) = { h t (d t )P + σ + h s (d s )P - σ } : C 0 b ( 30
)
k i (d t , d s ) = k 0 i exp(η k max(d t , d s )), i = lq or gs ( 31
)
where the drained elastic stiffness matrix C 0 b , the liquid permeability k 0 lq , and gas permeability k 0 gs are used to represent the undamaged material. The evolution rate of permeability due to damage effect is denoted by η k , and it is assumed to have the same value for the evolution of both liquid and gas permeability in this work.

Numerical implementation in FEM

The biphasic poroelastic problem, which incorporates damage, can be solved numerically using the finite element method. By applying Equations ( 12), ( 27) and ( 22), the weak forms for five unknown fields can be obtained as:

∫ Ω δε : C b (d t , d s ) : εdΩ - ∫ Ω δε : (b lq δp lq )IdΩ - ∫ Ω δε : (b gs δp gs )IdΩ = ∫ S T t.
∫ Ω {(g s c /l d + 2H s -)d s δd s + g s c ∇d s • ∇(∂d s )}dV = ∫ Ω 2H s -∂d s dV ( 36
)
for shear damage field. After discretization using the shape function of FEM, we obtain the weak formulations of these five fields as follows:

           C ϵϵ ∆U
K d t d t = F d t K d s d s = F d s ( 37
)
As a result, the five fields are coupled with each other. To improve computational efficiency, we first solve for the liquid and gas fields using a fully coupled calculation. Then, we use stress-fixed iteration to solve the coupling problem between the displacement field and the pore pressure field (liquid and gas fields). Finally, we employ the Alternate Minimization solver to handle the coupling between the two phase fields and the displacement field.

Numerical simulations

Geometry and mesh

The simulation study focuses on two gas injection in-situ tests: PGZ1002 and PGZ1003. These tests were conducted in boreholes excavated from the GEX gallery at a depth of 450 m, as shown in Figure 1. Since the experiments were performed at the main level of the MHM URL, the initial pore pressure is assumed to be uniform: p 0 lq =4.7 MPa and p 0 mg =0.1 MPa. The initial stress state was determined by the major, minor, and intermediate principal total stresses, which were set to σ H = -16.1 MPa (parallel to PGZ1003), σ h = -12.4 MPa (parallel to PGZ1002), and σ v = -12.7 MPa, respectively. The temperature was kept constant at 20 • C throughout the experiment. In this in-situ test, the borehole is equipped with a multiple packer system to monitor the liquid and gas pressure in five isolated intervals (e.g. #01 and #02) as shown in figure 2. The borehole annulus between these intervals is filled with resin and packers. Since the gas injection is applied to the interval, the two packers isolating the injection interval are connected to a high-pressure vessel to stabilize the pressure fluctuations induced by the tests.

•mol -1 M ol N 2 0.028 Ideal gas constant J•mol -1 •K -1 R 8.314 Nitrogen-vapour molecular diffusion coefficient 10 -5 m 2 •s -1 D vp N 2
2.12 Nitrogen-water molecular diffusion coefficient It is worth noting that the elastic parameters presented in the table are homogenized values on a macroscopic scale. To consider the effect of material heterogeneity, the Weibull distribution is used to obtain the non-uniform volume fraction of pores and inclusions.

10 -9 m 2 • s -1 D w N 2 2.0 Nitrogen henry's coefficient 10 -2 H N 2 1.5
φ α = m α β α ( f α β α ) mα-1 exp[-( f α β α ) m α ] , α = p, i (38)
where f α represents the volume fraction of pores (α = p) and mineral inclusions (α = i).

The scale parameters β p = 0.16 and β i = 0.4 represent the mean values of the corresponding random variables, and m α = 1.5 is the homogeneity index of the material. The homogenization-based models allow the macroscopic elastic properties become functions of the non-uniform porosity and inclusion volume fractions at microscopic (pore) and mesoscopic (inclusion) scales. The resulting equation is:

C mp = (1 -f p )C m : [(1 -f p )I + f p (I -P p : C m ) -1 ] -1 ( 39
)
and

C hom = C mp + [f i (C in -C mp ) : D in ] : [I + f i (C in -I)] -1 (40)
A detailed description can be found in our previous work [START_REF] Zhao | A constitutive model for anisotropic clay-rich rocks considering micro-structural composition[END_REF], which also refers to the original work [START_REF] Giraud | Application of results on eshelby tensor to the determination of effective poroelastic properties of anisotropic rocks-like composites[END_REF]. Based on this, the reference values of elastic properties for the solid clay matrix are as follows: In order to save on computation costs, the aforementioned material heterogeneity is only considered in the elements near the borehole. To demonstrate the heterogeneity of the study domain, the distribution of Reuss equivalent bulk modulus has been chosen and is displayed as shown in Figure 5. In the framework of the phase-field method, the scale length l d controls the width of localized damage bands, which is typically related to the mesh size. In this study, l d is set to 0.02 m. Following previous work on COx claystone [START_REF] Yu | Numerical modeling of deformation and damage around underground excavation by phase-field method with hydromechanical coupling[END_REF], we define the parameters of the phase-field model as follows: tensile toughness g t c = 200 N/m, shear toughness g s c = 200 N/m, friction angle φ = 15 • , cohesion c = 0.1 MPa, and the rate of permeability due to damage η k = 20.

E s ∥ = 3.6 GPa, E s ⊥ = 2.4 GPa, ν ∥ = 0.3, ν ∥⊥ = 0.
In the 3D model studied, there are six external faces excluding the borehole annulus. Among these six faces, the back and top faces are considered as the far field, while the remaining faces are symmetric. For the boundary conditions of the liquid and gas field, p lq = 4.7 MPa and p mg = 0.1 MPa are set on the far field faces, and no flux is imposed on the symmetric faces. As for the boundary conditions of the displacement field, the initial stresses σ H , σ h , and σ v are applied to the far field in accordance with the directions of PGZ1002 and PGZ1003, respectively. Moreover, the normal displacements are fixed on the symmetric faces. These boundary conditions on the six external faces remain unchanged throughout the simulation.

Table 3 presents the different stages and their corresponding dates for PGZ1002 and PGZ1003 tests. The different stages can be characterized by the variations in the boundary conditions that occur on the borehole wall and the interior annulus face of the interval (referred to as the interval wall in subsequent sections), which are presented in table 4. After the instantaneous excavation in S1, the initial stress is released from the borehole wall, and the pressure of the liquid and gas is maintained at atmospheric pressure.

After the excavation (S1), the boundary conditions of the displacement field on the borehole wall are modified due to the mechanical pressure exerted by the packers, which is measured and taken into account in the simulation. The liquid and gas field becomes impermeable due to the installation of resin. The interior annulus of the interval is fixed after the installation of the interval. In S3, a fast gas injection is performed to obtain a saturation less than 1 in the interval, and in order to obtain a similar saturation value in both tests, the liquid pressure is adjusted accordingly. In S4, nitrogen gas flow is applied on the interval wall according to the time variation provided by the experiment. The details of the stages and their corresponding dates can be found in Table 3, and the changes in boundary conditions on the borehole wall and interval wall are summarized in Table 4. The time variation of the packer pressure and gas flow rate are displayed in Figure 6.

Stage

Simulation results of PGZ1003 test

The damage effect is a significant factor in the PGZ test. The experiment has observed that the damage is primarily induced by two operations: borehole excavation and gas injection. During the excavation, shear damage occurs, while tensile damage does not appear due to the discharge of initial stress, see Figure 7. In the case of PGZ1003, where the borehole was excavated parallel to the direction of major initial stress, the unload of initial stress is almost isotropic since the horizontal and vertical directions of initial stress have almost the same value (σ h = -12.4 MPa and σ v = -12.7 MPa). As a result, the damage zone mainly appears to the horizontal direction of the borehole excavated due to the anisotropy of COx claystone [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of callovo-oxfordian claystone: from experimental studies to model calibration and validation[END_REF]. During the water flow towards the excavated surface, the low liquid pressure affects the deeper area in the horizontal direction because COx claystone has a higher horizontal permeability. Based on the excavation-induced damage, the subsequent gas injection further aggravates the damage, as shown in Figure 8. As the shear damage zone seems to remain unchanged compared to the state after excavation, it is reasonable to conclude that gas injection mainly induces tensile damage. Therefore, the heavily damaged zone induced by gas injection appears in the vertical direction of the borehole. We attribute this to the fact that the vertical Young's modulus (E hom v ) is smaller than the horizontal Young's modulus (E hom h ), implying that the high pore pressure induced by gas injection can lead to greater displacement in the vertical direction of the borehole. Based on Equation [START_REF] De Borst | Gradient-dependent plasticity: formulation and algorithmic aspects[END_REF], damage results in an increase in permeability, which happens when the pores are connected by fractures. To facilitate readability, the distribution of the multiplier for increasing permeability induced by damage is shown in Figure 9. As the maximum value of shear damage (d s = 0.72) is greater than that of tensile damage (d t = 0.23), the shape of the zone with increased permeability is similar to that of the shear damage zone. Therefore, the area with the highest multiplier appears near the borehole wall in the horizontal direction. Figure 10 depicts the evolution of gas pressure at the studied points, namely PAy-PDy and PAz-PDz, over a period of gas injection and one day after injection phases. The location of these points is illustrated in Figure 11, and their coordinates are provided in the appendix. Except for PAy, which is defined on the internal surface of the interval, the other points are located at different depths within the rock. The experimental data is measured by the sensor PRE04, which is placed on the interval, as shown in Figure 3. In the PGZ1003 test, as the gas injection is composed by two phases, as shown in Figure 6, the gas pressure measurements exhibit two peaks corresponding to the moments of the two injection phases. It is apparent that the numerical results reproduce this feature, and the peak values are also captured by the evolution of PAy, which is on the interval. However, the numerical results show a greater reduction after injection compared with the experimental data. We attribute this discrepancy to the loss of gas flow towards the gallery.

As there is a lack of experimental data measured in the rock, it is not possible to compare the pressure evolutions at the other points. However, these pressure evolutions can still help us to understand the coupling process with the aid of the distributions shown in Figure 11. At the start of gas injection (t=0), there is almost 5 MPa gas pressure close to the borehole due to the water/gas exchange operation. The pressure reaches its peaks at 90 min and 173 min and then decreases. The high gas pressure can spread to the deeper area after injection, and the form of high pressure is slightly anisotropic. We attribute this anisotropy to the permeability anisotropy and the position of the excavation damage zone. For ease of reading, the simulation results of liquid pressure are displayed in the same way as before. Figure 12 shows the pressure evolution at different points, and Figure 13 shows the pressure distribution at different moments. Figure 14 illustrates the decrease in saturation caused by the gas injection phase. As shown in Figure 14-left, the interval is not fully saturated due to water/gas exchange even at the beginning of the gas injection. On the other hand, the rock component remains fully saturated for a certain period of time at the start of the injection, as depicted in Figure 14-right. Therefore, since the gas flow pushes the liquid pressure towards the deeper area, the liquid pressure also sharply increases due to gas injection. Similar to the distribution for gas pressure, the high-pressure zone for liquid pressure is also anisotropic, but the overpressure of the liquid can affect the deeper area. In this way, the entire process from excavation to gas injection has been simulated and reproduced using the proposed biphasic hydromechanical model that considers the effect of damage. The damage zones induced by excavation and gas injection can be described and distinguished using the proposed double-phase field that considers both tensile and shear damage. Since damage seems to play an important role in this coupling problem, further analysis is carried out in the next subsection.

Study of damage effect

To further study the role of damage in this gas injection test, sensitivity analyses are conducted to compare with the previous simulation results, which can serve as reference results.

In this sensitivity study, the authors compared the simulation results without considering the damage effect (denoted as II') with the reference results (denoted as I') that included the damage effect. Figure 15 shows the comparison for the evolution of gas and liquid pressure at different points. The simulation without considering damage used a poro-elastic model, and the peak of gas and liquid pressure reached almost 25 MPa at the point on the interval (PAy). For the points in the rock, the pore pressure (both gas and liquid) at the closest point (PBy) had a higher increase induced by injection. Figure 16 shows the distribution of liquid pressure at different moments without considering damage. However, due to the low permeability of the impact rock, the gas and liquid were relatively hard to flow towards deeper areas, resulting in less effect on pore pressure at deeper points. Another sensitivity study is performed on the parameter η k as shown in Equation [START_REF] De Borst | Gradient-dependent plasticity: formulation and algorithmic aspects[END_REF]. Since η k controls the rate of permeability increase induced by damage, calibrating this parameter can affect the evolution of gas and liquid pressure, as shown in Figure 19. Specifically, smaller values of η k result in a weaker effect of permeability induced by damage, leading to higher gas and liquid pressure. Conversely, greater values of η k result in smaller gas and liquid pressure. After analyzing the previous sensitivity studies, it is clear that excavation induced shear damage is the key factor in reproducing the reasonable behavior of gas and liquid pressure during the injection phase. The proposed model able to reproduce this coupling process due to the reasonable calibration of phase field parameters. However, since it is difficult to quantify g s c and η k based solely on experimental data, the choice of these parameters is more likely based on experience. It is important to acknowledge a significant limitation in this study: the intended fracturing process, induced by gas injection, does not manifest as expected. Instead of developing actual crack patterns, the area predicted to fracture is represented as a damage zone. This contrasts with our recent work [START_REF] Yu | Numerical analysis of hydro-thermal fracturing in saturated rocks by considering material anisotropy and micro-structural heterogeneity[END_REF], where we successfully replicated the fracturing process in a 2D heating benchmark problem, considering material heterogeneity in the same manner. The lack of fracturing in this study can be attributed to the chosen element size and the length scale parameter l d . The finer element size may help to capture the fracturing process in this study. Given that a finer mesh would substantially increase the number of elements in 3D calculations, our in-house code's computational limitations prevent us from handling such an increase. This issue could potentially be resolved by enhancing the computational efficiency of our code. Alternatively, employing 2D simulations with much finer meshes might also enable us to capture the fracturing process more effectively.

Simulation of PGZ1002 test

To investigate the effect of initial stress state on gas injection tests, a further study is conducted on the PGZ1002 test. Unlike the PGZ1003 test, the borehole excavation in PGZ1002 is oriented parallel to the minor principal initial stress σ h , resulting in an anisotropic stress state on the plane perpendicular to the borehole: σ H = -16.1 MPa and σ v = -12.7 MPa, as shown in Figure 1. This anisotropic stress state leads to a different excavation damage zone due to the unload of initial stress during excavation. As shown in Figure 20-right, the excavation induced shear damage zone is located vertically close to the borehole, consistent with experimental observations. Similar to PGZ1003, this shear damage zone does not develop during gas injection. The primary damage induced by gas injection is tensile damage, as shown in Figure 20 Using the same way of displaying, Figures 21 and22 show the evolution and distribution of gas pressure and liquid pressure. When compared with the experimental data, the evolution of gas and liquid pressure is well reproduced during the gas injection phase. The increasing trend and peak values are accurately captured. However, the numerical results show lower pressures than the experimental data during the post-injection phase, which is attributed to the same reason explained in PGZ1003. Moreover, comparing with the results of PGZ1003, the form of the overpressure zone changes for both gas and liquid. Unlike the high gas pressure zone, which affects a deeper area in the horizontal direction for PGZ1003, the high gas pressure zone in PGZ1002 appears to have a deeper effect in the vertical direction. Additionally, the form of the high liquid pressure has changed from an anisotropic one, with deeper effect in the horizontal direction for PGZ1003, to an almost isotropic one for PGZ1002. We attribute this change to the in the location of the excavation damage zone. The excavation damage zone, oriented in the vertical direction, results in the high gas pressure affecting a deeper area. On the other hand, while the vertical permeability is smaller than the horizontal one, this may offset the effect of anisotropic damage and leads to a isotropic high liquid pressure form. PGZ1003 and PGZ1002 tests provide a good way to study the effect of initial stress on the hydromechanical behavior of rock induced by excavation and gas injection. By comparing these tests, it becomes clear that the initial stress can directly affect the location of the excavation damage zone, which is well reproduced in the simulation study using the proposed model. Additionally, the two different locations of excavation damage zones lead to different rock behavior during the gas injection phase, particularly in terms of the distribution of gas and liquid pressure. This comparison highlights the crucial role played by damage behavior in this biphasic hydromechanical coupling problem, especially when considering high gas injection.

Conclusion

In this study, the double phase field method coupled with biphasic hydromechanical problems is being used to consider both the tensile and shear cracks in rock-like material, and to model the transition from diffuse micro-cracks to localized macroscopic cracks by taking into account material heterogeneities. The proposed model is then applied to analyze the unsaturated hydromechanical response in two in-situ 3D gas injection experiments: PGZ1003 and PGZ1002 tests. In general, the numerical results successfully replicate the main trends observed in the in-situ measurements, including the evolution of gas pressure, liquid pressure, and induced damaged zones. The simulation result emphasized the importance of the damage effect during these in-situ tests of gas injection. It indicates the shear damage induced by excavation plays more significant role in this application of biphasic flow coupling problem. However, the numerical study has limitations as it does not capture macroscopic cracks in the simulation, despite considering material heterogeneity. Utilizing a finer mesh and a smaller length scale parameter may be effective in reproducing the fracturing process. We plan to explore this in future work, either by applying 2D simulations or by enhancing the computational efficiency of our code. Additionally, understanding the mechanisms governing biphasic flow with high gas flow rates in porous media is crucial for verifying our assumptions. Therefore, conducting laboratory-scale tests is a necessary step to deepen our understanding and validate our theoretical models.

Appendix A: Complete terms of the diffusion equation for liquid and gas fields

The expressions of the terms in diffusion equation [START_REF] Yue | Strong discontinuity fe analysis for heterogeneous materials: The role of crack closure mechanism[END_REF] are presented as:

                                         C lqlq = S lq (
By using these shape functions, the matrix of the terms in Equation [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]

F d t = ∫ Ω 2H t N T p dV F d s = ∫ Ω 2H s -N T p dV (51)

Figure 1 :

 1 Figure 1: 3D view of the boreholes PGZ1002 and PGZ1003

Figure 2 :

 2 Figure 2: 3D view of the boreholes PGZ1002 and PGZ1003

Figure 3 :

 3 Figure 3: Geometry of PGZ test: 3D model and two cutting planes

Figure 4 :

 4 Figure 4: Finite element mesh of PGZ test

3 ,

 3 and G ∥⊥ = 1.02 GPa. The elastic behavior of mineral inclusions is assumed to be isotropic, with E = 98 GPa and ν = 0.15. The mean values of the inclusion volume fraction and porosity are β p = 0.16 and β i = 0.4, respectively.

Figure 5 :

 5 Figure 5: Distribution of Reuss equivalent bulk modulus

Figure 7 :

 7 Figure 7: Distribution of shear damage and liquid pressure induced by excavation

Figure 8 :

 8 Figure 8: Distribution of tensile (top) and shear (bottom) damage after gas injection

Figure 9 :

 9 Figure 9: Distribution of multiplier of increasing permeability after gas injection

Figure 10 :

 10 Figure 10: Evolution of gas pressure on the studied points, compare with the experimental data

Figure 11 :

 11 Figure 11: Distributions of gas pressure on A-A' cutting plan at different time of gas injection

Figure 12 :

 12 Figure 12: Evolution of liquid pressure on the studied points, compare with the experimental data

Figure 13 :Figure 14 :

 1314 Figure 13: Distributions of liquid pressure on A-A' cutting plan at different time of gas injection

Figure 15 :

 15 Figure 15: Evolution of gas (left) and liquid (right) pressure on the studied points, comparing with the simulation with and without considering damage effect

Figure 16 :

 16 Figure 16: Distributions of gas (top) and liquid (bottom) pressure on A-A' cutting plan at different time of gas injection, simulated without considering damage effect

FigureFigure 18 :

 18 Figure Evolution of gas (left) and liquid (right) pressure on the studied points, comparing with the simulations using different g s c : 180, 200 and 220 N/m

Figure 19 :

 19 Figure 19: Evolution of gas (left) and liquid (right) pressure on the studied points, comparing with the simulations using different η k : 18, 20 and 22

  -left. Comparing the simulation results of PGZ1002 with those of PGZ1003, the form and location of gas injection induced damage zone does not change significantly.

Figure 20 :

 20 Figure 20: Distributions of tensile (left) and shear (right) damage on A-A' cutting plan after gas injection for simulation of PGZ1002 test

Figure 21 :

 21 Figure 21: Evolution of gas (left) and liquid (right) pressure on the studied points for simulation of PGZ1002 test

Figure 22 :

 22 Figure 22: Distributions of gas (top) and liquid (bottom) pressure on A-A' cutting plan at different time of gas injection for simulation of PGZ1002 test

NρR lqlq = ρ lq k lq µ lq k r lq + ρ vp D vp N 2

 2 ) + (ρ vp -ρ lq )(ϕ ∂S lq ∂pc + Sgs N S lq ) + ρ vp ϕM ol vp ρ lq RT (1 -S lq ) C lqgs = (ρ lq -ρ vp )(ϕ ∂S lq ∂pc + Sgs N S lq ) -ρ vp ϕM ol vp ρ lq RT S gs + ρvpSgs N C lqε = b(ρ lq S lq + ρ vp S gs )C gslq =lq ϕ(1 -(1 -H H2 )S lq ) + ρ H2 S lq N (1 -(1 -H H2 )S lq ) + ρ H2 ϕ ∂S lq ∂pc (1 -H H2 ) C gsgs = M ol H2 RT (1 + ρvp ρ lq )ϕ(1 -(1 -H H2 )S lq ) -ρ H2 ϕ ∂S lq ∂pc (1 -H H2 ) + ρ H2 Sgs N (1 -(1 -H H2 )S lq ) C gsε = bρ H2 (1 -(1 -H H2 )S lq )Sgsτ ϕM ol vp RT ρ lq R lqgs = ρ vp kgs µgs k r gs -ρ vp D vp N 2 S gs τ ϕ(

Table 1 :

 1 Table1and 2 show the parameters of solid material and liquid/gas, respectively. The parameters for COx claystone 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[END_REF][START_REF] Zhao | A constitutive model for anisotropic clay-rich rocks considering micro-structural composition[END_REF]. COx claystone is considered as a transversely isotropic material, and the interval is considered as an isotropic material. Material parameters for the gas flow injection benchmark.

	Parameters	Unit	Notation COx Interval
	Young's modulus(parallel to bedding) Young's modulus(perpendicular to bedding)	10 9 Pa 10 9 Pa	E hom h E hom v	6.08 4.42	0.001
	Poisson(parallel to bedding) Poisson(perpendicular to bedding)	--	v hom hh v hom vh	0.28 0.33	0.3
	Shear modulus(perpendicular to bedding)	10 9 Pa	G hom vh	2.7	-
	Density of solid grains	kg•m -3	ρ s		2770	1631
	Biot's coefficient	-	b		0.8	1
	Porosity	-	ϕ		0.16	0.3
	Water permeability(parallel to bedding) Water permeability(perpendicular to bedding) 10 -20 m 2 10 -20 m 2 Gas permeability(parallel to bedding) 10 -20 m 2 Gas permeability(perpendicular to bedding) 10 -20 m 2	k lq h k lq v k gs h k gs v		4.0 1.33 400 133	10 8 10 8
	Van Genuchten gas entry pressure	10 6 Pa	p r		15.0	0.05
	Van Genuchten parameter	-	n		1.49	1.6
	Residual degree of saturation	-	S res		0.01	0.01
	Maximum degree of saturation	-	S max	0.99	0.99
	Parameters	Unit	Notation Value
	Water compressibility	10 -10 Pa -1		c lq	4.5
	Water density	10 3 kg•m -3		ρ lq	1.0
	Nitrogen density	10 3 kg•m -3		ρ N 2	1.25
	Water dynamic viscosity	10 -3 Pa•s		µ lq	1.0
	Gas vapor dynamic viscosity	10 -5 Pa•s		µ gs	1.0
	Water Molar mass	kg•mol -1		M ol vp	0.018
	Nitrogen Molar mass	kg

Table 2 :

 2 Water and nitrogen properties.

Table 3 :

 3 Timeline of the gas injection test: PGZ1002 and PGZ1003

	Stage	Borehole wall Mechanical Liquid	Gas	Interval wall Mechanical Liquid	Gas
	S1	f (σ 0 )	0.1 MPa 0.1 MPa	-	-	-
	S2	f (t)	Imperm. Imperm. Displ. fixed Imperm. Imperm.
	S3	f (t)	Imperm. Imperm. Displ. fixed <4.7 MPa 4.7 MPa
	S4	f (t)	Imperm. Imperm. Displ. fixed Imperm.	q N 2 (t)

Table 4 :

 4 Boundary condition of borehole wall and interval wall

	16							600	14							600
	14							500	12							500
	12							400								400
								300	10							300
	10							200								200
	8							100	8							100
	6	0	0.2	0.4	0.6	0.8	1	1.2 0	6	0	0.2	0.4	0.6	0.8	1	1.2 0
	Figure 6: Evolution of packer pressure and gas injection flow rate for PGZ1003 test (left) and PGZ1002 test
	(right)

  : C b (d t , d s ) : [B u ]dΩ C ϵlq = ∫ Ω [B u ] t : b lq : [N p ] t C ϵgs = ∫ Ω [B u ] t : b gs : [N p ] t : C lqlq : [N p ]dΩ C lqgs = ∫ Ω [N p ] t : C lqgs : [N p ]dΩ C lqϵ = ∫ Ω [N p ] t : C lqϵ : [B u ]dΩ R lqlq = ∫ Ω [B p ] t : R lqlq : [B p ]dΩ R lqgs = ∫ Ω [B p ] t : R lqgs : [B p ]dΩ [N p ] t : C gsgs : [N p ]dΩ C gslq = ∫ Ω [N p ] t : C gslq : [N p ]dΩ C gsϵ = ∫ Ω [N p ] t : C gsϵ : [B u ]dΩ R gsgs = ∫ Ω [B p ] t :R gsgs : [B p ]dΩ R gslq = ∫ Ω [B p ] t : R gslq : [B p ]dΩ

						can be expressed
	as:	  	C ϵϵ =	∫	Ω [B (48)
		          	C lqlq =	∫	Ω [N (49)
		     	C gsgs =	∫
		    			(50)
	   	K d t = K d s =	∫ ∫ Ω {(g t c /l d + 2H t )N T p N p + g t c l d B T p B p }dV
	  

u ] t p ] t Ω Ω {(g s c /l d + 2H s -)N T p N p + g s c l d B T p B p }dV

Acknowledgment:

This work is jointly supported by the French National Agency for radioactive waste management (ANDRA), the DECOVALEX project and the National Natural Science Foundation of China (No. 12202099). 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; DECOVALEX2023 is the current phase of the project. The authors appreciate and thank the DECOVALEX-2023 Funding Organisations ANDRA, BASE, BGE, BGR, CAS, CNSC, COVRA, US DOE, ENRESA, ENSI, JAEA, KAERI, NWMQ, NWS, SURAO, SSM and Taipower 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 Organisations.

Appendix B: The detail discrete approximation of FEM

The discrete approximations to the mechanical field and liquid, gas field are (i = lq, gs):

Appendix C: Coordinates of study points

The detail coordinates of the study points for the simulation of PGZ1003 and PGZ1002 are given as:

Point Coordinates (PGZ1003) Coordinates (PGZ1002) PAy (0, 0.033, 0) (0, 0.040, 0) PBy (0, 0.076, 0) (0, 0.055, 0) PCy (0, 0.114, 0) (0, 0.114, 0) PDy (0, 0.190, 0) (0, 0.190, 0) PAz (0, 0, 0.040) (0, 0, 0.033) PBz (0, 0, 0.076) (0, 0, 0.075) PCz (0, 0, 0.114) (0, 0, 0.114) PDz (0, 0, 0.190) (0, 0, 0.190)


domain: sdv sde shs

text:
Marie-Pierre Guinchard

Christophe Robin

Armand Guckert

Claude Bazard

Do carbohydrates contribute to winter survival of white clover?