Modelling an in-situ ventilation test in the Andra Underground Research Facilities

Wastes resulting from the nuclear electricity production have to be isolated from the biosphere for a very long period of time. For this purpose, deep underground repository in weak permeable geological layers is considered as a reliable solution for the nuclear waste storage. It is however well established that during excavation, the underground drilling process engenders cracks and eventually fractures [1] that deteriorate the hydro-mechanical properties of the surrounding host material in the so-called Excavation Damaged Zone (EDZ). The EDZ behaviour is a major issue because it may constitute a preferential flow path for radionuclide migration. Consequently, the characterisation of the material transport properties and of the transfer kinetics that occur around galleries still need to be investigated. The EDZ properties may be also affected by host rock-gallery air interactions. Ventilation induced drying may also provoke additional cracking, which potentially alters the transport properties of the damaged zone. Large-scale air ventilation experiments are performed in Underground Research Laboratories (URL) that have been constructed to check the feasibility of the repository. A numerical modelling of the SDZ air ventilation test (Andra URL) performed in a low permeability rock is proposed in order to both predict the development of the EDZ during excavation and study the air interaction with the host formation during maintenance phases.


Introduction
Wastes resulting from the nuclear electricity production have to be isolated from the biosphere for a very long period of time.For this purpose, deep underground repository in low permeable geological layers is considered as a reliable solution for the nuclear waste storage.Networks of galleries and access shafts will be constructed in these geological formations.The design of these underground structures must ensure that the good confining properties of the host rock will not be affected during the different steps of the storage.It is however well established that during excavation, the underground drilling process engenders cracks and eventually fractures [1] that deteriorates the hydro-mechanical properties of the surrounding host material.These modifications take place in a zone called Excavation Damaged Zone (EDZ), located around the galleries, which is affected by important modifications of the material flow characteristics such as permeability increase [2].Since a low hydraulic conductivity is required to ensure a safe long-term disposal, the EDZ behaviour is a major issue because it may constitute a preferential flow path for radionuclide migration.
Consequently, the characterisation of the material transport properties and of the transfer kinetics that occur around galleries still need to be investigated.
The EDZ properties may be also affected by host rock-gallery air interactions.Ventilation during excavation and maintenance phases has a direct influence on relative humidity and temperature of the gallery atmosphere.Ventilation induced drying may also provoke additional cracking, which potentially modifies the transport properties of the damaged zone.
In order to better understand the impact of these processes, large-scale air ventilation experiments are performed in Underground Research Laboratories (URL) that have been constructed to check the feasibility of the repository.The objective of the in-situ experiment (SDZ) performed in the Andra URL (Meuse Haute Marne region in France) is the characterisation of the structure and the size of the excavated damaged zone (EDZ), potentially influenced by the desaturation and possible resaturation due to controlled gallery ventilation conditions.
A numerical modelling of the SDZ air ventilation test performed in a low permeability rock is proposed in order to both predict the development of the EDZ during excavation and study the air interaction with the host formation during maintenance phases.The studied rock is the Callovo-Oxfordian claystone (COx), which is envisaged for deep underground repository of nuclear wastes in France [3].Numerous modelling of the EDZ behaviour have been realised in the past decades.Theses modelling often include separately hydro-mechanical coupling [4,5], permeability variation [6], flow transfers [7], or strain localisation [8].In the present study, the objective is to predict, in a fully coupled approach, the shear band pattern in the EDZ, the permeability increase within the fractures and the impact of ventilation on the two latter processes.The numerical developments were performed in the non-linear finite element code Lagamine [9].
EDZ evidences and the considered large-scale air ventilation experiment are firstly detailed in section 2.Then, the hydro-mechanical constitutive model is developed in section 3 for an unsaturated porous media.Finally, the gallery excavation and the air ventilation experiment are numerically modelled in the section 4. The first numerical phase is the gallery excavation modelling during which the development of the fractures that compose the EDZ is reproduced by shear banding [8].Considering the fractured rock as a continuous medium at the macroscale, the intrinsic hydraulic permeability evolution is reproduced through a straindependent relation [10].Succeeding to the excavation, the air ventilation experiment is reproduced in the underground gallery.In this second numerical phase, the exchanges at gallery wall between the claystone and the cavity air are characterised with a non-classical mixed hydraulic boundary condition [11,12] that involves both water vapour and liquid water transfers.

SDZ experiment
The in-situ experiment (SDZ), devoted to the characterisation of the structure and the size of the excavated damaged zone (EDZ), is located in an experimental zone isolated from the GED experimental horizontal gallery by a 2 m long airlock.The test zone contains a first part without covering (with a length of 7.2 m) and, at the end, another 5 m long zone with a concrete covering with a thickness of 20 cm and a thin impervious geotextile (Figure 1).

Figure 1. SDZ experiment zone in the GED gallery
A series of drilling in the SDZ zone and in the GED gallery allow for a description of the fracturing around the galleries as well as permeability and water content measurement.During the test, the pore pressure evolution in different directions and radial distances is monitored.Furthermore, a geological survey of the fractured zone as well as measurements of the gallery convergence is also performed by Andra after the gallery excavation.All these measurements allow for a better observation of the coupled hydro mechanical couplings in the EDZ.

Excavation damaged zone
The EDZ has been carefully investigated in the Andra URL through, for instance, fracture measurements, permeability analysis, and flow transfers.For the Callovo-Oxfordian claystone, induced extension and shear fractures are detected in the proximity of the galleries, with shearing as the principal failure mechanism because of the high in situ stress environment [13].The SDZ ventilation experiment is located in the GED experimental gallery, which is oriented parallel to the minor horizontal principal stress.Around this gallery, the permeability increase has been highlighted by measurements performed under saturated conditions in boreholes that are drilled in different orientations.Three zones can be defined: an undisturbed zone with a permeability lower than 10 −19 m 2 far from the gallery, a slightly disturbed zone with a permeability ranging from 10 −19 m 2 to 10 −17 m 2 , and a highly disturbed zone close to the gallery with a permeability higher than 10 −17 m 2 (increase higher than 2 orders of magnitude).The extents of the zones are detailed in Fig. 2.

Ventilation experiment
Ventilation can lead to rock desaturation, stress modification, as well as modification of the fracturing structure close to the drifts.The major objective of the SDZ ventilation test is to characterise the influence of controlled ventilation on the hydraulic transfers taking place around the galleries, especially in the EDZ.In the experimental zone isolated from the rest of the laboratory, different levels of hygrometry are applied in order to highlight the effect of drainage and wetting.The ventilation can be sequenced in different phases (Fig. 3): (1) the first phase corresponds to the global laboratory ventilation before the airlock closure, (2) once the airlock is closed the ventilation is stopped in the experimental zone, later a controlled ventilation is imposed with hygrometric conditions of (3) 30 % of relative humidity and then of (4) 60 % of relative humidity.The air temperature T and relative humidity RH evolutions are monitored in the experimental zone in different gallery sections.
The pore water pressure evolution is monitored inside the rock mass, around the experimental zone without covering.After equilibrium of the sensors with the rock mass, the progressive drainage of the rock is monitored with a decrease of the pore water pressure.One can observe that the drainage is not important in the far field, where the pore water pressures are close to the initial rock water pressure.It is more important around the gallery where the measurements are about the atmospheric pressure (in a 3 m ring).In the excavation damaged zone, we can expect that the increase of permeability and the drainage lead to rock desaturation.In the long term, the two controlled ventilation phases (phases ( 3) and ( 4) in Fig. 3) only have a minor effect on the measured pore water pressure.To characterise the desaturation, water content measurements have also been performed on core samples coming from horizontal boreholes drilled in the experimental zone without covering.They indicate, firstly, a strong desaturation close to the gallery wall with 3 % ≤ w ≤ 5 % that can extend up to 0.5 m depth, secondly, a moderated desaturation deeper in the rock, and thirdly, a limited desaturation or no desaturation deeper than 2 m with 6.7 % ≤ w ≤ 8.2 % for the saturated claystone.

Framework for an unsaturated media with microstructure
The observations and the measurements of the SDZ experiment clearly evidence the existence of the excavation damaged zone, where modifications of mechanical and flow properties appear.In this work, the numerical modelling of the EDZ aims at representing the complexity of this zone in a fully coupled approach.The goal is to be able to predict the localisation of the deformation in shear bands, which are the precursor of fractures [14].Inside these localisation bands, the permeability is deeply increased.In fact, the coupling between mechanical deformation, permeability, and hydraulic flow is indubitably of paramount importance to obtain an adequate representation of the EDZ hydromechanical behaviour.
The first step is to define a correct and robust modelling of the strain localisation leading finally to the rupture.However, it is well known that classical finite element models suffer of a pathological problem of mesh dependency for strain localisation modelling.This is due to properties of the underlying mathematical problem.There is a need of an internal length scale for a correct modelling of the post peak behaviour.Several approaches have been proposed to introduce such internal length scale in the problem.Among the different regularization techniques, we propose in the following to use the second gradient model [15], which is a particular model of a microstructure medium.The second grade model has been further extended to saturated biphasic medium [16] and unsaturated conditions [8].
A more pronounced permeability increase inside the fractures (shear bands) due to strain localisation effect is expected.Different approaches exist in the literature to predict the evolution of the permeability with the mechanical loading (micromechanics, FE 2 computation, semi-empirical formulation …).A strain-dependent evolution of the material intrinsic permeability will be used in this work, in order to introduce this hydromechanical coupling effect.
The second step is the modelling of the gallery air ventilation.The SDZ experiment in the experimental zone without covering will be reproduced after the gallery excavation.The air-rock interaction has also to be characterised to reproduce the exchanges at gallery wall correctly.
In the following numerical modelling, these two steps will be investigated with an elaborated description of the hydro-mechanical coupling, especially within the EDZ.
To describe the numerical model, the balance equations are firstly recalled.Then the constitutive equations of the mechanical and the fluid transfer problems are presented, with special emphasis on the different coupling existing between the mechanical and the hydraulic parts.

Balance of momentum
Let us recall first the governing equations of a monophasic medium with microstructure.This kind of models can be traced back to the pioneering work of the Cosserat brothers [17], Mindlin [18] and Germain [19].In the framework of microstructure continuum theory, a microkinematic gradient field v ij is introduced to describe strain and rotation at the microscale.With respect to classical continuum mechanics, additional terms are then added in the internal virtual work of a given body [19].
In the local second gradient model used in the following, an assumption is added: the micro kinematic gradient is equal to the macro-deformation gradient     = Finally, the principle of virtual work can be rewritten as follows:  is the double stress dual to the second gradient of the macro displacement.
Assuming that the boundary Ω t is regular, the external virtual work can be defined as follows: where  ̅  is the external (classical) forces per unit area and  ̅  an additional external (double) force per unit area, both applied on a part Γ  of the boundary of Ω t .D denotes the normal derivative of any quantity q, (Dq = (q/x k )n k ).
The previous virtual work equation of second gradient models can of course be used in a finite element code.However, this equation needs the use of C 1 functions for the displacement field as second derivatives of the displacement are involved.In order to avoid such functions in the virtual work principle (Eq.2), the equalities between   and   and between   * and   * are introduced through a field of Lagrange multipliers  ij related to a weak form of the constraint (Eq. 1) (see [15]):

Fluid mass balance equation
The water mass balance equation reads in a weak form: where   * is the virtual pore water pressure field, Q t is a sink term and Γ  is the part of the boundary where the input water mass per unit area  ̅ is prescribed.M t and    are respectively the mass of the water inside the current configuration of the skeleton Ω t and the mass flow.
The mass flow    is defined as follows: where κ ij is the anisotropic intrinsic permeability,  ,  is water relative permeability and μ w is the water dynamic viscosity.
The compressible fluid is assumed to respect the following relationship [20].This predicts an increase of water density as a function of the pore water pressure, defining χ w as the water bulk modulus: The time derivative of the water mass is obtained directly by considering the influence of compressibility of the fluid phase and the solid grains, variation of water saturation and modification of the porous volume.

Constitutive law
Geomaterials exhibit a porous structure where an assembly of solid grains forms the solid skeleton and fluids can fill in the porous space between these grains.Generally, two fluids are considered: a wetting and a nonwetting fluid corresponding respectively to a liquid and a gaseous phase.The effect of fluid pressures on the Cauchy total stress field  ij is defined assuming Biot's postulate of effective stress [21].This definition accounts for the relative deformability of the solid grains by introducing its compressibility through the Biot's tensor: where   ′ is the Biot's effective stress field defined under soil mechanics convention in which compressive stress is positive, b ij is Biot's tensor and p w is the pore water pressure, negative if suction.An isotropic transverse linear elastic-perfectly plastic model (with a Van Eekelen yield surface) is assumed to reproduce the mechanical behaviour of the porous media [22].The elastic and plastic theories are defined hereafter for anisotropic materials, which means that the elastoplastic properties depend on the coordinate axes to which the properties are referred and therefore depend on the orientation.
By separating the elastic and plastic components of the strain tensor, a linear elastic relationship links the elastic strain rate    to the effective stress rate   ′ through the Hooke's law.   is the Hooke anisotropic elastic tensor [23] that must be symmetric because of thermodynamic requirement [24].In case of orthotropy and cross-anisotropy, the properties are related to the orthotropic (and orthogonal) axes e 1 , e 2 and e 3 and the (e 1 ; e 2 ; e 3 ) space represents the orthotropic configuration [25].For cross-anisotropy, five independent parameters are needed to describe the material elasticity: E  ; E // ;  //// ;  // ; G // , where the subscripts // and  indicate the directions parallel and perpendicular to the bedding, respectively.We consider that (e 1 ; e 2 ) is the isotropic planes orientation (bedding planes for sedimentary materials) and e 3 is the normal to these planes.The Hooke's law being formulate in the orthotropic axes, a change of the reference system has to be computed to obtain the stress state.
Another important aspect of the anisotropic elasticity is the compressibility of the solid grains skeleton.It is expressed through Biot's symmetric tensor which is defined for anisotropic material as follows [26]: [10] where   is the Kronecker symbol and K s is the bulk modulus of the solid phase.In addition to the anisotropic elastic behaviour, the plastic behaviour and the strength parameters of the material can also be anisotropic and therefore depend on direction and orientation.Various failure criteria for anisotropic materials have been developed [27] and hereafter, we focus on the material cohesion anisotropy which is defined with a second order microstructure fabric tensor a ij describing the cohesion spatial distribution.This tensor is a measure of the material fabric and its eigenvectors correspond to the principal material microstructure axes, i.e. the orthotropic axes, e 1 ; e 2 ; e 3 .The cohesion corresponds to the projection of this tensor on a generalised unit loading vector l i that characterises the loading direction relative to the material microstructure orientations [28].Each component l i corresponds to the stress resultant acting on facets of normal e i and is defined as follows: #′2   #′2 [11] where   #′ is expressed in reference to the material axes.
Therefore, the cohesion specifies the effect of load orientation relative to the material axes [29]:

Results and comparison
When considering deep underground repository of nuclear waste, there is a need for a more accurate modelling of the EDZ and of the gallery air ventilation effect on the flow transfers taking place around the galleries.Therefore, a hydro-mechanical modelling of gallery excavation and of the SDZ ventilation test is realised in two-dimensional plane strain state.

Numerical model
The considered rock is the Callovo-Oxfordian claystone which is a transversely isotropic material with horizontal bedding planes and an initial anisotropic stress state.For a gallery oriented along the minor horizontal principal stress, the initial stress state and pore water pressure are: where p w,0 is the initial pore water pressure and σ x,0 , σ y,0 , σ z,0 are the initial total stresses.
The modelling considers the SDZ experimental zone without covering and only one gallery quarter for symmetry reason.The mesh geometry is a gallery with a 2.3 m radius and is discretized with 2520 elements and 9802 nodes.To establish the symmetry, the symmetry xand y-axes are considered as impervious (no water flow) and the normal displacement, as well as the normal derivative of the radial displacement, cancel along these axes [8].At the mesh external boundaries, the normal total stress and the pore water pressure are imposed constant.

Excavation phase
As a first step, the drilling process and the EDZ development are modelled.During the excavation, the major objectives are to characterise the development of fractures induced by rock deconfinement and the modification of the hydraulic properties.The fractures are represented with shear banding and the permeability evolution is addressed with a strain-dependent relation.
The gallery excavation is performed with a convergence-confinement method [30] and decreasing stresses at gallery wall with the excavation front crossing the studied section after one week.The pore water pressure in the gallery is assumed to decrease rapidly to one atmospheric pressure when the excavation front crosses the studied section.It is assumed that the gallery air remains fully saturated with water vapour because the gallery is not ventilated yet.It means that the air relative humidity is of 100 % and the water pressure in the cavity corresponds to the atmospheric pressure.The evolution of the intrinsic hydraulic permeability is reproduced with a strain-dependent formulation.The evolution parameters for the permeability are calibrated based on the shear banding results and on in situ experimental measurements of permeability around the GED gallery (Fig. 2).An appropriate solution considers a formulation mixing the deformation and the plastic indicator (yield index).The following equation uses the equivalent deviatoric strain ̂  and the yield index YI, leading to: The results for β = 1 × 10 10 and YI thr = 0.95 illustrated in Fig. 4 show a permeability increase in a quite good agreement with the experimental measurements for all the directions.The permeability evolution during the excavation and the comparison to experimental data are DOI: 10.1051/ 04003 (2016) , 9 E 2016 consistent.In the far field, the permeability corresponds to the undamaged rock and is well captured in the modelling.In the near field, the permeability corresponds to the fractured rock (EDZ).

Ventilation phase
The second part of the numerical modelling consists in applying, after the excavation, an air ventilation in the underground gallery that corresponds to the SDZ ventilation experiment.After the first step, the fracture description as well as the permeability evolution and inhomogeneity in the EDZ have been taken into account.The following challenge is to analyse whether the introduced hydro-mechanical coupling and EDZ description allow reproducing the drainage and desaturation provoked in the surrounding media by the air ventilation.The principal objectives of the SDZ test modelling are to characterise, firstly, the air-rock transfers at gallery wall, and secondly, the influence of a controlled gallery ventilation on the claystone behaviour.After the drilling, the ventilation is progressively applied in the underground gallery with a drying of the air which is numerically reproduced by a ventilation initiation phase.During this initiation phase (between 21 and 50 days of simulation), the water pressure in the gallery air decreases from the atmospheric pressure to a value corresponding to the first measurement of air humidity.Then, the air ventilation is applied in the gallery by means of the hygrometric measurements.During the ventilation, the pore water pressure and the water vapour density in the cavity correspond to the relative humidity of the cavity air by Kelvin's law.
The first result that can be analysed is the claystone progressive drainage which occurs due to the ventilation.In Fig. 5 is illustrated the comparison between numerical results and experimental measurements of pore water pressure evolution in horizontal boreholes drilled around the drift.It is obvious that the equilibrium phase of the measure (increase of p w ) is not reproduced by the numerical modelling which starts with the initial pore water pressure value of the clayey rock, p w,0 = 4.5 MPa.Only the drainage is reproduced and has to be compared to experimental values after about 275 days of modelling.A close correspondence is observed between the numerical results and the experimental measurements in the different directions.The pore overpressures observed after the excavation in the vertical direction result of hydro-mechanical couplings mainly related to the anisotropy of the initial stress state.In this direction, the matching with the experimental data is slightly less good which may be related to the lower vertical permeability (k⊥⊥ < k  ) and to the positions of the shear bands that condition the intrinsic permeability evolution.
The vapour mass transfer coefficient used to obtain the results of Fig. 5 is of α = 10 −3 m/s.This transfer coefficient is introduced through a non-classical mixed hydraulic boundary condition [11,12] to represent the long-term exchange (evaporation) at gallery wall.It has indeed an influence on the drainage and its value is chosen to obtain a good reproduction of the experimental data.In the long term, the two controlled ventilation phases (phases (3) and (4) in Fig. 3) do not have a visible effect on the numerical results.

Conclusions
The presented study is devoted to the numerical modelling of failure process and hydraulic transfer occurring around underground galleries subjected to air ventilation, in a cross-anisotropic argillaceous rock.The fractures induced by the excavation process in the vicinity of the gallery are reproduced by shear banding, and the strain localisation is properly modelled by means of a regularisation method.The evolutions of flow transfers and transport properties have been characterised around the underground drift, especially in the excavation damaged zone.The evolution of intrinsic permeability is addressed by using a strain-dependent relation that involves plastic strain as well as a restricted part of the elastic strain.At the end of excavation, the permeability increase is of several orders of magnitude in the shear bands as measured experimentally.
The flow transfers in the damaged zone and in the undisturbed material are studied through the modelling of a large-scale gallery air ventilation experiment (SDZ).The air-rock interaction at gallery wall involves both water vapour and liquid water transfers that are implemented with a mixed hydraulic boundary condition.Such condition is relevant for the reproduction of the transfers because an instantaneous equilibrium between the air and the rock may not be assumed beforehand.Equilibrium is reached in the long term by evaporation process during the gallery air ventilation.Regarding the different ventilation phases, it has been observed that the progressive rock drainage is mostly affected by the longterm ventilation, not by the constant ventilation phases of the SDZ experiment.
The numerical results which are presented indicate that a strain localisation approach for the EDZ modelling, linked to an intrinsic permeability evolution, allows to

Figure 2 .
Figure 2. SDZ experiment zone in the GED gallery

Figure 3 .
Figure 3. Temperature and relative humidity evolution in the SDZ gallery.

Figure 4 .
Figure 4. Intrinsic permeability at the end of excavation in vertical (black), oblique (grey) and horizontal (dotted line) direction.

Figure 5 .
Figure 5. Pore water pressure evolution in horizontal boreholes with comparison to experimental measurements for different distances from gallery wall transfer in the Callovo-Oxfordian claystone.