# Gas-entry pressure impact on the evaluation of hydrogen migration at different scales of a deep geological disposal of radioactive waste

Zakaria Saâdi

► To cite this version:

Zakaria Saâdi. Gas-entry pressure impact on the evaluation of hydrogen migration at different scales of a deep geological disposal of radioactive waste. Scientific Reports, 2024, 14 (1), pp.6221. 10.1038/s41598-024-56454-y. irsn-04582711

HAL Id: irsn-04582711

https://irsn.hal.science/irsn-04582711

Submitted on 22 May 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Gas-entry pressure impact on the evaluation of hydrogen migration at different scales of a deep geological disposal of radioactive waste

Zakaria Saâdi

► To cite this version:

Zakaria Saâdi. Gas-entry pressure impact on the evaluation of hydrogen migration at different scales of a deep geological disposal of radioactive waste. Scientific Reports, 2024, 14 (1), pp.6221. 10.1038/s41598-024-56454-y. irsn-04582711

HAL Id: irsn-04582711

https://asnr.hal.science/irsn-04582711v1

Submitted on 22 May 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Gas-entry pressure impact on the evaluation of hydrogen migration at different scales of a deep geological disposal of radioactive waste

Zakaria Saâdi

Although the importance of gas-entry pressure in simulating two-phase liquid-gas flow in porous media has been studied at the column and borehole scales, its impact on the simulation of transient hydraulic-gas at different scales of a deep geological repository of radioactive waste (DGR) in low permeability clay rock during the post-closure phase has not yet been studied. The purpose of this work is to show that neglecting this phenomenon can lead to underestimation of the maximum gas pressure and water-gas fluxes simulated within the host rock and backfilled drift network. This could impact the performance of the engineered barrier system of a DGR. Simulations performed for a high-level waste disposal cell and for a simplified repository composed of hundreds of disposal cells situated in a clay host rock, show that gas preferentially migrates through the DGR components with low capillary entry pressures, such as the excavation damaged zone (Refers to the zone where fractures develop due to failure of the rock mass around galleries after tunneling) (EDZ), the engineered barriers materials (backfill, bentonite-plug...) and interfaces between the EDZ and these materials. Such a result could have significant consequences on the performance of a repository, due to the accumulation of gas in the drift network and high increase of gas pressure, which could lead to the host rock hydraulic fracturing.

Understanding gas transport processes is one of the key issues in the assessment of radioactive waste repository performance and is the focus of many research programs (e.g., Fate Of Repository GasEs: FORGE, 2009–2013; European Joint Programme on Radioactive Waste: EURAD, 2019–2024). The purpose of these research programs is (i) to develop and calibrate models on experimental set-ups under controlled flow-rate and/or gas pressure injection, and to carry out a series of experiments tests on clay rock samples to study the conditions under which gas breakthrough processes occur, and to analyze rock discontinuities, hydro-mechanical variables, induced desaturation, inflow and outflow gas, and preferential paths created, and (ii) to study upscaling of hydrodynamic properties of the clay material derived from these experiments to develop large spatial scale numerical models of gas migration in the repository.

Many in-situ (boreholes) and laboratory gas breakthrough experiments showed the importance of gas-entry pressure (P<sub>g,e</sub>) in two-phase flow through clay-rich materials (e.g.<sup>1-6</sup>). In laboratory experiments, gas injection pressure is usually applied on one side of a clay sample, where gas pressure is increased until the gas entry pressure is exceeded, and gas ingress in the saturated sample is observed. Once the gas migrates through the clay sample, it changes pore-water pressure, swelling pressure and flowrate observed by the monitoring instruments placed around the sample. When gas breakthrough occurs, gas outflow is recorded from the changes in pressure of the backpressure. Since the clay structure at the point of injection is minimally disturbed, lower injection pressures can be used and hence the risk for unwanted fracturing and reflux is decimated. Due to the lower injection pressures, large injection volumes can be injected what results in a larger radius of influence of the injection.

Constant gas injection rate experiments are usually carried out to induce high gas pressures at the boundary of the clay material to estimate specifically its gas breakthrough pressure. However, the change in gas pressure within the sample is usually abrupt and can lead to the creation of fractures especially at the edges of the clay

Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSE-ENV, SPDR, UEMIS, 92262 Fontenay-aux-Roses, France. email: zakaria.saadi@irsn.fr
<figure><img src="image_2.png" /></figure>

OPEN

# Gas-entry pressure impact on the evaluation of hydrogen migration at different scales of a deep geological disposal of radioactive waste

Zakaria Saâdi

Although the importance of gas-entry pressure in simulating two-phase liquid-gas flow in porous media has been studied at the column and borehole scales, its impact on the simulation of transient hydraulic-gas at different scales of a deep geological repository of radioactive waste (DGR) in low permeability clay rock during the post-closure phase has not yet been studied. The purpose of this work is to show that neglecting this phenomenon can lead to underestimation of the maximum gas pressure and water-gas fluxes simulated within the host rock and backfilled drift network. This could impact the performance of the engineered barrier system of a DGR. Simulations performed for a high-level waste disposal cell and for a simplified repository composed of hundreds of disposal cells situated in a clay host rock, show that gas preferentially migrates through the DGR components with low capillary entry pressures, such as the excavation damaged zone (Refers to the zone where fractures develop due to failure of the rock mass around galleries after tunneling) (EDZ), the engineered barriers materials (backfill, bentonite-plug...) and interfaces between the EDZ and these materials. Such a result could have significant consequences on the performance of a repository, due to the accumulation of gas in the drift network and high increase of gas pressure, which could lead to the host rock hydraulic fracturing.

Understanding gas transport processes is one of the key issues in the assessment of radioactive waste repository performance and is the focus of many research programs (e.g., Fate Of Repository GasEs: FORGE, 2009–2013; European Joint Programme on Radioactive Waste: EURAD, 2019–2024). The purpose of these research programs is (i) to develop and calibrate models on experimental set-ups under controlled flow-rate and/or gas pressure injection, and to carry out a series of experiments tests on clay rock samples to study the conditions under which gas breakthrough processes occur, and to analyze rock discontinuities, hydro-mechanical variables, induced desaturation, inflow and outflow gas, and preferential paths created, and (ii) to study upscaling of hydrodynamic properties of the clay material derived from these experiments to develop large spatial scale numerical models of gas migration in the repository.

Many in-situ (boreholes) and laboratory gas breakthrough experiments showed the importance of gas-entry pressure (P<sub>g,e</sub>) in two-phase flow through clay-rich materials (e.g.<sup>1-6</sup>). In laboratory experiments, gas injection pressure is usually applied on one side of a clay sample, where gas pressure is increased until the gas entry pressure is exceeded, and gas ingress in the saturated sample is observed. Once the gas migrates through the clay sample, it changes pore-water pressure, swelling pressure and flowrate observed by the monitoring instruments placed around the sample. When gas breakthrough occurs, gas outflow is recorded from the changes in pressure of the backpressure. Since the clay structure at the point of injection is minimally disturbed, lower injection pressures can be used and hence the risk for unwanted fracturing and reflux is decimated. Due to the lower injection pressures, large injection volumes can be injected what results in a larger radius of influence of the injection.

Constant gas injection rate experiments are usually carried out to induce high gas pressures at the boundary of the clay material to estimate specifically its gas breakthrough pressure. However, the change in gas pressure within the sample is usually abrupt and can lead to the creation of fractures especially at the edges of the clay

Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSE-ENV, SPDR, UEMIS, 92262 Fontenay-aux-Roses, France. email: zakaria.saadi@irsn.fr
sample (interfaces between the clay and the cell of confinement), and thus to a preferential gas transport through these pathways (fractures). Coupled hydro-gas-mechanical 3D numerical models are difficult to apply to such experiments, that why it is preferred to inject small gas fluxes during longer times to slow down the kinetics of gas migration in the sample and to avoid earlier such mechanical damages in the clay structure.

For example, for the case of the Callovo-Oxfordian (COx) claystone, the rock proposed for hosting the French Deep Geological Repository (DGR) project Cigéo⁷, the measured capillary entry pressure ($P_{c,e} = P_{g,e} - P_{w}$) ranges between 0,3 and 3 MPa for damaged (fissured) COx samples, and can be greater than 5 MPa for intact COx samples⁴.

Two in-situ experiments of gas injection in boreholes inside the Opalines clay have been carried out by Marshall et al.² to characterize gas entry pressure and gas permeability of the Opalines clay. These experiments have been implemented with two events as constant flow rate ("gas threshold pressure test") and as multi-step constant gas pressure. Simulation of these two-phase flow experiments was based on iTOUGH2 code⁸ with the classical VGM model for hydrodynamic properties of the Opalines clay. However, their analysis was not relevant because neither capillary gas entry pressure $P_{c,e}$ nor hysteresis were modelled.

A recent paper⁹ showed that $N_2$ gas pressure inside a chamber, where $N_2$ pulses are injected through a borehole inside the COx (PGZ1 experiment,¹⁰⁻¹²), can be highly underestimated when $P_{c,e}$ is neglected in the COx hydraulic properties, especially during long time periods when capillary pressure in the COx becomes higher than its capillary gas entry pressure ($P_c = P_w - P_g \ge P_{c,e}$). During these periods, the $N_2$ gas accumulates more in the chamber, and its pressure increases more significantly due to the high capillary gas entry pressure of the COx. Indeed, the COx behaves as a capillary barrier for gas entry from the chamber deeply in the COx leading to the delay of the desaturation of EDZ and COx, but in a more significative way and in a smallest extend in space during time (by comparison to a simulation by parametrization $P_{c,e} = 0$). The desaturation of the COx is limited to a thin layer near the EDZ/COx interface, and the kinetic of the diffusion of dissolved $N_2$ in the water saturated COx is reduced in time.

In the last past two decades, although many studies dealt with the numerical modelling of gas migration through the different components of a DGR (e.g.¹³⁻¹⁵), none of them studied the impact of a non-zero $P_{c,e}$ in the hydraulic properties of the porous materials representing the host rock, the EDZ, and the engineered barriers materials. The main reason for neglecting $P_{c,e}$ is related to the difficulty of the implementation and verification of the thermodynamic conditions necessary for the switch from single-liquid water to two-phase liquid-gas conditions and vice-versa. Materials in a DGR typically show strong differences in capillary entry pressures. Consideration of $P_{c,e}$ makes it more difficult to achieve a stable and convergent numerical solution of a mathematical problem that involves materials with very contrasting properties.

The purpose of this work is to review the models of gas migration in a DGR at the scale of a HLW disposal cell and of a small repository ("module") of hundreds of HLW cells in order to assess uncertainties on gas pressure and two-phase flow due to non-negligible $P_{c,e}$. Therefore, we re-simulate the benchmarks models that were proposed in the Euratom 7th Framework Program project FORGE (Fate Of Repository GasEs¹⁴,¹⁶,¹⁷) using the modified Van Genuchten-Mualem capillary model VGMPE proposed by⁹ with $P_{c,e} \neq 0$ for all materials composing the HLW cells and the drift network. This capillary model has been implemented in new versions of iTOUGH2/ EOS5⁸ and TOUGH2-MP¹⁸ codes by accounting for the thermodynamic changes from single to two-phase conditions with $P_{c,e} \neq 0$.

# The studied benchmark models

## The cell scale

The benchmark model represents an axisymmetric HLW-cell embedded in the COx clay host rock, as shown in Fig. 1a. The interfaces Canisters-EDZ and Plug (bentonite)-COx, considered as a centimeter-thick region, and EDZ around the borehole for canisters emplacement, are also represented. A hydrogen source term due to anoxic corrosion, corresponding to a gas flux rate of 100 mol/year/cell lasting 10 000 years is assumed. It is localized on the waste canister external surfaces (purple continuous line in Fig. 1a).

Gravity effect is assumed negligible compared to capillary forces, because of the very low permeability of all the clay-rich materials composing the cell. Time-dependent boundary conditions are considered in the backfilled access drift (green discontinuous line in Fig. 1a) by assuming time-varying gas pressure and gas saturation, as shown in Fig. 1b¹⁶. The simulation time is limited to 100 000 years.

In this work, the classical VGM-model describing the hydraulic properties of the different materials composing the cell has been changed to introduce a non-zero capillary gas entry pressure $P_{c,e}$ (Pa). The new VGMPE model proposed by⁹ has been used. It is a generalization of¹⁹-model for two-phase flow with a non-zero capillary entry pressure $P_{c,e}$ in the constitutive VGM²⁰⁻²² relationships for the water retention curve $P_c(S_l)$ and the relative permeability to liquid and gas curves $k_{r,\beta=l,g}(S_l)$. This model is described by the following equations:

$$P_c(S_l) = \begin{cases} -\frac{1}{\alpha} \left[ (S_e^* S_e)^{-\frac{1}{m}} - 1 \right]^{\frac{1}{n}} & \text{if } S_e \le 1 \\ 0 & \text{elsewhere} \end{cases} \quad (1)$$

$$k_{r,l}(S_l) = \frac{k_l(S_l)}{k_{0,l}} = \begin{cases} [S_e]^\tau \left[ \frac{1 - (1 - (S_e^* S_e)^{\frac{1}{m}})^m}{1 - (1 - (S_e^* S_e)^{\frac{1}{m}})^m} \right]^2 & \text{if } S_e < 1 \\ 1 & \text{if } S_e = 1 \end{cases} \quad (2)$$
sample (interfaces between the clay and the cell of confinement), and thus to a preferential gas transport through these pathways (fractures). Coupled hydro-gas-mechanical 3D numerical models are difficult to apply to such experiments, that why it is preferred to inject small gas fluxes during longer times to slow down the kinetics of gas migration in the sample and to avoid earlier such mechanical damages in the clay structure.

For example, for the case of the Callovo-Oxfordian (COx) claystone, the rock proposed for hosting the French Deep Geological Repository (DGR) project Cigéo⁷, the measured capillary entry pressure ($P_{c,e} = P_{g,e} - P_{w}$) ranges between 0,3 and 3 MPa for damaged (fissured) COx samples, and can be greater than 5 MPa for intact COx samples⁴.

Two in-situ experiments of gas injection in boreholes inside the Opalines clay have been carried out by Marshall et al.² to characterize gas entry pressure and gas permeability of the Opalines clay. These experiments have been implemented with two events as constant flow rate ("gas threshold pressure test") and as multi-step constant gas pressure. Simulation of these two-phase flow experiments was based on iTOUGH2 code⁸ with the classical VGM model for hydrodynamic properties of the Opalines clay. However, their analysis was not relevant because neither capillary gas entry pressure $P_{c,e}$ nor hysteresis were modelled.

A recent paper⁹ showed that $N_2$ gas pressure inside a chamber, where $N_2$ pulses are injected through a borehole inside the COx (PGZ1 experiment,¹⁰⁻¹²), can be highly underestimated when $P_{c,e}$ is neglected in the COx hydraulic properties, especially during long time periods when capillary pressure in the COx becomes higher than its capillary gas entry pressure ($P_c = P_w - P_g \ge P_{c,e}$). During these periods, the $N_2$ gas accumulates more in the chamber, and its pressure increases more significantly due to the high capillary gas entry pressure of the COx. Indeed, the COx behaves as a capillary barrier for gas entry from the chamber deeply in the COx leading to the delay of the desaturation of EDZ and COx, but in a more significative way and in a smallest extend in space during time (by comparison to a simulation by parametrization $P_{c,e} = 0$). The desaturation of the COx is limited to a thin layer near the EDZ/COx interface, and the kinetic of the diffusion of dissolved $N_2$ in the water saturated COx is reduced in time.

In the last past two decades, although many studies dealt with the numerical modelling of gas migration through the different components of a DGR (e.g.¹³⁻¹⁵), none of them studied the impact of a non-zero $P_{c,e}$ in the hydraulic properties of the porous materials representing the host rock, the EDZ, and the engineered barriers materials. The main reason for neglecting $P_{c,e}$ is related to the difficulty of the implementation and verification of the thermodynamic conditions necessary for the switch from single-liquid water to two-phase liquid-gas conditions and vice-versa. Materials in a DGR typically show strong differences in capillary entry pressures. Consideration of $P_{c,e}$ makes it more difficult to achieve a stable and convergent numerical solution of a mathematical problem that involves materials with very contrasting properties.

The purpose of this work is to review the models of gas migration in a DGR at the scale of a HLW disposal cell and of a small repository ("module") of hundreds of HLW cells in order to assess uncertainties on gas pressure and two-phase flow due to non-negligible $P_{c,e}$. Therefore, we re-simulate the benchmarks models that were proposed in the Euratom 7th Framework Program project FORGE (Fate Of Repository GasEs¹⁴,¹⁶,¹⁷) using the modified Van Genuchten-Mualem capillary model VGMPE proposed by⁹ with $P_{c,e} \neq 0$ for all materials composing the HLW cells and the drift network. This capillary model has been implemented in new versions of iTOUGH2/ EOS5⁸ and TOUGH2-MP¹⁸ codes by accounting for the thermodynamic changes from single to two-phase conditions with $P_{c,e} \neq 0$.

## The studied benchmark models

### The cell scale

The benchmark model represents an axisymmetric HLW-cell embedded in the COx clay host rock, as shown in Fig. 1a. The interfaces Canisters-EDZ and Plug (bentonite)-COx, considered as a centimeter-thick region, and EDZ around the borehole for canisters emplacement, are also represented. A hydrogen source term due to anoxic corrosion, corresponding to a gas flux rate of 100 mol/year/cell lasting 10 000 years is assumed. It is localized on the waste canister external surfaces (purple continuous line in Fig. 1a).

Gravity effect is assumed negligible compared to capillary forces, because of the very low permeability of all the clay-rich materials composing the cell. Time-dependent boundary conditions are considered in the backfilled access drift (green discontinuous line in Fig. 1a) by assuming time-varying gas pressure and gas saturation, as shown in Fig. 1b¹⁶. The simulation time is limited to 100 000 years.

In this work, the classical VGM-model describing the hydraulic properties of the different materials composing the cell has been changed to introduce a non-zero capillary gas entry pressure $P_{c,e}$ (Pa). The new VGMPE model proposed by⁹ has been used. It is a generalization of¹⁹-model for two-phase flow with a non-zero capillary entry pressure $P_{c,e}$ in the constitutive VGM²⁰⁻²² relationships for the water retention curve $P_c(S_l)$ and the relative permeability to liquid and gas curves $k_{r,\beta=l,g}(S_l)$. This model is described by the following equations:

$$P_c(S_l) = \begin{cases} -\frac{1}{\alpha} \left[ (S_e^* S_e)^{-\frac{1}{m}} - 1 \right]^{\frac{1}{n}} & \text{if } S_e \le 1 \\ 0 & \text{elsewhere} \end{cases} \quad (1)$$

$$k_{r,l}(S_l) = \frac{k_l(S_l)}{k_{0,l}} = \begin{cases} [S_e]^\tau \left[ \frac{1 - (1 - (S_e^* S_e)^{\frac{1}{m}})^m}{1 - (1 - (S_e^* S_e)^{\frac{1}{m}})^m} \right]^2 & \text{if } S_e < 1 \\ 1 & \text{if } S_e = 1 \end{cases} \quad (2)$$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Schematic representation of the radial cross section dimensions of the cylindrical waste cell model (assuming plane symmetry along its axial axis), as well as the boundary conditions applied at the surface of the domain. Dark points are those where results should be given (adapted from Wendling et al. 2013a). (b) Time-dependent gas saturation and gas pressure applied at the backfilled access drift.</figcaption></figure>

$$k_{r,g}(S_l) = \frac{k_g(S_l)}{k_{0,g}} = \begin{cases} [1 - S_e]^{\tau'} \left[ 1 - \frac{1 - (1 - (S_e^* S_e)^{\frac{1}{m}})^m}{1 - (1 - (S_e^* S_e)^{\frac{1}{m}})^m} \right]^2 & \text{if } S_e < 1 \\ 0 & \text{if } S_e = 1 \end{cases} \quad (3)$$

With:

$$S_e^* = \frac{S_{ls} - S_{lr}}{S_{ls}^* - S_{lr}} = [1 + (\alpha P_{c,e})^{n'}]^{-m} \quad S_e = \frac{S_l - S_{lr}}{S_{ls} - S_{lr}} \quad (4)$$

In the equations above:

m and n: are dimensionless shape parameters related by the relation $m = 1 - \frac{1}{n}$.

α = 1/P*: in (Pa⁻¹), with P_r is a scaling pressure parameter (Pa) (α > 0; if P_r is positive).

τ and τ': are dimensionless parameters representing tortuosity for relative permeability to liquid and gas respectively. In this study τ = τ' = 0.5.

S_{lr} and S_{ls}: are residual and maximum (or full) liquid saturation, respectively (-). In this study, S_{lr} = 0 and S_{ls} = 1.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Schematic representation of the radial cross section dimensions of the cylindrical waste cell model (assuming plane symmetry along its axial axis), as well as the boundary conditions applied at the surface of the domain. Dark points are those where results should be given (adapted from Wendling et al. 2013a). (b) Time-dependent gas saturation and gas pressure applied at the backfilled access drift.</figcaption></figure>

$$k_{r,g}(S_l) = \frac{k_g(S_l)}{k_{0,g}} = \begin{cases} [1 - S_e]^{\tau'} \left[ 1 - \frac{1 - (1 - (S_e^* S_e)^{\frac{1}{m}})^m}{1 - (1 - (S_e^* S_e)^{\frac{1}{m}})^m} \right]^2 & \text{if } S_e < 1 \\ 0 & \text{if } S_e = 1 \end{cases} \quad (3)$$

With:

$$S_e^* = \frac{S_{ls} - S_{lr}}{S_{ls}^* - S_{lr}} = [1 + (\alpha P_{c,e})^{n'}]^{-m} \quad S_e = \frac{S_l - S_{lr}}{S_{ls} - S_{lr}} \quad (4)$$

In the equations above:

m and n: are dimensionless shape parameters related by the relation $m = 1 - \frac{1}{n}$.

α = 1/P: in (Pa⁻¹), with P<sub>r</sub> is a scaling pressure parameter (Pa) (α > 0; if P<sub>r</sub> is positive).

τ and τ': are dimensionless parameters representing tortuosity for relative permeability to liquid and gas respectively. In this study τ = τ' = 0.5.

S<sub>lr</sub> and S<sub>ls</sub>: are residual and maximum (or full) liquid saturation, respectively (−). In this study, S<sub>lr</sub> = 0 and S<sub>ls</sub> = 1.
$S_e$: is the effective saturation of liquid phase (-).

$S_e^*$: is the new effective liquid saturation, which reaches its maximal value when $S_I = S_{Is} = 1$ (-).

$k_{0,l}$ and $k_{0,g}$: are scaling permeability parameters, defined as intrinsic permeabilities to liquid and gas, respectively (m²).

In Eq. (1) the water retention curve is extended beyond the full liquid saturation ($S_{Is} = 1$) by introducing the parameter $S_{Is}^* \ge 1$ in the capillary VGM model. As reported in⁹, relative permeability to gas and liquid is more sensitive to $P_{c,e}$ than the water retention curve near full liquid saturation. Notice that when $P_{c,e} = 0$, Eq. (4) implies that $S_e^* = 1$ and $S_{Is}^* = S_{Is}$, then Eqs. (1), (2) and (3) are reduced to the classical VGM-model equations proposed in²².

The physical properties (i.e., porosity, intrinsic permeability, pore compressibility...) of the different materials composing the disposal cell, and parameters of their hydrodynamic properties in Eqs. (1), (2) and (3), are given in Table 1¹⁶. The $P_{c,e}$-value 2 MPa for the COx has been optimized through best fit of $P_c(S_l)$, $k_{r,l}(P_c)$, $k_{r,g}(P_c)$ data and has been used for the PGZ1-experiment model⁹. The $P_{c,e}$-values of bentonite and EDZ are deduced in % of COx $P_{c,e}$-value optimized (2 MPa) by analogy to those measured on average for the three materials. This % is evaluated independently, from the ratio of measured $P_{c,e}$-values of bentonite and EDZ (~ 4 and 2 MPa, respectively,²³,²⁴) to measured $P_{c,e}$-value of the COx (~ 6 MPa,²³,²⁴). Therefore, $P_{c,e}$-values of bentonite and EDZ are equal to 4/6 (67%) and 2/6 (33%) times 2 MPa, respectively (Table 1). The $P_{c,e}$-value 0.1 MPa for the backfill is an approximation of the value 0.09875 MPa which has been evaluated through re-fitting simultaneously $P_c(S_l)$, $k_{r,l}(P_c)$, and $k_{r,g}(P_c)$ numerical data calculated by Eqs. (1), (2) and (3) with $P_{c,e} = 0$. This fitting, using the optimization algorithms developed by⁹, did not highly modify values of the shape parameter 'm' and the scale parameter 'P_r' in Eqs. (1), (2) and (3), given in Table 1. Like the EDZ material, the backfill has a porous structure composed of clay blocks with high capillarity (high $P_{c,e}$), and big inter-block voids without capillarity ($P_{c,e} = 0$). Although $P_r$-value of the backfill is higher than that of the EDZ, its $P_{c,e}$ value is smaller than that of the EDZ. This is in accordance with its higher permeability and porosity compared to those of the EDZ (Table 1).

As shown in Table 1, interface facing canister is constituted by voids (porosity of 100%, like a fracture with a thin layer between two parallel plates, whose permeability is calculated by Poiseuille law for the fluid flow). The interface facing bentonite plug is a mixture of big voids with the porous materials bentonite and EDZ (COx matrix and fractures), but it is largely made up of voids with important roughness. This explains its higher permeability, which is 100 times greater than that of bentonite and close to that of the EDZ. Interface facing backfill differs from that facing the canister by a different mixture of backfill, EDZ, and big voids. Its permeability is higher than the interface facing the plug because of the high permeability of the backfill compared to that of bentonite in the mixture, and to a small roughness. The three interfaces are largely made up of voids and flow by capillarity in their porous structures is negligible. Their scale parameter $P_r$ in VGMPE model is chosen very small (Table 1), and therefore have negligible $P_{c,e}$ ($P_{c,e} \approx 0$).

Transport properties are such that: the diffusion coefficients of water vapor and hydrogen in unsaturated porous materials are deduced from their diffusion coefficients in free gas and liquid through the proportionality constant (tortuosity factor) calculated by the Millington and Quirk formula²⁵. Diffusion coefficients of hydrogen in free liquid and gas phases are equal to $6.0 \times 10^{-9}$ and $9.5 \times 10^{-5} \text{ m}^2/\text{s}^{17}$, respectively, and the diffusion coefficient of water vapor in free gas phase is equal to $2.1 \times 10^{-5} \text{ m}^2/\text{s}^{15}$. The inverse Henry's constant of hydrogen is equal to $1.379 \times 10^{-10} \text{ Pa}^{-1}$.

The initial gas saturation is 30% for bentonite seals (cell and main drift) and for backfill (drifts). An initial gas saturation of 95% for all interfaces is considered. For these materials initially unsaturated, the gas pressure is equal to 1 atmosphere (~ 1 bar = 0.1 MPa), and the water pressure is deduced from the gas pressure and water saturation. Therefore, liquid pressure in capillary equilibrium with gas pressure is deduced from gas pressure and liquid saturation by applying the VGM-model of the water retention curve defined by Eq. (1). The clay host rock, and cell and drift EDZs are initially at full liquid saturation with an initial water pressure set to 5 MPa.

The mesh of the domain used for the simulations is like that recently presented in¹⁵. It is radial and rectilinear; it consists of $39 \times 122 = 4758$ elements in the plane (r, z). This number is reduced to 4358 active elements because elements representing the waste canister are not considered in the calculation. Two extra-cells with infinite volumes were added to handle boundary conditions (BC). The first one is for fixing the time-dependent gas pressure

<table><caption>**Table 1.** Physical properties and parameters of VGMPE hydrodynamic properties of the materials used in the simulations (¹⁶,¹⁷). aRadial intrinsic permeability. Calculated from horizontal and vertical intrinsic permeability values using geometric mean: ($K_x \times K_y \times K_z$)¹/³.</caption><thead><tr><th rowspan="2">Material</th><th rowspan="2">Porosity (-)</th><th colspan="2">Intrinsic permeability (m²)</th><th rowspan="2">Pore compressibility (Pa⁻¹)</th><th rowspan="2">Capillary entry pressure, P<sub>c,e</sub> (Pa)</th><th rowspan="2">VG pressure scale parameter, P<sub>r</sub> (Pa)</th><th rowspan="2">VG shape parameter, n (-)</th></tr><tr><th>Waste cell scale<sup>a</sup></th><th>Module scale</th></tr></thead><tbody><tr><td>Geological Medium (clay host rock)</td><td>0.15</td><td>7.94 × 10<sup>-21</sup></td><td>1.0 × 10<sup>-20</sup></td><td>1.11 × 10<sup>-9</sup></td><td>2.00 × 10<sup>6</sup></td><td>1.5 × 10<sup>7</sup></td><td>1.5</td></tr><tr><td>EDZ</td><td>0.15</td><td>7.94 × 10<sup>-18</sup></td><td>5.0 × 10<sup>-18</sup></td><td>1.11 × 10<sup>-9</sup></td><td>6.67 × 10<sup>5</sup></td><td>1.5 × 10<sup>6</sup></td><td>1.5</td></tr><tr><td>Bentonite</td><td>0.35</td><td>1.0 × 10<sup>-20</sup></td><td>1.0 × 10<sup>-20</sup></td><td>8.25 × 10<sup>-10</sup></td><td>1.33 × 10<sup>6</sup></td><td>1.6 × 10<sup>7</sup></td><td>1.6</td></tr><tr><td>Backfill</td><td>0.40</td><td>5.0 × 10<sup>-17</sup></td><td>5.0 × 10<sup>-17</sup></td><td>2.09 × 10<sup>-9</sup></td><td>1.00 × 10<sup>5</sup></td><td>2.0 × 10<sup>6</sup></td><td>1.5</td></tr><tr><td>Interface facing plug</td><td>0.30</td><td>7.94 × 10<sup>-18</sup></td><td>5.0 × 10<sup>-18</sup></td><td>1.11 × 10<sup>-9</sup></td><td>0</td><td>1.0 × 10<sup>4</sup></td><td>4</td></tr><tr><td>Interface facing canister</td><td>1.00</td><td>1.0 × 10<sup>-12</sup></td><td>1.0 × 10<sup>-12</sup></td><td>1.09 × 10<sup>-11</sup></td><td>0</td><td>1.0 × 10<sup>4</sup></td><td>4</td></tr><tr><td>Interface facing backfill</td><td>0.40</td><td>1.0 × 10<sup>-15</sup></td><td>1.0 × 10<sup>-15</sup></td><td>2.09 × 10<sup>-9</sup></td><td>0</td><td>1.0 × 10<sup>4</sup></td><td>4</td></tr></tbody></table>
$S_e$: is the effective saturation of liquid phase (-).

$S_e^*$: is the new effective liquid saturation, which reaches its maximal value when $S_I = S_{Is} = 1$ (-).

$k_{0,l}$ and $k_{0,g}$: are scaling permeability parameters, defined as intrinsic permeabilities to liquid and gas, respectively (m²).

In Eq. (1) the water retention curve is extended beyond the full liquid saturation ($S_{Is} = 1$) by introducing the parameter $S_{Is}^* \ge 1$ in the capillary VGM model. As reported in⁹, relative permeability to gas and liquid is more sensitive to $P_{c,e}$ than the water retention curve near full liquid saturation. Notice that when $P_{c,e} = 0$, Eq. (4) implies that $S_e^* = 1$ and $S_{Is}^* = S_{Is}$, then Eqs. (1), (2) and (3) are reduced to the classical VGM-model equations proposed in²².

The physical properties (i.e., porosity, intrinsic permeability, pore compressibility...) of the different materials composing the disposal cell, and parameters of their hydrodynamic properties in Eqs. (1), (2) and (3), are given in Table 1¹⁶. The $P_{c,e}$-value 2 MPa for the COx has been optimized through best fit of $P_c(S_l)$, $k_{r,l}(P_c)$, $k_{r,g}(P_c)$ data and has been used for the PGZ1-experiment model⁹. The $P_{c,e}$-values of bentonite and EDZ are deduced in % of COx $P_{c,e}$-value optimized (2 MPa) by analogy to those measured on average for the three materials. This % is evaluated independently, from the ratio of measured $P_{c,e}$-values of bentonite and EDZ (~ 4 and 2 MPa, respectively,²³,²⁴) to measured $P_{c,e}$-value of the COx (~ 6 MPa,²³,²⁴). Therefore, $P_{c,e}$-values of bentonite and EDZ are equal to 4/6 (67%) and 2/6 (33%) times 2 MPa, respectively (Table 1). The $P_{c,e}$-value 0.1 MPa for the backfill is an approximation of the value 0.09875 MPa which has been evaluated through re-fitting simultaneously $P_c(S_l)$, $k_{r,l}(P_c)$, and $k_{r,g}(P_c)$ numerical data calculated by Eqs. (1), (2) and (3) with $P_{c,e} = 0$. This fitting, using the optimization algorithms developed by⁹, did not highly modify values of the shape parameter 'm' and the scale parameter 'P_r' in Eqs. (1), (2) and (3), given in Table 1. Like the EDZ material, the backfill has a porous structure composed of clay blocks with high capillarity (high $P_{c,e}$), and big inter-block voids without capillarity ($P_{c,e} = 0$). Although $P_r$-value of the backfill is higher than that of the EDZ, its $P_{c,e}$ value is smaller than that of the EDZ. This is in accordance with its higher permeability and porosity compared to those of the EDZ (Table 1).

As shown in Table 1, interface facing canister is constituted by voids (porosity of 100%, like a fracture with a thin layer between two parallel plates, whose permeability is calculated by Poiseuille law for the fluid flow). The interface facing bentonite plug is a mixture of big voids with the porous materials bentonite and EDZ (COx matrix and fractures), but it is largely made up of voids with important roughness. This explains its higher permeability, which is 100 times greater than that of bentonite and close to that of the EDZ. Interface facing backfill differs from that facing the canister by a different mixture of backfill, EDZ, and big voids. Its permeability is higher than the interface facing the plug because of the high permeability of the backfill compared to that of bentonite in the mixture, and to a small roughness. The three interfaces are largely made up of voids and flow by capillarity in their porous structures is negligible. Their scale parameter $P_r$ in VGMPE model is chosen very small (Table 1), and therefore have negligible $P_{c,e}$ ($P_{c,e} \approx 0$).

Transport properties are such that: the diffusion coefficients of water vapor and hydrogen in unsaturated porous materials are deduced from their diffusion coefficients in free gas and liquid through the proportionality constant (tortuosity factor) calculated by the Millington and Quirk formula²⁵. Diffusion coefficients of hydrogen in free liquid and gas phases are equal to $6.0 \times 10^{-9}$ and $9.5 \times 10^{-5} \text{ m}^2/\text{s}^{17}$, respectively, and the diffusion coefficient of water vapor in free gas phase is equal to $2.1 \times 10^{-5} \text{ m}^2/\text{s}^{15}$. The inverse Henry's constant of hydrogen is equal to $1.379 \times 10^{-10} \text{ Pa}^{-1}$.

The initial gas saturation is 30% for bentonite seals (cell and main drift) and for backfill (drifts). An initial gas saturation of 95% for all interfaces is considered. For these materials initially unsaturated, the gas pressure is equal to 1 atmosphere (~ 1 bar = 0.1 MPa), and the water pressure is deduced from the gas pressure and water saturation. Therefore, liquid pressure in capillary equilibrium with gas pressure is deduced from gas pressure and liquid saturation by applying the VGM-model of the water retention curve defined by Eq. (1). The clay host rock, and cell and drift EDZs are initially at full liquid saturation with an initial water pressure set to 5 MPa.

The mesh of the domain used for the simulations is like that recently presented in¹⁵. It is radial and rectilinear; it consists of $39 \times 122 = 4758$ elements in the plane (r, z). This number is reduced to 4358 active elements because elements representing the waste canister are not considered in the calculation. Two extra-cells with infinite volumes were added to handle boundary conditions (BC). The first one is for fixing the time-dependent gas pressure

<table><caption>**Table 1.** Physical properties and parameters of VGMPE hydrodynamic properties of the materials used in the simulations (¹⁶,¹⁷). aRadial intrinsic permeability. Calculated from horizontal and vertical intrinsic permeability values using geometric mean: ($K_x \times K_y \times K_z$)¹/³.</caption><thead><tr><th rowspan="2">Material</th><th rowspan="2">Porosity (-)</th><th colspan="2">Intrinsic permeability (m²)</th><th rowspan="2">Pore compressibility (Pa⁻¹)</th><th rowspan="2">Capillary entry pressure, P<sub>c,e</sub> (Pa)</th><th rowspan="2">VG pressure scale parameter, P<sub>r</sub> (Pa)</th><th rowspan="2">VG shape parameter, n (-)</th></tr><tr><th>Waste cell scale<sup>a</sup></th><th>Module scale</th></tr></thead><tbody><tr><td>Geological Medium (clay host rock)</td><td>0.15</td><td>7.94 × 10<sup>-21</sup></td><td>1.0 × 10<sup>-20</sup></td><td>1.11 × 10<sup>-9</sup></td><td>2.00 × 10<sup>6</sup></td><td>1.5 × 10<sup>7</sup></td><td>1.5</td></tr><tr><td>EDZ</td><td>0.15</td><td>7.94 × 10<sup>-18</sup></td><td>5.0 × 10<sup>-18</sup></td><td>1.11 × 10<sup>-9</sup></td><td>6.67 × 10<sup>5</sup></td><td>1.5 × 10<sup>6</sup></td><td>1.5</td></tr><tr><td>Bentonite</td><td>0.35</td><td>1.0 × 10<sup>-20</sup></td><td>1.0 × 10<sup>-20</sup></td><td>8.25 × 10<sup>-10</sup></td><td>1.33 × 10<sup>6</sup></td><td>1.6 × 10<sup>7</sup></td><td>1.6</td></tr><tr><td>Backfill</td><td>0.40</td><td>5.0 × 10<sup>-17</sup></td><td>5.0 × 10<sup>-17</sup></td><td>2.09 × 10<sup>-9</sup></td><td>1.00 × 10<sup>5</sup></td><td>2.0 × 10<sup>6</sup></td><td>1.5</td></tr><tr><td>Interface facing plug</td><td>0.30</td><td>7.94 × 10<sup>-18</sup></td><td>5.0 × 10<sup>-18</sup></td><td>1.11 × 10<sup>-9</sup></td><td>0</td><td>1.0 × 10<sup>4</sup></td><td>4</td></tr><tr><td>Interface facing canister</td><td>1.00</td><td>1.0 × 10<sup>-12</sup></td><td>1.0 × 10<sup>-12</sup></td><td>1.09 × 10<sup>-11</sup></td><td>0</td><td>1.0 × 10<sup>4</sup></td><td>4</td></tr><tr><td>Interface facing backfill</td><td>0.40</td><td>1.0 × 10<sup>-15</sup></td><td>1.0 × 10<sup>-15</sup></td><td>2.09 × 10<sup>-9</sup></td><td>0</td><td>1.0 × 10<sup>4</sup></td><td>4</td></tr></tbody></table>
PG and gas saturation SG (Fig. 1b) at the drift exit boundary Z=0 (green discontinuous line in Fig. 1a), and the second one is for fixing a constant liquid pressure PL of 5 MPa at r = 20 m (blue discontinuous line in Fig. 1a).

The purpose is to follow during 100 000 y the state variables gas saturation (SG), and gas and liquid pressures (PG, PL) at points P1-P18, as well as convective and diffusive fluxes of H₂ in liquid and gas phases at sections S<sub>EDZ1</sub>, S<sub>drift</sub> and S<sub>COx</sub>, as shown in Fig. 1a.

**The repository scale "module"**

In this subsection, we aim to develop a 3D model of H₂ transport at the scale of a small repository (module) of many HLW disposal cells (which is only part of a repository of several parallel drifts and modules) by simulating simultaneously all H₂ sources (located in cells), and H₂ transfer pathways in the underground drifts network (access and main drifts) during the post-closure phase of a DGR, i.e. for 100 000 y. The studied module is that presented and documented in FORGE benchmark¹⁴,¹⁷.

The purpose of the study of gas migration within this module is two-fold:

(i) to verify our numerical results simulated for a parametrization P<sub>c,e</sub> = 0 by confronting them to those obtained by other teams that participated to this benchmark. This step is essential for the credibility of our model, which we intend to use in a future work to simulate gas migration scenarios at the scale of a DGR composed of several modules. These scenarios could include different physical sub-models and associated parameters representing the source term of H₂ and its solubility, and thermo-hydraulic-gas properties of the clay host rock, EDZ, engineered barriers materials (backfill and bentonite-plug), and interfaces between EDZ and these materials.

(ii) To investigate the impact of the new parameterization P<sub>c,e</sub> ≠ 0 (for all materials) on gas pressure and water-gas fluxes within the module by using the new VGMPE model for hydraulic properties (Eqs. (1)–(4)). This can help to draw new conclusions about the real situation of gas entry within the host rock and the engineered barriers materials.

The main difficulty in this study is to simulate the physical model with less CPU time. Therefore, the math-
ematical problem should describe properly the physics of two-phase flow in the porous materials composing the
module (disposal waste cell and drifts network). The “conceptual model” includes mathematical descriptions of
physical processes, as well as model geometry, and spatial distribution of materials with different properties, as
well as their spatial discretization (mesh).

The modelled domain

The studied module is in the middle of a raw of 5 modules along Y-axis (having the same geometry) in the whole repository (10 modules, with a raw of 5 modules at each side of the main drift; see Supplementary Fig. S1 online). It consists of 50 cells at each side of the access drift (Fig. 2). Because of the module position within the repository, time-dependent gas pressure PG(t) and gas saturation SG(t) at downstream and upstream BCs (Fig. 2 above) are very similar. These time-dependent PG(t) and gas saturation SG(t), like those shown in Fig. 1b, were calculated using the full 3D repository model. Therefore, the fact that the central XZ-plane of the access drift is a symmetry axis for the disposal cells, only the downstream BC has been modelled (Fig. 2). So, the size of the modelled domain is reduced to the volume 714 × 63 × 150 m along X, Y, and Z-axis, respectively. Figure 3 shows a 3D view of the modelled geometry of the module which has been embedded within the COx. Each HLW disposal cell shown in Fig. 2 has the same characteristics as described in Fig. 1a. The vertical extend of the domain is limited to the layer thickness of COx (150 m). The horizontal extend in the other directions (OX, OY) is representative to the distance between cells (10 m). The porous materials considered in the model are the backfill in drifts, bentonite for plugs (seals) in cells and main drift, the COx clay host rock, the EDZ, and interfaces EDZ-backfill, EDZ-bentonite and COx-bentonite. The so-called “cell” consists of the waste canisters and the bentonite plug of the micro-tunnel (Fig. 1a).

Geometries of the disposal cells and drifts composing the module within the COx are detailed in Wendling et al.¹⁷. However, some approximations have been considered in our simulations:

* Each circular cross section of the disposal cell has been equated to a square cross section (same area) in order to simplify introduction of disposal cells centimetric interfaces in a rectilinear mesh. This approximation is the same as that adopted by other teams that contributed to FORGE benchmark (e.g.²⁶,³²).

* The waste canisters are not modelled and assumed impermeable to water and gas. The elements that represent the canisters are excluded from the calculations. The boundaries to neighboring elements are no-flux boundaries. A constant hydrogen source term of 100 mol/y/cell lasting 10 000 y is attributed to each disposal cell, it is applied at the waste cylinder surface.

* To avoid the severe numerical convergence problems caused by the infinitesimal volume of the centimet-
ric interfaces located near coarser mesh elements, the space-discretization near interfaces is simplified by
homogenizing them with a fraction of one of their adjacent materials, as proposed by²⁶ and applied by²⁷.
In fact, the combination interface-ad adjacent material has been carried out as follows: for cells, the interfaces
EDZ-canister and EDZ-bentonite seal are combined with 0.24 m of EDZ; for drifts, the EDZ-backfill is com-
bined with 0.49 m of backfill while the COx-bentonite interface is combined with 0.49 m of bentonite. This
approach leads to changing transport and hydraulic properties, and initial conditions of the homogenized
layer. A recent analysis¹⁴ showed that differences between results simulated with and without homogenization
of interfaces in the mesh were not substantial.
PG and gas saturation SG (Fig. 1b) at the drift exit boundary Z=0 (green discontinuous line in Fig. 1a), and the second one is for fixing a constant liquid pressure PL of 5 MPa at r = 20 m (blue discontinuous line in Fig. 1a).

The purpose is to follow during 100 000 y the state variables gas saturation (SG), and gas and liquid pressures (PG, PL) at points P1-P18, as well as convective and diffusive fluxes of H₂ in liquid and gas phases at sections S<sub>EDZ1</sub>, S<sub>drift</sub> and S<sub>COx</sub>, as shown in Fig. 1a.

**The repository scale "module"**

In this subsection, we aim to develop a 3D model of H₂ transport at the scale of a small repository (module) of many HLW disposal cells (which is only part of a repository of several parallel drifts and modules) by simulating simultaneously all H₂ sources (located in cells), and H₂ transfer pathways in the underground drifts network (access and main drifts) during the post-closure phase of a DGR, i.e. for 100 000 y. The studied module is that presented and documented in FORGE benchmark¹⁴,¹⁷.

The purpose of the study of gas migration within this module is two-fold:

(i) to verify our numerical results simulated for a parametrization P<sub>c,e</sub> = 0 by confronting them to those obtained by other teams that participated to this benchmark. This step is essential for the credibility of our model, which we intend to use in a future work to simulate gas migration scenarios at the scale of a DGR composed of several modules. These scenarios could include different physical sub-models and associated parameters representing the source term of H₂ and its solubility, and thermo-hydraulic-gas properties of the clay host rock, EDZ, engineered barriers materials (backfill and bentonite-plug), and interfaces between EDZ and these materials.

(ii) To investigate the impact of the new parameterization P<sub>c,e</sub> ≠ 0 (for all materials) on gas pressure and water-gas fluxes within the module by using the new VGMPE model for hydraulic properties (Eqs. (1)–(4)). This can help to draw new conclusions about the real situation of gas entry within the host rock and the engineered barriers materials.

The main difficulty in this study is to simulate the physical model with less CPU time. Therefore, the math-
ematical problem should describe properly the physics of two-phase flow in the porous materials composing the
module (disposal waste cell and drifts network). The “conceptual model” includes mathematical descriptions of
physical processes, as well as model geometry, and spatial distribution of materials with different properties, as
well as their spatial discretization (mesh).

The modelled domain

The studied module is in the middle of a raw of 5 modules along Y-axis (having the same geometry) in the whole repository (10 modules, with a raw of 5 modules at each side of the main drift; see Supplementary Fig. S1 online). It consists of 50 cells at each side of the access drift (Fig. 2). Because of the module position within the repository, time-dependent gas pressure PG(t) and gas saturation SG(t) at downstream and upstream BCs (Fig. 2 above) are very similar. These time-dependent PG(t) and gas saturation SG(t), like those shown in Fig. 1b, were calculated using the full 3D repository model. Therefore, the fact that the central XZ-plane of the access drift is a symmetry axis for the disposal cells, only the downstream BC has been modelled (Fig. 2). So, the size of the modelled domain is reduced to the volume 714 × 63 × 150 m along X, Y, and Z-axis, respectively. Figure 3 shows a 3D view of the modelled geometry of the module which has been embedded within the COx. Each HLW disposal cell shown in Fig. 2 has the same characteristics as described in Fig. 1a. The vertical extend of the domain is limited to the layer thickness of COx (150 m). The horizontal extend in the other directions (OX, OY) is representative to the distance between cells (10 m). The porous materials considered in the model are the backfill in drifts, bentonite for plugs (seals) in cells and main drift, the COx clay host rock, the EDZ, and interfaces EDZ-backfill, EDZ-bentonite and COx-bentonite. The so-called “cell” consists of the waste canisters and the bentonite plug of the micro-tunnel (Fig. 1a).

Geometries of the disposal cells and drifts composing the module within the COx are detailed in Wendling et al.¹⁷. However, some approximations have been considered in our simulations:

* Each circular cross section of the disposal cell has been equated to a square cross section (same area) in order to simplify introduction of disposal cells centimetric interfaces in a rectilinear mesh. This approximation is the same as that adopted by other teams that contributed to FORGE benchmark (e.g.²⁶,³²).

* The waste canisters are not modelled and assumed impermeable to water and gas. The elements that represent the canisters are excluded from the calculations. The boundaries to neighboring elements are no-flux boundaries. A constant hydrogen source term of 100 mol/y/cell lasting 10 000 y is attributed to each disposal cell, it is applied at the waste cylinder surface.

* To avoid the severe numerical convergence problems caused by the infinitesimal volume of the centimet-
ric interfaces located near coarser mesh elements, the space-discretization near interfaces is simplified by
homogenizing them with a fraction of one of their adjacent materials, as proposed by²⁶ and applied by²⁷.
In fact, the combination interface-ad adjacent material has been carried out as follows: for cells, the interfaces
EDZ-canister and EDZ-bentonite seal are combined with 0.24 m of EDZ; for drifts, the EDZ-backfill is com-
bined with 0.49 m of backfill while the COx-bentonite interface is combined with 0.49 m of bentonite. This
approach leads to changing transport and hydraulic properties, and initial conditions of the homogenized
layer. A recent analysis¹⁴ showed that differences between results simulated with and without homogenization
of interfaces in the mesh were not substantial.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Above: 2D view of the module in the plane XY at Z=75 m. Below: 2D view in the XZ plane (cross section AA') of the main drift, and 2D view in the YZ plane (cross section BB') of the access drift and of the disposal waste cell. Access and main drifts have the same squared cross Sect. (6×6 m²). The access drift is completely backfilled (Maya blue color). The main drift is partly backfilled: presence of a bentonite seal with a cross section of 8×8 m² in the plane XZ (orange color). Interfaces of 1 cm thickness are shown by red lines (Adapted from Wendling et al.¹⁷).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The 3D view of the module embedded in the COx (preprocessing under PetraSim).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Above: 2D view of the module in the plane XY at Z=75 m. Below: 2D view in the XZ plane (cross section AA') of the main drift, and 2D view in the YZ plane (cross section BB') of the access drift and of the disposal waste cell. Access and main drifts have the same squared cross Sect. (6×6 m²). The access drift is completely backfilled (Maya blue color). The main drift is partly backfilled: presence of a bentonite seal with a cross section of 8×8 m² in the plane XZ (orange color). Interfaces of 1 cm thickness are shown by red lines (Adapted from Wendling et al.¹⁷).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The 3D view of the module embedded in the COx (preprocessing under PetraSim).</figcaption></figure>
## *The model implementation*

Initial and boundary conditions of the model have been implemented as described in the benchmark specifications¹⁷:

* A Newman type BC with zero-flux is imposed at the lateral boundaries, except at the exit boundary of the main drift, where a Dirichlet BC is imposed with time-dependent gas saturation and gas pressure (idem curves in Fig. 1a).

* At the bottom ($Z=0$) and top ($Z=150$ m) boundaries, Dirichlet BCs are imposed with constant hydrostatic pressures of 4 and 6 MPa, respectively.

* The initial unsaturated conditions in the materials composing the module are the same as those composing the disposal waste cell. However, for the saturated materials initially at full water saturation, clay host rock and EDZs, the initial water pressure is not constant but calculated from hydrostatic equilibrium between $Z=0$ and $Z=150$ m before drifts excavation.

The physical properties, VGMPE-hydrodynamic parameters (Eqs. (1)–(4)), and the two-phase hydrogen diffusion coefficients of the different materials composing the module are the same as those composing the disposal cell (subsection “The cell scale”, Table 1).

Pre-processing of the model mesh and its visualization are carried out via PetraSim²⁸, as shown in Fig. 4a. Figure 4b,c,d show a zoom on the mesh over all the materials which constitute the disposal cells and drifts of the module within the COx.

The mesh of the module is also rectilinear, as for the cell mesh. It is composed of about $412 \times 11 \times 21 = 95172$ elements (including 150 elements representing the mesh of waste canisters that were deactivated, not used for the calculations). Three extra-elements with infinite volumes are added for handling Dirichlet BC's, two of which for imposing constant hydrostatic liquid pressure at bottom and top vertical boundaries ($Z=0, 150$ m), and one of which for imposing time-dependent gas pressure and gas saturation (Fig. 1b) at the main drift exit ($Y=63$ m).

In this study, the state variables gas saturation (SG), and gas and liquid pressure (PG, PL) are followed during 100 000 y at points elements located near cells C1, C25 and C50, and the two points elements P-Pd-1 et P-Pd-2 located in the bentonite-plug of the main drift, Fig. 2. The points elements P-C1-1, P-C25-1, P-C50-1 are situated in the COx; the points elements P-C1-2, P-C25-2, P-C50-2 are situated in the homogenized EDZ-interface (upper part) connected to the bentonite plug; and the points elements P-C1-3, P-C25-3, P-C50-3 are situated in the backfilled access drift. In addition, convective and diffusive $H_2$ fluxes in liquid and gas phases at sections FC-1, F-AD/MD and F-Pd (green, dashed lines in Fig. 2) are also followed during 100 000 y. F-C1: vertical plane passing through the end (close to the waste) of the cell bentonite plug, cell EDZ, and homogenized layer interface-cell EDZ; F-AD/MD: vertical plane at the access drift-main drift interface, passing through backfill, drift EDZ,

<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Mesh of the studied domain. (b) Zoom on the drifts mesh. (c) Mesh of the access drift in plane XZ. (d) Mesh near cells in plane XY.</figcaption></figure>
## *The model implementation*

Initial and boundary conditions of the model have been implemented as described in the benchmark specifications¹⁷:

* A Newman type BC with zero-flux is imposed at the lateral boundaries, except at the exit boundary of the main drift, where a Dirichlet BC is imposed with time-dependent gas saturation and gas pressure (idem curves in Fig. 1a).

* At the bottom ($Z=0$) and top ($Z=150$ m) boundaries, Dirichlet BCs are imposed with constant hydrostatic pressures of 4 and 6 MPa, respectively.

* The initial unsaturated conditions in the materials composing the module are the same as those composing the disposal waste cell. However, for the saturated materials initially at full water saturation, clay host rock and EDZs, the initial water pressure is not constant but calculated from hydrostatic equilibrium between $Z=0$ and $Z=150$ m before drifts excavation.

The physical properties, VGMPE-hydrodynamic parameters (Eqs. (1)-(4)), and the two-phase hydrogen diffusion coefficients of the different materials composing the module are the same as those composing the disposal cell (subsection “The cell scale”, Table 1).

Pre-processing of the model mesh and its visualization are carried out via PetraSim²⁸, as shown in Fig. 4a. Figure 4b,c,d show a zoom on the mesh over all the materials which constitute the disposal cells and drifts of the module within the COx.

The mesh of the module is also rectilinear, as for the cell mesh. It is composed of about 412 × 11 × 21 = 95 172 elements (including 150 elements representing the mesh of waste canisters that were deactivated, not used for the calculations). Three extra-elements with infinite volumes are added for handling Dirichlet BC's, two of which for imposing constant hydrostatic liquid pressure at bottom and top vertical boundaries ($Z=0, 150$ m), and one of which for imposing time-dependent gas pressure and gas saturation (Fig. 1b) at the main drift exit ($Y=63$ m).

In this study, the state variables gas saturation (SG), and gas and liquid pressure (PG, PL) are followed during 100 000 y at points elements located near cells C1, C25 and C50, and the two points elements P-Pd-1 et P-Pd-2 located in the bentonite-plug of the main drift, Fig. 2. The points elements P-C1-1, P-C25-1, P-C50-1 are situated in the COx; the points elements P-C1-2, P-C25-2, P-C50-2 are situated in the homogenized EDZ-interface (upper part) connected to the bentonite plug; and the points elements P-C1-3, P-C25-3, P-C50-3 are situated in the backfilled access drift. In addition, convective and diffusive $H_2$ fluxes in liquid and gas phases at sections FC-1, F-AD/MD and F-Pd (green, dashed lines in Fig. 2) are also followed during 100 000 y. F-C1: vertical plane passing through the end (close to the waste) of the cell bentonite plug, cell EDZ, and homogenized layer interface-cell EDZ; F-AD/MD: vertical plane at the access drift-main drift interface, passing through backfill, drift EDZ,

<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Mesh of the studied domain. (b) Zoom on the drifts mesh. (c) Mesh of the access drift in plane XZ. (d) Mesh near cells in plane XY.</figcaption></figure>
homogenized layer interface-drift EDZ; and F-Pd: vertical plane passing through the middle of the downstream main drift plug and homogenized layer interface-bentonite.

## The numerical solution

Simulations are performed using the TOUGH-suite of codes iTOUGH2⁸ and TOUGH2-MP¹⁸. Both codes have been used to solve many numerical problems of non-isothermal two-phase flow and two components (water-hydrogen) transport in both liquid and gas phases in porous media, especially hydrogen migration at different scales of a DGR (e.g.¹⁴,¹⁵). The advantage of using the massively parallel code TOUGH2-MP is that it uses the domain decomposition method, which is designed for parallel simulation in multi-CPU platforms (multicores processing). The equation of state 5 (EOS5)²⁹ has been chosen to describe thermodynamics of H₂ gas in presence of water.

The numerical solution for each of the two specified models presented above is based on a modified version of TOUGH2-MP/EOS5 by including the VGMPE model as implemented in iTOUGH2/EOS5⁹. This modification led to the development of two numerical approaches:

* The P<sub>c</sub>-method (PCM) which consists of changing the primary variable S<sub>g</sub> by the capillary pressure variable P<sub>c</sub> for the switch from single-phase liquid to two-phase liquid-gas conditions, and vice-versa.

* The S<sub>g</sub>-method (SGM), which introduces thermodynamic conditions for this transition.

Both numerical methods give accurate results by comparison (benchmarking) to an analytical solution and two numerical codes for single-phase liquid unsaturated flow⁹. Moreover, solutions of both numerical schemes were in a good agreement with N₂-gas pressure measured in a borehole of the COx during periods of six pulses of N₂ injection in the PGZ1-experiment carried out at the URL of Bure¹⁰⁻¹². As shown in³⁰, the numerical scheme SGM is better in terms of CPU-time for PGZ1-experiment simulation by comparison to PCM scheme, because it drastically reduces the number of Newton-Raphson iterations and time steps when P<sub>c,e</sub> ≠ 0, and this for both coarse and fine meshes tested. Therefore, the SGM numerical scheme has been chosen for this study.

## Results and discussions

Because of isothermal conditions at 20 °C, the water vapor partial pressure is constant in time and equal to 2337 Pa. This value is very small compared to that of H₂ (between 0.1 and 9 MPa). Idem, the water vapor flux is very negligible compared to that of H₂. Therefore, in what follows, the gas phase is essentially composed of hydrogen.

### The cell scale: P<sub>c,e</sub> = 0 vs P<sub>c,e</sub> ≠ 0 simulations

From the comparison of time-varying liquid and gas pressure (PL and PG), and gas saturation (SG) simulated at points elements P1-P12 (Fig. 1a) during 100 000 y by both parametrizations, P<sub>c,e</sub> = 0 (Fig. 5a,b,c) and P<sub>c,e</sub> ≠ 0 (Fig. 5d,e,f), we can conclude:

- COx: although points P1 and P8-12 remain saturated during the whole simulation period (100 000 y), their liquid pressure PL decreases more during transient phases for the parametrization P<sub>c,e</sub> ≠ 0 (Fig. 5d). This is essentially due to the increase of dissolved H₂ mass generated by canisters metal corrosion that is diffusing in liquid water towards points P1 and P8-12. The gas flow from EDZ (lower gas entry pressure) into intact COx (higher gas entry pressure) is only possible if gas pressure is above gas entry pressure of the COx (development of 2-phase liquid-gas thin zone in the near EDZ-COx interface). This leads to a temporary blockage of H₂ gas at the EDZ-COx interface. The point P12 is much more impacted by the substantial dissolved H₂ mass (negative PL during the period ~10-240 y) due to its presence near the access drift BC (source of H₂, Fig. 1b).

- Interfaces: a fast displacement of gas near the access drift at points P3 and P4 is simulated for a parametrization P<sub>c,e</sub> ≠ 0 (Fig. 5e); the gas saturation SG quickly reaches 1.0. Far from the access drift (point P2), SG reaches 1.0 later after ~100 y (~2000 y for the parametrization P<sub>c,e</sub> = 0).

- EDZ: for the parametrization P<sub>c,e</sub> ≠ 0, there is a significant and fast water-desaturation of the EDZ in the nearest points to the access drift (points P6, P7); there is also a late arrival of H₂ at the farthest point to the access drift, P5 (~400 y, Fig. 5e).

- For points P2 and P3 in the interface and points P5 and P6 in the EDZ located far from the access-drift, the gas pressure increases by about 0.6 MPa at 10 000 y (from 6.47 to 7.1 MPa) for a parametrization P<sub>c,e</sub> ≠ 0 (Fig. 5f). This is due to an important accumulation of H₂ at interfaces and backfill of the access drift for this parametrization.

Looking at Fig. 6a,b,c, one can see, for both parametrizations (P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0), that the transport of dissolved H₂ in liquid-water phase through section SEDZ1 (Fig. 1a) is mainly by diffusion from the EDZ to the clay host rock (COx). Both parametrizations simulate a small convective flux for the transport of dissolved H₂ in liquid phase from the COx towards the EDZ, but with an important delay of ~6000 y for parametrization P<sub>c,e</sub> ≠ 0 before reaching negligible value after 10 000 y. This convective flux is essentially produced by a capillary pressure gradient between the COx (rich of dissolved hydrogen due to diffusion in liquid phase), and interface canister-EDZ which dries faster (liquid saturation of 0%) during the gas production period in the case of parametrization P<sub>c,e</sub> ≠ 0 (see also points P2, P3, P4 in Fig. 5e).
homogenized layer interface-drift EDZ; and F-Pd: vertical plane passing through the middle of the downstream main drift plug and homogenized layer interface-bentonite.

## The numerical solution

Simulations are performed using the TOUGH-suite of codes iTOUGH2⁸ and TOUGH2-MP¹⁸. Both codes have been used to solve many numerical problems of non-isothermal two-phase flow and two components (water-hydrogen) transport in both liquid and gas phases in porous media, especially hydrogen migration at different scales of a DGR (e.g.¹⁴,¹⁵). The advantage of using the massively parallel code TOUGH2-MP is that it uses the domain decomposition method, which is designed for parallel simulation in multi-CPU platforms (multicores processing). The equation of state 5 (EOS5)²⁹ has been chosen to describe thermodynamics of H₂ gas in presence of water.

The numerical solution for each of the two specified models presented above is based on a modified version of TOUGH2-MP/EOS5 by including the VGMPE model as implemented in iTOUGH2/EOS5⁹. This modification led to the development of two numerical approaches:

* The P_c-method (PCM) which consists of changing the primary variable S_g by the capillary pressure variable P_c for the switch from single-phase liquid to two-phase liquid-gas conditions, and vice-versa.

* The S_g-method (SGM), which introduces thermodynamic conditions for this transition.

Both numerical methods give accurate results by comparison (benchmarking) to an analytical solution and two numerical codes for single-phase liquid unsaturated flow⁹. Moreover, solutions of both numerical schemes were in a good agreement with N₂-gas pressure measured in a borehole of the COx during periods of six pulses of N₂ injection in the PGZ1-experiment carried out at the URL of Bure¹⁰⁻¹². As shown in³⁰, the numerical scheme SGM is better in terms of CPU-time for PGZ1-experiment simulation by comparison to PCM scheme, because it drastically reduces the number of Newton-Raphson iterations and time steps when P_c,e ≠ 0, and this for both coarse and fine meshes tested. Therefore, the SGM numerical scheme has been chosen for this study.

## Results and discussions

Because of isothermal conditions at 20 °C, the water vapor partial pressure is constant in time and equal to 2337 Pa. This value is very small compared to that of H₂ (between 0.1 and 9 MPa). Idem, the water vapor flux is very negligible compared to that of H₂. Therefore, in what follows, the gas phase is essentially composed of hydrogen.

### The cell scale: P<sub>c,e</sub> = 0 vs P<sub>c,e</sub> ≠ 0 simulations

From the comparison of time-varying liquid and gas pressure (PL and PG), and gas saturation (SG) simulated at points elements P1-P12 (Fig. 1a) during 100 000 y by both parametrizations, P<sub>c,e</sub> = 0 (Fig. 5a,b,c) and P<sub>c,e</sub> ≠ 0 (Fig. 5d,e,f), we can conclude:

- COx: although points P1 and P8-12 remain saturated during the whole simulation period (100 000 y), their liquid pressure PL decreases more during transient phases for the parametrization P<sub>c,e</sub> ≠ 0 (Fig. 5d). This is essentially due to the increase of dissolved H₂ mass generated by canisters metal corrosion that is diffusing in liquid water towards points P1 and P8-12. The gas flow from EDZ (lower gas entry pressure) into intact COx (higher gas entry pressure) is only possible if gas pressure is above gas entry pressure of the COx (development of 2-phase liquid-gas thin zone in the near EDZ-COx interface). This leads to a temporary blockage of H₂ gas at the EDZ-COx interface. The point P12 is much more impacted by the substantial dissolved H₂ mass (negative PL during the period ~10-240 y) due to its presence near the access drift BC (source of H₂, Fig. 1b).

- Interfaces: a fast displacement of gas near the access drift at points P3 and P4 is simulated for a parametrization P<sub>c,e</sub> ≠ 0 (Fig. 5e); the gas saturation SG quickly reaches 1.0. Far from the access drift (point P2), SG reaches 1.0 later after ~100 y (~2000 y for the parametrization P<sub>c,e</sub> = 0).

- EDZ: for the parametrization P<sub>c,e</sub> ≠ 0, there is a significant and fast water-desaturation of the EDZ in the nearest points to the access drift (points P6, P7); there is also a late arrival of H₂ at the farthest point to the access drift, P5 (~400 y, Fig. 5e).

- For points P2 and P3 in the interface and points P5 and P6 in the EDZ located far from the access-drift, the gas pressure increases by about 0.6 MPa at 10 000 y (from 6.47 to 7.1 MPa) for a parametrization P<sub>c,e</sub> ≠ 0 (Fig. 5f). This is due to an important accumulation of H₂ at interfaces and backfill of the access drift for this parametrization.

Looking at Fig. 6a,b,c, one can see, for both parametrizations (P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0), that the transport of dissolved H₂ in liquid-water phase through section SEDZ1 (Fig. 1a) is mainly by diffusion from the EDZ to the clay host rock (COx). Both parametrizations simulate a small convective flux for the transport of dissolved H₂ in liquid phase from the COx towards the EDZ, but with an important delay of ~6000 y for parametrization P<sub>c,e</sub> ≠ 0 before reaching negligible value after 10 000 y. This convective flux is essentially produced by a capillary pressure gradient between the COx (rich of dissolved hydrogen due to diffusion in liquid phase), and interface canister-EDZ which dries faster (liquid saturation of 0%) during the gas production period in the case of parametrization P<sub>c,e</sub> ≠ 0 (see also points P2, P3, P4 in Fig. 5e).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Evolutions in time of state variables simulated by SGM numerical scheme and the VGMPE model for hydraulic properties. (a, d) Liquid pressure at points P1 and P8-P12; (b, e) gas saturation degree at points P2-P7; (c, f) gas pressure at points P2-P7. (a, b, c) Parametrization $P_{c,e}=0$ (left); (d, e, f) Parametrization $P_{c,e} \neq 0$ (right). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

The most surprising result is that shown in Fig. 6a,b, the diffusion of dissolved H₂ in liquid-water towards the COx remains higher during all the period of gas production (10 000 y) when parametrization $P_{c,e} \neq 0$ is considered. This is essentially due to negligible desaturation of the COx around the EDZ, and the very small diffusive-convective flux of H₂ in the gas phase towards the COx (short, dashed lines in Fig. 6d,e,f). However, for parametrization $P_{c,e} = 0$, both diffusion and convection processes in the gas phase towards the EDZ are substantial (solid lines in Fig. 6e, f). Because of the high $P_{c,e}$-value of the COx, H₂ in the gas phase of EDZ cannot
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Evolutions in time of state variables simulated by SGM numerical scheme and the VGMPE model for hydraulic properties. (a, d) Liquid pressure at points P1 and P8-P12; (b, e) gas saturation degree at points P2-P7; (c, f) gas pressure at points P2-P7. (a, b, c) Parametrization $P_{c,e}=0$ (left); (d, e, f) Parametrization $P_{c,e} \neq 0$ (right). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

The most surprising result is that shown in Fig. 6a,b, the diffusion of dissolved H₂ in liquid-water towards the COx remains higher during all the period of gas production (10 000 y) when parametrization $P_{c,e} \neq 0$ is considered. This is essentially due to negligible desaturation of the COx around the EDZ, and the very small diffusive-convective flux of H₂ in the gas phase towards the COx (short, dashed lines in Fig. 6d,e,f). However, for parametrization $P_{c,e} = 0$, both diffusion and convection processes in the gas phase towards the EDZ are substantial (solid lines in Fig. 6e, f). Because of the high $P_{c,e}$-value of the COx, H₂ in the gas phase of EDZ cannot
<figure><img src="image_8.png" /></figure>

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

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

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

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

<figure><img src="image_13.png" /><figcaption>**Figure 6.** Evolutions in time of H₂ flux simulated at cross section S<sub>EDZ1</sub> (Fig. 1) by using both parametrizations P<sub>c,e</sub>=0 and P<sub>c,e</sub>≠0. (a, d) Convection–Diffusion (CD); (b, e) Diffusion (D); (c, f) Convection (C). (a, b, c) Liquid phase; (d, e, f) Gas phase. INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

deeply penetrate the COx only if the difference between gas pressure in the EDZ and liquid pressure in the COx is greater than the capillary entry pressure of the COx, P<sub>c,e</sub> = 2 MPa. Therefore, H<sub>2</sub> migrates fast through the gas phase of centimetric interfaces and EDZ with low P<sub>c,e</sub>-values towards the backfilled access drift, as shown by the convection dominant transport of H<sub>2</sub> at section SDrift in Fig. 7d,e,f, during the first ~300 y (short, dashed lines, parametrization P<sub>c,e</sub> ≠ 0). Beyond ~300 y, the convective H<sub>2</sub> flux in gas phase simulated by both parametrizations (P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0) is reduced due to the increase of water saturation in the access drift according to the time-dependent gas saturation (Fig. 1b) applied at the access drift boundary (Fig. 1a). Notice, however, that such a reduction is faster for simulations by parametrization P<sub>c,e</sub> ≠ 0. For this parametrization, the convective H<sub>2</sub> flux in gas phase vanishes drastically once the gas production stops at time 10 000 y, whereas a reversible convective H<sub>2</sub> flux in the gas phase from the access drift to the drift EDZ for parametrization P<sub>c,e</sub> = 0 is simulated for a period of about 20 000 y until the access drift becomes fully saturated at ~30 000 y (corresponding to a gas saturation of ~0% in Fig. 1b). This explains the extend of the unsaturated zone in the COx at a radial distance of about 2 m
<figure><img src="image_8.png" /></figure>

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

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

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

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

<figure><img src="image_13.png" /><figcaption>**Figure 6.** Evolutions in time of H₂ flux simulated at cross section S<sub>EDZ1</sub> (Fig. 1) by using both parametrizations P<sub>c,e</sub>=0 and P<sub>c,e</sub>≠0. (a, d) Convection–Diffusion (CD); (b, e) Diffusion (D); (c, f) Convection (C). (a, b, c) Liquid phase; (d, e, f) Gas phase. INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

deeply penetrate the COx only if the difference between gas pressure in the EDZ and liquid pressure in the COx is greater than the capillary entry pressure of the COx, P<sub>c,e</sub> = 2 MPa. Therefore, H<sub>2</sub> migrates fast through the gas phase of centimetric interfaces and EDZ with low P<sub>c,e</sub>-values towards the backfilled access drift, as shown by the convection dominant transport of H<sub>2</sub> at section SDrift in Fig. 7d,e,f, during the first ~300 y (short, dashed lines, parametrization P<sub>c,e</sub> ≠ 0). Beyond ~300 y, the convective H<sub>2</sub> flux in gas phase simulated by both parametrizations (P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0) is reduced due to the increase of water saturation in the access drift according to the time-dependent gas saturation (Fig. 1b) applied at the access drift boundary (Fig. 1a). Notice, however, that such a reduction is faster for simulations by parametrization P<sub>c,e</sub> ≠ 0. For this parametrization, the convective H<sub>2</sub> flux in gas phase vanishes drastically once the gas production stops at time 10 000 y, whereas a reversible convective H<sub>2</sub> flux in the gas phase from the access drift to the drift EDZ for parametrization P<sub>c,e</sub> = 0 is simulated for a period of about 20 000 y until the access drift becomes fully saturated at ~30 000 y (corresponding to a gas saturation of ~0% in Fig. 1b). This explains the extend of the unsaturated zone in the COx at a radial distance of about 2 m
<figure><img src="image_14.png" /></figure>

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

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

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

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

<figure><img src="image_19.png" /><figcaption>**Figure 7.** Evolutions in time of H₂ flux simulated at cross section S<sub>Drift</sub> (Fig. 1) by using both parametrizations P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0. (a, d) Convection–Diffusion (CD); (b, e) Diffusion (D); (c, f) Convection (C). (a, b, c) Liquid phase; (d, e, f) Gas phase. INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

around the EDZ at 10 000 y, as shown in Fig. 8a. For both parametrizations (P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0), there is a small convective transport of dissolved H₂ in liquid-water towards the access drift (section S<sub>Drift</sub>) during almost the whole simulation period. At ~5000 y the transport becomes diffusion-dominated towards drift EDZ, interface plug-EDZ and the plug for parametrization P<sub>c,e</sub> ≠ 0, and later at ~10,000 y for parametrization P<sub>c,e</sub> = 0 (delay of about 5000 y). Notice that the diffusion process in liquid phase is more important for parametrization P<sub>c,e</sub> ≠ 0 due to H₂ confinement in interfaces and EDZ, as also shown in Fig. 6b for the diffusive flux at SEDZ1. Both parameterizations simulate a constant diffusion in liquid phase once the access drift becomes fully saturated at ~30 000 y (Fig. 1b).

The flux of dissolved H₂ in liquid phase and that of H₂ in the gas phase simulated at section SCOx at a radial distance r = 1.8 m (0.8 m away from the waste cell EDZ-COx interface, Fig. 1a), shown in Fig. 9, confirm again
<figure><img src="image_14.png" /></figure>

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

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

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

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

<figure><img src="image_19.png" /><figcaption>**Figure 7.** Evolutions in time of H₂ flux simulated at cross section S<sub>Drift</sub> (Fig. 1) by using both parametrizations P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0. (a, d) Convection–Diffusion (CD); (b, e) Diffusion (D); (c, f) Convection (C). (a, b, c) Liquid phase; (d, e, f) Gas phase. INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

around the EDZ at 10 000 y, as shown in Fig. 8a. For both parametrizations (P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0), there is a small convective transport of dissolved H₂ in liquid-water towards the access drift (section S<sub>Drift</sub>) during almost the whole simulation period. At ~5000 y the transport becomes diffusion-dominated towards drift EDZ, interface plug-EDZ and the plug for parametrization P<sub>c,e</sub> ≠ 0, and later at ~10,000 y for parametrization P<sub>c,e</sub> = 0 (delay of about 5000 y). Notice that the diffusion process in liquid phase is more important for parametrization P<sub>c,e</sub> ≠ 0 due to H₂ confinement in interfaces and EDZ, as also shown in Fig. 6b for the diffusive flux at SEDZ1. Both parameterizations simulate a constant diffusion in liquid phase once the access drift becomes fully saturated at ~30 000 y (Fig. 1b).

The flux of dissolved H₂ in liquid phase and that of H₂ in the gas phase simulated at section SCOx at a radial distance r = 1.8 m (0.8 m away from the waste cell EDZ-COx interface, Fig. 1a), shown in Fig. 9, confirm again
<figure><img src="image_20.png" /><figcaption>**Figure 8.** 2D profiles of gas saturation at time $t = 10\,000\,y$ simulated by our model, by iTOUGH2/EOS5 (numerical scheme SGM) and VGMPE-model for hydraulic properties: (a) Parametrization $P_{c,e} \neq 0$; (b) Parametrization $P_{c,e} = 0$.</figcaption></figure>

simulations observed at section SEDZ1 (Fig. 6). For parametrization $P_{c,e} \neq 0$, the fact that neither diffusion nor convection of $H_2$ in the gas phase are responsible for hydrogen migration to SCOx (Fig. 9d,e,f), all the hydrogen reaching SCOx is that transported by diffusion-only of dissolved $H_2$ in liquid-water from EDZs to the saturated COx, Fig. 9a,b,c. Therefore, section SCOx remains saturated, as also confirmed by gas saturation profile at 10 000 y shown in Fig. 8b. However, for the parametrization $P_{c,e} = 0$, a convection-dominated transport in the gas phase is simulated at section SCOx, from COx towards EDZs and interfaces (cell EDZ-canister, cell EDZ-plug), showing that a two-phase flow is developed from SCOx to the EDZ/COx interface all over the period of gas production 10 000 y (see gas saturation profile at 10 000 y in Fig. 8a, where SCOx at $r = 1.8\,m$ is unsaturated, as well as solid line for convective-diffusive flux at section SDrift in Fig. 7d). This convective flux is reduced in parallel to the increase of water saturation in the access drift, EDZs and interfaces. Also notice the rapid cancellation of flow once gas production stops at 10 000 y for this parametrization ($P_{c,e} = 0$).

Some numerical oscillations that are non-physical appear in the calculation of small liquid fluxes after 10 000 y, especially when parametrization $P_{c,e} = 0$ is considered (whether by convection or by diffusion) in all Figs. 6, 7 and 9. These oscillations, due to numerical dispersion, disappear almost totally for parametrization $P_{c,e} \neq 0$. Vogel et al.¹⁹ made the same observations in their numerical simulations of an experiment of capillary rise of water in a 1D clay column of 1 m length by solving Richards’ equation. Notice, however, the drastic change of convective transport of $H_2$ in liquid phase for parametrization $P_{c,e} \neq 0$ near discontinuity in time at 10 000 y due to small changes in the time step.

To summarize, from analysis of 2D-profiles of gas saturation at 10 000 y (Fig. 8), time-evolution of PL, PG and SG in points P1-P12 (Fig. 5), and diffusive and convective fluxes of $H_2$ in liquid and gas phases through sections SEDZ1, SDrift and SCOx (Figs. 6, 7 and 9), we can conclude that:

- For a parametrization $P_{c,e} = 0$, there is simulation of two-phase flow pathways showing an important advance of $H_2$ front in the COx in the radial direction of the disposal cell, leading to development of an unsaturated zone around EDZs (Fig. 8a). This zone is especially due to important $H_2$ diffusion and convection in liquid and gas phase (supported by Figs. 6, 7 and 9).

- For parametrization $P_{c,e} \neq 0$, the advance of gas front in the COx is limited to the near EDZ/COx interface (Fig. 8b). However, there is simulation of an important dissolved $H_2$ front in the saturated COx due to diffusion-only from the EDZs to the COx as shown in Figs. 6 and 9.

- For parametrization $P_{c,e} \neq 0$, there is a fast water re-saturation of the access drift and of the thin unsaturated zone developed near the disposal cell. This is the consequence of the discontinuity in the VGMPE model
<figure><img src="image_20.png" /><figcaption>**Figure 8.** 2D profiles of gas saturation at time $t = 10\,000\,y$ simulated by our model, by iTOUGH2/EOS5 (numerical scheme SGM) and VGMPE-model for hydraulic properties: (a) Parametrization $P_{c,e} \neq 0$; (b) Parametrization $P_{c,e} = 0$.</figcaption></figure>

simulations observed at section SEDZ1 (Fig. 6). For parametrization $P_{c,e} \neq 0$, the fact that neither diffusion nor convection of $H_2$ in the gas phase are responsible for hydrogen migration to SCOx (Fig. 9d,e,f), all the hydrogen reaching SCOx is that transported by diffusion-only of dissolved $H_2$ in liquid-water from EDZs to the saturated COx, Fig. 9a,b,c. Therefore, section SCOx remains saturated, as also confirmed by gas saturation profile at 10 000 y shown in Fig. 8b. However, for the parametrization $P_{c,e} = 0$, a convection-dominated transport in the gas phase is simulated at section SCOx, from COx towards EDZs and interfaces (cell EDZ-canister, cell EDZ-plug), showing that a two-phase flow is developed from SCOx to the EDZ/COx interface all over the period of gas production 10 000 y (see gas saturation profile at 10 000 y in Fig. 8a, where SCOx at $r = 1.8\,m$ is unsaturated, as well as solid line for convective-diffusive flux at section SDrift in Fig. 7d). This convective flux is reduced in parallel to the increase of water saturation in the access drift, EDZs and interfaces. Also notice the rapid cancellation of flow once gas production stops at 10 000 y for this parametrization ($P_{c,e} = 0$).

Some numerical oscillations that are non-physical appear in the calculation of small liquid fluxes after 10 000 y, especially when parametrization $P_{c,e} = 0$ is considered (whether by convection or by diffusion) in all Figs. 6, 7 and 9. These oscillations, due to numerical dispersion, disappear almost totally for parametrization $P_{c,e} \neq 0$. Vogel et al.¹⁹ made the same observations in their numerical simulations of an experiment of capillary rise of water in a 1D clay column of 1 m length by solving Richards’ equation. Notice, however, the drastic change of convective transport of $H_2$ in liquid phase for parametrization $P_{c,e} \neq 0$ near discontinuity in time at 10 000 y due to small changes in the time step.

To summarize, from analysis of 2D-profiles of gas saturation at 10 000 y (Fig. 8), time-evolution of PL, PG and SG in points P1-P12 (Fig. 5), and diffusive and convective fluxes of $H_2$ in liquid and gas phases through sections SEDZ1, SDrift and SCOx (Figs. 6, 7 and 9), we can conclude that:

- For a parametrization $P_{c,e} = 0$, there is simulation of two-phase flow pathways showing an important advance of $H_2$ front in the COx in the radial direction of the disposal cell, leading to development of an unsaturated zone around EDZs (Fig. 8a). This zone is especially due to important $H_2$ diffusion and convection in liquid and gas phase (supported by Figs. 6, 7 and 9).

- For parametrization $P_{c,e} \neq 0$, the advance of gas front in the COx is limited to the near EDZ/COx interface (Fig. 8b). However, there is simulation of an important dissolved $H_2$ front in the saturated COx due to diffusion-only from the EDZs to the COx as shown in Figs. 6 and 9.

- For parametrization $P_{c,e} \neq 0$, there is a fast water re-saturation of the access drift and of the thin unsaturated zone developed near the disposal cell. This is the consequence of the discontinuity in the VGMPE model
<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" /></figure>

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

<figure><img src="image_26.png" /><figcaption>**Figure 9.** Evolutions in time of H₂ flux simulated at cross section S<sub>COx</sub> (Fig. 1) by using both parametrizations P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0. (a, d) Convection–Diffusion (CD); (b, e) Diffusion (D); (c, f) Convection (C). (a, b, c) Liquid phase; (d, e, f) Gas phase. INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

hydraulic properties near full water saturation (Eqs. (1), (2) and (3)) that are well implemented in the SGM-scheme of TOUGH2/EOS5 code.

- Finally, for parametrization P<sub>c,e</sub> ≠ 0, the bentonite-plug begins to re-saturate earlier, which explains also the abrupt decrease of S<sub>g</sub> from 1 to 0 in point P4 (Fig. 5e).

For the given H₂ generation rate produced by corrosion of waste canisters in the disposal cell during 10 000 y, it is clear from simulations by parametrization P<sub>c,e</sub> ≠ 0 that the COx remains saturated and behaves as a capillary barrier. Therefore, H₂ cannot deeply migrate in the COx through capillary two-phase flow pathways by convection-dominated transport as long as the difference between gas pressure in the EDZ and liquid-water pressure in the COx is less than the P<sub>c,e</sub>-value of the COx. This is due to the highest P<sub>c,e</sub>-value of the COx with regards to the other materials in the modelled domain (Table 1). Because the COx desaturates only in the near interface EDZ-COx, H₂ remains confined in the interfaces and EDZ, and therefore preferentially migrates only
<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" /></figure>

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

<figure><img src="image_26.png" /><figcaption>**Figure 9.** Evolutions in time of H₂ flux simulated at cross section S<sub>COx</sub> (Fig. 1) by using both parametrizations P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0. (a, d) Convection–Diffusion (CD); (b, e) Diffusion (D); (c, f) Convection (C). (a, b, c) Liquid phase; (d, e, f) Gas phase. INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

hydraulic properties near full water saturation (Eqs. (1), (2) and (3)) that are well implemented in the SGM-scheme of TOUGH2/EOS5 code.

- Finally, for parametrization P<sub>c,e</sub> ≠ 0, the bentonite-plug begins to re-saturate earlier, which explains also the abrupt decrease of S<sub>g</sub> from 1 to 0 in point P4 (Fig. 5e).

For the given H₂ generation rate produced by corrosion of waste canisters in the disposal cell during 10 000 y, it is clear from simulations by parametrization P<sub>c,e</sub> ≠ 0 that the COx remains saturated and behaves as a capillary barrier. Therefore, H₂ cannot deeply migrate in the COx through capillary two-phase flow pathways by convection-dominated transport as long as the difference between gas pressure in the EDZ and liquid-water pressure in the COx is less than the P<sub>c,e</sub>-value of the COx. This is due to the highest P<sub>c,e</sub>-value of the COx with regards to the other materials in the modelled domain (Table 1). Because the COx desaturates only in the near interface EDZ-COx, H₂ remains confined in the interfaces and EDZ, and therefore preferentially migrates only
by diffusion-dominated transport in liquid phase deeply in the COx, and by convection-dominated transport in gas phase towards the access drift.

Moreover, on the one hand there is a fast water drainage in the materials with the lowest P<sub>c,e</sub>-values (interfaces, EDZ, backfilled access drift) during the gas accumulation in these materials; on the other hand, there is a quasi-instantaneous re-saturation of these materials after the period of H<sub>2</sub> production (10 000 y). These results confirm again simulation results previously obtained for the PGZ1-experiment<sup>9</sup>.

## The module scale

### *Simulations with P<sub>c,e</sub>=0 and benchmarking*

Our numerical solution converges for the whole simulation period of 100 000 y. Figure 10a,b,c shows time-varying gas pressure at the three points P-C1-3, P-C25-3 and P-C50-3 located in the access drift, along axes of cells #1, #25, #50 (Fig. 2), respectively, which were initially unsaturated ($S_{g,init}$ = 0,3). The gas pressure increases significantly in the access drift and achieves its maximum (~7.57 MPa) at the downtime of gas production by cells (10 000 y). The gas saturation evolves in the same way at each of these points (Fig. 11a,b,c), showing a partial water re-saturation of the drift during ~ 1000 y ($S_g$ = 0,1), then its water desaturation following the gas pressure build-up to achieve $S_g$ = 0,15 at 10 000 y, and finally its total water re-saturation at about ~ 30 000 y ($S_g$ = 0) which corresponds to the transit time to the hydrostatic equilibrium in the drifts network (blue curves BC in Figs. 10 and 11: Time-dependent Dirichlet boundary condition at the main drift exit).

For the points P-C1-2, P-C25-2, and P-C50-2 situated in the homogenized EDZ-interface, which were initially unsaturated, gas pressure curves (Fig. 10a,b,c) follow the same trend in time as in points located in the backfilled access drift (P-C#-3, # = 1, 25, 50), except at early times, during the first 250 y. However, due to the retention properties of the homogenized layer EDZ-interface (like COx), gas saturation at these points elements (P-C#-2, # = 1, 25, 50) is much smaller than those in the backfilled access drift, and on the other it is like that at points P-C1-1, P-C25-1, and P-C50-1 located in the COx (5 m behind the cells, Fig. 2), initially fully water saturated ($S_g$ = 0). Points in the COx begin to desaturate slightly after ~ 1000 y, due to gas advance in the COx (Fig. 11a,b,c) during the period of significant H<sub>2</sub> production by cells. This time evolution is the result of the significant gas pressure

<figure><img src="image_27.png" /><figcaption>**Figure 10.** Evolutions in time of gas pressure simulated by a parametrization P<sub>c,e</sub>=0. (a, b, c) Points P-C#-1,2,3 (#: n° of the cell). (d) Points P-Pd-1,2 in the main drift plug (Fig. 2). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>
by diffusion-dominated transport in liquid phase deeply in the COx, and by convection-dominated transport in gas phase towards the access drift.

Moreover, on the one hand there is a fast water drainage in the materials with the lowest P<sub>c,e</sub>-values (interfaces, EDZ, backfilled access drift) during the gas accumulation in these materials; on the other hand, there is a quasi-instantaneous re-saturation of these materials after the period of H<sub>2</sub> production (10 000 y). These results confirm again simulation results previously obtained for the PGZ1-experiment<sup>9</sup>.

## The module scale

### *Simulations with P<sub>c,e</sub>=0 and benchmarking*

Our numerical solution converges for the whole simulation period of 100 000 y. Figure 10a,b,c shows time-varying gas pressure at the three points P-C1-3, P-C25-3 and P-C50-3 located in the access drift, along axes of cells #1, #25, #50 (Fig. 2), respectively, which were initially unsaturated ($S_{g,init} = 0,3$). The gas pressure increases significantly in the access drift and achieves its maximum (~7.57 MPa) at the downtime of gas production by cells (10 000 y). The gas saturation evolves in the same way at each of these points (Fig. 11a,b,c), showing a partial water re-saturation of the drift during ~1000 y ($S_g=0,1$), then its water desaturation following the gas pressure build-up to achieve $S_g=0,15$ at 10 000 y, and finally its total water re-saturation at about ~30 000 y ($S_g=0$) which corresponds to the transit time to the hydrostatic equilibrium in the drifts network (blue curves BC in Figs. 10 and 11: Time-dependent Dirichlet boundary condition at the main drift exit).

For the points P-C1-2, P-C25-2, and P-C50-2 situated in the homogenized EDZ-interface, which were initially unsaturated, gas pressure curves (Fig. 10a,b,c) follow the same trend in time as in points located in the backfilled access drift (P-C#-3, # = 1, 25, 50), except at early times, during the first 250 y. However, due to the retention properties of the homogenized layer EDZ-interface (like COx), gas saturation at these points elements (P-C#-2, # = 1, 25, 50) is much smaller than those in the backfilled access drift, and on the other it is like that at points P-C1-1, P-C25-1, and P-C50-1 located in the COx (5 m behind the cells, Fig. 2), initially fully water saturated ($S_g=0$). Points in the COx begin to desaturate slightly after ~1000 y, due to gas advance in the COx (Fig. 11a,b,c) during the period of significant H<sub>2</sub> production by cells. This time evolution is the result of the significant gas pressure

<figure><img src="image_27.png" /><figcaption>**Figure 10.** Evolutions in time of gas pressure simulated by a parametrization P<sub>c,e</sub>=0. (a, b, c) Points P-C#-1,2,3 (#: n° of the cell). (d) Points P-Pd-1,2 in the main drift plug (Fig. 2). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>**Figure 11.** Evolutions in time of gas saturation simulated by a parametrization $P_{c,e} = 0$. (a, b, c) Points P–C#-1,2,3 (#: n° of the cell). (d) Points P-Pd-1,2 in the main drift plug (Fig. 2). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

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

evolution at these points (Fig. 10a,b,c), where their maximum values achieve ~7.0 to 7.2 MPa (according to their positions from the access drift exit).

The gas pressure evolves in time in the same way at points P-Pd-1 et P-Pd-2 (overlap of grey and red curves in Fig. 10d) initially unsaturated ($S_g,ini = 0.3$), located in the bentonite-plug of the main drift (Fig. 2). This time-evolution, due to H₂ gas front advance in the main drift from both sides of the bentonite-plug, does not highly differ from that imposed at the exit of the main drift (blue curve: BC). The gas originating from the exit BC slightly increases the maximum gas pressure simulated at 10 000 y (~6.3 MPa). However, the gas saturation degree evolves differently at these two points (Fig. 11d). The point P-Pd-2 saturates prematurely, prior to point P-Pd-1, due to its proximity to the saturated COx. Ultimately, these evolutions show that the bentonite-plug in the main drift has undergone significant and rapid water re-saturation after ~1000 y ($S_g = 0.03$), then a low and very slow water re-saturation after the arrival of gas front, reaching then the minimum value $S_g = 1 \times 10^{-5}$ at time $t \sim 28$ 000 y. The full water re-saturation ($S_g = 0$) of the plug is achieved prior to the transit time to hydrostatic equilibrium (30 000 y), which explains the key role that plays this material in the reduction of gas flux that could reach the exit of the main drift.

The 2D profile of gas saturation degree (Fig. 12a) at slice plane Z = 75 m (vertical mid-height of the COx layer), simulated at time $t = 10$ 000 y, shows a significant desaturation of backfilled drifts and COx around cells. A zoom on some cells near the main drift (see Supplementary Fig. S2b online) shows that an unsaturated zone with high gas saturation (red zone with $S_g \sim 0,2$) is developed in the EDZs of backfilled drifts and cells. Notice, however, the full water re-saturation of the main drift plug with respect to the other backfilled parts of the drifts network (see also 3D profile of $S_g$ in Fig. 12d, and curves $S_g(t)$ in points P-Pd-1 and P-Pd-2 in Fig. 11d). The 2D profile of $S_g$ at slice plane y = 3 m (exit of cells plugs onto the access drift), shown in the Supplementary Fig. S2c online, confirms this result. Within the access drift at slice plane y = 2.5 m (see Supplementary Fig. S2a online), $S_g$ is much higher.

The 2D profile of $S_g$ along the cell axis at slice plane x = 690.5 m, shown in the Supplementary Fig. S2d online, explains well the low desaturation of the COx zone behind the cells, as also shown in Fig. 11a (point P-C1-1).
<figure><img src="image_28.png" /><figcaption>**Figure 11.** Evolutions in time of gas saturation simulated by a parametrization $P_{c,e} = 0$. (a, b, c) Points P–C#-1,2,3 (#: n° of the cell). (d) Points P-Pd-1,2 in the main drift plug (Fig. 2). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

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

evolution at these points (Fig. 10a,b,c), where their maximum values achieve ~7.0 to 7.2 MPa (according to their positions from the access drift exit).

The gas pressure evolves in time in the same way at points P-Pd-1 et P-Pd-2 (overlap of grey and red curves in Fig. 10d) initially unsaturated ($S_g,ini = 0.3$), located in the bentonite-plug of the main drift (Fig. 2). This time-evolution, due to $H_2$ gas front advance in the main drift from both sides of the bentonite-plug, does not highly differ from that imposed at the exit of the main drift (blue curve: BC). The gas originating from the exit BC slightly increases the maximum gas pressure simulated at 10 000 y (~6.3 MPa). However, the gas saturation degree evolves differently at these two points (Fig. 11d). The point P-Pd-2 saturates prematurely, prior to point P-Pd-1, due to its proximity to the saturated COx. Ultimately, these evolutions show that the bentonite-plug in the main drift has undergone significant and rapid water re-saturation after ~1000 y ($S_g = 0.03$), then a low and very slow water re-saturation after the arrival of gas front, reaching then the minimum value $S_g = 1 \times 10^{-5}$ at time $t \sim 28$ 000 y. The full water re-saturation ($S_g = 0$) of the plug is achieved prior to the transit time to hydrostatic equilibrium (30 000 y), which explains the key role that plays this material in the reduction of gas flux that could reach the exit of the main drift.

The 2D profile of gas saturation degree (Fig. 12a) at slice plane Z = 75 m (vertical mid-height of the COx layer), simulated at time $t = 10$ 000 y, shows a significant desaturation of backfilled drifts and COx around cells. A zoom on some cells near the main drift (see Supplementary Fig. S2b online) shows that an unsaturated zone with high gas saturation (red zone with $S_g \sim 0,2$) is developed in the EDZs of backfilled drifts and cells. Notice, however, the full water re-saturation of the main drift plug with respect to the other backfilled parts of the drifts network (see also 3D profile of $S_g$ in Fig. 12d, and curves $S_g(t)$ in points P-Pd-1 and P-Pd-2 in Fig. 11d). The 2D profile of $S_g$ at slice plane y = 3 m (exit of cells plugs onto the access drift), shown in the Supplementary Fig. S2c online, confirms this result. Within the access drift at slice plane y = 2.5 m (see Supplementary Fig. S2a online), $S_g$ is much higher.

The 2D profile of $S_g$ along the cell axis at slice plane x = 690.5 m, shown in the Supplementary Fig. S2d online, explains well the low desaturation of the COx zone behind the cells, as also shown in Fig. 11a (point P-C1-1).
<figure><img src="image_30.png" /><figcaption>**Figure 12.** Profiles of the state variables simulated at time $t=10\,000\,y$ in the module inside the COx with a parametrization $P_{c,e}=0$. 2D view in plane XY at slice plane $Z=75\,m$ of: (a) Gas saturation degree, (b) Pressure and (c) Mass fraction of dissolved $H_2$ in water. (d) 3D view of the gas saturation degree in the whole domain.</figcaption></figure>

The 2D profile of pressure (Fig. 12b) at slice plane $Z=75\,m$, simulated at time $t=10\,000\,y$, shows development of a red zone surrounding all the cells and which is at a pressure of ~7 MPa. Beyond that zone area, an advance of the front of water in both directions OX and OY is simulated, but without any substantial desaturation of the COx. This advance along the X-axis of about 100 m is due to the important solubility of hydrogen in the COx, confirmed as well by the 2D profile of the mass fraction of $H_2$ dissolved in water, simulated at slice plane $Z=75\,m$ at $10\,000\,y$ (Fig. 12c).

The maximum gas pressure obtained by our model at point P-C25-3 at $t \sim 10\,000\,y$ (7.57 MPa, Fig. 10b) doesn't differ from that obtained by²⁷ (see Fig. 11 of²⁷) who used TOUGH2-MP/EOS5 with LaGriT mesh³¹ and homogenization of interfaces with fractions of their adjacent materials as well. Moreover, our model simulates a maximum gas pressure of about 6.3 MPa at point P-Pd-2 (Fig. 10d) which is slightly higher than the one simulated by these authors (5.8 MPa, see Fig. 12 of²⁷). However, some discrepancies were simulated during the gas pressure evolution at the beginning of ~3000 y before maximum gas pressure achievement. In fact, the gas pressure increase from this time is much slower with our model. These discrepancies can be explained by the differences between meshes, between the models of the diffusion coefficient of $H_2$ in a two-phase porous medium, and between methods of the homogenization of interfaces with adjacent materials, used by each team modelers.

Also notice that³² simulated a maximum gas pressure of ~7.8 MPa at point P-C25-3 and of ~6.8 MPa at points P-PU-1 and P-PU-2 (upstream points rather than downstream points P-Pd-1 and P-Pd-2), which are approximately the same as those simulated by our model. Recall here that his³² simulations were carried out using a rectilinear and structured mesh by considering explicitly the centimetric interfaces.

### *Simulations with $P_{c,e} \neq 0*

Numerical simulations with a parametrization $P_{c,e} \neq 0$ did not achieve 100 000 y because of a numerical divergence of the solution at time $t \sim 4122\,y$. When this numerical divergence occurred, the gas pressure build-up in the backfilled access drift (points P-C1-3, P-C25-3, P-C50-3) was very high (more than 18 MPa). The gas saturation degree in the backfill achieved ~0,42. These results are confirmed by the profile of gas saturation degree $S_g$ (Fig. 13a,d), and profile of liquid or gas pressure $P_l$ or $P_g$ at slice plane $Z=75\,m$ (Fig. 13b), simulated at time $t=4122\,y$. The $S_g$-profile at slice plane $Z=75\,m$ (Fig. 13a) also show an unsaturated zone of high gas saturation degree, which envelops all the cells, and which induces a significant advance of an unsaturated front of water in the COx zone over a distance of ~50 m along X-axis for X between 150 and 200 m (gas 'piston effect',¹⁵). This desaturation of the COx is essentially due to the high pressure of $H_2$ which propagates from the backfilled drifts to the COx
<figure><img src="image_30.png" /><figcaption>**Figure 12.** Profiles of the state variables simulated at time $t=10\,000\,y$ in the module inside the COx with a parametrization $P_{c,e}=0$. 2D view in plane XY at slice plane $Z=75\,m$ of: (a) Gas saturation degree, (b) Pressure and (c) Mass fraction of dissolved $H_2$ in water. (d) 3D view of the gas saturation degree in the whole domain.</figcaption></figure>

The 2D profile of pressure (Fig. 12b) at slice plane $Z=75\,m$, simulated at time $t=10\,000\,y$, shows development of a red zone surrounding all the cells and which is at a pressure of ~7 MPa. Beyond that zone area, an advance of the front of water in both directions OX and OY is simulated, but without any substantial desaturation of the COx. This advance along the X-axis of about 100 m is due to the important solubility of hydrogen in the COx, confirmed as well by the 2D profile of the mass fraction of $H_2$ dissolved in water, simulated at slice plane $Z=75\,m$ at $10\,000\,y$ (Fig. 12c).

The maximum gas pressure obtained by our model at point P-C25-3 at $t \sim 10\,000\,y$ (7.57 MPa, Fig. 10b) doesn't differ from that obtained by²⁷ (see Fig. 11 of²⁷) who used TOUGH2-MP/EOS5 with LaGriT mesh³¹ and homogenization of interfaces with fractions of their adjacent materials as well. Moreover, our model simulates a maximum gas pressure of about 6.3 MPa at point P-Pd-2 (Fig. 10d) which is slightly higher than the one simulated by these authors (5.8 MPa, see Fig. 12 of²⁷). However, some discrepancies were simulated during the gas pressure evolution at the beginning of ~3000 y before maximum gas pressure achievement. In fact, the gas pressure increase from this time is much slower with our model. These discrepancies can be explained by the differences between meshes, between the models of the diffusion coefficient of $H_2$ in a two-phase porous medium, and between methods of the homogenization of interfaces with adjacent materials, used by each team modelers.

Also notice that³² simulated a maximum gas pressure of ~7.8 MPa at point P-C25-3 and of ~6.8 MPa at points P-PU-1 and P-PU-2 (upstream points rather than downstream points P-Pd-1 and P-Pd-2), which are approximately the same as those simulated by our model. Recall here that his³² simulations were carried out using a rectilinear and structured mesh by considering explicitly the centimetric interfaces.

### *Simulations with $P_{c,e} \neq 0*

Numerical simulations with a parametrization $P_{c,e} \neq 0$ did not achieve 100 000 y because of a numerical divergence of the solution at time $t \sim 4122\,y$. When this numerical divergence occurred, the gas pressure build-up in the backfilled access drift (points P-C1-3, P-C25-3, P-C50-3) was very high (more than 18 MPa). The gas saturation degree in the backfill achieved ~0,42. These results are confirmed by the profile of gas saturation degree $S_g$ (Fig. 13a,d), and profile of liquid or gas pressure $P_l$ or $P_g$ at slice plane $Z=75\,m$ (Fig. 13b), simulated at time $t=4122\,y$. The $S_g$-profile at slice plane $Z=75\,m$ (Fig. 13a) also show an unsaturated zone of high gas saturation degree, which envelops all the cells, and which induces a significant advance of an unsaturated front of water in the COx zone over a distance of ~50 m along X-axis for X between 150 and 200 m (gas 'piston effect',¹⁵). This desaturation of the COx is essentially due to the high pressure of $H_2$ which propagates from the backfilled drifts to the COx
<figure><img src="image_31.png" /><figcaption>**Figure 13.** Profiles of state variables simulated at time $t=4122$ y in the module inside the COx with a parametrization $P_{c,e} \neq 0$. 2D view in plane XY at slice plane $Z=75$ m of: (a) Gas saturation degree, (b) Pressure, and (c) Mass fraction of dissolved $H_2$ in water. (d) 3D view of the gas saturation degree in the whole domain.</figcaption></figure>

(yellow zone in Fig. 13b) leading to an important advance of the dissolved $H_2$ front in water (Fig. 13c). A zoom on pressure profiles near cells (see Supplementary Fig. S3b,d online) and in the access drift (see Supplementary Fig. S3a,c online) confirm these results. These results can be explained by:

* The mesh, initially developed for parametrization $P_{c,e} = 0$ (Fig. 4), is perhaps not adequate for simulating high pore gas velocities in the access drift, which can be higher than $10^{-5}$ m/s in some zones of the domain. Therefore, new meshes must be developed and tested to avoid the high gas pressure gradients in these regions for a parametrization $P_{c,e} \neq 0$. This new mesh should comply to a maximal Peclet number for all mesh elements with these high pore gas velocities to reduce the numerical dispersion. Moreover, these high gas velocities violate the validity domain of Darcy's law which is initially developed for laminar flows with slow gas velocities.

* The significant $H_2$ pressure calculated in the backfilled drifts (~18 MPa) induces a significant desaturation of the COx with an increase of gas pressure within it (see Supplementary Fig. S3 online). This $H_2$ pressure far exceeds the lithostatic (fracturing) pressure $11.3 \pm 0.3$ MPa of the COx indicated by Andra⁷. This pressure build-up can trigger a micro-fissuring of the COx due to pores dilation, as well as an opening of preferential flow paths which are not accounted by the current model. These micro-fissures and/or fractures can present an escape route for gas in the COx in high quantities, but at the expense of the COx structure changing.

The important accumulation of $H_2$ in the drifts is due to its preferential transport in the homogenized interfaces, backfill, and EDZ, having the smallest $P_{c,e}$-values. The COx, having the highest $P_{c,e}$-value, plays the role of a capillary barrier which does not allow at the beginning creation of two-phase flow pathways in the COx and blocks temporarily the gas phase in the drift until its pressure builds up outer in the COx and exceeds the fracturing pressure of the COx ($11.3 \pm 0.3$ MPa). This mechanism is likely accentuated by the presence of the bentonite-plug in the main drift.

It is worth mentioning, that attempts exist to control the unwanted pressure buildup in form of an Engineered Gas Transport System (EGTS) by modifying essentially the seals and the backfill materials (e.g. bentonite-plug in the main drift, Fig. 2) as proposed by Nagra³³.

As an alternative to SGM-solution which is difficult to converge for no more than 4122 y at the module scale for parametrization $P_{c,e} \neq 0$, simulations have been performed with an approximative solution using a
<figure><img src="image_31.png" /><figcaption>**Figure 13.** Profiles of state variables simulated at time $t=4122$ y in the module inside the COx with a parametrization $P_{c,e} \neq 0$. 2D view in plane XY at slice plane $Z=75$ m of: (a) Gas saturation degree, (b) Pressure, and (c) Mass fraction of dissolved $H_2$ in water. (d) 3D view of the gas saturation degree in the whole domain.</figcaption></figure>

(yellow zone in Fig. 13b) leading to an important advance of the dissolved $H_2$ front in water (Fig. 13c). A zoom on pressure profiles near cells (see Supplementary Fig. S3b,d online) and in the access drift (see Supplementary Fig. S3a,c online) confirm these results. These results can be explained by:

* The mesh, initially developed for parametrization $P_{c,e} = 0$ (Fig. 4), is perhaps not adequate for simulating high pore gas velocities in the access drift, which can be higher than $10^{-5}$ m/s in some zones of the domain. Therefore, new meshes must be developed and tested to avoid the high gas pressure gradients in these regions for a parametrization $P_{c,e} \neq 0$. This new mesh should comply to a maximal Peclet number for all mesh elements with these high pore gas velocities to reduce the numerical dispersion. Moreover, these high gas velocities violate the validity domain of Darcy's law which is initially developed for laminar flows with slow gas velocities.

* The significant $H_2$ pressure calculated in the backfilled drifts (~18 MPa) induces a significant desaturation of the COx with an increase of gas pressure within it (see Supplementary Fig. S3 online). This $H_2$ pressure far exceeds the lithostatic (fracturing) pressure $11.3 \pm 0.3$ MPa of the COx indicated by Andra⁷. This pressure build-up can trigger a micro-fissuring of the COx due to pores dilation, as well as an opening of preferential flow paths which are not accounted by the current model. These micro-fissures and/or fractures can present an escape route for gas in the COx in high quantities, but at the expense of the COx structure changing.

The important accumulation of $H_2$ in the drifts is due to its preferential transport in the homogenized interfaces, backfill, and EDZ, having the smallest $P_{c,e}$-values. The COx, having the highest $P_{c,e}$-value, plays the role of a capillary barrier which does not allow at the beginning creation of two-phase flow pathways in the COx and blocks temporarily the gas phase in the drift until its pressure builds up outer in the COx and exceeds the fracturing pressure of the COx ($11.3 \pm 0.3$ MPa). This mechanism is likely accentuated by the presence of the bentonite-plug in the main drift.

It is worth mentioning, that attempts exist to control the unwanted pressure buildup in form of an Engineered Gas Transport System (EGTS) by modifying essentially the seals and the backfill materials (e.g. bentonite-plug in the main drift, Fig. 2) as proposed by Nagra³³.

As an alternative to SGM-solution which is difficult to converge for no more than 4122 y at the module scale for parametrization $P_{c,e} \neq 0$, simulations have been performed with an approximative solution using a
regularization technique for P<sub>c</sub>(S<sub>i</sub>) (Eq. (1)) linearization near full liquid saturation. Although this approximation can underestimate the real maximum gas pressure within the drift, it can help to explain how P<sub>c,e</sub> can impact, through the VGMPE model, the different transport pathways of H<sub>2</sub> by diffusion and by convection in the module during 100 000 y. By using this approximation, the gas pressure simulated by parametrization P<sub>c,e</sub> ≠ 0 increases in all points elements by comparison to that simulated with parametrization P<sub>c,e</sub> = 0. Figure 14 shows this comparison for the example of points P-C1-1, P-C1-2, P-C1-3 near Cell#1 and for points P-Pd-1 and P-Pd-2 in the main drift plug. The maximum gas pressure (peak achieved at ~ 10 000 y) at the three points elements P-C1-1, P-C1-2, P-C1-3 increases by about 1.5 MPa (Fig. 14a,b,c), while that at points elements P-Pd-1 and P-Pd-2 increases by about 1 MPa (Fig. 14d). This increase in gas pressure is also accompanied by a substantial increase of water drainage and gas saturation in the EDZ (point P-C1-2_Pe, Fig. 15b), and access drift (point P-C1-3_Pe, Fig. 15c).

For the point P-C1-1 located in the clay host rock (COx), Fig. 15a, there is a premature desaturation of the COx during time when parametrization P<sub>c,e</sub> = 0 is considered. However, for parametrization P<sub>c,e</sub> ≠ 0, although desaturation of the COx is delayed, its magnitude is higher than that is simulated by the parametrization P<sub>c,e</sub> = 0. This result is like that obtained by<sup>9</sup> in their study of nitrogen gas migration in the COx by using both parametrizations. However, this result differs from that is simulated in the clay host rock at the cell scale (Fig. 9) for parametrization P<sub>c,e</sub> ≠ 0. This is explained by the high gas pressure that is involved at module scale in the EDZ, and access and main drifts due to higher P<sub>c,e</sub>-values of the COx and bentonite-plug that leaded to the accumulation of gas and its confinement in the materials with low P<sub>c,e</sub> values (EDZ, backfill).

Notice also in Fig. 15 a delay of about 10 000 y before reaching hydrostatic equilibrium (re-saturation) at all points in the EDZ and backfilled drifts after time 10 000 y when parametrization P<sub>c,e</sub> ≠ 0 is considered.

There is a fast re-saturation of the main drift bentonite plug (fast decrease of gas saturation, Fig. 15d) due to the increase of water flux in the plug (Fig. 16f) in the case of parametrization P<sub>c,e</sub> ≠ 0. Moreover, the fact that bentonite-plug has an important P<sub>c,e</sub>-value greater than that of the backfill of the drift, explains this increase.

The analysis of the results for points elements P-C25-# (# = 1,2,3) and P-C50-# (# = 1,2,3) (not shown here, just for avoiding many figures in the text) is like to that for points elements P-C1-# (# = 1,2,3).

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

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

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

<figure><img src="image_35.png" /><figcaption>**Figure 14.** Evolutions in time of gas pressure simulated by parametrizations P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0. (a, b, c) Points P-C#-1,2,3 (#: n° of the cell). (d) Points P-Pd-1,2 in the main drift plug (Fig. 2). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>
regularization technique for P<sub>c</sub>(S<sub>i</sub>) (Eq. (1)) linearization near full liquid saturation. Although this approximation can underestimate the real maximum gas pressure within the drift, it can help to explain how P<sub>c,e</sub> can impact, through the VGMPE model, the different transport pathways of H<sub>2</sub> by diffusion and by convection in the module during 100 000 y. By using this approximation, the gas pressure simulated by parametrization P<sub>c,e</sub> ≠ 0 increases in all points elements by comparison to that simulated with parametrization P<sub>c,e</sub> = 0. Figure 14 shows this comparison for the example of points P-C1-1, P-C1-2, P-C1-3 near Cell#1 and for points P-Pd-1 and P-Pd-2 in the main drift plug. The maximum gas pressure (peak achieved at ~ 10 000 y) at the three points elements P-C1-1, P-C1-2, P-C1-3 increases by about 1.5 MPa (Fig. 14a,b,c), while that at points elements P-Pd-1 and P-Pd-2 increases by about 1 MPa (Fig. 14d). This increase in gas pressure is also accompanied by a substantial increase of water drainage and gas saturation in the EDZ (point P-C1-2_Pe, Fig. 15b), and access drift (point P-C1-3_Pe, Fig. 15c).

For the point P-C1-1 located in the clay host rock (COx), Fig. 15a, there is a premature desaturation of the COx during time when parametrization P<sub>c,e</sub> = 0 is considered. However, for parametrization P<sub>c,e</sub> ≠ 0, although desaturation of the COx is delayed, its magnitude is higher than that is simulated by the parametrization P<sub>c,e</sub> = 0. This result is like that obtained by<sup>9</sup> in their study of nitrogen gas migration in the COx by using both parametrizations. However, this result differs from that is simulated in the clay host rock at the cell scale (Fig. 9) for parametrization P<sub>c,e</sub> ≠ 0. This is explained by the high gas pressure that is involved at module scale in the EDZ, and access and main drifts due to higher P<sub>c,e</sub>-values of the COx and bentonite-plug that leaded to the accumulation of gas and its confinement in the materials with low P<sub>c,e</sub> values (EDZ, backfill).

Notice also in Fig. 15 a delay of about 10 000 y before reaching hydrostatic equilibrium (re-saturation) at all points in the EDZ and backfilled drifts after time 10 000 y when parametrization P<sub>c,e</sub> ≠ 0 is considered.

There is a fast re-saturation of the main drift bentonite plug (fast decrease of gas saturation, Fig. 15d) due to the increase of water flux in the plug (Fig. 16f) in the case of parametrization P<sub>c,e</sub> ≠ 0. Moreover, the fact that bentonite-plug has an important P<sub>c,e</sub>-value greater than that of the backfill of the drift, explains this increase.

The analysis of the results for points elements P-C25-# (# = 1,2,3) and P-C50-# (# = 1,2,3) (not shown here, just for avoiding many figures in the text) is like to that for points elements P-C1-# (# = 1,2,3).

<figure><img src="image_32.png" /><figcaption>**Figure 14.** Evolutions in time of gas pressure simulated by parametrizations P<sub>c,e</sub> = 0 and P<sub>c,e</sub> ≠ 0. (a, b, c) Points P-C#-1,2,3 (#: n° of the cell). (d) Points P-Pd-1,2 in the main drift plug (Fig. 2). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>**Figure 15.** Evolutions in time of gas saturation simulated by parametrizations $P_{c,e}=0$ and $P_{c,e} \neq 0$. (a, b, c) Points P-C#-1,2,3 (#: n° of the cell). (d) Points P-Pd-1,2 in the main drift plug (Fig. 2). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

As shown by curves in Fig. 17a,b,c, there is an important convective transport of dissolved H₂ in liquid water from the homogenized EDZ-interface material to the bentonite plug in the disposal waste cell (FC-C1-L and FC-C1-LPe, Fig. 17c) during the first ~600 y and first 1000 y, for both parametrizations, $P_{c,e} \neq 0$ and $P_{c,e} = 0$, respectively. However, beyond these short periods, a return convective transport of dissolved H₂ towards the homogenized EDZ-interface material is simulated until the gas production ceases at 10 000 y and the cell plug begins to re-saturate. This return flow is more important in the case of parametrization $P_{c,e} \neq 0$. Beyond 10 000 y, this transport is again directed towards the cell plug during a period of about 20 000 y for parametrization $P_{c,e} = 0$, and 10 000 y later for parametrization $P_{c,e} \neq 0$ before vanishing at hydrostatic equilibrium.

For both parametrizations ($P_{c,e} \neq 0$ and $P_{c,e} = 0$), Fig. 17d,e,f, there is a convection-dominated transport in the gas phase from the homogenized EDZ-interface material towards the bentonite plug in the disposal waste cell (FC-C1-G and FC-C1-GPe, Fig. 17f). However, a return transport is simulated earlier for parametrization $P_{c,e} = 0$ (FC-C1-G at ~400 y) than that is simulated for parametrization $P_{c,e} \neq 0$ (FC-C1-GPe at ~3830 y). Both fluxes FC-C1-G and FC-C1-GPe become smaller after ~4500 y. The flux FC-C1-G vanishes earlier at ~30 000 y when a full water saturation is achieved at the main drift exit, whereas it continues for an additional duration of ~10 000 y (till time ~40 000 y) for the flux FC-C1-GPe, before reaching hydrostatic equilibrium, Fig. 15c.

For both parametrizations ($P_{c,e} = 0$ and $P_{c,e} \neq 0$), the simulated transport of dissolved H₂ in liquid water by diffusion towards the AD-exit (fluxes FD-AD/MD-L and FD-AD/MD-LPe, Fig. 18b) is negligible compared to that simulated by convection (fluxes FC-AD/MD-L and FC-AD/MD-LPe, Fig. 18c). There is a parallel increase of both fluxes, FC-AD/MD-L and FC-AD/MD-LPe, towards the main drift during the first ~700 y (with a very small difference, even between peaks). Subsequently, these fluxes decreased differently until a return transport towards the access drift is produced. This return is more significant and much larger in time for FC-AD/MD-LPe, which started earlier at ~1500 y than for FC-AD/MD-L at ~5000 y, until the gas production stops at 10 000 y. Beyond these times, both fluxes FC-AD/MD-L and FC-AD/MD-LPe become small, positive (flux towards the main drift) and practically constant in time (~0.2 to $0.3 \times 10^{-11}$ kg·s⁻¹) before vanishing at ~30650 and ~40770 y, respectively (delay of ~10 120 y). This delay is like that simulated by FC-C1-LPe (Fig. 17c),
<figure><img src="image_33.png" /><figcaption>**Figure 15.** Evolutions in time of gas saturation simulated by parametrizations $P_{c,e}=0$ and $P_{c,e} \neq 0$. (a, b, c) Points P-C#-1,2,3 (#: n° of the cell). (d) Points P-Pd-1,2 in the main drift plug (Fig. 2). BC: boundary condition; INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.</figcaption></figure>

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

As shown by curves in Fig. 17a,b,c, there is an important convective transport of dissolved H₂ in liquid water from the homogenized EDZ-interface material to the bentonite plug in the disposal waste cell (FC-C1-L and FC-C1-LPe, Fig. 17c) during the first ~600 y and first 1000 y, for both parametrizations, $P_{c,e} \neq 0$ and $P_{c,e} = 0$, respectively. However, beyond these short periods, a return convective transport of dissolved H₂ towards the homogenized EDZ-interface material is simulated until the gas production ceases at 10 000 y and the cell plug begins to re-saturate. This return flow is more important in the case of parametrization $P_{c,e} \neq 0$. Beyond 10 000 y, this transport is again directed towards the cell plug during a period of about 20 000 y for parametrization $P_{c,e} = 0$, and 10 000 y later for parametrization $P_{c,e} \neq 0$ before vanishing at hydrostatic equilibrium.

For both parametrizations ($P_{c,e} \neq 0$ and $P_{c,e} = 0$), Fig. 17d,e,f, there is a convection-dominated transport in the gas phase from the homogenized EDZ-interface material towards the bentonite plug in the disposal waste cell (FC-C1-G and FC-C1-GPe, Fig. 17f). However, a return transport is simulated earlier for parametrization $P_{c,e} = 0$ (FC-C1-G at ~400 y) than that is simulated for parametrization $P_{c,e} \neq 0$ (FC-C1-GPe at ~3830 y). Both fluxes FC-C1-G and FC-C1-GPe become smaller after ~4500 y. The flux FC-C1-G vanishes earlier at ~30 000 y when a full water saturation is achieved at the main drift exit, whereas it continues for an additional duration of ~10 000 y (till time ~40 000 y) for the flux FC-C1-GPe, before reaching hydrostatic equilibrium, Fig. 15c.

For both parametrizations ($P_{c,e} = 0$ and $P_{c,e} \neq 0$), the simulated transport of dissolved H₂ in liquid water by diffusion towards the AD-exit (fluxes FD-AD/MD-L and FD-AD/MD-LPe, Fig. 18b) is negligible compared to that simulated by convection (fluxes FC-AD/MD-L and FC-AD/MD-LPe, Fig. 18c). There is a parallel increase of both fluxes, FC-AD/MD-L and FC-AD/MD-LPe, towards the main drift during the first ~700 y (with a very small difference, even between peaks). Subsequently, these fluxes decreased differently until a return transport towards the access drift is produced. This return is more significant and much larger in time for FC-AD/MD-LPe, which started earlier at ~1500 y than for FC-AD/MD-L at ~5000 y, until the gas production stops at 10 000 y. Beyond these times, both fluxes FC-AD/MD-L and FC-AD/MD-LPe become small, positive (flux towards the main drift) and practically constant in time (~0.2 to $0.3 \times 10^{-11}$ kg·s⁻¹) before vanishing at ~30650 and ~40770 y, respectively (delay of ~10 120 y). This delay is like that simulated by FC-C1-LPe (Fig. 17c),
<figure><img src="image_37.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(f)</figcaption></figure>

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

**Figure 18.** Evolutions in time of H₂ flux simulated at cross section F-AD/MD (Fig. 2) by using both parametrizations P<sub>c,e</sub>=0 and P<sub>c,e</sub>≠0. (a, d) Convection–Diffusion (CD); (b, e) Diffusion (D); (c, f) Convection (C). (a, b, c) Liquid phase; (d, e, f) Gas phase. INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.

first directed towards the main drift exit BC during the periods ~25 000 y and ~40 000 y for parametrizations P<sub>c,e</sub>=0 and P<sub>c,e</sub>≠0, respectively, before a return diffusive transport towards the access drift begins when the exit boundary begins to be fully water saturated.

Figure 19c shows that the convective transport of dissolved H₂ in liquid phase simulated by both parametrizations, although negligible compared to that of the diffusion (Fig. 19b), is directed towards the exit BC during almost all the period of gas generation (10 000 y). Notice also that FC-Pd-LPe is higher than FC-Pd-L during this period, but both fluxes decrease by increasing water saturation at the exit boundary until reaching small return flow rates from the plug towards the access drift before vanishing entirely at ~32 000 y and ~42 000 y for parametrizations P<sub>c,e</sub>=0 and P<sub>c,e</sub>≠0, respectively (delay of about 10 000 y).
<caption>(a)</caption>

F-AD/MD

<figure><img src="image_35.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(e)</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(f)</figcaption></figure>

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

**Figure 18.** Evolutions in time of H₂ flux simulated at cross section F-AD/MD (Fig. 2) by using both parametrizations P<sub>c,e</sub>=0 and P<sub>c,e</sub>≠0. (a, d) Convection–Diffusion (CD); (b, e) Diffusion (D); (c, f) Convection (C). (a, b, c) Liquid phase; (d, e, f) Gas phase. INJ: injected mass of hydrogen per unit time, generated by the corrosion of canisters of the disposal waste cell during 10 000 y.

first directed towards the main drift exit BC during the periods ~25 000 y and ~40 000 y for parametrizations P<sub>c,e</sub>=0 and P<sub>c,e</sub>≠0, respectively, before a return diffusive transport towards the access drift begins when the exit boundary begins to be fully water saturated.

Figure 19c shows that the convective transport of dissolved H₂ in liquid phase simulated by both parametrizations, although negligible compared to that of the diffusion (Fig. 19b), is directed towards the exit BC during almost all the period of gas generation (10 000 y). Notice also that FC-Pd-LPe is higher than FC-Pd-L during this period, but both fluxes decrease by increasing water saturation at the exit boundary until reaching small return flow rates from the plug towards the access drift before vanishing entirely at ~32 000 y and ~42 000 y for parametrizations P<sub>c,e</sub>=0 and P<sub>c,e</sub>≠0, respectively (delay of about 10 000 y).
BC becomes saturated. This shows that mass of hydrogen reaching the drift plug, and thus the exit BC, through gas phase, by using parametrization $P_{c,e} \neq 0$ is more significant than that is simulated by using parametrization $P_{c,e} = 0$. This delay, as well as significant $H_2$ mass reaching the plug simulated by parametrization $P_{c,e} \neq 0$ are the consequences of higher $P_{c,e}$ values of clay host rock (COx) and drift plug bentonite (2 and 1.33 MPa, respectively, Table 1) which do not ease the escape of the cumulated gas in the backfilled drift network, and therefore lead to the increase of gas pressure and gas flow within this network. Gas entry in the drift plug is delayed until the difference between gas pressure in the drift backfill and liquid pressure in the drift plug becomes higher than the $P_{c,e}$-value of the drift bentonite plug.

The significant convective gas transport towards the main drift plug and the exit BC for parametrization $P_{c,e} \neq 0$, explains the importance of gas piston effect on water flow shown in Fig. 16f for this parametrization.

## Conclusion

This study is the first one to show the impact of a non-zero gas-entry pressure on the results of gas migration modeling in the various components of a DGR (cell and module). Illustration of this impact through FORGE benchmark models for both HLW cell and module shows that gas accumulates mainly in materials having the lowest $P_{c,e}$ values (interfaces, EDZ and backfill), which can substantially increase the gas pressure within the backfilled drift network and its neighborhood (host rock). In the case of FORGE benchmarks, the gas pressure can reach very high values, thus taking the model out of its range of validity and requiring additional complex hydromechanical processes to be considered, such as pathway dilations, opening of interfaces or cracking/fracturing of the host rock.

Enhancement of the numerical modelling of gas migration within the different components of a DGR may also require consideration of processes other than the gas-entry pressure, such as hysteresis, thermal, chemical, and hydro-gas-mechanical processes. These processes might influence porosity and pore size distribution, which would change parameterization of capillary-pressure—saturation—permeability—diffusivity relations.

Finally, all these enhancements could help in future studies to evaluate carefully advective and diffusive transport of radionuclides within the backfilled drifts.

## Data availability

All data that support the findings of this study are available from the corresponding author on request.

Received: 31 July 2023; Accepted: 5 March 2024

Published online: 14 March 2024

## References

1. Hildenbrand, A., Schloemer, S. & Krooss, B. M. Gas breakthrough experiments on fine-grained sedimentary rocks. *Geofluids* **2**, 3–23 (2002).

2. Marschall, P., Horseman, S. & Gimmi, T. Characterisation of gas transport properties of the Opalinus clay, a potential host rock formation for radioactive waste disposal. *Oil Gas Sci. Technol.* **60**, 121–139 (2005).

3. Egermann, P., Lombard, J.-M. & Bretonnier, P. A fast and accurate method to measure threshold capillary pressure of caprocks under representative conditions. Society of core analysts. In *International symposium Society of Core Analysts, Trondheim, Norway, September 2006 Vol. SCA2006-07* (2006).

4. Boulin, P., Bretonnier, P., Vassil, V., Samouillet, A., Fleury, M. & Lombard, J.-M. Entry pressure measurements using three unconventional experimental methods. In *International Symposium of the Society of Core Analysts, Austin, Texas, USA, 18–21 September 2011* (2011).

5. Cuss, R. J., Harrington, J. F., Noy, D. J., Graham, C. C. & Sellin, P. Evidence of localised gas propagation pathways in a field-scale bentonite engineered barrier system; results from three gas injection tests in the large scale gas injection test (Lasgit). *Appl. Clay Sci.* **102**, 81–92 (2014).

6. Harrington, J. F., Graham, C. C., Cuss, R. J. & Norris, S. Gas network development in a precompacted bentonite experiment: Evidence of generation and evolution. *Appl. Clay Sci.* **147**, 80–89 (2017).

7. Andra, S. Evaluation de la faisabilité du stockage géologique en formation argileuse. Dossier 2005 Argile, https://www.andra.fr/sites/default/files/2017-12/266.pdf (2005).

8. Finsterle, S. iTOUGH2 user's guide. LBNL-40040, Earth Sciences Division, LBNL, University of California, Berkeley, CA 94720 (2016).

9. Amri, A., Saâdi, Z. & Ababou, R. Parametric sensitivity to gas-entry pressure in two-phase flow models for a deep geologic repository of radioactive waste. *Transp. Porous Med.* **145**, 13–43 (2022).

10. de La Vaissière, R. Interpretation of in-situ gas experiment PGZ1. FORGE Report D5.9 & D5.12, 6 December 2011, 1–64 (2011).

11. de La Vaissière, R. et al. Gas injection test in the Callovo-Oxfordian claystone: Data analysis and numerical modelling. *Geol. Soc. Lond. Special Publ.* **400**, 427–441 (2014).

12. de La Vaissière, R., Talandier, J., Armand, G., Vu, M.-N. & Cornet F. H. From two-phase flow to gas fracturing into Callovo-Oxfordian claystone. In *53rd US Rock Mechanics/Geomechanics Symposium, American Rock Mechanics Association, New York City, New York, June 2019* (2019).

13. Enssle, C. P., Brommundt, J., Kaempfer, Th. U., Mayer, G. & Wendling, J. Full-scale 3D modelling of a nuclear waste repository in the Callovo-Oxfordian clay. Part 2: thermo-hydraulic two phase transport of water, hydrogen, 14C and 129I. *Geol. Soc. Lond. Special Publ.* **400**, 469–481 (2014).

14. Wendling, J. et al. Gas transport modelling at different spatial scales of a geological repository in clay host rock. *Environ. Earth Sci.* **78**, 221 (2019).

15. Saâdi, Z., Dymitrowska, M., Deleruyelle, F. & Marsal, F. An evaluation model of the impact of hydrogen 'piston effect' on water displacement in a deep geological repository for radioactive waste. *Environ. Earth Sci.* **79**, 434 (2020).

16. Wendling, J. et al. Final report on benchmark studies on repository-scale numerical simulations of gas migration. Part 1: Cell scale benchmark. European Commission FORGE Deliverable D1.6-R (2013).

17. Wendling, J. et al. Final report on benchmark studies on repository-scale numerical simulations of gas migration, Part 2: Module scale benchmark. FORGE Deliverable Report D1.6-R. Luxembourg: European Commission (2013).

18. Zhang, K., Wu, Y.-S. & Pruess, K. User's guide for TOUGH2-MP—A massively parallel version of the TOUGH2 code, report LBNL-315E, Lawrence Berkeley National Laboratory, Berkeley, California, USA (2008).
BC becomes saturated. This shows that mass of hydrogen reaching the drift plug, and thus the exit BC, through gas phase, by using parametrization $P_{c,e} \neq 0$ is more significant than that is simulated by using parametrization $P_{c,e} = 0$. This delay, as well as significant $H_2$ mass reaching the plug simulated by parametrization $P_{c,e} \neq 0$ are the consequences of higher $P_{c,e}$ values of clay host rock (COx) and drift plug bentonite (2 and 1.33 MPa, respectively, Table 1) which do not ease the escape of the cumulated gas in the backfilled drift network, and therefore lead to the increase of gas pressure and gas flow within this network. Gas entry in the drift plug is delayed until the difference between gas pressure in the drift backfill and liquid pressure in the drift plug becomes higher than the $P_{c,e}$-value of the drift bentonite plug.

The significant convective gas transport towards the main drift plug and the exit BC for parametrization $P_{c,e} \neq 0$, explains the importance of gas piston effect on water flow shown in Fig. 16f for this parametrization.

## Conclusion

This study is the first one to show the impact of a non-zero gas-entry pressure on the results of gas migration modeling in the various components of a DGR (cell and module). Illustration of this impact through FORGE benchmark models for both HLW cell and module shows that gas accumulates mainly in materials having the lowest $P_{c,e}$ values (interfaces, EDZ and backfill), which can substantially increase the gas pressure within the backfilled drift network and its neighborhood (host rock). In the case of FORGE benchmarks, the gas pressure can reach very high values, thus taking the model out of its range of validity and requiring additional complex hydromechanical processes to be considered, such as pathway dilations, opening of interfaces or cracking/fracturing of the host rock.

Enhancement of the numerical modelling of gas migration within the different components of a DGR may also require consideration of processes other than the gas-entry pressure, such as hysteresis, thermal, chemical, and hydro-gas-mechanical processes. These processes might influence porosity and pore size distribution, which would change parameterization of capillary-pressure-saturation-permeability-diffusivity relations.

Finally, all these enhancements could help in future studies to evaluate carefully advective and diffusive transport of radionuclides within the backfilled drifts.

## Data availability

All data that support the findings of this study are available from the corresponding author on request.

Received: 31 July 2023; Accepted: 5 March 2024

Published online: 14 March 2024

## References

1. Hildenbrand, A., Schloemer, S. & Krooss, B. M. Gas breakthrough experiments on fine-grained sedimentary rocks. *Geofluids* **2**, 3-23 (2002).

2. Marschall, P., Horseman, S. & Gimmi, T. Characterisation of gas transport properties of the Opalinus clay, a potential host rock formation for radioactive waste disposal. *Oil Gas Sci. Technol.* **60**, 121-139 (2005).

3. Egermann, P., Lombard, J.-M. & Bretonnier, P. A fast and accurate method to measure threshold capillary pressure of caprocks under representative conditions. Society of core analysts. In *International symposium Society of Core Analysts, Trondheim, Norway, September 2006 Vol. SCA2006-07* (2006).

4. Boulin, P., Bretonnier, P., Vassil, V., Samouillet, A., Fleury, M. & Lombard, J.-M. Entry pressure measurements using three unconventional experimental methods. In *International Symposium of the Society of Core Analysts, Austin, Texas, USA, 18-21 September 2011* (2011).

5. Cuss, R. J., Harrington, J. F., Noy, D. J., Graham, C. C. & Sellin, P. Evidence of localised gas propagation pathways in a field-scale bentonite engineered barrier system; results from three gas injection tests in the large scale gas injection test (Lasgit). *Appl. Clay Sci.* **102**, 81-92 (2014).

6. Harrington, J. F., Graham, C. C., Cuss, R. J. & Norris, S. Gas network development in a precompacted bentonite experiment: Evidence of generation and evolution. *Appl. Clay Sci.* **147**, 80-89 (2017).

7. Andra, S. Evaluation de la faisabilité du stockage géologique en formation argileuse. Dossier 2005 Argile, https://www.andra.fr/sites/default/files/2017-12/266.pdf (2005)

8. Finsterle, S. iTOUGH2 user's guide. LBNL-40040, Earth Sciences Division, LBNL, University of California, Berkeley, CA 94720 (2016).

9. Amri, A., Saâdi, Z. & Ababou, R. Parametric sensitivity to gas-entry pressure in two-phase flow models for a deep geologic repository of radioactive waste. *Transp. Porous Med.* **145**, 13-43 (2022).

10. de La Vaissière, R. Interpretation of in-situ gas experiment PGZ1. FORGE Report D5.9 & D5.12, 6 December 2011, 1-64 (2011).

11. de La Vaissière, R. et al. Gas injection test in the Callovo-Oxfordian claystone: Data analysis and numerical modelling. *Geol. Soc. Lond. Special Publ.* **400**, 427-441 (2014).

12. de La Vaissière, R., Talandier, J., Armand, G., Vu, M.-N. & Cornet F. H. From two-phase flow to gas fracturing into Callovo-Oxfordian claystone. In *53rd US Rock Mechanics/Geomechanics Symposium, American Rock Mechanics Association, New York City, New York, June 2019* (2019).

13. Enssle, C. P., Brommundt, J., Kaempfer, Th. U., Mayer, G. & Wendling, J. Full-scale 3D modelling of a nuclear waste repository in the Callovo-Oxfordian clay. Part 2: thermo-hydraulic two phase transport of water, hydrogen, 14C and 129I. *Geol. Soc. Lond. Special Publ.* **400**, 469-481 (2014).

14. Wendling, J. et al. Gas transport modelling at different spatial scales of a geological repository in clay host rock. *Environ. Earth Sci.* **78**, 221 (2019).

15. Saâdi, Z., Dymitrowska, M., Deleruyelle, F. & Marsal, F. An evaluation model of the impact of hydrogen 'piston effect' on water displacement in a deep geological repository for radioactive waste. *Environ. Earth Sci.* **79**, 434 (2020).

16. Wendling, J. et al. Final report on benchmark studies on repository-scale numerical simulations of gas migration. Part 1: Cell scale benchmark. European Commission FORGE Deliverable D1.6-R (2013).

17. Wendling, J. et al. Final report on benchmark studies on repository-scale numerical simulations of gas migration, Part 2: Module scale benchmark. FORGE Deliverable Report D1.6-R. Luxembourg: European Commission (2013).

18. Zhang, K., Wu, Y.-S. & Pruess, K. User's guide for TOUGH2-MP-A massively parallel version of the TOUGH2 code, report LBNL-315E, Lawrence Berkeley National Laboratory, Berkeley, California, USA (2008).
19. Vogel, T., Van Genuchten, MTh. & Cislerova, M. Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions. *Adv. Water Resour.* **24**, 133–144 (2001).

20. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. *Water Resour. Res.* **12**, 513–522 (1976).

21. van Genuchten, MTh. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Sci. Soc. Am. J.* **44**, 892–898 (1980).

22. Parker, J., Lenhard, R. & Kuppusamy, T. A parametric model for constitutive properties governing multiphase flow in porous media. *Water Resour. Res.* **23**, 618–624 (1987).

23. Andra. Référentiel du comportement THM des formations sur le site de Meuse/Haute Marne, Document technique D.RP. AMFS.12.0024, 2012, Andra (2012).

24. Andra. Mesure de la pression de percée de gaz et importance de la cinétique de montée en pression. (Document technique externe ANDRA, CG.RP.AMFS.12.0030). Agence Nationale pour la Gestion des Déchets Radioactifs (ANDRA), France, rapp. tech., 2012 (Juillet), 1–28 (2012).

25. Millington, R. J. & Quirk, J. M. Permeability of porous solids. *Trans. Faraday Soc.* **57**, 1200–1207 (1961).

26. Ahusborde, E., Amaziane, B. & Jurak, M. Three-dimensional numerical simulation by upscaling of gas migration through engineered and geological barriers for a deep repository for radioactive waste. *Geol. Soc. Lond. Special Publ.* **415**, 123–141 (2015).

27. Sentis, M. L. & Gable, C. W. Coupling LaGrit unstructured mesh generation and model setup with TOUGH2 flow and transport: A case study. *Comput. Geosci.* **108**, 42–49 (2017).

28. Thunderhead engineering. PetraSim user manual, Version 5.0 (2010).

29. Pruess, K., Oldenburg, C. & Moridis, G. TOUGH2 user's guide, version 2.1, report LBNL-43134, Lawrence Berkeley National Laboratory, Berkeley, California, USA (2012).

30. Amri, A. Réévaluation du modèle physique de transfert de l'hydrogène pour l'étude du transitoire hydraulique-gaz dans un stockage profond de déchets radioactifs. PhD Thesis, University of Toulouse, Institut National Polytechnique de Toulouse (INP Toulouse), 213 (2021).

31. LaGriT—Los Alamos Grid Toolbox. Los Alamos National Laboratory, http://lagrit.lanl.gov (2016).

32. Justinavicius, D. Investigation of gas migration in a geological repository. Summary of Doctoral Dissertation. Technological Sciences, Energetics and Power Engineering (06T), Lithuanian Energy Institute, Nuclear Engineering Laboratory, Kaunas, Lithuania, 36 (2014)

33. Nagra. Effects of post-disposal gas generation in a repository for low- and intermediate-level waste sited in the Opalinus Clay of Northern Switzerland. Nagra Technical Report NTB 08-07 No. 08-07. Nagra: Wettingen, Switzerland, 175 (2008).

## Acknowledgements

The author gratefully acknowledges the funding provided by EURAD, the European Joint Programme on Radioactive Waste Management for DONUT (Grant Agreement No 847593). I also thank both reviewers for their insightful comments, which significantly improved the paper.

## Author contributions

S.Z. preprocessed geometries, meshes, initial and boundary conditions (input files). He did the numerical simulations by TOUGH codes and post-processed the results (output files). He also wrote the main manuscript text and prepared all the figures.

## Competing interests

The author declares no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-024-56454-y.

**Correspondence** and requests for materials should be addressed to Z.S.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024
19. Vogel, T., Van Genuchten, MTh. & Cislerova, M. Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions. *Adv. Water Resour.* **24**, 133–144 (2001).

20. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. *Water Resour. Res.* **12**, 513–522 (1976).

21. van Genuchten, MTh. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Sci. Soc. Am. J.* **44**, 892–898 (1980).

22. Parker, J., Lenhard, R. & Kuppusamy, T. A parametric model for constitutive properties governing multiphase flow in porous media. *Water Resour. Res.* **23**, 618–624 (1987).

23. Andra. Référentiel du comportement THM des formations sur le site de Meuse/Haute Marne, Document technique D.RP. AMFS.12.0024, 2012, Andra (2012).

24. Andra. Mesure de la pression de percée de gaz et importance de la cinétique de montée en pression. (Document technique externe ANDRA, CG.RP.AMFS.12.0030). Agence Nationale pour la Gestion des Déchets Radioactifs (ANDRA), France, rapp. tech., 2012 (Juillet), 1–28 (2012).

25. Millington, R. J. & Quirk, J. M. Permeability of porous solids. *Trans. Faraday Soc.* **57**, 1200–1207 (1961).

26. Ahusborde, E., Amaziane, B. & Jurak, M. Three-dimensional numerical simulation by upscaling of gas migration through engineered and geological barriers for a deep repository for radioactive waste. *Geol. Soc. Lond. Special Publ.* **415**, 123–141 (2015).

27. Sentis, M. L. & Gable, C. W. Coupling LaGrit unstructured mesh generation and model setup with TOUGH2 flow and transport: A case study. *Comput. Geosci.* **108**, 42–49 (2017).

28. Thunderhead engineering. PetraSim user manual, Version 5.0 (2010).

29. Pruess, K., Oldenburg, C. & Moridis, G. TOUGH2 user's guide, version 2.1, report LBNL-43134, Lawrence Berkeley National Laboratory, Berkeley, California, USA (2012).

30. Amri, A. Réévaluation du modèle physique de transfert de l'hydrogène pour l'étude du transitoire hydraulique-gaz dans un stockage profond de déchets radioactifs. PhD Thesis, University of Toulouse, Institut National Polytechnique de Toulouse (INP Toulouse), 213 (2021).

31. LaGriT—Los Alamos Grid Toolbox. Los Alamos National Laboratory, http://lagrit.lanl.gov (2016).

32. Justinavicius, D. Investigation of gas migration in a geological repository. Summary of Doctoral Dissertation. Technological Sciences, Energetics and Power Engineering (06T), Lithuanian Energy Institute, Nuclear Engineering Laboratory, Kaunas, Lithuania, 36 (2014)

33. Nagra. Effects of post-disposal gas generation in a repository for low- and intermediate-level waste sited in the Opalinus Clay of Northern Switzerland. Nagra Technical Report NTB 08-07 No. 08-07. Nagra: Wettingen, Switzerland, 175 (2008).

## Acknowledgements

The author gratefully acknowledges the funding provided by EURAD, the European Joint Programme on Radioactive Waste Management for DONUT (Grant Agreement No 847593). I also thank both reviewers for their insightful comments, which significantly improved the paper.

## Author contributions

S.Z. preprocessed geometries, meshes, initial and boundary conditions (input files). He did the numerical simulations by TOUGH codes and post-processed the results (output files). He also wrote the main manuscript text and prepared all the figures.

## Competing interests

The author declares no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-024-56454-y.

**Correspondence** and requests for materials should be addressed to Z.S.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024