Development of a mathematical model for gas migration (twophase flow) in natural and engineered barriers for radioactive waste disposal
Correspondence: jinfante@uottawa.ca

Published:January 01, 2019
 Standard View
 Chapter PDF

CiteCitation
E. E. Dagher, T. S. Nguyen, J. A. Infante Sedano, 2019. "Development of a mathematical model for gas migration (twophase flow) in natural and engineered barriers for radioactive waste disposal", Multiple Roles of Clays in Radioactive Waste Confinement, S. Norris, E.A.C. Neeft, M. Van Geet
Download citation file:
 Share

Tools
Abstract
In a deep geological repository (DGR) for the longterm 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 bentonitebased materials, are currently the preferred choice of seal materials. Understanding the longterm performance of these seals as barriers against gas migration is an important component in the design and longterm safety assessment of a DGR. This study proposes a hydromechanical linear poroelastic viscocapillary mathematical model for advectivediffusive controlled twophase flow through a lowpermeability 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 lowpermeability 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 (LowLevel 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. Todate, this waste, in the form of low, intermediate and highlevel radioactive waste has been primarily stored onsite at nuclear power plants or in belowsurface radioactive waste management facilities.
The LowLevel 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 longterm management in Canada, the bulk of which is lowlevel radioactive waste. Table 1 provides waste inventory estimates by the end of 2011 and by the end of 2050.
Waste category  Waste inventory to the end of 2011  Waste inventory to the end of 2050 

Highlevel radioactive waste (nuclear fuel waste)  9400 m^{3}  20 000 m^{3} 
Intermediatelevel radioactive waste  33 400 m^{3}  67 000 m^{3} 
Lowlevel radioactive waste  2 343 000 m^{3}  2 594 000 m^{3} 
Waste category  Waste inventory to the end of 2011  Waste inventory to the end of 2050 

Highlevel radioactive waste (nuclear fuel waste)  9400 m^{3}  20 000 m^{3} 
Intermediatelevel radioactive waste  33 400 m^{3}  67 000 m^{3} 
Lowlevel radioactive waste  2 343 000 m^{3}  2 594 000 m^{3} 
In identifying a longterm 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 longterm 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 (Horseman et al. 1999; 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 buildup 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 longterm performance of these seals as barriers against gas migration is an important component in the design and longterm safety assessment of a DGR.
An investigation of the gas transport processes in lowpermeability 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 hydromechanical 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 advectivediffusive 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 viscocapillary twophase 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 poresize distribution and soil compressibility. The airentry value (AEV) or gasentry 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 twophase 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 fieldscale 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 pressureinduced 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 (Harrington & Horseman 1999; Horseman et al. 1999, 2003; Graham et al. 2012; Harrington et al. 2012, 2017; Cuss et al. 2014; Sellin 2014; Bennett et al. 2015).
Marschall et al. (2005) described this primary transport mechanism as dilatancycontrolled gas flow or ‘pathway dilation’ and identified it as an important transport mechanism in clay soil or clayrich rock with low tensilestrength. 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 stressstate and state of deformation of the soil.
Finally, gas transport along macroscopic fractures, also known as hydro and gasfracturing, 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 lowtensile strength material when gas pressure buildup 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 inter and intrastack 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 gasfilled 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 microcracks (Fall et al. 2014). Their model was limited in part because they did not consider the elastoplastic 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 elastoplastic 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 415 – entitled, 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 twophase flow, particularly dilatancycontrolled gas flow through lowpermeability swelling soils over the lifetime 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 longterm safety assessment of a DGR.
Stemming from these international initiatives to understand gas migration and multiphase 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 (DECOVALEX2019) was established. This task, led by the British Geological Survey (BGS), further attempts to identify the physical hydromechanical (HM) mechanisms required to adequately model dilatancycontrolled gas migration. This research is in part a contribution to Task A of DECOVALEX2019.
This paper presents a fullycoupled, hydromechanical linearelastic mathematical model for advectivediffusive viscocapillary 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 1dimensional 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 lowpermeable 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 airentry 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 stressstate 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 poroelasticity
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 BiotWillis 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)
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 saturation, S_{w,r}
From the fitted SWCC, the model can be used to predict the degree of saturation of water, S_{w}, from measuring the stressstate 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_{aev} – e_{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 rearranging equation (16) we can express this relationship in terms of porosity,
Darcy’s Law for twophase 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 particlesize. The KozenyCarman 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 KozenyCarmen Equation is given by,
where S_{0} is the specific surface (surface exposed to fluid per unit solid volume) (cm^{2} cm^{−3}).
For finegrained expansive clays, which have very small grainsizes less than 2 µm, and are nonuniform, the KozenyCarman equation may not be appropriate (Pall & Moshenin 1980; ValdesParada 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 volumesurface mean diameter (cm).
If the porosity and intrinsic permeability were known, the volumesurface 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 volumesurface 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},
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) by,
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)
Constitutive relations of a damage model
A damage model has been introduced based on that developed by Tang et al. (2002) and adopted by Fall et al. (2014), whereby an element under tension or compression will experience an elastic damage, D (0 < D < 1), when the stress of the element satisfies the strength criterion and the element begins to fail. 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
D is the damage (unitless)
E_{0} is the Young’s modulus before damage (Pa)
E is the Young’s modulus after damage (Pa)
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 multiaxial tension is expressed by equation (37), while the criterion for D_{c} in multiaxial compression is expressed by equation (38)
where
ε is the strain (unitless)
ε_{tu} is the tensile strain corresponding to a complete material failure
ε_{t}_{0} is the strain corresponding to the point of tensile strength
ε_{c}_{0} 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 poroelasticity 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})
∂/∂x_{i}(ρ_{w}n(v_{iw}−v_{is})) is the net advective flux of water for the soil element (kg m^{−3} s^{−1}).
Substituting in Darcy’s Law (18), and solving the righthand 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,
where V_{v} is the volume of voids (m^{−3}).
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})
∂/∂x_{i}(ρ_{g}n((v_{ig}−v_{is}) + H(v_{iw} − v_{is}))) 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 righthand 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 2Daxisymmetric and 3D timedependent (i.e. transient) HMcoupled multiphysics models were developed. The models simulate the simultaneous migration of gas and liquid (twophase flow) in porous media, which are coupled to the linear poroelastic 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 nearsaturated 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 MX80 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).
Sensor name  Measurement type  Axial distance from injection face (mm) 

Injection load cell  Total stress  0 
Radial load cell 1  Total stress  15.2 
Radial load cell 2  Total stress  60 
Radial load cell 3  Total stress  104.8 
Backpressure load cell  Total stress  120 
Radial porewater array 1  Porefluid pressure  38.6 
Radial porewater array 2  Porefluid pressure  60 
Radial porewater array 3  Porefluid pressure  81.4 
Central filter (Middle)  Porefluid pressure  60 
Sensor name  Measurement type  Axial distance from injection face (mm) 

Injection load cell  Total stress  0 
Radial load cell 1  Total stress  15.2 
Radial load cell 2  Total stress  60 
Radial load cell 3  Total stress  104.8 
Backpressure load cell  Total stress  120 
Radial porewater array 1  Porefluid pressure  38.6 
Radial porewater array 2  Porefluid pressure  60 
Radial porewater array 3  Porefluid pressure  81.4 
Central filter (Middle)  Porefluid pressure  60 
Modelling approach
The modelling approach included the simulation of a number of study scenarios. A validation study scenario S1 was run with the HM poroelastic model for both 2Daxisymmetric 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.
Scenario number  Model  Description  Porosity  Intrinsic permeability (m^{2})  AEV (Pa)  SWCC fitting parameters  

a′ (1/m)  n′  m′  
S1Validation  HM  Advection, dissolution and diffusion linear poroelasticity  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Parametric studies  
S2  H  Advection, no dissolution, no diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S3  H  Advection, dissolution, no diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S4  H  Advection, dissolution, diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Sensitivity analysis  
S5  HM  ↑ por  0.48  5.1E−21  6.2E + 06  0.0008  2.04  0.51 
S6  HM  ↓ por  0.30  6.9E−22  1.6E + 07  0.0003  2.04  0.51 
S7  HM  ↑ AEV  0.44  3.4E−21  2.0E + 07  0.00025  2.04  0.51 
S8  HM  ↓ AEV  0.44  3.4E−21  5.0E + 06  0.0006  2.04  0.51 
S9  HM  ↑ k_{ij}  0.44  1.0E−19  8.0E + 06  0.0006  2.04  0.51 
S10  HM  ↓ k_{ij}  0.44  1.0E−22  8.0E + 06  0.0006  2.04  0.51 
S11  HM  Inclusion of a damage model  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S12  HM  Effect of mesh size  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Scenario number  Model  Description  Porosity  Intrinsic permeability (m^{2})  AEV (Pa)  SWCC fitting parameters  

a′ (1/m)  n′  m′  
S1Validation  HM  Advection, dissolution and diffusion linear poroelasticity  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Parametric studies  
S2  H  Advection, no dissolution, no diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S3  H  Advection, dissolution, no diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S4  H  Advection, dissolution, diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Sensitivity analysis  
S5  HM  ↑ por  0.48  5.1E−21  6.2E + 06  0.0008  2.04  0.51 
S6  HM  ↓ por  0.30  6.9E−22  1.6E + 07  0.0003  2.04  0.51 
S7  HM  ↑ AEV  0.44  3.4E−21  2.0E + 07  0.00025  2.04  0.51 
S8  HM  ↓ AEV  0.44  3.4E−21  5.0E + 06  0.0006  2.04  0.51 
S9  HM  ↑ k_{ij}  0.44  1.0E−19  8.0E + 06  0.0006  2.04  0.51 
S10  HM  ↓ k_{ij}  0.44  1.0E−22  8.0E + 06  0.0006  2.04  0.51 
S11  HM  Inclusion of a damage model  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S12  HM  Effect of mesh size  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Model geometry and mesh
The 2Daxisymmetric 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 2Daxisymmetric 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 setup. The number of elements and degrees of freedom for each geometry are provided in Table 4.
Model geometry  Mesh qualifier  No. of elements  Degrees of freedom 

2D axisymmetric  Extremely fine  8886  117 132 
3D  Normal  10 698  161 741 
Model geometry  Mesh qualifier  No. of elements  Degrees of freedom 

2D axisymmetric  Extremely fine  8886  117 132 
3D  Normal  10 698  161 741 
Material properties
Material properties for the solid bentonite MX80 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 volumesurface mean diameter, D_{vs}, was calculated from equation (22), and the initial tortuosity, τ, was calculated from equation (26).
Material  Parameter  Value 

Solid soil matrix (Bentonite)  Saturated/intrinsic permeability, k_{ij} (m^{2})  3.40E−21 
Dry density, ρ_{d} (kg m^{−3})  1560  
Maximum degree of saturation, S_{ws}  1  
Residual degree of Saturation, S_{wr}  0.05  
Porosity, n  4.40E−01  
Young’s Modulus, E (MPa)  307  
Poisson ratio, ν  0.4  
Bulk Modulus of bentonite, K_{s} (Pa)  2.0E + 08  
Compression Index, 1/K_{s}(Pa^{−1})  5.00E−09  
Volumesurface mean diameter, D_{vs} (m) calculated based on equation (23)  1.50E−09  
Initial tortuosity, τ  0.69  
Gas (Helium)  Density, ρ_{g} (kg m^{−3})  $pgMRT$ 
Dynamic Viscosity @ 293.15 K (Pa s)  2.0E−05  
Bulk Modulus of helium, K_{g} (Pa)  p_{g}  
Henry’s Coefficient, H  0.0091  
Diffusivity of helium in water (m^{2} s^{−1})  6.29E−09  
Diameter of helium gas particle, d_{he} (m)  1.0E−10  
Liquid (Water)  Density, ρ_{w} (kg m^{−3})  1000 
Dynamic Viscosity @ 293.15 K (Pa s)  0.001  
Bulk Modulus of water, K_{w} (Pa)  2.2E + 09  
Compression Index, 1/K_{w} (Pa^{−1})  4.55E−10 
Material  Parameter  Value 

Solid soil matrix (Bentonite)  Saturated/intrinsic permeability, k_{ij} (m^{2})  3.40E−21 
Dry density, ρ_{d} (kg m^{−3})  1560  
Maximum degree of saturation, S_{ws}  1  
Residual degree of Saturation, S_{wr}  0.05  
Porosity, n  4.40E−01  
Young’s Modulus, E (MPa)  307  
Poisson ratio, ν  0.4  
Bulk Modulus of bentonite, K_{s} (Pa)  2.0E + 08  
Compression Index, 1/K_{s}(Pa^{−1})  5.00E−09  
Volumesurface mean diameter, D_{vs} (m) calculated based on equation (23)  1.50E−09  
Initial tortuosity, τ  0.69  
Gas (Helium)  Density, ρ_{g} (kg m^{−3})  $pgMRT$ 
Dynamic Viscosity @ 293.15 K (Pa s)  2.0E−05  
Bulk Modulus of helium, K_{g} (Pa)  p_{g}  
Henry’s Coefficient, H  0.0091  
Diffusivity of helium in water (m^{2} s^{−1})  6.29E−09  
Diameter of helium gas particle, d_{he} (m)  1.0E−10  
Liquid (Water)  Density, ρ_{w} (kg m^{−3})  1000 
Dynamic Viscosity @ 293.15 K (Pa s)  0.001  
Bulk Modulus of water, K_{w} (Pa)  2.2E + 09  
Compression Index, 1/K_{w} (Pa^{−1})  4.55E−10 
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, $Cg,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
Rearranging (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}
where
$Cg,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 setup.
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 data provided by BGS (Daniels & Harrington 2017). 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 gasphase throughout the specimen and does not account for entrapment of gas within the pore space.
Parameter name  Parameter  Initial value condition 

Initial poregas pressure  p_{ginitial}  1.01E + 05 Pa 
Initial degree of saturation  S_{winitial}  0.96 
Initial suction (from SWCC)  s_{initial}  5.8E + 06 Pa 
Initial porewater pressure  p_{winitial}  −5.7E + 06 Pa 
Initial displacement field  u_{i}  0 m 
Initial stress  σ_{0xx} = σ_{0yy} = σ_{0zz}  5.0E + 05 Pa 
Initial gas concentration in porewater @STP  $Cg,H2O$  0.073 mol m^{−3} 
Parameter name  Parameter  Initial value condition 

Initial poregas pressure  p_{ginitial}  1.01E + 05 Pa 
Initial degree of saturation  S_{winitial}  0.96 
Initial suction (from SWCC)  s_{initial}  5.8E + 06 Pa 
Initial porewater pressure  p_{winitial}  −5.7E + 06 Pa 
Initial displacement field  u_{i}  0 m 
Initial stress  σ_{0xx} = σ_{0yy} = σ_{0zz}  5.0E + 05 Pa 
Initial gas concentration in porewater @STP  $Cg,H2O$  0.073 mol m^{−3} 
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 GenuchtenMualem 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 MX80 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.
Parameter  Value 

SWCC fitting parameter, a′(1m)  0.00065^{*} 
SWCC fitting parameter, n′  2.04 (Villar 2004) 
SWCC fitting parameter, m′  0.51 (Villar 2004) 
Relative permeability function fitting parameter, L′  0.51^{†} 
Parameter  Value 

SWCC fitting parameter, a′(1m)  0.00065^{*} 
SWCC fitting parameter, n′  2.04 (Villar 2004) 
SWCC fitting parameter, m′  0.51 (Villar 2004) 
Relative permeability function fitting parameter, L′  0.51^{†} 
calculated from the calibrated AEV using equations (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.
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 twophase flow through the MX80 bentonite specimen. Gas breakthrough (i.e. into the sample) was observed at 61 days when the gasinjection pressure reached 8 MPa. Over the 120day duration of the experiment gas migrated as a gas front through the first 5% of the sample.
The 2Daxisymmetric 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 120day 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 dilationcontrolled 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 dilationcontrolled 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 inline 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 porespace 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 buildup 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 2Daxisymetric 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.
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 reestablish 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 reestablishing 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 linearelastic 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 linearelastic 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.
S7–S8: Effect of AEV
The effect of increasing and decreasing the AEV, independent of porosity, on the porefluid pressure profiles are depicted in Figure 18(a) and on gas inflow in 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.
Material  Parameter  Value 

Solid soil matrix (Bentonite)  Compressive strength, f_{c} (MPa)  1 
Residual compressive strength, f_{cr} (MPa)  0.09  
Strain at compressive strength, ε_{c}_{0}  0.0003  
Tensile strength, f_{t} (MPa)  −1  
Residual tensile strength, f_{tr}(MPa)  −0.09  
Strain at tensile strength, ε_{t}_{0}  −0.0003  
Maximum permeability at maximum damage, k_{max} (m^{2})  1E−18  
c, smoothing coefficient  2 
Material  Parameter  Value 

Solid soil matrix (Bentonite)  Compressive strength, f_{c} (MPa)  1 
Residual compressive strength, f_{cr} (MPa)  0.09  
Strain at compressive strength, ε_{c}_{0}  0.0003  
Tensile strength, f_{t} (MPa)  −1  
Residual tensile strength, f_{tr}(MPa)  −0.09  
Strain at tensile strength, ε_{t}_{0}  −0.0003  
Maximum permeability at maximum damage, k_{max} (m^{2})  1E−18  
c, smoothing coefficient  2 
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 1E18 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 selfhealing were not reproduced. There is, therefore, still a need to investigate additional mechanisms to model dilatancycontrolled gas flow. Such mechanisms will be investigated in future studies. These will include the consideration of spatially weighted heterogeneity for porosity and a poroelastoplastic 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 selfhealing 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 2Daxisymetric 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.
Mesh qualifier  No. of elements  Degrees of freedom 

2D axisymmetric HM  
Extremely fine  6336  83 632 
Extra fine  1778  22 094 
Finer  602  7254 
Normal  214  2450 
3D HM  
Extra fine  81 067  1 278 368 
Finer  21 553  333 479 
Normal  4009  59 086 
Mesh qualifier  No. of elements  Degrees of freedom 

2D axisymmetric HM  
Extremely fine  6336  83 632 
Extra fine  1778  22 094 
Finer  602  7254 
Normal  214  2450 
3D HM  
Extra fine  81 067  1 278 368 
Finer  21 553  333 479 
Normal  4009  59 086 
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 2Daxisymmetric 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 2Daxisymmetric geometry, as the microfractures cannot intersect the axis of symmetry. Therefore, in order to model dilatancycontrolled 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 longterm safety assessment of a DGR is the longterm 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 hydromechanical linear poroelastic viscocapillary mathematical model for advectivediffusive controlled twophase 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 dilationcontrolled 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 dilatancycontrolled gas flow. The modelled results were compared to key features of the experimental data, including the gas inflow and outflow behaviour, porepressure fluid behaviour and evolution of radial stresses over the duration of the experiment. Although the model was not able to reproduce dilationcontrolled gas outflow from the sample, it was able to simulate twophase 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 linearelastic 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 lowpermeability swelling soils and was able to simulate twophase flow with dissolution of gas. Further development of the model to better mimic dilationcontrolled gasflow 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 dilatancycontrolled 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 poroelastoplastic 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 selfhealing behaviour of swelling soils will be addressed.
Acknowledgments
Decovalex (http://www.decovalex.org) is an international research project comprising participants from industry, government and academia focusing on development of understanding, models and codes in complex coupled problems in subsurface geological and engineering applications. DECOVALEX2019 is the current phase of the project. The authors appreciate and thank the DECOVALEX2019 funding organizations, Andra, BGR/UFZ, CNSC, US DOE, ENSI, JAEA, IRSN, KAERI, NWMO, RWM, SÚRAO, SSM and Taipower, for their financial and technical support of the work described in this paper. The statements made in the paper are, however, solely those of the authors and do not necessarily reflect those of the funding organizations.
Funding
This work was supported by the Canadian Nuclear Safety Commission.
Figures & Tables
Waste category  Waste inventory to the end of 2011  Waste inventory to the end of 2050 

Highlevel radioactive waste (nuclear fuel waste)  9400 m^{3}  20 000 m^{3} 
Intermediatelevel radioactive waste  33 400 m^{3}  67 000 m^{3} 
Lowlevel radioactive waste  2 343 000 m^{3}  2 594 000 m^{3} 
Waste category  Waste inventory to the end of 2011  Waste inventory to the end of 2050 

Highlevel radioactive waste (nuclear fuel waste)  9400 m^{3}  20 000 m^{3} 
Intermediatelevel radioactive waste  33 400 m^{3}  67 000 m^{3} 
Lowlevel radioactive waste  2 343 000 m^{3}  2 594 000 m^{3} 
Sensor name  Measurement type  Axial distance from injection face (mm) 

Injection load cell  Total stress  0 
Radial load cell 1  Total stress  15.2 
Radial load cell 2  Total stress  60 
Radial load cell 3  Total stress  104.8 
Backpressure load cell  Total stress  120 
Radial porewater array 1  Porefluid pressure  38.6 
Radial porewater array 2  Porefluid pressure  60 
Radial porewater array 3  Porefluid pressure  81.4 
Central filter (Middle)  Porefluid pressure  60 
Sensor name  Measurement type  Axial distance from injection face (mm) 

Injection load cell  Total stress  0 
Radial load cell 1  Total stress  15.2 
Radial load cell 2  Total stress  60 
Radial load cell 3  Total stress  104.8 
Backpressure load cell  Total stress  120 
Radial porewater array 1  Porefluid pressure  38.6 
Radial porewater array 2  Porefluid pressure  60 
Radial porewater array 3  Porefluid pressure  81.4 
Central filter (Middle)  Porefluid pressure  60 
Scenario number  Model  Description  Porosity  Intrinsic permeability (m^{2})  AEV (Pa)  SWCC fitting parameters  

a′ (1/m)  n′  m′  
S1Validation  HM  Advection, dissolution and diffusion linear poroelasticity  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Parametric studies  
S2  H  Advection, no dissolution, no diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S3  H  Advection, dissolution, no diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S4  H  Advection, dissolution, diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Sensitivity analysis  
S5  HM  ↑ por  0.48  5.1E−21  6.2E + 06  0.0008  2.04  0.51 
S6  HM  ↓ por  0.30  6.9E−22  1.6E + 07  0.0003  2.04  0.51 
S7  HM  ↑ AEV  0.44  3.4E−21  2.0E + 07  0.00025  2.04  0.51 
S8  HM  ↓ AEV  0.44  3.4E−21  5.0E + 06  0.0006  2.04  0.51 
S9  HM  ↑ k_{ij}  0.44  1.0E−19  8.0E + 06  0.0006  2.04  0.51 
S10  HM  ↓ k_{ij}  0.44  1.0E−22  8.0E + 06  0.0006  2.04  0.51 
S11  HM  Inclusion of a damage model  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S12  HM  Effect of mesh size  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Scenario number  Model  Description  Porosity  Intrinsic permeability (m^{2})  AEV (Pa)  SWCC fitting parameters  

a′ (1/m)  n′  m′  
S1Validation  HM  Advection, dissolution and diffusion linear poroelasticity  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Parametric studies  
S2  H  Advection, no dissolution, no diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S3  H  Advection, dissolution, no diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S4  H  Advection, dissolution, diffusion  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Sensitivity analysis  
S5  HM  ↑ por  0.48  5.1E−21  6.2E + 06  0.0008  2.04  0.51 
S6  HM  ↓ por  0.30  6.9E−22  1.6E + 07  0.0003  2.04  0.51 
S7  HM  ↑ AEV  0.44  3.4E−21  2.0E + 07  0.00025  2.04  0.51 
S8  HM  ↓ AEV  0.44  3.4E−21  5.0E + 06  0.0006  2.04  0.51 
S9  HM  ↑ k_{ij}  0.44  1.0E−19  8.0E + 06  0.0006  2.04  0.51 
S10  HM  ↓ k_{ij}  0.44  1.0E−22  8.0E + 06  0.0006  2.04  0.51 
S11  HM  Inclusion of a damage model  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
S12  HM  Effect of mesh size  0.44  3.4E−21  8.0E + 06  0.0006  2.04  0.51 
Model geometry  Mesh qualifier  No. of elements  Degrees of freedom 

2D axisymmetric  Extremely fine  8886  117 132 
3D  Normal  10 698  161 741 
Model geometry  Mesh qualifier  No. of elements  Degrees of freedom 

2D axisymmetric  Extremely fine  8886  117 132 
3D  Normal  10 698  161 741 
Material  Parameter  Value 

Solid soil matrix (Bentonite)  Saturated/intrinsic permeability, k_{ij} (m^{2})  3.40E−21 
Dry density, ρ_{d} (kg m^{−3})  1560  
Maximum degree of saturation, S_{ws}  1  
Residual degree of Saturation, S_{wr}  0.05  
Porosity, n  4.40E−01  
Young’s Modulus, E (MPa)  307  
Poisson ratio, ν  0.4  
Bulk Modulus of bentonite, K_{s} (Pa)  2.0E + 08  
Compression Index, 1/K_{s}(Pa^{−1})  5.00E−09  
Volumesurface mean diameter, D_{vs} (m) calculated based on equation (23)  1.50E−09  
Initial tortuosity, τ  0.69  
Gas (Helium)  Density, ρ_{g} (kg m^{−3})  $pgMRT$ 
Dynamic Viscosity @ 293.15 K (Pa s)  2.0E−05  
Bulk Modulus of helium, K_{g} (Pa)  p_{g}  
Henry’s Coefficient, H  0.0091  
Diffusivity of helium in water (m^{2} s^{−1})  6.29E−09  
Diameter of helium gas particle, d_{he} (m)  1.0E−10  
Liquid (Water)  Density, ρ_{w} (kg m^{−3})  1000 
Dynamic Viscosity @ 293.15 K (Pa s)  0.001  
Bulk Modulus of water, K_{w} (Pa)  2.2E + 09  
Compression Index, 1/K_{w} (Pa^{−1})  4.55E−10 
Material  Parameter  Value 

Solid soil matrix (Bentonite)  Saturated/intrinsic permeability, k_{ij} (m^{2})  3.40E−21 
Dry density, ρ_{d} (kg m^{−3})  1560  
Maximum degree of saturation, S_{ws}  1  
Residual degree of Saturation, S_{wr}  0.05  
Porosity, n  4.40E−01  
Young’s Modulus, E (MPa)  307  
Poisson ratio, ν  0.4  
Bulk Modulus of bentonite, K_{s} (Pa)  2.0E + 08  
Compression Index, 1/K_{s}(Pa^{−1})  5.00E−09  
Volumesurface mean diameter, D_{vs} (m) calculated based on equation (23)  1.50E−09  
Initial tortuosity, τ  0.69  
Gas (Helium)  Density, ρ_{g} (kg m^{−3})  $pgMRT$ 
Dynamic Viscosity @ 293.15 K (Pa s)  2.0E−05  
Bulk Modulus of helium, K_{g} (Pa)  p_{g}  
Henry’s Coefficient, H  0.0091  
Diffusivity of helium in water (m^{2} s^{−1})  6.29E−09  
Diameter of helium gas particle, d_{he} (m)  1.0E−10  
Liquid (Water)  Density, ρ_{w} (kg m^{−3})  1000 
Dynamic Viscosity @ 293.15 K (Pa s)  0.001  
Bulk Modulus of water, K_{w} (Pa)  2.2E + 09  
Compression Index, 1/K_{w} (Pa^{−1})  4.55E−10 
Parameter name  Parameter  Initial value condition 

Initial poregas pressure  p_{ginitial}  1.01E + 05 Pa 
Initial degree of saturation  S_{winitial}  0.96 
Initial suction (from SWCC)  s_{initial}  5.8E + 06 Pa 
Initial porewater pressure  p_{winitial}  −5.7E + 06 Pa 
Initial displacement field  u_{i}  0 m 
Initial stress  σ_{0xx} = σ_{0yy} = σ_{0zz}  5.0E + 05 Pa 
Initial gas concentration in porewater @STP  $Cg,H2O$  0.073 mol m^{−3} 
Parameter name  Parameter  Initial value condition 

Initial poregas pressure  p_{ginitial}  1.01E + 05 Pa 
Initial degree of saturation  S_{winitial}  0.96 
Initial suction (from SWCC)  s_{initial}  5.8E + 06 Pa 
Initial porewater pressure  p_{winitial}  −5.7E + 06 Pa 
Initial displacement field  u_{i}  0 m 
Initial stress  σ_{0xx} = σ_{0yy} = σ_{0zz}  5.0E + 05 Pa 
Initial gas concentration in porewater @STP  $Cg,H2O$  0.073 mol m^{−3} 
Parameter  Value 

SWCC fitting parameter, a′(1m)  0.00065^{*} 
SWCC fitting parameter, n′  2.04 (Villar 2004) 
SWCC fitting parameter, m′  0.51 (Villar 2004) 
Relative permeability function fitting parameter, L′  0.51^{†} 
Parameter  Value 

SWCC fitting parameter, a′(1m)  0.00065^{*} 
SWCC fitting parameter, n′  2.04 (Villar 2004) 
SWCC fitting parameter, m′  0.51 (Villar 2004) 
Relative permeability function fitting parameter, L′  0.51^{†} 
calculated from the calibrated AEV using equations (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.
Material  Parameter  Value 

Solid soil matrix (Bentonite)  Compressive strength, f_{c} (MPa)  1 
Residual compressive strength, f_{cr} (MPa)  0.09  
Strain at compressive strength, ε_{c}_{0}  0.0003  
Tensile strength, f_{t} (MPa)  −1  
Residual tensile strength, f_{tr}(MPa)  −0.09  
Strain at tensile strength, ε_{t}_{0}  −0.0003  
Maximum permeability at maximum damage, k_{max} (m^{2})  1E−18  
c, smoothing coefficient  2 
Material  Parameter  Value 

Solid soil matrix (Bentonite)  Compressive strength, f_{c} (MPa)  1 
Residual compressive strength, f_{cr} (MPa)  0.09  
Strain at compressive strength, ε_{c}_{0}  0.0003  
Tensile strength, f_{t} (MPa)  −1  
Residual tensile strength, f_{tr}(MPa)  −0.09  
Strain at tensile strength, ε_{t}_{0}  −0.0003  
Maximum permeability at maximum damage, k_{max} (m^{2})  1E−18  
c, smoothing coefficient  2 
Mesh qualifier  No. of elements  Degrees of freedom 

2D axisymmetric HM  
Extremely fine  6336  83 632 
Extra fine  1778  22 094 
Finer  602  7254 
Normal  214  2450 
3D HM  
Extra fine  81 067  1 278 368 
Finer  21 553  333 479 
Normal  4009  59 086 
Mesh qualifier  No. of elements  Degrees of freedom 

2D axisymmetric HM  
Extremely fine  6336  83 632 
Extra fine  1778  22 094 
Finer  602  7254 
Normal  214  2450 
3D HM  
Extra fine  81 067  1 278 368 
Finer  21 553  333 479 
Normal  4009  59 086 