Phenomenological Model of Nonlinear Dynamics and Deterministic Chaotic Gas Migration in Bentonite: Experimental Evidence and Diagnostic Parameters

Understanding gas migration in compacted clay materials, e.g., bentonite and claystone, is important for the design and performance assessment of an engineered barrier system of a radioactive waste repository system, as well as many practical applications. Existing field and laboratory data on gas migration processes in low-permeability clay materials demonstrate the complexity of flow and transport processes, including various types of instabilities, caused by nonlinear dynamics of coupled processes of liquid–gas exchange, dilation, fracturing, fracture healing, etc., which cannot be described by classical models of fluid dynamics in porous media. We here show that the complexity of gas migration processes can be explained using a phenomenological concept of nonlinear dynamics and deterministic chaos theory. To do so, we analyzed gas pressure and gas influx (i.e., input) and outflux (i.e., output), recorded during the gas injection experiment in the compact Mx80-D bentonite sample, and calculated a set of the diagnostic parameters of nonlinear dynamics and chaos, such a global embedding dimension, a correlation dimension, an information dimension, and a spectrum of Lyapunov exponents, as well as plotted 2D and 3D pseudo-phase-space strange attractors, based on the univariate influx and outflux time series data. These results indicate the presence of phenomena of low-dimensional deterministic chaotic behavior of gas migration in bentonite. In particular, during the onset of gas influx in the bentonite core, before the breakthrough, the development of gas flow pathways is characterized by the process of chaotic gas diffusion. After the breakthrough, with inlet-to-outlet movement of gas, the prevailing process is chaotic advection. During the final phase of the experiment, with no influx to the sample, the relaxation pattern of gas outflux is resumed back to a process of chaotic diffusion. The types of data analysis and a proposed phenomenological model can be used to establish the basic principles of experimental data-gathering, modeling predictions, and a research design.


Introduction
Investigations of gas migration in compacted clay materials, e.g., bentonite and claystone, are important for the design and performance assessment of nuclear waste disposal facilities and engineered barrier systems of a radioactive waste repository system. Special attention in these investigations is given to understanding and predicting the mechanisms controlling gas entry and flow through bentonite, which may impact the repository's barrier performance and host rock. The sources of gas (such as hydrogen) are anoxic corrosion of metallic structural components (Ortiz et al. 2002), radioactive decay of the disposed radioactive waste, radiolysis and the microbial breakdown (Pekala et al. 2019) of organic materials. The production of these gases may continue for over 100,000 years after the waste emplacement (Norris 2014).
There are three principal clay-based repository projects ongoing in Europe, using the Opalinus clay rock (Early Jurassic) in Switzerland, the Callovian-Oxfordian clay rocks of the Middle to Late Jurassic in France, and the Boom clay rock (Early Oligocene) in Belgium (Grambow 2016). Over the past 30 years, a series of projects have been performed to understand gas migration-MEGAS (Volckaert et al. 1995;Ortiz et al. 1997Ortiz et al. , 1998, PROGRESS (Rodwell 1999(Rodwell , 2000, and FORGE (Rueedi et al., 2011). The purpose of Task A (modeling Gas Injection ExpERiments-ENGINEER) of the International Project DECOVALEX-2019 (https:// decov alex. org/D-2023/ task-a. html) is to better understand the processes governing the movement of gas in two types of low-permeability materials: engineered compacted Mx80 bentonite, and a natural repository host clay rock. As part of these projects, various experimental field and laboratory studies, along with numerical simulations, have been conducted to investigate the mechanisms and controlling factors, such as gas entry and flow, pathway stability and sealing, which are expected to impact barrier performance and clay-based host rock behavior. It was found, for example, that if the rate of gas production exceeds the rate of gas diffusion within the barrier or host rock, discrete gas migration pathways are expected to form (Weetjens and Sillen 2006;Ortiz et al. 2002;Wikramaratna et al. 1993). Under these conditions, gas will continue to accumulate until its pressure becomes sufficiently large to generate a breakthrough to enter the surrounding material.
Understanding the complex mechanisms of gas migration are fundamentally important for the development, testing, and verification of numerical models for the quantitative simulations and prediction of gas fluxes, which could be used to support future repository design, and would be important in the assessment of gas pressure evolution as part of safety case development. However, existing field and laboratory data on gas migration processes in low-permeability clay materials demonstrate the complexity of flow and transport processes, including various types of instabilities, caused by nonlinear dynamics of coupled processes of liquid-gas exchange, dilation, fracturing, fracture healing, etc., which cannot be described by classical models of fluid dynamics in porous media. For example, Berkowitz et al. (2006) demonstrated that even in a relatively homogeneous porous medium, the dispersive flux cannot be expressed as a Fickian-type law. Bolshov et al. (2007) reviewed field and laboratory observations and basic physical concepts of nonclassical transport processes in geologic media. Cushman and O'Malley (2015) stated that "…the occurrence of Fickian dispersion in geophysical settings is a rare event and consequently should be labeled as anomalous." Misztal et al. (2015) simulated anomalous dispersion in porous media using the unstructured lattice Boltzmann method. Ruspini et al. (2014) reported the results of experiments demonstrating two-phase flow instabilities. Ortiz et al. (2002) conducted studies of the potential effects of the mechanical deformation of clay on radionuclide transport, a long-term in situ gas migration experiment through an artificially created pathway in the Boom Clay. Podgorney et al. (2000) and Faybishenko et al. (2001) described the results of field infiltration tests in fractured basalt, using the deterministic chaotic analysis. Faybishenko (2002Faybishenko ( , 2004Faybishenko ( , 2005 introduced a concept of chaotic liquid-gas flow and mixing in fractured rock based on the results of field infiltration experiments and laboratory experiments conducted by Persoff and Pruess (1995). Faybishenko (2017) identified dynamic causal inference involved in flow and transport processes in complex fractured-porous media. Yang et al. (2020) described a stochastic gas-liquid two-phase flow seepage model based on numerical simulations of a horizontal Hele-Shaw experiment.
The results of a series of laboratory experiments in Mx80 bentonite Graham et al. 2012Graham et al. , 2015Harrington et al. 2019) and Callovo-Oxfordian claystone (COx) (Cuss et al. 2011;Harrington et al. 2017a, b) showed that the advective transport of repository gas in bentonite causes the development of new porosity and propagation of dilatant pathways, resulting in the generation and evolution of a gas network. It was also found a strong coupling between stress, gas pressure and flow in bentonite, which takes place during multiple discrete propagation events, demonstrating spatial variability and timedependency of bentonite permeability. Harrington et al. (2019) concluded that characterization of the gas network distribution is of fundamental importance in predicting gas dissipation rates and understanding the long-term fate of gas in radioactive waste repositories.
The objective of this paper is to demonstrate the presence of phenomena of nonlinear dynamics and chaos of gas migration through the bentonite, based on the results of computations of diagnostic parameters of nonlinear dynamics and deterministic chaos, using the gas influx and outflux time series data obtained from the gas injection experiment in the Mx80 bentonite.
The paper is structured as follows. Sect. 2 includes a general description of the gas injection experiment in the Mx80 bentonite core, which was conducted at the British Geological Survey; Sect. 3 describes methods of the evaluation of diagnostic parameters of nonlinear dynamics and deterministic chaos of gas influx and outflux time series, using the time-domain and phase-space analyses; Sect. 4 presents the results of calculations of these parameters; Sect. 5 includes a phenomenological model of nonlinear dynamics and deterministic chaos of gas transport in bentonite; and Sect. 6 contains conclusions.
The types of QA/QC analysis of raw time series datasets and a decomposition of time series to remove high-frequency random data are presented in the Appendix.

Experimental Studies
Controlled gas injection experiments have been performed by the British Geological Survey on pre-compacted samples of Mx80 buffer bentonite using the injection of helium as a safe replacement for hydrogen gas. Basic physical properties of a bentonite sample are summarized in Table 1. A schematic diagram in Fig. 1 illustrates the design of the constant volume and radial flow gas migration apparatus, used by the British Geological survey. The schematic shows filters and total stress sensor locations. Each radial flow array comprises four filters set at 90° to each other; five total stress sensors are located around the periphery of the sample, two axial and three radial sensors. The schematic also shows the location of the central gas injection filter and radial drainage arrays. The pressure and flow to each radial array (comprised of four separate filters) was controlled by a separate syringe pump. Sample dimensions are 120 mm in length and 60 mm in diameter. Technical details of the testing assembly can be found in Harrington et al. (2017a, b). Once the samples have been installed and the end-closures attached, a backpressure was applied to fully water-saturate each bentonite sample, resulting in bentonite hydration due to the uptake of water and generation of a swelling pressure. The tests were conducted under ambient temperature (20°C), and inlet and outlet gas fluxes and pressure were monitored. Time series graphs of injection and measured gas pressure are shown in Fig. 2. Nonlinear dynamics calculations were performed for the periods before the breakthrough (segments A and B) and after the breakthrough (periods E and H).
The sample was initially hydrated by applying a constant pore-water pressure of 1 MPa (applied through the large filters embedded in each end-closure and those of the radial Time series of injection and measured gas pressures recorded during the Mx80 gas injection experiment conducted by the British Geological Survey. Vertical dashed lines indicate the test stages A through H, which are identified based on the gas injection pressure and pore pressure arrays) and allowed to re-equilibrate. It was observed that both axial and radial stresses increased because of sample swelling. The results of observations on the Mx80 core analyzed in this paper are typical for many other experiments conducted by the British Geological Survey.

Methods of Calculations of Diagnostic Parameters of Nonlinear Dynamics and Chaos
The time lag (time delay) is the time between the discrete time series points when the relationship between the points almost vanishes (Abarbanel 1996). The time lag is determined in the time domain using the method of averaged mutual information (AMI). The AMI function provides the information learned about one observation from another one on the average over all measurements (Fraser and Swinney 1986). The AMI is a measure of information (predictability) that x(t) can provide about x(t + τ). The first minimum of the AMI function vs. the time lag is used to determine an optimal time lag. Global embedding dimension (D GED ) is the dimension of the phase space needed to unfold the attractor of a nonlinear system from the observed scalar signals in such a way that the trajectories of the D GED -dimensional attractor do not cross each other. The global embedding dimension is determined using a method called false nearest neighbors, FNN, which is used to determine whether a nearest neighbor of a certain point in phase space is a true neighbor characterizing a dynamic system or a false projection caused by a too low dimensional phase space (Kennel et al. 1992;Abarbanel 1996;Abarbanel and Tsimring 1998). A minimum value of FNN function corresponds to the optimal D GED .
Local embedding dimension (D L ) characterizes how the dynamic system evolves on a local scale, and indicates the number of degrees of freedom governing the system dynamics, i.e., how many dimensions should be used to predict the system dynamics (Abarbanel and Tsimrin 1998). It also identifies the number of the Lyapunov exponents to be used in evaluating the system. The local embedding dimension, D L , is less or equals to the global embedding dimension, D GED . In the current paper, we adopt that D L = 3, and then use the D L values for calculations of the spectrum of Lyapunov exponents.
Correlation dimension, D cor , is a scaling exponent characterizing a cloud of points in an n-dimensional phase space and is calculated from (Grassberger and Procaccia 1983).
or as the slope of log(C(r)) versus log r: where C(r) is the number of pairs separated by distances less than r.
Information dimension, D inf , is a measure of the fractal dimension of a probability distribution of the time series. The information dimension is a particular case of the generalized correlation dimension when setting the order q = 1. The information dimension D inf may be defined as: D inf = lim{r→0} < ln p(r) > /ln(r), where p(r) is the probability of finding a neighbor in a neighborhood of size r, and < > indicates the mean value. The information dimension is computed for different embedding dimensions for checking if D inf saturates (as with calculations of D cor ).
Analysis of a spectrum of Lyapunov exponents is one of the most valuable diagnostic criteria of a deterministic chaotic system. The Lyapunov exponents are used to evaluate the sensitivity of a nonlinear dynamical system to small changes in initial conditions, and, in other words, the stability and an asymptotic behavior of nonstationary solutions of the numerical model describing the system's behavior. The spectrum of Lyapunov exponents is determined by measuring the mean exponential growth (positive Lyapunov exponents) or shrinking (negative Lyapunov exponents) of perturbations with respect to a nominal trajectory of the attractor pathway.
The number of Lyapunov exponents is equal to the Local Embedding Dimension (D L ), which is either less or equal to the D GED . Each Lyapunov exponent gives the rate of expansion or contraction along the coordinates of the attractor's phase space. For a chaotic system, the largest Lyapunov exponent must be positive, which implies that a line segment in the phase space grows exponentially with time -exp(λ i T), and a negative sum of the Lyapunov exponents (Σλ i < 0) indicates that the phase-space attractor does not expand indefinitely, but occupies a limited phase-space volume. A random non-deterministic chaotic system is characterized by all non-positive Lyapunov exponents. Generally, a Lyapunov exponent indicates how two nearby points in a phase space of a dissipative nonlinear dynamical system move exponentially apart in time. The exponentially rapid separation of initially close points leads to sensitivity of the system to initial conditions.
An attractor of a dynamical system is a set of numerical values toward which a system tends to evolve from initial (starting) conditions of the system. The dimension of an attractor is one of the most important parameters necessary to characterize the system's properties. Generally, the attractor dimension is equal to the value of D GED , which indicates a number of essential variables needed to model the system's dynamics (Farmer et al. 1983;Abarbanel 1996). A dynamical system with a chaotic attractor can be locally unstable, but globally stable, i.e., it behaves within an n-dimensional basin of attraction (Grebogi et al. 1987). An attractor is generally a space used to represent all possible states of a system, within which each possible state corresponds to one unique point in the phase space. Phase-space attractors are the graphs of the relationship between the physical parameters of the system. However, dynamical system attractors can be presented graphically only as 2D or 3D graphical objects. A concept of a pseudo-phase attractor is used to characterize the behavior of a dynamical system based on a single (i.e., univariate) time series. 2D or 3D pseudo-phase-space attractors are plotted using one-dimensional vectors with points separated by a time lag: for a 2D attractor-as scalar vectors X i (t) and X i (t + Δt), and for a 3D attractor-as X i (t), X i (t + Δt), and X i (t + 2Δt). One of the main features of the attractor of a deterministic chaotic system is that the size of the attractor remains practically the same even if initial conditions are slightly disturbed. However, the geometry of chaotic attractors can be complex and difficult to describe (Ott 2002). Such attractors are also called strange attractors, and they have a fractal structure.

Results of Time Series Analysis
The plots of time variations of the injection pressure and pressure measured at different locations in the core are shown in Fig. 2. The time series was subdivided into the segments according to the injection pressure pattern shown in Fig. 2. Segments A and B were selected to analyze the time series of influx oscillations, and Segments E and H were selected to analyze the outflux oscillations.  Fig. 3b shows the plots of the cumulative influx and outflux during the experiment. These plots reveal that the breakthrough of gas through the outlet occurred after 64.5 days of gas injection, during Segment D, when the injection pressure reached about 10,000 kPa. The total volume of the injected gas was 810 ml, while the total outflux was only 260 ml, indicating the overall accumulation of gas in the tested bentonite core.
To assess whether the flux fluctuations were caused by temporal variations of the injection pressure, Figs. A2 (Segment A) and A3 (Segment B) depict the graphs of the influx and expanded graphs of the injection pressure, and Fig. A4 shows the box-and-whisker plots of the injection pressure. The injection pressure varied in a very narrow range (from 990 to 1027 kPa for Segment A, and from 2994 to 3023 kPa for Segment B), which could have been caused partially by minor thermal noise. Bottom plots of Figs. A2 and A3 show no specific relationships between the influx and the injection pressure, indicating that the flux fluctuations are unlikely caused by an external driving force. Statistical analysis of the spread of fluctuations of the pressure indicates that 95.9% of data points were within the two standard deviations around the mean for Segment A, and 98.6%-for Segment B. Figure A5 shows the results of a decomposition of the original time series datasets into a trend (middle graphs in Fig. A5) and a random component (lower graphs in Fig. A5) for Segments A, B, E, and H. Figure 4 shows a comparison of the original (gray color) and the trend (i.e., noise-reduced) time series (red color). The trend of time series data was then used for calculations of diagnostic parameters of nonlinear dynamics and chaos. Figure 5 shows the results of calculations of the time lags for Segments A, B, E, and H, which are Δt = 7, 5, 9, and 7, correspondingly. The calculated time lags were then used for calculations of other nonlinear dynamics parameters (called calculations in the phase space)-global embedding dimension, correlation and information dimensions, a spectrum of Lyapunov exponents, and plotting 2D and 3D pseudo-phase attractors of influx and outflux. Fig. 6 shows that the FNN curves drop to practically zero at the values of the embedding dimension D GED = 5, 4, 4, and 4, correspondingly for Segments A, B, E, and H. The range of D GED from 4 to 5 indicates phenomena of low-dimensional deterministic chaos. Note that for a purely deterministic chaotic system, such as the Lorenz model, the percentage of false nearest points drops to zero at D GED = 3 and remains zero for higher dimensions, implying that the attractor is fully unfolded at D GED = 3 and remains unfolded at greater embedding dimensions.
Calculated values of the correlation dimension D cor vary from 3.21 to 3.77, and the values of D inf are from 3.26 to 4.14, which are summarized in Table 2. The results shown in   Table 2 indicate that D inf values are greater than the D cor (except for Segment H), which is typical for a deterministic chaotic system (Tsonis 1992). Figure 7 shows the calculated Lyapunov exponents of the time series for Segments A, B, E, and H. For all segments, calculations were performed using D L = 3. (Note that taking D L = D GED will not change eventually the conclusions.) For all segments, there is one positive Lyapunov exponent (red line in Fig. 7), one exponent close to zero (green line in Fig. 7), one negative Lyapunov exponent (blue line in Fig. 7), and the sum of Lyapunov exponents is non-positive (Table 2). This spectrum of Lyapunov exponents is typical for the deterministic chaotic behavior, indicating the converging attractor. Figures 8 and 9 present 2D and 3D pseudo-phase-space attractors for gas influx and outflux. For Segments A and B, the basins of attraction are close to the beginning of coordinates, attracting nearby oscillations. For Segments E and H, the attractors reflect the pattern of bounded oscillations, taking place without loss of stability, and which are called hidden oscillations. These hidden oscillations are filling in a compact subset of the phase space, attracting nearby oscillations, and which are called hidden attractors.

Phenomenological Model of Gas Flow in Bentonite
A phenomenological model can be defined as a theory that expresses the behavior of observed phenomena without paying detailed attention to their fundamental significance (Thewlis 1973). A phenomenological model can also be defined as a scientific model that is not derived from first principles and describes the empirical relationship of phenomena, which is consistent with fundamental theory, but is derived from observational data. For example, a phenomenological description of gas transport processes in Opalinus Clay was given by Marshall et al. (2005), based on the microstructural model concept and basic transport mechanisms.
According to the results summarized in Horseman et al. (1999) and Harrington and Horseman (2003), the transient patterns of influx and outflux and pressure measurements provide clear evidence of stimulation of the microfracture network, discrete gas pathway 'sealing' events, and localization of gas flow within multiple pathways resulting in a nonuniform gas discharge rate pattern at the outlet. Also, all tests exhibited a threshold pressure for breakthrough and showed a post-peak negative transient leading to steady-state gas flow. At the beginning of the test, with no pressure-induced pathways, the water-saturated bentonite was impermeable to gas. Harrington and Horseman (2003) hypothesized that gas inlet flux and gas breakthrough were affected by the development of pathways which propagated through the clay from gas inlet to outlet, and the flow exhibited nonlinear behavior, although two-phase flow was impossible.
In the current paper, based on the results of time series analysis of gas influx and outflux we extended a conceptual model of the gas transport given by Marshall et al. (2005) and Cuss et al. (2014). Fig. 10 (upper row) demonstrates phenomenological considerations of the basic gas transport mechanisms. The onset of gas influx in the bentonite core, which is developed before the gas breakthrough through the core outlet (Segments A, B, and C), is characterized by the process of chaotic diffusion, which lasted 63.5 days. We postulate that, in chaotic diffusion, the injected gas attempts to make its way through the water saturated bentonite with frequent gas bulb formation and annealing, leading to high-frequency influx fluctuations. After the breakthrough, when gas reached the backpressure end of the sample, with end-to-end movement of gas, the prevailing process is chaotic advection taking place along microscopic and macroscopic tensile fractures (Segment E). The chaotic flux behavior at this stage is probably caused by the instability of gas permeation pathways. The green lines on the schematics of the gas transport in Fig. 10 indicate schematically the randomly generated gas flow paths within porous media; the yellow dashed lines and red arrows indicate that the overall volume of affected porous media is expanding prior to and postbreakthrough, and then the volume shrinks during the recovery. Black lines for the period of the chaotic advection indicate dilatant inlet-to-outlet pathways, which could intersect, leading to mixing phenomena (as it is observed for chaotic advective liquid flow in fractured rock Faybishenko 2005). The 3D pseudo-phase-space attractors corresponding to the initial chaotic diffusive transport prior to the breakthrough are shown in Figs. 9a,b, and corresponding to the advective chaotic transport-in Fig. 9c. A 3D pseudo-phase-space attractor for the final relaxation/recovery stage of gas outflux, with no influx to the bentonite sample is shown in Fig. 9d.

Conclusions
Understanding gas migration in compacted clay materials (e.g., bentonite or claystone) is important for a performance assessment of an engineered barrier system of a repository system. Existing field and laboratory experimental data on gas migration processes in lowpermeability clay materials demonstrate the complexity of flow and transport processes, such as various instabilities, caused by nonlinear dynamics of liquid-gas exchange, dilation, fracturing, fracture healing, etc., which, however, cannot be described by classical models of fluid dynamics. In the current paper, we applied methods of nonlinear dynamics and chaos to the analysis of the results of observations of gas pressure and gas influx and outflux, recorded during the gas injection experiment in the Mx80 bentonite sample, and calculated a set of the diagnostic parameters of nonlinear dynamics and chaos. At the onset of gas influx in the bentonite core, before the breakthrough, the development of dilatant pathways (Harrington et al. 2013) can be characterized as the process of chaotic gas diffusion. After the breakthrough, with inlet-to-outlet movement of gas, the prevailing process is chaotic advection. During the final phase of the experiment, with no influx to the sample, the Mx80 core recovery is again characterized by a pattern of gas chaotic diffusion. The types of data analysis and a proposed phenomenological model, not derived from first principles, can be used to establish the basic principles of a phenomenological experimental Fig. 9 3D pseudo-phase-space strange attractors plotted using time lags determined by means of the AMI approach data-gathering and a modeling research design. Calculations of diagnostic parameters of nonlinear dynamics and deterministic chaos provide important constraints on physical and numerical models development, which will be a subject of future investigations. The proposed concept provides a new perspective for modeling gas migration in low-permeability materials and is of direct relevance to other clay-based scientific and engineering issues associated with advective gas flow, such as shale gas, hydrocarbon migration, carbon capture and storage, gas storage and landfill design.  Cuss et al. 2014), which are representative for different stages of gas transport through bentonite. The green lines on the schematics of the gas transport indicate the randomly generated gas flow paths within porous media; the yellow dashed lines and red arrows indicate that the overall volume of affected porous media is expanding prior to and post-breakthrough, and then the volume shrinks during the recovery. Black lines for the period of the chaotic advection indicate dilatant pathways originating at the inlet and ending up at the outlet The experimentally measured time series records usually contain a certain noisy, i.e., high-frequency random, component caused by uncertainty or technical errors of measurements, along with the trend and a periodic component. Removing high-frequency noise, i.e., the random component, by means of a low-pass filter is commonly used for more accurate analysis of nonlinear dynamical processes. The random component of the original time series was determined using a function 'decompose" of the R library "stats," and then it was extracted from the original dataset. (The minimal frequency value of f=2 was used to convert the original timeseries datasets into the "ts" object to perform calculations with the "decompose" function. Note that the use of a minimal value of frequency allowed for removing the minimal amount of noise.) Diagnostic parameters of chaos were calculated using R libraries "fractal," fNonlinear," "nonlinearTseries," "tseriesChaos," and 3D attractors were plotted using the library "scatterplot3d."    The results of the decomposition of the measured records of influx and outflux, showing the original data, time trend, and a random component that was extracted to determine the general trend to be used for data analysis