Skip to Main Content
Skip Nav Destination

Correspondence: jinfante@uottawa.ca

This is an Open Access article distributed under the terms of the Creative Commons Attribution 3.0 License (http://creativecommons.org/licenses/by/3.0/).

Abstract

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

In Canada, nuclear waste has been generating and accumulating since the 1930s when the Port Radium radium mine began operating in the Northwest Territories (Low-Level Radioactive Waste Management Office 2012). Since then Canada has become a sustainable nuclear country with operating nuclear facilities across the nuclear fuel cycle, all producing various forms of radioactive waste. To-date, this waste, in the form of low-, intermediate- and 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 m3 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.

Table 1.

Canadian Waste Inventory Projections to 2011 and 2050 (adapted fromLLRWMO 2012)

Waste categoryWaste inventory to the end of 2011Waste inventory to the end of 2050
High-level radioactive waste (nuclear fuel waste)9400 m320 000 m3
Intermediate-level radioactive waste33 400 m367 000 m3
Low-level radioactive waste2 343 000 m32 594 000 m3
Waste categoryWaste inventory to the end of 2011Waste inventory to the end of 2050
High-level radioactive waste (nuclear fuel waste)9400 m320 000 m3
Intermediate-level radioactive waste33 400 m367 000 m3
Low-level radioactive waste2 343 000 m32 594 000 m3

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 (H2) and the radiolysis of water producing H2 (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.

Fig. 1.

Four main processes of gas migration in clays (Cuss et al. 2014).

Fig. 1.

Four main processes of gas migration in clays (Cuss et al. 2014).

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.

Mathematical model

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, kij, as a function of changing porosity; and

  • (5) consideration of a damage model.

The applicable constitutive relations and governing equations for conservation of momentum, water mass and gas mass are presented below.

Constitutive relations for the mechanical behaviour

Bishop’s modified effective stress principle

Many effective stress equations have been proposed to characterize the stress-state of an unsaturated soil or porous media (Nguyen & Le 2015). This paper proposes the use of Bishop’s effective stress principle, which is dependent on both net normal stress and matric suction, and may be more suitable for expansive clay’s, as described by equation (1).

 

σ=(σpg)+χ(pgpw)
(1)

where

  • σ′ is the effective stress (Pa)

  • σ is the total normal stress (Pa)

  • pg is the poregas pressure (Pa)

  • pw is the porewater pressure (Pa)

  • (σpg) is the net normal stress (Pa)

  • (pgpw) is the matric suction (Pa)

  • χ is a parameter related to the degree of saturation of the soil (unitless).

Expanding and rearranging for σ

 

σ=σ+(1χ)pg+χpw
(2)

The porefluid pressure can be defined as,

 

p¯=(1χ)pg+χpw
(3)

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,

 

χ={[(pgpw)(pgpw)b]0.55,if(pgpw)>(pgpw)b1,otherwise
(4)

where (pgpw)b is the AEV of the soil and only applies when the matric suction >AEV (Khalili & Khabbaz 1998).

Equation of poro-elasticity

The general form for the equation of poromechanics can be expressed by,

 

σij=Cijklεkl+αBp¯δij
(5)

where

  • σij total normal stress tensor acting on the soil element (Pa)

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

 

σij=2Gεij+λεkkδij+αBp¯δij
(6)

where

 

εij=12(ui,j+uj,i)=12(uixj+ujxi)
(7)

 

εkk=uk,k=ukxk
(8)

where

  • G is the shear modulus (Pa)

  • λ is the Lamé constant (Pa)

  • εkk (=uk,k = trεij) is the volumetric strain (unitless)

  • uij is the displacement tensor (m).

substituting into equation (6)

 

σij=G(uixj+ujxi)+(λukxk+αBp¯)δij
(9)

or in terms of Young’s modulus, E, and Poisson’s ratio, ν

 

σij=E2(1+ν)(uixj+ujxi)+(νE(1+ν)(12ν)ukxk+αBp¯)δij
(10)

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

 

Se,w={1[1+|a(sρwg)|n]m,s>01,s0
(11)

where

  • Se,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 (pgpw) (Pa).

The fitting parameter, a′, is related to the AEV and can be expressed by equation (12)

 

a=0.5ρwgAEV
(12)

The fitting parameter m′ is related to fitting parameter n′ by equation (13),

 

m=11n
(13)

Se,w can be calculated from scaling the water saturation, Sw, with respect to the maximum degree of saturation, Sw,s, and the residual degree of saturation, Sw,r

 

Se,w=SwSw,rSw,sSw,r
(14)

From the fitted SWCC, the model can be used to predict the degree of saturation of water, Sw, from measuring the stress-state variable matric suction (or capillary pressure), s.

The degree of saturation of the gas phase, Sg, can be related to Sw via the following relationship

 

Sg=1Sw
(15)

It should be noted that the model does not consider entrapped gases which would be present if Sw,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.

 

log(saev)=log(saev)+m(eaeveaev)
(16)

where

  • saev is the AEV (Pa) for a void ratio, eaev

  • s′aev is the arbitrary reference AEV (Pa) for a void ratio e′aev

  • m is the slope of the saeveaev 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,

 

log(saev)=log(saev)+m(naev1naevnaev1naev)
(17)

Darcy’s Law for two-phase flow

Darcy’s Law for the water and gas phases are expressed in equations (18) and (19).

 

viwvis=kijkrwnμw(pwxj+ρwg)
(18)

 

vigvis=kijkrgnμg(pgxj+ρgg)
(19)

where

  • pw is the porewater pressure (Pa)

  • pg is the poregas pressure (Pa)

  • (viw−vis) is the velocity of water in the pores relative to the velocity of the soil matrix (m s−1)

  • (vig−vis) is the velocity of gas in the pores relative to the velocity of the soil matrix (m s−1)

  • n porosity (m3 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)

  • kij intrinsic permeability tensor of the porous medium (m2)

  • krw relative permeability of the water phase (unitless)

  • krg relative permeability of the gas phase (unitless)

  • ρw density of water phase (kg m−3)

  • ρg density of the gas phase (kg m−3).

Intrinsic permeability as a function of porosity

Numerous equations have been proposed to describe the relationship between intrinsic permeability and measurable properties of soils such as the porosity and particle-size. The Kozeny-Carman equation can be used to express the intrinsic permeability as a function of porosity using the capillary tube (hydraulic radius) model and is sufficient for laminar flow and low porefluid velocity (Carman 1956). The Kozeny-Carmen Equation is given by,

 

kij=15S02n3(1n)2
(20)

where S0 is the specific surface (surface exposed to fluid per unit solid volume) (cm2 cm−3).

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.

 

kij=Dvs2180n3(1n)2
(21)

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

 

Dvs=180kij(1n)2n3
(22)

In this study, the experimentally determined initial porosity and intrinsic permeability were used to calculate the volume-surface mean diameter of the soil sample. Any elastic deformations to the soil structure will influence the porosity, subsequently influencing the intrinsic permeability across the soil specimen.

Relative permeability of water and gas phases

The coefficient of permeability is strongly dependent on the degree of saturation, Sw. Using Mualem’s Model (Nguyen & Le 2015), krw and krg, can be calculated using the following expressions:

 

krw={Se,wL[1(1Se,w1/m)m]2,s>01,s0
(23)

 

krg=(1Se,w)L[1Se,w1/m]2m
(24)

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

 

De=DnSwτ
(25)

where

  • De is the effective diffusivity of gas dissolved in water through porous media (m2 s−1)

  • D is the diffusivity of gas in water (m2 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,

 

τ=n13Sw73
(26)

Storage due to suction of water and gas (storage coefficient)

Storage terms are applied in the governing equations when differentiating dSw/ds. The storage coefficient of water, Cw, can be calculated by differentiating the van Genuchten SWCC,

 

Cw=dSwds={1ρwgam1m(Sw,sSw,r)Se,w1/m(1Se,w1/m)m,s>00,s0
(27)

The storage coefficient of gas, Cg, can be related to the storage coefficient of water by,

 

Cg=dSgds=d(1Sw)ds=d(Sw)ds=Cw
(28)

Bulk modulus of water and gas

The bulk modulus of water, Kw, can be calculated as follows

 

Kw=VwdpwdVw=ρwdpwdρw
(29)

The bulk modulus of gas, Kg, can be calculated as follows

 

Kg=VgdpgdVg=ρgdpgdρg
(30)

From the ideal gas law,

 

pgρg=RTM
(31)

From Boyle’s Law for isothermal conditions,

 

dpgdρg=RTM
(32)

substituting equation (32) back into equation (30)

 

Kg=ρgRTM=pg
(33)

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,

 

E=E0(1Dc)
(34)

where

  • D is the damage (unitless)

  • E0 is the Young’s modulus before damage (Pa)

  • E is the Young’s modulus after damage (Pa)

 

kij=kUD+kD
(35)

where

  • kij is the permeability at any given time (m2)

  • kUD is the permeability before damage (m2) and is expressed by equation (21)

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

 

kD=Dc(kmaxkUD)
(36)

where

  • kmax(m2) 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

  • kmax and c, are determined through model calibration.

The criterion for Dt in multi-axial tension is expressed by equation (37), while the criterion for Dc in multi-axial compression is expressed by equation (38)

 

Dt={0εεt01ftrE0εεt0εεtu1εεtu
(37)

 

Dc={0εc0ε1fcrE0εεεc0
(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

  • ftr is the residual tensile strength (Pa)

  • fcr 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)

 

σijxj+Fv,i=0
(39)

where,

  • σij is the stress tensor (Pa)

  • Fv,i is the volumetric body force tensor (kg m−2 s−2)

  • σij/∂xj represents the change in normal and shear stresses across the soil element (kg m−2 s−2).

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

 

G2uixjxj+(G+λ)2ujxjxi+αB(1χ)pgxi+χpwxi+Fv,i=0
(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).

Fig. 2.

Control volume and phase relations.

Fig. 2.

Control volume and phase relations.

 

xi(ρwn(viwvis))=(ρwnSw)t
(41)

where,

  • (∂(ρwnS))/∂t is the rate of change of mass of water in the soil element over time (kg m−3 s−1)

  • ∂/∂xi(ρwn(viw−vis)) 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 right-hand side of equation (41) results in the following governing equations for the conservation of water mass,

 

xi(ρwkijkrwμw(pwxj+ρwgj))=ρw[n(dSwds)st+nSwKwpwt+Sw(1n)t(ukxk)]
(42)

where Kw is the bulk modulus of water (Pa).

It should be noted that the RHS of equation (42) is derived from the differentiation of ∂(ρwnSw)∂t. The final term in equation (42) is multiplied by (1−n) which is derived from differentiating the equation for porosity,

 

n=VvV
(43)

where Vv is the volume of voids (m−3).

Using the Quotient Rule to differentiate,

 

dn=d(VvV)=dVvVVvdVV2=dVvVVvVdVV=d(VVs)VndVV=dVV(1n)dVsV
(44)

where, dV/V = kk, the volumetric strain, and where it is assumed that dVs/V ≪ kk(1−n),

 

dn=d(VvV)=dεkk(1n)=(1n)(ukxk)
(45)

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

 

(ρgn((vigvis)+H(viwvis)))xi+xi(Dexk(ρgn(HSw)))=(ρgn(1Sw+HSw))t
(46)

where,

  • ρg, is the gas density (kg gas m−3)

  • ∂/∂xi(ρgn((vig−vis) + H(viw − vis))) is the net advective flux of gas for the soil element (kg m−3 s−1)

  • ∂/∂xi(−De(∂/∂xk)(ρgn(HSw))) is the net diffusive flux of gas for the soil element (kg m−3 s−1)

  • ∂(ρgn(1−Sw + 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 m3 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

 

(ρg(kijkrgμg(pgxj+ρgg)+Hkijkrwμw(pwxj+ρwg)))xi+xi(Dexk(ρgn(HSw)))=ρg[n(H1)(dSwds)st+n(1Sw+HSw)Kgpgt+(1Sw+HSw)(1n)t(ukxk)]
(47)

where Kg is the bulk modulus of the gas phase (Pa).

Numerical model description and modelling approach

Overview of the numerical model

A 2D-axisymmetric and 3D time-dependent (i.e. transient) HM-coupled multiphysics models were developed. The models simulate the simultaneous migration of gas and liquid (two-phase flow) in porous media, which are coupled to the linear poro-elastic behaviour of the solid matrix using the FEM. The commercially available code COMSOL Multiphysics® was used to numerically solve the governing equations of the model.

The numerical model was based on an experimental setup conducted by the BGS to assess the 1D flow through a saturated bentonite sample under a constant volume boundary stress condition (Daniels & Harrington 2017). In this experiment, a confined cylindrical sample of near-saturated bentonite was injected on one side with helium gas with increasing gas pressures over a period of 120 days. The other side was left at a constant water backpressure throughout the duration of the experiment. The experiment was conducted under isothermal conditions at a temperature of 293.15 K. During the experiment, a number of parameters where measured including the gas inflow and outflow, the porefluid pressure at defined porefluid arrays, and the total radial stresses at radial load cell arrays. The cylindrical specimen of MX-80 bentonite had a diameter of 60 mm and a length of 120 mm. Table 2 provides the location of the monitoring sensors within the specimen. The BGS provided the experimental data that was used in the model validation (Daniels & Harrington 2017).

Table 2.

Model dimensions and location of the porefluid arrays

Sensor nameMeasurement typeAxial distance from injection face (mm)
Injection load cellTotal stress0
Radial load cell 1Total stress15.2
Radial load cell 2Total stress60
Radial load cell 3Total stress104.8
Backpressure load cellTotal stress120
Radial porewater array 1Porefluid pressure38.6
Radial porewater array 2Porefluid pressure60
Radial porewater array 3Porefluid pressure81.4
Central filter (Middle)Porefluid pressure60
Sensor nameMeasurement typeAxial distance from injection face (mm)
Injection load cellTotal stress0
Radial load cell 1Total stress15.2
Radial load cell 2Total stress60
Radial load cell 3Total stress104.8
Backpressure load cellTotal stress120
Radial porewater array 1Porefluid pressure38.6
Radial porewater array 2Porefluid pressure60
Radial porewater array 3Porefluid pressure81.4
Central filter (Middle)Porefluid pressure60

Modelling approach

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

Table 3.

Base study and parametric study model scenarios

Scenario numberModelDescriptionPorosityIntrinsic permeability (m2)AEV (Pa)SWCC fitting parameters
a′ (1/m)nm
S1-ValidationHMAdvection, dissolution and diffusion linear poro-elasticity0.443.4E−218.0E + 060.00062.040.51
Parametric studies
S2HAdvection, no dissolution, no diffusion0.443.4E−218.0E + 060.00062.040.51
S3HAdvection, dissolution, no diffusion0.443.4E−218.0E + 060.00062.040.51
S4HAdvection, dissolution, diffusion0.443.4E−218.0E + 060.00062.040.51
Sensitivity analysis
S5HM↑ por0.485.1E−216.2E + 060.00082.040.51
S6HM↓ por0.306.9E−221.6E + 070.00032.040.51
S7HM↑ AEV0.443.4E−212.0E + 070.000252.040.51
S8HM↓ AEV0.443.4E−215.0E + 060.00062.040.51
S9HMkij0.441.0E−198.0E + 060.00062.040.51
S10HMkij0.441.0E−228.0E + 060.00062.040.51
S11HMInclusion of a damage model0.443.4E−218.0E + 060.00062.040.51
S12HMEffect of mesh size0.443.4E−218.0E + 060.00062.040.51
Scenario numberModelDescriptionPorosityIntrinsic permeability (m2)AEV (Pa)SWCC fitting parameters
a′ (1/m)nm
S1-ValidationHMAdvection, dissolution and diffusion linear poro-elasticity0.443.4E−218.0E + 060.00062.040.51
Parametric studies
S2HAdvection, no dissolution, no diffusion0.443.4E−218.0E + 060.00062.040.51
S3HAdvection, dissolution, no diffusion0.443.4E−218.0E + 060.00062.040.51
S4HAdvection, dissolution, diffusion0.443.4E−218.0E + 060.00062.040.51
Sensitivity analysis
S5HM↑ por0.485.1E−216.2E + 060.00082.040.51
S6HM↓ por0.306.9E−221.6E + 070.00032.040.51
S7HM↑ AEV0.443.4E−212.0E + 070.000252.040.51
S8HM↓ AEV0.443.4E−215.0E + 060.00062.040.51
S9HMkij0.441.0E−198.0E + 060.00062.040.51
S10HMkij0.441.0E−228.0E + 060.00062.040.51
S11HMInclusion of a damage model0.443.4E−218.0E + 060.00062.040.51
S12HMEffect of mesh size0.443.4E−218.0E + 060.00062.040.51

Model geometry and mesh

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

Fig. 3.

2D-axisymmetric geometry and meshing.

Fig. 3.

2D-axisymmetric geometry and meshing.

Fig. 4.

3D geometry and meshing.

Fig. 4.

3D geometry and meshing.

Table 4.

Number of elements and degrees of freedom

Model geometryMesh qualifierNo. of elementsDegrees of freedom
2D axisymmetricExtremely fine8886117 132
3DNormal10 698161 741
Model geometryMesh qualifierNo. of elementsDegrees of freedom
2D axisymmetricExtremely fine8886117 132
3DNormal10 698161 741

Material properties

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

Table 5.

Material properties of the solid soil matrix and gas and liquid phases

MaterialParameterValue
Solid soil matrix (Bentonite)Saturated/intrinsic permeability, kij (m2)3.40E−21
Dry density, ρd (kg m−3)1560
Maximum degree of saturation, Sws1
Residual degree of Saturation, Swr0.05
Porosity, n4.40E−01
Young’s Modulus, E (MPa)307
Poisson ratio, ν0.4
Bulk Modulus of bentonite, Ks (Pa)2.0E + 08
Compression Index, 1/Ks(Pa−1)5.00E−09
Volume-surface mean diameter, Dvs (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, Kg (Pa)pg
Henry’s Coefficient, H0.0091
Diffusivity of helium in water (m2 s−1)6.29E−09
Diameter of helium gas particle, dhe (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, Kw (Pa)2.2E + 09
Compression Index, 1/Kw (Pa−1)4.55E−10
MaterialParameterValue
Solid soil matrix (Bentonite)Saturated/intrinsic permeability, kij (m2)3.40E−21
Dry density, ρd (kg m−3)1560
Maximum degree of saturation, Sws1
Residual degree of Saturation, Swr0.05
Porosity, n4.40E−01
Young’s Modulus, E (MPa)307
Poisson ratio, ν0.4
Bulk Modulus of bentonite, Ks (Pa)2.0E + 08
Compression Index, 1/Ks(Pa−1)5.00E−09
Volume-surface mean diameter, Dvs (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, Kg (Pa)pg
Henry’s Coefficient, H0.0091
Diffusivity of helium in water (m2 s−1)6.29E−09
Diameter of helium gas particle, dhe (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, Kw (Pa)2.2E + 09
Compression Index, 1/Kw (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, ∂pg/∂xi = 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

 

pgV=nRT
(48)

 

pgV=mgMRT
(49)

Re-arranging (50)

 

mgV=ρg=pgMRT
(50)

where pg 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, nSw

 

Cg,H2O=ρgH(nSw)
(51)

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 m3 in gas phase).

Water conservation BCs

For water transport, a no flow Neumann BC, ∂pw/∂xi = 0, was set for the lower and radial boundaries. A Dirichlet BC was set as the water backpressure for the upper boundary.

Momentum conservation BCs

For the moment conservation equation, a roller constraint was applied along the upper, lower and radial boundaries to represent a condition whereby the boundary is free to move in the tangential direction, but is fixed in the normal direction, simulating a constant volume condition, as per the experimental set-up.

The gas injection pressure and water backpressure Dirichlet BCs were imported from the experimental data provided by BGS and have been plotted in Figure 5(a), while the concentration of dissolved gas BCs have been plotted in Figure 5(b).

Fig. 5.

Dirichlet BCs for the (a) lower gas injection pressure and upper water backpressure (b) lower and upper gas concentrations in porewater.

Fig. 5.

Dirichlet BCs for the (a) lower gas injection pressure and upper water backpressure (b) lower and upper gas concentrations in porewater.

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

Table 6.

Initial value conditions

Parameter nameParameterInitial value condition
Initial poregas pressurepginitial1.01E + 05 Pa
Initial degree of saturationSwinitial0.96
Initial suction (from SWCC)sinitial5.8E + 06 Pa
Initial porewater pressurepwinitial−5.7E + 06 Pa
Initial displacement fieldui0 m
Initial stressσ0xx = σ0yy = σ0zz5.0E + 05 Pa
Initial gas concentration in porewater @STPCg,H2O0.073 mol m−3
Parameter nameParameterInitial value condition
Initial poregas pressurepginitial1.01E + 05 Pa
Initial degree of saturationSwinitial0.96
Initial suction (from SWCC)sinitial5.8E + 06 Pa
Initial porewater pressurepwinitial−5.7E + 06 Pa
Initial displacement fieldui0 m
Initial stressσ0xx = σ0yy = σ0zz5.0E + 05 Pa
Initial gas concentration in porewater @STPCg,H2O0.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.

Fig. 6.

Intrinsic permeability (a), effectivity diffusivity (b) and AEV (c) curves as a function of porosity.

Fig. 6.

Intrinsic permeability (a), effectivity diffusivity (b) and AEV (c) curves as a function of porosity.

Fig. 7.

χ as a function of both AEV and suction.

Fig. 7.

χ as a function of both AEV and suction.

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.

Fig. 8.

Soil water and soil gas characteristic curves for bentonite at n = 0.44.

Fig. 8.

Soil water and soil gas characteristic curves for bentonite at n = 0.44.

Fig. 9.

Relative permeability curves (water and gas) in bentonite v. suction at n = 0.44.

Fig. 9.

Relative permeability curves (water and gas) in bentonite v. suction at n = 0.44.

Table 7.

van Genuchten-Mualem SWCC and relative permeability function fitting parameters

ParameterValue
SWCC fitting parameter, a′(1m)0.00065*
SWCC fitting parameter, n2.04 (Villar 2004)
SWCC fitting parameter, m0.51 (Villar 2004)
Relative permeability function fitting parameter, L0.51
ParameterValue
SWCC fitting parameter, a′(1m)0.00065*
SWCC fitting parameter, n2.04 (Villar 2004)
SWCC fitting parameter, m0.51 (Villar 2004)
Relative permeability function fitting parameter, L0.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 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.

Fig. 10.

S1 results: gas migration through bentonite MX-80 sample. For the 3D model, at times (a) 0 days, (b) 65 days and (c) 120 days; and for the 2D-axisymmetric model, at times (d) 0 days, (e) 65 days and (f) 120 days.

Fig. 10.

S1 results: gas migration through bentonite MX-80 sample. For the 3D model, at times (a) 0 days, (b) 65 days and (c) 120 days; and for the 2D-axisymmetric model, at times (d) 0 days, (e) 65 days and (f) 120 days.

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.

Fig. 11.

S1 results: gas (a) inflow and (b) outflow profiles over time.

Fig. 11.

S1 results: gas (a) inflow and (b) outflow profiles over time.

Fig. 12.

S1 results: porefluid pressure profiles over (a) time and (b) axial length of the specimen.

Fig. 12.

S1 results: porefluid pressure profiles over (a) time and (b) axial length of the specimen.

Fig. 13.

S1 results – Total radial stresses over time.

Fig. 13.

S1 results – Total radial stresses over time.

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.

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.

Fig. 14.

S2–S4 results: (a) contribution of advection, dissolution and diffusion to flow behaviour and (b) diffusion gradient over axial length of the soil specimen.

Fig. 14.

S2–S4 results: (a) contribution of advection, dissolution and diffusion to flow behaviour and (b) diffusion gradient over axial length of the soil specimen.

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.

Fig. 15.

S2–S4 results: contribution of (a) advection, dissolution and diffusion and (b) elastic deformation to porefluid pressure profiles.

Fig. 15.

S2–S4 results: contribution of (a) advection, dissolution and diffusion and (b) elastic deformation to porefluid pressure profiles.

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.

Fig. 16.

Change in porosity over time as a result of elastic deformation.

Fig. 16.

Change in porosity over time as a result of elastic deformation.

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.

Fig. 17.

S5–S6 results: effect of changes in initial porosity on (a) porefluid pressure and (b) gas inflow.

Fig. 17.

S5–S6 results: effect of changes in initial porosity on (a) porefluid pressure and (b) gas inflow.

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.

Fig. 18.

S7–S8 results: effect of changes in AEV on (a) porefluid pressure and (b) gas inflow.

Fig. 18.

S7–S8 results: effect of changes in AEV on (a) porefluid pressure and (b) gas inflow.

S9–S10: Effect of intrinsic permeability

The effect of increasing and decreasing the intrinsic permeability, kij, 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 kij by two orders of magnitude to 10−19 m2 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.

Fig. 19.

S9–S10 results: effect of changes in kij on (a) porefluid pressure and (b) gas inflow.

Fig. 19.

S9–S10 results: effect of changes in kij on (a) porefluid pressure and (b) gas inflow.

For the S10 scenario, a decrease of kij by over an order of magnitude to 10−22 m2 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 m2 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.

Table 8.

Material properties for the damage model

MaterialParameterValue
Solid soil matrix (Bentonite)Compressive strength, fc (MPa)1
Residual compressive strength, fcr (MPa)0.09
Strain at compressive strength, εc00.0003
Tensile strength, ft (MPa)−1
Residual tensile strength, ftr(MPa)−0.09
Strain at tensile strength, εt0−0.0003
Maximum permeability at maximum damage, kmax (m2) 1E−18
c, smoothing coefficient2
MaterialParameterValue
Solid soil matrix (Bentonite)Compressive strength, fc (MPa)1
Residual compressive strength, fcr (MPa)0.09
Strain at compressive strength, εc00.0003
Tensile strength, ft (MPa)−1
Residual tensile strength, ftr(MPa)−0.09
Strain at tensile strength, εt0−0.0003
Maximum permeability at maximum damage, kmax (m2) 1E−18
c, smoothing coefficient2

The results of Scenario S11, investigating the inclusion of a damage model, are presented below. Figures 2022 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, kmax, which was calibrated to a value of 1E-18 m2 to match the peak experimental inflow.

Fig. 20.

Effect of damage on gas (a) inflow and (b) outflow.

Fig. 20.

Effect of damage on gas (a) inflow and (b) outflow.

Fig. 21.

Effect of damage on the porefluid pressure profiles.

Fig. 21.

Effect of damage on the porefluid pressure profiles.

Fig. 22.

Effect of damage on total radial stresses.

Fig. 22.

Effect of damage on total radial stresses.

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

Fig. 23.

S11 results: gas migration through bentonite MX-80 sample for the 2D-axisymmetric model at times (a) 0 days, (b) 65 days and (c) 120 days.

Fig. 23.

S11 results: gas migration through bentonite MX-80 sample for the 2D-axisymmetric model at times (a) 0 days, (b) 65 days and (c) 120 days.

Although inclusion of damage to the model provides much better agreement with the experimental results, key features of the experimental gas outflow results depicting dilatancy and self-healing were not reproduced. There is, therefore, still a need to investigate additional mechanisms to model dilatancy-controlled gas flow. Such mechanisms will be investigated in future studies. These will include the consideration of spatially weighted heterogeneity for porosity and a poro-elasto-plastic model to simulate flow through preferential pathways. Inclusion of a swelling stress term will also be added to investigate the swelling behaviour in expansive soils. Finally, an investigation into constitutive relations to simulate the self-healing behaviour of swelling soils will be included.

S12: Effect of mesh size

A major finding of the study was the effect of mesh size on the shape of the porefluid pressure profile curves and the gas inflow and outflow curves. This was initially observed when moving from a 2D-axisymetric geometry to a 3D geometry, where the number of degrees of freedom needed to solve the model significantly increases. As a result of the increased computational expense, the 3D model was not able to complete a run at higher mesh sizes, using the fully coupled, or segregated solvers in COMSOL®. To exemplify this, the number of elements and degrees of freedom for each geometry are provided in Table 9.

Table 9.

Number of elements and degrees of freedom for 2D-axisymmetric and 3D models

Mesh qualifierNo. of elementsDegrees of freedom
2D axisymmetric HM
Extremely fine633683 632
Extra fine177822 094
Finer6027254
Normal2142450
3D HM
Extra fine81 0671 278 368
Finer21 553333 479
Normal400959 086
Mesh qualifierNo. of elementsDegrees of freedom
2D axisymmetric HM
Extremely fine633683 632
Extra fine177822 094
Finer6027254
Normal2142450
3D HM
Extra fine81 0671 278 368
Finer21 553333 479
Normal400959 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.

Fig. 24.

Effect of mesh size on porefluid pressure profiles.

Fig. 24.

Effect of mesh size on porefluid pressure profiles.

A similar comparison can be made between the 2D-axisymmetric solution with an extremely fine mesh size and the 3D solution with a finer mesh size as depicted in Figure 12(a). Due to the computational cost of modelling in 3D this poses a number of challenges. For one, preferential flow pathways cannot be adequately modelled using a 2D-axisymmetric geometry, as the microfractures cannot intersect the axis of symmetry. Therefore, in order to model dilatancy-controlled gas flow in heterogeneous specimens and propagate preferential flow pathways, it will be crucial to model in 3D. To address this current limitation, future studies will assess the use of different solvers and methods in COMSOL®, which may include the incorporation of discrete models, and identify ways to optimize the model for improved computational efficiency.

Conclusions

An important component in the design and long-term safety assessment of a DGR is the long-term performance of bentonite seals as barriers against gas migration in order to ensure the future safety to human health and the environment as a result of the exposure to these radionuclides. This paper proposes a hydro-mechanical linear poro-elastic visco-capillary mathematical model for 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.

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

Funding

This work was supported by the Canadian Nuclear Safety Commission.

References

1.
Ahusborde
,
E.
,
Amaziane
,
B.
&
Jurak
,
M.
2015
.
Three-dimensional numerical simulation by upscaling of gas migration through engineered and geological barriers for a deep repository for radioactive waste
. In:
Shaw
,
R.P.
(ed.)
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
123
141
, .
2.
Bennett
,
D.P.
,
Cuss
,
R.J.
,
Vardon
,
P.J.
,
Harrington
,
J.F.
,
Sedighi
,
M.
and
Thomas
,
H.R.
2015
.
Exploratory data analysis of the Large Scale Gas Injection Test (Lasgit) dataset, focusing on ‘second-order’ events around macro-scale gas flows
. In:
Shaw
,
R.P.
(ed.)
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
225
239
, .
3.
Carman
,
P.C.
1956
.
Flow of Gases Through Porous Media
 .
Academic Press Inc
,
New York
.
4.
CNSC
2006
.
Regulatory Guide G-320, Assessing the Long Term Safety of Radioactive Waste Management
 .
Canadian Nuclear Safety Commission
,
Ottawa
.
5.
Cuss
,
R.
,
Harrington
,
J.
,
Giot
,
R.
&
Auvray
,
C.
2014
.
Experimental observations of mechanical dilation at the onset of gas flow in Callovo-Oxfordian claystone
. In:
Norris
,
S., Bruno, J.
et al
(eds)
Clays in Natural and Engineered Barriers for Radioactive Waste Confinement
 .
Geological Society, London, Special Publications
400
,
507
519
, .
6.
Daniels
,
K.A.
&
Harrington
,
J.F.
2017
.
The response of compact bentonite during a 1D gas flow test
 .
British Geological Survey Open Report, OR/17/067
,
British Geological Survey
,
s.l
.
7.
Davis
,
R.O.
&
Selvadurai
,
A.P.
2002
.
Plasticity and Geomechanics
 .
Cambridge University Press
,
Cambridge
.
8.
Fall
,
M.
,
Nasir
,
O.
&
Nguyen
,
T.
2014
.
A coupled hydro-mechanical model for simulation of gas migration in host sedimentary rocks for nuclear waste repositories
.
Engineering Geology
 ,
176
,
24
44
.
9.
Graham
,
C.C.
,
Harrington
,
J.F.
,
Cuss
,
R.J.
&
Sellin
,
P.
2012
.
Gas migration experiments in bentonite: implications for numerical modelling
.
Mineralogical Magazine
 ,
76
,
3279
3292
.
10.
Harrington
,
J.F.
&
Horseman
,
S.T.
1999
.
Gas transport properties of clays and mudrocks
. In:
Aplin
,
A.C.
,
Fleet
,
A.J.
&
MacQuaker
,
J.H.
(eds)
Muds and Mudstones: Physical and Fluid Flow Properties
 .
Geological Society of London, London, Special Publications
,
158
,
107
124
, .
11.
Harrington
,
J.F.
,
de la Vaissière
,
R.
,
Noy
,
D.J.
,
Cuss
,
R.J.
&
Talandier
,
J.
2012
.
Gas flow in Callovo-Oxfordian claystone (COx): results from laboratory and field-scale measurements
.
Mineralogical Magazine
 ,
76
,
3303
3318
.
12.
Harrington
,
J.F.
,
Graham
,
C.C.
,
Cuss
,
R.J.
&
Norris
,
S.
2017
.
Gas network development in a precompacted bentonite experiment: evidence of generation and evolution
.
Applied Clay Science
 ,
147
,
80
89
.
13.
Ho
,
C.K.
&
Webb
,
S.W.
(eds)
2006
.
Gas Transport in Porous Media
 .
Springer
,
Dordrecht
.
14.
Hoch
,
A.R.
,
Cliffe
,
K.A.
,
Swift
,
B.T.
&
Rodwell
,
W.R.
2004
.
Modelling Gas Migrationin Compacted Bentonite: GAMBIT Club Phase 3
 .
Final Report
,
POSIVA
,
Olkiluoto, Finland
.
15.
Horseman
,
S.T.
,
Harrington
,
J.F.
&
Sellin
,
P.
1999
.
Gas migration in clay barriers
.
Engineering Geology
 ,
54
,
139
149
.
16.
Horseman
,
S.T.
,
Harrington
,
J.F.
&
Sellin
,
P.
2003
.
Water and Gas Movement in Mx80 Bentonite Buffer Clay
.
Symposium on the Scientific Basis for Nuclear Waste Management XXVII (Kalmar), Materials and Research Society
,
15–19 June 2003
,
Kalmar, Sweden
,
807
,
715
720
.
17.
Huang
,
S.
,
Barbour
,
S.L.
&
Fredlund
,
D.G.
1998
.
Development and verification of a coefficient of permeability function for a deformable unsaturated soil
.
Canadian Geotechnical Journal
 ,
35
,
411
425
.
18.
Khalili
,
N.
&
Khabbaz
,
M.H.
1998
.
A unique relatiship for χ for the determination of the shear strength of unsaturated soils
.
Géotechnique
 ,
48
,
681
687
.
19.
Laliberte
,
G.E.
,
Corey
,
A.T.
&
Brooks
,
R.H.
1966
.
Properties of Unsaturated Porous Media, Hydrology Paper No. 17
 .
Colorado State University
,
Fort Collins, CO
.
20.
Low-Level Radioactive Waste Management Office
.
2012
.
Inventory of Radioactive Waste in Canada
 .
Low-Level Radioactive Waste Management Office
,
Ottawa
.
21.
Man
,
A.
&
Martino
,
J.B.
2009
.
NWMO TR-2009-20: Thermal, Hydraulic and Mechanical Properties of Sealing Materials
 .
Atomic Energy of Canada Limited
,
Ottawa
.
22.
Marschall
,
P.
,
Horseman
,
S.
&
Gimmi
,
T.
2005
.
Characterisation of gas transport properties of the opalinus clay, a potential host rock formation for radioactive waste disposal
.
Oil and Gas Science and Technology – Rev. IFP
 ,
60
,
121
139
.
23.
Nash
,
P.J.
,
Swift
,
B.T.
,
Goodfield
,
M.
&
Rodwell
,
W.R.
1998
.
Modelling Gas Migration in Compacted Bentonite: A report produced for the GAMBIT Club
 .
POSIVA
,
Helsinki
.
24.
Nguyen
,
T.
&
Le
,
A.
2015
.
Simultaneous gas and water flow in a damage-susceptible bedded argillaceous rock
.
Canadian Geotechnical Journal
 ,
52
,
18
32
.
25.
Norris
,
S.
2015
.
EC FORGE project: updated consideration of gas generation and migration in the safety case
. In:
Shaw
,
R.P.
(ed.)
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
241
258
, .
26.
Pall
,
R.
&
Moshenin
,
N.N.
1980
.
Permeability of porous media as a function of porosity and particle size distribution
.
TRANSACTIONS of the ASAE
 ,
23
,
742
745
.
27.
Sellin
,
P.
(ed.)
2014
.
Experiments and Modelling on the behaviour of EBS
 .
FORGE Reprot, D3.38
. :
s.n
,
s.l
.
28.
Sellin
,
P.
&
Leupin
,
O.X.
2013
.
The use of clay as an engineered barrier in radioactive-waste management – a review
.
Clays and Clay Minerals
 ,
61
,
477
498
.
29.
Shaw
,
R.P.
(ed.)
2015a
.
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
30.
Shaw
,
R.P.
2015b
.
The Fate of Repository Gases (FORGE) project
. In:
Shaw
,
R.P.
(ed.)
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
1
7
, .
31.
Tang
,
C.A.
,
Tham
,
L.G.
,
Lee
,
P.K.K.
,
Yang
,
T.H.
and
Li
,
L.C.
et al
2002
.
Coupled analysis of flow, stress and damage (FSD) in rock failure
.
International Journal of Rock Mechanics & Mining Sciences
 ,
39
,
477
489
.
32.
Tripathy
,
S.
,
Tadza
,
M.Y.
&
Thomas
,
H.R.
2014
.
Soil-water characteristic curves of clays
.
Canadian Geotechnical Journal
 ,
51
,
869
883
.
33.
Valdes-Parada
,
F.J.
,
Ochoa-Tapia
,
J.A.
&
Alvarez-Ramirez
,
J.
2009
.
Validity of the permeability Carman-Kozeny equation: a volume averaging approach
.
Physica A
 ,
388
,
789
798
.
34.
van Genuchten
,
M.T.
1980
.
A closed-form equation for predicting the hydraulic conductivity of unsaturated soils
.
Soil Science Society of America Journal
 ,
44
,
892
898
.
35.
Villar
,
M.
2004
.
MX-80 Bentonite: Thermo-Hydro-Mechanical Characterisation Performed at CIEMAT in the Context of the Prototype Project
.
36.
Villar
,
M.V.
&
Lloret
,
A.
2008
.
Influence of dry density and water content on the swelling of a compacted bentonite
.
Applied Clay Science
 ,
39
,
38
49
.

Figures & Tables

Fig. 1.

Four main processes of gas migration in clays (Cuss et al. 2014).

Fig. 1.

Four main processes of gas migration in clays (Cuss et al. 2014).

Fig. 2.

Control volume and phase relations.

Fig. 2.

Control volume and phase relations.

Fig. 3.

2D-axisymmetric geometry and meshing.

Fig. 3.

2D-axisymmetric geometry and meshing.

Fig. 4.

3D geometry and meshing.

Fig. 4.

3D geometry and meshing.

Fig. 5.

Dirichlet BCs for the (a) lower gas injection pressure and upper water backpressure (b) lower and upper gas concentrations in porewater.

Fig. 5.

Dirichlet BCs for the (a) lower gas injection pressure and upper water backpressure (b) lower and upper gas concentrations in porewater.

Fig. 6.

Intrinsic permeability (a), effectivity diffusivity (b) and AEV (c) curves as a function of porosity.

Fig. 6.

Intrinsic permeability (a), effectivity diffusivity (b) and AEV (c) curves as a function of porosity.

Fig. 7.

χ as a function of both AEV and suction.

Fig. 7.

χ as a function of both AEV and suction.

Fig. 8.

Soil water and soil gas characteristic curves for bentonite at n = 0.44.

Fig. 8.

Soil water and soil gas characteristic curves for bentonite at n = 0.44.

Fig. 9.

Relative permeability curves (water and gas) in bentonite v. suction at n = 0.44.

Fig. 9.

Relative permeability curves (water and gas) in bentonite v. suction at n = 0.44.

Fig. 10.

S1 results: gas migration through bentonite MX-80 sample. For the 3D model, at times (a) 0 days, (b) 65 days and (c) 120 days; and for the 2D-axisymmetric model, at times (d) 0 days, (e) 65 days and (f) 120 days.

Fig. 10.

S1 results: gas migration through bentonite MX-80 sample. For the 3D model, at times (a) 0 days, (b) 65 days and (c) 120 days; and for the 2D-axisymmetric model, at times (d) 0 days, (e) 65 days and (f) 120 days.

Fig. 11.

S1 results: gas (a) inflow and (b) outflow profiles over time.

Fig. 11.

S1 results: gas (a) inflow and (b) outflow profiles over time.

Fig. 12.

S1 results: porefluid pressure profiles over (a) time and (b) axial length of the specimen.

Fig. 12.

S1 results: porefluid pressure profiles over (a) time and (b) axial length of the specimen.

Fig. 13.

S1 results – Total radial stresses over time.

Fig. 13.

S1 results – Total radial stresses over time.

Fig. 14.

S2–S4 results: (a) contribution of advection, dissolution and diffusion to flow behaviour and (b) diffusion gradient over axial length of the soil specimen.

Fig. 14.

S2–S4 results: (a) contribution of advection, dissolution and diffusion to flow behaviour and (b) diffusion gradient over axial length of the soil specimen.

Fig. 15.

S2–S4 results: contribution of (a) advection, dissolution and diffusion and (b) elastic deformation to porefluid pressure profiles.

Fig. 15.

S2–S4 results: contribution of (a) advection, dissolution and diffusion and (b) elastic deformation to porefluid pressure profiles.

Fig. 16.

Change in porosity over time as a result of elastic deformation.

Fig. 16.

Change in porosity over time as a result of elastic deformation.

Fig. 17.

S5–S6 results: effect of changes in initial porosity on (a) porefluid pressure and (b) gas inflow.

Fig. 17.

S5–S6 results: effect of changes in initial porosity on (a) porefluid pressure and (b) gas inflow.

Fig. 18.

S7–S8 results: effect of changes in AEV on (a) porefluid pressure and (b) gas inflow.

Fig. 18.

S7–S8 results: effect of changes in AEV on (a) porefluid pressure and (b) gas inflow.

Fig. 19.

S9–S10 results: effect of changes in kij on (a) porefluid pressure and (b) gas inflow.

Fig. 19.

S9–S10 results: effect of changes in kij on (a) porefluid pressure and (b) gas inflow.

Fig. 20.

Effect of damage on gas (a) inflow and (b) outflow.

Fig. 20.

Effect of damage on gas (a) inflow and (b) outflow.

Fig. 21.

Effect of damage on the porefluid pressure profiles.

Fig. 21.

Effect of damage on the porefluid pressure profiles.

Fig. 22.

Effect of damage on total radial stresses.

Fig. 22.

Effect of damage on total radial stresses.

Fig. 23.

S11 results: gas migration through bentonite MX-80 sample for the 2D-axisymmetric model at times (a) 0 days, (b) 65 days and (c) 120 days.

Fig. 23.

S11 results: gas migration through bentonite MX-80 sample for the 2D-axisymmetric model at times (a) 0 days, (b) 65 days and (c) 120 days.

Fig. 24.

Effect of mesh size on porefluid pressure profiles.

Fig. 24.

Effect of mesh size on porefluid pressure profiles.

Table 1.

Canadian Waste Inventory Projections to 2011 and 2050 (adapted fromLLRWMO 2012)

Waste categoryWaste inventory to the end of 2011Waste inventory to the end of 2050
High-level radioactive waste (nuclear fuel waste)9400 m320 000 m3
Intermediate-level radioactive waste33 400 m367 000 m3
Low-level radioactive waste2 343 000 m32 594 000 m3
Waste categoryWaste inventory to the end of 2011Waste inventory to the end of 2050
High-level radioactive waste (nuclear fuel waste)9400 m320 000 m3
Intermediate-level radioactive waste33 400 m367 000 m3
Low-level radioactive waste2 343 000 m32 594 000 m3
Table 2.

Model dimensions and location of the porefluid arrays

Sensor nameMeasurement typeAxial distance from injection face (mm)
Injection load cellTotal stress0
Radial load cell 1Total stress15.2
Radial load cell 2Total stress60
Radial load cell 3Total stress104.8
Backpressure load cellTotal stress120
Radial porewater array 1Porefluid pressure38.6
Radial porewater array 2Porefluid pressure60
Radial porewater array 3Porefluid pressure81.4
Central filter (Middle)Porefluid pressure60
Sensor nameMeasurement typeAxial distance from injection face (mm)
Injection load cellTotal stress0
Radial load cell 1Total stress15.2
Radial load cell 2Total stress60
Radial load cell 3Total stress104.8
Backpressure load cellTotal stress120
Radial porewater array 1Porefluid pressure38.6
Radial porewater array 2Porefluid pressure60
Radial porewater array 3Porefluid pressure81.4
Central filter (Middle)Porefluid pressure60
Table 3.

Base study and parametric study model scenarios

Scenario numberModelDescriptionPorosityIntrinsic permeability (m2)AEV (Pa)SWCC fitting parameters
a′ (1/m)nm
S1-ValidationHMAdvection, dissolution and diffusion linear poro-elasticity0.443.4E−218.0E + 060.00062.040.51
Parametric studies
S2HAdvection, no dissolution, no diffusion0.443.4E−218.0E + 060.00062.040.51
S3HAdvection, dissolution, no diffusion0.443.4E−218.0E + 060.00062.040.51
S4HAdvection, dissolution, diffusion0.443.4E−218.0E + 060.00062.040.51
Sensitivity analysis
S5HM↑ por0.485.1E−216.2E + 060.00082.040.51
S6HM↓ por0.306.9E−221.6E + 070.00032.040.51
S7HM↑ AEV0.443.4E−212.0E + 070.000252.040.51
S8HM↓ AEV0.443.4E−215.0E + 060.00062.040.51
S9HMkij0.441.0E−198.0E + 060.00062.040.51
S10HMkij0.441.0E−228.0E + 060.00062.040.51
S11HMInclusion of a damage model0.443.4E−218.0E + 060.00062.040.51
S12HMEffect of mesh size0.443.4E−218.0E + 060.00062.040.51
Scenario numberModelDescriptionPorosityIntrinsic permeability (m2)AEV (Pa)SWCC fitting parameters
a′ (1/m)nm
S1-ValidationHMAdvection, dissolution and diffusion linear poro-elasticity0.443.4E−218.0E + 060.00062.040.51
Parametric studies
S2HAdvection, no dissolution, no diffusion0.443.4E−218.0E + 060.00062.040.51
S3HAdvection, dissolution, no diffusion0.443.4E−218.0E + 060.00062.040.51
S4HAdvection, dissolution, diffusion0.443.4E−218.0E + 060.00062.040.51
Sensitivity analysis
S5HM↑ por0.485.1E−216.2E + 060.00082.040.51
S6HM↓ por0.306.9E−221.6E + 070.00032.040.51
S7HM↑ AEV0.443.4E−212.0E + 070.000252.040.51
S8HM↓ AEV0.443.4E−215.0E + 060.00062.040.51
S9HMkij0.441.0E−198.0E + 060.00062.040.51
S10HMkij0.441.0E−228.0E + 060.00062.040.51
S11HMInclusion of a damage model0.443.4E−218.0E + 060.00062.040.51
S12HMEffect of mesh size0.443.4E−218.0E + 060.00062.040.51
Table 4.

Number of elements and degrees of freedom

Model geometryMesh qualifierNo. of elementsDegrees of freedom
2D axisymmetricExtremely fine8886117 132
3DNormal10 698161 741
Model geometryMesh qualifierNo. of elementsDegrees of freedom
2D axisymmetricExtremely fine8886117 132
3DNormal10 698161 741
Table 5.

Material properties of the solid soil matrix and gas and liquid phases

MaterialParameterValue
Solid soil matrix (Bentonite)Saturated/intrinsic permeability, kij (m2)3.40E−21
Dry density, ρd (kg m−3)1560
Maximum degree of saturation, Sws1
Residual degree of Saturation, Swr0.05
Porosity, n4.40E−01
Young’s Modulus, E (MPa)307
Poisson ratio, ν0.4
Bulk Modulus of bentonite, Ks (Pa)2.0E + 08
Compression Index, 1/Ks(Pa−1)5.00E−09
Volume-surface mean diameter, Dvs (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, Kg (Pa)pg
Henry’s Coefficient, H0.0091
Diffusivity of helium in water (m2 s−1)6.29E−09
Diameter of helium gas particle, dhe (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, Kw (Pa)2.2E + 09
Compression Index, 1/Kw (Pa−1)4.55E−10
MaterialParameterValue
Solid soil matrix (Bentonite)Saturated/intrinsic permeability, kij (m2)3.40E−21
Dry density, ρd (kg m−3)1560
Maximum degree of saturation, Sws1
Residual degree of Saturation, Swr0.05
Porosity, n4.40E−01
Young’s Modulus, E (MPa)307
Poisson ratio, ν0.4
Bulk Modulus of bentonite, Ks (Pa)2.0E + 08
Compression Index, 1/Ks(Pa−1)5.00E−09
Volume-surface mean diameter, Dvs (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, Kg (Pa)pg
Henry’s Coefficient, H0.0091
Diffusivity of helium in water (m2 s−1)6.29E−09
Diameter of helium gas particle, dhe (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, Kw (Pa)2.2E + 09
Compression Index, 1/Kw (Pa−1)4.55E−10
Table 6.

Initial value conditions

Parameter nameParameterInitial value condition
Initial poregas pressurepginitial1.01E + 05 Pa
Initial degree of saturationSwinitial0.96
Initial suction (from SWCC)sinitial5.8E + 06 Pa
Initial porewater pressurepwinitial−5.7E + 06 Pa
Initial displacement fieldui0 m
Initial stressσ0xx = σ0yy = σ0zz5.0E + 05 Pa
Initial gas concentration in porewater @STPCg,H2O0.073 mol m−3
Parameter nameParameterInitial value condition
Initial poregas pressurepginitial1.01E + 05 Pa
Initial degree of saturationSwinitial0.96
Initial suction (from SWCC)sinitial5.8E + 06 Pa
Initial porewater pressurepwinitial−5.7E + 06 Pa
Initial displacement fieldui0 m
Initial stressσ0xx = σ0yy = σ0zz5.0E + 05 Pa
Initial gas concentration in porewater @STPCg,H2O0.073 mol m−3
Table 7.

van Genuchten-Mualem SWCC and relative permeability function fitting parameters

ParameterValue
SWCC fitting parameter, a′(1m)0.00065*
SWCC fitting parameter, n2.04 (Villar 2004)
SWCC fitting parameter, m0.51 (Villar 2004)
Relative permeability function fitting parameter, L0.51
ParameterValue
SWCC fitting parameter, a′(1m)0.00065*
SWCC fitting parameter, n2.04 (Villar 2004)
SWCC fitting parameter, m0.51 (Villar 2004)
Relative permeability function fitting parameter, L0.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.

Table 8.

Material properties for the damage model

MaterialParameterValue
Solid soil matrix (Bentonite)Compressive strength, fc (MPa)1
Residual compressive strength, fcr (MPa)0.09
Strain at compressive strength, εc00.0003
Tensile strength, ft (MPa)−1
Residual tensile strength, ftr(MPa)−0.09
Strain at tensile strength, εt0−0.0003
Maximum permeability at maximum damage, kmax (m2) 1E−18
c, smoothing coefficient2
MaterialParameterValue
Solid soil matrix (Bentonite)Compressive strength, fc (MPa)1
Residual compressive strength, fcr (MPa)0.09
Strain at compressive strength, εc00.0003
Tensile strength, ft (MPa)−1
Residual tensile strength, ftr(MPa)−0.09
Strain at tensile strength, εt0−0.0003
Maximum permeability at maximum damage, kmax (m2) 1E−18
c, smoothing coefficient2
Table 9.

Number of elements and degrees of freedom for 2D-axisymmetric and 3D models

Mesh qualifierNo. of elementsDegrees of freedom
2D axisymmetric HM
Extremely fine633683 632
Extra fine177822 094
Finer6027254
Normal2142450
3D HM
Extra fine81 0671 278 368
Finer21 553333 479
Normal400959 086
Mesh qualifierNo. of elementsDegrees of freedom
2D axisymmetric HM
Extremely fine633683 632
Extra fine177822 094
Finer6027254
Normal2142450
3D HM
Extra fine81 0671 278 368
Finer21 553333 479
Normal400959 086

Contents

References

References

1.
Ahusborde
,
E.
,
Amaziane
,
B.
&
Jurak
,
M.
2015
.
Three-dimensional numerical simulation by upscaling of gas migration through engineered and geological barriers for a deep repository for radioactive waste
. In:
Shaw
,
R.P.
(ed.)
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
123
141
, .
2.
Bennett
,
D.P.
,
Cuss
,
R.J.
,
Vardon
,
P.J.
,
Harrington
,
J.F.
,
Sedighi
,
M.
and
Thomas
,
H.R.
2015
.
Exploratory data analysis of the Large Scale Gas Injection Test (Lasgit) dataset, focusing on ‘second-order’ events around macro-scale gas flows
. In:
Shaw
,
R.P.
(ed.)
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
225
239
, .
3.
Carman
,
P.C.
1956
.
Flow of Gases Through Porous Media
 .
Academic Press Inc
,
New York
.
4.
CNSC
2006
.
Regulatory Guide G-320, Assessing the Long Term Safety of Radioactive Waste Management
 .
Canadian Nuclear Safety Commission
,
Ottawa
.
5.
Cuss
,
R.
,
Harrington
,
J.
,
Giot
,
R.
&
Auvray
,
C.
2014
.
Experimental observations of mechanical dilation at the onset of gas flow in Callovo-Oxfordian claystone
. In:
Norris
,
S., Bruno, J.
et al
(eds)
Clays in Natural and Engineered Barriers for Radioactive Waste Confinement
 .
Geological Society, London, Special Publications
400
,
507
519
, .
6.
Daniels
,
K.A.
&
Harrington
,
J.F.
2017
.
The response of compact bentonite during a 1D gas flow test
 .
British Geological Survey Open Report, OR/17/067
,
British Geological Survey
,
s.l
.
7.
Davis
,
R.O.
&
Selvadurai
,
A.P.
2002
.
Plasticity and Geomechanics
 .
Cambridge University Press
,
Cambridge
.
8.
Fall
,
M.
,
Nasir
,
O.
&
Nguyen
,
T.
2014
.
A coupled hydro-mechanical model for simulation of gas migration in host sedimentary rocks for nuclear waste repositories
.
Engineering Geology
 ,
176
,
24
44
.
9.
Graham
,
C.C.
,
Harrington
,
J.F.
,
Cuss
,
R.J.
&
Sellin
,
P.
2012
.
Gas migration experiments in bentonite: implications for numerical modelling
.
Mineralogical Magazine
 ,
76
,
3279
3292
.
10.
Harrington
,
J.F.
&
Horseman
,
S.T.
1999
.
Gas transport properties of clays and mudrocks
. In:
Aplin
,
A.C.
,
Fleet
,
A.J.
&
MacQuaker
,
J.H.
(eds)
Muds and Mudstones: Physical and Fluid Flow Properties
 .
Geological Society of London, London, Special Publications
,
158
,
107
124
, .
11.
Harrington
,
J.F.
,
de la Vaissière
,
R.
,
Noy
,
D.J.
,
Cuss
,
R.J.
&
Talandier
,
J.
2012
.
Gas flow in Callovo-Oxfordian claystone (COx): results from laboratory and field-scale measurements
.
Mineralogical Magazine
 ,
76
,
3303
3318
.
12.
Harrington
,
J.F.
,
Graham
,
C.C.
,
Cuss
,
R.J.
&
Norris
,
S.
2017
.
Gas network development in a precompacted bentonite experiment: evidence of generation and evolution
.
Applied Clay Science
 ,
147
,
80
89
.
13.
Ho
,
C.K.
&
Webb
,
S.W.
(eds)
2006
.
Gas Transport in Porous Media
 .
Springer
,
Dordrecht
.
14.
Hoch
,
A.R.
,
Cliffe
,
K.A.
,
Swift
,
B.T.
&
Rodwell
,
W.R.
2004
.
Modelling Gas Migrationin Compacted Bentonite: GAMBIT Club Phase 3
 .
Final Report
,
POSIVA
,
Olkiluoto, Finland
.
15.
Horseman
,
S.T.
,
Harrington
,
J.F.
&
Sellin
,
P.
1999
.
Gas migration in clay barriers
.
Engineering Geology
 ,
54
,
139
149
.
16.
Horseman
,
S.T.
,
Harrington
,
J.F.
&
Sellin
,
P.
2003
.
Water and Gas Movement in Mx80 Bentonite Buffer Clay
.
Symposium on the Scientific Basis for Nuclear Waste Management XXVII (Kalmar), Materials and Research Society
,
15–19 June 2003
,
Kalmar, Sweden
,
807
,
715
720
.
17.
Huang
,
S.
,
Barbour
,
S.L.
&
Fredlund
,
D.G.
1998
.
Development and verification of a coefficient of permeability function for a deformable unsaturated soil
.
Canadian Geotechnical Journal
 ,
35
,
411
425
.
18.
Khalili
,
N.
&
Khabbaz
,
M.H.
1998
.
A unique relatiship for χ for the determination of the shear strength of unsaturated soils
.
Géotechnique
 ,
48
,
681
687
.
19.
Laliberte
,
G.E.
,
Corey
,
A.T.
&
Brooks
,
R.H.
1966
.
Properties of Unsaturated Porous Media, Hydrology Paper No. 17
 .
Colorado State University
,
Fort Collins, CO
.
20.
Low-Level Radioactive Waste Management Office
.
2012
.
Inventory of Radioactive Waste in Canada
 .
Low-Level Radioactive Waste Management Office
,
Ottawa
.
21.
Man
,
A.
&
Martino
,
J.B.
2009
.
NWMO TR-2009-20: Thermal, Hydraulic and Mechanical Properties of Sealing Materials
 .
Atomic Energy of Canada Limited
,
Ottawa
.
22.
Marschall
,
P.
,
Horseman
,
S.
&
Gimmi
,
T.
2005
.
Characterisation of gas transport properties of the opalinus clay, a potential host rock formation for radioactive waste disposal
.
Oil and Gas Science and Technology – Rev. IFP
 ,
60
,
121
139
.
23.
Nash
,
P.J.
,
Swift
,
B.T.
,
Goodfield
,
M.
&
Rodwell
,
W.R.
1998
.
Modelling Gas Migration in Compacted Bentonite: A report produced for the GAMBIT Club
 .
POSIVA
,
Helsinki
.
24.
Nguyen
,
T.
&
Le
,
A.
2015
.
Simultaneous gas and water flow in a damage-susceptible bedded argillaceous rock
.
Canadian Geotechnical Journal
 ,
52
,
18
32
.
25.
Norris
,
S.
2015
.
EC FORGE project: updated consideration of gas generation and migration in the safety case
. In:
Shaw
,
R.P.
(ed.)
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
241
258
, .
26.
Pall
,
R.
&
Moshenin
,
N.N.
1980
.
Permeability of porous media as a function of porosity and particle size distribution
.
TRANSACTIONS of the ASAE
 ,
23
,
742
745
.
27.
Sellin
,
P.
(ed.)
2014
.
Experiments and Modelling on the behaviour of EBS
 .
FORGE Reprot, D3.38
. :
s.n
,
s.l
.
28.
Sellin
,
P.
&
Leupin
,
O.X.
2013
.
The use of clay as an engineered barrier in radioactive-waste management – a review
.
Clays and Clay Minerals
 ,
61
,
477
498
.
29.
Shaw
,
R.P.
(ed.)
2015a
.
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
30.
Shaw
,
R.P.
2015b
.
The Fate of Repository Gases (FORGE) project
. In:
Shaw
,
R.P.
(ed.)
Gas Generation and Migration in Deep Geological Radioactive Waste Repositories
 .
Geological Society, London, Special Publications
,
415
,
1
7
, .
31.
Tang
,
C.A.
,
Tham
,
L.G.
,
Lee
,
P.K.K.
,
Yang
,
T.H.
and
Li
,
L.C.
et al
2002
.
Coupled analysis of flow, stress and damage (FSD) in rock failure
.
International Journal of Rock Mechanics & Mining Sciences
 ,
39
,
477
489
.
32.
Tripathy
,
S.
,
Tadza
,
M.Y.
&
Thomas
,
H.R.
2014
.
Soil-water characteristic curves of clays
.
Canadian Geotechnical Journal
 ,
51
,
869
883
.
33.
Valdes-Parada
,
F.J.
,
Ochoa-Tapia
,
J.A.
&
Alvarez-Ramirez
,
J.
2009
.
Validity of the permeability Carman-Kozeny equation: a volume averaging approach
.
Physica A
 ,
388
,
789
798
.
34.
van Genuchten
,
M.T.
1980
.
A closed-form equation for predicting the hydraulic conductivity of unsaturated soils
.
Soil Science Society of America Journal
 ,
44
,
892
898
.
35.
Villar
,
M.
2004
.
MX-80 Bentonite: Thermo-Hydro-Mechanical Characterisation Performed at CIEMAT in the Context of the Prototype Project
.
36.
Villar
,
M.V.
&
Lloret
,
A.
2008
.
Influence of dry density and water content on the swelling of a compacted bentonite
.
Applied Clay Science
 ,
39
,
38
49
.

Related

Citing Books via

Close Modal

or Create an Account

Close Modal
Close Modal