# Assessment of a high-level spent nuclear fuel disposal model

A. Rodriguez-Dono, S. Olivella, N. Mokni

<figure><img src="image_1.png" /></figure>

## To cite this version:

A. Rodriguez-Dono, S. Olivella, N. Mokni. Assessment of a high-level spent nuclear fuel disposal model. Environmental Geotechnics, 2018, 7 (1), pp.42-58. 10.1680/jenge.18.00017. hal-02881791

HAL Id: hal-02881791

https://hal.science/hal-02881791

Submitted on 19 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 **Assessment of a high-level spent nuclear fuel disposal model**

2

3 Alfonso Rodriguez-Dono<sup>a,*</sup>, Sebastia Olivella<sup>a</sup>, Nadia Mokni<sup>c</sup>

4 <sup>a</sup> Department of Civil and Environmental Engineering, Polytechnic University of

5 Catalonia (UPC), Barcelona, Spain

6 <sup>b</sup> Institute of Environmental Assessment and Water Research (IDAEA), Spanish

7 National Research Council (CSIC), Spain

8 <sup>c</sup> Institute for Radiation Protection and Nuclear Safety (IRSN), France

9 * Corresponding author. Address: D2-306/2 UPC Campus Nord. Jordi Girona 1-3,

10 08034 Barcelona, Spain. E-mail: alfonso.rodriguez@upc.edu

11

12 **Abstract**: This paper assesses a model for the full-scale engineered barrier experiment

13 (Febex) in situ test, simulating the disposal of heat-emitting, high-level spent nuclear fuel.

14 The model has been developed using Code_Bright and is based on previous modelling

15 efforts. The model focuses on the thermo-hydro-mechanical behaviour of a bentonite

16 barrier surrounding a heater that simulates the heat production of radioactive spent fuel.

17 The new model incorporates the new dismantling operations and has been compared with

18 new in situ measurements. Once calibrated, the long-term response of the model has been

19 analysed. Furthermore, the implementation of additional processes and their impact on

20 hydration and stress development have been assessed. In this assessment, it has been

21 found that the bentonite barrier does not fully saturate near the canister in any case, which

22 is an important goal for the isolation of the spent nuclear fuel. The intrinsic permeability

23 has been found to give results closer to measurements when its value is doubled. The

24 double-structured model has given interesting results, but the authors found that it does

25 not have a great influence in the general behaviour of the model.

26 **Keywords**: FEBEX, CODE_BRIGHT, numerical modeling, nuclear waste disposal.

27

28 **1. Introduction**

29 Deep geological storage or disposal remains the preferred option for waste management

30 of heat-emitting, high-level radioactive nuclear waste (HLW) in several countries.

31 Excavation of deep underground repositories is limited to rock units that are reasonably

32 stable and without major groundwater flow, and to depths of a few hundred metres

33 below the surface.

34 The required degree of waste isolation needed for HLW is provided by a combination of

35 engineered and natural barriers placed between the potentially harmful radionuclides

36 and the biosphere. The natural barrier is fundamentally the host rock, and the artificial

37 barriers are the solid matrix of the waste itself, the metallic canister enclosing the waste,

38 and the backfill—the sealing material placed around the canisters. This backfill, also

39 called an 'engineered barrier', is often constructed using compacted expansive clay.

40 Bentonite has generally been chosen because of its high swelling capacity, low

41 permeability and favorable retardation properties (Gens et al., 2009).
The bentonite barrier fulfills several important functions. In the first instance, a very low hydraulic conductivity restricts water penetration and retards significantly solute transport due to its low diffusion coefficient and to additional sorption effects. It should also provide a favorable chemical environment and be able to self-heal if subjected to physical perturbation such as cracking and fissuring events (Gens, 2003). The engineered clay barrier and adjacent host rock (usually called the 'near field') will be subjected to the heating effect of the nuclear waste, and also to various associated hydraulic and mechanical phenomena that interact in a complex way.

In addition, compacted bentonite is initially unsaturated, and will therefore be subjected to hydration from the surrounding rock, triggering further coupled thermo-hydro-mechanical (THM) phenomena. In order to achieve a safe and robust repository design, it is necessary to have a good understanding of the processes that occur in the near field, and of their evolution over time (Gens et al., 2009).

Given the importance of the problem, it has proved useful to perform large-scale and medium-scale heating tests simulating repository conditions in underground laboratories around the world (Dixon et al., 2002; Pusch et al., 1985; Selvadurai, 1997; Volckaert et al., 1996). Because of the low permeability of the materials involved, such experiments usually require long testing times, measured in years, in order to obtain meaningful results. The complexity of the phenomena involved, and of their interactions, implies that process understanding and interpretation of results are limited unless supported by suitable numerical models able to reproduce the main features of the test. Owing to the coupled nature of the THM interactions, coupled THM formulations are inevitably required.

The present paper contains the description of an updated model for the large scale in situ heating test FEBEX. In this test, after five years of heating, one of the heaters was switched off and the experiment was partially dismantled, allowing the final state of the barrier to be observed directly. In this way, very valuable information on the state of the bentonite at the end of the test was obtained. The test has received attention during the initial (Gens et al., 1998) and intermediate stages (Alonso & Alcoverro, 2005). Finally, Gens et al. (2009) discussed the thermal, hydraulic and mechanical observations in the bentonite barrier and in the host rock, paying special attention to the progress of hydration in the barrier, the effects of heating and vapor transport, and the development of swelling pressures in the barrier.

In the FEBEX experiment, heaters are emplaced in the axis of a tunnel excavated in granite to simulate the heat production of radioactive waste. The test is fully instrumented, and attention is focused on the thermo-hydro-mechanical (THM) behavior of the near-field region constituted by the compacted bentonite barrier surrounding the heater and the immediately adjacent rock. Interpretation of the test is assisted by the performance of a coupled numerical analysis based on a formulation that incorporates the relevant THM phenomena. Further description of the FEBEX in situ test can be found in Gens et al. (2009).

The model of this test has been developed using CODE_BRIGHT (Olivella et al., 1996), a Finite Element Method program that enables coupled thermo-hydro-mechanical analysis in geological media. This is a code version that incorporates customized pre- and post- process interfaces so that simulation models that use CODE_BRIGHT –available from: https://deca.upc.edu/en/projects/code_bright- can be developed more easily. It allows easy modification of parameters, boundary conditions, excavation protocols, meshing and organization of calculation intervals. The reason is
that it has been developed using the GiD interface (http://www.gidhome.com) that
permits pre- and post- process of data in a user-friendly way.

The updated model for FEBEX is based on the geometry and material properties of the
previous modeling efforts (Gens et al., 2009). The main objective of this work is to
develop a new model of the well-known FEBEX experiment that can be further
exploited in terms of two main goals: (1) comparison with new measurements and (2)
implementation of additional processes such as coupled flows and double structure
models to investigate their impact of hydration and stress development.

This updated model is intended to be complementary with the work done by Sanchez et
al. (2012) who also made long term predictions about the FEBEX model but didn't
contemplate the introduction and assessment of the additional processes that are
included in this article.

## 2. The FEBEX updated model

Figure 1 shows the model geometry and materials considered in this model. The host
rock, the bentonite buffer, and the concrete plug are the geomaterials considered. In
addition, a steel liner confines the bentonite and plug. The model assumes axisymmetry
along the tunnel axis. The geometry and material distribution are prepared to simulate
the partial dismantling of the test which was carried out during the performance of the
test in such a way that one of the heaters has continued operating. Boundary conditions
are modified accordingly for the excavation of the plug and part of the bentonite is
removed from the model thus simulating excavation.

This model takes advantage of the unstructured mesh development for the discretization
of the host rock. However, the buffer is discretized with structured mesh as this is more
adequate in the zone of interest to reduce numerical errors and to facilitate the
monitoring of calculated variables that should be compared with measurements.

Table 1 shows the intervals considered in the analysis up to 2500 days, i.e. before
dismantling operations. In order to match the measured evolution of temperature during
the experiment (see more details about this in Gens et al., 2009), the heat power was
progressively modified (see table 1) until the 100°C target temperature was reached.
From that point on, the temperature was prescribed at the value of 100°C in both
heaters. Later, one of the heaters was removed but the other continued heating.

A mechanical boundary condition restraining the normal displacements has been
applied all along the external boundary. Also, a constant pressure of 0.9 MPa is imposed
on the external boundary. As indicated above, during intervals 3 to 5 the power of the
heaters is prescribed and therefore the temperature increases. Once the temperature
reaches the target value on the surface canister (100°C) the boundary condition on the
heater is changed from constant power to constant temperature (more information about
boundary condition implementation in the Code_Bright User's Guide, downloadable
from the Code_Bright web page – https://deca.upc.edu/en/projects/code_bright).

Tables 2 and 3 show the constitutive laws used in the model and the properties
corresponding to the FEBEX bentonite as considered with single porosity and
considering that it can be represented with the Barcelona Basic Model or BBM (Alonso
et al., 1990), van Genuchten model, Darcy's law, power relative permeability and
geometric mean for thermal conductivity. The parameters are divided in hydraulic and
mechanical parameters and correspond to the Code_Bright materials input window
136 (more information about this parameters in Code_Bright User's Guide). More precise
137 information may be obtained from the model, which is available by request to the
138 authors.

<figure><img src="image_3.png" /><figcaption>139
140 Figure 1. Model geometry, materials and mesh considered for FEBEX updated
141 modeling (some materials are differentiated to simulate excavation only).</figcaption></figure>

<table><caption>142
143 Table 1. Time intervals considered for the simulation up to 2500 days</caption><thead><tr><th>Time interval</th><th>Initial Time (days)</th><th>Final Time (days)</th><th>Prescribed heat flow (W/heater)</th><th>Prescribed Temperature (°C)</th></tr></thead><tbody><tr><td>1 (stress equilibrium)</td><td>0</td><td>20</td><td>-</td><td>12</td></tr><tr><td>2 (construction)</td><td>20</td><td>155</td><td>-</td><td>12</td></tr><tr><td>3</td><td>155</td><td>176</td><td>1200</td><td>-</td></tr><tr><td>4</td><td>176</td><td>197</td><td>2000</td><td>-</td></tr><tr><td>5</td><td>197</td><td>207</td><td>2400</td><td>-</td></tr><tr><td>6</td><td>207</td><td>2500</td><td>-</td><td>100</td></tr></tbody></table>

144
145 Figure 2a shows the distribution of temperature around the 2 heaters at 21 days, 1.5
146 years and 5.5 years, i.e. during normal operation of the test and thus before dismantling.
147 It can be observed that the temperature in the vicinity of the heaters rises till it reaches a
148 value close to 100°C and the temperature of the surrounding host rock continues
149 increasing after that.

150 Figure 2b shows the distribution of degree of saturation. In general, the rock remains
151 saturated except at early times. As the bentonite is installed in the tunnel with a high
152 initial suction, some desaturation of the rock takes place. This zone is re-saturated very
153 soon as the water from the host rock is mobilized. At 5.5 years, the bentonite buffer is
154 not fully saturated, but the drying near the canister is less intense than at 1.5 years,
155 indicating that the saturation process continues.

156 Figure 2c shows the distribution of porosity in the bentonite. Clearly, swelling of the
157 bentonite takes place near the hydration boundary (i.e. the contact with the rock) while
158 contraction takes place in the zone of drying (i.e. the contact with the heaters). As it can
159 be seen in the figure, the range of porosity variation after 5.5 years is 0.378 to 0.413,
160 with an average very close to the initial value of 0.4.

161

162 Table 2. Hydraulic parameters for bentonite. See Code_Bright User's Guide for further
details.

Retention curve

<table>
  <tr>
    <td>P1: P<sub>o</sub> (MPa)</td>
    <td>20</td>
  </tr>
  <tr>
    <td>P2: &sigma;<sub>o</sub> (Nm<sup>-1</sup>)</td>
    <td>0.072</td>
  </tr>
  <tr>
    <td>P3: &lambda;</td>
    <td>0.18</td>
  </tr>
  <tr>
    <td>P4: S<sub>rl</sub></td>
    <td>0.01</td>
  </tr>
</table>

Van Genuchten model:

$$
S_e = \frac{S_l - S_{rl}}{S_{ls} - S_{rl}} = \left[ 1 + \left( \frac{P_g - P_l}{P} \right)^{\frac{1}{1-\lambda}} \right]^{-\lambda} \quad P = P_o \frac{\sigma}{\sigma_o}
$$

Intrinsic permeability

<table>
  <tr>
    <td>P1: (k<sub>11</sub>)<sub>o</sub> (m<sup>2</sup>)</td>
    <td>1.9E-21</td>
  </tr>
  <tr>
    <td>P2: (k<sub>22</sub>)<sub>o</sub> (m<sup>2</sup>)</td>
    <td>1.9E-21</td>
  </tr>
  <tr>
    <td>P3: (k<sub>33</sub>)<sub>o</sub> (m<sup>2</sup>)</td>
    <td>1.9E-21</td>
  </tr>
  <tr>
    <td>P4: &phi;<sub>o</sub></td>
    <td>0.4</td>
  </tr>
  <tr>
    <td>P5: &phi;<sub>min</sub></td>
    <td>0.001</td>
  </tr>
</table>

Darcy's law: $\mathbf{q}_l = - \frac{\mathbf{k} k_{rl}}{\mu_l} (\nabla P_l - \rho_l \mathbf{g})$

Kozeny's model: $\mathbf{k} = \mathbf{k}_o \frac{\phi^3}{(1-\phi)^2} \frac{(1-\phi_o)^2}{\phi_o^3}$

Liquid phase relative permeability

<table>
  <tr>
    <td>ITYCL</td>
    <td>6</td>
  </tr>
  <tr>
    <td>P2: A</td>
    <td>1</td>
  </tr>
  <tr>
    <td>P3: λ</td>
    <td>3</td>
  </tr>
</table>

$k_{rl} = A S_e^{\lambda}$

Diffusive flux of vapor

<table>
  <tr>
    <td>ITYCL</td>
    <td>1</td>
  </tr>
  <tr>
    <td>P1: D (m<sup>2</sup>s<sup>-1</sup>K<sup>-n</sup>Pa)</td>
    <td>5.9E-6</td>
  </tr>
  <tr>
    <td>P2: n</td>
    <td>2.3</td>
  </tr>
  <tr>
    <td>P3: &tau;<sub>o</sub></td>
    <td>0.8</td>
  </tr>
</table>

Fick's law for molecular diffusion:

$$
\mathbf{i}_g^w = -(\tau \phi \rho_g S_g D_m^w \mathbf{I}) \nabla \omega_g^w D_m^{\text{vapor}} = D \left[ \frac{(273.15+T)^n}{P_g} \right]
$$

$\tau$ = constant = $\tau_o$

Conductive flux of heat 1

<table>
  <tr>
    <td>ITYCL</td>
    <td>1</td>
  </tr>
  <tr>
    <td>P1: λ<sub>dry</sub> (WmK<sup>-1</sup>)</td>
    <td>0.47</td>
  </tr>
  <tr>
    <td>P2: λ<sub>sat</sub> (WmK<sup>-1</sup>)</td>
    <td>1.15</td>
  </tr>
</table>

Fourier's law: $\mathbf{i}_c = -\lambda \nabla T$

$$
\lambda = \lambda_{sat}^{S_l} \lambda_{dry}^{(1-S_l)}
$$

164

165 Figure 3 shows the variation of temperature at some points within the bentonite and the
166 comparison between the Operational Base Case (OBC) model calculations –
167 corresponding to the analysis reported in Gens et al. (2009)– and the current model. It
168 can be seen that the temperature is reproduced in the same way by the model as it was
169 obtained with the OBC model, so this new model is equivalent to that one. It is
170 important to recall that one of the objectives of this contribution was to prepare a new
171 geometry and mesh (GiD environment) to be able to consider additional cases.

172 A similar response is obtained in terms of relative humidity (Figure 4) and the
173 correspondence of the present model with the OBC model shows that no significant
174 variation has occurred. Actually, small variations on the parameters or shapes of the
175 constitutive equations do not produce a significant variation of relative humidity. This
176 shows that it is possible that some additional mechanisms of water flow (in addition to
177 Darcy and Fick equations) may play a role and may help to explain the different
178 hydration rates.

179 The evolution of stresses in the bentonite is represented in Figure 5. The response shows
180 the same limitations as observed in the OBC especially concerning the overestimation
181 of stress development at early stages of hydration. The low measured stresses can be
182 interpreted qualitatively in different ways, one of them being the possible hydration of
183 the bentonite aggregates controlled by double structure effects.

184 Table 3. Mechanical parameters for bentonite. Thermo-elasto-plastic (TEP) model. See
185 Code_Bright User's Guide for further details.

Elastic parameters

<table>
  <tr>
    <td>ITYCL</td>
    <td>1</td>
  </tr>
  <tr>
    <td>P1: k<sub>io</sub></td>
    <td>-0.05</td>
  </tr>
  <tr>
    <td>P2: k<sub>so</sub></td>
    <td>0.25</td>
  </tr>
  <tr>
    <td>P3: K<sub>min</sub> (MPa)</td>
    <td>0.1</td>
  </tr>
  <tr>
    <td>P5: v</td>
    <td>0.4</td>
  </tr>
  <tr>
    <td>P8: α<sub>i</sub></td>
    <td>-0.003</td>
  </tr>
  <tr>
    <td>P9: α<sub>sp</sub></td>
    <td>-0.161</td>
  </tr>
  <tr>
    <td>P10: p<sub>ref</sub> (MPa)</td>
    <td>0.01</td>
  </tr>
</table>

$$
d\epsilon_v^e = \frac{k_i(s) dp'}{1+e p'} + \frac{k_s(p', s) ds}{1+e (s+0.1)} + (\alpha_o) dT
$$

Thermal parameters

<table>
  <tr>
    <td>ITYCL</td>
    <td>1</td>
  </tr>
  <tr>
    <td>P1: α<sub>o</sub></td>
    <td>1.5E-4</td>
  </tr>
  <tr>
    <td>P5: T<sub>ref</sub> (°C)</td>
    <td>20</td>
  </tr>
</table>

where:

$$
k_i(s) = k_{io} (1 + \alpha_i s)
$$

$$
k_s(p', s) = k_{so} (1 + \alpha_{sp} \ln p' / p_{ref})
$$

Plastic parameters 1

<table>
  <tr>
    <td>ITYCL</td>
    <td>1</td>
  </tr>
  <tr>
    <td>P1: λ(0)</td>
    <td>0.15</td>
  </tr>
  <tr>
    <td>P2: r</td>
    <td>0.925</td>
  </tr>
  <tr>
    <td>P3: β (MPa<sup>-1</sup>)</td>
    <td>0.05</td>
  </tr>
  <tr>
    <td>P4: ρ (°C<sup>-1</sup>)</td>
    <td>0.2</td>
  </tr>
  <tr>
    <td>P5: k</td>
    <td>0.1</td>
  </tr>
  <tr>
    <td>P6: p<sub>so</sub> (MPa)</td>
    <td>0.1</td>
  </tr>
</table>

$$
p_o = p^c \left( \frac{p_o^*(T)}{p^c} \right)^{\frac{\lambda(o)-kio}{\lambda(s)-kio}}
$$

$$
p_o^*(T) = p_o^*
$$

$$
\lambda(s) = \lambda(o) \left[ (1 - r) \exp(-\beta s) + r \right]
$$

Plastic parameters 2

<table>
  <tr>
    <td>ITYCL</td>
    <td>1</td>
  </tr>
  <tr>
    <td>P1: p<sup>c</sup> (MPa)</td>
    <td>0.5</td>
  </tr>
  <tr>
    <td>P2: M</td>
    <td>1</td>
  </tr>
  <tr>
    <td>P3: α</td>
    <td>0.53</td>
  </tr>
  <tr>
    <td>P4: e<sub>o</sub></td>
    <td>0.6</td>
  </tr>
  <tr>
    <td>P5: p<sub>o</sub><sup>*</sup> (MPa)</td>
    <td>12</td>
  </tr>
</table>

$$
p_s = p_{so} + k_s \exp(-\rho\Delta T)
$$

$$
\Delta T = T - T_{ref}
$$

187 Finally, the evolution of water pressure in the rock is represented in Figure 6. This
188 model uses 0.9 MPa as a boundary condition in the outer boundary. From the
189 observation of the results (overestimated pressure in the rock roughly by 0.05 MPa), it
190 is probable that 0.85 MPa would be more appropriate. This will be modified in future
191 models but the effect on the other variables is expected to be very small taking into
192 account the large suction gradients in the bentonite.
193

194
195 **3. Long term results**

196 Once the presented model was calibrated with real measurements, a long term
197 prediction has been performed. This section contains some results corresponding to the
198 model considering a maximum calculation time of 20 years and the incorporation of the
199 operations of shut down and dismantling of one canister, concrete abutment and
200 bentonite. The distribution of temperatures is shown in Figure 7a. It can be observed
201 that the operations of dismantling affect the canister close to the access drift together
202 with part of the barrier, while the second canister heating continues. Figure 7b shows
203 the distribution of degree of saturation and Figure 7c the distribution of mean stresses.
204 Degree of saturation indicates that at 20 years, full saturation has not been achieved yet.

205 The evolution of temperatures is presented in Figure 8. Temperatures near the
dismantled heater show that the cooling takes place quite fast and that it has only a
206 small influence on section D1. In section F2 and G, the effect is very important and the
207 model captures well the rapid fall of the temperature rapid down to the measured lower
208 values due to the new configuration of the test after partial dismantling.
209

210 Pore pressures are represented in Figure 9. The negative pore pressures correspond
211 actually to suction. The strong evaporation induced by heating near the canister
212 produces a large increase of suction. The maintained heating in section F2 delays the
213 pore pressure increase near the canister. This is because the evaporation is strong and
214 vapor migration compensates the water flow from the rock. In section G, which
215 undergoes cooling caused by the canister removal, the water pressure recovers more
216 rapidly at the same time that cooling takes place. The rate of water pressure recovery
217 increases rapidly.

218 Relative humidity (HR) depends on pressure and temperature and can be calculated
219 from these variables at any point. This is necessary in order to compare the calculated
220 variables with the measurements of relative humidity at the sensors. The comparison is
221 shown in Figure 10. Note that there is not real data for r = 0.52 m at section F2. While
222 in section F2, near the heater, relative humidity does not reach 100%, which would
223 correspond to full saturation, the 100% relative humidity is reached after 20 years at
224 section H, indicating that when the heating stops, full saturation is reached more
225 rapidly.

226 Figure 10 also shows the evolution of stresses at sections F2 and G. Note that sections H
227 and G can be considered equivalent sections, since they are very close to each other (see
228 Figure 3) and since neither there is data for stress at section H nor there is data for
229 relative humidity at section G. The stress development shows also the effect of the
dismantling but it is less pronounced in the central zone of the prevailing heater.
230 Stresses reach a value in the range of 4 to 5 MPa according to the model calculations.
231 The correspondence with the measurements is quite good at F2. In contrast, at section G
232 in which the measurements show a delay in the response, the sensor seems to have been
233 reactivated when dismantling took place. This model will be considered the base case
234 for the additional processes assessment that will be described next.
235

236

237 **4. Additional processes assessment**

238 The model presented above will be considered from now on as a base case for the
239 comparison with the results from alternative models including additional processes such
240 as the incorporation of a coupled flow induced by a thermal gradient, the variation of
241 the intrinsic permeability or the introduction of a double porosity model, for which the
242 impact on hydration and stress development has been assessed.

243

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

244
<figure><img src="image_7.png" /><figcaption>Figure 2. Distributions at 21 days, 1.5 years and 5.5 years of: a) temperature; b) degree of saturation; c) porosity.</figcaption></figure>

a) Results from OBC case as compared with measurements

<figure><img src="image_8.png" /></figure>

b) Results from the current new model

<figure><img src="image_9.png" /><figcaption>Figure 3. Temperature evolution at different points.</figcaption></figure>
## 4.1. Thermal gradient induced water flow

The water flow induced by a thermal gradient can be expressed as an additional flux to the usual Darcy's flux. The generalized advective flux including thermal osmosis is expressed as in Eq. (1):

$$ \mathbf{q}_l = -\frac{\mathbf{k}_{rl}}{\mu_l} (\nabla p + \rho_l g \nabla z) - k_{rl} k_T \nabla T \quad (1) $$

Where $\mathbf{k}$ [m²] denotes the intrinsic permeability tensor, $k_{rl}$ [-] denotes the relative permeability of the liquid, $g$ is the gravity vector, $\mu_l$ is the dynamic viscosity, $p$ [MPa] is the hydraulic pressure, $T$ [K] is temperature, $\rho_l$ is liquid density, and $k_T$ [m²/K/s] is the thermo-osmotic permeability.

This is the flux of liquid phase induced by both pressure and temperature gradients. It is an advective flux in the sense that it drags the water. The total flux of water is then written as in Eq. (2):

$$ \mathbf{j}_l^w = -D \rho_l \nabla w_l^w + \rho_l \mathbf{q}_l \quad (2) $$

Where $D$ [m²/s] is the diffusion coefficient, which actually corresponds to the dissolved species, in this case the air, which is dissolved in the water.

a) Results from OBC case as compared with measurements

<figure><img src="image_10.png" /></figure>

b) Results from the current new model

<figure><img src="image_11.png" /><figcaption>Figure 4. Relative humidity evolution.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

a) Results from OBC case as compared with measurements

<figure><img src="image_13.png" /></figure>

b) Results from the current new model

<figure><img src="image_14.png" /><figcaption>Figure 5. Evolution of stresses</figcaption></figure>

a) Results from OBC case as compared with measurements

<figure><img src="image_15.png" /></figure>

b) Results from the current new model

<figure><img src="image_16.png" /><figcaption>Figure 6. Evolution of pressure in the rock.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 7. Distributions for 2500 days and for 20 years of: a) temperature; b) distribution of degree of saturation; c) distribution of mean stress.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 8. Evolution of temperature at sections D1, F2 and G.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 9. Evolution of pressure at sections F2 and G.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 10. Evolution of relative humidity (up) and stress (down) at different sections for the base case.</figcaption></figure>

Soler (2001) describes the variations of thermo-osmotic permeability in the following way: “Reported values of k for Na-saturated kaolinite and Na bentonite (Dirksen, 1969),
293 at an average temperature of 25 °C and various temperature gradients and porosities, 294 range between 10⁻¹⁴ and 3·10⁻¹³ m²/K/s. In another study [11] thermal osmosis across a 295 kaolinite membrane was investigated. An estimate of k from the data of this study 296 (Srivastava, 1975), corresponding to conditions similar to the ones giving the high k 297 values (Dirksen, 1969), yields a value of 2.6·10⁻¹⁰ m²/K/s, which is three orders of 298 magnitude greater. As a first approximation, the values of k mentioned above can be 299 used to define a range of values kT = 10⁻¹⁴ – 10⁻¹⁰ m²/K/s”.

300 In order to compare this value with the value corresponding to the pressure gradient (i.e. 301 the dominant effect), it is necessary to calculate the ratio between intrinsic permeability 302 and viscosity. For the intrinsic permeability value of the Febex bentonite a value of 303 1.9·10⁻²¹ m² which divided by the viscosity of water 0.001 Pa·s, leads to: 1.9·10⁻¹² 304 m²/MPa/s. These two values are comparable.

305 In order to discriminate the relative dominance of fluxes in this case, it is also necessary 306 to introduce the gradients of pressure and temperature. While the temperature gradient 307 can be estimated easily and does not change very much with time (although temperature 308 and water content change), the pressure gradient varies more, in space and time. 309 Therefore, depending on the zone in the bentonite and the time, the advective flux may 310 be influenced by thermal osmosis and this should be investigated.

311 In order to implement in an easy way the thermal gradient term, Equation (1) can be 312 written as:

$$
313 \mathbf{q}_l = - \frac{k_{r,l}}{\mu_l} (\nabla p + \alpha \nabla T + \rho_l g \nabla z) \qquad (3)
$$

314 In FEBEX, the pressure and the temperature gradients are opposite in sign. Actually the 315 thermo-osmotic flux delays hydration and therefore it may explain some delayed 316 effects. The parameter α in this equation is the ratio between the thermo osmotic 317 conductivity and the hydraulic conductivity following Equation (3). As we have seen 318 that the order of magnitude can be similar, it follows that α=1 gives a condition in 319 which the thermal osmosis induced flux is comparable to the hydraulic flux. This 320 situation occurs because the bentonite has very low hydraulic conductivity and therefore 321 it is possible to find a value of the conductivity for thermal osmosis which leads to α=1 322 (i.e. 1.9·10⁻¹² m²/K/s, the value calculated above but in the units of kT).

323 The evolution of relative humidity in the case in which the thermal gradient induced 324 flux is incorporated is shown in Figure 11 (upper left) and Figure 12 (upper left), for 325 section F2. The water flow induced by the temperature gradient delays hydration. 326 Actually, for the values of the parameters considered, it manages to compensate the 327 water inflow up to a point that for steady state temperature conditions, the relative 328 humidity stabilizes and does not tend to 100%. This means that the bentonite barrier is 329 not fully saturated and hence the isolation would be achieved.

330 The evolution of the stress development is also shown in Figure 11 (lower left) and 331 Figure 12 (lower left). Since the hydration reaches a point that does not progress, and 332 relative humidity stabilizes, the stress also stabilizes below the swelling pressure.

333 Compared to the base case (see Figure 10), we can observe that the results of the base 334 case are closer to the real data than those obtained considering this thermal gradient 335 induced flux. In fact, considering this induced flux, the relative humidity seems to stay 336 constant with time, at least at points not far from the heater, which is not the result that 337 the real measurements are given.

338
## 4.2. Sensitivity to permeability

An analysis about the sensitivity of the model to permeability changes has also been performed, since permeability is the variable subjected to bigger variations due to changes in density, structure and construction of the barrier, measurement difficulties, etc. For the base case presented in this paper, permeability of the bentonite has been increased by a factor of 2, which is a reasonable variation. This accelerates somewhat the hydration as the permeability of the buffer controls the hydration rate in FEBEX. Moreover, as liquid pressure increases faster due to permeability increase, the stress developed also increases faster leading to values which are still comparable to the measurements (Figures 11 and 12).

349

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>Figure 11. Evolution of relative humidity (up) and stress (down) at section F2 for two additional cases: thermal gradient induced flux (left) and double permeability (right).</figcaption></figure>

352
353

Actually, the stress development in the case of double permeability is quite consistent with measurements, giving closer values to the real data than the base case. It can also be observed that even doubling the permeability, full saturation of the bentonite barrier is not achieved, even after 20 years.

354

355

356

357

<figure><img src="image_25.png" /><figcaption>Thermal gradient induced flux case</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Double permeability case</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Thermal gradient induced flux case</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Double permeability case</figcaption></figure>

358

Figure 12. Evolution of relative humidity (up) and stress (down) at section H for relative humidity and G for stresses for two additional cases: thermal gradient induced flux (left) and double permeability (right).

359

360

361

362

### 4.3. Double porosity model

363

The effect on the hydration and swelling including a double structure model is considered in this section. The double porosity model is an extension of the BBM approach presented above. The double structure effects are considered in the following way: two structural

364

365
366 levels are considered with different suctions. The two suctions are related by an
367 interaction term.

368 Equation 4 shows the simplified mass balance equations incorporating the microstructure.

369 It is assumed that porosity can be decomposed into two adding contributions.

$$
\frac{\partial}{\partial t} \left( \phi_M \rho_{IM}^w S_{IM} + \phi_M \rho_{gM}^w S_{gM} \right) + \nabla \cdot \left( \mathbf{j}_{IM}^w + \mathbf{j}_{gM}^w \right) = -\Gamma^w \quad (M) \quad (4)
$$

$$
\frac{\partial}{\partial t} \left( \phi_m \rho_{lm}^w S_{lm} \right) = \Gamma^w \qquad (m)
$$

370 Where an interchange term has been introduced. This term is calculated as a function of
371 the suctions in the micro pores and the macro pores (Eq. 5).

$$
\Gamma^w = \alpha_{\text{micro-macro}} (s_m - s_M) \tag{5}
$$

372 Where $s_m$ is the suction micro and $s_M$ corresponds to the suction macro. The macro and
373 micro porosities are defined following Eq. 6.

$$
\phi = \phi_M + (1 - \phi_M) \phi_{m'} = \phi_M + \phi_M \quad (6)
$$

374 Where the porosity of the microstructure may be defined in two different ways, i.e. with
375 respect to the volume of the microstructure or with respect to the total volume. The
376 backfill material in the FEBEX is composed by blocks. The macro porosity is assumed
377 to be 0.40 while the micro porosity is set to 0.15 in these model calculations.

378 The combination of the double structure and the BBM model consists –in the model used
379 here– in using the suction in the microstructure for the mechanical calculations. For
380 instance, for the swelling terms, the volumetric elastic deformation is calculated as in Eq.
381 7.

$$
d\varepsilon_v = \frac{\kappa_p}{1+e} \frac{dp'}{p'} + \frac{\kappa_s}{1+e} \frac{ds_m}{(s_m + p_{atm})} \quad (7)
$$

382 In addition, the micropores contribute to water storage and a specific retention curve is
383 considered, which essentially imposes that the micropores remain saturated (Eq. 8).

$$
S_{lm} = \exp\left(-\frac{(s_m)^{\beta_{\text{micro}}}}{P_{o \text{ micro}}}\right) \qquad (8)
$$

384 Where $\alpha_{\text{micro-macro}} = 5 \times 10^{-12} \text{ kg/MPa/m}^3/\text{s}$, $\beta_{\text{mic}} = 1.150$ and $P_{o \text{ micro}} = 700$ MPa.

385 The results of suction development are shown in figure 13. The micro suction shows a
386 smoothed response when sharp variations of macro suction take place. This is an expected
387 response as the water flow from macro to micro and vice-versa is delayed by the
388 interaction parameter $\alpha$.
389

390

391

392

393

394

395

396

Moreover, the results of the stress development when the double porosity is considered are quite similar to the results obtained in the base case, and thus are not shown in this paper for the sake of briefness. The swelling stress in the FEBEX experiment is an integrated response of expansion of the backfill in a heterogeneous way as the process of hydration is modified by the heating. The suction in the micro and macro pores is not equal but it may happen that the averaged values are not so different. Actually, this is controlled by the capacity to hydrate and this, in turn, is controlled by the permeability of the backfill.

397

<figure><img src="image_29.png" /></figure>

398

399

400

<caption>Figure 13. Evolution of macro and micro suction. Dashed line is considered for micro suction while the continuous line corresponds to the macro suction.</caption>

401

402

## 5. Conclusions

403

404

405

406

407

408

409

410

This paper presents a model that has been prepared based on previous modeling efforts performed for FEBEX. The model produces similar results as the previous models taking into account that some slight modifications have been introduced on the constitutive equations (for instance, the retention curve). Once calibrated, the long term response of the model predictions up to 20 years has also been performed, taking in consideration the operations of shut down and dismantling of the first canister, concrete abutment and bentonite. The evolution of temperature, degree of saturation, stress, water pressure, relative humidity have been analyzed.

411

412

413

414

415

416

417

The model has been used to analyze the effect of various processes, such as the incorporation of thermal osmosis, the variation of intrinsic permeability and the introduction of a double porosity model. The incorporation of thermal-osmosis leads to the possibility of developing a steady state situation in which the water flow induced by the thermal gradient is compensated by the water flow induced by the hydraulic gradient. This implies that the clay will not reach a full saturation even for long term conditions.

418

419

420

421

422

The vapor flux, although it can also be related to a thermal gradient is no able to reach this condition as it vanishes when the gas degree of saturation reduces. As the material progressively saturates, the vapor flux reduces and finally is zero at full saturation. In contrast, the liquid water flow induced by a thermal gradient increases with water saturation and reaches its maximum for a given thermal gradient, at full saturation.
The intrinsic permeability has been shown to be a parameter that permits a significant improvement simply by varying it by a factor of two. In fact, the results obtain for the double permeability case are closer to real data than the base case, meanwhile the thermal gradient induced flux case gives worse results than the base case, if we compare them with the real data.

Finally, the double structure model has given different results for micro and macro suction but we found that does not have great influence in the general behavior of the model, as the stress development obtained is quite similar to the base case.

432

433

## 434 References

Alonso E, Gens A, Josa A. A constitutive model for partially saturated soils. Geotech. 1990;40(3):405-30.

Alonso E, Alcoverro J. The FEBEX test as a benchmark case for THM modelling: historical perspective and lessons learnt. In: Alonso E, Ledesma A, editors. Advances in understanding engineered clay barriers, London: Taylor & Francis; 2005, p. 3-19.

Dirksen C. Thermo-osmosis through compacted saturated clay membranes. Soil Sci. Soc. Am. Proc. 1969;33:821-6.

Dixon D, Chandler N, Graham J, Gray MN. Two large-scale sealing tests conducted at Atomic Energy of Canada's underground research laboratory: the buffer-container experiment and the isothermal test. Can Geotech J. 2002;39(3):503-18.

Gens A, Garcia-Molina AJ, Olivella S, Alonso E, Huertas F. Analysis of a full scale in situ test simulating repository conditions. Int J Numer Anal Methods Geomech. 1998;22(7):515-48.

Gens A. The role of geotechnical engineering in nuclear energy utilisation: special lecture. In: Proc 13th Eur Conf Soil Mech Geotech Eng, Prague; 2003;3:25-67.

Gens A, Sanchez M, Guimaraes LDN, Alonso E, Lloret A, Olivella S et al. A full-scale in situ heating test for high-level nuclear waste disposal: observations, analysis and interpretation. Geotech. 2009;59(4):377-99. doi: 10.1680/geot.2009.59.4.377

Olivella S, Gens A, Carrera J, Alonso E. Numerical formulation for a simulator 'CODE_BRIGHT' for the coupled analysis of saline media. Engng Comput. 1996;13(7):87-112.

Pusch R, Borgesson L, Ramqvist G. Final report of the mass buffer test -volume II: Test results, Stripa Project 85/12. Stockholm: SKB; 1985.

Sanchez M, Gens A, Guimaraes L. Thermal-hydraulic-mechanical (THM) behaviour of a large-scale in situ heating experiment during cooling and dismantling. Canadian Geotechnical Journal, 2012, 49(10): 1169-1195, https://doi.org/10.1139/t2012-076

Selvadurai APS. Hydro-thermo-mechanics of engineered clay barriers and geological barriers. Eng Geol. 1997;47(4):311-2. doi: 10.1016/S0013-7952(97)00036-7

Soler JM. The effect of coupled transport phenomena in the Opalinus Clay and implications for radionuclide transport. J Contaminant Hydrology 2001;53(1-2):63-84.
465 Srivastava RC, Avasthi PK. Non-equilibrium thermodynamics of thermo-osmosis of
466 water through kaolinite. J. Hydrol. 1975;24:111-20.

467 Volckaert G, Bernier F, Alonso E, Gens A, Samper J, Villar MV, Martin PL, Cuevas J,
468 Campos R, Thomas HR, Imbert C, Zingarelli V. Thermal-hydraulic-mechanical and
469 geochemical behaviour of the clay barrier in radioactive waste repositories (model
470 development and validation), Nuclear Science and Technology, EUR 16744.

471 Luxembourg: Commission of the European Communities; 1996.

472