Modelling of Predictive Hydraulic Impacts of a Potential Radioactive Waste Geological Repository on the Meuse / Haute-Marne Multilayered Aquifer System ( France )

The French National Agency for Nuclear Waste Management (Andra) conducted a site investigations program within the project of a deep geological disposal of radioactive waste in the Meuse/ Haute-Marne region. The construction of the tunnel of 5 Km length and the shafts of about 500 m depth to access the repository located in the clay host formation of Callovo-Oxfordian age, will lead to the perturbations of the groundwater flow fields. The prediction of the behaviour of these perturbations is needed to support: 1) the engineering and monitoring operations, and 2) the assessment of the consequences on groundwater resources. A variably-saturated flow model of a local multi-layered aquifer system is developed. It integrates the Oxfordian aquifer (limestone), the Kimmeridgianaquitard (marl) and the Barrois limestone aquifer including the karst conduits network. The variably-saturated flow Richard’s equation is solved with a finite element simulator. Prior to the simulation of the predictive repository impacts, a transient flow model is calibrated with respect to Underground Research Laboratory (URL) construction data. The results are analysed and evaluated by the use of performance measures.


Introduction
Prior to the construction and the operation, by Andra (French Agency for Radioactive Waste), of the future deep geological repository for high and intermediate long lived radioactive waste, an analysis and evaluation of predictive hydraulic impacts which could be induced by the underground structures (shafts and tunnels) linking surface installations to the rad waste disposal located in the clay formation of Callovo-Oxfordian age (160 My) at depth of 500 m.The clay host formation of Callovo-Oxfordian age is found throughout the multilayered sedimentary fill of the Paris basin.The Callovo-Oxfordian layer is located between an overlying limestone of Oxfordian (Jurassic) age and an underlying limestone of Dogger (Jurassic) age.The potential repository emplacement is located in the Meuse/Haute-Marne sector which includes two nested zones: 1) transposition zone ZT and 2) the ZIRA.The transposition zone has been chosen as a suitable host domain based on 20 years of site investigations results; it extends over an area of approximately 250 km 2 .In the transposition zone, the Callovo-Oxfordian formation is encountered at a depth of roughly 500 m, with a minimum thickness of approximately 130 m and a very small hydraulic conductivity (less than 10 −12 m/s).Andra has carried out several investigation campaigns considering the chosen host domain and its surroundings and also carried out extensive site descriptive modelling.The aim of the site-descriptive modelling is to develop a discipline-integrated description of the past and present conditions at the site, by analysing, assessing, and modelling the data obtained from the investigation campaigns."ZIRA" is an area extending over 30 Km 2 selected for further investigations, will be emplaced, five shafts with a diameter of 5 m to access the repository.Two access tunnels represented in the model by one ramp of over 5 Km length, 13 m of diameter and with a slope of about 12.This ramp will join the shafts in the host formation of the Callovo-Oxfordian at about 500 m of depth.

Geological Model, Engineered Structures and Data Monitoring
The hydrogeological conceptual model consists of 27 layers from the Triassic to Tertiary at the scale of the Paris basin and the 28 layers of the Triassic to Neocomian/Berriasian on the sector scale.The geological conceptual model represents the multilayered system including regional local faults and fracturing zone mainly located south west of the "ZIRA" which is the potential emplacement for the future repository.The construction of this geological concept was based on over 3300 drill holes data and about 2800 Km of seismic profiles.
Modelled domain is extracted from an integrated hydrogeological region-sector model which consists of the 1) 27 layers from the Triassic to the Quaternary on the scale of the Paris basin and 2) the 28 layers from the Triassic to the Portlandian on the scale of the sector.Lateral extend covers an area of about 3000 Km 2 which large enough to contain the hydraulic perturbation of repository construction.Vertically, only the geological layers above the Callovo-Oxfordian host formation which will host the underground structures (shafts and access tunnels) are represented.The first formation overlying the Callovo-Oxfordian is the Oxfordian formation considered as an aquifer and consists of limestones of Oxfordian age and early Kimmeridgian age and with a mean thickness at the repository potential emplacement of about 290 m.In the area of interest, this aquifer is split into two aquifer units separated by a marl layer called "sériegrise".Macro-pores zones identified through drill holes data and seismic survey results were characterized as relatively water productive with hydraulic conductivity ranging from 10 −9 m/s to 10 −7 m/s.Marne and limestone hydraulic conductivity are of 10 −12 m/s and 10 −9 m/s to 10 −10 m/s respectively.
The 7 identified macropores zones are represented in the model by four relatively productive units: • Hp1-4: lower unit located in deep Oxfordian with a mean thickness of 55 m; • Hp5 confined by the marl; • Hp6 and Hp7 belong to the Upper Oxfordien.Hp5, Hp6 and Hp7 have a mean thickness of about 5 m.
Kimmeridgianmarl layer overlying Oxfordian has men thickness of about 100 m and hydraulic conductivity estimated to 10 −12 m/s.Measured specific storage values of the macro-pores zones ranges from 2.3 × 10 −6 m −1 to 3.0 × 10 −6 m −1  For the simulations of the hydraulic impact of the access shafts to the Underground Research Laboratory (URL), it has been supposed that the perturbations of the flows will not affect the host formation northe underlying layers.It has thus been decided to uniquely model the layers on top of the Callovo-Oxfordian as a surface that comprises the sedimentary structures of the Oxfordian (the porous horizons and gray marl series).

Mathematical Formulation and Numerical Model Variably Saturated Flow
Variably saturated groundwater flow is described by the modified form of Richard's equation which can be solved by Groundwater computer code in pressure head ( ) ψ form and in water content/hydraulic head form ( ) relative permeability of the medium with respect to degree of saturation [-]; The empirical law of dependence of the relative permeability and the pressure on the saturation is generally based on the experience achieved from the soil columns (granular porous media) on the decimeter to meter scale.A sensitivity of the variably saturated flow on the parameters of the model of Mualem (1976) [1]-Van Genuchten (1978) [2] allowed to test the parameters controlling the dependence of the relative permeability on the pressure or on the saturation with the aim to guarantee a rapid numerical convergence.The model of Mualem (1976) [1]-Van Genuchten (1978) [2] describes the dependence between the relative permeability and the saturation according to: and n [-] are fit parameters to experimental results; ν : pore connectivityparameter [-] equal to 0.5 for most soils [Mualem  1976].
The controlling coefficients of this model are α [1/L] and n [-], the coefficient υ being generally equal to 0.5.With w S , m S , r S and ( ) ( )

S S S S S = − −
: the saturation with water, the maximum saturation, the residual saturation and the effective saturation, respectively.The parameter α could be considered as being the inverse of the thickness of the capillary zone.A low value of α is synonymous with a weak gradient of satura- tion, favorable for easy numerical convergence.
Figure 1 represents the range of van Genuchten curves utilized for parameterization of the unsaturated zone in the model under development.A α of 0.1 has been used for the semi-permeable formations so that they are twice as high as the aquifer formations (e.g.: the porous horizons).

• Meshing
The multilayer conceptual model described above has been discretized in an unstructured 2D and 3D finiteelement mesh with a view to simulate flow and transport.The 3D mesh comprises 6'090'802 elements; the total number of nodes is 3'138'876.The major faults are discretized by means of 2D finite elements representing the fault planes.The diffuse fracturing is discretized as 3D elements and is represented as an equivalent porous medium.The flow solutions have been produced by means of the calculation code Ground Water (Cornaton, 2007) [3].

• Boundary conditions
In accordance with the hydrodynamics of regional hydrogeological system of the Paris Basin, boundary conditions for flow of the local model are applied to the geological medium or geosphere and to the engineered structures while they are in construction and later during the operational phase of the Underground Research Laboratory (URL) and the future Radwasterepository.Two types of boundary conditions are used to simulate first a steady state flow within the multilayered system before the construction of URL shafts in 2001 and prior to start the transient flow simulations.These boundary conditions are: 1) of Neumann type and correspond to the specified flux or recharge on the top of the model and 2) of Dirichlet type: constant hydraulic heads are assigned on the trace of the hydrographic network and at the lateral boundaries of the model.For the transient flow, the construction of the access shafts for the URL requires assignment of special boundary conditions that are variable in time.The boundary conditions applied at the level of the access shafts are seepage conditions: a hydraulic head equal to the elevation is attributed to the node (the pressure will be equal to the atmospheric pressure) while the flux is outward.Otherwise a null flux is assumed.The seepage conditions are variable in time in the sense that they are activated according to the calendar of advancement of the excavation.

Results of the Groundwater Flow Simulations
The analysis of the pressure records in multiple observation points shows that a quasi-steady state is obtained.This result has motivated an initial step of calibration for the local model which has consisted of calibrating the hydraulic conductivities of the different modeled units in the steady-state regime, supposing an instantaneous excavation of the main access shaft (PPA) to the URL and secondary shaft (PAX) where a seepage condition has been specified.The objective was to adjust the initial condition prior to excavation of the subterranean works but particularly calibration of the drainage flow rates for an equilibrium state.The values of hydraulic conductivity resulting from this calibration, along with the specific storage [coefficients] coming from the interpretation of hydraulic tests constitute the initial parameters for a transient-regime calibration.The calibration of the local model aims to better reproduce the measurements of hydraulic head in the monitoring drill holes of the perturbations and the flow rates of the water arrivals registered in the shafts PPA and PAX from October 2001 through August 2012, thus 11 years of continuous measurements.The method of calibration by "trial and error" is carried out by successive simulations and observation of the fit between the observed and modeled data.In each new simulation the user adjusts the hydrodynamic parameters to improve the fit of the measured values and those calculated.Computer code Ground Water (Cornaton, 2012) has been utilized to carry out the flow simula-tion for a time period of 3957 days.The march in time is managed by an automatic estimation method with optimal time steps (predictor-corrector time integration algorithm).The time steps thereby are increased or reduced according to the convergence between two time steps.The calibration seeks to define the values of the parameters K and S of the 7 hydro geological entities that better reproduces the temporal records of the heads and flow rates in the porous horizons.One establishes for this purpose an objective function that represents the variable of interest.This function is the relative error squared between the measured data and the numerically simulated data.14 records of head and 8 records of flow rates are used to calibrate the observed transient flow.The Objective Function is calculated for each simulation.The response surfaces method (RSM) allows establishing a mathematical approximation model from the values of the Objective Function.This polynomial model substitutes for the numerical model and thus allows evaluating values of the Objective Function for different combinations of the parameters.The RSM model is constructed by polynomial multiple regression based on the least squares criterion [4].

Concluding Remarks
The results of predictive hydraulic impact of the underground structures of the future rad waste repository seek to be improved by reducing the conceptual and numerical uncertainties.One of the scopesis to compute the transient evolution of the unsaturated zone around the tunnel and the shafts under construction and during the operational phase.This will be achieved by: 1) refining the hydro geological concept of Barrois limestone aquifer system, 2) coupling the surface and subsurface flow and finally 3) refining the mesh around the tunnel and shafts in order to better control the unsaturated zone expansion during the construction and exploitation of the repository and its retreat induced by the closure of the facilities.The hydro-mechanical effect on the hydraulic overall hydraulic impact of the ramp and shaft construction will be investigated.

Figure 1 .
Figure 1.Range of van Genuchten curves utilized in the parameterization of the unsaturated zone.In red for α =0.1 [1/m] and n = 1.8 for the the semi-permeable units, in blue for α = 0.2 [1/m] and n = 3 for the aquifer units.The residual saturation of 0.1 is also shown in gray.

Figure 2
represents the calibration of the transient flow based on the matching of the observed discharge and pressure/head time series.The predictive simulations of the hydraulic impacts rely on the quality of the transient flow calibration.

Figure 3 .
Figure 3. Meshing, layout of shafts and tunnels, and hydraulic impact at 100 years after the repository construction.