Modelling a hydration and heating column test on unsaturated bentonite using a double-porosity approach

. A controlled re-saturation of a heated 50 cm long column made of MX-80 bentonite pellets in the CIEMAT laboratory facilities (Madrid, Spain) has been running for 10 years (2011-2021). The experimental data on the saturation state and the swelling response of the bentonite material at high temperatures provided useful information for the verification of a double-porosity approach that considers the hydro-mechanical coupling between the micro-and macro-pore levels and the definition of retention curves for each structural domain. The numerical modelling with this double-porosity model has shown to be able to reproduce the main behaviour features observed in the lab test.


Introduction
The thermo-hydro-mechanical (THM) behaviour of bentonite-based materials have been extensively studied in the last decades due to their potential use as sealing or buffers in geological repositories for high-level radioactive wastes. In such environments, the engineered barrier system (EBS) is usually subjected to: i. changes in temperature due to the heat released from the radioactive decay and its propagation through the barrier; ii. changes in the moisture content due to the water uptake from the rock and the vapour generation and diffusion from the hotter towards the cooler zones; iii. changes in stresses due to the physical confinement of the bentonite barrier under non-isothermal hydration conditions.
In order to understand and quantify the main THM processes that take place in such an environment and to minimize the uncertainties during the transient phase, a hydration and heating test was carried out in the CIEMAT laboratory (Madrid, Spain). Such a laboratory test, performed in the framework of the Long-Term Performance of Engineered Barrier Systems (PEBS) Project [1], is an important contribution to the verification of mathematical formulations capable of reproducing the main constitutive response of expansive clays submitted to changes in their moisture content under isothermal and non-isothermal conditions.
The occurrence of irreversible changes in the fabric of active clays due to the physico-chemical processes taking place at particle level [2] restrict the use of singleporosity formulations in many geotechnical problems involving multi-structure porous media. In this respect, the development of mechanical constitutive models by [3] and [4] considering the existence of two dominant and overlapping pore structural levels and their * Corresponding author: ramon.barboza@upc.edu application to the modelling of bentonite buffers constituted a crucial enhancement of previous constitutive formulations for unsaturated clayey soils. In the present work, the modelling task is performed using an improved double-porosity model (DPM) that will be briefly described in the following sections.

Infiltration test at high temperatures
The hydration of a heated 50 cm-long column of pelletized MX-80 bentonite (termed as "Cell B" in [5]) (2011-2021) was aimed to the completion of the experimental database of the full-scale HE-E test in order to identify the THM processes occurring inside and around the EBS. Therefore, the granular sealing material in Cell B was subjected to heating and hydration from different ends simulating the same boundary conditions as in the Nagra Section of the in situ HE-E test [1].
The heating of Cell B was performed by means of a plane stainless steel heater placed at the bottom of the cell. Heater temperature was increased in two steps. In the first heating step, temperature was increased instantaneously from a value around 23°C (lab temperature) to 100°C. During the second heating step, the heater temperature was raised from 100°C to the target value of 140°C. Both heating episodes occurred without hydration of Cell B. An improvement of the external insulation took place between these two heating steps with the goal of reducing the lateral heat losses through the body of the cell.
The hydration phase started 5015 hours after the beginning of heating. A cooling chamber located at the top of the cell allowed the circulation and the injection of saline water (at lab temperature) into Cell B. A constant water pressure of 0.06 bar was applied until the end of the experiment.
The evolution of temperature and relative humidity (RH) inside Cell B during the test was monitored by thermo-couple sensors installed at 10 cm, 22 cm and 40 cm from the heater. Furthermore, a ring load cell was emplaced at the top of Cell B in order to record the axial pressure developed during the hydration of the mass of pellets under high thermal gradient. The schematic layout of the instrumented cell is shown in Fig. 1(a).
The dismantling of Cell B was carried out 10 years after the beginning of the test. Several samples were taken from the dismantled column of pellets in order to determine the water content and dry density distribution along the column. Detailed information about the experimental procedures carried out prior and during the construction of Cell B, the main observations and the data collected during the test as well as the dismantling operations and the final state of the bentonite can be found in [5]- [7].

Bentonite material characterization
The MX-80 granulate inside Cell B is the same bentonite material emplaced in the Nagra Section of the HE-E experiment. The initial water content and the average dry density of the bentonite column were 6.4% (RH≈40%) and 1.53 g/cm 3 , respectively, while the dry density of the solid grains is about 2.75 g/cm 3 [5]. Therefore, the initial average saturation and porosity of the column were 22% and 0.44, respectively. It was also reported that 69% of pores are larger than 200 nm [7], which could be considered as the initial macro-pore fraction for the double-porosity approach used in the numerical modelling of this test (see Section 3).
The thermal conductivity of the mass of pellets at the initial state conditions of the in situ test and at laboratory temperature varies between 0.12 W/m/K [5] and 0.45 W/m/K [8]. The hydraulic conductivity, estimated from the percolation of saline water through samples of powdered MX-80 bentonite compacted at different densities, varies between 7.5e-13 m/s and 2.5e-13 m/s for dry densities between 1.40 g/cm 3 and 1.60 g/cm 3 , respectively [9], [10]. Swelling pressure tests on small samples of powdered and pelletized MX-80 bentonite compacted at that dry density range yielded values ranging from 2 MPa to 8 MPa [11], [12].

Geometry and boundary conditions
The infiltration test at high temperature was modelled by means of the finite element (FE) code CODE_BRIGHT [13] as a 2-D axisymmetric geometry. Therefore, the bentonite column is simulated as a 3.5 cm radius and 50 cm high geometry. The lower and the upper horizontal boundaries of the column represent the heating element and the hydration surface, respectively. The body of Cell B (made of Teflon PTFE) and the layers of dense foam and wool that act as a thermal barrier are also included in the domain of analysis. The change of geometry between the two heating steps due to the improvement of the insulation system is also considered in the modelling of the test. The modelled geometries before and after the enhancement of the insulation system are shown in Fig. 1(b). All dimensions are taken in accordance to the test dimensions specified in [5].

Theoretical formulation of the DPM
In the current double-structure approach, the expansive soil is described by two distinct but overlapping structural domains: a micro-pore domain, composed by the internal structure of individual particles and aggregates of active clays and their vicinity and a macrostructural level corresponding to the arrangement of clay aggregates and other soil grains with macro-pores among them [2]. Hence, the micro-pore fraction includes both the intra-and inter-particle voids within a clay aggregate (micro-and meso-pore domains).
A three-phase diagram representing the soil element as a double-porosity medium is shown in Fig. 2. It can be noted that the volume of voids in the soil ( ) is the sum of the volume occupied by macro-pores ( ) and micro-pores ( ). The volume of microstructure ( ) includes the volume of micro-pores and the volume of the solid phase ( ). The total porosity for the active clay ( ) and the porosity for the clay aggregate ( ) are defined as To develop the mathematical formulation of the DPM, two additional porosity variables are introduced in Eq. (1) in order to describe the natural decomposition of total porosity as the sum of the pore-volume fractions in microstructure ( ̅ ) and macrostructure ( ̅ ). The saturation state in the porous domain α (« m » for micro and « M » for macro) by the γ phase (« l » for liquid and « g » for gas), , is defined as follows where denotes the fraction of the volume of voids in the α-structural level ( ) filled with the γ phase.

Additive decomposition of strain
The explicit consideration of two pore-structure levels in the mathematical formulation and the definition of micro-and macro-pore volume fractions in Eq. (1) lead to the additive partition of the strain tensor (̇) in a micro (̅̇m) and a macro (̅̇M) component. The strain tensor can be also split into elastic and plastic components. In addition, it is reasonable to assume that micro-structural strains are purely volumetric and fully reversible [2]. The elastic behaviour of expansive clays is described by the elastic components at micro (̅̇e ) and macro (̅̇M e ) levels, while the irreversible structural changes are given by the plastic strain component of macro-structure (̅̇M p ). Therefore, ̇= ̅̇m + ̅̇M = ̅̇m e + ̅̇M e + ̅̇M p (4) The micro-structural strain component of the active clay is related to the volumetric expansion or contraction of the clay particle/aggregate, ̇, according to This expression plays a crucial role in the mathematical development of the DPM since it allows to correlate the local microstructural behaviour with the global deformation response in terms of the fraction of the soil (1 + ) occupied by the microstructure (1 + ).

Main THM constitutive equations
The THM constitutive equations relate the independent variables (displacements, fluid pressures, temperature, etc.) with the dependent ones (stress tensor, degree of saturation, fluxes of liquid and gas phases, heat fluxes, etc.). A brief description of the general expressions for the THM constitutive laws used in the current doubleporosity formulation is presented below. A detailed explanation of such equations can be found in [14].
Heat conduction through the soil ( c ) is driven by temperature gradients ( ) according to Fourier's law c = − • (6) where the global thermal conductivity ( ) is a function of the thermal conductivity of each phase in the soil and the average degree of saturation ( ), defined as = ̅ + ̅ (7) The expression in Eq. (7) evidences the dependence of the thermal conductivity on the water content and the pore volume fraction in each pore-structure level.
The advective fluxes are formulated for the macrostructural domain, since the water and the gas mobility through micro-pores are neglected. The dependence of the water permeability on the pore structure is given by a modified Kozeny-Carman equation while its dependence on the effective saturation of macro-pores ( , ) is given by an exponential function , = ( , ) In such equations, M is the intrinsic permeability tensor for the current macrostructure; ̅ ,0 is the pore fraction for which the reference intrinsic permeability tensor ( M,0 ) is estimated; , is the liquid phase relative permeability; and , are model parameters. The assumption of non-saturation of micro-pores leads to the description of hydraulic constitutive laws for the micro-pore domain as well. Therefore, it is necessary to define a water retention curve (WRC) for each structural level ( = for micro; = for macro) as where is the matric suction; , are model parameters; is a parameter related to the pore-air entry value dependent on temperature ; and is the suction at fully dry conditions. The occurrence of unbalanced water potentials in micro-and macro-pores leads to a local exchange of water between both structural media ( ) proportional to the difference in micro suction ( ) and macro suction ( ) = ( − ) (11) where the leakage parameter, , is related to geometric characteristics of the porous media [15]. If the gas mobility between the two pore levels is considered high enough to assume local balance of gas pressure, the mass transfer is only driven by the difference between the liquid pressures in each pore domain.
The generalized stress-strain relationship for the αstructural domain can be expressed as ̇α = [ α ]̇α + lα̇+ gα̇+ Tα̇ (12) where ̇α is the engineering stress vector; ̇α is the strain vector; [ α ] is the constitutive stiffness matrix; lα , gα are constitutive vectors relating changes in liquid and gas pressures, respectively, to stress increments; and Tα is the vector relating stress to temperature changes.
In the present double-porosity approach, the elastic strain components at micro (̅̇m e ) and macro ( (14) where ̅ m e , ̅ M e are the elastic constitutive matrices for micro and macrostructure, respectively; ̇ is the net effective stress for the clayey soil; ̅̇m →M e is the elastic strain component due to the structural coupling between both pore domains; ̅ is the volumetric elastic modulus for the microstructure; , ̅ are the thermal bulk modulus of the porous medium and the bulk elastic modulus related to changes in macro suction, ̇, respectively; ℋ is a non-linear function of temperature, saturation and fluid pressures in micro-pores; and = [1 1 1 0 0 0] is an auxiliary vector. The mechanical coupling of the DPM formulation is achieved by the following expression relating the elastic tensors and by the assumption of equivalence between the elastic constitutive tensor for the soil, [ e ] −1 , and for the clay aggregate, [ m e ] −1 : (15) Irreversible changes in the clay fabric are attributed to the loading-collapse (LC) mechanism and to the micro-macro coupling mechanism (β-mechanism) [3], [4], [14]. The plastic LC mechanism of the DPM formulation is defined by the same set of equations and parameters that describe the single-porosity, thermoelastoplastic Barcelona Basic Model (BBM). However, in the current DPM the evolution of the yield surface in the stress space is dependent on the stress state, history variables and macro suction. The yield function for a general state, , in terms of the stress invariants ( , , ), is defined as (16) where is the slope of the the critical state line; is a function of Lode's angle ( ); expresses the dependence of the shear strength on macro suction and temperature; and 0 is the apparent non-saturated preconsolidation pressure that defines the LC locus in the : plane according to where 0 * is the pre-consolidation pressure at saturation; is a reference pressure; (0), ( ) are the macrostructural compressibility at saturated and non-saturated conditions, respectively; and ̅ is the elastic stiffness parameter at macro for changes in the mean stress, .
The generation of macro-structural plastic deformations due to the LC mechanism is given by the plastic flow rule ̅̇M ,LC p =̇ (18) where the magnitude of the plastic strain (̅̇M ,LC p ) is determined by the plastic multiplier, ̇. The plastic potential function, , is defined as ) 2 ( + )( 0 − ) = 0 (19) in which is the non-associativity parameter. The plastic β-mechanism is related to macro-structural irreversible changes induced by the volumetric deformation of microstructure. For the sake of simplicity, the magnitude of this component of plastic strain (̅̇M ,β p ) is taken proportional to the micro-structural elastic strains at isothermal conditions, that is, where defines a pair of coupling functions dependent on the degree of compactness of macro-structure ( ) and on the sign of the mean effective stress increment at the micro-structural domain (̇) , where ̅ , , is the volumetric plastic strain due to the LC mechanism and ̅ , , is the volumetric component of the plastic strain related to the coupling mechanism.

Model parameters
The thermal parameters required in the numerical modelling were fitted in accordance to the experimental data in the literature [5], [8], [16] and by backcalculating the results of the in situ HE-E test [17]. The WRC for the macro and micro media were calibrated from data provided by oedometer tests on small samples of granular MX-80 bentonite submitted to controlled wetting or drying paths at different temperatures [18] and covering a wide range of densities (see [19]). The hydraulic parameters for the water flux through the bentonite material were estimated from back-calculating the in situ data from RH sensors installed inside the clay barriers of the HE-E test (see [17]). The leakage coefficient for the local water mass transfer between the two families of pores was estimated from the evolution of the axial pressure at the top of Cell B together with the experimental data for the saturation state inside the bentonite column after its dismantling. The mechanical parameters for the bentonite pellets were taken from experimental data on the mechanical characterization of pelletized and granular bentonite mixtures with analogous composition and/or with similar grain and dry density ranges. However, some of those model parameters had to be calibrated during the modelling of the HE-E experiment [17] in order to fit the model predictions to the data recorded by the sensors installed in the in situ and in the lab heating tests. The main input parameters used to simulate the artificial hydration of Cell B are listed in Table 1.

Numerical results and discussion
The ability of the fully coupled DPM to reproduce the evolution of temperature and RH in the early stages of the test can be observed in Fig. 3. Distances are taken from the bottom of Cell B. The good agreement between the numerical and the measured data is more evident at the sensors closer to the heater. Temperature oscillations recorded by the sensors were associated with the changes in the lab temperature [10] and were not reproduced by the model.
The changes in temperature are related to the increments in temperature caused by the heating process itself and by the improvement of the external insulation system. In addition, the high thermal gradient induced by the heating element at the bottom and the cooling system at the top of the column [5] was also reproduced in the modelling. The thermal gradient was steeper in the 10 cm closest to the heater, which reflected the saturation state along the height of the column. It can also be noted that temperature recorded by the sensor located at 40 cm from the heater remained practically unchanged and equal to the average room temperature.
The two heating episodes generated an increase in the RH at the sensor closer to the heater due to the differential expansion of macro-pore water in comparison to the thermal expansion of the solid species and the soil skeleton. The magnitude of the thermal load and the improvement of the insulation system induced a marked drying (a reduction of RH) of the bentonite material near the heater. The condensation of water vapour coming from the bottom (by diffusion) and the   downward flux of liquid water during the hydration phase led to a progressive increase of saturation (an increase of RH) in the cooler zones of the column.
The saturation state along Cell B at the end of the infiltration test is plotted in Fig. 4. It can be noted that the upper portion of the bentonite column is almost fully saturated by the water supplied during hydration. However, the model overestimated the saturation of the middle portion of the column. This could be a consequence of the higher thermal gradient observed in the lab test together with the reduction in porosity due to the compression exerted by the upper layers of the column as the material swells under restricted volume conditions. It is important to note that the evolution of the saturation of macro-pores is governed by the water permeability while the evolution of the water content in micro-pores is due to the local water mass transfer controlled by the leakage parameter. The evolution of the axial pressure at the top of Cell B is mostly due to the expansion of the bentonite pellets during the hydration phase, as shown in Fig. 5. Although the model predictions for the swelling behaviour in the early stages of hydration fit the laboratory data quite well, the model was unable to reproduce the magnitude of the axial stress measured in the laboratory at the end of the test. This difference could be a consequence of the lateral friction between the bentonite column and the body of the cell, which is not considered in the calculation. It is important to mention that the performance of constitutive analyses with the same model parameters used in the modelling of Cell B leads to swelling pressures of around 4 MPa. This value is in agreement with the experimental data obtained for this type of bentonite material with similar dry density values (see [11], [12]).
Despite some noted discrepancies between the experimental data and the model results, it can be concluded that the current double-porosity approach is able to satisfactorily simulate the main behaviour features observed during the infiltration of a heated column of MX-80 pellets. Additional calibration of some model parameters that control the evolution of the retention curves and the clay fabric may be required to achieve a closer representation of the observations.