At the post-closure stage of a geological disposal facility for higher activity radioactive waste several species of gas are likely to be generated in the near-field environment. These could alter the sealing and chemical properties of the bentonite buffer and the local geochemical environment significantly. The authors' attempt to simulate multicomponent gas flow through variably saturated porous media is presented. Governing equations have been developed for a reactive gas-flow model to simulate the thermo-hydro-gas-chemical-mechanical behaviour, with specific reference to the performance of highly compacted bentonite buffer subjected to repository gas generation and migration. The developed equations have been included in the bespoke numerical model COMPASS and some generic simulations are also presented. The model presented extends current capability to assess buffer performance.


In a deep geological disposal facility of higher activity radioactive waste various species of gases are likely to be generated due to mechanisms including water vaporization and radiolysis, corrosion of the metal canister and overpack, biodegradation and microbial activity (e.g. Bonin et al., 2000; Ortiz et al., 2002). The generation, accumulation and release of gases from the disposal facility may significantly influence the long term radiological safety of the repository via a number of processes (Nuclear Energy Agency, 2001).

The anaerobic corrosion of ferrous elements causes the generation of hydrogen which is considered to be the major contributor of repository gases (e.g. Rodwell et al., 1999). Other gas species which are of potential interest include methane (CH4), carbon dioxide (CO2) and water vapour.

In the context of higher activity radioactive waste disposal, most of the research activities on gas migration have been carried out based on a fully saturated buffer. It is considered that prior to gas generation, uniform saturation of the buffer took place and all the voids have been occupied by porewater, causing the buffer to be nearly impermeable to the generated gases. Several research teams have carried out both laboratory based experiments and numerical simulations to understand the gas breakthrough and flow phenomena in fully saturated clay buffers (e.g. Pusch and Forsberg, 1983; Horseman and Harrington, 1997; Harrington and Horseman, 2003). Experiments were conducted with various boundary conditions, including variation of gas injection rates or pressures, resulting in the build-up of large pressures and consequential fracture/dilation processes that could occur, although no models to explicitly include fractures and fracture propagation have been developed.

To advance the knowledge and understanding of gas related issues in the context of higher activity radioactive waste disposal, more research activity focussing on gas flow through a partially saturated buffer is essential. Anaerobic corrosion of metallic canisters and consequent generation of hydrogen gas could initiate before the buffer becomes fully saturated. Senger (2008) reported that the corrosion of metallic canisters starts at a relative humidity (RH) of 60% and increases dramatically at a RH of 90% which corresponds to a degree of saturation of ~50–90% . Therefore, the generated hydrogen together with other repository gases might flow through the unsaturated portion of the buffer. Also, the prediction of uniform resaturation is questionable due to complex nature of clay–water interactions. In reality, possible non-uniform resaturation could lead the porewater to reach various spots of the canisters and set off corrosion due to local anaerobic condition.

The flow of repository gases through a variably saturated buffer could undergo various geochemical processes. As bentonite clay buffers contain significant amounts of chemicals and minerals, gas-chemical/geochemical interactions are of great importance. Therefore, to ensure the long term safety and performance efficiency of bentonite buffers, research activities focussing on both experimental and numerical modelling of multicomponent reactive gas transport through compacted clays are essential. The aim of this research is to address the important gas migration issues by developing a multicomponent reactive gas transport model under a coupled thermo-hydro-chemical-mechanical (THCM) framework. The completed model will be an advanced tool to simulate the transport and fate of gases in unsaturated soil and to predict the long term behaviour of physical, chemical and mechanical properties of the soil, particularly in compacted bentonite buffer.

In this paper, the authors' attempt to demonstrate some of the developments of multicomponent reactive gas transport model. The coupled model consists of a gas transport module and reaction module. The numerical code COMPASS (code for modelling of partially saturated soil) has been used to develop the transport module whereas the geochemical code PHREEQC, version 2.0 (Parkhurst and Appelo, 1999) has been used to model the geochemical reactions. The aim of this paper is to present the development of multicomponent gas transport processes in COMPASS, such as multicomponent advective and diffusive gas flow. The implementation accuracy of the sink/source term (or the linking between COMPASS-PHREEQC) is illustrated by a simple example. The detailed gas-chemical/geochemical processes are beyond the scope of this study and are not presented in this paper.

Proposed model

The newly developed coupled model is capable of simulating moisture flow, multicomponent gas flow, multicomponent dissolved chemical flow, heat flow and mechanical deformation. The existing THCM model has been extended by introducing multicomponent gas transport processes. The previous developments of THCM model and its modelling capabilities have been presented in the work of Thomas and He (1998) and Seetharam et al. (2007). Therefore, detailed discussions of moisture flow, heat flow, dissolved chemical flow and mechanical deformation are not repeated here. Instead some developments in multicomponent reactive gas transport are briefly described.

The governing equation for multicomponent gas transport is based on the law of mass conservation.


[the net rate ofchange of moles]±[net rate of moles gained/lostdue to the source/sink]=[net rate ofmolar flux]

In terms of primary variables, e.g. porewater pressure (ul), temperature (T), gas concentrations (cg) deformation (u) the governing equation yields:  
where, C and K represent the storage and flux terms, respectively. The subscripts, c, g and w stand for storage or flux of a gas due to flow of moisture (l), temperature (T), other gas components (g) in the system and mechanical deformation (u); the subscript and superscript i, j is a gas component counter and N is the total number of gas components in the system; Ssi represents the sink/source term of ith component; t and ∇ are time and gradient operator, respectively.
The transport processes in a gas phase are dominated by the mechanisms of advection and diffusion. Darcy's Law has been adopted to obtain multicomponent advective flow of gas components in a gas phase. In a mixture of gases, the advective flux of individual species can be calculated from the proportion of that species in the bulk flow since the mixture is non-separative (e.g. Mason and Malinauskas, 1964). The advective flow of ith component in terms of molar flux, Jadvi can be expressed as:  
where, xgi is mole the fraction of gas component i and Jadv represents the total flux of gas mixture. By using Darcy's law the total gas flux yields:  
where c is the total mole concentration, R is universal gas constant, T represents absolute temperature at isothermal condition and kg is unsaturated gas conductivity.
Fick's law of binary diffusion has been extended to develop the governing equation of multicomponent gas diffusion through unsaturated soil. The diffusive flux of ith component Jdiffi can be expressed as:  
where, n is the porosity of the medium, Sg represents the degree of gas saturation, τ is the tortuosity factor due to flow paths of porous medium. In this case, τ = nS7/3g; Dij represents the multicomponent diffusion coefficients. The calculation of the multicomponent diffusion matrix has been described in detail by Taylor and Krishna (1993).

The sink/source term of the governing gas transport equation (1) is determined based on the equilibrium chemistry of geochemical reactions. Important gas-chemical/geochemical processes, such as ion exchange, mineral precipitation/dissolution, surface complexation, gas phase equilibrium and redox reactions can be calculated using geochemical reaction modelling via PHREEQC2. The transport model and geochemical reaction model has been linked using sequential non-iterative approach (SNIA). After convergence with the transport model, the program proceeds to the geochemical model and during a single time step the chemical equations are solved only once. Therefore, it has been termed as sequential non-iterative approach.

The higher order partial differential equation has been discretized using finite element method (Galerkin weighted residual method) for spatial discretization and finite difference method for temporal discretization.

Problem setup

Multicomponent reactive gas transport research has received very little attention and therefore there is a significant lack of available information to validate the developed model. Verification of individual processes implemented in the model has been presented in this study. Three sets of simulations have been carried out in this context to model multicomponent gas transport through a variably saturated compacted bentonite buffer. The simulation scenarios have been discussed briefly below.

Multicomponent advective gas flow through a dry compacted bentonite buffer has been modelled in first simulation. Once the corrosion and consequent hydrogen generation has begun, the flow of other gas components present in the system is affected by the inflow of hydrogen gas. As the developed model is aiming to model multicomponent diffusive gas flow, the calculation of a multicomponent diffusion matrix has been included following the approach suggested by Taylor and Krishna (1993). The implementation of a multicomponent diffusion matrix in general mass conservation equation is complex, as it requires the information about n + 1 components to define an n component system. To address this issue of multicomponent diffusion, a comparison study has been carried out to investigate the flow behaviour of gases under multicomponent diffusion and binary diffusion approaches. The aim of the simulation is to check whether the multicomponent diffusion coefficients can be replaced with the binary diffusion coefficients of gases through air. More details have been provided later in this paper.

The implementation accuracy of the sink/source term in the numerical model has been verified in a third simulation. The sink/source term of the governing equation (1) has been calculated using PHREEQC2. The amount of gases may vary due to possible geochemical reactions, such as ion-exchange, dissolution or precipitation, surface complexation and gas–liquid phase equilibrium. One of the simplest forms of verification is by calculating the sink/source term which represents the exchange of the phase e.g. gas phase to liquid phase equilibrium using Henry's law. Theoretically implemented Henry's law as a sink/source term in the developed COMPASS model has been tested against the coupled COMPASS-PHREEQC2 model. The database of geochemical model PHREEQC2 contains gas–liquid phase equilibrium reactions which has been used without any modification.


It has been assumed that the generated hydrogen flows through the unsaturated portion of the buffer at a rate similar the generation source. The simulations have been carried out under isothermal conditions. The gas mixture and individual components have been assumed to behave as ideal gases and the mixture is diluted.

Model geometry and discretization

A representative domain of a bentonite buffer of 30 cm by 3 cm has been considered in this study. The domain has been discretized into 80 elements of 4-noded quadrilaterals. Time steps of 100 s and rate of convergence 0.02% has been considered for the simulations.

Data and model parameters

Only hydrogen gas generation has been considered at the canister-buffer interface in this study. NAGRA (2004) suggested a range of hydrogen generation rate from a realistic rate of (2.938 × 10−11 kg s−1) to a conservative rate of (1.049 × 10−9 kg s−1) based on anaerobic corrosion of metallic canisters in a higher activity radioactive waste disposal. In this paper, the realistic generation rate is in terms of molar flux.

The buffer material considered in this study is bentonite Mx80 compacted at a dry density of 1600 kg m−3. Considering a porosity of approximately 0.4, an intrinsic permeability of 1.0 × 10−21 m2 has been reported by (Horseman and Harrington, 1997). The unsaturated gas conductivity, kg can be obtained using relative gas permeability, krg intrinsic permeability, Kint and dynamic viscosity, μg of the gas mixture as:  
Relative gas permeability has been developed using the van Genuchten-Mualem relationship as:  
where, m is a fitting parameter and in this case the value is close to 0.43. The effective saturation, Se can be obtained from the following equation:  
where, Sr,l is residual degree of saturation. In this case the value is equal to 0.40.

The diffusion coefficient data used in the multicomponent diffusive flow simulation has been presented in Table 1.


The initial and boundary conditions of multicomponent advective gas flow simulation is presented in Fig. 1. Hydrogen (H2) and carbon dioxide (CO2) are considered in this simulation as two representative gases. The model domain is fully gas saturated and hydrogen inflows from the left edge of the buffer.

Multicomponent diffusive gas flow has been demonstrated in second simulation. Two sets of simulation have been carried out in this case. Initially, simulation with multicomponent diffusion coefficients has been undertaken. The model calculates the elements of a multicomponent diffusion matrix at every time step. Using the initial concentration presented in Fig. 2 the calculated multicomponent diffusion matrix takes the form:  

In the second simulation, only the binary diffusion coefficients of gases through air have been used (data: first three rows of Table 1). The conditions of the simulation have been presented in Fig. 2. At the right face of model domain, x = 0.3 m, concentrations of oxygen, hydrogen and carbon dioxide have been fixed at 50.0 mol m−3; the left edge is impermeable.

The third simulation considers hydrogen flow through an unsaturated buffer. The amount of gas that is lost from the gas phase to liquid phase can be calculated using Henry's law. Two simulations have been carried out. Henry's law has been adopted theoretically in the COMPASS model at the first simulation. Later COMPASS-PHREEQC2 linked program has been used for the same simulation condition. A schematic of the simulation condition is presented in Fig. 3.

Results and discussions

The results of the multicomponent advective gas flow simulation are presented in Fig. 4. The gas pressure distribution profile along the buffer is at the end of simulation period of 10,000 years. The constant inflow of hydrogen pushed the carbon dioxide front away from the injection face due to advective flow. The calculated mole fraction of hydrogen near the injection face is equals one whereas for carbon dioxide it is zero. As the sum of mole fractions is always equal to unity; the results suggest good implementation accuracy of the model in simulating multicomponent advective gas flow.

The evolution of gas concentrations due to multicomponent diffusion (simulation 2) is presented in Fig. 5. The graphs have been plotted at the left edge of the buffer, x = 0. Gases with higher diffusion (hydrogen) coefficients attain steady state faster than that of (oxygen and carbon dioxide) lower diffusions. By comparing the results of Figs 5 and 6, it has been found that the multicomponent diffusive gas flow is dominated by the diagonal elements of the multicomponent diffusion coefficient matrix, [D] and these values are similar to the binary diffusion coefficients. The simulation considering binary diffusion has shown a similar flow pattern to the multicomponent diffusive flow pattern except the initial slight dip in oxygen and carbon dioxide curve which arise due to the cross diagonal elements of the multicomponent diffusion matrix. Therefore, the complex application of multicomponent diffusion can be substituted with simple binary diffusion in case of diluted gas mixtures.

The results of a third simulation are presented in Fig. 7. In the presence of a liquid phase, some of the hydrogen dissolves in porewater following Henry's law to maintain the equilibrium. Simulations have been carried out for 1000 years and the evolution of gas concentration (or pressure) using the COMPASS analysis follows a similar trend with those obtained by modelling using the linked COMPASS-PHREEQC2 analysis. It can be concluded that the approach adopted for linking COMPASS model (transport model) and geochemical model (PHREEQC2) via the sequential non iterative approach has been correctly implemented in the code in terms of gas phase exchange.


An investigation of multicomponent pore gas flow through a variably saturated buffer has been carried out in this work. In the case of multicomponent advective and diffusive flow the modelling capabilities have been demonstrated. Implementation accuracy of the sink/source term and the coupling technique have been verified.

P1Freely available online through the publisher-supported open access option.