Corrosion of steel canisters, disposed of in a repository for high-level waste (HLW), leads to generation of hydrogen gas for a long period after the repository's closure. The accumulation of hydrogen gas may lead to significant desaturation and unacceptable build-up of pressure in the backfilled disposal tunnels if the gas cannot escape through the low-permeability host rock. Consequently, the investigation of gas migration is of high relevance in the assessment of the repository's performance.

In this paper, the results of numerical investigations on gas migration performed using the computer code TOUGH2 (USA) are presented. The objective was to investigate migration of gas generated in a single disposal tunnel of a conceptual geological repository in a clay formation, which was suggested for benchmark studies in the European Commission project FORGE (Fate Of Repository GasEs). The analysis was focused on evaluation of the impact of an initial temperature in the repository and of different tortuosity models on gas migration. It was revealed that gas migration results were dependent on tortuosity model, while temperature variation in the repository had minor impact.

Gas influence on the long-term safety of geological repositories has been investigated for more than fifteen years. EVAGAS (Manai, 1997), MEGAS (Ortiz et al., 1997), PEGASUS (Haijtink and Rodwell, 1998), PROGRESS (Rodwell, 2000), GASNET (Rodwell et al., 2003), PAMINA (Grupa and Schröder, 2009) and FORGE (Norris, 2014) are international projects which were undertaken to investigate the gas problem in the context of the safety of the disposal of radioactive waste. During these investigations it was established that radioactive waste and their disposal packages (steel canisters), as well as reinforcement steel used in the construction of the repository, are the principal sources of gas generation in the case of high-level waste and spent nuclear fuel disposal. Among four possible ways of gas generation (corrosion of metal components, microbial degradation, radiolysis and radioactive decay), anaerobic corrosion of steel components in the engineered barrier system (EBS) is the dominant source of hydrogen gas. The impact on the geological repository's long-term safety could be attributed to a certain category based on the type of the effect caused by gas generation: (1) the mechanical effects on the repository and rock structures as a result of pressure build-up that may happen if dissolved gas cannot escape through the low-permeability host rock (disruption of the EBS, fracturing of the host rock), (2) the direct effect of gas on groundwater flow around the repository (induced groundwater flow, forcing the contaminated water from the repository or canisters, transport associated with bubbles) and (3) release of volatile radionuclides and toxic gas into the biosphere.

The recent EU-FP7 project FORGE was concerned with examining aspects of gas transport around and within potential geological repositories. The numerical modelling section of this project was focused on an understanding of the role of gas from a safety perspective. Investigations that were carried out in this section concentrated on the analysis of the generated hydrogen behaviour in a conceptual geological repository in clay formation taking into account a network of thin interfaces (i.e. contact surfaces within the EBS). Three benchmark exercises at three different spatial scales (at the cell scale (single disposal tunnel), at the module scale (100 interconnected disposal tunnels) and at the whole repository scale (10 modules and the access infrastructure)) were analysed (Wendling et al., 2013a,18b,19c). The Lithuanian Energy Institute (LEI) participated in the benchmark exercises at the cell and the module scales using the computer code TOUGH2 (USA). The LEI results of reference and some sensitivity cases defined by Wendling (2009) are presented in Justinavicius et al. (2012). Results of two additional independent analyses, not defined in the FORGE benchmark specification, are presented in this paper. First, the impact of the initial temperature in the disposal tunnel was investigated because the temperature in the geological repository depends on its depth, clay type and could vary significantly in different countries (Dalage et al., 2010). Two isothermal analyses assuming a lower (10°C) and a higher (30°C) initial temperature in the disposal tunnel compared to a reference case (20°C) in the FORGE benchmark were performed. Second, the impact of the tortuosity model was analysed because the saturation dependence on tortuosity is not well known at present (Ghanbarian et al., 2013). The numerical code TOUGH2 has three alternative formulations for a tortuosity model in a porous medium: Millington-Quirk, relative permeability-based and constant-tortuosity. Because the constant-tortuosity model is very conservative and unrealistic, only the impact on gas migration of the first two models was analysed.

In the generic repository concept (Fig. 1), HLW canisters are placed into the horizontal disposal cells situated in a clay formation at a depth of ~500 m. Shaft, main and access tunnels are used for HLW transportation. The geological repository is divided into modules on both sides of the main tunnel. Each module has 100 interconnected disposal cells that are symmetrically arranged on both sides of the access tunnel (drift).

A detailed description of a single disposal cell is presented in Fig. 2. HLW canisters (2) are transported by shafts down to the main tunnel at 500 m below ground level, from which they are transported by access drifts (7) to the disposal cell situated in the clay formation (1). The HLW canisters are placed in the cells and finally they are sealed with a bentonite plug (5). Because the contact between waste canisters and the engineered disturbed zone (EDZ) (6) as well as between the bentonite plug and the EDZ is not perfect, it is assumed that these materials are separated by thin centimetre-thick layers of a porous medium. These layers are called interfaces (3,4) and have different characteristics (Table 1). During the operational phase of the repository, the shaft, the main tunnel and the drifts are ventilated. The undisturbed clay and the EDZ are depressurized and desaturated. After the emplacement of the HLW canisters, the access drift, the main tunnel and the shaft are backfilled with the mixture of crushed clay and compacted bentonite. Then resaturation of the clay and backfill materials starts. At the same time, corrosion of the steel canisters causes hydrogen to be generated in the repository.

Hydrogen gas transport modelling was performed with the original version of the TOUGH2 code (Pruess et al., 1999) using PETRASIM (Thunderhead Engineering, 2010) as the graphical interface. TOUGH2 is a numerical simulator for non-isothermal flows of multicomponent, multiphase fluids in one, two and three-dimensional porous and fractured mediums. The original TOUGH2 code has several fluid (e.g. water-air, water-CO2, water-hydrogen, water-air-radionuclides) property or EOS (equation-of-state) modules where thermophysical properties of different fluid mixtures are described. In this study the EOS5 (water-hydrogen) module was used to analyse the behaviour of the groundwater system in which hydrogen gas release occurs. Groundwater and gas flow in the repository involves two-phase flow process. The mathematical model of gas migration in TOUGH2 was implemented taking into account hydrogen dissolution in groundwater, advective and diffusive flows.

The amount of gas that can dissolve in a liquid follows Henry's law, which states that at a constant temperature, the mass of the gas that dissolves in a liquid is directly proportional to the absolute pressure of the gas:
where, c is concentration of the solute (mol/m3), Pp is partial pressure of the solute (Pa) and KH is Henry's constant [mol/(m3·Pa)], which depends on the solute, the solvent and the temperature. In TOUGH2, Henry's constant for hydrogen solubility in water is a linear function of temperature up to 25°C. For higher temperatures a constant value for 25°C is assumed:
Advective flow density of gas and liquid is governed by a multiphase extension of Darcy's law:
where, subscript β denotes the fluid phase (liquid or gas), F.βA is advective flow density of fluid (kg/(m2·s)), k is intrinsic permeability (m2), kr is relative permeability for phase β (–), ρβ is fluid density (kg/m3), μβ is coefficient of dynamic viscosity of the fluid (Pa·s), ▿Pβ is pressure gradient of the fluid (Pa/m) and g is acceleration of gravity (m/s2).
Diffusion process in the gas and liquid phases is described by Fick's Law and can be written as:
where, F.βD is diffusive flow density of fluid (kg/(m2·s)), ϕ is effective porosity of porous medium (–), τ0τβ is tortuosity (m/m), Dβ is molecular diffusion coefficient (m2/s) and ▾Xβ is gradient of fluid mass fraction (1/m). Tortuosity represents the resistance to diffusion due to the influence of pore structure. Tortuosity is essentially a multiplier on the diffusion coefficient with a range 0–1. Smaller values indicate a more tortuous and longer route, while a value of 1 specifies no tortuosity and that the diffusive transport is controlled only by the mass fraction gradient and a diffusion coefficient for a particular phase. Tortuosity has a porous-medium-dependent factor τ0 and a coefficient that depends on phase saturation τβ. Two alternative models for tortuosity are available in TOUGH2:
τ0τβ={τ0krβrelative permeability (RP) modelφ13Sβ103Millington - Quirk (MQ) model}
where, Sβ is saturation of the fluid (—).
A molecular diffusion coefficient of dissolved hydrogen in water depends on temperature:
where T is temperature (K) and μwater(T) is the dynamic viscosity of water (Pa·s) at a particular temperature.

TOUGH2 solves mass balance equations by the integrated finite difference method. The equations are set up by combining mass balance, advective and diffusive fluxes, and solved by supplementing a number of equations which express interrelationships of physical processes, variables and parameters. The Van Genuchten-Mualem model (Mualem, 1976) for relative permeability and the Van Genuchten (1980) model for capillary pressure were selected to express dependence on effective saturation of a porous medium.

A two-dimensional disposal cell model in a cylindrical r-z system was developed and analysed. As the scale is quite small (few tens of metres), it was possible to represent fine geometric features and especially the interfaces. The modelling domain was meshed with 4836 rectangular grid elements and was refined in engineered materials because of the higher gradients of the physical variables and finer result resolution required. Both interfaces were represented with one layer of grid elements.

Because the HLW canisters are made of a material impermeable to liquid and gas flow, they were not explicitly represented in the model. The hydrogen gas source term (100 mol/year/disposal cell) was represented by a simple step function lasting 10,000 years. Hydrogen was injected directly into the interface between the HLW canisters and the EDZ. The constant pressure (5 MPa) condition at the outer radial boundary and time-varying boundary conditions (gas pressure and gas saturation) at the outer boundary of the access drift were assumed in the model. All other boundary conditions were designated as no-flow as indicated in Fig. 3.

Gas transport analysis was considered to be isothermal, i.e. residual heat removal from the canisters was neglected. Initial gas saturation was set to be 0.3 in the bentonite plug and the access drift and 0.95 in both interfaces. The remaining part of the modelling domain was assumed to be fully water saturated. The initial groundwater pressure was set to be 5 MPa in all water-saturated parts of the model (according to a continuous formulation in TOUGH2, the gas and water pressures are the same if only dissolved gas appears). Gas pressure in other parts of the model was set to be 0.1 MPa.

Analysed cases

Four independent analyses were performed in order to obtain information on the disposal system behaviour using different tortuosity models and different initial temperatures in the disposal tunnel. Descriptions of these cases are presented in Table 2.

The key data allowing evaluation of gas migration regularities in the geological repository are: gas saturation in porous materials, peak pressure in the analysed system and gas release (in gaseous and dissolved state) from the repository (Towler et al., 2012).

Reference case

Modelling results revealed that significant levels of gas saturation were reached in both interfaces, while in the EDZ the free gas phase only exists for particular periods of time and never exceeds 0.12. Only dissolved hydrogen was detected in the undisturbed clay formation. Full resaturation of the system was reached ~20,000 years after the repository closure. The peak pressure of ~5.7 MPa was observed in the interfaces and the EDZ. Gas-induced peak over-pressure was fixed at the end of the corrosion process and exceeded the hydrostatic pressure in the disposal cell by ~15%. However, such pressure was much lower than the lithostatic pressure (~10 MPa) in the selected geological repository concept; therefore mechanical effects on the engineered and natural clay structures of the repository caused by the pressure build-up were not assumed.

The graph of gaseous hydrogen flux through different surfaces in a single disposal cell is presented in Fig. 4. A positive flow rate indicates H2 flow towards the access drift.

The results showed that hydrogen was transported from the disposal cell towards the access drift during the period of gas generation (for up to 10,000 years). Most of the generated hydrogen (~85%) flows towards the access drift in gaseous form and the major part of this hydrogen (~81%) flows through the very thick waste interface (surface S-int1). The remaining part (~4%) of the hydrogen flows through the EDZ (surface S-EDZ). As the hydrogen accesses the plug interface, with a much lower permeability (equal to permeability of the EDZ), part of it (~49%) flows through the plug interface (surface S-int2) and the remaining part (~24%) flows towards the drift through the EDZ. Overall, ~73% of the generated hydrogen is transported from the disposal cell towards the access drift in a gaseous form, while the remaining part dissolves in porewater and diffuses in the undisturbed clay formation. Gaseous hydrogen that enters the access drift (surface S-drift) leaves it instantly through the lateral side of the drift boundary (surface S-drift(out)).

Temperature effect

Three isothermal analyses assuming different initial temperatures in the disposal tunnel were performed. Equations 2 and 6 indicate the temperature effect on hydrogen solubility and a molecular diffusion coefficient of dissolved hydrogen. Values of these parameters are summarized in Table 3.

Results obtained in cases of a lower (S1) and higher (S2) initial temperature were compared to the results of the reference case and are summarized in Table 4. As can be seen, the temperature variation has minor effect on the peak pressure in the disposal tunnel, but it has an impact on hydrogen flux through different surfaces.

In the case of a temperature of 10°C in the analysed system, the total mass flux of gaseous hydrogen toward the access drift (surface S-drift in Fig. 4) increased by 5.5% and total mass flux of dissolved hydrogen towards the undisturbed clay formation (surface S-clay in Fig. 4) reduced by 15.4% compared to the reference case. This was mainly because of a decreased molecular diffusion coefficient of dissolved H2. In the 30°C temperature case, the molecular diffusion coefficient of dissolved H2 increased and the total mass flux of gaseous hydrogen toward the access drift decreased by 7.5% and the total mass flux of dissolved hydrogen towards the undisturbed clay formation increased by 43.2% compared to the reference case. In total, increasing the initial temperature in the disposal tunnel from 10°C to 30°C gives ~9% lower mass transport of gaseous hydrogen from the disposal tunnel towards the access drift and higher mass transport of dissolved hydrogen to the undisturbed clay formation compared to the total mass of the generated hydrogen. The temperature influence on the molecular diffusion coefficient of dissolved hydrogen is more important compared to the temperature influence on the solubility of the hydrogen in all analysed cases.

Tortuosity effect

As indicated above, two analyses using relative permeability (RP) and Millington-Quirk (MQ) models for tortuosity (equation 5) were performed. It can be seen in Fig. 5 that dependence of tortuosity on liquid saturation for materials (in which diffusion process is relevant) is different using both models. Model MQ gives higher tortuosity values in an (almost) fully saturated porous medium and also increased diffusion in these materials compared to the RP model.

The MQ model for tortuosity had a significant impact on the results compared to the RP model. The main differences were determined in mass flux of the hydrogen through the analysed surfaces. As presented in Fig. 6, a much higher (35%) mass flux of gaseous hydrogen toward the access drift was obtained in the RP model. In the case of dissolved hydrogen flow toward undisturbed clay formation, the situation was opposite; much higher mass flux was determined in the MQ model. The peak pressure in this case also reduced by 1.8% compared to the RP model.

After assessing the influence of temperature variation in the disposal tunnel and different tortuosity models on gas migration using the computer code TOUGH2, the following conclusions can be made:

  1. If the Millington-Quirk tortuosity model is applied, total mass of gaseous hydrogen transported from the disposal tunnel towards the access drift decreases by ~35% because of the more effective diffusion of dissolved hydrogen compared to the relative permeability model.

  2. Increasing the initial temperature in the disposal tunnel from 10°C to 30°C gives ~9% lower mass transport of gaseous hydrogen from the disposal tunnel towards the access drift and higher diffusion of dissolved hydrogen to the undisturbed clay formation compared to total mass of the generated hydrogen.

  3. Temperature variations in the disposal tunnel and different tortuosity models have a minor impact on peak pressure in the analysed system.

This research was performed within the frame of EC Project on ‘Fate of Repository Gases (FORGE)’ under the 7th EURATOM FP (FP7-230357) and co-funded by Lithuanian Agency for Science, Innovation and Technology (MITA).

The publication of this research has been funded by the European Union's European Atomic Energy Community's (Euratom) Seventh Framework programme FP7 (2007–2013) under grant agreements n°249396, SecIGD, and n°323260, SecIGD2.
P1Freely available online through the publisher-supported open access option.