A Mathematical Model of Gas and Water Flow in a Swelling Geomaterial—Part 2. Process Simulation

: Gases can potentially generate in a deep geological repository (DGR) for the long-term containment of radioactive waste. Natural and engineered barriers provide containment of the waste by mitigating contaminant migration. However, if gas pressures exceed the mechanical strength of these barriers, preferential ﬂow pathways for both the gases and the porewater could form, providing a source of potential exposure to people and the environment. Expansive soils, such as bentonite-based materials, are widely considered as sealing materials. Understanding the long-term performance of these seals as barriers against gas migration is an important component in the design and the long-term safety assessment of a DGR. This study proposes a hydro-mechanical mathematical model for migration of gas through a low-permeable swelling geomaterial based on the theoretical framework of poromechanics. Using the ﬁnite element method, the model is used to simulate 1D ﬂow through a conﬁned cylindrical sample of near-saturated low-permeable soil under a constant volume boundary stress condition. The study expands upon previous work by the authors by assessing the inﬂuence of heterogeneity, the Klinkenberg “slip ﬂow” e ﬀ ect, and a swelling stress on ﬂow behavior. Based on the results, this study provides fundamental insight into a number of factors that may inﬂuence two-phase ﬂow. permeability of the gas increases significantly with a decrease in poregas pressure, and as poregas pressures increase, the intrinsic permeability of the gas approaches that of water. As there is little information in literature regarding the slip factor for helium in bentonite clays, in order to assess whether the Klinkenberg effect could provide a significant role in matching key features of the experimental results, a constant, c = 1.0 × 10 8 , was used for the Klinkenberg effect. A value of c = 1.0 × 10 4 , which matches that of the empirical equation, would not provide the significant change in permeability needed to obtain complete breakthrough into the sample. helium in bentonite clays, in order to assess whether the Klinkenberg effect could provide a significant role in matching key features of the experimental results, a constant, c = 1.0 × 10 8 , was used for the Klinkenberg effect. A value of c = 1.0 × 10 4 , which matches that of the empirical equation, would not provide the significant change in permeability needed to obtain complete breakthrough into the sample. Klinkenberg e ﬀ ect; intrinsic permeability of gas as a function poregas pressure, k g . results reinforce the conclusion that additional mechanisms of mechanical deformation are required to initiate dilation-controlled gas ﬂow and the formation of discrete preferential ﬂow pathways.


Introduction
In Canada, nuclear waste has been generated and accumulated since the 1930s when the Port Radium radium mine began operating in the Northwest Territories [1]. Over the past century, 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.
In identifying a long-term solution to manage global radioactive waste, the international community, including Canada, has been evaluating the use 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 to the environment and radiological exposure to people.
In developing a safety case for a DGR to provide supportive arguments that the management of radioactive waste will be protective of human health and the environment over the long term, relevant features, events, and processes (FEPs) must be evaluated [2,3]. One such process with the potential means for radiological exposure to the biosphere is the generation of radioactive gas that may migrate to the surface [4]. Gas could be generated through a number of processes, including degradation of organic matter, radioactive decay of the waste, corrosion of metals producing hydrogen gas (H 2 ), Several works have identified that a number of factors contribute to the effect of gas fingering. These include the displacement rate of the gas with the liquid (i.e., differences in mobility of the individual phases), the extent of heterogeneity of the porous medium and spatial variations of porosity and permeability, the viscosity ratio between gas and liquid, the three-dimensional dispersion within the porous media, and the aspect ratio and boundary conditions representing the narrowness of the porous media [21][22][23]. Additionally, in order to numerically simulate fingering effects accurately, consideration of mesh size also plays a critical role [22]. Moortgat et al. also identified that other factors may help stabilize gas fingering, which include consideration of gravity, diffusion, and capillary pressure [23].
Dagher et al. [11] developed a fully-coupled hydro-mechanical (HM) linear-elastic mathematical model for advective-diffusive visco-capillary controlled two-phase flow through geomaterials in order to model the first two transport mechanisms (i.e., advection and diffusion of dissolved gas and visco-capillary "two-phase" flow) proposed by Marschall et al. [10]. Results from a constant volume 1D flow experiment was used to validate the model. A number of parametric studies were investigated to assess the contribution of advection of poregas, diffusion of dissolved gas in porewater, advection of dissolved gas in porewater, and inclusion of mechanical deformation (linear elasticity) on flow behavior with increasing gas pressures over time. Additionally, sensitivity analyses were conducted to gain an understanding of the influence of a number of soil properties on flow behavior, such as the effect of modifying the air-entry value (AEV), intrinsic permeability, and initial porosity of the soil specimen. Finally, the study investigated the use of a linear elastic damage model to better represent the experimental results. Although the model results reproduced some of the general features noted in the experimental results, the model was unable to simulate dilatancy-controlled gas flow. The study concluded that, in order to simulate dilatancy-controlled gas flow, additional mechanisms need to be considered within the model. These include the use of advanced mechanism of mechanical deformation to be coupled to gas transport, consideration of heterogeneity within the soil sample to help induce preferential flow, inclusion of a swelling stress term to incorporate the swelling behavior of expansive soils, and the incorporation of a self-healing mechanism to represent observed phenomena of experimental studies [11].
This study builds upon the authors' previously published work by investigating the use of advanced physical relationships on flow behavior, specifically the influence of: (i) heterogeneity within the geomaterial to promote preferential flow, (ii) the Klinkenberg "slip flow" effect on gas permeability, and (iii) a linear swelling stress to promote swelling of the geomaterial. Finally, the authors investigate the conditions that promote fingering in two-phase flow. The models are validated against experimental data from a 1D flow test. The results of this study will be used to inform the development of a more sophisticated HM model that can be used to simulate dilation-controlled gas flow in a swelling geomaterial.
This research continues to be, in part, a contribution to 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). This task, led by the British Geological Survey (BGS), further attempts to identify the physical HM mechanisms required to adequately model dilatancy-controlled gas migration.

Study Overview
This study expands upon the work performed by Dagher et al. [11] on the development of a mathematical model for gas migration (two-phase flow) through a low-permeable swelling geomaterial. Using the theoretical framework of poromechanics, this study provides a process simulation and analysis of enhanced processes for two-phase flow related to the inclusion of the Minerals 2020, 10, 32 4 of 33 following advanced physical relations to the HM model to investigate their effect on two-phase (water and gas) flow behavior, specifically: i.
heterogeneity through the use of a random normal distribution of porosity; ii. Klinkenberg "slip flow" effect for the intrinsic permeability of the gas in soil; iii. a swelling strain to take into account volume changes of the geomaterial with suction.
The mathematical model presented in this paper follows the general formulation by Dagher et al. [11] and Nguyen and Le [24] and includes the constitutive relations for the soil water characteristic curves (SWCCs) through the application of the van Genuchten equation and the relative permeability-saturation relationships using Mualem's model. The applicable constitutive relations and governing equations for conservation of momentum, water mass, and gas mass are also discussed in the companion paper to this one [25]. For sake of consistency, the final governing equations are also provided below along with a mathematical description of the Klinkenberg effect and swelling strain.

Conservation of Water Mass
The governing equation for the conservation of water mass can be expressed by Equation (1), where ρ w is the density of water phase (kg m −3 ) p w is the porewater pressure (Pa) k ij,w is the intrinsic permeability tensor of water flow in the porous medium (m 2 ) k r,w is the relative permeability of the water phase (adimensional) µ w is the dynamic viscosity of the water phase (Pa s or kg m −1 s −1 ) g is the acceleration due to gravity (m s −2 ) n is the porosity (m 3 voids · m −3 total) S w is the degree of saturation of water (adimensional) s is the matric suction p g − p w (Pa) K w is the bulk modulus of the water phase (Pa s or kg m −1 s −1 ) ∂u k ∂x k = ε vol is the volumetric strain (adimensional) u is the displacement of the solid skeleton (m) t is time (s) Note that, in this study, the permeability is assumed to be isotropic, therefore k ij = k; however, we keep the tensorial notation in the governing equation for the sake of generalization.

Conservation of Gas Mass
The governing equation for the conservation of gas mass can be expressed by Equation (2), where ρ g is the density of the gas phase (kg m −3 ) p g is the poregas pressure (Pa) k ij,g is the intrinsic permeability tensor of gas flow in the porous medium (m 2 ) k r,g is the relative permeability of the gas phase (adimensional) µ g is the dynamic viscosity of the gas phase (Pa s or kg m −1 s −1 ) H is Henry's coefficient (kg species A m −3 in aqueous phase kg −1 species A m 3 in gas phase) D e is the effective diffusivity of gas dissolved in water through porous media (m 2 s −1 ) K g is the bulk modulus of the gas phase (Pa).

Conservation of Momentum (Quasi-Static Equilibrium)
For a linear poro-elastic swelling geomaterial, the governing equation for the conservation of momentum can be expressed by Equation (3), where G is the shear modulus (Pa) λ is the Lamé constant (Pa) β sw is a swelling coefficient (1/Pa), as discussed in Section 2.4 K is the bulk modulus of the solid skeleton (Pa) X is a parameter related to the degree of saturation of the soil [26] F v,i is the volumetric body force tensor (kg m −2 s −2 )

Intrinsic Permeability and the Klinkenberg Effect "Slip Flow"
In addition to the constitutive relations for the hydraulic behavior provided in Dagher et al. [11] and the companion paper to this one [25], this paper investigates the use of the Klinkenberg effect to represent a value of the intrinsic permeability of gas that is different from that of water. In this study, the intrinsic permeability tensor is assumed to be isotropic. Therefore, off-diagonal components are nil and diagonal components are equal to a value k, where δ ij is the Kronecker delta (identity tensor) (adimensional). For fine-grained expansive clays, which have very small grain-sizes less than 2 µm and are non-uniform, Pall and Moshenin [27] proposed the following equation based on a volume-surface mean diameter and changing porosity to account for non-uniformity of the soil particles. This equation was applied to the author's previous model [11].
where D vs is the volume-surface mean diameter (cm). The porosity is calculated by, n = n 0 + ε vol (6) where n 0 is the initial porosity of the porous media (m 3 voids · m −3 total). As the intrinsic permeability is not a property of the porefluid but rather a property of the geometry of the porous media (i.e., porosity, pore shape, and pore size distribution), it is reasonable to assume that the permeability of porewater and poregas in the same soil specimen would be the same. However, a number of studies have shown that the intrinsic permeability of gas is often higher than that of Minerals 2020, 10, 32 6 of 33 water and may be related to the poregas pressure [28][29][30]. The Klinkenberg effect, also described as "slip flow", has been used to conceptualize this phenomenon, whereby slip occurs between the gas molecules and the solid wall, resulting in an increase in the intrinsic permeability of gas as poregas pressure decreases [31].
Equations (7) and (8) provide the relationship between the intrinsic permeability value of water and the intrinsic permeability value of gas a function of pressure.
where k g is the intrinsic permeability value of gas (m 2 ). k w is the intrinsic permeability value of water (m 2 ) b ke is the Klinkenberg slip factor (Pa) c ke is a constant (adimensional) κ B is Boltzmann's constant (J K −1 ) T is the temperature (K) r is the pore radius (m).

Swelling Strain Component
Using the elastic framework, the total strain can be divided into an elastic and a swelling strain component as follows: ε kl = ε el ij + ε sw ij (9) where ε kl is the total strain tensor (adimensional) ε el ij is the elastic strain component (adimensional) ε sw ij is the swelling strain component (adimensional) For an expansive/swelling soil, the swelling strain increment, dε sw ij is a function of suction, s, or the degree of saturation, S w , whereby an increase in suction (or likewise a decrease in degree of saturation) leads to shrinkage of the geomaterial and an increase in porosity due to the increase in volumetric strain provided by the swelling strain [32]. Assuming isotropic swelling, this can be expressed by the following equation, where β sw is a swelling coefficient (1/Pa) s = p g − p w is the suction (Pa).
To simplify the relationship between swelling strain and suction, Nasir et al. [32] proposed a linear swelling strain with a constant swelling coefficient, β sw , for a bentonite-sand mixture with a relatively good fit to experimental swelling stress results.

Overview of the Numerical Model
For the process simulation study, a 3D HM coupled multiphysics numerical model was built using COMSOL ® to numerically solve the governing equations and the constitutive relations described in Section 2. Using the FEM, the model is used to simulate the simultaneous migration of gas and liquid (two-phase flow) in porous media, which are coupled to the mechanical behavior of the solid matrix.
For process testing and enhanced flow analysis study, the numerical model setup is based on laboratory experiments conducted by the BGS to simulate 1D gas flow through a saturated bentonite sample under constant volume boundary stress conditions [33]. As per the 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. Key phases of the experiment include the hydration phase (t = 7.3 days to 39 days) where the sample was left to saturate, the gas injection phase (t = 39 days to 71 days) where gas was continuously injected into the system, resulting in increasing gas pressures, and the gas shut-off (t = 71 days) where the injection pump was stopped, and gas pressures were allowed to dissipate naturally. The experiment was conducted under isothermal conditions at a temperature of 293.15 K. During the experiment, a number of parameters were 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. The BGS provided the experimental data that were used in the process testing and enhanced flow analysis study [33].

Model Geometry and Mesh
The model geometry and the finite element mesh are presented in Figure 1a,b, respectively. Although the experiment attempts to simulate 1D conditions, the actual conditions are truly 3D due to heterogeneity of the sample and the presence of a central injection rod inserted at the outlet face. This central injection rod is not used in this experiment but in subsequent ones to be simulated in another paper. Therefore, the FE mesh, as shown in Figure 1b, is 3D, consisting of approximately 11,000 elements.
Minerals 2019, 9, x FOR PEER REVIEW 7 of 32 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. The BGS provided the experimental data that were used in the process testing and enhanced flow analysis study [33].

Model Geometry and Mesh
The model geometry and the finite element mesh are presented in Figure 1a,b, respectively. Although the experiment attempts to simulate 1D conditions, the actual conditions are truly 3D due to heterogeneity of the sample and the presence of a central injection rod inserted at the outlet face. This central injection rod is not used in this experiment but in subsequent ones to be simulated in another paper. Therefore, the FE mesh, as shown in Figure 1b, is 3D, consisting of approximately 11,000 elements.

Material Properties
Material properties for solid bentonite MX-80 soil matrix, helium gas, and water are provided in Table 1 and were adopted from Dagher et al. [11]. In the proposed numerical model, the expression for the density of gas, ρ , and bulk modulus, K , provided in Table 1, were derived from the ideal gas law, as described in Dagher et al. [11]. As the behavior of gases at high pressures deviate from that of an ideal gas, this presents an over-simplification. In light of this, future work will consider applying a correction factor to these expressions.

Material
Parameter Name Symbol (Units) Value

Material Properties
Material properties for solid bentonite MX-80 soil matrix, helium gas, and water are provided in Table 1 and were adopted from Dagher et al. [11]. In the proposed numerical model, the expression for the density of gas, ρ g , and bulk modulus, K g , provided in Table 1, were derived from the ideal gas law, as described in Dagher et al. [11]. As the behavior of gases at high pressures deviate from that of an ideal gas, this presents an over-simplification. In light of this, future work will consider applying a correction factor to these expressions. (1) As a simplification to the model, the density of gas and bulk modulus are derived from the ideal gas law, as described in Dagher et al. [11].

Initial Value Condition
The initial conditions at t = 0 across the domain are provided in Table 2. The initial conditions assume an initial pore gas pressure equal to atmospheric within the bentonite sample and a 96% degree of water saturation. Table 2. 1D flow case: initial value conditions.

Boundary Conditions
The hydraulic and the mechanical boundary conditions (BC) for gas transport, water transport, and momentum transport are provided below.

Gas Conservation BCs
For gas transport, a no flow Neumann BC, ∂p g ∂x i = 0, was set at the radial boundaries. Dirichlet BCs were set at a specified gas injection pressure, P inj , for the front-face boundary and atmospheric pressure, P atm , at the back-face boundary.
For the concentration of gas in porewater, C g, H 2 O , the ideal gas law was assumed, where p g is the poregas pressure (Pa). Assuming instantaneous dissolution, the concentration of gas in the porewater can be calculated by multiplying Equation (13) by Henry's coefficient, H, and the portion of water in a unit volume, nS w where C g, H 2 O is the concentration of gas in water (kg gas m −3 water).

Water Conservation BCs
For water transport, Dirichlet BCs were set at a value equal to the water backpressure at the back-face boundary as well as at the radial porefluid arrays. This is a change from the BCs applied at the radial porefluid arrays in the previous study by the authors [11], which assumed no water flow through the radial porefluid arrays. This change in radial porefluid BC allowed for increased hydration during the hydration period, which ran from t = 7.3 days to 39 days and was more representative of the experiment [11,33]. A no flow Neumann BC, ∂p w ∂x i = 0, was set at the front face boundary and along the radial boundaries.

Momentum Conservation BCs
For the momentum conservation equation, a roller constraint was applied along the upper, the lower, and the radial boundaries to represent a condition whereby the boundary was free to move in the tangential direction but was fixed in the normal direction, simulating an overall constant volume condition, as per the experimental set-up.
The gas injection pressure and the water backpressure Dirichlet BCs were imported from the experimental data provided by BGS and were calculated from the previously described theory plotted in Figure 2a, while the concentration of dissolved gas BCs was plotted in Figure 2b.

Momentum Conservation BCs
For the momentum conservation equation, a roller constraint was applied along the upper, the lower, and the radial boundaries to represent a condition whereby the boundary was free to move in the tangential direction but was fixed in the normal direction, simulating an overall constant volume condition, as per the experimental set-up.
The gas injection pressure and the water backpressure Dirichlet BCs were imported from the experimental data provided by BGS and were calculated from the previously described theory plotted in Figure 2a, while the concentration of dissolved gas BCs was plotted in Figure 2b.

Modelling Approach
In this study, simulation of a number of study scenarios was performed with increasing model complexity in order to gain an understanding of the contribution of heterogeneity, the Klinkenberg "slip flow" effect, and the inclusion of a linear swelling strain on two-phase flow. The results of each simulation were compared to the experimental results and the effect of each on flow behavior analyzed [33]. Table 3 summarizes the simulations run and the values of key model parameters. It should be noted that these model parameters were selected based on calibrating the model to best match the experimental results while maintaining numerical stability and convergence.

Modelling Approach
In this study, simulation of a number of study scenarios was performed with increasing model complexity in order to gain an understanding of the contribution of heterogeneity, the Klinkenberg "slip flow" effect, and the inclusion of a linear swelling strain on two-phase flow. The results of each simulation were compared to the experimental results and the effect of each on flow behavior analyzed [33]. Table 3 summarizes the simulations run and the values of key model parameters. It should be noted that these model parameters were selected based on calibrating the model to best match the experimental results while maintaining numerical stability and convergence.
Additionally, a special study, S0, utilizes the simplified single-phase flow model that was used in the verification of the numerical model [25]. This scenario was included to assess the difference in increasing model complexity from a single-phase model to a two-phase flow model. This model is identical to that of the base case, S1, but is simplified to be consistent with the verification study through the following assumptions:  For this process simulation study, the simplified single-phase flow model applies the same material properties used in the verification study [25] but now applies boundary conditions and initial conditions that were described above. As per the assumptions, the material properties for the parameters used for scenario S0 remain constant through the model run and are presented in Table 4. The van Genuchten equations for the SWCC and the AEV are provided in the authors' original paper [11].

Introduction of Heterogeneity
Heterogeneity was introduced into the simulated sample by spatially applying a random normal distribution to the initial porosity material property with a mean porosity set at the experimentally determined value of 0.44 and applying a standard deviation of 0.06. There was no spatial correlation in heterogeneity intended with this approach. The application of this method to generate heterogeneity is not based on any experimental literature and was applied in this study in an attempt to trigger the generation of preferential flow pathways by the given models. A standard deviation of 0.06 was selected, as it provided the largest degree of heterogeneity while maintaining model stability. The initial porosity distribution within the bentonite sample is depicted in Figure 3, while Figure 4 provides the XZ-plane cross section of the initial porosity distribution. The initial porosity distr Based on the experimental intrinsic permeability of water in this study, k (3.4 × 10 −21 m 2 ), Figure 5 shows the relationship between poregas pressure and intrinsic permeability using a Klinkenberg slip factor based on Equations (7) and (8) with a constant of, c = 1.0 × 10 8 and c = 1.0 × 10 4 , as well as using Equation (15) to empirically solve for the slip factor. As depicted by the figure, the intrinsic permeability of the gas increases significantly with a decrease in poregas pressure, and as poregas pressures increase, the intrinsic permeability of the gas approaches that of water. As there is little information in literature regarding the slip factor for helium in bentonite clays, in order to assess whether the Klinkenberg effect could provide a significant role in matching key features of the experimental results, a constant, c = 1.0 × 10 8 , was used for the Klinkenberg effect. A value of c = 1.0 × 10 4 , which matches that of the empirical equation, would not provide the significant change in permeability needed to obtain complete breakthrough into the sample. The initial porosity distr Based on the experimental intrinsic permeability of water in this study, k (3.4 × 10 −21 m 2 ), Figure 5 shows the relationship between poregas pressure and intrinsic permeability using a Klinkenberg slip factor based on Equations (7) and (8) with a constant of, c = 1.0 × 10 8 and c = 1.0 × 10 4 , as well as using Equation (15) to empirically solve for the slip factor. As depicted by the figure, the intrinsic permeability of the gas increases significantly with a decrease in poregas pressure, and as poregas pressures increase, the intrinsic permeability of the gas approaches that of water. As there is little information in literature regarding the slip factor for helium in bentonite clays, in order to assess whether the Klinkenberg effect could provide a significant role in matching key features of the experimental results, a constant, c = 1.0 × 10 8 , was used for the Klinkenberg effect. A value of c = 1.0 × 10 4 , which matches that of the empirical equation, would not provide the significant change in permeability needed to obtain complete breakthrough into the sample.
Based on the experimental intrinsic permeability of water in this study, k w (3.4 × 10 −21 m 2 ), Figure 5 shows the relationship between poregas pressure and intrinsic permeability using a Klinkenberg slip factor based on Equations (7) and (8) with a constant of, c = 1.0 × 10 8 and c = 1.0 × 10 4 , as well as using Equation (15) to empirically solve for the slip factor. As depicted by the figure, the intrinsic permeability of the gas increases significantly with a decrease in poregas pressure, and as poregas pressures increase, the intrinsic permeability of the gas approaches that of water. As there is little information in literature regarding the slip factor for helium in bentonite clays, in order to assess whether the Klinkenberg effect could provide a significant role in matching key features of the experimental results, a constant, c = 1.0 × 10 8 , was used for the Klinkenberg effect. A value of c = 1.0 × 10 4 , which matches that of the empirical equation, would not provide the significant change in permeability needed to obtain complete breakthrough into the sample.

Coefficient of Swelling Strain
This study applies a swelling coefficient β (1/Pa) of 1 × 10 −9 . This coefficient was initially informed by that obtained by Nasir et al. [32] and calibration to best match the experimental data while maintaining numerical stability and convergence.

Results and Discussion
In this section, the compression-positive sign convention is used for stress components and pressures.
The results of the 1D flow process simulation and enhanced two-phase flow analysis study are presented below. For each study scenario, the following results are discussed: 1. poregas pressure evolution through the specimen 2. gas inflow/outflow 3. gas storage in the system 4. total stress evolution The experimental results provided for comparison were discussed in Daniels and Harrington [33] and in Dagher et al. [11]and are not discussed in detail in this paper. This paper focuses on the contribution of each enhanced flow mechanism to flow behavior.

Poregas Pressure Evolution through the Specimen
Poregas pressure profiles over time for the central cross section (XZ-plane) are depicted in Figure  6 for S0 and Figure 7a-g for S1 to S6. For the single-phase flow case depicted by Figure 6, the poregas pressure migration results follow a uniform gradient. This differs from the transient pore-pressure results provided in the companion paper [25], which demonstrate a non-linear kinked shape following a sudden increase in injection pressure until steady-state is reached. However, in the companion paper, at these same poregas pressures (10 MPa), steady-state is reached in approximately 1.5 h. In this study, our timescale is in days. Therefore, it is expected that transient response is very fast, and equilibrium is reached quickly.

Coefficient of Swelling Strain
This study applies a swelling coefficient β sw (1/Pa) of 1 × 10 −9 . This coefficient was initially informed by that obtained by Nasir et al. [32] and calibration to best match the experimental data while maintaining numerical stability and convergence.

Results and Discussion
In this section, the compression-positive sign convention is used for stress components and pressures.
The results of the 1D flow process simulation and enhanced two-phase flow analysis study are presented below. For each study scenario, the following results are discussed: 1.
poregas pressure evolution through the specimen 2.
gas storage in the system 4.
total stress evolution The experimental results provided for comparison were discussed in Daniels and Harrington [33] and in Dagher et al. [11] and are not discussed in detail in this paper. This paper focuses on the contribution of each enhanced flow mechanism to flow behavior.

Poregas Pressure Evolution through the Specimen
Poregas pressure profiles over time for the central cross section (XZ-plane) are depicted in Figure 6 for S0 and Figure 7a-g for S1 to S6. For the single-phase flow case depicted by Figure 6, the poregas pressure migration results follow a uniform gradient. This differs from the transient pore-pressure results provided in the companion paper [25], which demonstrate a non-linear kinked shape following a sudden increase in injection pressure until steady-state is reached. However, in the companion paper, at these same poregas pressures (10 MPa), steady-state is reached in approximately 1.5 h. In this study, our timescale is in days. Therefore, it is expected that transient response is very fast, and equilibrium is reached quickly.  For the base case scenario S1 representing the linear elastic advective-diffusive visco-capillary two-phase flow model represented in Figure 7a, the results match those first published in Dagher et al. [11], whereby once the AEV is exceeded (at 69 days), there is gas breakthrough into the sample followed by a slow moving uniform gas flow through the specimen. The advective poregas front only travels roughly 5% through the sample.
Results of S2 are provided by Figure 7b and show that, once heterogeneity is introduced, the poregas front travels slightly further than the base case and is no longer uniform with the formation of small gas fingers. Figure 7c shows the results of S3 once the Klinkenberg "slip flow" effect has been introduced. With a fairly large coefficient of the Klinkenberg effect of 1 × 10 8 , significantly increased gas migration into the sample and complete gas breakthrough are experienced. Up to 68 days, a relatively uniform gas front advances into the sample; however, as the front approaches the central rod, its shape becomes concave. As poregas moves forward, there is increased resistance to gas migration as a result of increased porewater pressure buildup around the central rod. This results in a lag in the poregas migration front, which is more pronounced closer to the center of the specimen. Once past the central rod, the poregas front returns to uniformity. It should be noted that, with the introduction of slip flow, there is some migration of poregas at low gas injection pressures; however, the predominant breakthrough of gas into the sample occurs once the AEV has been reached. One notable observation in the evolution of poregas migration is that the Klinkenberg effect tends to saturate the gas in the bentonite specimen and does not aid in the formation of distinct preferential flow pathways.
The results of S4, which introduce both the Klinkenberg effect and heterogeneity, are presented in Figure 7d. Characteristics of the results for both S2 and S3 can be seen, whereby the introduction of heterogeneity provides the formation of small fingers with the poregas now traveling quickly through the specimen. In fact, with the inclusion of heterogeneity, the poregas front migrates further relative to S3. One particular item to note is that, even with increased flow through the system relative to S2, the size of the fingers remains small and subtle, not large and distinct as the authors originally anticipated. Any fingers that do form eventually shrink and disappear. This may be due to a number of factors that could be physical, such as the presence of suction and diffusion, which aid in impeding the fingering capacity, or numerical and related to the method applied to introducing heterogeneity within the sample and to the mesh size.
In order to test whether Klinkenberg slip flow plays a role in mitigating the formation of gas fingers, the model is run using a much higher initial permeability of 5 × 10 −19 m 2 compared to that provided experimentally of 3.4 × 10 −21 m 2 . These results are provided for S5 in Figure 7e. The results do not show an increase in the number or the extent of gas fingers in comparison to S4, thereby demonstrating that the Klinkenberg effect may not play a large role in stabilizing the formation of gas fingers.
Finally, the effect of introducing a linear swelling strain is provided by S6 as depicted in Figure 7f. In comparison to S4, the presence of swelling does not result in a noticeable difference in the time it takes for the poregas migration front to travel through the sample. However, once gas shut-off occurs (at day 71), the presence of swelling in S6 appears to stabilize or reduce the effect of fingering in comparison to S4. Additionally, at day 75, an odd point of low gas pressure is exhibited. This is not a numerical figment but a result of fingering that occurs around a spot of low permeability and corresponds to that same point of low permeability presented in Figure 4.

Gas Inflow and Outflow
Total gas inflow and outflow profiles over time for the experimental data as well as for each study scenario are depicted in Figures 8-10. The experimental gas inflow shows a steep rise in gas into the sample when the injection pressure exceeds the AEV of 1 × 10 7 Pa (at 63 days). This is followed by chaotic inflow behavior and then a sudden drop in inflow, which corresponds to shut-off of the experimental injection pump (at 71 days). A similar behavior is exhibited for the experimental gas outflow, whereby gas flow out of the sample occurs almost immediately following gas flow into the sample, followed by a period of chaotic flow behavior and the occurrence of several sudden bursts

Gas Inflow and Outflow
Total gas inflow and outflow profiles over time for the experimental data as well as for each study scenario are depicted in Figures 8-10. The experimental gas inflow shows a steep rise in gas into the sample when the injection pressure exceeds the AEV of 1 × 10 7 Pa (at 63 days). This is followed by chaotic inflow behavior and then a sudden drop in inflow, which corresponds to shut-off of the experimental injection pump (at 71 days). A similar behavior is exhibited for the experimental gas outflow, whereby gas flow out of the sample occurs almost immediately following gas flow into the sample, followed by a period of chaotic flow behavior and the occurrence of several sudden bursts of outflow near the end of the experimental run. The authors interpreted this as the almost immediate formation of microfractures and the propagation of preferential flow pathways. Figure 8 depicts the modeled inflow and outflow curves for the single-phase flow scenario, S0. Figure 8a shows that the modeled gas inflow is not able to match the shape, the timing, or the magnitude of the experimental inflow curve. However, the inflow curve responds almost immediately to the injection pressure boundary condition with no rebound curve at each large change in injection pressure. This similar behavior is observed at the outflow with a magnitude several orders lower compared to the inflow, corresponding to the poregas pressure at outflow. These results are expected and are due to the very low constant gas permeability applied to the model (of the order of 10 −22 m 2 ). Figure 9a depicts the modeled inflow curves for scenarios S1, S2, and S3. The curve representing the base case is capable of capturing the general shape and timing of gas breakthrough into the bentonite specimen but is unable to capture the peak magnitude of the inflow or the sporadic inflow behavior during the period where gas is flowing into the sample. For S2, introducing heterogeneity results in slightly early breakthrough and slightly higher inflow compared to the base case. This is likely a result of the presence of areas within the injection face with larger pore sizes and consequently localized areas with lower AEVs, resulting in earlier breakthrough into the sample. For S3, when the Klinkenberg effect is introduced, there is a chaotic inflow behavior observed, which is similar to that of the experimental data. However, the model is unable to match the inflow magnitude and shows that inflow continues to occur even after the injection pump is stopped (i.e., day 71). The reason for this can be attributed to the Klinkenberg "slip flow" factors being applied. At approximately 63 days, once the AEV is reached and airflow into the sample begins, the gas permeability corresponding to a poregas pressure of 10 MPa is approximately 3.1 × 10 −17 m 2 (see Figure 5). Once gas injection is shut-off at day 71, the residual poregas pressure still remains above the AEV and increases exponentially as the poregas injection pressures decrease. Figure 9b depicts gas outflow curves. No outflow is observed for the base case or when applying heterogeneity to the sample. However, when the Klinkenberg effect is applied, a rapid increase in gas outflow is observed at approximately 88 days and continues with a steady declining outflow. This is due to the lag time it takes for the poregas migration front to completely flow through the sample (i.e., complete breakthrough). This modeled outflow behavior is not observed experimentally and is more representative of a plug flow behavior as opposed to dilation-controlled outflow as indicated by the experimental data.
as the poregas injection pressures decrease. Figure 9b depicts gas outflow curves. No outflow is observed for the base case or when applying heterogeneit   Figure 9. Results of (a) gas inflow and (b) gas outflow for S1, S2, and S3. Figure 10a depicts the modeled inflow curves for scenarios S4, S5, and S6. For S4, where both heterogeneity and the Klinkenberg effect are applied, the results show significantly higher gas inflow and a similar chaotic behavior when compared to the results of S3. A noticeable difference between the modeled results and the experimental results is that, following gas shut-off, there is an even larger increase in gas inflow between 78 and 95 days. As with S3, this is likely attributed to increasing gas permeability, as the residual injection gas pressures decrease while remaining above the AEV. The effect is even more pronounced as slip flow is now occurring along more localized pathways due to the presence of heterogeneity. When assessing the S4 outflow curve depicted in Figure 10b, outflow now occurs at around 69 days and matches the same shape and timing of the modeled inflow curve without the chaotic behavior. This may be due to flow stabilization occurring during migration through the sample as a result of suction and diffusion. Another unusual observation is that the Figure 9. Results of (a) gas inflow and (b) gas outflow for S1, S2, and S3. Figure 10a depicts the modeled inflow curves for scenarios S4, S5, and S6. For S4, where both heterogeneity and the Klinkenberg effect are applied, the results show significantly higher gas inflow and a similar chaotic behavior when compared to the results of S3. A noticeable difference between the modeled results and the experimental results is that, following gas shut-off, there is an even larger increase in gas inflow between 78 and 95 days. As with S3, this is likely attributed to increasing gas permeability, as the residual injection gas pressures decrease while remaining above the AEV. The effect is even more pronounced as slip flow is now occurring along more localized pathways due to the presence of heterogeneity. When assessing the S4 outflow curve depicted in Figure 10b, outflow now occurs at around 69 days and matches the same shape and timing of the modeled inflow curve without the chaotic behavior. This may be due to flow stabilization occurring during migration through the sample as a result of suction and diffusion. Another unusual observation is that the magnitude of outflow under the presence of heterogeneity is roughly 2.5 times higher than that of the inflow curve. This could only be attributed to a rapid increase in gas permeability throughout the sample that occurs as a result of a decrease in poregas pressures during shut-off, leading to rapid outflow and desaturation of the sample. None of the modeled simulations are able to reproduce the experimental inflow and outflow well. This suggests that the role of AEV and the presence of heterogeneity, slip flow, and reduction in porosity with increasing suction as result of a swelling strain cannot fully describe breakthrough through the system. Although the presence of heterogeneity and slip flow contribute to some preferential flow and chaotic inflow behavior, in order to obtain a sharp and rapid increase in flow For S5, when the Klinkenberg effect is restricted in the model and a much higher average initial intrinsic permeability of water is applied, the inflow profile more closely resembles that of the experimental data with some chaotic behavior, although not as sporadic as when Klinkenberg slip flow is present. However, following gas shut-off, the results of S5 do not result in a sharp decrease as observed experimentally, indicating that other mechanisms controlling the gas entry into the sample are occurring, such as damage or plasticity. For S5, complete gas breakthrough occurs at 90 days, but the outflow behavior is not representative of the experimental data.
For scenario S6, when a linear swelling strain is present, the inflow curves match the general shape and timing of those of S4 (without swelling) but reach a slightly higher peak inflow at 68 days when breakthrough into the sample first occurs. In comparison to S4, the peak inflow following gas shut-off achieved numerically at 104 days is reduced in the instance of linear swelling. This is expected by the model, as a decrease in suction at the injection interface due to the decrease in poregas injection pressure would result in a decrease in the volumetric strain and the corresponding porosity and gas permeability. This results in less flow into the sample. This effect is also observed in the outflow.
None of the modeled simulations are able to reproduce the experimental inflow and outflow well. This suggests that the role of AEV and the presence of heterogeneity, slip flow, and reduction in porosity with increasing suction as result of a swelling strain cannot fully describe breakthrough through the system. Although the presence of heterogeneity and slip flow contribute to some preferential flow and chaotic inflow behavior, in order to obtain a sharp and rapid increase in flow as observed experimentally, some other HM-coupled mechanism is needed to promote the rapid increase and decrease of flow into the system.

Gas Storage in the System
An important part of this process simulation study is a comparison of the modeled results to the experimental volume of gas stored in the system over time (i.e., inflow minus outflow). It is expected that dilation-controlled gas flow and the formation of preferential flow pathways would result in minimal gas storage within the system. This is observed by the experimental inflow and outflow results whereby the majority of gas entering the system almost immediately exits the system with several delayed bursts of outflow (at~80 days and~116 days) attributed to pathway sealing and the creation of new pathways. This is a key feature of the experimental inflow and outflow results originally noted by the experimenters [33]. This feature of the experimental results can be observed by a comparison of the experimental inflow and outflow (black lines) presented in Figure 9a,b, respectively. Figure 11 displays the experimental volume of gas storage in the system and the volume of gas stored for each model scenario. The experimental results show little gas storage within the soil specimen with approximately 0.002 m 3 of gas stored in the system at 120 days. The results of S0 show a large increase in volume of gas stored. Although flow into the sample is low due to the high-pressure gradient and the constant permeability, continuous gas migration occurs through the sample, resulting in a large volume of gas stored in comparison to the experimental data.
As for the two-phase flow models, the results of S1 and S2 match the experimental values quite well; however, this is likely only a coincidence, as in both cases (S1 and S2), gas flows into the sample with no advective gas outflow occurring.
For S3, the Klinkenberg effect results in a small volume of gas entering the system at lower gas pressures. Once the AEV is reached, the amount of gas stored in the system increases significantly. Once gas shut-off occurs, there is a desaturation of the gas as it flows out of the system. This behavior is in line with what is visually depicted in Figure 7c-h. When both heterogeneity and slip flow are present, there is an even earlier increase in gas stored within the system, followed by a longer period of gas desaturation. Introduction of a swelling strain in S6 and S7 does not significantly affect the volume of gas stored. The inclusion of the Klinkenberg effect provides a large deviation from the behavior observed experimentally, supporting the notion that inclusion of the Klinkenberg effect results in more plug flow-like behavior as opposed to the formation of more discrete preferential flow pathways. These results reinforce the conclusion that additional mechanisms of mechanical deformation are required to initiate dilation-controlled gas flow and the formation of discrete preferential flow pathways. present, there is an even earlier increase in gas stored within the system, followed by a longer period of gas desaturation. Introduction of a swelling strain in S6 and S7 does not significantly affect the volume of gas stored. The inclusion of the Klinkenberg effect provides a large deviation from the behavior observed experimentally, supporting the notion that inclusion of the Klinkenberg effect results in more plug flow-like behavior as opposed to the formation of more discrete preferential flow pathways. These results reinforce the conclusion that additional mechanisms of mechanical deformation are required to initiate dilation-controlled gas flow and the formation of discrete preferential flow pathways. Figure 11. Volume of gas stored for S1-S6.

Evolution of Total Stresses
The total stress evolution as measured by injection load cell, radial load cells, and backpressure load cell is presented in Figure 12. The experimental results depict an initial large increase in total stress during the hydration phase (t = 7.3 days to 39 days) followed by stabilization of the stresses. Figure 11. Volume of gas stored for S1-S6.

Evolution of Total Stresses
The total stress evolution as measured by injection load cell, radial load cells, and backpressure load cell is presented in Figure 12. The experimental results depict an initial large increase in total stress during the hydration phase (t = 7.3 days to 39 days) followed by stabilization of the stresses. Once the injection pressure exceeds the AEV (t = 63 days), there is another rapid increase in total stresses until gas injection ceases. In their previous work, the authors were not able to replicate the magnitude of the total stress evolution within the hydration phase and attributed this, in part, to neglecting the effect of swelling pressure in the mathematical model, a key behavior of expansive clays. Three key attributes are assessed within this study to attempt to replicate the total stress evolution observed experimentally: (i) modified BCs to allow for additional hydration of the sample from the radial porefluid arrays, (ii) Klinkenberg effect leading to increased gas flow through the system, and (iii) the application of a swelling strain to simulate swelling of the clay. Figure 12a provides the total stress evolution for the single-phase flow case, S0. The results show a minor stress response with increasing injection gas pressures when compared to the experimental data. This is to be expected, as the single-phase flow model only considers changes in porefluid pressure associated with an increase in poregas pressure at a constant χ = 0.9 and ignores the porefluid pressure associated with the immobile water phase. Figure 12b provides the total stress evolution for the base case, S1. The results show much better agreement of the total stress evolution during the hydration phase in comparison to the authors' previous work [11], as it is able to capture the shape, the timing, and the magnitude of the stress evolution. This is evidently a direct result of the change in BCs applied in this study. However, S1 is not capable of simulating the total stress response observed following breakthrough into the sample. Figure 12c provides the total stress evolution for S2. The results are similar to those of S1 but with a larger spread between the total stress curves due to the introduced heterogeneity. During the hydration phase, the total stress profiles at the axial load cells (i.e., injection load cell and backpressure load cell) do not increase to the same magnitude as the experimental data, although the total radial stresses are captured fairly well. This may be due to increased heterogeneity at the injection and backpressure fronts during installation of the specimen. Again, S2 is not capable of simulating the total stress response observed following breakthrough.
The results of introducing "slip flow" are provided in Figure 12d. For S3, when the Klinkenberg effect is introduced, the modeled results following breakthrough capture the experimental behavior very well, with a steep increase in total stress at both the axial and the radial load cells, followed by a subtle tail off once injection stops. This increase in total stress is due to achieving complete breakthrough through the sample. The model does have difficulties reaching the maximum total stresses recorded experimentally at the injection and the backpressure load cells; however, there is reasonably good agreement. Figure 12e shows the stress evolution for S4. When heterogeneity and Klinkenberg effect are both introduced, the model is able to achieve a higher separation of stress curves and match the range of peak total stresses recorded experimentally, although it is still not quite able to match the total axial stresses exhibited during the initial hydration period. Furthermore, a deeper analysis of the stress evolution at each individual load cell array shows an overprediction of total stress at radial load cell 1 (RLC1) and an underprediction at the backpressure load cell (BLC). However, this may be due to the specific random distribution of heterogeneity applied in our model. If multiple iterations of the model are run using different randomly generated porosity distributions, the total stress evolution behavior exhibited experimentally by the individual load cell arrays may be realized by the model. Figure 12f shows the stress evolution for S5, whereby the initial intrinsic permeability is set two orders of magnitude higher than that measured experimentally in order to observe the effect on flow behavior. For this scenario, the results show a much greater increase in total stresses from t = 0 days to t = 7.3 days, followed by fairly stable stress curves throughout the remainder of the simulation. Of particular note is that there is minimal stress response during breakthrough. These results, when compared to the experimental results, support the notion that the Klinkenberg effect may play an important role in the gradual stress response of the system.
The results of S6 are provided in Figure 12g. The introduction of a swelling strain provides a larger separation between the total stress curves. It also reduces the total stresses on the system in comparison to S4, with the exception of a peak total stress observed once gas breakthrough into the sample occurs, followed by a rapid decrease. This behavior in total stress under the presence of swelling is expected, as an increase in gas pressure results in an increase in suction and a decrease in volumetric strain. In our process simulation study, this is most evident at the peak at RLC1, whereby once the AEV is reached, there is a rapid increase in total stress in the system. Following gas shut-off, there is an even steeper decline. It should be noted that this behavior is not observed experimentally, although use of a non-linear swelling strain may lead to a better fit. sample occurs, followed by a rapid decrease. This behavior in total stress under the presence of swelling is expected, as an increase in gas pressure results in an increase in suction and a decrease in volumetric strain. In our process simulation study, this is most evident at the peak at RLC1, whereby once the AEV is reached, there is a rapid increase in total stress in the system. Following gas shut-off, there is an even steeper decline. It should be noted that this behavior is not observed experimentally, although use of a non-linear swelling strain may lead to a better fit. In light of these results, more work is needed to better fit the total stress evolution during the hydration period. This may be through re-evaluating the initial degree of saturation of the specimen, the initial hydration provided by the experiment, and additional calibration of the non-linear swelling stress parameters.

Mesh Dependency, Time Stepping, and Relative Tolerance
Based on the results of this study, it was identified that a number of numerical factors play a significant role in the numerical modelling of two-phase flow. These factors include the output time step selected in COMSOL ® , the mesh size, and the relative tolerance.
During the numerical analysis, an output time step of 0.1 day was set in COMSOL ® in order to capture the chaotic peak behavior observed by the experimental results. During the model runs, when a higher output time step of 1 day was applied, the model results were not able to capture the full oscillations well. Output time steps of 0.01 day were also performed, and although the results provided even greater resolution of the chaos observed, in order to balance computational time and the resolution needed to capture main features of the experimental data, it was determined that an output time step of 0.1 day was adequate. It should be noted that COMSOL ® applies an effective time step smaller than this but based on the user defined time step.
For each model run, a user defined optimum relative tolerance had to be set in order to reach convergence. The selected relative tolerance had to be set low enough to minimize the relative error that could be generated by the numerical model while simultaneously allowing for enough buffer to converge in areas with dynamic changes in poregas pressure.

Extent of Heterogeneity
It should be noted that the authors originally modeled small heterogeneities within the bentonite sample considering a mean porosity distribution of 0.44 with a standard deviation of 0.01. However, In light of these results, more work is needed to better fit the total stress evolution during the hydration period. This may be through re-evaluating the initial degree of saturation of the specimen, the initial hydration provided by the experiment, and additional calibration of the non-linear swelling stress parameters.

Mesh Dependency, Time Stepping, and Relative Tolerance
Based on the results of this study, it was identified that a number of numerical factors play a significant role in the numerical modelling of two-phase flow. These factors include the output time step selected in COMSOL ® , the mesh size, and the relative tolerance.
During the numerical analysis, an output time step of 0.1 day was set in COMSOL ® in order to capture the chaotic peak behavior observed by the experimental results. During the model runs, when a higher output time step of 1 day was applied, the model results were not able to capture the full oscillations well. Output time steps of 0.01 day were also performed, and although the results provided even greater resolution of the chaos observed, in order to balance computational time and the resolution needed to capture main features of the experimental data, it was determined that an output time step of 0.1 day was adequate. It should be noted that COMSOL ® applies an effective time step smaller than this but based on the user defined time step.
For each model run, a user defined optimum relative tolerance had to be set in order to reach convergence. The selected relative tolerance had to be set low enough to minimize the relative error that could be generated by the numerical model while simultaneously allowing for enough buffer to converge in areas with dynamic changes in poregas pressure.
With respect to the investigation into the formation of gas fingers, as discussed above, a number of physical conditions are required in order to simulate viscous fingering. These include differences in viscosity of both fluids, heterogeneity in the porosity and the permeability of the geomaterial, and a very fine mesh. Our modeling results were able to produce gas fingers. This was a result of differences in viscosities of each fluid (helium and water) considered within their own transport equations (Darcy's law), the inclusion of heterogeneity in pore size, and the selected mesh size. The results showed that the extent of gas fingering was more strongly related to changes in permeability within the soil specimen than size of the mesh. Additionally, when fingering was present, the size of the gas fingers was short and their formation was short lived, likely because the effects of gravity, diffusion, and suction did in fact hinder formation of long discrete fingers, as was identified by previous studies on gas fingering [23].
The results of this paper support the conclusion that other highly coupled HM mechanisms must be at play. In light of this, future studies will assess the contribution of several types of stress-strain relationships on flow behavior, including damage and poro-elastoplasticity.
The results of this study conclude that, in order to mimic dilation and dilatancy-controlled gas flow, additional considerations to the stress-strain behavior are required. Future studies will look at the effects of different stress-strain constitutive models on flow behavior. These will include a detailed assessment of the coupling between flow and stress state triggered by mechanical damage and plasticity. Future studies will also investigate the application of strain-localization and channeling in an attempt to simulate dilatancy-controlled gas flow and exasperate dilation and creation of preferential flow paths. Finally, multiple model iterations of randomly normally-distributed initial porosities will be assessed in an attempt to define, through a probabilistic assessment, the chaotic nature of dilation pathways.