Anisotropic hydromechanical modelling of a deep drift at the Meuse/Haute-Marne URL

The Callovo-Oxfordian claystone is considered as a potential geological host formation for high-level and intermediate-level long-lived radioactive waste in France. Since 2000, the French National Radioactive Waste Management Agency (Andra) has been constructing an Underground Research Laboratory (URL) at Bure with intent to demonstrate the feasibility of a geological repository in the Callovo-Oxfordian claystone (COx) formation. The excavation of galleries at the main level (490 m depth) of the URL showed a significant fracturing induced by the excavation (Armand et al. 2014). Indeed, the authors have reported the existence of two adjacent zones around the laboratory galleries: a "connected" fracturing zone and a more extended "discrete" fracturing zone overhanging the first zone. The "connected" fracturing zone would correspond to the well-known EDZ (Excavation Damaged Zone) and would represent a domain for which the fracturing is particularly developed, whereas the "discrete" fracturing area located between the EDZ and the intact rock mass would represent the EdZ (Excavation Disturbed Zone). The in situ observations show an anisotropic extent of the induced fractured zone even for the drifts excavated following the in situ major horizontal stress for which, the initial total stress is quasi-isotropic in the drift section. Different factors contribute probably to this anisotropic response of the COx to the excavation operation. Based on the experimental results, various failure criteria for anisotropic materials exhibiting a visible inherent anisotropy have been proposed and constitutive models were developed in the literature. The developed models can be split mainly in (a) the empirical models based on the theory of variational cohesion and / or friction (e.g., Wang et Yu 2014); (b) the models built on the concept of ubiquitous joints with several planes of weakness (Sainsbury et al. 2008), (c) the models where damage and/or plasticity are incorporated and formulated in the framework of irreversible thermodynamics (Pietruszczak et al. 2002). For this purpose a phenomenological macroscopic anisotropic model (including both the elastic anisotropy and the induced anisotropic plasticity) has been proposed in this paper. The basic assumption is that the failure of an anisotropic material is due to either fracturing of bedding planes (or weakness planes) or the failure of the rock matrix. The formalism is therefore based on approaches with weakness planes commonly called “Discontinuous weakness plane” or “Ubquituous joints” aimed to account for rock mass strength and its anisotropy within large-scale continuum approaches. In the proposed model, any predefined weakness plane is not considered. The joints are introduced following an estimated orientation of the induced fractures to reproduce the induced anisotropy. In this framework, it is assumed that the rock is composed of a matrix and of potential planes of weakness (reflecting the induced damage approached by the theory of plasticity) as observed on the results of biaxial tests under plane strain conditions with X-ray micro-tomography (Besuelle and Lanata 2014) aiming to characterize the deformation mechanism. Then the matrix is assumed to be linear, transversely isotropic and the plasticity is described by an isotropic non linear yield function derived from the laboratory characterization. A non-associated flow rule is used with a distinction between compression and extensional stress paths, as well as the absence of volumetric strain beyond large plastic distortion. Induced fractures are modelled as the planes of weakness represented by joints. Experiments on brittle failure reveal two fundamental types of fractures: tensile/extension (mode I) and shear (modes II) which individually induces different orientation of the facture failure plane relative to the principal stresses. Mode I, fractures are associated to (a) tensile fractures occurring when the minimum principal stress s3 reaches the tensile stress, or (b) longitudinal splitting when the minimum principal stress is close to zero as in uniaxial compression. The fracture orientation is perpendicular to the tensile stress (tension case) or parallel to the major compression stress (uniaxial compression). Mode II, shear fractures occur under triaxial compression configuration with fracture plane angles θ less than 45° with the maximum compressive stress. In the latter case, we assumed that the rocks fail along conjugate fractures with θ=±(45° - ϕw/2) where ϕw is the friction angle along the weakness plane. This is in accordance with the experimental results of failure plane orientation reported by Zhang (2016). Finally, a perfectly plastic behaviour according to the Mohr-Coulomb criterion is assumed along these weakness planes; while the elastic part is considered as linear and transversely 396 isotropic. The constitutive equations were implemented in the three-dimensional explicit finite-difference code, FLAC3D, where fully coupled hydromechanical modelling can be performed under fully saturated conditions for both isotropic and anisotropic rock masses. The ability of the proposed model to reproduce the plastic zones around a drift excavated following in situ major horizontal stress (GCS drift) at the main level of the Meuse/Haute-Marne (M/H-M) URL, is firstly successfully tested from a purely mechanical (monophasic) approach. Comparison between simulation results and the in situ measurements around the GCS drift provides new insights on the understanding of the deformation mechanisms observed around the structures of the M/H-M URL. More precisely, numerical simulation predicted convergence ratios (horizontal / vertical) at different measurement section positions along the length of the gallery, which remain between 1.5 and 2 as illustrated in Fig. 1. Those obtained with in situ measurements for the instantaneous response are ranged between 1.3 and 2.1 (Guayacan-Carillo et al. 2016). The extension of the predicted plastic zones with the proposed model varies from 0.7xD to 1xD (Fig.2), which well corresponds to the in situ observations. Indeed, according to Armand et al. (2014), the extent of EDZ is 0.15xD and 0.5xD in the roof/floor and drift sides, respectively, against 0.5xD and 1xD for EdZ. Hydromechanical coupling was also investigated with the proposed model on the same GCS drift. For simplicity, the coupled hydromechanical simulation was carried out on a 2D geometry (a cross section of the GCS drift, excavated with a drainage condition along its wall). The initial pore pressure at the main level (-490 m) is about 4.7 MPa. The instantaneous (0+) response leads to interstitial over-pressures with a maximal amplitude of 1.5 MPa in the direction of the initial minor stress (horizontal) and underpressures in the direction of the initial major stress (vertical) with respect to the GCS section (Fig. 3ab), which is in agreement with a poroelastic analysis presented by Guayacan-Carillo et al. (2016). With time and in the direction of the initial minor stress, the peak of pressure moves towards the interior of the rock mass while increasing in intensity at first and then decreasing the value of the peak of pore pressure with time. This remains in agreement with the in situ measurements (Seyedi et al. 2016). The apparition of overpressure is explained by elastic anisotropic behaviour of rock as described by Guayacan-Carillo et al. (2016) or Fig. 3a. However, the authors showed the overpressure reaches the maximum value at the instantaneous (0+) response of excavation in the vicinity of the drift wall and this overpressure disappears between 1 to 6 months. When the mechanical behaviour of the COx is elastoplastic (as the case of the present study), the peak of pore pressure is reached at 15 days after the excavation and the overpressures vanish after 2 years. This is closer to the in situ observations