Development of a mathematical model for gas migration (two-phase flow) in natural and engineered barriers for radioactive waste disposal

Abstract In a deep geological repository (DGR) for the long-term containment of radioactive waste, gases could be generated through a number of processes. If gas production exceeds the containment capacity of the engineered barriers or host rock, these gases could migrate through these barriers and potentially expose people and the environment to radioactivity. Expansive soils, such as bentonite-based materials, are currently the preferred choice of seal materials. Understanding the long-term performance of these seals as barriers against gas migration is an important component in the design and long-term safety assessment of a DGR. This study proposes a hydro-mechanical linear poro-elastic visco-capillary mathematical model for advective-diffusive controlled two-phase flow through a low-permeability expansive soil. It is based on the theoretical framework of poromechanics, incorporates Darcy's Law for both the porewater and poregas, and a modified Bishop's effective stress principle. Using the finite element method (FEM), the model was used to numerically simulate 1D flow through a low-permeability expansive soil. The results were verified against experimental results found in the current literature. Parametric studies were performed to determine the influence on the flow behaviour. Based on the results, the mathematical model looks promising and will be improved to model flow through preferential pathways.

In Canada, nuclear waste has been generating and accumulating since the 1930s when the Port Radium radium mine began operating in the Northwest Territories (Low-Level Radioactive Waste Management Office 2012). Since then Canada has become a sustainable nuclear country with operating nuclear facilities across the nuclear fuel cycle, all producing various forms of radioactive waste. To-date, this waste, in the form of low-, intermediate-and highlevel radioactive waste has been primarily stored on-site at nuclear power plants or in below-surface radioactive waste management facilities.
The Low-Level Radioactive Waste Management Office (LLRWMO) estimated that by the end of 2050 there would be more than 2.6 million m 3 of radioactive waste requiring long-term management in Canada, the bulk of which is low-level radioactive waste. Table 1 provides waste inventory estimates by the end of 2011 and by the end of 2050.
In identifying a long-term solution to manage the global radioactive waste, the international community, including Canada, have been assessing the technical practicability of a deep geological repository (DGR) for its long-term management. The primary purpose of a DGR is to contain and isolate wastes to minimize impact on the environment and radiological exposure to people.
To adequately design a safety case for a DGR, consideration must be given to the relevant features, events and processes (FEPs), and their impact on the primary objectives of a DGR (CNSC 2006;Norris 2015). One such process with the potential means for radiological exposure to the biosphere is the generation of radioactive gas which may migrate to the surface (Shaw 2015a). Gas could be generated through a number of processes including the degradation of organic matter, radioactive decay of the waste, corrosion of metals producing hydrogen gas (H 2 ) and the radiolysis of water producing H 2 (Cuss et al. 2014;Ahusborde et al. 2015;Shaw 2015b). If production exceeds the containment capacity of the engineered barriers or host rock, these gases could migrate through the engineered barriers and/or the host rock Sellin & Leupin 2013). The preferential migration pathway of these radioactive gases, to potentially expose people and the environment to radioactivity, might be through the access and ventilation shafts; as these components are typically part of the repository design.
In recent years, a number of international projects have focused on the topics of gas generation and migration. In Europe, the Fate Of Repository GasEs (FORGE) project was undertaken to address key research areas on gas generation and migration, in order to better ascertain the mechanisms of each, and to support the safety case for DGR (Sellin 2014;Shaw 2015b). As part of the FORGE project, the large scale gas injection test (Lasgit) was designed to study the impact of gas build-up and migration through an engineered barrier system (EBS) (Sellin 2014;Shaw 2015b).
Expansive or swelling soils, such as bentonitebased materials, are currently the preferred choice of seal materials used for an EBS. Understanding the long-term performance of these seals as barriers against gas migration is an important component in the design and long-term safety assessment of a DGR.
An investigation of the gas transport processes in low-permeability clay material was described by Marschall et al. (2005), and further investigated by Cuss et al. (2014). Marschall et al. (2005) recognized that gas transport through porous media is controlled by a number of the media's hydraulic and mechanical characteristics such as the intrinsic permeability, porosity and material strength. They identified the importance of the hydro-mechanical state of the rock or soil media (i.e. water saturation, porewater pressure and stress state) and the gas pressure at focal points which could lead to microfracturing (Marschall et al. 2005). Marschall et al. (2005) divided the basic transport mechanisms into four processes conceptualized in Figure 1.
In the first process, the advective-diffusive transport of gas dissolved in the porewater is governed by Darcy's law for advective groundwater flow, Fick's law for the diffusion of dissolved gas due to concentration gradients in the porewater and Henry's law which describes the solubility of gas in the porewater.
The second process of visco-capillary two-phase flow is related to the viscous and capillary forces (i.e. matric suction) which describe the underlying principles of unsaturated soil mechanics. Soil properties influencing this transport process are the pore-size distribution and soil compressibility. The air-entry value (AEV) or gas-entry pressure, which is the value of matric suction that must be exceeded before air recedes into the soil pores, is a controlling factor in this transport mechanism of two-phase flow. It is also a measure of the maximum pore size in a soil, as water will migrate from and gas will enter through the largest pores first.
A wealth of laboratory and field-scale experimental studies have investigated gas transport processes through natural (host rock) and engineered barriers. These studies have provided considerable evidence suggesting that gas flow at gas pressures above a reference level is accompanied by the creation of pressure-induced preferential pathways and dilation of the clay, yet have not been able to determine the exact mechanisms which control gas entry, flow and pathway sealing Horseman et al. 1999Horseman et al. , 2003Graham et al. 2012;Harrington et al. 2012Harrington et al. , 2017Cuss et al. 2014;Sellin 2014;Bennett et al. 2015). Marschall et al. (2005) described this primary transport mechanism as dilatancy-controlled gas  flow or 'pathway dilation' and identified it as an important transport mechanism in clay soil or clay-rich rock with low tensile-strength. Dilatancy will occur when the gas pressure reaches a reference stress level acting on the clay medium, forcing the clay particles to align in a dispersed orientation and resulting in the formation of microfractures accompanied by an increase in plastic deformation (Davis & Selvadurai 2002). As a result, gas flow along these microfractures will be promoted due to the increased pore space. The presence of microfractures, therefore, leads to an increase in the intrinsic permeability of the material and in turn results in changes in the relationship between the capillary pressure (i.e. matric suction) and degree of saturation (i.e. SWCC of the material). As a result of dilatencycontrolled gas flow, transport properties are now dependent on the stress-state and state of deformation of the soil. Finally, gas transport along macroscopic fractures, also known as hydro-and gas-fracturing, becomes a significant transport mechanism when a macroscopic tensile fracture occurs. Unlike the formation of microfractures that result in dilatancy, these macroscopic fractures develop in low-tensile strength material when gas pressure build-up is rapid and becomes larger than the sum of the minimum principal stress and the tensile strength of the material (i.e. porewater displacement and formation of microfractures can no longer counter balance the gas production rate) (Marschall et al. 2005).
A number of studies have attempted to model gas migration in natural or engineered barrier systems. In 1998 the results of Phase 1 of the GAMBIT Club programme were published (Nash et al. 1998). The main objective of the GAMBIT Club programme was to establish a computational model that could represent the principle features of gas migration through compacted bentonite observed through experiments (Nash et al. 1998). Nash et al. (1998 applied the theory of linear elastic fracture mechanics to develop a model of crack propagation through clay, whereby providing a continuous gas flow in a sample, following gas breakthrough, gas flow would be observed across the sample and dilation pathways occur. Phase 1 saw difficulties obtaining agreement between simulation and experimental results over the whole test sequence (Nash et al. 1998).
Phases 2 and 3 of the GAMBIT Club programme further explored the refinement and extension of their earlier models, including the investigation of clay fabric effects (i.e. distinction between interand intra-stack pore space) and how to treat hysteresis effects (Hoch et al. 2004). Hoch et al. (2004) examined results with three basic mechanisms that could describe flow, (1) gas flow governed by conventional concepts of capillary pressure and relative permeability, (2) microfissuring (dilation) of the clay and (3) macroscopic fracturing of the clay, but could not provide a unique explanation of the results in terms of these mechanisms. They concluded that the current uncertainties associated with understanding the transport mechanisms made it considerably difficult to establish an adequate model for gas migration (Hoch et al. 2004). The report also concluded that any further development of models of gas migration in bentonite will depend on obtaining better characterization of gas pathways as well as the couplings that exist in the clay between stress and strain, gas and water fluid pressures, and gas-filled porosity (Hoch et al. 2004).
Fall et al. (2014) developed a coupled hydromechanical model for simulating gas migration in host sedimentary rocks for nuclear waste repositories. Their study identified that high gas pressure can lead to mechanical damage and the formation of micro-cracks (Fall et al. 2014). Their model was limited in part because they did not consider the elasto-plastic mechanical behaviour of the rock. Nguyen & Le (2015) developed a mathematical model for gas migration in Opalinus clay. The inherent anisotropy due to bedding, and the elasto-plastic behaviour of the clay were considered in the model. The model was successfully validated against laboratory experiments and field gas injection tests, especially up to the dilatancy controlled flow phase. However, it was recognized that the last phase of gas flow controlled by macroscopic fractures need further understanding.
In 2015, the Geological Society of London, published Special Publication 415entitled, Gas Generation and Migration in Deep Geological Radioactive Waste Repositories. This publication presented current research to date on experimental studies on gas migration and the development of models to describe gas behaviour in several systems of a DGR (Shaw 2015b). However, there is still much work to be done to understand the basic mechanisms and processes of two-phase flow, particularly dilatancy-controlled gas flow through low-permeability swelling soils over the life-time of the repository. In doing so, the development of mathematical relations describing these processes is imperative, along with the development of suitable numerical models which can be used to support the design and long-term safety assessment of a DGR.
Stemming from these international initiatives to understand gas migration and multi-phase flow in low permeability geomaterials, Task A of the current project phase of the international working group for the DEvelopment of COupled models and their VALidation against EXperiments (DECOVALEX-2019) was established. This task, led by the British Geological Survey (BGS), further attempts to identify the physical hydro-mechanical (HM) mechanisms required to adequately model dilatancycontrolled gas migration. This research is in part a contribution to Task A of DECOVALEX-2019.
This paper presents a fully-coupled, hydromechanical linear-elastic mathematical model for advective-diffusive visco-capillary controlled twophase flow through geomaterials. The objective of the study is to gain an intricate understanding of the model processes and model parameters which can influence flow behaviour. A validation study is performed by numerically simulating experimental data from a 1-dimensional flow constant volume boundary condition (BC) experiment using the finite element method (FEM). This is followed by parametric studies and sensitivity analysis that are used to assess the effects of several features of the numerical model on the flow behaviour.

Mathematical model
An HM model to describe the migration of gas (twophase flow) through a low-permeable expansive soil was developed using the theoretical framework of poromechanics. The model follows the general formulation from Nguyen & Le (2015), with the addition of the following features: (1) a Bishop's effective stress, with a χ parameter generalized from the work of Khalili & Khabbaz (1998); (2) gas dissolution into the liquid phase and subsequent gas migration from (a) the advection of water, and (b) diffusion of gas through the liquid phase; (3) a relation for the air-entry value (AEV), and corresponding soil water characteristic curves (SWCCs), as a function of changing porosity; (4) a relation of the intrinsic permeability, k ij , as a function of changing porosity; and (5) consideration of a damage model. The applicable constitutive relations and governing equations for conservation of momentum, water mass and gas mass are presented below.

Constitutive relations for the mechanical behaviour
Bishop's modified effective stress principle. Many effective stress equations have been proposed to characterize the stress-state of an unsaturated soil or porous media (Nguyen & Le 2015). This paper proposes the use of Bishop's effective stress principle, which is dependent on both net normal stress and matric suction, and may be more suitable for expansive clay's, as described by equation (1).
where σ′ is the effective stress (Pa) σ is the total normal stress (Pa) p g is the poregas pressure (Pa) p w is the porewater pressure (Pa) (σ−p g ) is the net normal stress (Pa) ( p g −p w ) is the matric suction (Pa) χ is a parameter related to the degree of saturation of the soil (unitless). Expanding and rearranging for σ The porefluid pressure can be defined as, Khalili & Khabbaz (1998), proposed the following unique relationship for the determination of χ based on the ratio of suction over the AEV, also termed the suction ratio, where ( p g −p w ) b is the AEV of the soil and only applies when the matric suction >AEV (Khalili & Khabbaz 1998).
Equation of poro-elasticity. The general form for the equation of poromechanics can be expressed by, where σ ij total normal stress tensor acting on the soil element (Pa) C ijkl is the fourth order stiffness tensor (Pa) ε kl is the strain tensor (unitless) α B is the Biot-Willis coefficient (unitless) δ ij is the Kronecker delta (identity tensor) (unitless).
Applying Hooke's Law for linear elasticity where where G is the shear modulus (Pa) λ is the Lamé constant (Pa) ε kk (=u k,k = trε ij ) is the volumetric strain (unitless) u ij is the displacement tensor (m). substituting into equation (6) or in terms of Young's modulus, E, and Poisson's ratio, ν

Constitutive relations for the hydraulic behaviour
Equation for the SWCC. The SWCC describes the relationship between the degree of saturation, S or water content (θ or w) in the soil and the matric suction, s. A number of mathematical models can be used to model the behaviour of the SWCC. This paper utilizes the van Genuchten equation to model the SWCC; a specific soil fitted to model parameters a, n′ and m′ (van Genuchten 1980).
where S e,w is the effective degree of saturation of water a′, n′ and m′ are fitting parameters for the SWCC ρ w is the density of water (kg m −3 ) g is the acceleration due to gravity (m s −2 ) s is the matric suction ( p g −p w ) (Pa). The fitting parameter, a′, is related to the AEV and can be expressed by equation (12) a ′ = 0.5 r w g AEV The fitting parameter m′ is related to fitting parameter n′ by equation (13), S e,w can be calculated from scaling the water saturation, S w , with respect to the maximum degree of saturation, S w,s , and the residual degree of From the fitted SWCC, the model can be used to predict the degree of saturation of water, S w , from measuring the stress-state variable matric suction (or capillary pressure), s.
The degree of saturation of the gas phase, S g , can be related to S w via the following relationship It should be noted that the model does not consider entrapped gases which would be present if S w,s < 1.
Equation for the AEV as a function of pore size. Huang et al. (1998), using data published by Laliberte et al. (1966), showed that for sand, sandy loam and silt loam, the logarithm of the AEV is linearly proportional to the void ratio or porosity at the AEV, as expressed by equation (16), and this linear relationship can be assumed for deformable unsaturated porous media.
where s aev is the AEV (Pa) for a void ratio, e aev s′ aev is the arbitrary reference AEV (Pa) for a void ratio e′ aev m is the slope of the s aeve aev curve. It should be noted that the slope, m, is negative as an increase in the void ratio (i.e. pore size) will result in a decrease in the AEV. In the current model, the swelling potential of expansive clays is not considered and therefore this relationship is considered reasonable for this study.
By re-arranging equation (16) we can express this relationship in terms of porosity, Darcy's Law for two-phase flow. Darcy's Law for the water and gas phases are expressed in equations (18) and (19).
where p w is the porewater pressure (Pa) p g is the poregas pressure (Pa) (v iw −v is ) is the velocity of water in the pores relative to the velocity of the soil matrix (m s −1 ) (v ig −v is ) is the velocity of gas in the pores relative to the velocity of the soil matrix (m s −1 ) n porosity (m 3 voids · m −3 total) μ w dynamic viscosity of the water phase (Pa s or kg m −1 s −1 ) μ g dynamic viscosity of the gas phase (Pa s or kg m −1 s −1 ) k ij intrinsic permeability tensor of the porous medium (m 2 ) k rw relative permeability of the water phase (unitless) k rg relative permeability of the gas phase (unitless) ρ w density of water phase (kg m −3 ) ρ g density of the gas phase (kg m −3 ).

Intrinsic permeability as a function of porosity.
Numerous equations have been proposed to describe the relationship between intrinsic permeability and measurable properties of soils such as the porosity and particle-size. The Kozeny-Carman equation can be used to express the intrinsic permeability as a function of porosity using the capillary tube (hydraulic radius) model and is sufficient for laminar flow and low porefluid velocity (Carman 1956). The Kozeny-Carmen Equation is given by, where S 0 is the specific surface (surface exposed to fluid per unit solid volume) (cm 2 cm −3 ). For fine-grained expansive clays, which have very small grain-sizes less than 2 µm, and are nonuniform, the Kozeny-Carman equation may not be appropriate (Pall & Moshenin 1980;Valdes-Parada et al. 2009). Pall & Moshenin (1980), proposed a modification of equation (20) by taking the volumesurface mean diameter to better account for the nonuniformity of soil particles.
where D vs is the volume-surface mean diameter (cm). If the porosity and intrinsic permeability were known, the volume-surface mean diameter could be obtained by rearranging equation (21).
In this study, the experimentally determined initial porosity and intrinsic permeability were used to calculate the volume-surface mean diameter of the soil sample. Any elastic deformations to the soil structure will influence the porosity, subsequently influencing the intrinsic permeability across the soil specimen.
Relative permeability of water and gas phases. The coefficient of permeability is strongly dependent on the degree of saturation, S w . Using Mualem's Model (Nguyen & Le 2015), k rw and k rg , can be calculated using the following expressions: where L′ is the relative permeability function fitting parameter (unitless).

Diffusivity of gas in water through porous media.
Diffusion of gas through the liquid phase of a porous media can be expressed as an effective diffusivity, D e , D e = DnS w t (25) where D e is the effective diffusivity of gas dissolved in water through porous media (m 2 s −1 ) D is the diffusivity of gas in water (m 2 s −1 ) n is the porosity (unitless) τ is the tortuosity (unitless).
The tortuosity for gas diffusion through porous water can be calculated using the Millington and Quirk model (Ho & Webb 2006) Storage due to suction of water and gas (storage coefficient). Storage terms are applied in the governing equations when differentiating dS w /ds. The storage coefficient of water, C w , can be calculated by differentiating the van Genuchten SWCC, The storage coefficient of gas, C g , can be related to the storage coefficient of water by, Bulk modulus of water and gas. The bulk modulus of water, K w , can be calculated as follows The bulk modulus of gas, K g , can be calculated as follows From the ideal gas law, From Boyle's Law for isothermal conditions, substituting equation (32) back into equation (30) K Constitutive relations of a damage model Failure results in a decrease in the elastic modulus and an increase in the intrinsic permeability, as expressed in this paper by equations (34) and (35), respectively, where k ij is the permeability at any given time (m 2 ) k UD is the permeability before damage (m 2 ) and is expressed by equation (21) k D is the increase in permeability as a result of damage (m 2 ) and can be expressed by equation (36) where k max (m 2 ) is the maximum permeability obtained when D = 1 c is a smoothing coefficient (unitless) and may be calibrated anywhere from 2 to 20 depending on the mesh size and the rate of increase in damage k max and c, are determined through model calibration.
The criterion for D t in multi-axial tension is expressed by equation (37), while the criterion for D c in multi-axial compression is expressed by equation (38) where ε is the strain (unitless) ε tu is the tensile strain corresponding to a complete material failure ε t0 is the strain corresponding to the point of tensile strength ε c0 is the strain corresponding to the compressive strength f tr is the residual tensile strength (Pa) f cr is the residual compressive strength (Pa).
Note for equations (37) and (38) the following sign convention is used: (1) strain components and volumetric deformations are positive for expansion and negative for contraction; and (2) stress components are positive for tension and negative for compression.

Governing equations
Conservation of momentum (quasistatic equilibrium). The total equilibrium equation for a cubical soil element can be expressed in tensor notation by equation (39) where, σ ij is the stress tensor (Pa) F v,i is the volumetric body force tensor (kg m −2 s −2 ) ∂σ ij /∂x j represents the change in normal and shear stresses across the soil element (kg m −2 s −2 ).
Substituting equation (8) and differentiating, the governing equation for the conservation of momentum using a linear poro-elasticity model is expressed by equation (40).
Conservation of water mass. For a soil element depicted in Figure 2, the general equation for the conservation of water mass in porous media can be expressed by equation (41).
where, (∂(ρ w nS))/∂t is the rate of change of mass of water in the soil element over time (kg m −3 s −1 ) Substituting in Darcy's Law (18), and solving the right-hand side of equation (41) results in the following governing equations for the conservation of water mass, where K w is the bulk modulus of water (Pa). It should be noted that the RHS of equation (42) is derived from the differentiation of ∂(ρ w nS w )∂t.
The final term in equation (42) is multiplied by (1−n) which is derived from differentiating the equation for porosity, Using the Quotient Rule to differentiate, where, dV/V = dε kk , the volumetric strain, and where it is assumed that dVs/V ≪ dε kk (1−n), Conservation of gas mass. For a soil element depicted in Figure 2, the general equation for the conservation of gas mass in porous media can be expressed by equation (46), respectively where, ρ g , is the gas density (kg gas m −3 ) ) is the net advective flux of gas for the soil element (kg m −3 s −1 ) ∂/∂x i (−D e (∂/∂x k )(ρ g n(HS w ))) is the net diffusive flux of gas for the soil element (kg m −3 s −1 ) ∂(ρ g n(1−S w + HS))/∂t is the rate of change of mass of gas in the soil element over time (kg m −3 s −1 ) H is Henry's coefficient (kg species A m −3 in aqueous phase kg −1 species A m 3 in gas phase). Substituting in Darcy's Law (19) and solving the right-hand side of equation (46), results in the following governing equation for the conservation of gas mass where K g is the bulk modulus of the gas phase (Pa).

Numerical model description and modelling approach
Overview of the numerical model A 2D-axisymmetric and 3D time-dependent (i.e. transient) HM-coupled multiphysics models were developed. The models simulate the simultaneous migration of gas and liquid (two-phase flow) in porous media, which are coupled to the linear poro-elastic behaviour of the solid matrix using the FEM. The commercially available code COMSOL Multiphysics ® was used to numerically solve the governing equations of the model. The numerical model was based on an experimental setup conducted by the BGS to assess the 1D flow through a saturated bentonite sample under a constant volume boundary stress condition (Daniels & Harrington 2017). In this experiment, a confined cylindrical sample of near-saturated bentonite was injected on one side with helium gas with increasing gas pressures over a period of 120 days. The other side was left at a constant water backpressure throughout the duration of the experiment. The experiment was conducted under isothermal conditions at a temperature of 293.15 K. During the experiment, a number of parameters where measured including the gas inflow and outflow, the porefluid pressure at defined porefluid arrays, and the total radial stresses at radial load cell arrays. The cylindrical specimen of MX-80 bentonite had a diameter of 60 mm and a length of 120 mm. Table 2 provides the location of the monitoring sensors within the specimen. The BGS provided the experimental data that was used in the model validation (Daniels & Harrington 2017).

Modelling approach
The modelling approach included the simulation of a number of study scenarios. A validation study scenario S1 was run with the HM poro-elastic model for both 2D-axisymmetric and 3D models. Several parametric studies were performed to assess the contribution of advection, diffusion and linear elastic deformation on flow behaviour (S2-S4). Similarly, sensitivity analyses were performed to identify the influence of changing the initial porosity, AEV and intrinsic permeability on flow behaviour (S5-S10). Scenario S11 was run to determine whether better agreement with the experimental results could be obtained by inclusion of a damage model. Lastly, the effect of mesh size on the numerical solution was investigated by running S1 using a number of mesh sizes (S12). Table 3 summarizes the simulations which were run.

Model geometry and mesh
The 2D-axisymmetric and 3D geometry and mesh sizes are presented in Figures 3 and 4, respectively. Due to the nature of axisymmetric models (i.e. rotates along a vertical axis), the 2D-axisymmetric  model was configured vertically. As a result, the effects of gravity have been turned off for this specific configuration in COMSOL ® , as it differs from the actual experimental set-up. The number of elements and degrees of freedom for each geometry are provided in Table 4.

Material properties
Material properties for the solid bentonite MX-80 soil matrix, helium gas and water are provided in Table 5. Material properties were provided by BGS for the initial porosity, intrinsic permeability,  Young's Modulus and Poisson ratio. The dry density of bentonite was calculated from the dry weight and volume of the bentonite sample. The volume-surface mean diameter, D vs , was calculated from equation (22), and the initial tortuosity, τ, was calculated from equation (26).

BCs
The hydraulic and mechanical BCs for gas transport, water transport and momentum transport are provided below.
Gas conservation BCs. For gas transport, a no flow Neumann BC, ∂p g /∂x i = 0, was set at the radial boundaries. Dirichlet BCs were set at the gas injection pressure as a function of time for the lower boundary. The upper boundary was set at atmospheric pressure to allow for gas to flow out of the sample. For the concentration of gas in porewater, C g,H2O , Dirichlet BCs were applied to the lower boundary as a function of the gas injection pressure and to the upper boundary based on an atmospheric pressure condition as follows; from the Ideal Gas Law Re-arranging (50) where p g is the poregas pressure (Pa). Assuming instantaneous dissolution, the concentration of gas in the porewater can be calculated by multiplying equation (50) by Henry's coefficient, H, and the portion of water in a unit cell, nS w C g, H2O = r g H(nS w ) (51) where C g,H2O is the concentration of gas in water (kg gas m −3 water)  H is Henry's coefficient (kg species A m −3 in aqueous phase kg −1 species A m 3 in gas phase).
Water conservation BCs. For water transport, a no flow Neumann BC, ∂p w /∂x i = 0, was set for the lower and radial boundaries. A Dirichlet BC was set as the water backpressure for the upper boundary.
Momentum conservation BCs. For the moment conservation equation, a roller constraint was applied along the upper, lower and radial boundaries to represent a condition whereby the boundary is free to move in the tangential direction, but is fixed in the normal direction, simulating a constant volume condition, as per the experimental set-up. The gas injection pressure and water backpressure Dirichlet BCs were imported from the experimental data provided by BGS and have been plotted in Figure 5(a), while the concentration of dissolved gas BCs have been plotted in Figure 5(b).

Initial value conditions
The initial conditions at t = 0 s across the domain are provided in Table 6, and are based on experimental It should be noted that the bentonite sample was not fully saturated upon the start of the experiment but underwent a hydration phase which began on t = 7.3 days and ended on t = 39 days. Thereby an initial degree of saturation, S w,inital , measured at 0.96 was used as the initial condition. This is important as the model assumes a continuous gas-phase  Initial poregas pressure p g initial 1.01E + 05 Pa Initial degree of saturation S w initial 0.96 Initial suction (from SWCC) s initial 5.8E + 06 Pa Initial porewater pressure p w initial −5.7E + 06 Pa Initial displacement field u i 0 m Initial stress σ 0 xx = σ 0 yy = σ 0 zz 5.0E + 05 Pa Initial gas concentration in porewater @STP C g,H2O 0.073 mol m −3 throughout the specimen and does not account for entrapment of gas within the pore space.
Intrinsic permeability. AEV, effective diffusivity, chi parameter, SWCCs and relative permeability functions The intrinsic permeability, AEV and effective diffusivity as a function of porosity are provided in Figure 6. Khalili and Khabbaz's χ-value is displayed as a function of suction and AEV in Figure 7. Table 7 lists the fitting parameters for the van Genuchten-Mualem model SWCC and relative permeability function for the S1 validation study. Values for the van Genuchten SWCC fitting parameters were assessed from literature taking the calculated dry density of MX-80 bentonite (Villar 2004;Villar & Lloret 2008;Man & Martino 2009;Tripathy et al. 2014). Figures 8 and 9 depict the SWCC in the form of the effective degree of saturation and relative permeability of gas and water through bentonite as a function of suction.

Results and discussion
In this section, the following sign convention is used: (1) strain components and volumetric deformations are positive for expansion and negative for contraction; and (2) stress components and pressures are positive for compression and negative for tension.

S1: Validation study
The results of the validation study are presented below. As shown in Figure 10, the proposed HM model was capable of simulating two-phase flow through the MX-80 bentonite specimen. Gas breakthrough (i.e. into the sample) was observed at 61   (12) and (18) where AEV = 8.0E 6Pa and s′ aev = 3.46E7 Pa. † A value of L′ = 1 was inputted by the user and required for COMSOL ® to converge, and does not significantly change the relative permeability curves.
days when the gas-injection pressure reached 8 MPa.
Over the 120-day duration of the experiment gas migrated as a gas front through the first 5% of the sample. The 2D-axisymmetric and 3D simulation results were compared against (1) the experimental inflow and outflow profiles (Fig. 11), (2) the experimental porefluid pressure profiles at the radial porewater arrays (Fig. 12), and (3) the evolution of total stresses at the radial load cell arrays (Fig. 13) over the 120-day study period.
Total gas flow into the sample is depicted in Figure 11(a). From the figure it is clearly observed that the onset of gas flow into the sample (i.e. gas breakthrough) occurs when the injection pressure reaches the AEV (around 61 days). The AEV of 8 MPa was the only model parameter calibrated to match the experimental results, with the corresponding van Genuchten SWCC fitting parameter a′ calculated from equation (12). The model was also able to reproduce the point in time when the gas injection pump was stopped, resulting in a rapid decline in  gas inflow (71 days). The shape and timing of the modelled gas inflow curves matches the experimental data reasonably well, however, the peak magnitude of the modelled gas inflow was approximately half the magnitude of the experimental results. The fact that the simulation results of gas inflow were significantly lower than the experimental results, may be due to a lower modelled intrinsic permeability at the breakthrough pressure. As the intrinsic permeability applied to the model was based on those measured by the BGS, it may be that some damage occurs at the breakthrough poregas pressure, which allow for a significant increase in poregas pressure. The 2D and 3D models were able to reproduce similar results to each other. Differences are likely due to mesh size and will be discussed further in scenario S12. Figure 11(b) shows the experimental and modelled gas outflow. The experimental gas outflow results align with those of the inflow results, providing evidence that dilation-controlled gas flow is occurring, as the majority of gas moving into the system also leaves the system, with very little gas being stored within the system. Using the proposed model, no gas outflow was observed from the specimen, indicating that the model was not able to reproduce dilation-controlled gas flow, as could be observed via the peaks of the experimental dataset. This further supports the reasoning that inclusion of a damage model may be required to simulate the increase in gas permeability within the specimen in order to obtain gas breakthrough flow in-line with the experimental results.
Similar results, though not all, are observed when comparing the porefluid pressure profiles obtained between the experimental and modelled results, as depicted in Figure 12(a). Following hydration of the sample which occurred between day 7.3 and day 39, gas testing began. At day 46 there is a sudden decrease in porefluid pressure observed in the experimental results when the gas pressure into the system is ramped up. This is followed by an increase in porewater pressure in Porewater Array 1 and 2, while Porewater Array 3 continues to decrease. This initial decrease in porewater pressure could have occurred due to several factors. Supporting  the notion that some damage or dilation may be occurring, an increase in the pore-space due to damage would result in a decrease in porewater pressure. Additionally, this could be further explained by the effect of suction, where the system attempts to maintain a capillary pressure close to zero. Unfortunately, these features were not able to be reproduced by the model and will be further investigated in scenario S11.
Once the gas injection pressure reaches the AEV however, increases in porefluid pressures at the porefluid arrays were observed in response to increasing gas injection pressures. At this point the model was able to adequately reproduce key features of the experimental results including both the shape and timing of changes in porefluid pressure. However, the extent of porefluid pressure build-up is not as large or as sharp as that depicted by the experimental results. Regarding the lower peak of the porefluid pressure profiles, this is to be expected as the model was not able to simulate gas flow through the porefluid pressure arrays, and therefore there is little contribution of poregas pressure to the porefluid pressure. As a result, the increase in porefluid pressure observed is dominated by the displacement of porewater by the poregas as it migrates through the front end of the bentonite specimen, creating a porewater pressure gradient.
This leads to another distinguishable feature; in the modelled results, following gas breakthrough into the sample, there is a clear separation in the porefluid pressure profiles between the three porewater arrays. This lag can be attributed to the time it takes for the displaced water to push through the soil column, and results in the increase in porefluid pressures. This can be seen in the experimental results as well, although less pronounced due to the use of a log scale in Figure 12(a).
The difference in the sharpness or steepness of the porefluid pressure profiles in comparison to the experimental results is also expected and can be partially attributed to the use of van Genuchten SWCC equation, which provides a continuous smooth function at the interface between the saturation zone and the desaturation zone (i.e. AEV). In addition to this, during analysis of the modelling results obtained using the 2D-axisymetric model v. 3D model, it was identified that mesh size plays a large role in the shape and steepness of the porefluid pressure profile at the AEV. This is further investigated in scenario S12. Fig. 13. S1 results -Total radial stresses over time.
In Figure 12(b), the results of the porefluid pressure profiles along the axial length of the soil specimen are presented. At the beginning of the simulation (t = 0 days), the modelled results show a negative porefluid pressure along the soil specimen. This is expected as the soil specimen was not fully saturated. During the hydration phase (day 7.3 to day 39), the porefluid pressures stabilize at 1 MPa, representing equilibrium with the defined BCs. Following a ramp up of gas injection pressures front from day 39 to day 63, a small increase in porefluid pressures at the gas injection front is observed (at t = 60 days). Once the gas injection pressure exceeds the AEV at day 63, a continuous gas phase begins entering the sample, and a large increase in the porefluid pressure gradient across  the sample is observed as depicted by the curve t = 80 days. Once gas injection is stopped, the porefluid pressure gradients begin to decrease and equilibrate back towards the backpressure BC of 1 MPa as depicted by the decrease in porefluid pressure gradients along the length of the soil specimen at day 100 and day 120. The behaviour depicting the change in porefluid pressure gradients over time depicted by the model results in Figure 12(b) agree with our conceptualization of the experiment.
For the 3D case, the total radial stresses, which are a measure of both the net normal stress acting on the solid matrix and porefluid pressure is presented in Figure 13. Again, we see the model able to reproduce both the shape and timing of the stress evolution acting on the radial load cell arrays. However, the magnitude could not be reproduced and was an order lower. The increase in the total radial stresses, which leads to the deviation between the experimental results and the modelled results, occurs during the initial hydration phase. This is, in part, likely a result of neglecting the effect of swelling pressure in the numerical model, a key behaviour of expansive clays. If a swelling pressure is considered in the model, then an increase in the total radial stress response being modelled during the hydration phase would be expected as well and could account for this difference.
From the investigation of the validation scenario, S1, it was identified that improved results could be obtained if at some critical pressure, a large change in porosity, resulting in an increase in the gas permeability of between two to three orders of magnitude was introduced into the model as will be demonstrated by S11.

S2-S4: Parametric studies
The 3D results of the parametric studies S2, S3 and S4 are provided in Figure 14. These results assess the contribution of advection, diffusion of gas within porewater and the advection of dissolved gas in porewater to the flow behaviour, given the benchmark dataset available.
In Figure 14(a), the total gas inflow is broken down to the contributions from the individual mechanisms of gas migration being modelled: (1) advective gas flow; (2) diffusive gas flow (diffusion of gas dissolved in porewater); and (3) advection of dissolved gas within porewater. As can be seen, gas migration is diffusion dominated until the AEV is reached, and gas migration into the sample occurs, whereby gas advection becomes the predominant transport mechanism. The results also show that advection of dissolved gas in porewater is minimal, which is to be expected as there is little movement of porewater through the system.
To demonstrate that diffusion is being adequately presented by the model, Figure 14(b) depicts the gas concentration along the length of the specimen over time. There are four distinct phases identified in the figure, and each is separated by a jump corresponding to periods when the gas injection pressure is rapidly increased requiring the system to re-establish equilibrium. The first section represents the initiation of the experiment at t = 0 days. The second phase    represents the hydration period (day 7.3 to day 39), where the gas injection pressure is maintained with the water backpressure at 1 MPa. The third phase begins at the onset of gas testing on day 39, and as gas injection pressure is ramped up, so does the concentration of gas in porewater resulting in a new gradient and the system re-establishing equilibrium. The fourth phase can be identified by a decreasing concentration gradient which occurs once gas injection has stopped (day 73). These are the expected behaviours that would be observed based on the governing equations being applied in the model.
In Figure 15(a) a comparison of the porefluid profile curves for the S2, S3 and S4 scenarios are indistinguishable. This is to be expected as the porefluid pressure is a function of the poregas and porewater pressures via equation (3), and the porewater pressure is not dependent on the concentration of dissolved gas.
The results of S4 are compared to S1 in Figure 15 (b) to see the effect of linear-elastic deformation to the flow behaviour. Although there is little contribution to flow behaviour, the HM model has smoother curves than the hydraulic models alone and could be an indication of a more stable model.
The results of the parametric studies demonstrate that for the experiment under consideration, at high gas pressures, gas flow is controlled by advection, with diffusion and elastic deformation contributing little to flow behaviour. It is also evident by the porefluid profile curves, that the mechanical deformation resulting from the inclusion of a linear-elastic model is not sufficient to result in dilatancy and selfhealing of the bentonite pores. This is further illustrated by Figure 16, which shows small changes in porosity over time at each porefluid array. It should be noted, that as there is a constant volume BC applied, the total deformation throughout the sample is in equilibrium.
Sensitivity analysis S5-S6: Effect of porosity. The effect of increasing and decreasing the initial porosity on the porefluid pressure profiles are depicted in Figure 17(a) and on gas inflow in Figure 17(b). As described by equation (17) and equation (21), both the AEV (and corresponding SWCCs and relative permeability curves) and intrinsic permeability are a function of the porosity. An increase in the initial porosity to 0.48, results in both a decrease in the AEV, noted by earlier gas breakthrough, and an increase in the intrinsic permeability, as depicted by higher peak porefluid pressures. Similarly, by decreasing the initial porosity to 0.3, the required AEV for gas breakthrough is not reached by the experimental gas injection pressures. As a result, there is no gas migration into the specimen, and no change in the porefluid pressure profiles. This is as expected, providing confidence in the model.  Figure 18(b). Since the amount of elastic deformation is small (see Fig. 16), and therefore contributes little to changes in the AEV and intrinsic permeability, results of S7 and S8 are like those of S5 and S6. A decrease of the AEV to 5.0E + 06 Pa results in earlier gas breakthrough and higher peak porefluid pressures, and an increase of the AEV to 2.0E + 07 Pa, results is very little migration into the specimen.
S9-S10: Effect of intrinsic permeability. The effect of increasing and decreasing the intrinsic permeability, k ij , independent of porosity, on the porefluid pressure profiles are depicted in Figure 19(a) and on gas inflow in Figure 19(b). For the S9 simulation an increase in k ij by two orders of magnitude to 10 −19 m 2 shows an increase in the porefluid pressure profile once the AEV is reached, however the peak of the curve is now lower than the S1 scenario. This is to be expected as there will be a lower porefluid pressure gradient.
For the S10 scenario, a decrease of k ij by over an order of magnitude to 10 −22 m 2 resulted in no visible increase in the porefluid pressure profile, even when the AEV is exceeded. This is to be expected, as the experimental intrinsic permeability of 3.4E −21 m 2 is already low, and by decreasing it further, advective gas migration becomes so slow that once the gas pressure exceeds the AEV, there is little displacement of porewater within the experimental timeframe. S11: Effect of damage. The parameters for the damage model presented by equations (34) to (38) are provided in Table 8. The compressive and tensile strengths were estimated based on information provided by Man & Martino (2009), while the strains at compressive and tensile failure were calculated from the elastic modulus.
The results of Scenario S11, investigating the inclusion of a damage model, are presented below. Figures 20-22 show the results for the effect of damage on the gas inflow and outflow behaviour, porefluid pressure profile behaviour and total radial stresses behaviour, respectively. Overall, the modelled results demonstrate much better agreement with the experimental results. Of interest is the maximum intrinsic permeability, k max , which was calibrated to a value of 1E-18 m 2 to match the peak experimental inflow.
As depicted in Figure 20(a), when considering damage, the modelled inflow was able to be reproduced quite well, however as shown in Figure 20 (b), the model was still not able to reproduce gas outflow even with the inclusion of damage. If k max was set higher than it would be possible to obtain gas outflow, however, this would result in overshooting the modelled gas inflow well above the experimental inflow. In Figure 21, the magnitude, shape and timing of the porefluid pressure profile curves are closer to the experimental results, however the full magnitude was still not realized. The total radial stress profiles are depicted in Figure 22 and have only improved slightly as a result of damage. Figure 23 depicts the point of gas entry into the specimen and the migration of gas over time. Significantly increased gas flow is observed with the inclusion of the damage model, however still no preferential flow pathways representing dilationcontrolled gas flow were observed. Instead there is a moving gas front along the soil specimen, with a minimum modelled S w of approximately 0.9 obtained at the point of gas injection and 0.94 at the gas front (approximately 20 mm from the injection point).
Although inclusion of damage to the model provides much better agreement with the experimental results, key features of the experimental gas outflow results depicting dilatancy and self-healing were not reproduced. There is, therefore, still a need to investigate additional mechanisms to model dilatancy-controlled gas flow. Such mechanisms will be investigated in future studies. These will include the consideration of spatially weighted heterogeneity for porosity and a poro-elasto-plastic model to simulate flow through preferential pathways. Inclusion of a swelling stress term will also be added to investigate the swelling behaviour in  expansive soils. Finally, an investigation into constitutive relations to simulate the self-healing behaviour of swelling soils will be included. S12: Effect of mesh size. A major finding of the study was the effect of mesh size on the shape of the porefluid pressure profile curves and the gas inflow and outflow curves. This was initially observed when moving from a 2D-axisymetric geometry to a 3D geometry, where the number of degrees of freedom needed to solve the model significantly increases. As a result of the increased computational expense, the 3D model was not able to complete a run at higher mesh sizes, using the fully coupled, or segregated solvers in COMSOL ® . To exemplify this, the number of elements and degrees of freedom for each geometry are provided in Table 9. Figure 24 shows the effect that mesh size has on the shape of the porefluid pressure curve. In Figure 24 results of the porefluid pressure profiles for Array 1 (blue curves) and Array 2 (red curves) with an extremely fine, finer and normal mesh are compared. By applying a more refined mesh (i.e. solving more degrees of freedom), results in a steeper porefluid pressure curve at gas inflow and a higher peak porefluid pressure, leading to a higher resolution solution.
A similar comparison can be made between the 2D-axisymmetric solution with an extremely fine mesh size and the 3D solution with a finer mesh size as depicted in Figure 12(a). Due to the computational cost of modelling in 3D this poses a number of challenges. For one, preferential flow pathways cannot be adequately modelled using a 2D-axisymmetric geometry, as the microfractures cannot intersect the axis of symmetry. Therefore, in  order to model dilatancy-controlled gas flow in heterogeneous specimens and propagate preferential flow pathways, it will be crucial to model in 3D.
To address this current limitation, future studies will assess the use of different solvers and methods in COMSOL ® , which may include the incorporation of discrete models, and identify ways to optimize the model for improved computational efficiency.

Conclusions
An important component in the design and long-term safety assessment of a DGR is the long-term performance of bentonite seals as barriers against gas migration in order to ensure the future safety to human health and the environment as a result of the exposure to these radionuclides. This paper proposes a hydro-mechanical linear poro-elastic visco-capillary mathematical model for advectivediffusive controlled two-phase flow using a Bishop's effective stress. The objective of this study was to gain insight into the key parameters influencing gasflow behaviour in multiphase flow and to identify whether dilation-controlled gas flow could be represented by such a model. A validation study was conducted against experimental data from 1D flow through a saturated bentonite sample under a constant volume boundary stress condition. The experimental data demonstrated features of dilatancy-controlled gas flow. The modelled results were compared to key features of the experimental data, including the gas inflow and outflow behaviour, pore-pressure fluid behaviour and evolution of radial stresses over the duration of the experiment. Although the model was not able to reproduce dilation-controlled gas outflow from the sample, it was able to simulate two-phase flow of water and gas through the sample, once the AEV was reached, and gas flow into the specimen was observed.
The model was able to reproduce key features of the experimental results, including both the shape and timing of changes in gas inflow, porefluid pressure and radial stress evolution through the sample. However, the magnitude of such features was not able to be reproduced. It was identified that improved results could be obtained if at some critical pressure, a large change in porosity, resulting in an increase in the gas permeability of two to three orders of magnitude was introduced into the model.
In light of this, the effect of the inclusion of a damage model was investigated. The study confirmed that a damage mechanism calibrated to the maximum gas permeability required to match the peak experimental inflow could be successfully applied and results in much better agreement with the experimental results.
Parametric studies were performed to assess the contribution of advection, dissolution and diffusion of gas, as well as linear-elastic deformation on flow behaviour. The results of the parametric studies demonstrated that at high gas pressures, gas flow was controlled by advection, with diffusion and elastic deformation contributing little to flow behaviour. Sensitivity analysis were conducted by changing the initial porosity, permeability and AEV of the swelling soil to test the model behaviour. In all cases the model behaved as expected.
This study provides a preliminary model for migration of gas through low-permeability swelling soils and was able to simulate two-phase flow with dissolution of gas. Further development of the model to better mimic dilation-controlled gas-flow will be addressed in subsequent studies. This model provides a basis for developing more advanced models to describe the phenomena of dilatancycontrolled gas flow.
The results of this study conclude that in order to mimic dilation, and dilatancy-controlled gas flow, additional mechanisms need to be considered within the model. Such mechanisms will be investigated in future studies. These will include the consideration of spatially weighted heterogeneity for porosity and a poro-elasto-plastic model to simulate flow through preferential pathways. Inclusion of a swelling stress term will also be added to investigate the swelling behaviour in expansive soils. Finally, an investigation into constitutive relations to simulate the self-healing behaviour of swelling soils will be addressed.