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

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 high-level 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³ 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₂) and the radiolysis of water producing H₂ (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 build-up and migration through an engineered barrier system (EBS) (Sellin 2014; Shaw 2015b).

Expansive or swelling soils, such as bentonite-based 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 (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 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 dilatency-controlled 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 inter- and 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 hydro-mechanical 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 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 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 dilatancy-controlled gas migration. This research is in part a contribution to Task A of DECOVALEX-2019.

This paper presents a fully-coupled, hydro-mechanical linear-elastic mathematical model for advective-diffusive visco-capillary controlled two-phase 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.

An HM model to describe the migration of gas (two-phase 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.

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).

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, ν

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⁻³)

• g is the acceleration due to gravity (m s⁻²)

• 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 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.

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 re-arranging equation (16) we can express this relationship in terms of porosity,

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⁻¹)

• (v_ig−v_is) is the velocity of gas in the pores relative to the velocity of the soil matrix (m s⁻¹)

• n porosity (m³ voids · m⁻³ total)

• μ_w dynamic viscosity of the water phase (Pa s or kg m⁻¹ s⁻¹)

• μ_g dynamic viscosity of the gas phase (Pa s or kg m⁻¹ s⁻¹)

• k_ij intrinsic permeability tensor of the porous medium (m²)

• 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⁻³)

• ρ_g density of the gas phase (kg m⁻³).

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₀ is the specific surface (surface exposed to fluid per unit solid volume) (cm² cm⁻³).

For fine-grained expansive clays, which have very small grain-sizes less than 2 µm, and are non-uniform, 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 volume-surface mean diameter to better account for the non-uniformity 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.

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).

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² s⁻¹)

• D is the diffusivity of gas in water (m² s⁻¹)

• 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 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,

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)

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₀ 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²)

• k_UD is the permeability before damage (m²) and is expressed by equation (21)

• k_D is the increase in permeability as a result of damage (m²) and can be expressed by equation (36)

where

• k_max(m²) 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.

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⁻² s⁻²)

• ∂σ_ij/∂x_j represents the change in normal and shear stresses across the soil element (kg m⁻² s⁻²).

Substituting equation (8) and differentiating, the governing equation for the conservation of momentum using a linear poro-elasticity model is expressed by equation (40).

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,

• (∂(ρ_wnS))/∂t is the rate of change of mass of water in the soil element over time (kg m⁻³ s⁻¹)

• ∂/∂x_i(ρ_wn(v_iw−v_is)) is the net advective flux of water for the soil element (kg m⁻³ s⁻¹).

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 ∂(ρ_wnS_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⁻³).

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),

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⁻³)

• ∂/∂x_i(ρ_gn((v_ig−v_is) + H(v_iw − v_is))) is the net advective flux of gas for the soil element (kg m⁻³ s⁻¹)

• ∂/∂x_i(−D_e(∂/∂x_k)(ρ_gn(HS_w))) is the net diffusive flux of gas for the soil element (kg m⁻³ s⁻¹)

• ∂(ρ_gn(1−S_w + HS))/∂t is the rate of change of mass of gas in the soil element over time (kg m⁻³ s⁻¹)

• H is Henry’s coefficient (kg species A m⁻³ in aqueous phase kg⁻¹ species A m³ 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).

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).

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.

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 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).

The hydraulic and mechanical BCs for gas transport, water transport and momentum transport are provided below.

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

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

where

• $Cg,H2O$ is the concentration of gas in water (kg gas m⁻³ water)

• H is Henry’s coefficient (kg species A m⁻³ in aqueous phase kg⁻¹ species A m³ in gas phase).

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.

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).

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 gas-phase throughout the specimen and does not account for entrapment of gas within the pore space.

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.

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.

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 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.

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.

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 self-healing 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.

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.

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.

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⁻¹⁹ m² 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⁻²² m² 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² 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.

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² 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 dilation-controlled 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.

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.

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 advective-diffusive controlled two-phase flow using a Bishop’s effective stress. The objective of this study was to gain insight into the key parameters influencing gas-flow 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 dilatancy-controlled 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.

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 sub-surface geological and engineering applications. DECOVALEX-2019 is the current phase of the project. The authors appreciate and thank the DECOVALEX-2019 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.

This work was supported by the Canadian Nuclear Safety Commission.