## Abstract

This paper presents the results of a benchmark study that compares a number of numerical models applied to a specific problem in the context of hydrogen flow and transport in a nuclear waste repository. The processes modelled are two-phase (water and hydrogen) immiscible compressible two-component transient flow in a heterogeneous porous medium under isothermal conditions. The three-dimensional (3D) model represents a module of a repository for high-level waste in a clay host rock. An upscaling technique and a vertex-centred finite-volume method are employed to yield very accurate solutions. Since the full range of results required in the benchmark is too large to be displayed in this paper, we focus on the evolution of the pressures, the saturations, the fluxes and the comparison of the numerical results with the other participants. A homemade C++ upscaling code and the parallel multiphase flow simulator DuMu^{X} have been adopted for this study.

The long-term safety of the disposal of nuclear waste is an important issue in all countries with a significant nuclear programme. Repositories for the disposal of high-level and long-lived radioactive waste generally rely on a multibarrier system to isolate the waste from the biosphere. The multibarrier system typically comprises the natural geological barrier provided by the repository host rock and its surroundings, and an engineered barrier system (i.e. engineered materials placed within a repository, including the waste form, waste canisters, buffer materials, backfill and seals).

In the frame of designing nuclear waste geological repositories, a problem of possible two-phase flow of water and gas appears (for more details see, for instance, OECD/NEA 2006; Norris 2010). Multiple recent studies have established that in such installations important amounts of gases, mainly hydrogen, are expected to be produced, in particular due to the corrosion of metallic components used in the repository design. The creation and transport of a gas phase is an issue of concern with regard to the capability of the engineered and natural barriers to evacuate the gas phase and avoid overpressure, thus preventing mechanical damage. It has become necessary to carefully evaluate these issues while assessing the performance of a geological repository (see, e.g., Croisé *et al.* 2011; Senger *et al.* 2011; Enssle *et al.* 2012 and the references therein).

Significant effort has been, and continues to be, expended in numerous national and international programmes in an attempt to understand the potential impact of gas migration, and of two-phase gas–water processes, on the performance of underground radioactive waste repositories, and to provide modelling tools or approaches that will allow these impacts to be assessed. Verification of numerical models for two-phase immiscible compressible flow in porous media by the means of appropriate benchmark problems is a very important step in developing and using these models.

In 2007, the French Agency for the Management of Radioactive Waste (ANDRA 2005) and the French research group MoMaS (http://www.gdrmomas.org/) launched a two-phase flow-simulation benchmark, named Couplex-Gaz, to improve the simulation of the migration of hydrogen produced by the corrosion of nuclear waste packages in underground storage. This is a system of two-phase (water–hydrogen) flow in a porous medium. This benchmark generated some interest, and engineers encountered difficulties in the handling numerical simulations for this model. One of the purposes of this test case is to model the resaturation of a disposal cell for intermediate-level long-lived waste (known as Type-B waste) from as soon as the cell is closed until the end of the gas-production period (for more details see, for instance, Zhang *et al.* 2011).

Since 2010, three benchmarks have been studied within the framework of the European Project FORGE: Fate Of Repository Gases WP1.2 (http://www.bgs.ac.uk/forge/). To be able to understand water and gas flow at the repository scale, it is necessary to understand the water and gas flow associated with subrepository ‘cell’-scale components (e.g. on the scale of a vault or tunnel). Information is needed at an appropriate level of detail in order to upscale from the cell level to the full repository level, perhaps via intermediate modules of scale (e.g. several tens of cells). For this reason, the first benchmark was defined at cell level (Wendling 2013*b*). The type of cell considered is a somehow generic high-level waste cell in a clay host rock. As the scale is quite small, it is possible to represent fine geometric features, and, especially, the interface between the canister and clay rock, as well as interface between the bentonite seal and the clay rock (considered as a centimetre-thick region). As permeability is increased in the excavation damaged zone (EDZ), this zone is also represented. The aim of this test case is to better understand the mechanisms of the gas migration at the cell scale, and, in particular, to analyse the effect of the presence of different material and interfaces on such mechanisms. The results can be found in Amaziane *et al.* (2013) and Wendling (2013*b*).

The second benchmark was realized at a module scale (several tens of cells) (Wendling 2013*a*). This is the intermediate scale towards a full repository representation. Its main objective is to understand how gas is moving from a cell towards a drift and, finally, towards a drift plug. (Is the convection still the main process at this scale? Which part of the gas generated inside the cells is moving across the drift plugs? What is the characteristic time for this transfer? What pressure can be achieved? …). The mathematical model in this test case involves two-phase immiscible compressible (liquid and gas phases), transient flow processes in a porous medium under isothermal conditions, including gas dissolution and water evaporation. For gas flow, a mixture of gases is involved, with hydrogen being the chief constituent. In this study, all gases produced by the corrosion phenomena of metallic components are treated as hydrogen. Hydrogen dissolution in the porewater follows Henry's law; gas and water flow are governed by Darcy–Muskat's law. The main difficulties related to this benchmark are the size and heterogeneities of the three-dimensional (3D) module and the complex non-linear processes involved.

For the performance assessment, a good understanding of the behaviour of the system over a large timscale of several hundreds of thousands of years is compulsory. Furthermore, simulating a number of different geological realizations or scenarios is required to reduce uncertainty in performance assessment. Thousands of such runs may be required to cover the range of parameter variation.

Theoretically, the best way is to simulate every single detail of the module. It is not computationally feasible (or desirable) to perform these simulations on the fine-grid model because of the huge size and the complex structures of the calculation domain, and the large simulation period, as well as the limit of the computational power. Therefore, an upscaling method has been used in the benchmark to represent the domain on the coarse grid and to give simulation results comparable to the results obtained on the fine grid, in order to satisfy all of the requirements of the test-case simulation.

The upscaling or homogenization of multiphase flow through heterogeneous porous media has been an interesting problem for many years and many methods have been developed. There is extensive literature on this subject. We will not attempt a literature review here but merely mention a few references. Here we restrict ourselves to the mathematical upscaling or homogenization method as described in Hornung (1997) for flow and transport in porous media. The homogenization method has the advantage of mathematical rigorousness and allows us to establish explicit connections between processes at various scales. A recent review of the methods developed for incompressible immiscible two-phase flow in porous media can be viewed in Farmer (2002) and Durlofsky (2005). Let us also mention that few homogenization results were obtained in the case of fields with different rock types: porosity, absolute permeability, relative permeabilities and capillary pressure curves being different in each type of porous media (see, e.g., Bourgeat & Hidani 1995; Bourgeat *et al.* 1996; Bourgeat & Jurak 2010; Henning *et al.* 2013). We refer, for instance, to Bourgeat *et al.* (2004), Bourgeat & Marusic-Paloka (2005), Gipouloux & Smaï (2008), Bourgeat & Piatniski (2010) and Bourgeat *et al.* (2010*a*, *b*) for more information on the homogenization of one-phase flow within the framework of the geological disposal of radioactive waste.

However, the situation is quite different for immiscible compressible, such as water–gas, flow in porous media, where, only recently, a few results have been obtained for the the case of a domain made up of several zones with different characteristics (Amaziane *et al.* 2010*a*, 2014). A rigorous mathematical derivation, by means of the two-scale convergence, of the upscaled model used in the numerical simulation for the FORGE module-scale benchmark is presented in Amaziane *et al.* (2010*a*, 2014). Upscaling is used to transform a fine-grid model into a more practical, coarser one. Effective parameters are assigned to a group of neighbouring fine-grid blocks. This method allows the determination of these effective parameters from the knowledge of the geometrical structure of a representative elementary volume (REV) and its heterogeneities by numerically solving appropriate local problems.

The purpose of this paper is to describe the results obtained for the FORGE module-scale benchmark with the parallel multiphase flow simulator DuMu^{X} (http://www.dumux.org: Flemisch *et al.* 2011) and the upscaling method that has been used. The rest of the paper is organized as follows. First, we give a short description of the benchmark. Secondly, we describe the mathematical and physical model used in this study. Thirdly, we present the mathematical upscaling method to represent all model data on a coarser simulation grid. We have developed our C++ homemade code using a finite-element method to compute effective porosity, absolute permeability, capillary pressure, relative permeability and tortuosity curves obtained by the upscaling procedure. Two upscaled models are derived to simulate the benchmark. We exhibit the numerical results for the module-scale benchmark. Lastly, some concluding remarks are forwarded.

## Description of the benchmark

The complete definition of the benchmark can be found in Wendling (2013*a*). Here we only present a brief description of the configuration of the benchmark at the module scale. Figure 1 represents a horizontal cross-section of the whole repository model, which is composed of 10 modules connected by the main drift going to the well that connects the repository to the surface.

Each module contains 100 waste canisters, 50 on each side of the access drift that connects them to the main drift. One of the modules is shown in Figure 2. Each canister in the module is separated from the access drift by a bentonite plug (see Fig. 3). The bentonite plugs are also placed in the main drift in order to separate one module from another and, finally, there is a bentonite plug in the well that separates the whole repository from the surface.

All tunnels that make the repository are surrounded by a layer of the EDZ, which is a fractured medium and thus more permeable than the surrounding geological medium. The EDZ is found around the canisters, the access and the main drifts (Fig. 3), and is only not present around the main drift bentonite plugs that are in direct contact with the geological medium. The main and access drifts are filled with the material called the backfill in Figure 3. The contact between the waste canister and the EDZ is not perfect, and a thin space (of 1 cm) exists there called the interface, which is presented in the model as a very permeable porous medium with very low capillary pressure curve (compared to other material curves). This interface is present at all contacts of the bentonite plugs and the EDZ (or the geological medium for the main drift bentonite plugs), and at the contacts of the backfill and the EDZ. Its presence can be seen in Figure 3 where the vertical cross-section B–B' is depicted.

During the first 10 000 years, all of the canisters generate hydrogen gas by metal corrosion at a rate of 100 mol per year, per canister, which then migrates through the repository. The aim of the module-scale benchmark is to simulate the gas migration during the first 100 000 years after the closure of the repository. Since in this benchmark the whole repository is not considered, an appropriate boundary condition is given at the ends of the main drift where the chosen module is connected to the other modules. A complete description of the boundary conditions is given in Wendling (2013*a*).

The initial conditions are different in the different materials that make up the repository: interfaces, backfills, bentonite plugs, the EDZ and the geological medium. Table 1 exhibits the initial conditions for the phase pressures, *P*_{g} (gas) and *P*_{w} (liquid), and the liquid saturation, *S*_{w}, in different materials.

The geological medium is initially water saturated with imposed hydrostatic pressure. The whole computational domain is a box 150 m high with the repository in the middle of it (at a height of 75 m). The hydrostatic water pressure at the upper surface is 3 MPa and 5 MPa at the lower surface, and these pressures are maintained during the simulations. On the lateral sides of the box, no-flow conditions are imposed. We also note that the canisters do not make up part of the computational domain. Only their surfaces serve as the gas sources at the given rate.

Tables 2 and 3 exhibit the physical properties of the different materials; there are three different type of interfaces. The capillary pressures and relative permeability functions are given by van Genuchten–Mualem model (van Genuchten 1980). The residual water and gas saturations are equal to zero.

The main difficulty of the benchmark is its highly heterogeneous structure, shown in Tables 2 and 3, and the presence of structures such as the interfaces and the EDZ, which have very small dimensions compared to the dimensions of the computational domain.

## Mathematical model

The mathematical model described in the FORGE module-scale benchmark (Wendling 2013*a*) is an isothermal two-phase flow through heterogeneous porous media under local thermodynamic equilibrium, without chemical reactions or biological decomposition. Only two chemical components are present in the system (H_{2}O and H_{2}) and dissolution/evaporation processes are taken into consideration (see Table 4 for the notation used in the following equations).

The dependence between water and gas saturation is expressed by:
(1)
where *V*_{w} is the liquid volume [m^{3}], *V*_{g} is the gas volume [m^{3}] and *V*_{p} is the pore volume [m^{3}].

The phase pressures are connected by the capillary pressure law: where [Pa] is the capillary pressure curve.

The van Genuchten–Mualem model is used to express the capillary pressure function and the relative permeabilities. We have:
(2)
(3)
(4)
where *S*_{we} is the effective water saturation:
(5)
*S*_{wr} and *S*_{gr} are the residual water and gas saturations, *P*_{r} is the reference pressure [Pa] and *n*, *m*=1−1/*n* are coefficients for van Genuchten's law.

Liquid and gas flow are governed by Darcy–Muskat's law:

The mass conservation law for each phase writes as follows: (6)

The mass fraction of the component *c* in phase α is expressed by:
(7)
where *m*_{α}^{c} is the mass of component *c* in phase α [kg] and *m*_{α} is the total mass of phase α [kg]. Mass fractions satisfy:

The mass conservation law for each component *c*∈{H_{2}O, H_{2}} is given by:
(8)
where the sum goes over both phases. The diffusion process in the gas and the liquid phase are described by Fick's law, and the diffusive fluxes [kg m^{2} s^{−1}] in both phases α∈{w, g} and for both components *c*∈{H_{2}O, H_{2}} can be written as:
(9)

The mixture is supposed ideal and the diffusion coefficients are given by:
(10)
(11)
where *T*=*T*_{o}=293 K, *P*_{o}=10^{5} Pa and *D*_{o}=9.5×10^{−5} m^{2} s^{−1}, and we have taken μ_{water}(*T*)=10^{−5} Pa s.

The viscosity of the gas mixture (water vapour+hydrogen) can be estimated by a classical Wilke approximation or by a simplified formula as follows: (12) where and

The water density is given by:
(13)
where ρ_{atm} [kg m^{−3}] and *P*_{atm} [Pa] are atmospheric density and pressure, respectively, and *S*_{s} is the specific storage coefficient [Pa^{−1}].

The quantity of the hydrogen dissolved in the water is given by Henry's law, written in the following form:
(14)
where is the maximum molar concentration of the hydrogen in the water [mol m^{−3}], *H*_{H2} (*T* = 7.6 × 10^{−6} [mol Pa^{−1} m^{−3}] is the constant of the Henry's law at *T*=293 K, is the partial pressure of the hydrogen in the gas phase [Pa] and is the molar mass of the hydrogen [kg mol^{−1}].

The relationship between partial pressure of each gas present in the total gas phase and total gas pressure is given by Dalton's law, which states that for a binary mixture (gaseous H_{2} and water vapour):
(15)

Each of the gases is supposed perfect:
(16)
where *P*_{g}^{c} is the partial pressure of the component *c* in the gas phase [Pa], *M _{c}* is the molar mass of component

*c*[kg mol

^{−1}], ρ

_{g}

^{c}is the mass concentration of component

*c*in the gas phase [kg m

^{−3}],

*R*=8.314 J mol

^{−1}K

^{−1}is the constant of the perfect gas and

*T*is the temperature [K].

The saturation pressure for water vapour, *P*_{sat} [Pa], is only dependent on temperature and can be expressed as:
(17)
where *T*_{c} is the temperature [°C].

Kelvin's law gives a relationship between the saturation pressure for water vapour, the effective pressure for water vapour and the capillary pressure: (18)

**Remark** Given the actual hydrogen concentration in water, , it is possible to define a ‘pseudo-pressure’ of the dissolved hydrogen, , by inversion of Henry's law:
(19)

## Upscaling

A detailed 3D modelling of the FORGE module-scale benchmark would require a tremendous computational effort, even when using high-performance simulator codes. Because such a detailed model can be highly CPU-time consuming, there exist various attempts that aim at the optimum compromise between the quality of the simulation and the required amount of CPU time. Upscaling refers to the techniques used to transform a fine-grid model into a more practical, coarser one.

In an upscaled model, each coarse-grid block is composed of a number of fine-grid blocks, all having different physical properties. The upscaling replaces these heterogeneous properties of the porous material within the coarse-grid block with the equivalent homogeneous ones, which will be called the effective properties. The two grids, fine and coarse, represent distinct spatial scales: the fine one is used to represent the heterogeneities of the model; and the coarser one is used for representation of the gas-transport process.

The aim of the upscaling procedure is to incorporate the fine-scale information given on the fine grid into the coarse-scale medium properties in order to represent well the gas-transport process on the coarse-grid scale. The fine-scale details of the gas-transport process will be lost but the coarse-scale details must be adequately represented. The lifting of the fine-scale variations of the medium properties to the coarse scale is done by a simplified representation of the fine-scale transport process on each coarse-grid block. This representation is called the local cell problem, and has to be solved for each grid block in the coarse mesh. From these solutions, one calculates the mean fluxes by an averaging procedure and these fluxes are then used in the calculation of the effective properties of the medium on the coarse grid. The effective medium properties are homogeneous properties given to each coarse-grid block, which are defined by the requirement that they must produce the same mean fluxes as the local problems. These procedures are detailed in the next subsection. In that sense, an equivalent homogeneous medium replaces the heterogeneous one in each coarse-grid block. Simplifying assumptions used in the formulation of the local cell problems may be different but, in the case of capillary dominated two-phase flow, the capillary equilibrium hypothesis is very efficient. It produces local problems that can be easily solved, and the coarse-scale transport problem is governed by the same type of differential equations as the fine-scale one.

The mathematical upscaling or homogenization allows these effective parameters to be determined from the knowledge of the geometrical structure of a REV and its heterogeneities by solving appropriate local problems. The local problem equation is numerically solved and an effective parameter value is calculated based on the obtained solution. To solve this local equation we have to assume a set of boundary conditions. The most general case is to assume linear or periodic boundary conditions. Another alternative is the method known as the Pressure Solver Method where constant-pressure no-flow boundary conditions are used. However, in Amaziane & Koebbe (2006), it is shown that in many cases the differences between the various methods were small. Linear boundary conditions are convenient, particularly in complicated geometries with unstructured grids.

For the numerical simulation of the FORGE module-scale benchmark, to circumvent the heterogeneities and the difficulty in taking into account the thin interface surrounding the waste canisters and the bentonite plugs, we have used the above-described mathematical upscaling. In this subsection, we briefly recall its main properties and its application to the FORGE module-scale benchmark. More details about the method can be found in Amaziane *et al.* (2001, 2006), Amaziane & Koebbe (2006) and Bourgeat & Jurak (2010).

### Mathematical upscaling

In this subsection, we will give a brief description of how the effective properties have been obtained by the mathematical upscaling or homogenization method. More precisely, the upscaled model is obtained using the two-scale asymptotic expansions as described in many textbooks (see, e.g., Hornung 1997). The presentation of the whole of the calculations is out of the scope of this paper but can be found in Amaziane *et al.* (2010*b*). Here we will only present the basic idea of the method and the local problems to be solved in order to compute the effective properties.

The method is based on the following assumptions: the periodic porous medium has two characteristic lengths: the macroscopic length, *L*, and the microscopic length, *l*. Their ratio is the scale parameter, which should be small, ε=*l*/*L*≪1. There are two spatial scales: *x* being macroscopic and *y*=*x*/ε being microscopic.

Denoting the sequence (indexed by ε) of solutions of the considered partial differential equation with periodically oscillating coefficients by *u*_{ε}, a two-scale asymptotic expansion is an ansatz of the form:
(20)
where each function *u*_{i}(*x*, *x*/ε) in this series depends on two variables, *x*, the macroscopic (or slow) variable, and *y*, the microscopic (or fast) variable, and is *Y*-periodic in *y* (*Y* is the unit period).

Inserting the ansatz (equation 20) into the equation satisfied by *u*_{ε} and identifying powers of ε leads to a cascade of equations for each term *u*_{i}(*x*, *x*/ε). In general, averaging with respect to *y* yields the homogenized equation for *u*_{0} and the local problems to be solved in order to compute the effective parameters. The asymptotic analysis of the problem (i.e. the limit as ε tends to zero) yields a rigorous definition of the effective parameters. The results of this procedure can easily be extended to a non-periodic case by replacing the unit period *Y* by the representative elementary volume.

The complete benchmark model can be represented only on a very fine grid. By means of mathematical upscaling, we represent all model data on a coarser simulation grid, which cannot directly represent the data small-scale variations such as the interfaces and the EDZ. We will refer to upscaled data assigned to the cells of the coarser grid as the *effective* data. Our goal is to calculate the effective values for the porosity Φ(*x*), the absolute permeability *K*(*x*), the capillary pressure *P*_{c}(*x*, *S*), the relative permeability *kr*_{α}(*x*, *S*) and the tortuosity coefficient τ(*x*). The mathematical upscaling is a steady-state upscaling technique based on the local capillary equilibrium hypothesis.

Let us explain the calculation of the effective absolute permeability. In a chosen REV, denoted *V* (i.e. a cell of the coarse grid in which the absolute permeability *K*(*x*) is heterogeneous), this heterogeneous value will be replaced by an homogeneous effective absolute permeability, *K**, in the following way: by imposing suitable boundary conditions we generate heterogeneous fluxes in all three space dimensions. Then we take the mean value of these fluxes and we search for a constant effective absolute permeability, *K**, that will generate these mean fluxes when the same boundary conditions are imposed.

More precisely, we have to solve in *V* the following steady-state Darcy's flow problems for *i*=1, 2, 3:
(21)
where *x _{i}* is the

*i*-th coordinate. Then we calculate the mean flux: and we define the ‘effective absolute permeability’

*K** by

*K**, where is the

*i*-th vector of the canonical basis in

*IR*

^{3}. Now it is not difficult to see that with the same boundary conditions,

*P*=

_{i}*x*, the permeability

_{i}*K** will produce the Darcy flux

*K**, which is by our definition equal to .

Solving the local problems (equation 21) gives the following expression for the coefficients of the tensor *K**:
(22)

Assume that *K* is a symmetric positive definite tensor, then it is easy to see from equation (22) that *K** is symmetric and positive definite. But, in general, even with each porous medium being isotropic, we may have an effective tensor that is significantly anisotropic (i.e. *K** is typically a full tensor quantity). The procedure for computing *K** requires the solution of the fine-scale flow (equation 21) over the target coarse region *V*. A numerical code (Amaziane & Koebbe 2006), using a finite-element method, has been performed to solve the local problem (equation 21) and to compute the effective absolute permeability tensor for any given REV.

The *effective porosity*, designated Φ*, is computed such that pore volume is exactly conserved between the fine and coarse scales. Specifically, Φ* is computed by:

The idea of conservation of the mean fluxes for calculating the effective absolute permeability can be extended to upscale the relative permeability and the tortuosity.

Let the coarse grid cell *V* be composed of *n* different materials *V*_{1}, … , *V*_{n}, and to the *j*-th material we associate the capillary pressure function *P*_{c}^{j} (*S*), and the relative permeability functions *k*_{rw}^{j}(*S*) and *k*_{rg}^{j}(*S*). To upscale the relative permeabilities, we need a representative distribution of the liquid saturation in the coarse grid cell *V*. That distribution is given by the capillary equilibrium hypothesis that imposes that the capillary pressure must be constant in the whole volume *V*.

This saturation can be represented as:
(23)
where is the characteristic function of *V*_{j}, and *S*_{j} is given by the capillary pressure *u* which is imposed in *V*. Therefore, the values *S*_{j} are computed from the following equations:
(24)

Then we define the *effective saturation*, *S**, corresponding to the capillary pressure level *u*, by:
(25)

This definition of *S** preserves the quantity of water in the cell *V*. The *effective capillary pressure* function, , can be now naturally defined as . This function is calculated in a form of a table in the following way: for each chosen capillary pressure value *u*, we calculate *S** using equations (23)–(25). Then, the pair (*S**, *u*) represent a point on the graph of the function . Note that in equation (24), we do not solve non-linear equations, but simply, for given *u*, we use the inverse of capillary pressure functions to calculate the saturation values.

For a given effective saturation, *S**, and corresponding local saturation distribution, *S*^{0}(*x*), we can solve the following local problems for *i*=1, 2, 3:
calculate average fluxes:
and define the effective relative permeability functions (*k*_{rα}(*S**))_{t} in the space direction *i*, for *i*=1, 2, 3, by the formula:

Again, by this definition of the effective relative permeability, will produce the mean flux in the *i*-th direction generated by the heterogeneous permeability *K*(*x*)*k*_{rα}(*x*, *S*^{0}(*x*)).

Similarly, the effective tortuosity function , in the space direction *i*, is defined by preservation of the mean diffusive flux −Φ(*x*) and, therefore, will depend on the phase saturation. For each space direction *i*, we solve:
and we calculate the mean fluxes:

Then we define the effective tortuosity in the phase α=w, g and the space direction *i*=1, 2, 3, by:

Since the above procedure gives directional effective relative permeabilities and tortuosities, in each coarse-grid block *V* we choose dominant direction and we replace all of three directional curves with the one in the dominant direction.

We have developed our C++ homemade code using a finite-element method to compute effective porosity, absolute permeability, capillary pressure, relative permeability and tortuosity curves obtained by the upscaling procedure described in this subsection. Some numerical results will be presented in the following subsection for the FORGE module-scale benchmark.

### Upscaled models for the benchmark

The FORGE module-scale benchmark is limited to a portion of a repository containing an access drift serving a number of cells. The calculation domain is 3D, and includes several materials such as argillite (geological medium), the EDZ (excavation damaged zone), bentonite (cells and drift plugs), backfill of the drifts and interfaces between different materials. One of the major difficulties in the numerical simulation of this test case is to take into account simultaneously all of the sources terms (generally located on the boundary of the storage cells) and to represent, at the same time, the transfer network constituted by the different drifts and cells.

In order to use our upscaling procedure, we group several heterogeneous fine-grid blocks into one coarse-grid block in which we replace the heterogeneous properties by homogeneous, effective ones. This grouping can be done in different ways, and different levels of homogeneity of the upscaled model can be achieved. Our aim is to eliminate the interfaces and the EDZ from the model because of their small size. We have designed two levels of upscaling, denoted UM1 and UM2 (Upscaled Model), of increasing homogeneity in order to verify, by comparison, the capacity of the upscaled models to simulate global behaviour of the module. Figure 4 represents a vertical cross-section of these upscaling strategies near the canister (bottom) and in the drifts (top).

In UM1, the canisters and their plugs, EDZ, and the interface are substituted by a homogeneous block. In UM2, the homogeneous block is enlarged by the surrounding geological medium: that is, the EDZ and the interface, the canister, the plug and the surrounding geological media are substituted by one homogeneous block. This is illustrated at the bottom of Figure 4. We also note that after this upscaling the gas source, originally placed on the boundary of the canister, has to be distributed onto the whole homogeneous block containing the canister.

In the main and access drifts, the interface, the backfill, and the EDZ are substituted by one homogeneous block for both upscaled models, UM1 and UM2. This is illustrated at the top of Figure 4.

Each upscaled model is characterized by effective porosity, absolute permeability, relative permeability, capillary pressure and tortuosity curves, calculated by the procedure described previously. In order to save space, we will show only these effective properties for one effective homogeneous material, the one obtained by upscaling the backfill, the EDZ and the interface in the access drift.

The porosity and the absolute permeability of these materials are presented in Tables 2 and 3. The effective values are given by: The effective relative permeability in the direction of the drift (dominant direction) and the capillary pressure curves are shown in Figure 5, where they are compared to the backfill curves.

## Numerical results

In this section, we present the numerical solution of the benchmark problem defined previously. The simulations conducted in this benchmark are computationally challenging because of the size and heterogeneities of the 3D module, and the complex non-linear processes involved. The problem can be formulated as a highly coupled degenerate system of partial differential equations. Note that the benchmark deals with an immiscible compressible two-phase two-component flow in a porous medium under highly contrasted capillary pressures that lead to a degeneracy in the system. The main difficulty in numerical simulation of this system consists of dealing with the possible appearance and disappearance of a fluid phase. This process strongly reduces the time step during the simulation; because there is a final large simulation time of 100 000 years, this makes the simulation time extremely long. We need to upscale the interfaces and the EDZ in the benchmark in order to avoid problems of flow resolution in very permeable interface that can lead to convergence problems, and can result in small time steps.

Simulations are performed with the simulator DuMu^{X} (http://www.dumux.org: Flemisch *et al.* 2011). DuMu^{X}, DUNE for Multi-{Phase, Component, Scale, Physics, …} flow and transport in porous media, is a free and open-source simulator for flow and transport processes in porous media, based on the Distributed and Unified Numerics Environment DUNE (http://www.dune-project.org/). DuMu^{X} is developed at the University of Stuttgart and is part of the OPM (Open Porous Media) initiative.

DuMu^{X} has several mathematical flow models implemented in different modules. We use a 2p2c module that implements two-phase flow with two components. A coupled fully implicit approach is adopted. For the spatial discretization, a vertex-centred finite-volume approach is used, while the time discretization is performed thanks to an implicit Euler scheme. An automatic control and adaptive time-stepping is based on the number of iterations required by the Newton method. Finally, the software Gmsh (Geuzaine & Remacle 2009) is used for the mesh generation

A simulation with the complete model would require a very fine grid to take into account the small-scale details, such as the interfaces, and the duration of the computation would be prohibitively long. Therefore, some upscaling strategy is mandatory in order to reduce the size of the problem; hence, we have implemented the mathematical upscaling described previously. We will present the numerical results obtained for the two upscaled models UM1 and UM2 defined above. The CPU time on eight processors for the UM2 is about 4 h, while the UM1 takes approximately 1 month on 24 processors.

The FORGE benchmark requests three types of output data: evolution of the pressures; saturations and mass fractions in prescribed points and along prescribed lines; and evolution of the mass flows through prescribed surfaces. Figure 6 presents the points, lines and surfaces where the data have been extracted. In the following sections of this paper, we will present just a selection of these results owing to space constraints.

### Results for UM1

The UM1 simulation is performed on a grid composed of 188 298 elements. This is presented in Figure 7 without the surrounding geological medium, which is not shown in order to make the structure of the repository visible.

In Figures 8 and 9 we present the time evolution of the fluxes through the surfaces F-C25 and F-Pd. F-C25 is a vertical plane cross-section of the waste container in the 25th row, passing through the contact of the container and the plug. F-Pd is a vertical plane cross-section of the one main drift plug, passing through the middle of the plug (see Fig. 6). Figures 8 and 9 exhibit, respectively, the dissolved and gaseous hydrogen fluxes for surfaces F-C25 and F-Pd. Each total flux is the sum of corresponding advective flux and diffusive flux. First, we note that the gaseous hydrogen total fluxes (Fig. 9) are three orders of magnitude larger than the dissolved hydrogen total fluxes (Fig. 8). We also note that, in the gaseous phase, the advective fluxes dominate the diffusive fluxes by several orders of magnitudes; in the liquid phase, however, they are comparable.

Figure 10 displays the evolution of the gas pressure in points defined by the benchmark (see Fig. 6 for their positions). The pressure maximum is about 7 MPa and it is achieved before the extinction of the gas source at 10 000 years.

Finally, Figures 11 and 12 represent, respectively, the gas pressure and the gas saturation on the line L-MD that passes along the middle of the main drift (see Fig. 6). Figure 12 shows that the main drift plugs are almost fully resaturated after 1000 years. Figure 11 shows pressurization of the main drift, which does not exceed 7 MPa, and the return to equilibrium, which is not completely achieved even after 100 000 years.

### Comparison between UM1 and UM2

The UM2 simulation is performed on a grid composed of 27 776 elements, as shown in Figure 13. Note that the blocks surrounding the access drift in Figure 13 represents a homogeneous material obtained by upscaling of the canisters, interface, the EDZ and the bentonite plugs with surrounding geological medium. Again, the geological medium surrounding the upscaled module structure is not shown in order to make the structure visible.

The differences between the UM1 and UM2 numerical models are in the grid sizes (188 298 v. 27 776 elements, respectively, which is a reduction by a factor of 7 from UM1 to UM2), and in the effective properties of different media: porosity, absolute permeability, relative permeability and capillary pressure curves, but also in the tortuosity and the initial conditions.

Figure 14 shows several comparisons between the UM1 and the UM2 simulations.

The mathematical upscaling replaces the heterogeneous medium by a homogeneous one with the same global behaviour. Therefore, it cannot be expected that local evolution of the system will be close to the UM1 and UM2 simulations in all aspects, but only that their global behaviour (e.g. field pressurization, gas preferential paths, resaturation periods) will follow the same pattern. We also expect that the UM2 simulation will locally give results that are more accurate as it retains more heterogeneities. As seen in Figure 14, the largest differences can be found in the saturation variable as, in differently upscaled models, the saturation represents the mean saturation calculated over different volumes. Since the initial data must also be upscaled in order to keep the initially given amount of gas, the saturations in the two simulations have different initial conditions. However, the time evolution follows the same pattern. The gas and water pressures at various points show, as expected, more local consistency. Both simulations have the maximum gas pressure of about 7 MPa. Also, both simulations show that the main drift plugs are approximately resaturated after about 1000 years (see Fig. 12). The maximum desaturation in the access drift is about 18% and in the geological medium below 2% in both simulations, which is not shown in the figures of this paper.

The results show that, locally, the differences between the two simulations can be large but global behaviour is similar enough to confirm the same conclusions about the response of the system to the migration of generated gas. Finally, comparison between the duration of the computations for models UM1 and UM2 shows the efficiency of the upscaling procedure (reduction from 1 month to 4 h) while retaining the same overall behaviour.

### Comparison with other participants

Figure 15 compares the gas pressure at the points C50-3 (top) and Pd-1 (bottom) between our simulations and simulations of four other participants of the benchmark: the Lithuanian Energy Institute (LEI), the Nuclear Waste Management Organization (NWMO), the Agence Nationale pour la gestion des Déchets RAdioactifs (ANDRA) and the Nuclear Decommissioning Authority (NDA). The results are in good accordance and show the reliability of our approach.

## Conclusion

We have presented numerical simulations of the FORGE module-scale benchmark. Our approach for the numerical modelling includes the use of a mathematical upscaling procedure to reduce the heterogeneity of the model, thus reducing its complexity. First, we obtained effective properties for two upscaled models. Then a vertex-centred finite-volume method was used within the framework of DuMu^{X}, a parallel multiphase flow simulator, to simulate the proposed 3D modelling of hydrogen gas migration in a radioactive waste repository at a module scale.

The simulation results have been compared to corresponding results from other teams and have shown a reasonable match. Our simulations show that the maximum pressure in the module will be about 7 MPa and that, at the places where the fluxes were calculated, the convection in the gaseous phase will be the main method of hydrogen transport. The transport of hydrogen dissolved in water is about two–three orders of magnitude less significant than the transport of gaseous hydrogen. Hydrogen is transported from the cells to the access and main drifts, which represent preferential paths for the hydrogen migration. The geological medium shows only a slight desaturation of less than 2%. The main drift bentonite blocks are almost fully resaturated after 1000 years.

The above results illustrate that the proposed mathematical upscaling combined to a finite-volume method is capable of tackling, in a robust and accurate fashion, various physical phenomena relevant to hydrogen flow and transport in a nuclear waste repository.

## Acknowledgments

The research leading to these results has received funding from the European Atomic Energy Community's Seventh Framework Program (FP7/2009-2013) under Grant Agreement No. 230357, the FORGE project. This work was partially supported by the GnR MoMaS (PACEN/CNRS, ANDRA, BRGM,CEA, EDF, IRSN) France, all of their support is gratefully acknowledged. We also wish to thank all of the teams who took part in the test case presented herein for their very active participation: L. Yu (SCK-CEN, Belgium), E. Treille and J. Wendling (ANDRA, France), M. Dymitrowska and D. Pellegrini (IRSN, France), D. Justinavicius (LEI, Lithuania), M. Sentis (ENSI, Switzerland), S. Norris (NDA RWMD, UK), A. Bond (Quintessa, UK), H Leung (NWMO, Canada) and N. J. Calder (Geofirma, Canada).

- © The Geological Society of London 2015