Implementation and validation of pressure-dependent gas permeability model for bentonite in FEM code Thebes

In an Engineered Barrier System of a nuclear waste repository, gas migrates through: a) diffusion/advection of dissolved gases, b) two-phase continuum flow, c) dilatant pathway flow and d) single-phase gas flow through macro-fractures in the soil. The gas production rate and the corresponding gas pressure accumulation affect the clay material behaviour and its properties such as air entry value. For the safe design of the EBS system, computational models need to account for the identified transport mechanisms. This study presents an enhancement in the finite element code Thebes [1, 2] that replicates the observed increase in permeability at higher gas pressures, e.g. due to pore dilatancy and gas fracture as proposed by Xu et al. [3]. The formulation links permeability to gas pressure and threshold/critical pressure. For model validation, the study utilizes a gas injection experiment carried out in IfG (Institute for Rock Mechanics, Germany) on Opalinus Clay [4]. The results show a good fit against the measurements while giving insight into gas flow through clays.


Introduction
Geological nuclear waste repositories can produce multiple gases through anaerobic corrosion of metal structures or aerobic and organic activities [5,6]. Depending upon the production rate and Engineered Barrier System (EBS) material properties gas pressure may lead to fractures in the barrier layers that would impair the repository's ability to isolate hazardous waste by facilitating radioactive gas or radionuclide transport processes. Preventing such failure requires an efficient quantification of gaseous accumulation and flow process.  Figure 1 illustrates four gas flow phases in soils such as bentonite or claystone, a key component of the barrier layer system [7,8] at progressively higher pressures: a) Gas diffusion phase occurs under low gas pressures, b) Two-phase flow that is both diffusive and adjective occurs when gas entry pressure value is exceeded, c) Flow through dilatant pathway happens when gas pressure exceeds minimum in-situ stress, and d) Single phase gas flow initiates at soil fracturing due to high gas pressures. In continuum approaches like FEM modelling of the transitions between those gas flow phases changes remains challenging. Studies implement various approaches to simulate the preferential pathway generation in soil. Simple hydraulic or hydromechanical models use permeability enhancement formulation that is dependent on a threshold gas pressure or external/confining pressure component [3,9,10]. The more advanced models employ hydromechanical coupled approaches such as damage models [3,11,12].
Further, studies such as Olivella and Alonso [14] and Chittenden et al. [15] represent key experimental phenomena of fissures and fracture gas flow [16,17] by using mixed approaches like embedding tubes or plates in the continuum soil domain. This allows the models to achieve single-phase gas flow at onset conditions (critical pressure/stress state). However, modelling the interaction between the soil domain and the capillaries can be challenging. Additionally, the opening and closing conditions of the capillaries largely remain the same as in the previous models.
The present study aims at enhancing the FEM code Thebes to account for preferential gas flow. Originally the code is designed to model the response of unsaturated (expansive) materials to thermal, hydraulic and mechanical (THM) loading [1,2]. Now, Thebes includes a gas pressure-dependent permeability  [3]. In comparison to other approaches, the formulation replicates gas dilatant and fracture flow behaviour with fewer modifications. Moreover, for future work, its simplicity allows for easy incorporation of random field theory and stochastic analysis [18] to simulate more practical scenarios.
To validate the newly implemented model, this study employs a gas injection experiment from IfG (Institute for Rock Mechanics, Germany) on Opalinus Clay [4]. The experiment uses stepwise gas load increments that gradually lead to gas breakthrough and clay fracture. Thebes results match well with the experimental data and OpenGeoSys codes [3] while offering insight into gas flow through clays. We also consider different coupling scenarios (gas-only vs water-gas coupling), along with mesh and time step sensitivity.

Governing equations
The Finite Element Method code Thebes assumes a three-component representation of the material: a) air, b) soil, and c) water in 3 phases (gas, liquid and solid). To perform thermal-hydraulic and mechanical coupled analysis the framework employs heat conservation, compositional method for mass conservation, and balances of mechanical forces.
Thebes key aspects are as follows: a) Advection flow -Darcy's law, b) Diffusion -Philip and Vries [19] for vapours and Fick's law for other gaseous species, c) Heat migration -Fourier's heat law, d) water-vapour phase change, and e) Mechanical behaviour of unsaturated expansive clay by adapting modified Barcelona Basic Model (BBM) that additionally considers the temperature effects [20,21,22].
Herein, the paper shows only the relevant part of the framework that is active during the gas fracture test simulation (water-gas coupling). It includes mass balances of water and gas (Sect. 2.1 and 2.2), water retention curve, Darcy's Law and relative permeability functions (Sect. 2.3). For full details on the Thebes refer to Abed and Sołowski [1, 2].

Water mass balance
shows the mass balance of water component. Here, i w  (i = l, g) are the water density in the liquid phase (l) and vapour density in the gas phase (g). i S (i = l, g) is the degree of saturation for liquid (l) and gas (g), respectively. jT  (j = s, w) are the coefficients of volumetric thermal expansion of solids (s) and water (w). wp  is the coefficient of water compressibility, T is the temperature, w M is the molar mass of vapour, g is the gravity and R is the universal gas constant. i q (i = l, g) are the water (l) and vapour (g) fluxes. k h (k = w, g) are the water (w) and gas pressure (g) heads.  is the matric suction head and g w j represents vapour diffusion.
 is a function of the pressure pw and the temperature T [23]: Note that the present work ignores vapour contribution. Equation 3 shows the gas mass balance. Thebes accounts for gas solubility by employing Henry's volumetric coefficient (H) of solubility [24]. Additionally, the sum of all the diffusive species in a phase satisfies the condition 0 mi kk j = (k components and m phases), implying

Air mass balance
Where, a  is the dry air density (same in liquid and gas phase), which follows the ideal gas law (see, Eq. 4). a M signifies the molar mass of air/gas.

Darcy's Law, retention curve and permeability functions
The current work assumes only bulk pressure flow and ignores diffusion of fluids for the gas fracture test. The governing Darcy's Law equation in the direction of gravity for i th phase (gas or liquid phase) is written as: Several models exist for soil water retention behaviour. The current work employs a different suction and gas entry pressure dependent retention and relative permeability function (Eq. 7 -9) than in the early version of Thebes (Eq. 6) [25,26], to accommodate claystone gas fracture test simulation [3,25,27].
Where, e P is the gas entry pressure. m is the retention curve fitting parameter of van Genuchten model.
Finally, Xu et al. [3] proposes two approaches for gas intrinsic permeability enhancement to simulate dilatant and micro-fracture gas flow behaviour in soils. Thebes adopts its pressure-dependent permeability enhancement approach (See Eq. 10). The second method uses the concept of damage model (See, Eq. 11) where permeability is a function of volumetric () vol  and plastic strain () p  . The formulation assumes that the permeability increases rapidly during plasticity. Additionally, it considers two separate functions to define permeability enhancement under compaction or extension.
Where g p [MPa] is the gas pressure,

Gas fracture experiment and Thebes verification
Thebes replicates gas transport experiments on Opalinus clays from IfG [4]. The test uses a claystone sample of 150.45 mm in length and 73.59 mm in diameter. Further, it employs nitrogen gas injection from the bottom borehole in incremental steps from 1 to 3.5 MPa (see, Figures 2 and 3) while maintaining atmospheric conditions (0.1 MPa) at the outlet borehole. The experiment lasts 4 hrs, maintaining hydrostatical loading conditions and 3 MPa confining pressures (see, Figure 3) by surrounding the sample with metal plates and rubber sealing. The experiment simulation uses axisymmetric geometry (see, Figure 2) with four-noded quadrilateral elements (1st order). Further, along with mesh and time step sensitivity analysis, the study considers two scenarios a) single gas flow model at constant saturation, and b) coupled water -gas flow simulations. Mesh sensitivity analysis uses 4000 and 8000 elements for gas-only cases, while it is 8000 and 45000 elements for water-gas coupling. We perform the time step sensitivity analysis using 715 and 3887 steps only for the water-gas coupling case (45000 elements). Further, the simulations apply similar conditions as Xu et al. [3] in OpenGeosys code. The model keeps all the boundaries closed to fluid flow (gas and water) except at the inlet and outlet. It applies the gas pressure as previously specified and assumes 0.9 initial saturation. Xu et al. [3] further note that the saturation remains unobserved during the experiment, and no water flows out. Therefore, the saturation likely varies in the simulation. Since the hydraulic boundary remains unclear in the study, the gas-only simulation uses constant saturation of 0.9, while the water-gas coupled case maintains constant hydraulic pressure at the injection and extraction boreholes, such that it represents 0.9 saturation at atmospheric gas pressure. The results show that saturation levels in claystone remain between 0.7 and 0.9 during the simulation. The model also assumes soil anisotropy due to the horizontal orientation of the bedding plane in claystone Xu et al. [3] . Hence, it uses ten times higher permeability along the bedding plane than in the normal direction. See Table 1 for the list of modelling parameters. Table 1. Parameters for modelling Gas fracture experiment [3].

Parameters values
Initial  Figure 3 indicates that all the simulations replicate the measured gas flow qualitatively well and follow the same pattern of four gas flow stages: no outflow until gas pressure reaches 2.5 MPa, moderate flow with the gas pressure raising to 3 MPa, significant flow at 3.5 MPa of gas injection pressure, likely representing flow due to micro-fractures and decline after the gas injection pressure drops to 1 MPa. Xu et al. [3] calculation with the deformation predict a too rapid drop in the outflow when the experiment reduces the gas pressure, while the gas pressure dependent approach gives improved results. However, the mechanically uncoupled, gas pressure dependent approach leads to worse results when mechanical loading cycles on the gas flow are considered [3].
Thebes was successful in matching the flow rate in both approaches (i.e. single gas flow with 8000 finite elements and coupled water-gas flow with 45000 finite elements). The difference in the results is marginal everywhere except at the gas flow peak value, where the results are approximately 0.7ml/min less than the experiment. Fig. 3. Thebes gas flow rate results vs Gas fracture experiment [3,4]  In the case of gas-only analysis, increasing the mesh density twice (4000 to 8000 elements) leads to the peak gas flow increases by 1.9 ml/min. However, it slightly overestimates the results before the onset of gas fracture. It is reasonable to assume that another round of increase in mesh density could bring the peak values closer to the experimental data, but the differences before the gas fracture state will increase. Hydraulic-gas coupling requires many more elements to reach accurate results, with the results converge at 45000 elements producing the closest fit to the experiment. From the coarser mesh scenario, the peak gas flow value improves by 3 ml/min. Furthermore, changing the time step shows negligible effects on the results. Overall, examining all the cases suggests that the water-gas coupling increases the accuracy of prediction especially before the peak flow and hence is important to take into account.

Conclusion
The paper presents an extension of the framework of FEM code Thebes so it accounts for gas pressuredependent preferential gas flow paths. For this, a gas pressure-dependent permeability enhancement approach is adapted from Xu et al. [3]. The successful replication of the gas fracture experiment results on Oplanius clay [4] validates the proposed formulation. Furthermore, the study highlights the effect of water-gas coupled modelling over gas-only analysis of the experiment. The result suggests that the water-gas coupled analysis allows the soil to dry, providing a closer representation of reality, while the gas-only analysis slightly overestimates gas outflow rates. Mesh sensitivity analysis shows a significant increase in gas flow values, whereas the time step change does not affect the outcomes. Moreover, results show Thebes capability to replicate four stages of gas flow from negligible to two-phase and preferential micro-fracture flows.
In general, the permeability enhancement approach is easy to incorporate and proves effective in this study. However, the method is mechanically uncoupled and cannot replicate the change in gas flow patterns due to an increase in confining pressure [3]. Hence, the formulation needs modification to account for mechanical conditions. Further, the model requires more test cases of gas flow experiments in different soils. Even though the method is relatively more empirical than some other advanced models [14,15], its simplicity and low numerical overhead pave the way for more novel work, such as permeability function application with random field theory and stochastic analysis [18] applied to larger scale practical problems.