3D modelling of gas injection tests on FEBEX material: incorporating heterogeneity effects

A

relative permeability constant

$a$

spacing of the fractures (m)

$b0$

initial aperture (m)

$bmax$

maximum aperture (m)

C

parameter in $P0ϕ$

D

parameter in $λϕ$

$Dαi$

diffusion coefficient of the phase α (m²s⁻¹)

$Dα′$

hydrodynamic dispersion coefficient (m²s⁻¹)

$dl$

longitudinal dispersivity (m)

$dt$

transversal dispersivity (m)

$e$

void ratio

$G$

plastic potential

$g$

gravitational constant (ms⁻²)

$iαi$

non-advective mass flux of species $i$ in phase α (air or water vapour) (kgs⁻¹m⁻²)

$k$

intrinsic permeability of the phase α (m²)

$kCubic$

strain dependent gas permeability (m²)

$keffl, keffg$

effective permeability of gas and liquid phase (m²)

$kint$

intrinsic permeability of gas (m²)

$kKozeny$

porosity dependent permeability (m²)

$krl$⁠, $krg$⁠, $krα$

relative permeability in phase α (l:liquid, g:gas)

$M$

critical state line parameter

n

relative permeability power

$Pg$

gas pressure (MPa)

$Pl$

liquid pressure (MPa)

$P0$

air entry pressure (MPa)

$P0ε$

strain dependent air entry pressure (MPa)

$P0ϕ$

porosity dependent air entry pressure (MPa)

$p$

mean net stress (MPa)

$pc$

reference stress (MPa)

$ps$

tensile strength in saturated conditions (MPa)

$p′$

mean effective stress (MPa)

$p0$

preconsolidation stress

$p0*$

preconsolidation stress for saturated conditions (MPa)

$q$

deviatoric stress (MPa)

$qα$

advective flux in phase α (ms⁻¹)

$r$

parameter defining the maximum soil stiffness

$Sl$⁠, $Sg$

degree of liquid and gas saturation

$s$

suction (MPa)

$α$

non-associativity parameter

$β$

parameter controlling the rate of increase of the soil stiffness with suction (MPa⁻¹)

$ε0$

reference strain to calculate aperture variations

$εve$, $εvp$

elastic and plastic volumetric deformation

$κi$

volumetric stiffness parameter under changes of mean stress

$κs$

volumetric stiffness parameter under changes of suction

$λ0$

slope of void ratio – mean stress curve at zero suction

$μα$

dynamic viscosity (Pa s)

$ρα$

phase density (kgm⁻³)

$σ′$

effective stress (MPa)

$τ$

tortuosity coefficient

$ϕ$

porosity

$ψ$

shape parameter

$ψϕ$

porosity dependent shape parameter

$ωαi$

mass fraction of species $i$ in phase α

Geological disposal is regarded as a sustainable solution for the long-term management of radioactive spent fuel (Posiva, 2021). Both thermo-hydro-mechanical (THM) modelling and the design of a final nuclear waste repository – considering safety requirements and performance targets established for buffer and backfill materials in terms of long-term THM response – pose a significant challenge, as discussed in Toprak et al. (2024). Besides the THM processes, gas transport mechanisms will take place in a disposal system during the post-closure phase of a nuclear waste repository. Gas migration is a critical component within the safety assessment of deep geological repositories (DGRs) in low-permeability formations (Nagra, 2008).

Gas migration into clay barriers was investigated in various international projects, such as the FORGE project (2009–2013), the DECOVALEX project (2019–2023), and finally the EURAD-Gas project (2018–2024). In the FORGE project, gas migration issues in repository performance assessment were investigated (Norris, 2015).

In the DECOVALEX-2019 Project (Tamayo-Mas et al., 2021), various types of modelling approaches were developed. One of the biggest challenges, as explained by Tamayo-Mas et al. (2024), was to characterise and localise dilatancy-controlled flow in the models performed by several teams. As concluded by Tamayo-Mas et al. (2024), heterogeneity might provide one possible route to represent localisation of flow in continuum models, but the distribution functions could not be physically justified, as these functions were arbitrary.

EURAD-GAS (EURAD, 2024) was built on the outcomes of FORGE and focused on the mechanistic understanding of gas transport in clay material. The main objectives of this project were to (i) improve the mechanistic understanding of gas transport processes in engineered clay materials and their couplings with the mechanical behaviour and (ii) evaluate the gas transport regimes that can be active at the scale of a geological disposal system and their potential impact on barrier integrity and repository performance.

This paper shows hydro-mechanical and gas transport (HM-G) modelling results of successive gas injection tests (Gutiérrez-Rodrigo, 2018) that were part of an experimental campaign submitted to the EURAD-GAS project aimed at understanding the mechanisms controlling advective gas flow in the Spanish reference barrier material (FEBEX).

The objective of the experiment was to investigate gas advective flux in FEBEX material. The modelling work aimed to simulate the formation of dilatant pathways during advective gas flux in FEBEX bentonite. In addition, a key objective of this study was to simulate successive gas injection tests, incorporating hydration stages to capture the coupled hydro-mechanical-gas transport processes occurring throughout the experiment.

Figure 1 outlines the methodology employed in this study, which involved six major modelling steps. The calibration of hydraulic (Sections 3.1 and 3.2) and gas transport (Sections 3.3 and 3.4) model parameters (Step 1) was performed in accordance with Lloret et al. (2004), Lloret and Villar (2007), and Gutiérrez-Rodrigo et al. (2015, 2021). Following the calibration of hydro-gas model parameters, a simulation of the swelling pressure test (Step 2) was conducted to establish mechanical model parameters (Step 3) as detailed in Section 3.5.

Subsequently, HM-G model parameters were refined through the processing of test data and back-calculations from a gas injection test. Three-dimensional (3D) models incorporating preferential pathways (Models A and B, Section 2) and heterogeneous porosity distribution (Model C, Section 2) were developed (Step 4). In these models, the injected volume of gas was monitored (Step 5), and the relative gas permeability of both the filters and the sample was correlated. Finally, the back pressure measured during testing was compared with model predictions (Step 6), and further calibration of cubic law parameters for the dilatant pathways was undertaken. As described in Section 3.4, the primary gas flow mechanism is advective transport of gas in its free state. Therefore, special emphasis was placed on modelling gas advective flow. However, the model is also capable of reproducing other gas transport mechanisms, such as diffusion and the advection of dissolved gas into the liquid phase. In Section 4, a comparison of the proportions of total mass fluxes is presented.

There is a variety of numerical approaches to model gas injection tests on compacted clays. Two-dimensional (2D) modelling of gas injection tests on initially saturated clays, with the mechanical behaviour predominantly elastic, was reported by Gerard et al. (2014), Guo and Fall (2018), Radeisen et al. (2023), and Mo et al. (2024). 3D modelling (the mechanical part is mainly elastic) of gas injection tests on initially saturated clays was reported in Damians et al. (2020), Tamayo-Mas et al. (2024), and EURAD (2024).

While previous contributions include modelling of the gas breakthrough (BT) behaviour of compacted bentonites primarily under saturated conditions and in 2D geometries within an elastic regime, the underlying mechanisms, such as 3D heterogeneity and the successive gas injection procedure, remain inadequately represented. Under repository-like conditions, the processes of gas generation, accumulation, and BT typically occur in a sequential manner. In the DGR, the gas generated will migrate in a cyclic manner, regulated by the opening and closure of pathways. Consequently, the modelling of successive gas injection tests, encompassing all stages from saturation to gas injection, is crucial for accurately assessing the performance of DGRs. BT pressures are associated with the actual dry density achieved at the end of the saturation process, rather than the initial dry density (Radeisen et al., 2023). Therefore, the simulation of saturation and re-saturation stages preceding gas injection becomes increasingly important, as both dry density and, consequently, permeability vary during the saturation process. Furthermore, 2D geometries, whether plane or axisymmetric, are not optimal for simulating gas flow in compacted clays due to the asymmetric nature of heterogeneity and the development of dilatant pathways.

This paper demonstrates the feasibility of modelling successive gas injection tests that might correspond to repository-like conditions, encompassing all stages: saturation, gas injection, development and closure of dilatant pathways, and re-saturation.

Analysing the impact of heterogeneity and incorporating it into numerical models is a challenging task. Heterogeneity may be associated with a single parameter or multiple parameters, depending on the nature of the problem. In engineering, and particularly in geotechnical engineering, the influence of heterogeneity spans several parameters – for example, the effective thermal conductivity of heterogeneous materials (Wang et al., 2008), inter-well transmissivity in heterogeneous aquifers (Desbarats, 1993), and the effect of heterogeneity on the strength characterisation of rock (Tang et al., 2007).

Random finite-element methods (FEM) that account for heterogeneity have been applied in the analysis of settlements in geotechnical applications (Griffiths and Fenton, 2009) and in slope stability assessments (Li et al., 2016). Heterogeneity in dry density observed in full-scale experiments (Sakaki et al., 2023), as well as in coupled THM processes (Watanabe et al., 2010; Rodriguez-Dono et al., 2023), has been linked to flow-related problems.

In this study, a framework and methodology were developed for determining the HM-G properties of compacted clay materials based on the simulation of laboratory-scale tests. This approach allows for the configuration of multiscale heterogeneity (such as air entry pressure, gas BT pressure, swelling pressure, porosity, and gas permeability) within a 3D model domain. The methodology established in this study can be applied to configure and replicate similar laboratory-scale tests, particularly in the context of nuclear waste repository design, where accounting for multiscale heterogeneity is essential.

Although the models presented in this study incorporate validated HM-G parameters through the processing and simulation of laboratory-scale tests, the generated heterogeneous permeability and porosity zones exhibit an arbitrary (non-localised) distribution. Improved instrumentation, sensor data, and advanced imaging techniques (such as the use of transparent walls) could enhance the definition of heterogeneous permeability zones (Wiseall et al., 2015). Another limitation is that, in contrast to the hydro-mechanical (HM) model parameters, the cubic law parameters (Section 3.3) calibrated in this paper are test-specific parameters. In any new gas injection test performed in FEBEX bentonite, the directions and magnitude of micro-apertures will vary. Consequently, cubic law parameters must be recalibrated based on back-calculations of the new experiment.

The sample used in the test was FEBEX bentonite, with a diameter of 38 mm and an initial height of 20 mm. The initial degree of saturation was 81%, corresponding to an initial water content of 15.3%. Since this is a constant volume test (Figure 2(a)), all boundaries are confined, thereby preventing any displacement. The specifications and the details of the gas BT experiments on FEBEX material can be found in Gutiérrez-Rodrigo (2018) and Gutiérrez-Rodrigo et al. (2021).

Before the gas (nitrogen at room temperature) injection tests, the bentonite was fully saturated with deionised water in non-deformable steel cylindrical cells. To optimise computational efficiency, the saturation period was fixed at 60 days across all simulations. Full saturation was achieved within this period in all cases. After saturation, the filters used for the process were replaced with dry ones. The gas cylinders, equipped with pressure transducers, were then connected to both cells for the gas injection tests (Figure 2(b)). The swelling pressure developed during the hydration of FEBEX bentonite was deemed sufficient to prevent gas flow between the bentonite and the cell wall. An initial pressure of 400 kPa was applied in the upstream cylinder, while a vacuum of 1 kPa was set in the downstream cylinder. The pressure in the upstream cylinder was gradually increased by 200 kPa increments over 24-h periods until gas BT occurred. The back pressure was monitored in the downstream cylinder. The moment when gas instantaneously crosses the saturated material is referred to as ‘BT’. The pressure required to achieve this BT is termed the BT pressure. After the first BT, any remaining gas in the sample was drained through the downstream cylinder. Once steady-state conditions were reached (with constant pressure at the upper boundary), a second BT test was conducted on the same sample.

Table 1 summarises the phases in the test and the magnitude of gas BT and back pressures obtained in the corresponding phases. GID as a CAD system and Code_Bright as an FEM programme have been used to simulate gas BT tests on FEBEX bentonite material.

The injector system has been introduced to the model as a boundary condition. Prescribed gas pressure was applied at one end of the sample. Model calibration was conducted by comparing the resulting back pressures with experimental data.

A total of three models were developed, each incorporating distinct heterogeneity configurations (Figure 3). Two distinct levels of heterogeneity are observed during gas injection testing. The first level is linked to variations in dry density (or porosity) during the hydration phase. The second level of heterogeneity concerns the effective permeability of gas. Models A and B represent arbitrary realisations of heterogeneity in terms of gas permeability, with Model A incorporating preferential flow paths and Model B combining both preferential paths and random heterogeneous zones. Model C, constructed as a representative scenario (not derived from specific laboratory observations), incorporates an initial heterogeneous porosity distribution. The evolution of gas permeability in relation to strain development and desaturation, along with the magnitude of permeabilities (Figure 9) in preferential pathways, is detailed in Section 3.3.

In Models A and B, the initial porosity was kept constant, corresponding to a dry density of 1.8 g/cm³, consistent with the experiment. In contrast, Model C introduces three discrete levels of initial dry density, each occupying equal volumes. Furthermore, Model C adopts the same distribution of permeability zones – comprising both connected pathways and unconnected sections employed in Model B. The performance targets established for buffer materials in a future DGR (Posiva, 2021) may range from 1.6 to 1.8 g/cm³. In addition, laboratory test data for these dry densities are available (Gutiérrez-Rodrigo et al., 2021). Consequently, the modelling work encompassed this range of dry densities to ensure consistency with experimental data and performance requirements.

All models have a 3D full geometry where porous stones have been simulated as separate materials. The geometry has a structured mesh, and the number of hexahedral elements is 960. There are 1197 nodes.

The hydraulic boundary conditions, along with time stepping and the prescribed gas pressure at the upper and lower boundaries, are illustrated in Figure 4(a). The total inflow from the upper boundary, where gas injection occurs, and the corresponding accumulated gas volume in the model are presented in Figure 4(b). At the conclusion of the first gas BT, both in the experiment and the models (A, B, and C), the total injected gas volume was approximately 0.29 cm³. The modelling approach involves applying gas pressure from the upper boundary (Figure 4(a)) and comparing the resulting back pressures.

The boundary conditions and porous stone parameters were held consistent across all models shown in Figure 3. As gas injection occurs by way of the porous stone, neither the configuration of the model nor the bentonite material parameters significantly affect the resulting mass flow rates or the injected volume. Instead, the hydro-gas parameters of the porous stone (see Section 3) predominantly control the mass flow rate and the cumulative volume within the system.

In this section, material model calibrations based on laboratory experiments (water retention curve (WRC), permeability tests, and swelling pressure tests) are described. In intact clay, where free gas advection does not occur, gas permeability and the WRC are dependent on porosity. However, in dilatant pathways, gas permeability and WRC are strain-dependent, with the aperture sizes in these pathways varying significantly. Therefore, the model must incorporate at least two distinct gas flow pathways, each characterised by different aperture properties.

The equations governing the HM-G response of clay-based materials used in CODE_BRIGHT (Olivella and Vaunat, 2023) are provided, and corresponding material model parameters are listed for each process. Calibration of hydro-gas and mechanical parameters based on laboratory experiments provides important information to guide modelling and to prepare a proper modelling configuration leading to a more realistic design.

Laboratory experiments on FEBEX material (Villar, 2007; Lloret et al., 2004) show that the permeability of bentonite mixtures depends on the bentonite content and on the dry density. The WRC and permeability of bentonite-based materials (Toprak et al., 2017; Toprak et al., 2020) are crucial parameters in THM modelling of engineered barrier systems. The gas transport through buffer and backfill materials in nuclear waste repository design, as highlighted in EURAD (2024), requires special attention.

In CODE_BRIGHT, several model options are available for WRC (Equations 1–4), as well as water (Equations 5–7) and gas (Equations 8–12) permeability. Porosity- and strain-dependent hydro-gas transport models were developed by analysing test data from Lloret et al. (2004) and Gutiérrez-Rodrigo (2018, 2021) to more accurately simulate the HM-G transport behaviour of FEBEX bentonite. The calibrated model parameters are listed in Table 2. An additional list of symbols and their corresponding usage within the equations is provided in the  Appendix for clarity.

For the porous stone filters, a high permeability for both water and gas, set at 5 × 10⁻¹⁵ m², was employed. This elevated permeability, combined with a higher porosity (0.95), facilitates efficient water and gas flow into compacted bentonite during the hydration and gas injection stages. Concurrently, a high elastic modulus of 6300 MPa ensures that constant volume conditions are maintained. The gas permeability of the porous stones was calibrated based on test data, specifically targeting the injected gas volume during the first gas BT phase.

The WRC (Figure 5), on which the approach in the following is based, is the van Genuchten equation (van Genuchten, 1980):

With $P=P0σσ0$⁠, where P₀ is the air entry value at a certain temperature, σ₀ the water surface tension at that temperature, and σ the surface tension as a function of the temperature.

The parameters P and $ψ$ are determined for different porosities following the relations, where C and D are model parameters:

A porosity-dependent WRC was used in zones where there is no significant advection of gas in its free state (intact clay). A strain-dependent WRC (Figures 6(a) and 6(b)) is necessary to accurately simulate gas entry pressure in dilatant preferential pathways, where permeability follows the cubic law (Equation 8).

To accurately represent the reduction in air entry pressure during the advection of free gas through dilatant pathways, a strain-dependent WRC (Olivella and Alonso, 2008) is required.

In this framework, the air entry pressure (⁠$P0$⁠) is governed by $kint$ (Equation 8), which, in turn, is influenced by $kcubic$⁠. Since $kcubic$ exhibits strain-dependency (Equation 9), $P0$ also becomes strain-dependent:

Figure 6(a) presents a comparison between a representative strain-dependent WRC, calibrated with a mean strain rate of 0.1, for the dilatant pathways, and a porosity-dependent WRC for intact clay. The results show that air entry pressure decreases significantly when a strain-dependent WRC is considered. In all model configurations (A, B, and C), a total of three WRCs were employed. A porosity-dependent WRC was used to represent the intact clay, while two strain-dependent WRCs were applied along two distinct flow paths, as illustrated in Figure 6(a). During the gas transportation, different levels of apertures (Equation 9) develop. The magnitude of these apertures is independent of porosity. Two sets of $kCubic$ parameters, as outlined in Table 2, were proposed. Strain-dependent WRC was generated based on these parameter sets, as shown in Figure 6(b). Volumetric deformations concentrated along preferential pathways (logarithmic scale in Model A) at the conclusion of the first BT test are depicted in Figure 6(b). The air entry pressure is influenced by the strains induced during the successive gas injection test.

For a continuum medium, Kozeny’s model is defined as below, where $k0$ is reference intrinsic permeability (m²) and $ϕo$ is reference porosity.

The liquid phase permeability depends on the degree of saturation. $krl$ is the relative permeability for the liquid phase, and it is controlled by the degree of saturation (S_l) and parameters A and power n.

Finally, the effective permeability of the liquid phase $keffI$ (m²) is defined as a product of matrix permeability (⁠$kKozeny$⁠) and relative permeability (⁠$krl$⁠).

Test data from Lloret et al. (2004) were utilised to calibrate the permeability model parameters. As illustrated in Figure 7, liquid permeability is dependent on both porosity and the degree of saturation. For higher initial dry densities, the sample exhibits lower intrinsic permeability at the same level of saturation. Hydraulic conductivity reaches its maximum when full saturation is achieved, as all pores are filled with water, providing the maximum surface area for water flow. The initial state of the sample is also depicted in Figure 7, with the calibration following a single continuous path.

To model the propagation of preferential pathways during gas flow in a porous medium, a hydro-mechanical model was developed by Olivella and Alonso (2008). The cubic law is applied to establish a simple relationship between deformation and intrinsic permeability. Figure 8 illustrates the geometry of these apertures caused by gas flow. The mean aperture distance, denoted as ‘a’, represents the average distance between adjacent apertures. The aperture opening is represented by ‘b’, corresponding to the size of each aperture. The orientation of the apertures is characterised by ‘n’, while the height is denoted as ‘s’. In this study, a simplified permeability function for embedded fractures was employed (Olivella and Alonso, 2008), in which the strain-dependent permeability model remains independent of both ‘n’ and ‘s’.

The intrinsic gas permeability is defined as below, where $kKozeny$ depends on porosity, following Equation 5.

Permeability in the cubic law depends on aperture characteristics, where $b0$ corresponds to the initial aperture and $bmax$ is the maximum aperture.

The aperture, ‘b’ is variable with deformation (⁠$Δb$⁠), influenced by the mean aperture distance (a) and the difference in strain (⁠$ε-ε0$⁠) along the aperture’s normal direction (n).

The relative permeability for the gas phase, $krg$⁠, is given below, where $Sg$ is (⁠$1-Sl$⁠)

Finally, the gas effective permeability (⁠$keffg$⁠) is a product of intrinsic and relative gas permeability:

The test data for gas permeabilities were derived from Gutiérrez-Rodrigo et al. (2021). Different sets of $kCubic$ functions were prepared as shown in Figure 9.

In this study, some reference parameters were calibrated for two different magnitudes of $kCubic$ based on the back-calculations of the gas injection test.

In CODE_BRIGHT, four fluxes (Equations 13–16) are considered for gas flow under three distinct conditions, as depicted in Figure 10. Below a threshold known as the gas entry pressure, the primary transport mechanism is gas diffusion within pore water (Condition 1). When the gas entry pressure is exceeded, advection of free gas initiates (Condition 2). After gas BT, gas advection becomes the dominant transport mechanism (Condition 3). The liquid phase is modelled as a mixture of water and dissolved air, while the gas phase consists of dry air (free gas) and water vapour.

Diffusion of dissolved air and water vapour begins under Condition 2 (⁠$PgasBT>Pgas>Pentry$⁠) as shown in Figure 10. In this study, temperature was assumed to be constant. Vapour diffusion was neglected.

Diffusive (or non-advective) mass fluxes are usually described by Fick’s law (τ: tortuosity, ϕ: porosity, ρ_α: density of the related phase, S_α: degree of saturation of the related phase, D: diffusion coefficient of the related phase, and α: air or water vapour phase):

The objective of this study was to investigate advective gas flow in a compacted bentonite. Since the gas pressure increments were applied over short time intervals, gas diffusion was limited during the tests, and no direct measurements of gas diffusion were performed. Although the experiment itself was not designed to study gas diffusion, the numerical model is capable of simulating the diffusive component of gas transport. In this context, a gas diffusion coefficient (Table 2) was assumed, based on values of similar magnitude reported for comparable materials by Gonzalez-Blanco et al. (2016), Jacops et al. (2017), and EURAD (2024). As the test was dominated by advective flow, the assumed diffusion coefficient had only a minor impact on the overall results.

Under Condition 2 (⁠$PgasBT>Pgas>Pentry$⁠), hydrodynamic dispersion can be considered as non-advective flow. Hydrodynamic dispersion (⁠$dl$ is longitudinal dispersivity and $dt$ is transversal dispersivity) mass flux is computed by means of Fick’s law and written as:

As shown in Figure 10, gas accumulation lasts until reaching the gas BT pressure. Under Condition 3 (⁠$PgasBT≥Pgas$⁠), advection of free gas occurs as follows:

The intrinsic permeability $k$ (m²) for the related phase (liquid or gas) is calculated as above, where $μ$ (Pa.s) is the dynamic viscosity, $ρ$ (kg/m³) is the density of the phase, and $g$ (m/s²) is the gravitational constant. In the calculations of the injected gas volume into the sample (Figure 4(b)), the density of nitrogen gas was assumed to be 1.250 kg/m³. The calculations were performed at a constant temperature of 20°C. The dynamic viscosity of nitrogen was considered as 1.75 × 10⁻⁵Pa⋅s. Default values defined in CODE_BRIGHT were used for the density and viscosity of the water and air phases.

The Barcelona Basic Model (BBM) (Alonso et al., 1990) was used to simulate the swelling pressure development of FEBEX bentonite in the models (Section 4). BBM as an elasto-plastic model would not only capture the development of swelling pressures during the hydration phase but also simulate the potential irreversible strains induced by gas injection, particularly under a heterogeneous model configuration.

BBM is one of the most widely used constitutive models for the simulation of the compacted bentonite in the design of final spent nuclear fuel repositories (Gens et al., 2009; Sánchez et al., 2023; Toprak et al., 2024).

The calibrated loading-collapse (LC) curve in the suction (⁠$s$⁠) – mean effective stress (⁠$p′$⁠) –deviatoric stress (⁠$q$⁠) plane is presented in Figure 11. BBM parameters for FEBEX were previously reported in Lloret et al. (2004) and Gens et al. (2009). In this study, the BBM parameters were updated (Table 3) through simulations of a swelling pressure test, targeting a swelling pressure of 8.8 MPa. The swelling pressure test is a designed procedure aimed at calibrating material model parameters, taking into account the initial conditions of the sample, including initial suction and dry density. It should be noted that this test does not correspond to a physical test conducted on the sample itself.

In the updated BBM parameters, suction dependency on κ_i and κ_s was neglected due to computational cost and model limitations. The constraints of BBM parameters are discussed in Alcoverro et al. (2023).

Overall code specifications (basic features and governing equations, mass balance and energy equations) can be found in CODE_BRIGHT and Olivella et al. (1996).

BBM model approaches are summarised below. Here, a brief description of BBM associated with material model calibration is summarised (Equations 17–23). All the model specifications with details can be found in Alonso et al. (1990).

The volumetric compressive behaviour of the BBM is defined by the elastic functions implemented for considering highly expansive soils.

In BBM, $κ$ (logarithmic compliance with respect to changes in net pressure) and parameter $κs$ (logarithmic compliance with respect to changes in suction) are constant.

The yield surface (⁠$F$⁠) depends on stresses and suction and can be expressed using stress invariants. Here, $q$ represents the deviatoric stress, $p$ denotes the mean total stress, $α$ is the non-associativity parameter, and M is the slope of the critical state line. When $α=1$⁠, $F=G$ (associated plasticity). In this study, the associated plasticity was considered.

The increase in soil stiffness with suction is defined as:

where $pc$ is the reference stress; $p0*$ is the initial preconsolidation stress for saturated conditions; $λ(0)$ is the slope of void ratio in saturated conditions; $r$ defines the maximum soil stiffness; and $β$ controls the rate of increase of soil stiffness with suction.

In all simulations, full saturation was achieved during the 60-day hydration phase. Models A and B share the same initial dry density configuration, which is uniform, whereas Model C incorporates a spatially variable initial dry density, as shown in Figure 3(a). The impact of the predefined dry density configuration on the stress state at the end of the saturation phase is illustrated in Figure 12(a). The only difference between Models A and B lies in the distribution of gas permeability zones (Figure 3(b)). Since no dilatant pathways develop during the hydration phase, the distribution of the generated mean total stresses in Models A and B is identical due to their shared dry density configuration.

In contrast, Model C exhibits some local variations in the distribution of mean total stress, attributable to the initial heterogeneity in dry density. However, given the small sample size and constant volume conditions, these variations are not significant. In all models, lower stress magnitudes are observed near the porous stones. This occurs due to the interpolation method used by the code to calculate mean total stresses. Across all configurations, the mean total stress developed in the centre of the model domain was approximately 9 MPa.

The distribution of mean total stresses in Model A is presented in Figure 12(b). The mean total stresses achieved at the end of the saturation phase (60 days) and during the first BT (118 days) exhibit similar distributions and magnitudes. A reduction in stresses, due to the dismantling of the sample and its recovery during gas injection (90 days), is illustrated in Figure 12(b).

As illustrated in Figure 12(c), the evolution of computed total mean stresses in Model A exhibits three distinct peaks. The first peak (1) corresponds to the development of swelling pressure during the hydration phase, where full saturation was reached in the simulation. After hydration, gas injection commenced, leading up to the first gas BT pressure. The second peak (2) reflects the magnitude of the total mean stresses achieved at the end of the first BT. Following the first BT, the remaining gas was drained from the lower section, and the second phase of gas injection began, as described in Table 1.

There was no re-saturation of the sample following the first BT. The second BT (peak 3) occurred at a lower gas pressure than that required for the first BT. Gas flow was facilitated during the second BT due to increased permeability (Equation 9) and decreased air entry pressure (Equation 4) resulting from the first BT.

The gas pressure at the first BT was the highest, subsequently decreasing during the following gas BT tests (Table 1). This phenomenon can be attributed to the presence of voids filled with residual gas or damage to the fabric along the flow pathways (Horseman et al., 1999). As observed in this study, the follow-up BT pressures (Phases II and IV; Table 1) for FEBEX bentonite did not recover to their previous values, even after re-saturating the specimen (Phase III). Similar laboratory studies (Cui et al., 2022) indicate that gas BT processes may lead to permanent degradation of the sealing capacity.

Figure 13(a) compares the evolution of gas pressure in Model A for the upper and lower sections of the sample. Similar to Models B and C (Table 4), Model A slightly overestimated back pressures during Phase II. However, it efficiently reproduced the development of back pressure during Phase IV.

As illustrated in Figure 13(b), gas enters the compacted clay when the gas pressure (⁠$Pgas$⁠) exceeds the gas entry pressure (⁠$Pentry$⁠). As dilatant pathways developed during gas injection, the gas entry pressure significantly decreased (Equation 4), facilitating the accumulation of gas and the formation of new apertures within the compacted clay matrix. At the initial stages of gas injection, both dissolved advective flux (Equation 16) and non-advective fluxes (diffusion: Equation 13 and dispersion: Equation 14) were active. As shown in Figure 13(b), the diffusion of dissolved gas into the liquid phase was greater than the advection of free gas until the first BT occurred. Following the first BT, advection of free gas became the dominant flow mechanism, resulting in a slight desaturation at the considered location. The gas pressure dropped to 6.6 MPa after the first BT, which corresponds to the final gas pressure (⁠$PgasF$⁠) shown in Figure 13(b).

A comparison of back pressures between the experimental tests and the models is presented in Table 4. The evolution of back pressure for Models A, B, and C, along with the test data, is illustrated in Figure 14.

The evolution of back pressures is primarily influenced by the boundary conditions rather than the model configuration, resulting in similar trends and magnitudes across the three models, with some local differences. The models slightly overestimate back pressures during Phase II. However, they demonstrate improved accuracy in predicting back pressures during Phase IV, where gas BT pressures are lower.

The gas BT pressure is an intrinsic property of the material, directly associated with the magnitude of the swelling pressure. In contrast, back pressure is not a material property, and its magnitude can vary between tests depending on sample dimensions and the test procedure. As demonstrated by Gutiérrez-Rodrigo (2018), there is no direct correlation between gas BT pressure and back pressure. The presence of back pressure confirms end-to-end gas flow; however, its magnitude is not directly linked to the gas BT pressure. In gas injection tests conducted at varying dry densities – where different gas BT pressures were expected – similar back pressure magnitudes were observed. In other words, a higher gas BT pressure did not consistently correspond to a higher gas back pressure across all tests. A higher back pressure in Phase IV compared with Phase II might be associated with residual gas pressure. Although the gas was drained after the BT steps, some residual gas may have remained in the sample (Gonzalez-Blanco et al., 2024). This remaining gas could contribute to the increased back pressure observed in Phase IV. On the one hand, the residual gas pressure following re-saturation might decrease the gas BT pressure; on the contrary, it might increase back pressure. The models satisfactorily reproduce this behaviour.

Figure 15 compares the evolution of key parameters across the three models, including the degree of saturation (a), gas diffusion and dispersion (b), advection of free gas (c), and permeability (d) and porosity (e) at different intervals (Table 1). As depicted in Figure 15(a), the distribution of degree of saturation varied between the models at the end of the second gas BT event. Desaturation occurred along preferential pathways and/or specific sections due to gas-water exchange. The desaturation observed in the lower part of the models indicates that gas reached this section, resulting in completed gas BT. The magnitude of desaturation is governed by the WRC. The differing configurations of the three models result in variations in the spatial distribution of desaturated zones.

Gas diffusion and dispersion are active throughout the model domain. However, gas diffusion (Figure 15(b)) is less dominant in sections where gas advection is the primary flow mechanism, such as in the dilatant pathways. Since gas diffusion is porosity-dependent, its distribution in Model C is more heterogeneous compared with Models A and B. As the objective of the experiment was to investigate particularly gas advective flux in FEBEX, gas diffusion was limited by maintaining short duration for stress increment steps, as discussed in Section 2. While the injected volume of gas is comparable with test data (Figure 4), the proportion of gas advection to gas diffusion is merely representative (Figures 13(b), 15(a), and 15(b)). It is clear that the main transport path is through gas advection. However, under a different model parameter set, the proportion of gas fluxes might vary, while always maintaining the advective flux higher.

Figure 15(c) illustrates the distribution of gas advective fluxes across the three models, where the advective flux of free gas is concentrated along preferential pathways in Model A and corresponding sections in Models B and C. Permeability, following the cubic law (Equation 9), increased in sections where gas advection was the main flow mechanism. This increase in permeability, due to the formation of apertures during gas flow, is shown in Figure 15(d). Previous studies (e.g. Gutiérrez-Rodrigo et al., 2015; Horseman et al., 1999) indicate that the gas permeability of FEBEX and MX80 bentonites, after BT, may increase to as much as 6.4 × 10⁻¹⁹ m², which is two orders of magnitude higher than the water permeability of these bentonites. The modelling results presented in this study are consistent with these observations, showing gas permeability values along the gas advection path in the range of 10⁻¹⁹ m².

The desaturated zones, advective gas flux, and regions of increased permeability exhibit a direct correlation, indicating a coupled relationship between these phenomena. Notably, the presence of high-permeability zones in the lower part of the model domain confirms the advection of free gas through the lower section. This observation suggests that gas migration originated from the upper part of the sample by way of dilatant pathways, resulting in desaturation, particularly along the gas flow paths.

A summary of the impact of system heterogeneity on the modelling of gas BT tests for FEBEX material is provided in Table 5. This coupled process is based on several assumptions, considering varying levels of effective gas permeability, porosity, swelling pressure, gas BT pressure, and water retention and permeability functions.

Due to the hydration phase preceding the final gas injection, some sections of the FEBEX material (lower part) underwent swelling, while others (upper part) were compressed, leading to changes in porosity (Figure 15(e)). Gas BT pressures are more closely associated with the final dry density before gas injection rather than the initial dry density. The initial dry density plays an important role until full saturation is reached prior to gas injection. Gas BT pressure is influenced by both hydraulic and mechanical changes, including swelling, compression, and/or plasticity, occurring before gas injection. Nevertheless, this is a constant volume test, and while the initial porosity is set at 0.33, it can vary from 0.30 to 0.36 in different sections of the sample according to model response (Figure 15(e)). This variation must be considered to achieve a realistic design.

The configuration generated in this work demonstrates that the gas flow within compacted FEBEX material strongly depends on the specific conditions of the investigated system. Various HM-G aspects need to be taken into account to attain a comprehensive understanding of the localised gas flow in gas BT tests on compacted FEBEX material.

This paper presents a 3D HM-G simulation of a successive gas injection test conducted on FEBEX bentonite material. The model incorporates saturation processes to analyse the relationship between swelling pressure during hydration and gas BT pressure during the initial gas injection phase. A key feature of the model is the integration of both connected and unconnected dilatant pathways, which are employed to capture flow localisation and establish a robust numerical framework for hydro-mechanical coupling.

The study evaluates the impact of an initially heterogeneous porosity distribution on the HM-G response of compacted bentonite during successive gas injections. The experimental and modelling results reveal the complex behaviour of bentonite under gas transport, characterised by multiple interacting coupling mechanisms. The findings highlight the sensitivity of heterogeneous model configurations, particularly in terms of permeability and porosity, to water uptake and gas migration processes, affecting the distribution of gas permeability and fluxes.

The simulated results under various model configurations demonstrate strong alignment with experimental data on injected gas volume and back pressure (gas outflow rate). However, local variations in permeability, gas fluxes (both diffusive and advective), and saturation were observed across the models. Given the small sample volume, significant differences between configurations A (connected paths), B (connected paths + unconnected sections), and C (initial heterogenous porosity) were not evident.

This study successfully integrates BBM for mechanical behaviour with strain-dependent permeability (following the cubic law) and a strain-dependent WRC. It validates the feasibility of simulating all phases of a successive gas injection test, representative of repository-like conditions, within a 3D heterogeneous model domain that accounts for variations in dry density and gas permeability.

The calibration of cubic law parameters needs to be improved by incorporating, for instance, better instrumentation, higher-quality sensor data, advanced imaging techniques, and transparent walls. Despite the limitations of the localisation of gas flow pathways in the laboratory tests, the multiscale heterogeneity configuration (dry density and gas permeability) developed in this study can be utilised in laboratory or full-scale modelling of the gas injection experiments.

The modelling approach developed in this study for simulating 3D multiscale heterogeneity under successive hydration and gas injection conditions – incorporating advanced material models such as BBM, as well as porosity- and strain-dependent WRC and permeability functions – has been applied to one of the most comprehensive in-situ gas injection experiments: the gas permeable seal test (GAST). The theoretical framework and methodology established in this work served as the basis for the modelling of GAST, as presented in Toprak et al. (2025).

Financial support has been granted by EURAD. Tests were performed in the laboratories of CIEMAT.

Table 6 provides the nomenclature used in the study, detailing symbols, parameter definitions, and their associated equations.