The generation and migration of gas within and around proposed radioactive waste disposal facilities is potentially a safety critical process. A safety case for a facility that generates significant quantities of gas (e.g. through metal corrosion or radiolysis) will require demonstration that gas migration around and away from the waste is sufficiently understood and will not breach the safety case for the facility. Models can be used to understand the likely hydraulic evolution of such a disposal facility, but the models need to consider processes over a range of scales. Awhole repository may extend over kilometres, with individual disposal cells at the scale of tens of metres and features which provide pathways for gas migration on a centimetre scale. All of these features may be significant from a safety perspective and capturing the impact of all of these features in a single model is a significant challenge.

This paper presents an approach to tackling this multi-scale problem, which allows the whole repository to be modelled in a computationally efficient manner. The approach involves identifying areas within the modelled domain that show very similar behaviour, and representing these areas with sub-models, so that small-scale features are retained, but computational overhead is decreased by using the results in more than one location in the model domain. The approach allowed a model of a whole repository to be run on a single processor core, whilst maintaining the small-scale features of the system. The model results were compared against more conventional upscaling techniques and show the advantage of a more detailed representation of small-scale features. The model results reflect the conceptual understanding of how gas would migrate in a repository.


Migration of gas in a radioactive waste disposal facility is a multi-scale problem involving gas movement in features on the centimetre scale over distances of kilometres. Conventionally, representing such multi-scale systems would require a ‘homogenization’ step whereby features smaller than the desired scale of discretization are combined and the parameterization adjusted to compensate. In practice, for multi-phase flow applications, this involves taking the volumetric weighted mean for scalar properties such as porosity used in the volumetric conservation equations, and adopting various weighted means for the flux-based tensor properties (e.g. intrinsic permeability) depending on the geometry of the homogenized features relative to the volume boundary transfer. For linear problems the equivalent properties are often trivial to calculate; however for non-linear systems true equivalence is much more difficult to establish especially under strongly transient conditions. In all cases however, some information is ‘lost’ in the process of homogenization; the degree to which this is significant depends on the end application of the model. In this particular case, given the emphasis on the behaviour of small-scale structures, it was considered that a homogenization approach might give rise to unacceptable simplifications and so an alternative was sought.

Previous authors have addressed this problem in a number of ways including the use of simple nested models and more sophisticated ‘architectural elements’ (McDermott et al., 2006, and references therein) and a similar approach to this paper using finite volumes in the parallelized code TOUGH2-MP (Poller et al., 2011). In Poller et al. (2011) the modelled domain is subdivided into sections and these sections are represented by repeated units, but the repeat units are not fully coupled. However, when small-scale features are prevalent throughout large parts of the modelling domain, as is the case in he present work, these approaches would not be suitable for modelling the system. This paper presents a novel methodology that is based on a control volume approach, rather than a finite element approach, in which sub-models are embedded into a single model of the whole repository.

The work described in the paper was carried out for NDA RWMD (now Radioactive Waste Management Limited) as part of the EC Fate Of Repository GasEs (FORGE) project. The models herein were built as part of a benchmarking work package with models based on the French high level waste (HA-cell) disposal concept, which consists of High Level Waste (HLW) placed in steel canisters and installed in horizontal boreholes from the side walls of human-accessible tunnels. Further information on this work, including benchmarking of the different teams can be found in Wendling et al. (2013ac).

Model specification

The system modelled was specified in detail as part of the EC FORGE project (Wendling, 2013c). A repository containing 1000 HA cells, divided into groups of 100 cells that make up ‘modules’ was modelled (Fig. 1). The modules were linked by a main access drift with large bentonite seals emplaced along the drift. Between these seals, the module access drifts project out horizontally at 90° with the HA cells attached horizontally on both sides of the drift. The individual HA cells consist of a horizontal borehole 1 m in diameter containing a 5 m long bentonite plug at the access tunnel end and the remaining 40 m borehole length consisting of the emplaced HA waste packages. Excavation/Engineered Damage Zone (EDZ) structures were assumed to be present around all the excavations (zones of typically low capillary pressure and relatively high permeability compared with undisturbed rock, e.g. De La Vaissière, 2014) and narrow interface gap structures were assumed to be present on the upper half of all bentonite plugs, wastes and backfills. Around the main drift seals, it was specified that the EDZs have been removed and that the seal extends to undisturbed rock. The benchmarking specification (Wendling, 2013c) also required that the properties of the EDZ and interfaces were not to change with time. Connection to the surface was represented by a shaft with associated sealing structure (Fig. 2).

Bentonite plugs in the main drift separate each module, with a final drift plug in the drift that leads to the shaft. There is also a bentonite plug in the shaft, just below the aquifer layer.

A simplified geology was adopted with a low permeability (1 × 10−20 m2 horizontally, 5 × 10−21 m2 vertically), moderate porosity (15%), homogenous argillite (clay) host rock that was overlain by a more permeable ‘aquifer’ formation with a 1:100 hydraulic gradient running along the long axis of the facility, with the higher heads at the shaft end.

The following processes were specified (Wendling, 2013c): (1) Multi-phase Darcy flow of water, water vapour and hydrogen gas. (2) Dissolution and advective/diffusive transport of hydrogen in water. (3) Advective and diffusive transport of hydrogen in water vapour. (4) Specific models for the equations of state for water and gas, water vapour diffusivity and water and hydrogen gas viscosity. (5) Constant gas generation of 100 mol/y per HA cell for 10,000 years from closure – approximately 20 kg/yr per module. This rate of production is a simplification of the calculations of hydrogen production via corrosion of metals and radiolysis for the French disposal concept (Wendling, 2013a).

Initial conditions assumed that the interfaces, backfill and bentonite were at atmospheric gas pressure and partially saturated, while the EDZ and host rock were fully water saturated and at hydrostatic water pressure. For simplicity, the initial gas saturation was assumed to be hydrogen rather than air. The upstream and downstream boundaries of the aquifer were set to specified water heads and free gas was allowed to escape through the upper surface of the aquifer; all other external surfaces of the model were considered to be impermeable.

Numerical formulation

A conventional finite volume (or control volume) approximation was adopted (Versteeg and Malalasekra, 2007). In the finite volume method, volume integrals in a partial differential equation that contain a divergence term can be converted to surface integrals, using the divergence theorem, enabling the original partial differential equation to be replaced by ordinary differential equations for the volume-averaged variables. The advantage of using the finite volume method over the commonly used finite element approach (Zienkiewicz et al., 2005) in this work is that in the finite volume approach there is a uniquely defined flux between volumes, ensuring that the relevant quantity is conserved locally. In the finite element approach, there is no unique flux between elements, the normal gradients will generally differ across the element boundary, and so local conservation is lost. For this application the ability to guarantee local mass conservation and have a clear separation between geometry and physical process makes the implementation of the proposed multi-scale technique considerably simpler, while guaranteeing global and local mass conservation.

The Darcy multi-phase flow and vapour transport equations can be readily transferred into the control volume framework and a brief summary of the relevant relationships is provided below; full details are provided by Wendling et al. (2013c).

We specify that the fluids occupy the total available pore space such that the fractional volume occupancy of each phase (saturation, S) sums to unity, i.e.  
where the subscript w denotes the wetting phase (water) and nw is a non-wetting phase. For this analysis we assume there is only one wetting phase (water) and one non-wetting phase (hydrogen). The difference in pressure (p, MPa) between the wetting and non-wetting phase is normally termed the capillary pressure or suction (pc) and for a two-phase system is defined as a function of the saturation of the wetting phase:  
If we denote the fluid accessible porosity of the porous medium by θ(−) and the density of phase k by ρk (kg m −3), then θρkSk (kg m−3) is a conserved quantity in any representative volume of porous media, which inversely relates changes in density to changes in porosity or saturation. If the Darcy velocity of fluid k is denoted uk (m s−1), the fluid vapour flux or dissolved hydrogen flux is denoted by vk (kg m−2 s−1) then we can write the equation for the conserved quantity as  
where qk (kg m−3 s−1) is a source or sink rate of fluid per unit volume. This basic conservation equation can be related directly to equation (1).
The Darcy velocity of a specific fluid in a multiphase flow setting is given by a simple extension of Darcy's law for single phase flow (Chen et al., 2006).  
Here, g (m s −2) is the acceleration due to gravity, z (m) is the vertical coordinate and μk (Pa s) is the viscosity of fluid k. kk (m2) is the effective permeability tensor of the porous medium for fluid k which is a function of fluid saturation. The vapour flux vw and dissolved gas is flux vnw are defined as follows  
where Dv is the net vapour diffusivity (m2 s−1), ρv is the vapour density (directly related to the capillary pressure; Engel et al., 2003), Dnw is the net diffusivity of the non-wetting phase (hydrogen) in water (m2 s−1) and cnw is the concentration of the non-wetting phase in water (kg m−3).

For this study, a formulation of the above relationships that has been successfully used in a range of other applications was adopted (Benbow et al., 2014; Bond et al., 2013ac; Bond and Watson, 2012; Towler and Bond, 2011). In this formulation the primary conserved variables are the masses of fluid in each control volume, while the fluid pressures and capillary pressure are calculated as algebraic variables through the fluid compressibility and rock matrix elasticity terms while satisfying the saturation constraint equation (equation 1), giving rise to five unknowns per control volume.

The equations are solved using a monolithic, fully implicit technique, implemented in ‘QPAC’, a proprietary multi-physics code developed by Quintessa Ltd (Quintessa, 2012; www.quintessa.org/qpac). Non-linear iterations for each time step are solved using a Newton-Raphson scheme while the linearized equations are solved using a GMRES iterative solver. Time-stepping is managed using a polynomial history-derived predictor-corrector method based on the work of Byrne and Hindmarsh (1975). Derivatives are calculated symbolically to allow arbitrary, non-linear and mutually dependent definitions of equation properties to be specified without recompiling the code. In order to give flexibility in the physical process models being solved, ‘sub-systems’ can be specified which use different physical relationships that are then coupled together via ‘joiners’ to create a single set of equations to be solved (see Fig. 3). Such sub-systems are arbitrary in nature and can be joined to represent continuous sub-domains within a volume, associations between non-contiguous geometries or even separation of different physical process models within a single volume, depending on the application. The use of such sub-systems is critical to the modelling approach outlined in this paper.

In all cases the models were run on conventional desktop computers with no parallelization.

Approach and implementation

The geometry of the repository is such that the majority of HA cells are situated in similar geometrical positions, i.e. between two other HA cells, attached at one end to a drift and with host rock followed by another cell at the other end. Results from detailed models of a single module, where the HA cells were represented explicitly, indicated that HA cells in similar positions behave in a very similar way, with only the HA cells at either end of the module showing significantly different behaviour (Wendling, 2013a,b). This similarity in localized behaviour could be exploited by representing the majority of HA cells by a single detailed sub-model of one HA cell. HA cells that might behave differently, e.g. at the ends of modules and at the ends of the repository, could then be represented by different detailed submodels. Interaction between detailed sub-models of individual HA cells could then be fully coupled into the larger-scale model grid at numerous locations, making use of the ability to explicitly conserve fluxes between volumes in the finite volume method.

This approach allows the smaller-scale features of the HA cells to be retained, the compromise being that identical behaviour of the sub-modelled HA cells is enforced across sections of the repository. The challenge was to use a sufficient number of detailed sub-models to capture the variation across the repository, without creating an excessive numerical overhead.

As specified, the whole repository has bilateral symmetry down the main access drift, and thus the modelling domain could consider only half the facility with no loss of accuracy. During testing, interface structures were found to add a considerable numerical overhead due to their specification as high intrinsic permeability and very low air entry pressure features, so they were omitted for the base case model and included as a variant calculation.

The HA sub-models were fully coupled to the overall grid (Fig. 4) using scaled fluxes to allow a single HA sub-model to represent an average behaviour of multiple HA cells at multiple locations. The model was set up such that volumes of rock containing sets of HA cells with the same expected behaviour were represented by large compartments. All compartments representing volumes that contain HA cells were then discarded but the boundaries between the discarded compartments and the rest of the model remained. The sub-models were then used in place of the discarded volumes with appropriate scaling of fluxes. This flux scaling required that mass fluxes on a single transfer applied to the sub-model and the larger-scale model were different in order to give overall mass conservation. The plug of the HA cell in the sub-model connected directly to the tunnel backfill in the larger-scale model with the other sides of the sub-model connecting to the host rock. If a submodel is used in multiple locations within the overall model then multiple transfers are created on each sub-model boundary surface, one for each use in the larger-scale, whole repository model. On the basis of earlier modelling work, flow symmetry boundaries between HA cells was assumed, so submodels were not connected to each other directly.

The reference mass fluxes (qo,i,j) between the sub-model i and an associated whole repository discrete volume j were first calculated based on the pressure and concentration gradients between the sub-model compartment and the whole repository compartments, using the compartment face areas associated with the sub-model (equations 2, 3, 4, 5 and 6). These reference mass fluxes were then scaled such that the correct flux was seen at the whole repository scale (qw) and the sub-model scale (qs);  
where ni,j is the number of HA cells being represented by sub-model i at discrete volume j in the wider model. Ni is the total number of HA cells a sub-model is representing, i.e.  

The above relationships mean that for fluxes applied on the sub-model, the volume (total volume represented by the sub-model in the whole repository grid) weighted arithmetic mean of reference fluxes across all connections between a given boundary in the sub-model and the whole repository model was used. The net fluxes applied to the sub-model therefore represent an average behaviour of the connections to the whole repository model when the sub-model is used multiple times. The mass flux applied to the whole repository model side of each transfer is the reference flux multiplied by the number of submodels being represented in the discarded volume; the flux is simply scaled to reflect the number of HA cells that the sub-model represents. The flux scaling is illustrated in Fig. 5 for a particular simplified geometry. As implemented, each HA cell is a separate sub-system, thus making the flux scaling shown in Fig. 5 easier to implement. The approach has certain similarities to that adopted by McDermott et al. (2009) but is distinct in its use of using a single numerical method at all scales and control volumes to guarantee local mass conservation.

The initial distribution of sub-models used for the case is shown in Fig. 6. This was judged to be the sensible maximum averaging that could be adopted while retaining the general behaviour shown by modelling performed at the module scale (Wendling et al., 2013b). HA sub-model 1 represents only 2 HA cells in the whole model, compared to HA sub-model 6, which represents 384 HA cells. The consequence of selecting such a pattern was subject to extensive testing, as discussed in the following section.

Testing and variants

The model was built following ISO:9001 standards which in this case required that the model is tested to ensure that the outputs are self-consistent, physically plausible (e.g. mass conservation is retained globally and locally; no significant solution instabilities) and are not dependent on arbitrary model construction decisions. The standards as applied to this project also required sufficient documentation to allow the model to be re-implemented and re-used in the future. Initial testing was performed on ‘mini-models’ considering only one half of a module to allow by-hand checking of the fluxes between sub-systems and the larger-scale model, and to facilitate global and local mass conservation analysis. In addition, a grid convergence testing programme was conducted to ensure that the discretization scheme for the sub-models was sufficient to provide representative fluxes, pressures and fluid saturations. This was achieved through progressive and selective refinement of the sub-model and larger-scale model grids until the results remained essentially unchanged for a wide range of metrices. No flux stability issues were identified in the implementation of this approach even when using a single sub-model to represent multiple parts of the system.

Testing was also performed on the full model through consideration of three sets of variant calculations:

Sub-modelling: Changing and increasing the pattern of sub-models used in order to obtain a finer representation of the variation of HA cells. Cases involving 16 and 30 sub-models were implemented (see Fig. 7).

Interfaces: Inclusion of all interfaces on the upper half and side of engineered features to be fully consistent with the specification.

Upscaled modules: The use of simpler homogenization techniques was also investigated. This involved not removing those portions of the larger-scale model that would be represented using the submodels and instead assigning effective properties using simple homogenization techniques.

The key results from these and other variants are discussed later in the paper; however the process of investigating these cases showed no significant issues with the base model and that the implementation and sub-modelling assumptions were robust.

Key results

Base case

The model results for the full repository model show that migration of H2 away from the facility is almost entirely diffusion dominated, with peak free gas fluxes to the shaft of 0.06 kg/yr being small in comparison to 200 kg/yr generated in the whole facility. The model testing process showed that while relatively small, the shaft fluxes were stable versus solver tolerances and grid refinement. Significant free gas is retained in the model (Fig. 8), but this is largely due to the very small amount of desaturation of a large volume of host rock and EDZ immediately surrounding the canisters. Free gas fluxes through the top of the argillite host rock are also small (0.1 kg/yr maximum, only significant after ~9000 years) and almost entirely caused by exsolution of dissolved gas as the water pressure drops, rather than the migration of a discrete gas phase from the disposal facility. The dissolved and free H2 reaching the aquifer is diluted by the relatively high flow of fresh water above the argillite, which helps to maintain an upwards diffusive gradient. A plume of dissolved gas moves through the aquifer from around 15,000 to 50,000 years, but the concentration of H2 is significantly lower than the concentrations around the waste canisters.

Gas pressures in the facility exceed hydrostatic pressures after ~200 years, peaking at ~8.5 MPa, compared with a hydrostatic pressure of ~5 MPa (Fig. 9). Note that the vertical lithostatic pressure is expected to be 13 MPa. Pressures only start to drop after 10,000 years when gas generation in the HA cells ceases and return to hydrostatic conditions after ~55,000 years. In contrast, at the base of the shaft no pressurization is seen and the seals resaturate smoothly to reach equilibrium conditions after 5000 years (Fig. 10).

The pressurization and dissolution dominance is caused by the main drift seals resaturating sufficiently quickly that the effective gas permeability of the seals are low enough to prevent significant gas advection. In contrast, the bentonite plugs in the HA cells are found to be largely ineffective in retarding gas advection because of the presence of an EDZ around the plugs and the proximal generation of H2 slowing the EDZ resaturation, and hence keeping an effective gas permeability around the plugs (Fig. 11).

The results from the base case model demonstrate that implementation of the homogenization approach works well, giving sensible results in line with the design concept of the repository.

Variant results

The sub-modelling variant cases, using more HA sub-models to represent the system showed that the basic scheme of using only 6 HA sub-models was robust relative to the overall behaviour of the system. The use of additional sub-models only showed small differences in pressure and gas flux at key evaluation points and surfaces and the overall evolution of the system in all cases was functionally identical.

Inclusion of interface structures had only a minor impact on the gas pressures in the facility, but did increase the flux of gas around the main plugs to the base of the shaft by a factor of ten. However, this enhanced flux of 0.25 kg/yr for the half-model is still very small in comparison with the 100 kg/yr generated by the HA cells in the half-model. Even with the interfaces present, the vast majority of the H2 gas is transported in solution.

The most interesting results related to the comparison of the new sub-modelled approach versus the more conventional homgenization upscaling techniques. There are clear differences between the pressures in the modules depending on whether sub-models or an upscaled, homogenized representation is used, particularly in the timing of the evolution, but the overall form of the behaviour is similar (Fig. 12). The modelled reduction in gas pressure in the host rock during resaturation was smaller and later when using sub-models compared to the upscaled approach – this is consistent with the increased discretization and more accurate material representation in the sub-models. Pressures in the drift increase more rapidly in the sub-models, as the pathway for gas transport from the canister to the drift is explicitly represented and maintained. However, late time behaviours are very similar. Overall, the results between cases are coherent and illustrate the advantages of the sub-modelling approach in capturing some key behaviours as well as confirming the basic response of the sub-modelled response is correct.

While it might be argued that the differences for the overall system behaviour are relatively minor, unless there is a clear safety requirement to capture early transient behaviours, it is clear that the localized behaviour in the HA cells are much more realistically captured when using the submodelled representation. This is especially true when considering the immediate vicinity of the waste. Such differences between approaches will be magnified when non-linear coupling effects, such as using more complete models of gas generation where corrosion processes are linked closely to water availability and temperature. Indeed, such an approach allows for complex physical and chemical processes which are often very restricted in the spatial extent to be represented directly in a whole repository model, a key direction for future work in this area.

The results of the variant cases illustrated that the selection of sub-modelled volumes was appropriate and that key information on the early transient behaviour is retained in the sub-modelled representation that would otherwise be lost when using homogenization techniques.


The work has shown that it is practicable to model multi-phase flow physics of an entire radioactive waste disposal facility on a desktop computer while explicitly retaining key features down to the cm scale, something that was found to be difficult to achieve using more conventional techniques. This approach therefore provides another technique that can be used in conjunction with other homogenization/simplification approaches to assess the evolution of radioactive waste disposal systems.

The main feature of the model which enabled such a representation was the ‘sub-modelling’ approach, which enabled averaging over features that are repeated many times in the repository, in this case, the HA cells. This approach considerably reduced the discretization required to represent the modelled repository and hence reduced the computational cost of running the model, while at the same time retaining explicit representation of very small but potentially important features in the system (such as EDZ and interface structures).

The models have been designed to represent ANDRA's disposal concept, in which the seals are intended to prevent free gas migration to the shaft (which when sealed may provide a less robust barrier than the natural host rock, and hence permit relatively more rapid transport of radioactive gas to the surface), and ensure the dominant phase in which gas is transported away from the repository is as a dissolved phase in the host rock. We have demonstrated that, for the model parameterization used, very nearly all of the gas generated within the repository dissolves in water in the host rock, with very little advection of free gas outside the main disposal area. The main drift plugs are very effective in sealing the system, but if they were to fail or work less effectively, there would be significantly more advection of gas towards the shaft. The results of the models are also sensitive to the presence or absence of the very small interface structures during the early transient period.

Comparison of this approach to conventional upscaling techniques has demonstrated that key information may be lost when not considering small-scale features, but also demonstrated that the overall, large-scale system behaviours are similar, building confidence in the sub-modelling approach. Extending this approach to include thermal processes, coupled gas generation and other very localized non-linear processes, which are very difficult to capture accurately using conventional homogenization techniques, is possible without further code modification. Such work has been performed and will be published in due course.


The funding from NDA RWMD and technical review from NDA RWMD and the other FORGE partners in Work Package 1 are gratefully acknowledged.

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.