Journal of Applied Mathematics and Physics, 2014, 2, 1085-1090
Published Online November 2014 in SciRes. http://www.scirp.org/journal/jamp
http://dx.doi.org/10.4236/jamp.2014.212125
How to cite this paper: Hakim, B., Jaouher, K., Laurent, T., Gregory, D. and Pierre, P. (2014) Modelling of Predictive Hydrau-
lic Impacts of a Potential Radioactive Waste Geological Repository on the Meuse/Haute-Marne Multilayered Aquifer Sys-
tem (France). Journal of Applied Mathematics and Physics, 2, 1085-1090. http://dx.doi.org/10.4236/jamp.2014.212125
Modelling of Predictive Hydraulic Impacts
of a Potential Radioactive Waste
Geological Repository on the
Meuse/Haute-Marne Multilayered
Aquifer System (France)
Benabderrahmane Hakim1*, Kerrou Jaouher2, Tacher Laurent3, Deman Gregory2,
Perrochet Pierre2
1Research & Dev elop men t De par tment , ANDRA, 1-7 rue Jean Monnet, Châtenay-Malabry Cedex, France
2Centre of Hydrogeology and Geothermics (CHYN), University of Neuchâtel, Rue Emile Argand 11, Neuchâtel,
Switzerl and
3Laboratoire de Mécanique des Sols (EPFL-LMS), Ecole Polytechnique Fédérale de Lausanne, Lausanne,
Switzerl and
Received August 2014
Abstract
The French National Agency for Nuclear Waste Management (Andra) conducted a site investiga-
tions 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 as-
sessment of the consequences on groundwater resources. A variably-saturated flow model of a lo-
cal 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 ana-
lysed and evaluated by the use of performance measures.
Keywords
Numerical Modelling, Variably-Saturated Flow, Karst, Repository, Hydraulic Impacts
*
Corresponding author.

B. Hakim et al.
1086
1. 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 pre-
dictive 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 sedi-
mentary fill of the Paris basin. The Callovo-Oxfordian layer is located between an overlying limestone of Ox-
fordian (Jurassic) age and an underlying limestone of Dogger (Jurassic) age. The potential repository emplace-
ment 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 inves-
tigations results; it extends over an area of approximately 250 km2. In the transposition zone, the Call-
ovo-Oxfordian formation is encountered at a depth of roughly 500 m, with a minimum thickness of approxi-
mately 130 m and a very small hydraulic conductivity (less than 10−12 m/s). Andra has carried out several inves-
tigation 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. “ZIR A” is an area extending over 30 Km2 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 Cal lo vo -Oxfordian at about 500 m of depth.
2. Conceptual Model and Mathematical Formulation
2.1. 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 “ZI R A” which is the potential emplacement for the future repository. The construction of this geo-
logical 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 Trias-
sic to the Portlandian on the scale of the sector. Lateral extend covers an area of about 3000 Km2 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 tun-
nels) are represented. The first formation overlying the Callovo-Oxfordian is the Oxfordian formation consi-
dered 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 under-
lying 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).

B. Hakim et al.
1087
2.2. 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
( )
-H θ
.
qQ
t
θ ∂+∇⋅ =∆
∂
r
qkK H=−∇
( )()
r
ckKz Q
t
ψ
ψψ
∂=∇⋅∇+ +∆
∂
( )
w
ws
S
c SS
ψθ
ψ
∂
=+ ∂
( )
r
kKz Q
t
θψ ∂=∇⋅∇+ +∆
∂
[ ]
LHz ψ = +
,
ψ
: pressure head
[ ]
L
; z: elevation
[]
L
, q: Flux vector
1
LT
−


; t: time
[]
T
;
Q∆
: source and/or
sink term
1
T−


;
K
: Tensor of saturated hydraulic conductivity
1
LT
−


;
θ
:water content [-];
( )
r rw
k kS=
:
relative permeability of the medium with respect to degree of saturation [-];
ws
S θθ =
[-];
s
θ
: saturated water
conte nt [ -];
( )
c ψ
: nonlinear capacitive coefficient.
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:
( )
1 for 0
m
n
w rmr
SS SS αψ ψ
−

=++ +<

for 0
wm
SS ψ = ≥
( )
2
1
1
11, 1, 1
m
m
re e
kS Smn
n
ν

= −−=−>


wr
emr
SS
SSS
−
=−
r
S
: residualsaturation [-];
m
S
: maximum saturation [-];
e
S
: effective saturation [-];
α
[L−1] 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
() ()
ewr mr
SSSSS=−−
: 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).

B. Hakim et al.
1088
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.
• Meshing
The multilayer conceptual model described above has been discretized in an unstructured 2D and 3D fini-
teelement mesh with a view to simulate flow and transport. The 3D mesh comprises 6’090’802 elements; the to-
tal 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 con-
ditions 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 varia-
ble 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.
3. 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 conduc-
tivity 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 perturba-
tions 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 car-
ried 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-

B. Hakim et al.
1089
tion for a time period of 3957 days. The march in time is managed by an automatic estimation method with op-
timal 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 parame-
ters 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 da-
ta. 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 mathe-
matical 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]. 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.
4. Results of the Predictive Impact Simulations
Predictive impact simulations were conducted based on the transient flow model reproducing the hydraulic per-
turbations induced by the URL construction and the drained shafts during 12 first years of the operational phase.
Evaluation of the zone of hydraulic perturbation due to excavation of the access shafts has been carried out in the
steady-state regime as well as in transient regime. The results have shown that the perturbation might at maximum
exceed the extent of the transposition zone (steady-state regime). In the transient regime, the zone affected by the
work of excavating the underground openings is less extensive with a maximum perturbation at the level of the
lower Argovian around the tunnel/ramp [Figure 3]. After utilization of [the facility for] 100 years, it is expected to
close and completely seal the underground infrastructure of the future repository. The predictive simulations of
transient flow have been carried out to evaluate the time needed to return to a state of hydraulic equilibrium [as]
prior to the perturbation. Tunnel and shafts construction induce hydraulic disturbance which is added to the one
in place originated fromthe URL construction period and its present activity. Extensions and amplitudes of the
hydraulic impact gradually evolve according to the schedule of the works. After six months of excavation of the
first two shafts, it follows: 1) in the upper Oxfordian, a high amplitude perturbation localized around the tunnel
while the shafts do not yet generate significant disruption, and 2) in the middle Oxfordian, the homogenized
disturbance includes impacts of the tunnel and the shaft. At the end of construction, impacts: 1) in the upper
Oxfordian indicate a local homogeneous perturbation with an amplitude of 40 m around the tunnel, disconti-
nuous and very low extension and amplitude around the four shaft and 2) in the Middle Oxfordian a ho mo-
Figure 2. Measured and computed discharge (shafts) and pressure (monitoring borehole) time series.

B. Hakim et al.
1090
Figure 3. Meshing, layout of shafts and tunnels, and hydraulic impact at 100 years after the repository construction.
geneous disturbance of tunnel and shafts with maximum amplitude of about a hundred meters; at this stage the
disturbance remains separated from that induced by underground laboratory.
The weak hydraulic response to shaft in the Upper Oxfordian is due to the absence of the porous horizons
HP6 Hp7 in this area which is characterized by relatively low hydraulic diffusivity. During the operational phase:
1) in the Upper Oxfordian: disturbances of the tunnel and shafts as well as the Underground Laboratory are
merged 50 years approximately; at 100 years approximately, the overall disturbance covers the South-West part
of the transposition zone [Figure 3] and 2) in the Middle Oxfordian, the interference of the hydraulic distur-
bances occurs at the end of the construction phase. The overall impact covers the south-west part of the Trans-
position with maximum amplitude of 200 m.During the construction phase the maximum outflow rate of drained
groundwater into the shafts, is about 20 m3/day/shaft while the tunnel maximum water production is estimated to
150 m3/day. Tunnel length crossing Oxfordian aquifer is of 3800 m. At the end of the operational phase (100
years approximately) of the repository, the water discharge developed by the macro-pores (Hp1-4, Hp5, Hp6) is
about 4 m3/day while the tunnel discharge is one order of magnitude greater.
5. 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 tran-
sient evolution of the unsaturated zone around the tunnel and the shafts under construction and during the opera-
tional 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 reposi-
tory 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.
References
[1] Mu alem, Y. (1976) A New Model to Predict the Hydraulic Conductivity of Unsaturated Porous Media. Water Re-
sources R esearch , 12, 513-522. http://dx.doi.org/10.1029/WR012i003p00513
[2] Van Genuchten, M.Th. (1978) Mass Transport in Satur ated-Unsaturated Media: One-Dimensional Solutions. Research
Rep. No. 87-WR-11, Water Resources Program, Princeton University.
[3] Cornaton, F. (2007) Ground Water: A 3-D Ground Water Flow, Mass Transport and Heat Transfer Finite Element Si-
mulator. Reference Manual, 363 p. http://www1.unine.ch/chyn/php/softwares.php
[4] Kerrou, J., Deman, G. , Tach er, L., Ben abd errahmane, H. and P erro chet, P. (2015) Calibrating a 3D Numerical Model
of V ariab ly-Saturated Flow to Assess the Hydraulic Impacts of the Future Radwaste Repository on the Meuse/Haut e-
Marne S ite. The Clay Conference, Brussels.