Discrete Dilatant Pathway Modeling of Gas Migration Through Compacted Bentonite Clay

A coupled multiphase fluid flow and discrete fracturing model is applied to simulate bench-scale gas migration experiments on compacted bentonite. The numerical modeling is based on the linking of the multiphase fluid flow simulator TOUGH2 with a Rigid-Body-Spring Network model, which enables a discrete (lattice) represen- tation of elasticity and individual fractures. The evolution of a complex network of dilatant flow paths is modeled through opening and breakage of lattice interface bonds between porous-elastic matrix elements. To achieve a good match with the experimental results, including an abrupt gas breakthrough along with pressure and stress responses, it was necessary to calibrate model parameters for (1) air-entry pressure, (2) shear and tensile failure of lattice interface bonds, (3) moisture swelling/shrinkage effects on stress, and (4) aperture-dependent permeability of dilatant flow paths. Our best-fit conceptual model considers a pervasive network of discrete flow paths propagating from the gas injection point, whereas some of the experimental data indicate the potential for heterogeneous and unstable flow paths. approaches for the quantitative treatment of gas migration in clay-based repository systems is needed. 12 In this paper we present the development and testing of one


Introduction
Gas migration in clay-based buffer and host rock materials has been one of the primary subjects in the field of radioactive waste disposal. Gas generation in a geological repository, a result of anaerobic corrosion of metals and other processes, could pressurize low-permeability materials in the engineered barrier systems to develop new dilatant pathways for gas movement, which could degrade the barrier performance in the long term. Many international research programs have been conducted, involving both laboratory scale and in situ experiments, 1-3 and substantial insights have been gained on the phenomenology of gas transport processes in bentonite and claystone.
Four primary phenomenological models describing gas flow, shown in Fig. 1, can be defined as follows: (1) gas movement by diffusion and/ or solution within interstitial fluids along prevailing hydraulic gradients; (2) gas flow in the original porosity of the fabric, commonly referred to as visco-capillary (or two-phase) flow; (3) gas flow along localized dilatant pathways, which may or may not interact with the continuum stress field; and (4) gas fracturing of the rock similar to that performed during stimulation of hydrocarbon reservoirs. 4,5 These processes can occur simultaneously in different parts of a buffer, whereas local behavior would likely involve a transition of the processes with the gas pressure buildup, from diffusion dominated towards dilatancy controlled and potential fracturing. A discrete gas phase form once the rate of gas production exceeds the rate of gas diffusion within the pores of the barrier or host rock. 5 The pressure buildup is expected to take thousands of years to reach the conditions required for the formation of dilatant pathway gas flow. Thus, modeling will be necessary to predict gas migration over the long-term repository post-closure period.
Various modeling approaches have been proposed for the interpretation of the experimental results and for the analysis of gas release scenarios from geological repositories in the context of long-term performance. However, the predictive capabilities of gas transport models are limited, and basic mechanisms of gas transport in bentonite are not sufficiently understood to provide for robust conceptual and quantitative models. [6][7][8][9][10][11] It is clear that conventional porous media advection, diffusion and two-phase flow concepts are not sufficient. Indeed, recent studies have shown growing evidence of complex fluid path behavior that is not fully understood. 2,10,11 As such, development and testing of new and novel numerical approaches for the quantitative treatment of gas migration in clay-based repository systems is needed. 12 In this paper we present the development and testing of one approach for modeling the gas migration through compacted bentonite, adopting a discrete fracture modeling approach. The work was performed as part of DECOVALEX (Development of Coupled models and their VALidation against Experiments) project, an international research and model comparison collaboration for understanding and modeling of coupled thermo-hydro-mechanical-chemical processes in geological systems. The DECOVALEX-2019 phase was running from 2016 through 2019, and this study falls under Task A related to gas migration through low permeability clay barriers. 5,13 In DECOVALEX-2019, Task A, data from a series of flow tests performed on initially saturated samples of compacted bentonite were made available to participating modeling teams. As summarized in Tamayo-Mas et al., 5 a wide variety of conceptual approaches and numerical models were applied by the different modeling teams. This paper presents the development and results of one of these conceptual approaches and models, a discrete fracture model with the added capability of handling dilatant pathways, applied to model these experiments. The idea of applying the discrete fracture model for modeling dilatant flow path is in line with recent studies showing the growing evidence of such flow path behavior. 2,10,11,14 For example, Fig. 2 shows μ-CT from Gonzales-Blanco et al. 14 of MX-80 bentonite at different states, as compacted, after saturation, and after gas injection. As observed in the figure, the compacted state displayed a grain structure that is easily distinguishable. After saturation, the sample showed a more homogeneous structure although some boundaries between the original grains are still present (the original grain structure and weak grain boundaries did not completely vanish during the saturation process). It appeared that the gas injection process reactivated these low-density boundaries (clay gel and pores) by inducing some opening and desaturation of the gas pathways. The model approach investigated here, reactivation of weak boundaries can be simulated, with opening and associated gas entry.
The simulator applied in this study for modeling dilatant flow gas is the TOUGH-RBSN simulator that combines the TOUGH2 multiphase fluid flow simulator 16 with the rigid-body-spring network (RBSN) model. 17,18 RBSN enables a discrete (lattice) representation of elasticity, individual fractures and fracture networks that are mapped onto an unstructured, 3-D Voronoi tessellation of a spatially random set of points. 19 The TOUGH-RBSN has been applied in the past for modeling drying induced cracks in clay, 17 which exhibit similar complex pattern as shown in Fig. 2. Previous studies with the TOUGH-RBSN have demonstrated the potential of modeling discrete fracture propagation through heterogeneous geological media in 3D settings. 17,18,20 The TOUGH-RBSN simulator could be categorized within the broad group of discrete fracture models for simulating coupled hydro-mechanical behavior in porous fractured geological media at multiple scales. [21][22][23] It is one of several TOUGH-based geomechanics simulators in which the TOUGH2 multi-phase flow simulator is sequentially coupled to a geomechanics simulator. 24 The TOUGH-RBSN simulator was selected because of its capability for modeling development of complex and intersecting fracture-like flow paths. 20 In Section 2 we briefly introduce the gas injection experimental data that were part of the DECOVALEX-2019, Task A. This is followed in Section 3 by a more detailed presentation of the TOUGH-RBSN simulator and its adaptation for modeling flow along dilatant pathways these experiments. Thereafter, in Sections 4 and 5, the modeling of the actual experiments is presented, starting with a linear flow experiment in Section 4 and a spherical flow experiment in Section 5. We end with a discussion and conclusions on the modeling results with respect to discrete fracture modeling of these types of processes.

Gas injection experimental data
Task A of DECOVALEX-2019 focuses on modeling two different gas injection experiments undertaken by the British Geological Survey (BGS): (a) a one-dimensional (1D) and (b) a three-dimensional (3D) test. The two tests were conducted on initially water-saturated MX-80 bentonite cylindrical samples (120 × 60 mm) emplaced in the same test cell, but in different configurations to study 1D or 3D gas flow (Fig. 3). The test cell is a constant-volume vessel meaning that the sample cannot expand and that stress will not be constant during the injection. The pressure vessel was instrumented with (i) two axial and three radial load cells, (ii) three radial arrays which allowed the continuous monitoring of pore pressure or gas flow within each array, and (iii) a central filter mounted at the end of a 6.4 mm diameter steel tube inserted into the sample, where gas injection pressure and flow can be measured in the case of 3D spherical gas flow (Fig. 3b). The experimental setup and an overview of the measurements are given in Daniels and Harrington 15 and Harrington et al. 10 whereas the experimental data for the two experiments used in Task A of DECOVALEX-2019 are described in detail in Tamayo-Mas et al. 5 While hydrogen will be the primary gas generated in a high level waste repository, helium was selected as a safe substitute based on its inert nature and similar molecular diameter. 10 Fig . 4 presents data from the 1D gas flow experiment. The test includes an initial hydration phase (0-39 days) followed by the actual gas injection test. Hydration was carried out to saturate the sample and develop swelling stress to represent the in-situ conditions for a buffer around a waste canister in radioactive waste repository. The sample was compacted to a dry density of 1.56 kg/m 3 and the hydration resulted in a radial total stress of about 7.5 MPa while the axial stress was higher, at about 9.5 MPa (Fig. 4c). Under these initial conditions, helium gas injection commenced at 39 days to simulate what could happen as a result of gas production in situ with gas overpressure increasing until gas breakthrough into the saturated bentonite. In detail, at Day 39 helium was added to the injection system to increase gas pressure to 3 MPa. 5 This was then held constant for a period of 7 days. At Day 46, the injection pump was set to a constant gas compression rate of 500 μl/h and the injection pressure gradually increased for the next 8 days from 3 MPa to 5 MPa whilst the volume of gas in the injection system decreased from 235 ml to 139.7 ml. At this point (Day 54), the gas compression rate was reduced to 375 μl/h. At Day 61, 59.95 ml further helium was added to the interface vessel, whilst maintaining constant gas pressure in the system. 5 At about 64 days, gas started to enter the bentonite sample and made a sudden breakthrough. The injection rate determines how fast pressure can increase before breakthrough. In these experiments the pressure build up takes about 1 month, whereas at a real repository setting such pressure build up would be expected to take thousands of years.
The experimental data in Fig. 4 show a typical behavior observed in these types of gas flow experiments on low permeability clay material. The gas injection experiment started with injection of gas into the injection system at a small rate leading to a slow increase in gas pressure in the injection system. At Day 64, when the injection pressure reached about 8 MPa, a sudden gas breakthrough occurred and gas entered the sample from one end as shown in Fig. 3a. The pore pressure and stress measured on the outside of the sample increased almost at the same time as the gas breakthrough into the sample occurred. Following the initial gas breakthrough, the total stress closely followed the pore pressure within the sample. This behavior continued following the cessation of pumping at 71 days, when gas pressure, total stress and pore pressure began to decay. During this later phase of testing (71-121 days), gas outflow was observed sporadically suggesting new gas pathways continued to open and close. 5 The spherical gas flow test consisted of a 3D gas injection test performed on another pre-compacted MX-80 bentonite sample. Here, the gas injection was conducted at a point inside the sample from the central filter leading to a spherical gas flow (Fig. 3b). Gas outflow was measured at the three radial arrays located at different axial distances from the injection point. Fig. 5 shows how several outburst of gas flow occurs at two of the arrays and later outflow focuses to one of the arrays until the end of the experiment. Radial and axial stresses show a complex response during the episodic flow, but later becomes uniform (except for one axial stress), following the pressure decaying injection pressure. 5,10

Modeling approach and development
This section provides a more detailed description of TOUGH-RBSN and how it is adapted to enable simulation of dilatant gas flow in bentonite based on a conceptual model and key features described in Section 3.1. A description of how these key feature are implemented into the numerical modeling along with a more detailed description of TOUGH-RBSN are given in Sections 3.2 through 3.5.
Initial simulations of the experimental data on gas migration in bentonite were performed using both continuum-and discontinuumbased coupled multiphase and geomechanics models. 25 Here we envision that gas flow taking place through the macro-pores or interfaces between individual clay aggregates or between clusters of clay aggregates separated by flow channels. We found that modeling of these experiments requires the following key features and properties: 1) A multiphase flow model containing an air-entry pressure, because the experimental results show that no gas enters the sample until gas pressure has increased to about 8 MPa, when a sudden gas entry into the sample occurs. 2) A poro-elastic model to simulate the stress response to pore pressure changes, while considering matrix moisture swelling/shrinkage or some other way to reduce the poro-elastic stress, or otherwise total stress will be overestimated. 3) A criterion for the creation and opening of discrete flow channels considering reactivation and damage of interfaces due both tensile and shear failure to produce observed abrupt gas breakthrough. 4) A conceptual model for channel-like gas migration in pathways with permeability described as a function of their aperture, to simulate high-peak gas flow rate and subsequent sealing.
Using the TOUGH-RBSN, this can be realized mathematically by dividing the model domain into (1) fractures or flow channels with their unique channel-like properties and (2) porous-elastic solids with claylike properties including swelling/shrinkage. The next two sections describe how these features are implemented in TOUGH-RBSN.
The computational domain for both the TOUGH2 and RBSN calculations is tessellated using a Voronoi diagram. 26 Fractures that in this case represents dilatant flow channels are explicitly modeled within the Voronoi grid. Voronoi cells generally represent the clay matrix component and pre-existing or newly generated fractures or opened flow channels are placed on the Voronoi cell boundaries (Fig. 6).  Delaunay edge defines the nodal connection of the corresponding lattice element. Through the dual Voronoi tessellation, the spatial domain is collectively filled with discrete polyhedral cells that render the elemental volumes. More detailed procedure of the domain partitioning is presented elsewhere. 27,28 In the simulations using a grid structure of the ordinary Voronoi discretization, flow and mass transfer are enacted only through the connections of the neighboring Voronoi nodes (called cell-cell connections in Fig. 6). However, if a dilatant flow path is opened through activation of an interface, substantially enhanced flow may arise through the dilatant flow path. As shown in Fig. 6, an interface node with increased permeability is inserted where the cell-cell connection intersects the Voronoi cell boundary. The original cell-cell connection is divided into two cell-interface (and vice versa) connections by the interface node. In addition, the connections between the interface nodes are established to activate flow channels in discrete fractures.
Elasticity and fracture-damage of geomaterials are modeled using the RBSN approach. As a discrete modeling approach, the RBSN represents the mechanical system by a collection of simple lattice (two-node) elements. The lattice topology is defined by the Delaunay tessellation, and the dual Voronoi diagram is used to render the volume of a discretized domain (Fig. 7a). The elemental formulations are based on the rigid-body-spring concept. 19,29 A lattice element consists of a zero-size spring set located at the centroid of the common Voronoi cell boundary and two rigid arm constraints that relate the spring set to the nodes (Fig. 7b). For 3D modeling, a spring set is formed from three axial springs and three rotational springs in local n − s − t coordinates (Fig. 7c). The n− axis is normal, and the st plane is parallel to the Voronoi cell boundary.
consists of the six spring coefficients which are defined according to the geometrical properties of Voronoi diagram: where E is the elastic modulus, J p , I ss , and I tt are the polar and two principal moments of inertia of the Voronoi cell boundary with respect to the centroid, respectively. The spring coefficients are scaled by the   Poisson ratio can be represented by adjusting α 1 and α 2 . By setting α 1 = α 2 = 1, the models behave with elastic homogeneity under uniform straining, albeit with zero effective Poisson ratio. 28,30 The stress tensors is evaluated at Voronoi cell nodes by considering the equilibrium conditions of the spring forces. Sets of the spring forces are applied at the boundaries surrounding a Voronoi cell (Fig. 8a), and nodal force components Fnn, Fns, and Fnt can be calculated for an arbitrary section passing through the Voronoi cell node with its corresponding local n − s − t coordinates, which satisfy the equilibrium condition with all the forces acting on the remaining cell boundaries (Fig. 8b). Moment contributions to equilibrium are not considered here. By dividing these force components by the cut-face area, the corresponding stress components σ n , σ s , and σ t can be obtained. By repeating this process for three mutually perpendicular sections, the full stress tensor is obtained (Fig. 8c). Details are given elsewhere. 27 From the stress tensors at two neighboring nodes, the stress tensor of the inter-element is calculated according to where σ i and σ j are the stress tensors at the neighboring nodes i and j, respectively.
In RBSN, the fracturing process is represented by the damage/ breakage of the springs. For a damaged spring set, the local spring coefficients are degraded as where ω is a scalar damage index with a range from 0 (undamaged) to 1 (completely damaged). In the modeling of brittle fracturing, which is  applied to the cases presented in this paper, ω is either 0 or 1. Fracture may initiate within a lattice element when the stress state exceeds the given material strength. To determine the criticality of the stress state, a stress ratio is calculated for each lattice element: where σ e is the element stress state and σ is the critical stress defined by fracture criteria. The maximum principal tensile stress of σ evaluated by Equation (2) serves as σe in Equation (4).
During iterative calculations, only one element, with the most critical stress state (i.e., the largest R f ≥ 1), is allowed to break per iteration, and the fracture event entails a reduction of spring stiffnesses and a release of the associated elemental forces. Strength properties of each element are determined as critical stress components on a Mohr-Coulomb surface. Fig. 9 shows the fracture surface with a tension cutoff, where the envelop is defined by three parameters: the angle of internal friction ψ (surface inclination with respect to σ n -axis); cohesive strength c (surface intersection with the shear axes); and the tensile strength f n . Thus, in this study the interfaces between Vorinoi cells can be activated by both shear and tensile failure.
Model simulation was carried out using a fluid equation of state module in TOUGH2 that handles liquid and gas phase with components of water and a gas obeying ideal gas properties of the helium gas injected in the experiment. Important properties for gas migration are the water retention curve and the relative gas permeability, which are different for the matrix and fractures. Here, we adopt a van Genuchten capillary pressure model 31 and a Corey relative permeability model, 32 which are both functions of the degree of saturation. Pressure and flow responses of a gas injected into an initially saturated bentonite sample may be only realized at a certain level of injection pressure or above. This conditional gas penetration is implemented by introducing an air entry pressure P ae , corresponding to a residual gas saturation (1 − S lr ) in the capillary pressure function. The van Genuchten capillary pressure model is used to define the water retention curve as with S * = (S − S lr )/(1 − S lr ), P' 0 and λ are capillary scaling and shape parameters For fractured elements, the apparent air entry pressure P' 0 of the element is scaled by the function of permeability as The relative permeability-saturation relationships of liquid and gaseous phases are parameterized using Corey model as where Ŝ = (S − S lr )/(1 − S lr − S gr ). and m g is a multiplying factor for the enhanced gas permeability. The residual saturations S lr and S gr are provided to limit the mobility of the respective phase, i.e., both liquid and gaseous phases can vary their mobilities only in the range of S = [S lr , 1 − S gr ]. To avoid unphysical situation with P c = ∞, larger S lr for the relative permeability is usually chosen as compared to S lr for the capillary pressure. 32 S gr can be used to define the air entry pressure in relationship with the capillary pressure function. 33 In two-phase cases, the pore pressure can be regarded as the sum of partial contributions of liquid and gas. Using Bishop's effective stress relationship, an equivalent pore pressure P φ is calculated as a function of gas pressure, P g , and liquid pressure, P l , according to: where χ is the Bishop's coefficient dependent on the degree of saturation. Khalili and Khabbaz proposed a general formulation for χ based on the relationship of capillary pressure and air entry pressure 34 : for P g − P l < P ae ( P g − P l P ae ) − 0.55 for P g − P l ≥ P ae (9) The equivalent pore pressure P φ is applied when calculating the effective stress as part of the two-way coupled procedure described in the next section. Fig. 10 shows a schematic flow diagram of the coupling procedure between TOUGH2 and RBSN codes. Coupling modules are implemented in each side of the modeling codes, by which material properties and mechanical boundary conditions are updated with the outputs of primary variables of physical quantities. Pore pressure and degree of saturation are primary variables of the TOUGH2 analysis, which are involved in the RBSN simulation to determine the effective stress and swelling/shrinkage strain. In return, the primary variables of the RBSN model, stress/strain states, damage index, and fracture aperture, are used to evaluate hydrological properties and conditions for the TOUGH2 simulation.
Based on the linear poro-elasticity theory, the effective (grain-tograin) stress σ n ′ is calculated from the pore pressureP φ35 : σ n ′ = σ n + α p P φ where σ n is the total normal stress obtained from overall loading, including external loads, and α p is the Biot's effective stress parameter.
The Biot's effective stress parameter can vary between 0 and 1, and is normally close to 1 for soft clay. 35 Note that the tensile stress is taken to be positive for the sign convention. Also, shrinkage/swelling effects due to the local changes of liquid saturations S can be considered using a simple linear swelling model as 36 : where ε s is shrinkage/swelling strain; and α s is the hydraulic shrinkage coefficient, which would be an empirical parameter to calibrate against laboratory test data. If a poro-elastic geomaterial is under mechanically confined conditions, the change in effective stress due to swelling/ shrinkage can be calculated as where E is the Young's modulus. Next, the RBSN to TOUGH2 link supplies the effective stress and the strain calculated in the lattice element to update the hydrological properties of the corresponding TOUGH2 grid blocks on the left side of Fig. 10. Porosity, permeability, and capillary pressure are generally related to the effective stress and strain values. 37 The flow model configuration comprises two types of elements. The first type of element is associated with individual Voronoi cells, representing the matrix/grain bulk in the computational domain. The permeabilities of these elements depend on porosity which can be described based on Kozeny-Carman equation applied to bentonite as 38 : where k 0 is the intrinsic permeability, and φ 0 is the initial porosity. The other type of element represents potential fractures or preexisting discontinuities or planes of weakness embedded in a portion of the matrix volume. It is assumed that these elements are positioned at the common boundary I between two adjacent Voronoi cells i and j (see Fig. 10), and their permeabilities can be calculated as the sum of two components: If an element is yet to be fractured (i.e., ω = 0), k matrix is calculated as in Eq. (6) and k fracture is simply assumed to be zero. On the other hand, once the interface element is fractured (i.e., ω > 0), k fracture will be predominant to enhance the total permeability and k matrix is assumed to revert to the initial intrinsic permeability k 0 . This conditional expression of the permeability is written as where a is the TOUGH2 element width, and b h is the hydraulic conducting fracture aperture. The hydraulic aperture is coupled to the mechanical aperture b m 39 : where x = 1/2(x +|x|) represents Macaulay brackets, and b r is the residual hydraulic aperture. f ≤ 1.0 is a dimensionless factor that accounts for the slowdown of flow in a natural rough fracture in comparison to the ideal case of parallel smooth fracture surfaces. 39 In the TOUGH-RBSN sequential coupling, TOUGH2 is a main driver that controls the time stepping during the coupling procedure, while RBSN solves the mechanical response as a quasi-static process at each time step. The selection of small time steps is important to find stable solutions of the mechanical response, so the TOUGH2 input defines the upper limit for time step size with a small value to avoid any abrupt change of hydrological conditions over time steps. Time steps tends to be reduced to very small values when interfaces are activated either in shear or tensile mode. Such shortening of time steps is also a result of having abrupt changes in saturations when gas invades and newly created gas flow path. Similar to other TOUGH-based geomechanical simulators, the sequential coupling scheme is based on the so-called stress-fixed iteration scheme, with the sequential solution becoming unconditionally stable. 40

Results -modeling the 1D gas flow experiment
The 1D gas flow experiment was simulated using a two-dimensional 120 × 60 mm rectangular domain discretized with an unstructured Voronoi grid with 1401 cells and 3840 lattice elements (Fig. 11). An additional 2 mm-thick layer of elements was placed as padding on each side of the domain boundaries to provide a uniform mechanical confinement (fixed displacement boundary conditions). Colored marks in Fig. 11 indicate the measuring points of pressure evolutions during simulations, which coincide with the locations of porewater pressure sensors in the experiments. 5,15 The initial confining stress values are given from measurements with load cells, which are 9 MPa in the axial direction and 6 MPa in the radial direction. The initial pore pressure is given as 1 MPa, and the initial saturation of the model is set to 0.999 for fully saturated conditions. While the padding elements surrounding the specimen are mechanically fixed, gas (air) injection is applied at the left boundary, for which the injection pressure is controlled to match the experimental values. A backpressure of 1 MPa at the right boundary is maintained during simulations. Hydrologically, the left and right padding elements are A large number of calibration simulations were performed in trial and error mode to investigate the key parameters that needed to be varied to be able to match the experimental results. The final set of parameters are listed in Table 1 with markings on which of the parameters that were fixed and which were calibrated Some supposedly wellknown material properties for MX-80 bentonite were fixed within the DECOVALEX-2019, Task A. This included a saturated water permeability, k = 3.4 × 10 -21 m 2 , porosity, φ = 44%, Young's modulus, E = 307 MPa, and Poisson's ratio, ν = 0.4. 5 In addition, the parameters for the capillary pressure and relative permeability functions (Equations (5) and (7)) were adopted from Senger and Marschall, 33 who used data from MX-80 bentonite of similar dry density. With these parameters fixed, some of the other parameters were calibrated to capture some of the key features in the experiment as follows.
1) With permeability fixed at 3.4 × 10 -21 m 2 , gas had to be prevented from entering the sample during the first 65 days by applying a sufficiently high air-entry pressure (8 MPa). Here we calibrated the residual gas saturation, S gr = 0.0877, which corresponds to the air entry pressure of 8 MPa at a minimal desaturation. This implies that the injection overpressure has to be 8 MPa before gas enters the sample.   Fig. 12 shows the comparison of simulation results with the experimental data for pore pressure, outflow rate, as well as for axial and radial stresses. In the pore pressure and stress evolutions, thick solid lines denote simulation results and thin dotted lines represent experimental data. The simulation captures very well the abrupt changes that occur after 65 days as a result of the gas breakthrough. Noted discrepancies between the simulation and the experiment are the longer tail of gas flow in the simulation (red curve in Fig. 12b) and the somewhat overestimated total stress magnitude ( Fig. 12c and d). Nevertheless, there is an overall good agreement regarding these macroscopic hydromechanical responses that are measured on the outside of the sample. Fig. 13 presents simulation results of the dilatant flow paths (fracturing) and pressure contours during the gas breakthrough, at around 61-66 days. In the early time (61.33 days), some preferential fracturing forms distinct flow paths, but later, along with the higher-pressure advection into the sample, a uniform fracture pattern can be observed across the sample. The dominant failure mechanism to form the fracturing pattern is shear failure. The shear strength of these bonds was calibrated to a low value (using c = 0.1 MPa and β = 18 • ), which may reflect the granular texture of bentonite where the grain-to-grain strength is expected to be very low.
Fig. 14 presents one relevant example result from a large number of sensitivity studies that were performed. In this case, a high cohesive strength was applied so that tensile failure became a dominant failure mechanism. Results in Fig. 14 are reasonable with the formation of a preferential path of gas breakthrough. However, the pressure does not change at the boundary of the model and hence there would be no response at the pressure sensors. Thus, there is no possibility to match the observed pressure response on the outside of the sample when model simulations display a preferential flow path confined to the center of the sample.

Results -modeling the 3D gas flow experiment
For modeling of the 3D (spherical) gas flow experiment, we generated a 3D Voronoi mesh of a cylinder with dimensions of 120 mm height and 60 mm diameter, which is composed of 7856 Voronoi cells and 33,316 element connections (Fig. 15). Additional padding elements were placed on the axial ends and the circumferential surface to provide a constant volume boundary condition (zero-displacement constraints).
Initial mechanical confinement is assumed with the axial stress of 7.25 MPa and the radial stress of 7.75 MP in the sample. The initial pore pressure is set to 1 MPa representing a fully water-saturated state. While the padding elements surrounding the sample domain are hydrologically tight, a fully gas-saturated element with 168 ml of volume represents the injection system, which is directly connected to the center of the sample domain.
A constant injection rate of 2.75 × 10 -9 m 3 /s at standard temperature and pressure (STP) conditions is maintained throughout the simulation with fluid storage of the injection system calibrated to match the observed early time injection pressure evolution. As indicated by green marks in Fig. 15, 12 locations of radial filters on the circumferential surface of the sample are taken as outflow boundaries, where the corresponding elements are connected to the backpressure boundary element with a constant 1 MPa pressure. Outflow values are measured at the connections between those circumferential elements and the backpressure boundary element. Reaction forces are measured at the constrained padding elements to derive stress values, which are colored in red in Fig. 15b.
The 3D gas flow experiment was simulated starting with the set of parameters used in the previous 1D flow simulation, though some parameters had to be adjusted in order to match the 3D experimental data. This included a slight reduction in air entry pressure from 8 to 7 MPa. Moreover, alternative coupled hydromechanical models were applied considering either the effective stress in Equation (10) together with the linear swelling stress model in Equation (12), or the Bishop's effective stress model with equivalent pore pressure from Equation (8). This difference in the poro-elastic model as well as difference in model geometry (2D vs 3D) also resulted in a recalibration of strength parameters with lower cohesion and frictional angle than those used in the previous 1D flow simulation. Fig. 16 shows snapshots of fracture development during gas injection. Fractures are initiated at the center of the sample, where the gas injection occurs. Then, the main fracture cluster grows into a spherical shape toward the lateral surface of the sample. After the fracture cluster reaches the lateral surface, near-surface fractures propagate in the longitudinal direction. These non-uniform fracture patterns form a heterogeneous enhanced permeability domain within the sample, where preferential flow paths are generated. Fig. 17 presents the comparison of simulated and experimental results. Overall, the model reasonably well captures the pressure and flow response, while the stress responses tend to be overestimated. The pressure and flow responses are in good agreement regarding the timing of gas breakthrough and peak rate, but the model result tends to be smoother without some of the strong oscillations displayed in the experimental data ( Fig. 17a and b). Simulated stress is in good agreement until gas breakthrough, whereas after gas breakthrough, stresses are overestimated at all locations except for the radial stress at R1 (Fig. 17c and d). There is, however, another significant difference regarding flow behavior. In the experiment, only one radial array takes a dominant outflow, whereas the simulation exhibits more uniform outflow rates for all three radial arrays. This shows that the flow behavior in the experiment was more heterogenous in terms of channelized gas flow paths than the modeling results display.

Discussion
The model results presented in this article correspond to one of the several models that were applied to simulate the same set of experiments within the DECOVALEX-2019, Task A. 5 The model presented here is the only discrete fracture model and it enables more explicit modeling of the underlying physics associated with the evolution of a network of discrete gas fracturing and dilatant flow paths. Table 1 presents all the parameters, including those fixed and those calibrated. For the calibrated parameters, the simulations showed that shear failure was dominant although tensile failure could also locally occur. Shear failure became dominant in the modeling because the sample is confined and the total stress increases along with fluid pressure. Under such conditions, it is difficult to achieve the tensile effective stress required for inducing tensile failure. The findings of dominant shear failure may be impacted by the fact that homogenous strength properties were assumed, whereas heterogeneous properties might provide some room for local tensile failure.
The total stress was challenging to model as the total stress increases significantly with the uniformly increasing pore pressure in the confined sample. One simple way to better match the observed stress magnitude would be to assume a very low value of Biot's parameter, e.g. α p = 0.1-0.5, as was done by many modeling teams in the DECOVALEX-2019 model comparison. 5 Such a low α p might be motivated if gas flow takes place through narrow discrete flow paths or channels. Using TOUGH-RBSN, we instead tried to match the stress evolution by the However, the standard Bishop's model applied still led to an overestimate of the total stress after the gas breakthrough and further studies would be needed to resolve this issue.
We note that the strength parameters for the fracturing process were calibrated to different values for the 1D (linear) flow model and 3D (spherical) flow model simulations. The difference in the values for the 1D vs 3D flow models may not be surprising as changes were made going from the linear swelling model to the Bishop's equation and going from a 2D plane strain analysis to a full 3D analysis. More consistent values will require the use of consistent poro-elastic model, e.g. Bishop's approach in both cases, as well as to run the two experiment with the same 3D cylindrical model geometry and a consistent mesh discretization. Thus further sensitivity modeling and calibration would be required before achieving consistent input parameters for both the 1D and 3D modeling cases. It would require to run full 3D models for the two types of experiments. Here we are facing limitation by the efficiency of the numerical simulations with TOUGH-RBSN. As mentioned, the time steps becomes very small when new discrete flow paths are opening with invasion of the gas phase. The simulations for the 3D cases took days on a standard desktop computer meaning that model calibration is time consuming.
Finally, we like to point out that uniform strength values were applied in the model and this resulted in a spatially uniform fracturing from the injection point, whereas data from the 3D experiment shows a more heterogeneous flow behavior. Application of more heterogeneous bond strength properties may give rise to more localized flow paths.
However, to fully evaluate the accuracy of such a model simulation, it would be necessary to have images or information about the actual flow paths inside the sample. Development of the recent work by Harrington et al. 10,11 using stress perturbation analysis to identify flow events may provide such data in the future.

Conclusions
We have applied a coupled multiphase fluid flow and discrete fracturing model (TOUGH-RBSN model) to simulate gas migration experiments on compacted bentonite. We modeled the evolution of a complex network of dilatant flow paths through opening and breakage of lattice bonds between porous-elastic matrix elements. For matching experimental results, including an abrupt gas breakthrough along with pressure and stress responses, it was necessary to calibrate model parameters for (1) an air-entry pressure, (2) shear and tensile failure of lattice bonds, (3) moisture swelling/shrinkage effects on stress, and (4) aperturedependent permeability of dilatant flow paths. For a good match to the 1D linear flow experimental results, the simulations indicate a pervasive network of discrete flow paths propagating from the gas injection point, whereas some data of the 3D experiment indicate more heterogenous and unstable flow paths. Moreover, model simulations showed that shear failure was dominant although tensile failure could also locally occur. Shear failure became dominant in the modeling because the sample is confined and the total stress increases along with fluid pressure. Under such conditions, it is difficult to achieve the tensile effective stress required for inducing tensile failure. The model tool developed and applied in this study is promising for modeling such underlying processes and the application to this problem would benefit improved imaging of such processes.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.