Advective gas flow in bentonite: development and comparison of enhanced multi-phase numerical approaches

Geomechanics for Energy and the Environment 37 (2024) 100528 Available online 15 December 2023 2352-3808/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Advective gas flow in bentonite: Development and comparison of enhanced multi-phase numerical approaches E. Tamayo-Mas a , * , J.F. Harrington a , I.P. Damians b , S. Olivella b , E. Radeisen c , J. Rutqvist d , Y. Wang e a British Geological Survey, Keyworth, Nottingham, United Kingdom b Universitat Polit` ecnica de Catalunya⋅BarcelonaTech and International Centre for Numerical Methods in Engineering, Barcelona, Spain c Federal Institute for Geosciences and Natural Resources, Hanover, Germany d Lawrence Berkeley National Laboratory, Berkeley, CA, USA e Sandia National Laboratories, Albuquerque, NM, USA ARTICLE INFO Editors-in-Chief: Professor Lyesse Laloui and Professor Tomasz Hueckel Keywords: Dilatancy-controlled gas flow Enhanced numerical hydro-mechanical models Pathway flow Comparative study Bentonite gas tests DECOVALEX ABSTRACT Understanding the impact of repository gas, generated from degradation of waste and its interaction with the host rock, is essential when assessing the performance and safety function of long-term disposal systems for radioactive waste. Numerical models based on conventional multi-phase flow theory have historically been applied to predict the outcome and impact of gas flow on different repository components. However, they remain unable to describe the full complexity of the physical processes observed in water-saturated experiments (e.g., creation of dilatant pathways) and thus, the development of novel representations for their description is required when assessing fully saturated clay-based systems. This was the primary focus of Task A within the international cooperative project DECOVALEX-2019 (D-2019) and refinement of these approaches is the primary focus of this study (Task B in the current phase of DECOVALEX-2023). This paper summarises development of enhanced numerical representations of key processes and compares the performance of each model against high-quality laboratory test data. Experimental data reveals that gas percolation in water-saturated compacted bentonite is characterised by four key features: (i) a quiescence phase, followed by (ii) the gas breakthrough, which leads to a (iii) peak value, which is then followed by (iv) a negative decay. Three models based on the multiphase flow theory have been developed. These models can provide good initial values and reasonable responses for gas breakthrough (although some of them still predict a too-smooth response). Peak gas pressure values are in general reasonably well captured, although maximum radial stress differences are observed at 48 mm from the base of the sample. Here, numerical peak values of 12.8 MPa are predicted, whereas experimental values are about 11 MPa. These models are also capable of providing a reasonable representation of the negative pressure decay following peak pressure. However, other key specific features (such as the timing of gas breakthrough) still require a better representation. The model simulations and their comparison with experimental data show that these models need to be further improved with respect to model parameter calibration, the numerical representation of spatial heterogeneities in material properties and flow localisation, and the upscaling of the related physical processes and parameters. To further understand gas flow localisation, a new conceptual model has been developed, which shows that discrete channels can possibly be induced through the instability of gas-bentonite interface during gas injection, thus providing a new perspective for modeling gas percolation in low-permeability deformable media. 1. Introduction In all repository concepts for the geological disposal of radioactive waste, an engineered barrier system (EBS) is used to encapsulate the waste canister, or, to act as borehole or gallery seals. 1 These systems are often based on bentonite clays due to their low permeability 1–4 and high swelling capacity enabling the closure of engineering voids. 5–7 However, in all repository concepts gases will be generated through the corrosion * Corresponding author. E-mail address: [email protected] (E. Tamayo-Mas). Contents lists available at ScienceDirect Geomechanics for Energy and the Environment journal homepage: www.elsevier.com/locate/gete https://doi.org/10.1016/j.gete.2023.100528 Received 7 March 2023; Received in revised form 10 October 2023; Accepted 6 December 2023 Geomechanics for Energy and the Environment 37 (2024) 100528 2 of metallic materials (under anoxic conditions), the radioactive decay of waste and the radiolysis of water. 1, 8 These gases will move by the combined processes of molecular diffusion (governed by Fick’s law), solubility (governed by Henry’s law) and bulk advection (described by an extended form of Darcy’s law). 9 While the processes governing diffusion and solubility are well understood, the mechanisms controlling the advective movement of gas (as a discrete phase) in clay-based materials remain an active area of research. There is now a substantial body of evidence, spanning multiple decades, indicating that in these materials, gas migration occurs through the creation of dilatant pathways rather than by moving within the original porosity of the clay as conceptualised in Darcy’s law. 10–19 These pathways interact with the surrounding clay mass to introduce a complex hydro-mechanical coupling. 15, 17, 20 Understanding these processes is essential to accurately modelling repository gas migration and to assess its long-term impact on the safety function of the EBS. To this end, multiple researchers have recently developed new numerical algorithms within the framework of Biot’s consolidation theory, where features such as damage, plasticity, embedded fractures, etc., have been considered to include or explicitly represent dilatant gas flow in clay-based materials. 20–26 All these models are capable of capturing some of the main response features observed in laboratory studies, but they do not appear to effectively reproduce all of the main experimental observations that characterise dilatancy-controlled flow. Indeed, although most of them can satisfactorily simulate specific experimental features such as breakthrough time, maximum pressure/stress measurements and the dilation of the sample, there remain difficulties when trying to match the whole experimental history. Features such as the evolution of stress during and after breakthrough or the temporal variations in gas outflow are not well reproduced. 21 Other features such as the improvement of parameter calibration procedures, the further development of spatially- and temporally-dependent processes and the inclusion of pathways that behave highly dynamic and unstable, still remain uncertainties and thus, development of new and novel numerical representations for the quantitative treatment of gas in clay-based repository systems are still required. This was the primary focus of Task A in the DECOVALEX-2019 (D-2019) project, 21 in which 8 teams from around the world attempted to model advective gas flow in 1D and 3D experiments performed on compacted bentonite samples under controlled laboratory conditions. 15, 27 In Task A D-2019, four types of modelling approaches were developed: (i) two-phase flow models incorporating a range of different mechanical deformation behaviours, (ii) enhanced two-phase flow models in which fractures were embedded within a plastic material (continuous techniques) or incorporated into the model using a rigid-body-spring network (discrete approaches), (iii) a single-phase model incorporating a creep damage function in which only gas flow was considered, and (iv) a conceptual approach used to examine the chaotic nature of gas flow. In contrast to previous international gas projects such as EVEGAS 28–30 or GAMBIT, 31–33 where some model parameters were heuristically adjusted to overcome the poorly characterised couplings between the stress field and gas and water pressures, the D-2019 models featured more robust hydro-mechanical couplings based on pre-defined physical quantities. However, some of the important underlying physics (e.g., creation of dilatant pathways) associated with advective gas flow were poorly described and therefore unable to represent the full complexity of the processes in these low-permeability materials. Several concerns were raised in Task A D-2019 as some key features in the modelling of advective gas were still unclear: 1. Parameter calibration and model constraints: model complexity was significantly different among the proposed strategies and some models were clearly over-parameterised. Marked differences were also found in the calibration outcomes. Indeed, both the number of the calibrated parameters and the experimental outputs used to calibrate them were significantly different between teams. While, Fig. 1. Photos of apparatus and the test sample installed within the bore of the vessel. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 3 some of these parameters (e.g., tensile strength, swelling pressure etc.) had a clear measurable physical meaning, others (e.g., damage smoothing coefficients, capillary spacing etc.) were numerical constructs which had an indirect physical meaning and were often unmeasurable. Hence, their definition was complex and their extrapolation to other tests was found to be difficult thus hampering the use of models as predictive tools to assess gas movement. 2. Heterogeneity: two models included explicit representations of material heterogeneity, but the assumed distribution functions were arbitrary. Heterogeneity might provide one possible route to represent localisation of flow in continuum models but the distribution functions could not be physically justified and thus, teams concluded that further exploration was required. 3. Stochasticity: the experimental data from the 1D and the 3D gas injection tests exhibited a combination of stochastic and deterministic behaviours. Gas breakthrough (identified by the discharge of gas at the backpressure filter, see Fig. 1) occurred after a period of increasing gas pressure, when bulk gas flow was observed through the development of an emergent pathway. The instability and pathway switching observed in the 3D experiment before a main flow path was established, suggested that the precise timing of the gas breakthrough and associated gas flows could be semi-stochastic in nature. It is therefore important for a deterministic analysis to understand and distinguish between the key experimental features reproducible across all experiments and those that only occur in specific experiments. Therefore, being able to analyse and model similar high-quality experimental datasets was required to help give confidence in process understanding. 4. Upscaling: numerically upscaling (from laboratory to field scale) poses significant challenges, which go beyond the current study. However, a future companion paper will address these issues based on an examination of a full-scale test undertaken at the ¨ Asp¨ o Hard Rock Laboratory, Sweden. With these concerns in mind, it was concluded that enhanced numerical representations for the quantitative treatment of gas in claybased repository systems were required. This was the primary focus of Task B in DECOVALEX-2023 (D-2023). This paper summarises the outcomes of stages 1 and 2 with work conducted from May 2020 to January 2022 by the four modelling teams participating in the first part of the task: 1. BGR/UFZ (Germany): Federal Institute for Geosciences and Natural Resources and the Helmholtz Centre for Environmental Research. 2. LBNL (United States of America): Lawrence Berkeley National Laboratory. 3. SNL (United States of America): Sandia National Laboratories. 4. CIMNE-UPC/Andra (Spain/France): Universitat Polit` ecnica de Catalunya (International Centre for Numerical Methods in Engineering), funded by l ′ Agence nationale pour la gestion des des d´ echets radioactifs. First (Section 2), a comprehensive description of the apparatus used for the test is given. Details of the test history undertaken and the new, hitherto unpublished, results can be found. Second (Section 3), team models are collated and assessed against the above-mentioned test to Fig. 2. Schematic of apparatus. [A] shows the dimensions of the sample and filters. [B] shows the dimensions of the injection filter. [C] and [D] show the location of the axial sensors in contact with the base and top of sample respectively. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 4 analyse their performance (Section 4). Finally, a discussion of main findings and key conclusions are provided in Sections 5 and 6 respectively. It is not the intention of this paper to provide an exhaustive description of the individual contributions from each team, but rather give a technical overview and synthesis of key conclusions and results, highlighting their main differences and similarities (see Appendix A). 2. Experimental data As part of Task B, a new unpublished dataset was required to facilitate the comparison exercise. In previous stages of DECOVALEX-2019, teams had been given data at the outset with which to develop and calibrate their models. These experiments were based on data derived from tests performed on samples with a nominal length of 120 mm and nominal diameter of 60 mm. However, for the purposes of this numerical exercise (test FPR-21–004), the sample geometry was changed from a 2:1 ratio to a 1:1 ratio (sample of 60 mm in length and 60 mm diameter). As before, experiments were performed in a constant volume cell with material subject to rehydration first and then gas flow. 2.1. Experimental set-up In this test geometry, the specimen is volumetrically constrained, preventing bulk dilation of the clay in any direction. This BGS customdesigned apparatus, Fig. 1, has six main components: (i) a thickwalled, dual-closure pressure vessel; (ii) an injection pressure system; (iii) a backpressure system; (iv) 24 total stress gauges to measure radial and axial total stresses; (v) two filters for monitoring porewater pressure; and (6) a data acquisition system based around a National Instruments Compact Rio set-up. The pressure vessel comprised of a dual-closure tubular vessel manufactured from Invar (a nickel-iron alloy with a low thermal expansion coefficient, also known as Alloy 36) and was pressure-tested to 70 MPa. Each end-closure was secured by eight high tensile cap screws which could be used to apply a small pre-stress to the specimen if required. The vessel was mounted vertically with injection of gas through a rod mounted in the lower end-closure, Fig. 1, the dimensions of which are presented in Fig. 2. The 60 mm internal bore of the pressure vessel was honed and hardchromed to give a highly polished surface. Two pore pressure filters, labelled F1 and F2, Fig. 1, were mounted in the lower end-closure and used to provide local measurements for pore pressure during the course of the experiment. The total stress sensors were located in a geometric pattern, Table 1, with 4 radial arrays, each comprised of four sensors, spaced evenly along the sample. Each end-closure also contained four axial total stress sensors, mounted in a square configuration, visible in Fig. 1 and Fig. 2[C] and [D]. The central or injection filter was embedded at the end of a 6.4 mm diameter stainless steel tube and was used to inject helium (see Section 2.3) during gas testing. The end of the filter was profiled to match a standard twist drill to minimize voidage around the injection tip. Pressure and flow rate of test fluids were controlled using two ISCO- 260, Series D, syringe pumps, operating from a single digital control unit. Given the potential for gas leakage past the injection pump seal, a constant flow rate was developed by displacing gas from a pre-charged cylinder by injecting water. This also helped to ensure that the helium was saturated with water vapour prior to injection, reducing the potential for desiccation. A second pre-charge vessel was placed in the backpressure circuit to collect the gas as it was discharged from the sample. Test data was over-sampled at an acquisition rate of one scan per 20 s. Thereafter, data processing was based on every 10th datapoint, unless something of note occurred within the data. All pressure transducers were calibrated to an accredited laboratory standard and linear least-squares regression fits were applied to the data to ensure its comparability. 2.2. Geotechnical properties Geotechnical properties for the test sample FPR-21–004 were not available at the onset of the modelling exercise, as testing was still ongoing. As such, geotechnical data was provided by a second sample, FPR-20–030, manufactured in exactly the same way as that for FPR- 21–004, with data presented in Table 2. Data are calculated based on a grain density of 2.77 g/cm 3 . 2.3. Test history Details of each test stage are presented in Table 3. Sample FPR- 21–004 was subject to a series of test stages to promote swelling and resaturation (stages 1 and 2), equilibration (stage 3) gas injection (stage 4) and self-sealing (stage 5). 2.4. Results Following installation of the sample, a small backpressure of 250 kPa was simultaneously applied to the backpressure filter and both F1 and F2 filters, Fig. 1. On day 3.1 (test stage 2), pressure in all three filters was increased to 1.0 MPa, and the sample allowed to hydrate. The subsequent development of axial and radial stresses is shown in Fig. 3 and Table 1 Sensor locations. Sensor prefixed with the letter A=axial and R=radial. Radius equates to the centre line of the vessel running axially along its length. Height (z) is the distance from the base of the sample. Surface area relates to the circular size of the sensor/filter, see Fig. 1. Sensor name Rotation (degrees) Radius (mm) Height, z (mm) Surface area (mm 2 ) Sensor name Rotation (degrees) Radius (mm) Height, z (mm) Surface area (mm 2 ) A1 0 20 0 50.27 R10 90 30 36 50.27 A2 90 20 0 50.27 R11 180 30 36 50.27 A3 180 20 0 50.27 R12 270 30 36 50.27 A4 270 20 0 50.27 R13 45 30 48 50.27 R1 0 30 12 50.27 R14 135 30 48 50.27 R2 90 30 12 50.27 R15 225 30 48 50.27 R3 180 30 12 50.27 R16 315 30 48 50.27 R4 270 30 12 50.27 A5 0 20 60 50.27 R5 45 30 24 50.27 A6 90 20 60 50.27 R6 135 30 24 50.27 A7 180 20 60 50.27 R7 225 30 24 50.27 A8 270 20 60 50.27 R8 315 30 24 50.27 F1 135 24 0 28.27 R9 0 30 36 50.27 F2 315 24 0 28.27 Table 2 Geotechnical properties based on sample FPR-20–030. Moisture content Bulk density (kg/m 3 ) Dry density (kg/m 3 ) Void ratio Porosity Saturation 0.25 1989 1592 0.740 0.425 0.93 E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 5 Fig. 4. At the same time, water pressure was increased and the gas (helium) pressure in the injection filter was also incremented to match the change in water pressure. This was done in order to prevent accidental water flow into the filter and thus remove the possibility of slug flow (i.e., gas displacement of water from the filter ahead of gas entry into the sample) during gas injection (stage 4). Examination of the axial stress data, Fig. 3, shows a progressive increase in stress, axially throughout the sample. Unsurprisingly, stresses were systematically higher at the top of the core, adjacent to the large backpressure filter, where access to water was unlimited. At the base of the sample, access to water was limited to that available in the small filters F1 and F2, see average stresses A1-A4 in Fig. 3[B]. By the end of stage 2, this resulted in an average axial stress of 7.5 MPa, compared to 8.5 MPa at the top of the sample. Close inspection of the Table 3 Experimental stages and boundary conditions for test FPR-21–004. Stage number Stage type Start time (d) Injection pressure gas (kPa) Backpressure water (kPa) Filter F1 (kPa) Filter F2 (kPa) Gas displacement rate ( μ l/h) 1 Equilibration 0 0 250 250 250 - 2 Hydration 3.1 1000 1000 1000 1000 - 3 Equilibration 27.2 1000 1000 1000 1000 - 4 Gas ramp 28.4 2000 +1000 - - 180 5 Shut-in 84.9 - 1000 - - - Fig. 3. Test FPR-21–004 showing the development of axial stress during test stages 1 and 2. In [A] each sensor is shown and in [B] the average stress is shown: sensors A1 to A4 (located at the base of the sample) and A5 to A8 (positioned on the top face of the sample). Fig. 4. Test FPR-21–004 showing the development of radial stress during test stages 1 and 2. Sensors R1-R4 are closest to the base of the sample and R13-R16 are closest to the top of the sample, Table 1. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 6 data in Fig. 3 indicates that swelling was ongoing at the end of stage 2. However, based on the final 5 days of data from this stage, the rates of change were relatively small: 11.3 and 3.4 kPa/d for A1-A4 and A5-A8 respectively. The higher rate of change in the base of the sample relates to the aforementioned availability (or rather lack) of water. It is also noteworthy that, by the end of the stage, the variation in axial stress measured across each plane was very similar (0.69 MPa and 0.63 MPa for A1-A4 and A5-A8 respectively), see Table 4. The development of radial stress evolves in a similar manner to that of the axial stress. When the clay was exposed to water, radial stresses also rapidly increased, Fig. 4, approaching a well-defined asymptote by the end of the stage. Inspection of the data shows some variability, even across the same plane of measurement, such as the difference in values between sensors R9 to R12, or R13 to R16. Analysis of the data in Table 4, indicates larger variations in pressure exist as the distance to the points of measurement increase from the three sources of hydration. This results in some degree of heterogeneity in the resultant stress field at the end of stage 2, probably linked to a combination of incomplete homogenisation of the bentonite and subtle variations in the microstructure of the sample. The data in Fig. 5 indicates that average radial stress does not change along the length of the sample. This suggests that friction between the sample and vessel wall is likely to play a small (if any) role in the development of the stress on the scale of the experiment. The difference between axial stress at the top and bottom of the sample (Fig. 3) therefore relates to the availability of water (which is less at the base of the sample) suggesting incomplete homogenisation at this region of the sample. However, the length of time required to fully homogenise bentonite is unclear and remains a research priority within the European Community. As time was limited and the average values of radial stress, Fig. 5, were relatively similar across each measurement, the sample was deemed ready for gas testing. Following a brief period of further equilibration, gas testing (stage 4) began on day 28.4, with the injection pump set at a constant displacement rate of 180 μ l/h. Gas pressure gradually increased for the following 49.4 days, reaching a peak value of 12.36 MPa at day 77.8, see Fig. 6. This was followed by a spontaneous negative pressure transient leading to a quasi-steady state around day 84.9. At this point, the injection pump was stopped, test stage 5, and the pressure allowed to slowly decay. Gas pressure continued to decline until the test was stopped at day 117. Closer inspection of the axial stress data, see Fig. 7[A], suggests a small gas entry (i.e., onset of flow into the clay) event occurred around day 70.1, shortly after gas pressure exceeded the lowest value of axial stress. However, inspection of the outflow data, Fig. 7[B], shows no obvious sign of discharge. In the absence of a suitable sink, gas pressure therefore continued to increase, finally resulting in a major gas entry event around day 74.3 at a gas pressure of around 9.84 MPa. In the intervening time between the precursor and major entry events, the rate of change in all stress traces began to increase, providing clear evidence of a hydrodynamic effect caused by the coupling between gas pressure and axial stress. Further examination of the outflow data following the major gas entry event at day 74.3, Fig. 7[B], suggests a very small increase in outflow occurred around this time, possibly caused by compaction of the clay following major gas entry. However, gas breakthrough did not occur until day 75.5, from which it can be inferred that it took 1.2 days for sufficient pathways to develop and locate the backpressure filter. At this point, outflow rapidly increased and was associated with a series of spontaneous changes in the stress field, Fig. 7[A]. However, the initial development of permeability was short-lived, as conductive pathways began to spontaneously close and outflow quickly decreased. From day 76.1 to 76.5, outflow reduced to pre-major breakthrough levels and was accompanied by further increases in gas pressure and stresses within the sample, as gas pathways developed to reconnect with the backpressure filter. Towards the right margin of Fig. 7[B], outflow from the sample spontaneously increased for a second time at day 76.5 and was again associated with complex changes in the stress field. Fig. 8 shows the evolution of axial stress in response to the development of flow out of the sample. Gradual and abrupt changes in gas pressure and stress were observed that are interpreted as linking to the development of flow paths within the sample. This behaviour is consistent to previous experimental observations (e.g., 15 and 17 ). Following major gas breakthrough outflow, stress and gas pressure were highly coupled. In general, gas pressure remained slightly above that of axial stress, suggesting some degree of pressure drop along gas pathways. However, the orientation of the pathways to the direction of stress measurement strongly impacts the measured value of stress. Therefore, a detailed analysis of the stress response (axial and radial) is required and will be undertaken at a future date. Table 4 Minimum and maximum pressures at the end of stage 2 for test FPR-21–004. Sensors Max value at end of stage 2 (MPa) Min. value at end of stage 2 (MPa) Difference in pressure (MPa) A1-A4 7.86 7.17 0.69 R16- R13 8.71 7.79 0.92 R12-R9 8.90 7.32 1.58 R8-R5 8.83 7.82 1.01 R4-R1 8.75 7.88 0.87 A5-A8 8.70 8.07 0.63 Fig. 5. Test FPR-21–004 showing the development of radial stress at each plane of measurement during test stages 1 and 2. The position of each ‘array’ is in Table 1. Fig. 6. Test FPR-21–004 showing the development of gas pressure and axial stress, stage 5. The gradual increase in most sensors from day 28 to 70 can be attributed to continued hydration of the bentonite. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 7 While the peak gas pressure response was more rounded in shape compared to the that of outflow, peak gas pressure occurred prior to peak outflow. This suggests a degree of time-dependent behaviour in the development of gas pathways and thus permeability within the sample. However, by the end of the test stage at day 84.9, Fig. 8[B], inflow and outflow were converging, signifying the test was approaching a near steady state condition. However, close examination of the data shows small discrepancy between inflow and outflow data at day 84.9. As the system was not in true steady-state, it is not possible to identify the true origin of this discrepancy, but might stem from a small background leak. As soon as the injection pump was stopped, gas pressure, axial stress and outflow from the sample rapidly decreased, Fig. 9. For the first time, gas pressure dropped below axial stresses A2 and A5 through A8, suggesting depressurisation of some pathways and trapping of residual gas Fig. 7. Data from part of stage 4 of test FPR-21–004 with dotted lines denoting significant events. [A] expanded view of gas pressure and axial stress data showing precursor gas entry around day 70.1 followed by major gas entry at day 74.3. Line colours are the same as those on Fig. 6. [B] shows time averaged outflow data under experimental conditions with major gas outflow occurring at day 75.5. Fig. 8. Data from part of stage 4 of test FPR-21–004. [A] development of axial stress and [B] evolution of inflow and outflow data at STP, from day 75–85. Combined, the data illustrates the complex coupling between development of flow and accompanying changes in stress within the bentonite. Fig. 9. Shut-in, test stage 5, for sample FPR-21–004. [A] change in gas pressure and axial stress and [B] evolution of outflow (at STP) following the cessation of pumping at day 84.9 (the missing data at day 103 was caused by a problem with the data logger). E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 8 in others. From day 87, a more gradual decline in gas pressure and axial stress was observed. By the end of the stage, gas pressure was over 1.0 MPa lower than axial stress at sensors A5 and A6. Examination of Fig. 9 indicates the variation in total stress increases as the shut-in stage progresses. This is likely caused by the localised closure of gas pathways, resulting in the entrapment of gas within the sample, perturbating the stress field. The disconnect between outflow and changes in axial stress and gas pressure suggests complex patterns of flow drainage, which spatially evolved within the clay. This indicates that pathway closure can be rapid, as in the early stages of the shut-in response, or slow, as illustrated by the subsequent gradual decline in axial stresses and gas pressure. The early, rapid, pathway closure may be driven by the elastic (compressibility) of the material and the subsequent gradual response by the slow diffusion and egress of gas trapped along now partially sealed pathways. Fig. 10 shows the evolution of radial stress during gas injection, test stages 4 and 5. These sensors showed qualitatively similar behaviour to that previously discussed for axial stress. Considerable variations in stress values were observed across each plane, Fig. 11, ranging from 140 kPa to 2650 kPa. The data clearly shows the evolutionary and complex nature of the stress field induced during advective gas flow and is likely linked to the localised internal development of gas pathways within the clay. A plot of average stress for each plane of measurement, Fig. 12, shows that once major gas entry occurred, stress increased fairly consistently throughout the sample, with the exception of R13-R16. Here peak stress was substantially lower compared to all other arrays, suggesting fewer pathways propagated towards the top of the sample. In Fig. 10. Radial stress data for FPR-21–004. Radial stresses R13 through R16 register the lowest values as these sensors were located below the injection filter. Fig. 11. Difference between maximum and minimum stress at each axial and radial plane of measurement plotted as a function of time. Fig. 12. Average stress at each plane of measurement along the axis of sample FPR-21–004. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 9 contrast, during the major gas entry and breakthrough events, average stress values for A1-A4 were similar to those in array A5-A8, suggesting some pathways must have been horizontally orientated across the vertical axis of the sample. The highest stress values occurred in array R1- R4, closest to the backpressure filter, indicating that (i) the gas pressure gradient induced presumably small pathways and (ii) most of the pressure drop between the gas phase and the water within the backpressure filter occurred in a narrow zone close, or at, the edge of the sample. The evolution in pore pressure at filters F1 and F2, within sample FPR-21–004, is shown in Figure 13[A]. Following the cessation of hydration at day 27.2, a gradual decrease in pressure was observed, illustrating that the sample was not in full hydraulic equilibrium at the onset of gas testing. As hydration continued from the backpressure filter, a gradual increase in water pressure was seen from around day 38 and 49 in filters F1 and F2 respectively. Figure 13[B] shows the development of pore pressure during gas entry and subsequent breakthrough events. As major gas entry occurred at day 74.3, both filters showed an increase in pressure, probably related to a hydrodynamic effect during the displacement and compaction of the clay at gas entry. As injection gas breakthrough occurred at day 75.5, filter F2 showed a rapid increase in pressure, most likely due to the arrival of gas at the filter. Gas pressure peaked at 6.64 MPa and then began to gradually decrease. A second breakthrough event in the same filter, occurred at day 84 when pressure rapidly increased. Filter pressure then appeared correlated to the gas pressure with an offset of around 0.7 MPa by the end of stage 5. This data also demonstrates that the pressure drop along the connected gas pathways, from the injection point to the filter, was relatively small. In contrast, the development of pressure within filter F1 showed no abrupt changes in pressure as with filter F2. This suggests that filter F1 was likely recording local hydrodynamic changes in porewater pressure due to displacement and localised compaction of the clay as gas pathways formed. However, the increase in pressure of filter F1 from around day 87, as gas pressure declined, might be linked to the slow discharge of gas into the filter as the pathways began to drain. 3. Enriched multi-phase modelling approaches Three different numerical approaches have been adopted by BGR/ UFZ, LBNL and CIMNE-UPC/Andra (see also Appendix A for details about their key similarities and differences). These three strategies are based on the general theory of multi-phase flow modelling and are combined with additional features to describe some of the key aspects observed in gas-laboratory tests that distinguish clay-rich media from other rock-types (e.g., deformation of the porous media, creation of dilatant preferential pathways that open and eventually self-seal). Indeed, the bentonite sample is modelled as a deformable porous medium that behaves as either an elastic (LBNL, CIMNE-UPC/Andra) or an elasto-plastic (BGR/UFZ) solid. Intrinsic permeability is treated as a function of other properties during the gas injection test, representing the dilatant pathways. In particular, intrinsic permeability is assumed to depend either on the strain tensor (BGR/UFZ), on the effective minimum compressive stress (LBNL) or on the embedded fractures aperture and spacing (CIMNE-UPC/Andra). The water retention curve is also assumed to be a function of embedded fractures aperture and spacing in the model developed by CIMNE-UPC/Andra. 3.1. Model developed by BGR/UFZ 3.1.1. Conceptual model The model expands upon the work performed by BGR/UFZ within the Task A DECOVALEX-2019 (D-2019), where the team developed a hydro-mechanical model for migration of gas through a low-permeable linear elastic geomaterial that included a pressure-dependent intrinsic permeability. Indeed, for the current phase of DECOVALEX-2023, BGR/ UFZ developed a fully coupled, hydro-mechanical model based on multiphase flow theory, 23 whose key features are: •Constitutive relations for the hydraulic behaviour: the Mualem model is adopted for the description of the relative permeabilities of gas and water whereas the relationship between water saturation and capillary pressure is based on the van Genuchten formulation, 34 . •Constitutive relations for the mechanical behaviour: bentonite is assumed to behave as an elasto-plastic porous medium. Hooke’s law is assumed to describe the stress-strain relationship in the elastic regime of the material whereas the Drucker-Prager failure criterion is assumed to describe the plastic deformation (perfect plasticity with non-associated flow). This is enhanced with a tension cut-off parameter to limit the load carrying capacity of the model near the tensile region. •Hydro-mechanical coupling: Biot’s theory 35 is assumed to describe the hydro-mechanical coupling. That is, the effective stress tensor σ ′ (Pa) is calculated from the pore pressure p and the total stress tensor σ as σ ′ = σ − α p(1) (where α (-) is the Biot’s coefficient) and used to define the linear momentum balance equation of the porous medium ∇[ σ ′ − α (pg−Swpc)I]+ ρ g=0(2) where pg (Pa) is the gas pressure, Sw (-) is the water saturation, pc (Pa) is the capillary pressure, I is the identity tensor, ρ (kg/m 3 ) is the total density and g (m/s 2 ) is the gravitational acceleration. Material heterogeneity is included in the model by means of the gas entry pressure and Young’s modulus, which are described with a heterogeneous distribution. Both can be derived from a non-uniform distribution of the dry density, which has been observed in similar experiments. 36 The approach does not attempt to represent the scale of Fig. 13. [A] pore pressure development in filters F1 and F2 during test FPR-21–004. [B] expanded data from day 65 to day 77 showing the filter response during initial gas entry, major gas entry, initial gas breakthrough and major gas breakthrough respectively from left to right. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 16 Fig. 22. Average gas saturation profiles at each layer, at the base and at the top of the sample obtained with the three numerical models. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 17 stress values. However, the rapid transition phase seen at breakthrough is still too smoothly and too early predicted (Figs. 18 and 19). In terms of stresses, a similar behaviour is obtained when CIMNE-UPC/Andra’s approach - a coupled hydro-gas-mechanical 3D numerical model that assumes a heterogeneous initial permeability field and embedded fractures - is employed. In terms of outflows, an earlier (BGR/UFZ) and a later (CIMNE-UPC/Andra) peak are obtained. LBNL applied a continuum approach, with a direct relation between permeability and least principal effective stress, that models dilatant permeability increase once the local fluid pressure approaches the least principal stress, i.e. the least principal effective stress approaches zero. This approach along with a gas entry pressure that is scaled with the dilation-induced permeability can be used to simulate the macroscopic gas flow behaviour (Figs. 18 and 19), although the model is not fine enough to capture microscopic dilatancy and heterogeneous gas channelling. The observed shape of the post-breakthrough curve, which shows a shut-in behaviour and is directly linked to the stress state, is poorly reflected by all the models thus suggesting a better understanding or representation of the hydro-mechanical coupling is still required. New conceptual models may therefore need to be developed to explain this complexity (Appendix B). However, their development is still in a very preliminary phase and at the current stage of the work, they are not able to reproduce the experiments. Thus, from a pragmatic perspective, enhanced flow models remain the only viable tool to represent these systems at this time. In summary, this study illustrates further analyses are needed for a better understanding before models can be used, with confidence, as a predictive tool to assess advective gas movement: •Parameter calibration and model constraints: developed models need calibration of fitted parameters. This poses a computational challenge as a large number of simulations are required. In addition, their adjustment is a complex task: indeed, although it is possible to evaluate the main effect of a single parameter on the global response, several parameter correlations exist due to the complex model couplings. A great effort from the modelling teams has been made to tackle this issue and detailed sensitivity analysis have been carried out to characterise the one-dimensional gas injection test. However, their extrapolation to other tests can be difficult. This was already observed in Task A D-2019, 21 where different parameter values for the 1D and the 3D tests were sometimes arbitrarily prescribed. Thus, more numerical analyses and a better understanding of the complex couplings within the codes are needed before the models could be used with confidence as a predictive tool to assess gas movement. •Heterogeneity: the role of material heterogeneity needs to be further explored, as it has shown that it might provide a possible way to represent flow localisation. Indeed, some teams have assumed heterogeneous fields (e.g., permeability, gas entry pressure), as the introduction of heterogeneity at the capillary scale is not practicable. However, the assumed distribution functions are usually arbitrary prescribed and lack a physical justification. In the absence of data (e. g., characterisation of pore morphology, variation in material properties), this issue remains a challenge for all teams. •Stochasticity: experimental data exhibits a combination of stochastic and deterministic behaviours. Despite the growing awareness of the need to reflect non-deterministic process features, developed approaches in this study are purely deterministic and hence, it is therefore important for numerical modelling to consider how best to account for key stochastic features such as the number of gas pathways, the precise timing of the gas breakthrough and the associated gas flows. •Upscaling: although this paper presents a summary of the modelling work performed when simulating laboratory-scale experiments, models that are tractable at the engineered barrier and repository scales are needed. Due to the need to accurately and efficiently include small-scale heterogeneities in a field-scale model, this poses a major challenge and will constitute the main aim of the next stage of the D-2023 project. In this forthcoming stage, the modelling of a fullscale in situ test (Lasgit) will be undertaken. This is a highly instrumented test based on a mock canister, encapsulated in bentonite/ pellets, and placed within a deposition hole at the ¨ Asp¨ o underground research laboratory (Sweden). At the current stage of the modelling study, substantial future work is needed to account for some of the foreseen complexities (e.g., representation of the heterogeneous porous field-scale domain and the inclusion of the interfaces, since they play an important role). However, the expansion of the numerical approaches to a full-scale test will provide an invaluable insight to help to inform the modelling approaches and assess the impact of upscaling of gas flow on repository layout and therefore the design of any future facility. •Application: the validity of these models for the prediction of gas flow in repository systems remains limited. Upscaling the detailed processes represented in the current models would be computationally challenging, requiring compromises in both mesh density and simplification, or, averaging of less important processes. However, it is important to remember that models need only to be fit for their application, so long as end-users clearly understand their deficiencies. Any changes/simplifications to the models would require experimental validation where possible and be derived from a detailed understanding of the importance/relevance of key phenomenological features. 6. Conclusions Different numerical representations for the quantitative description of advective gas flow in clay-based repository systems have been developed. This study shows that models that were used to simulate particular controlled-laboratory experiments (Task A in the D-2019 project) can be employed to reproduce other tests, even with different sample geometries. This study has also shown that these numerical approaches can successfully represent some of the main deterministic experimental features typically observed in dilatancy-controlled gas flow (e.g., initial and peak stress values). However, there are other features that are not correctly captured (e.g., rapid transition phase seen at breakthrough). This suggests that numerical models still lack a correct description of the full complexity of the physical processes observed in water-saturated experiments. This study has identified three key learning points that need to be born in mind when numerically modelling gas flow through water saturated low permeable clay samples: 1. Model calibration: models need calibration of fitted parameters. Each model needs a specific calibration process and thus, their adjustment requires detailed sensitivity analysis. Since gas flow is semi-stochastic by nature, these analyses should distinguish between Fig. 24. Buoyant rising of a less viscous fluid within a more viscous fluid. Image is taken from [B.4]. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 18 the key experimental features reproducible across all experiments and those that only occur in specific experiments. 2. Heterogeneity: preferential pathways are not explicitly represented and hence, fracture-specific properties such as the evolving fracture front are not included into the models. Instead, heterogeneous material distributions can be used to implicitly represent preferential pathways. 3. Model coarsity: developed models are not fine enough to include microscopic dilatancy or heterogeneous gas channelling. Although this might represent a limitation, the use of fine meshes is not recommended, as approaches that can potentially be tractable at engineered barrier and repository scales are ultimately needed. CRediT authorship contribution statement E. Tamayo-Mas: supervision, validation, writing – original draft. J. F. Harrington: supervision, validation, data acquisition, writing – review and editing, funding acquisition. I. P. Damians, S. Olivella, E. Radeisen, J. Rutqvist, Y. Wang: model development, writing – review and editing. 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. Data Availability Approaches have been assessed against a new and hitherto unpublished gas test. A brief overview of the experimental data is given in the manuscript. Data might be made available on request. Acknowledgements DECOVALEX is an international research project comprising participants from industry, government and academia, focusing on development of understanding, models and codes in complex coupled problems in sub-surface geological and engineering applications; DECOVALEX- 2023 is the current phase of the project. The authors appreciate and thank the DECOVALEX-2023 Funding Organisations Andra, BASE, BGE, BGR, CAS, CNSC, COVRA, US DOE, ENRESA, ENSI, JAEA, KAERI, NWMO, NWS, SÚRAO, SSM and Taipower for their financial and technical support of the work described in this paper. The statements made in the paper are, however, solely those of the authors and do not necessarily reflect those of the Funding Organisations. Appendix A. model comparison Differences between the proposed numerical strategies lie in conceptual features, the software used by the teams, in the assumed geometry to represent the saturated bentonite, in the initial/boundary conditions prescribed for the tests and in the material parameters. Here, these differences are reported. Conceptual differences between the numerical strategies Table 5 Brief description of the three numerical models developed by the participating teams. BGR/UFZ LBNL CIMNE-UPC/Andra Hydro-mechanical (HM) coupling Fully coupled (via the Biot’s effective stress tensor) HM model Sequentially coupled (via the Biot’s effective stress tensor) HM model Fully coupled (via the Biot’s effective stress tensor) HM model Key hydraulic features Van Genuchten- Mualem model Van Genuchten- Corey model Fracture-dependent van Genuchten model +fracture-dependent relative permeabilities Mechanical deformation Elasto-plasticity (Drucker-Prager with a tension cut-off parameter) Linear elasticity (with swelling stress) Elasticity (with a dilatancy term in the deviatoric component of the volumetric strains) Dilatant pathways description Strain-dependent intrinsic permeability triggered by a heterogeneous Young’s modulus Permeability is assumed to depend on pressure and the effective minimum compressive stress Intrinsic and relative permeabilities are assumed to be decomposed into matrix and fracture terms Codes used by the teams Table 6 Software employed by the participating teams. BGR/UFZ LBNL CIMNE-UPC/Andra Software OpenGeoSys TOUGH2 + FLAC 3D CODE_BRIGHT Version 5.8 TOUGH2 V2.1 FLAC3D V5 8.6 Reference 38 45 + 46 48 E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 19 Test geometries used by the teams Table 7 Test geometries employed by the teams. BGR/UFZ LBNL CIMNE-UPC/Andra Discretisation method Finite element Integral finite difference Finite element Geometry 2D triangular axisymmetric mesh 3D quadrilateral axisymmetric mesh 3D hexahedral mesh Number of elements 2716 857 3430 Number of nodes 1448 1846 3909 Explicit description of the injector Via a source term Via a representative volume of the injection chamber With all the relevant components (i.e., injection rod and injection filter tip) Element order 2nd order (quadratic triangular) Linear Linear Prescribed initial conditions Table 8 Initial conditions prescribed by the teams. BGR/UFZ LBNL CIMNE-UPC/Andra Pore-water pressure (MPa) 1.00 1.00 0.25 Pore-gas pressure (MPa) 3.00 1.00 0.10 Water saturation 97% 100% 100% Swelling pressure Compressive stress: no real swelling 8.0 - Capillary pressure (MPa) 2.0 0.0 0.0 Displacements along x axis - 0.0 0.0 Displacement along y axis - 0.0 0.0 Displacements along z axis (if appropriate) - 0.0 0.0 Stress (MPa) 7.0 (effective compressive stress) 8.0 8.0 Temperature 20◦20◦20◦ Prescribed boundary conditions Table 9 Boundary conditions prescribed by the teams. BGR/UFZ LBNL CIMNE-UPC/Andra Pore-water pressure - 1 (at outlet) Variable (as per given test specifications) Pore-gas pressure (MPa) 1 (top of sample, outflow) 1 (at outlet) Variable (as per given test specifications) Flowrate injector ˙ m(t) = nV0 ΔVIn,t Prescribed until 85 days to match injection pressure rise Prescribed from 22 to 85 days as per test Flow rate at STP plot given Capillary pressure - 0 at outlet/inlet Variable (as per given test specifications) Displacements along x axis Right boundary =0 Left boundary =symmetry axis 0 normal to outer surface 0 at all outer surfaces Displacement along y axis Top and bottom boundary =0 0 normal to outer surface 0 at all outer surfaces Displacements along z axis (if appropriate) - 0 normal to outer surface 0 at all outer surfaces Temperature 20◦20◦20◦ Final parameter values employed by the teams The basic material parameters employed by the teams are listed and compared in Table 10. Specific model parameters that need to be calibrated (or that are considered known) are reported in following sections. Table 10 Basic parameters used by the teams (* indicates values beyond specified parameters). BGR/UFZ LBNL CIMNE-UPC/Andra Elastic modulus (MPa) Heterogeneous* 307 307 Poisson’s ratio (-) 0.4 0.4 0.4 (continued on next page) E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 20 Table 10 (continued) BGR/UFZ LBNL CIMNE-UPC/Andra Porosity (-) 0.43 * 0.44 0.4327 (initial porosity) Biot’s coefficient (-) 0.9 1 0.5 Intrinsic permeability of water (m 2 ) 3.0 ×10 −20 * 3.4 ×10 −21 Heterogeneous* Parameters employed by BGR/UFZ Table 11 Calibrated parameters employed by BGR/UFZ. Parameter Symbol Units Calibrated value Sensitivity analysis Effect with respect to the global response Minimum value Maximum value Plastic strain multiplier b1 - 4000 1000 4000 Impacts the intrinsic permeability Volumetric strain multiplier b3 - 250 50 250 Impacts the intrinsic permeability Intrinsic permeability kin m 2 7.0 ×10 −21 3.0 ×10 −21 6.0 ×10 −20 Impacts the breakthrough time and the total flow Mean of the Young’s modulus distribution E Pa 3.5 ×10 9 4.5 ×10 8 3.5 ×10 9 Impacts the breakthrough time Minimum value of the Young’s modulus distribution Emin Pa 5.0 ×10 8 - - Impacts the breakthrough time Table 12 Other parameters (assumed known) employed by BGR/UFZ to characterise the bentonite. Parameter Symbol Units Calibrated value Taken from Density of solid grains ρ s kg/m 3 1570 Decovalex 2019 21 Friction angle φ ◦32 - Dilatancy angle ψ ◦20 - Cohesion c MPa 3 - Tortuosity τ - 1 - Apparent gas entry pressure (van Genuchten model) P0 MPa 10 [A.1] Measure of the pore size distribution (van Genuchten model) n - 1.49 [A.2] Residual degree of saturation (van Genuchten model) Sres - 0.01 [A.1] Maximum degree of saturation (van Genuchten model) Smax - 0.99 [A.1] Parameters employed by LBNL Table 13 Calibrated parameters employed by LBNL. Parameter Symbol Units Calibrated value Sensitivity analysis Effect with respect to the global response Minimum value Maximum value Swelling coefficient β sw - 0.015 0.01 0.03 Impacts the total stress level Max aperture for stress-permeability bh0 m 4.9 ×10 −6 1.0 ×10 −7 1.0 ×10 −5 Impacts gas flow rate after break through Reference stress normal to the fracture σ n,ref Pa 0.1 0.05 1.0 Impacts gas flow rate after break through Residual saturation of gas (Corey model) Sgr - 0.13 0.05 0.2 Impacts gas entry (pressure) Multiplying factor for the enhanced gas permeability (Corey model) mg - 375 100 1000 Impacts gas flow rate Table 14 Other parameters (assumed fixed) employed by LBNL. Parameter Symbol Units Calibrated value Taken from Apparent gas entry pressure (van Genuchten model) P0 MPa 18 41 Shape factor (van Genuchten model) λ - 0.45 41 Residual liquid saturation (van Genuchten model) Slr - 0.01 41 E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 21 Parameters employed by CIMNE-UPC/Andra Table 15 Other parameters (assumed fixed) employed by CIMNE-UPC/Andra. Parameter Symbol Units Calibrated value Taken from Homogeneous part (same parameter value for the entire sample modelled) 21 Dilatancy angle ψ ◦24 Tortuosity for dissolved gas (Fick’s law) τ - 0.5 Transverse dispersion coefficient (Fick’s law) DT - 0.001 Longitudinal dispersion coefficient (Fick’s law) DL - 0.01 Power for liquid state for both matrix and fractures (Relative permeability) nmatrix/fractures liquid - 3 Maximum liquid saturation (Relative permeability) Smax liquid - 1 Minimum liquid saturation (Relative permeability) Smin liquid - 0 Power for gas state for matrix (Relative permeability) nmatrix gas - 2 Power for gas state for fractures (Relative permeability) nfractures gas - 1 Maximum gas saturation (Relative permeability) Smax gas - 0.3 Minimum gas saturation (Relative permeability) Smin gas - 0 Shape function (van Genuchten model λVG - 0.45 Initial porosity ϕ0 - 0.44 Dry density ρ dry kg/m 3 1512 Molar mass of Helium M kg/mol 0.004 Henry’s constant H MPa 10000 Heterogeneous part (three different sample zones randomly distributed). M1 (2/3) - M2 (1/6) - M3 (1/6), seeFig. 16. Initial capillary pressure (van Genuchten model) P0 MPa M1: 48.6 21 M2: 22.5 M3: 10.8 Finite air entry value P00 MPa M1: 5.4 M2: 2.5 M3: 1.2 Reference permeability (matrix intrinsic permeability) k0 m 2 M1: 1.0 ×10 −21 M2: 1.0 ×10 −20 M3: 1.0 ×10 −19 Internal associated width for each fracture a m M1: 5.0 ×10 −6 M2: 5.0 ×10 −5 M3: 5.0 ×10 −4 Initial aperture of the fractures b0 m M1: 1.5 ×10 −9 M2: 5.0 ×10 −9 M3: 9.5 ×10 −9 Maximum aperture of the fractures bmax m M1: 1.5 ×10 −7 M2: 3.5 ×10 −7 M3: 7.5 ×10 −7 Initial strain ε 0 % M1: 0.05 M2: 0.03 M3: 0.01 A. 7 Appendix references. [A.1] R. Senger, E. Romero, A. Ferrari and P. Marschall (2014). Characterization of gas flow through low-permeability claystone: laboratory experiments and two-phase flow analyses. Geological Society, London, Special Publications, 400, 1, pp. 531–543. DOI: 10.1144/SP400.15. [A.2] Z. Dai, J. Samper, A. Wolfsberg and D. Levitt (2008). Identification of relative conductivity models for water flow and solute transport in unsaturated bentonite. Physics and Chemistry of the Earth, Parts A/B/C, 33, 1, pp. S177-S185. DOI: 10.1016/j.pce.2008.10.012. Appendix B. new theoretical model developed by SNL SNL developed an alternative approach to describe the key features observed in gas-laboratory tests that distinguish clay-rich media from other rock-types. The new model expands the work performed in [B.1], where a phenomenological concept of nonlinear dynamics and deterministic chaos theory was employed to analyse the gas pressure and the gas influx and outflux obtained from a one-dimensional test performed on a pre-compacted Mx80 bentonite sample at the British Geological Survey. Indeed, the computation of a set of diagnostic parameters [e.g., global embedding dimension, correlation dimension, information dimension, spectrum of Lyapunov exponents (B.1. and references therein)] suggested that the prevailing processes during the gas test are (i) chaotic diffusion both at the onset of gas influx and during the final phase of the experiment and (ii) chaotic advection after the breakthrough. The time series analyses of the data indicate that the observed chaotic behaviour of the system can be described by 3–5 independent variables. Here, focus is placed on the development of a conceptual model for such chaotic systems. A detailed mathematical formulation and analysis will be provided elsewhere. The primary research focus of the SNL team within the current D-2023 project is to understand the actual mechanism(s) for the emergence of the observed complex behaviour of gas migration in water saturated compact bentonite. In particular: •Channelling postulate: although gas migration in a porous medium is generally treated as an immiscible displacement process within a “rigid” solid framework, this is unlikely to be the case for hydrated compacted bentonite, where the pore size is generally extremely small (on a scale of nanometers). Indeed, the capillary pressure pc in such a medium can be estimated by E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 22 pc=2 σ cos(θ) r(19) where σ (N/m) is the surface tension between gas and water, θ (◦) is the contact angle water on the solid; and r (m) is the radius of the pore necks. Typical values of these parameters ( σ =~70 mN/m, θ =~40◦, r =1 - 10 nm, see [B.2] and [B.3]) lead to a significantly higher pc (~10 - 100 MPa) than the gas pressure generally used in experiments and therefore, given the low tensile strength of bentonite, the only possible way for gas to move through such a medium is by channelling. •Channelling description: this channelling process is mathematically analogous to the buoyant rising of a viscous fluid within another viscous fluid, [B.4]. This allows the characterisation of both the spacing of channels λ=2 π L1 31/3( μ 2 μ 1)1/3 (20) and the perturbation growth of the interface ξ=2ΔPgL1 3(L1+L2) μ 2( μ 2 24 μ 1)1/3 (21) where μ 1 and μ 2 are the viscosities of gas and water-saturated compacted bentonite (Pa⋅s) respectively, ΔPg (Pa) is an externally imposed gas pressure gradient and L1, L2 (m) are the lengths defined in Fig. 23[A], characterising the compressibility of the bentonite, which is in turn related to the swelling pressure of bentonite as a function of the dry density of the material. As seen in Eq. (20), the spacing of channels only depends on the viscosity ratio μ 2/ μ 1 and the length L1. Hence, the number of channels to be developed in the system is determined by the viscosity ratio only, see Fig. 23[C]. Moreover, it is interesting to note that the spacing of channels is linearly proportional to the length L1, implying that the channeling mechanism is linearly scalable. In other words, the complex fluid flow behaviour in gas migration in bentonite observed in small-scale experiments can probably occur in a large-scale experiment in the field, as highlighted by the multiple gas tests experiments performed within the Lasgit project. •Gas breakthrough postulate: SNL postulates that the breakthrough point of gas migration corresponds to the time when a fully percolating channel network is established across the material domain. One important constraint on the establishment of such networks is that the gas pressure must be high enough to make enough room for gas percolation by compressing the bentonite matrix. The room available for gas migration is directly determined by the swelling capacity of the material and the dry density used in the testing, and thus the breakthrough pressure is around the swelling pressure. Since this room is usually a small fraction of the total material volume in a test, the apparent gas saturation degree is expected to be small. •Gas movement postulate (after the channel network is created): as experimentally and theoretically shown in [B.4], if the upwelling material has greater viscosity than the surrounding material, the structure is a long vertical column with gradually decreasing diameter and, if the upwelling material has less viscosity than the surrounding material, the structure envelops a rim syncline and a pronounced overhang and eventually ascends as a spherical pocket of fluid fed by a pipe, see Fig. 24. Analogously, and due to the large contrast in viscosity between the bentonite matrix and the gas injected, SNL postulates that gas percolates through a water saturated bentonite matrix in the form of individual gas bubbles. The size of bubbles is determined by the apparent surface tension of the gas-bentonite interface. •Gas movement description (after the channel network is created): previous work carried out by SNL in [B.5] uses a logistic map model to describe the movement of a gas bubble in a deformable clay matrix which can display a complex dynamic behaviour that can be explained with a bifurcation and chaos concept. Indeed, the dynamic behaviour of the system has been shown closely related to clay matrix dilation, fracturing and fracture healing as induced by gas bubble movement. Therefore, these components should be included in any model describing the advective movement of gas through low-permeability materials. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 23 Fig. 23. Stability analysis of gas-water saturated bentonite interface: [A] modelling system, [B] perturbation growth of the interface and [C] spacing of channels determined by the viscosity ratio. . B.1. Appendix references [B.1] B. Faybishenko, Y. Wang, J.F. Harrington, E. Tamayo-Mas, J. Birkholzer and C. Jov´ e -Col´ on (2022). Phenomenological Model of Nonlinear Dynamics and Deterministic Chaotic Gas Migration in Bentonite: Experimental Evidence and Diagnostic Parameters. Transport in Porous Media, 141, p. 585–606. DOI: 10.1007/s11242–021-01733–9. [B.2] Y. Wang, C. Bryan, T. Dewers, J.E. Heath and C. Jove-Colon (2012). Ganglion Dynamics and Its Implications to Geologic Carbon Dioxide Storage. Environmental Science and Technology, 47, 1, p. 219–226. DOI: 10.1021/es301208k. [B.3] Y. Wang (2014). Nanogeochemistry: Nanostructures, emergent properties and their control on geochemical reactions and mass transfers. Chemical Geology, 378–379, pp. 1–23. DOI: 10.1016/j.chemgeo.2014.04.007. [B.4] J.A. Whitehead Jr and D.S. Luther (1975). Dynamics of laboratory diapir and plume models. Journal of Geophysical Research, 80, 5, pp. 705–717. DOI: 10.1029/JB080i005p00705. [B.5] Y. Wang, T. Hadgu, E. Kalinina, C. Jove-Colon and B. Faybishenko (2018). International Collaboration on Spent Fuel Disposition in Crystalline Media: FY18 Progress Report. DOI: 10.2172/1469442. References 1. Sellin P, Leupin OX. The use of clay as an engineered barrier in radioactive-waste management–a review. Clays Clay Miner. 2013;61(6):477–498. https://doi.org/ 10.1346/CCMN.2013.0610601. 2. Shirazi SM, Kazama H, Salman FA, Othman F, Akib S. Permeability and swelling characteristics of bentonite. Int J Phys Sci. 2010;5(11):1647–1659. 3. Daniels KA, Harrington JF, Zihms SG, Wiseall AC. Bentonite permeability at elevated temperature. Geosciences. 2017;7(1):3. https://doi.org/10.3390/ geosciences7010003. 4. Harrington JF, Volckaert G, Noy DJ. Long-term impact of temperature on the hydraulic permeability of bentonite. Lond: Geol Soc Lond. 2014;400:589–601. https://doi.org/10.1144/SP400.31. 5. Dixon DA, Stone J, Birch K, Kim CS. Sealing materials for a deep geological repository: Evaluation of swelling pressure and hydraulic conductivity data for bentonite-based sealing materials proposed for use in placement rooms. Can Geotech J. 2023;60(7):951–965. https://doi.org/10.1139/cgj-2021-0609. 6. Li KP, Chen YG, Ye WM, Wang Q. Self-sealing behavior of bentonite-based materials in high-level radioactive waste disposal: A systematic review, 106873 Appl Clay Sci. 2023;235(106873):2023. https://doi.org/10.1016/j.clay.2023.106873. 7. Harrington JF, Daniels KA, Wiseall AC, Sellin P. Bentonite homogenisation during the closure of void spaces. Int J Rock Mech Min Sci. 2020;136, 104535. https://doi. org/10.1016/j.ijrmms.2020.104569. 8. International Atomic Energy Agency (1997). Experience in Selection and Characterization of Sites for Geological Disposal of Radioactive Waste, IAEA TECDOC-991, IAEA, Vienna. 9. W.R. Rodwell, A.W. Harris, S.T. Horseman, P. Lalieux, W. Müller, L. Ortiz Amaya, K. Pruess (1999). Gas Migration and Two-phase Flow through Engineered and Geological Barriers for a Deep Repository for Radioactive Waste: a Joint EC/NEA Status Report (Nuclear Science and Technology). EC Report EUR19122. European Commission, Luxembourg. 10. Horseman ST, Harrington JF, Sellin P. Gas migration in MX80 buffer bentonite. MRS Online Proc Libr. 1996;465:1003–1010. https://doi.org/10.1557/PROC-465-1003. 11. J.F. Harrington, S.T. Horseman (2003). Gas migration in KBS-3 buffer bentonite. Sensitivity of test parameters to experimental boundary conditions. Technical Report TR-03–02. Svensk K¨ arbr¨ anslehantering AB (SKB), Stockholm, Sweden. 12. Harrington JF, Milodowski AE, Graham CC, Rushton JC, Cuss RJ. Evidence for gasinduced pathways in clay using a nanoparticle injection technique. Mineral Mag. 2012;76(8):3327–3336. https://doi.org/10.1180/minmag.2012.076.8.45. E. Tamayo-Mas et al. Geomechanics for Energy and the Environment 37 (2024) 100528 24 13. Cuss RJ, Harrington JF, Giot R, Auvray C. Experimental observations of mechanical dilation at the onset of gas flow in Callovo-Oxfordian claystone. Lond: Geol Soc Lond. 2014;400:507–519. https://doi.org/10.1144/SP400.26. 14. Cuss RJ, Harrington JF, Noy DJ, Graham CC, Sellin P. Evidence of localised gas propagation pathways in a field-scale bentonite engineered barrier system; results from three gas injection tests in the Large scale gas injection test (Lasgit). Appl Clay Sci. 2014;102:81–92. https://doi.org/10.1016/j.clay.2014.10.014. 15. Harrington JF, Graham CC, Cuss RJ, Norris S. Gas network development in a precompacted bentonite experiment: Evidence of generation and evolution. Appl Clay Sci. 2017;147:80–89. https://doi.org/10.1016/j.clay.2017.07.005. 16. Graham CC, Harrington JF, Sellin P. Gas migration in pre-compacted bentonite under elevated pore-water pressure conditions. Appl Clay Sci. 2016;132:353–365. https://doi.org/10.1016/j.clay.2016.06.029. 17. Harrington JF, Graham CC, Cuss RJ, Norris S. Gas network development in compact bentonite: key controls on the stability of flow pathways. Geofluids. 2019;2019:1–19. https://doi.org/10.1155/2019/3815095. 18. Xu L, Ye WM, Ye B, Chen B, Chen YG, Cui YJ. Investigation on gas migration in saturated materials with low permeability. Eng Geol. 2015;197:94–102. https://doi. org/10.1016/j.enggeo.2015.08.019. 19. Cui LY, Ye WM, Wang Q, Chen YG, Chen B. Gas migration behavior in saturated bentonite under flexible conditions with consideration of temperature effects. Acta Geotech. 2023;18(2):971–984. https://doi.org/10.1007/s11440-022-01624-3. 20. Senger R, Romero E, Marschall P. Modeling of gas migration through lowpermeability clay rock using information on pressure and deformation from fast air injection tests. Transp Porous Media. 2018;123:563–579. https://doi.org/10.1007/ s11242-017-0962-5. 21. Tamayo-Mas E, Harrington JF, Brüning T, et al. Modelling advective gas flow in compact bentonite: lessons learnt from different numerical approaches. Int J Rock Mech Min Sci. 2021;139, 104580. https://doi.org/10.1016/j.ijrmms.2020.104580. 22. Guo G, Fall M. Modelling of dilatancy-controlled gas flow in saturated bentonite with double porosity and double effective stress concepts. Eng Geol. 2018;243: 253–271. https://doi.org/10.1016/j.enggeo.2018.07.002. 23. Radeisen E, Shao H, Hesser J, Kolditz O, Xu W, Wang W. Simulation of dilatancycontrolled gas migration processes in saturated bentonite using a coupled multiphase flow and elastoplastic H2M model. J Rock Mech Geotech Eng. 2023;15(4): 803–813. https://doi.org/10.1016/j.jrmge.2022.05.011. 24. Chittenden N, Benbow S, Bond A, Norris S. Development of an upscalable HM model for representing advective gas migration through saturated bentonite. Int J Rock Mech Min Sci. 2020;133, 104415. https://doi.org/10.1016/j.ijrmms.2020.104415. 25. Yang J, Fall M. Coupled hydro-mechanical modelling of dilatancy controlled gas flow and gas induced fracturing in saturated claystone. Int J Rock Mech Min Sci. 2021;138, 104584. https://doi.org/10.1016/j.ijrmms.2020.104584. 26. Damians IP, Olivella S, Gens A. Modelling gas flow in clay materials incorporating material heterogeneity and embedded fractures. Int J Rock Mech Min Sci. 2020;136, 104524. https://doi.org/10.1016/j.ijrmms.2020.104524. 27. K.A. Daniels, J.F. Harrington (2017). The response of compact bentonite during a 1D gas flow test. British Geological Survey Open Report. OR/17/067. 28. T. Manai (1995). EVEGAS: European validation exercise of GAS migration models through geological media: Phase 1 - Final Report. European Commission Report EUR 16639 EN. European Commission, Nuclear Science and Technology Series, Directorate-General, Science, Research and Development, Luxembourg. 29. T. Manai (1997). EVEGAS: European validation exercise of gas migration models through geological media: Phase 2 - Final Report. European Commission Report EUR 17556 EN. European Commission, Nuclear Science and Technology Series, Directorate- General, Science, Research and Development, Luxembourg. 30. T. Manai (1997). EVEGAS: European validation exercise of gas migration models through geological media: Phase 3 - Final Report. European Commission Report EUR 17557 EN. European Commission, Nuclear Science and Technology Series Directorate- General, Science, Research and Development, Luxembourg. 31. P.J. Nash, B.T. Swift, M. Goodfield, W.R. Rodwell (1998). Modelling Gas Migration in Compacted Bentonite: A report produced for the GAMBIT Club. POSIVA Report 98–08. Posiva Oy, Helsinki. 32. B.T. Swift, A.R. Hoch, W.R. Rodwell (2001). Modelling Gas Migration in Compacted Bentonite: GAMBIT Club Phase 2 Final Report. POSIVA Report 2001–02. Posiva Oy, Olkiluoto. 33. A.R. Hoch, K.A. Cliffe, B.T. Swift, W.R. Rodwell (2004). Modelling Gas Migration in Compacted Bentonite: GAMBIT Club Phase 3 Final Report. POSIVA Report 2004–2. Posiva Oy, Olkiluoto. 34. van Genuchten MT. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J. 1980;44(5):892–898. https:// doi.org/10.2136/sssaj1980.03615995004400050002x. 35. Biot MA, Willis DG. The elastic coefficients of the theory of consolidation. J Appl Mech. 1957;24(4):594–601. https://doi.org/10.1115/1.4011606. 36. Villar MV, Iglesias RJ, García-Si˜ neriz JL. State of the in situ Febex test (GTS, Switzerland) after 18 years: a heterogeneous bentonite barrier. Environ Geotech. 2018;7(2):147–159. https://doi.org/10.1680/jenge.17.00093. 37. W.J. Xu, H. Shao, H. Hesser, W.Q. Wang, O. Kolditz, T. Popp (2011). Simulation of dilatancy-controlled gas migration process in saturated argillaceous rock. Proceedings of the 2nd International Symposium on Computational Geomechanics (Comgeo II), Cavtat-Dubrovnik. 38. Kolditz O, Bauer S, Bilke L, et al. OpenGeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media. Environ Earth Sci. 2012;67(2):589–599. https://doi.org/10.1007/ s12665-012-1546-x. 39. Gray WG, Hassanizadeh SM. Unsaturated flow theory including interfacial phenomena. Water Resour Res. 1991;27(8):1855–1863. https://doi.org/10.1029/ 91WR01260. 40. Kunhwi K, Rutqvist J, Harrington JF, Tamayo-Mas E, Birkholzer JT. Discrete dilatant pathway modeling of gas migration through compacted bentonite clay. Int J Rock Mech Min Sci. 2021;137, 104569. https://doi.org/10.1016/j.ijrmms.2020.104569. 41. R.K. Senger, P. Marschall (2008). Task Force on EBS/Gas Transport in Buffer Material. NAGRA NAB 08–24, Wettingen, Switzerland. 42. Rutqvist J, Ijiri Y, Yamamoto H. Implementation of the Barcelona Basic Model into TOUGH-FLAC for simulations of the geomechanical behavior of unsaturated soils. Comput Geosci. 2011;37:751–762. https://doi.org/10.1016/j.cageo.2010.10.011. 43. Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. Int J Rock Mech Min Sci Geomech Abstr. 1983;20(6):249–268. https://doi.org/10.1016/ 0148-9062(83)90595-8. 44. Rutqvist J. An overview of TOUGH-based geomechanics models. Comput Geosci. 2017;108:56–63. https://doi.org/10.1016/j.cageo.2016.09.007. 45. Pruess K, Oldenburg C, Moridis G. TOUGH2 user’s guide, version 2.1, LBNL-43134 (revised). Berkeley, USA: Lawrence Berkeley Natl Lab. 2012. 46. Itasca (2018). FLAC3D VERSION 6.0. Fast Lagrangian Analysis of Continua in 3 Dimensions. ITASCA Consulting Group Inc. 47. Olivella S, Alonso EE. Gas flow through clay barriers. G´ eotechnique. 2008;58(3): 157–176. https://doi.org/10.1680/geot.2008.58.3.157. 48. Olivella S, Gens A, Carrera J, Alonso EE. Numerical formulation for simulator (CODE_BRIGHT) for coupled analysis of saline media. Eng Comput. 1996;13(7): 87–112. https://doi.org/10.1108/02644409610151575. 49. I.P. Damians, B.S. Noghretab, S. Olivella, A. Gens (2022). Numerical Mesh Sensitivity Works Performed to Analyze Advective Gas Flow in a Compact Clay. Paper presented at the 56th U.S. Rock Mechanics/Geomechanics Symposium, Santa Fe, New Mexico, USA, June 2022. 〈doi:10.56952/ARMA-2022–0679〉. E. Tamayo-Mas et al.