THC Modelling of Bentonite Barrier of Geological Repository in Granite and Its Impact on Long-Term Safety

Narkuniene, Asta; Grigaliuniene, Dalia; Poskas, Gintautas

doi:10.3390/app14177851

Open AccessArticle

THC Modelling of Bentonite Barrier of Geological Repository in Granite and Its Impact on Long-Term Safety

by

Asta Narkuniene

^*

,

Dalia Grigaliuniene

and

Gintautas Poskas

Nuclear Engineering Laboratory, Lithuanian Energy Institute, 3 Breslaujos Str., 44403 Kaunas, Lithuania

^*

Author to whom correspondence should be addressed.

Appl. Sci. 2024, 14(17), 7851; https://doi.org/10.3390/app14177851

Submission received: 7 August 2024 / Revised: 28 August 2024 / Accepted: 28 August 2024 / Published: 4 September 2024

Download

Browse Figures

Versions Notes

Abstract

:

As in any other industry, nuclear energy results in the accumulation of some waste, which needs to be managed safely and responsibly due to its radiotoxicity. In the case of highly radioactive waste, geological disposal in stable rock is considered a broadly accepted solution. For the evaluation of the long-term safety of a geological repository, the assessment of radionuclide transport needs to be carried out. Radionuclide transport through engineered and natural barriers of the repository will highly depend on the barriers’ transport-related properties, which will be determined by coupled thermal, hydraulic, chemical, mechanical, biological, and radiation processes taking place in those barriers. In this study, the thermo-hydro-chemical (THC) state of bentonite was analysed considering CO₂ gas diffusion and temperature-dependent solubility in water. Reactive transport modelling of bentonite under non-isothermal conditions was performed with the COMSOL Multiphysics software (v6.0), coupled with the geochemical solver Phreeqc via the iCP interface. The modelling demonstrated that the consideration of chemical processes in bentonite had no significant influence on non-reactive Cl⁻ transport; however, it would be important for other radionuclides whose sorption in porous media depends on the porewater pH. Based on the modelling results, changes in the bentonite mineralogical composition and, subsequently, porosity depend on the partial CO₂ pressure at the bentonite–granite boundary. In the case of low CO₂ partial pressure at the bentonite–granite interface, the calcite dissolution led to a slight porosity increase, while higher CO₂ partial pressure led to decreased porosity near the interface.

Keywords:

geological disposal; coupled thermo-hydro-chemical processes; modelling; COMSOL Multiphysics; iCP; Phreeqc; reactive transport

1. Introduction

In any industry, including nuclear energy, huge benefits inevitably bring some less beneficial residual effects. These might be in the form of emissions (like the coal industry) or in the form of solid materials (microplastics, toxic waste, radioactive waste, etc.). As a society thinking about the future, we need to take care of those materials in a responsible way. In the case of highly radioactive waste, it is internationally demonstrated that reprocessing and transmutation of specific radionuclides do not eliminate the need for permanent disposal underground but would help reduce its footprint. Therefore, a number of disposal concepts have been developed, and much research has been initiated over the last decades. Geological disposal (repository) is realized by emplacing high-level radioactive waste in a geological environment surrounded by carefully selected and designed engineered barriers. Such a system of barriers compatible with each other (multibarrier system) ensures the containment and isolation of radionuclide-containing material over a hundred thousand years. In a carefully selected stable geological environment, the characteristics of the host rock, aquifers, aquitards up to biosphere, and prevailing conditions (stresses, groundwater pressures, temperature, etc.) are determined by processes and conditions billions of years earlier, and the engineered barrier properties are those that could be adjusted to be compatible with the surrounding geological environment. Thus, the processes taking place in engineered barriers and their state at different times in the geological repository lifetime is crucial.

Bentonite material is among other materials considered for the construction of engineered barriers in geological repositories. Due to its favourable characteristics, such as low hydraulic conductivity, diffusion dominated transport, high sorption capacity, etc., it is often foreseen as having a use in spent nuclear fuel or high-level waste disposal in granite in such countries as Sweden, Finland, Canada, and Japan [1]. Under geological disposal conditions, a number of processes will take place in the bentonite barrier such as heat transfer, groundwater flow, vapour transport, transport of diluted species, chemical interaction with accessory minerals in the bentonite and bentonite structure itself, etc.

From the perspective of long-term safety of the repository, it is important to assess the possibilities and extent of radionuclide release and transport from the waste matrix. Released radionuclides could dissolve in porewater and be transported further with it or could be transported in gaseous phase. In low permeability porous media, radionuclides are transported by diffusion mainly. Therefore, an effective diffusivity is a key parameter describing this transport mechanism. It is well known that diffusivity in a porous medium is lower than in free water due to pore structure. The latter could be affected by a number of geochemical processes, which, in turn, depend on thermal, hydraulic, mechanical, etc., conditions. Apart from the porosity change, geochemical processes are important from the perspective of sorption of radionuclides (pH [2,3], redox conditions, availability of complexing agents, etc.).

In order to evaluate the long-term safety of a disposal facility, a reliable prediction should be made about the state of the repository barriers (engineered and natural). The complexity of the repository and disposal conditions requires a multidisciplinary approach. In addition to a number of in situ heating experiments in the underground research laboratories [4] with the focus on the thermo-hydraulic (TH), thermo-hydro-mechanical (THM) state of the host rock (such as ALC experiment in the Callovo-Oxfordian claystone [5], FE experiment in the Opalinus clay [6] or the PRACLAY test in the Boom clay [7]), there are also experiments focused on the bentonite barrier. The FEBEX (full-scale engineered barrier experiment in crystalline host rock) in situ test is among these. The test had been carried out for 18 years and has been dismantled, providing a large set of chemical data that can be used for studying chemical state in the FEBEX bentonite. The other bentonites are being investigated as well. For example, the Wyoming-type bentonite and the Czech bentonite of type BCV in the HotBENT test started in 2021 (Grimsel Test Site, Switzeland), but the data from dismantling will be available after 5 years [8]. In line with experimental studies and knowledge from various disciplines, numerical models are being developed for complex analysis and long-term predictions. Numerical models differ as they were developed with varying states of interest, i.e., some for the bentonite hydration state and some for dry density distribution [9,10]. Other studies focused on the reactive transport of dissolved species under non-isothermal conditions [11,12,13,14], taking into account the effects of microstructure evolution during hydration and dehydration [15]. Zheng et al. [11,14] developed a THMC model for the FEBEX in situ test, and their numerical analysis showed that water content and the concentrations of dissolved species were strongly dependent on the changes in intrinsic permeability, thermal osmotic permeability, and initial dissolved concentration. Samper et al. [13] revisited the THMC model developed by Zheng et al. [11] by updating boundary conditions near the heater, refining the finite element grid, and increasing the vapour tortuosity factor and analysed species concentrations sensitivity to the bentonite retention curve, smectite dissolution, vapour tortuosity, and diffusion coefficients. However, these studies did not consider the impact of reactive gases such as CO₂ [11,14] on the bentonite performance. The groundwater coming into contact with the bentonite barrier contains a number of dissolved chemical species and dissolved CO₂, which will affect calcite dissolution/precipitation in the bentonite. Thus, it is important to assess the CO₂ potential transport. The CO₂ transport in the gas phase was considered in [13,16], but the study did not take into account the CO₂ degassing due to elevated temperature. CO_2g released near the heater could then be transported in the gas phase and redissolve in the condensation zone, and it could affect the pH in the bentonite. The impact of separate processes under isothermal geological repository conditions—such as variable flow conditions, CO₂ gas diffusion, mineral dissolution/precipitation, and cation exchange—was analysed in [17], where six cases of increasing complexity were benchmarked with different numerical simulators. The sensitivity of pH and CO₂ profiles to chemical activity calculation and mesh discretisation were reported. In the study, different codes were used and reasonable agreement between results was obtained.

This study considers non isothermal conditions and additional phenomena (temperature-dependent CO₂ solubility in water, surface complexation) for the evaluation of the thermo-hydro-chemical (THC) state of the bentonite. The reactive transport modelling of bentonite under non-isothermal conditions was performed with COMSOL Multiphysics, coupled with the geochemical solver Phreeqc via the iCP interface [18]. This required performance of the modelling of the two-phase flow of a miscible fluid (water and air) considering important phenomena, such as heat transfer; gas dissolution; advective–diffusive transport in the aqueous and gaseous phase, including CO₂ transport; and chemical interactions (aqueous complexation, mineral dissolution/precipitation, ion exchange, and surface complexation). In order to evaluate the effect of the CO₂ gas diffusion on the bentonite barrier performance, additional cases of different complexity were simulated, and the results between the cases were compared. In addition, this paper presents some insights about possible CO₂ degassing and flux of different magnitude imposed by surrounding rock.

2. Materials and Methods

The importance of the thermo-hydraulic-chemical state of engineered barriers of a geological repository and of understanding key phenomena responsible for it is acknowledged by the international community. The system to be analysed was inspired by the FEBEX in situ test [19] (Figure 1).

The FEBEX in situ test was installed and carried out in a tunnel excavated in granite in the Grimsel underground research Laboratory (Switzerland). There were a couple of heaters (diameter of 0.9 m) surrounded by a bentonite barrier of 0.69 m thickness. The tunnel diameter was 2.28 m. The test was designed to keep the temperature of 100 °C on the heater–bentonite interface. Due to heating, the temperature gradient was established across the bentonite layer. In such a system, a number of coupled processes took place: groundwater flow into the bentonite barrier; vapour generation and transport by advection and diffusion; advective and diffusive gas flow; heat transfer in gas, liquid, and solid phases; reactive advective–dispersive transport of solute dissolved in porewater; and chemical interactions between porewater and material.

2.1. Conceptual Model

Considering the cylindrical geometry of the bentonite layer around the heater and striving to save computational resources, a 1D model has been selected for the analysis of coupled THC processes in the bentonite (Figure 2).

The bentonite layer of 0.69 m thickness starts at the coordinates corresponding to the heater surface (z = 0.45 m) and ends at the coordinates corresponding to the bentonite–granite interface (z = 1.14 m). The bentonite in the model was discretized into finite elements to solve the coupled PDEs with the finite element method. Non-uniform discretisation was applied following [18], as indicated in Table 1. The duration of numerical simulation was set 18 years in line with the FEBEX in situ test duration.

Bentonite is a complex porous medium. As it contains large quantities of montmorillonite (for example, 61–95 wt.% in the MX-80 bentonite, 92 ± 4 wt.% in the Febex bentonite (smectite–illite mixed-layer) [20]), a mineral that has a significant impact on bentonite properties. The structure of montmorillonite is a unit made of an alumina octahedral sheet surrounded by two silica tetrahedral sheets [21]. Due to substitution of trivalent Al for tetravalent Si and divalent cations (e.g., Mg, Fe) for trivalent Al in the montmorillonite structure, a permanent negative lattice charge (the cation exchange capacity, CEC) persists [21]. The cations located in spaces between adjacent particles (interlayer spaces) compensate this charge [22]. The other type of reactive sites relevant in montmorillonite are surface hydroxyl groups existing along the edges of the clay platelets (“edge” or “broken bond” sites). According to [23], these sites provide ~10% of the CEC and can protonate and deprotonate so that the concentrations of neutral, protonated and deprotonated edge sites change as a function of pH. Saturation of bentonite with porewater of a particular composition will impose the re-distribution of dissolved species between the porewater and those bonded on montmorillonite. Apart from this main mineral, bentonite contains small amounts of other minerals as impurities, such as calcite, gypsum, and quartz. If impurities are not chemically inert, they can influence the chemical state of bentonite by mineral precipitation/dissolution processes and in turn affect its physical properties. In the long-term perspective, the precipitation of illite and dissolution of montmorillonite might also take place. Zheng et al. [12] analysed the illitisation in the FEBEX in situ test and concluded that illitisation would continue for 50 years but would not proceed. The major controlling factor in this case is chemical conditions, while temperature plays a secondary role. Samper et al. [13] concluded that consideration of smectite dissolution via kinetic law was significant in terms of the increase in aluminium concentration. The concentrations of other dissolved cations increased slightly, and concentrations of anions and neutral species showed no sensitivity to smectite dissolution [13]. Considering the outcomes from previous studies [12,13], illitisation and smectite dissolution were not considered in the current study. Apart from dissolved species of K, Na, Ca, Mg, and C, the inflowing granite porewater will contain dissolved gases, among which CO₂ is important. It could, depending on its level, impose CO_2g transport into or out of bentonite and affect the pH, concentration of dissolved HCO₃⁻, and calcite precipitation/dissolution.

Therefore, in order to improve the understanding of the THC state of the bentonite barrier under disposal conditions, the coupled thermo-hydro-chemical processes and tracer transport were analysed under different cases (Table 2). Case A was defined as a reference case where no chemical processes, such as mineral precipitation/dissolution, ion exchange, or surface complexation, were taken into account. Case B includes mineral precipitation/dissolution, ion exchange, and surface complexation but disregards CO₂ gas transport. Considering the CO₂ transport and possible degassing, Case C was developed, and the system behaviour was analysed with fixed CO₂ pressure on the bentonite–granite interface: 10⁻⁴ bar based on dissolved CO₂ concentration in the Spanish granite (Case C1) and with much higher flux of CO₂ from the host rock (fixed CO₂ partial pressure 10⁻¹ bar, Case C2).

2.2. Mathematical Model

In this study, the THC state of bentonite was analysed with COMSOL Multiphysics (v6.0), coupled with the geochemical solver Phreeqc (v. 3.7.3) via the iCP interface [19]. The main advantage of the interface was that it uses the strengths and capabilities of the mentioned codes, providing a tool for numerical simulation of relevant processes coupled with geochemistry. As stated in [19], iCP is written in Java and uses the IPhreeqcC++ dynamic library and the COMSOL Java-API. The geochemical reactions are solved in parallel by balancing the computational load over multiple threads.

For the analysis of the THC state of bentonite, it was decided to perform modelling of the two-phase flow of a miscible fluid (water and air) considering important phenomena such as heat transfer, gas dissolution, advective–diffusive transport in the gaseous phase, CO₂ transport, and chemical interactions (aqueous complexation, mineral dissolution/precipitation, ion exchange, and surface complexation).

2.2.1. Transport Processes

The transport of species and CO₂ in porous material could be advective/dispersive and diffusive in both (liquid and gaseous) phases. It was assumed that during the ventilation phase, there will be no gas pressure gradient imposed, and thus the CO₂ transport in the gaseous form will be governed by diffusion and the dissolved CO₂ will be transported by advection, molecular diffusion, and dispersion. No other volatile species were considered in the analysis. Mathematically, species transport could be expressed by a partial differential equation describing i species mass (

m_{i}

) conservation over time in space [24]:

\frac{\partial m_{i}}{\partial t} + q_{L} \cdot \nabla c_{i} = \nabla \cdot [(D_{D, i} + D_{e, i}) \nabla c_{i}] + R_{i},

(1)

where

c_{i}

is the dissolved species concentration (mol/m³),

c_{i, g}

is the species concentration in the gaseous form (mol/m³),

q_{L}

is Darcy velocity (m/s),

D_{D, i}

is the dispersion tensor (m²/s), and

D_{e, i}

is the effective diffusion coefficient (m²/s).

R_{i}

is a term describing the source and/or sink of the species (mol/m³s) due to chemical reactions, decay, etc.

Chemical processes in the bentonite were simulated using mass conservation law (mass action law), and all the reactions were assumed at chemical equilibrium, where activities of the species are related though the equilibrium constant K [25].

K = \frac{\prod_{p r o d u c t s} a_{p}^{ν_{p}}}{\prod_{r e a c t a n t s} a_{r}^{ν_{r}}},

(2)

where

a_{p}

and

a_{r}

are activities and

ν_{p}

and

ν_{r}

are stoichiometric coefficients of the products and reactants, accordingly.

With a chemical solver, the evaluation of mineral volume changes could be conducted. Changes in mineral volumes were converted into porosity changes:

n = 1 - V_{m i n e r a l s} - V_{i n e r t},

(3)

where

V_{m i n e r a l s}

is the volume of minerals 1 m³ of the porous medium and

V_{i n e r t}

is the volume of inert materials in 1 m³ of the porous medium.

Components of the dispersion tensor are [24]:

D_{D i i} = α_{L} \frac{q_{L i}^{2}}{|q_{L}|} + α_{T} \frac{q_{L j}^{2}}{|q_{L}|}

(4)

D_{D i j} = (α_{L} - α_{T}) \frac{q_{L i} \cdot q_{L j}}{|q_{L}|}

where

D_{D i i}

is the diagonal components of the dispersion tensor (m²/s);

D_{D i j}

is off-diagonal components of the dispersion tensor (m²/s); and

α_{L}

,

α_{T}

are longitudinal and transverse dispersivities (m).

The effective diffusion coefficient

D_{e, i}

for solutes in porous media was described as a function of liquid volume fraction and tortuosity [24]:

D_{e, i} = D_{L, i} = θ_{l} τ D_{o},

(5)

where

τ

is the medium tortuosity and

θ_{l}

is volumetric water content (dimensionless). For partially saturated porous media, tortuosity is related to water content through relationships such as [26]:

τ = \frac{θ_{l}^{7 / 3}}{n^{2}},

(6)

In a porous medium in which pores are fully water-filled, the liquid volume fraction

θ_{l}

is equal to porosity n. Meanwhile for an unsaturated medium, it is related to water saturation

S_{w}

(

θ_{l} = n S_{w}

). At the same time, the gas volume fraction in the case of an unsaturated porous medium is:

θ_{g} = n - θ_{l} = (1 - S_{w}) \cdot n

(7)

In case of CO₂ gas diffusion, mass is expressed as:

\frac{\partial m_{C O 2}}{\partial t} = \frac{\partial (θ_{l} c_{C O 2, d})}{\partial t} + \frac{\partial (θ_{g} c_{C O 2, g})}{\partial t}

(8)

According to Henry’s law, the concentration of CO₂ in gas phase (

c_{C O 2, g}

(mol/m³)) is linearly related to dissolved CO₂ concentration (

c_{C O 2, d}

(mol/m³)):

c_{C O 2, g} = \frac{c_{C O 2, d}}{H},

(9)

where H is Henry’s constant for CO₂ (dimensionless).

Then, effective diffusion coefficient could be expressed as follows:

D_{e, C O 2} = D_{L, C O 2} + \frac{1}{H} D_{G, C O 2},

(10)

where

D_{L, C O 2}

is diffusivity of dissolved CO₂ (m²/s),

D_{G, C O 2}

is diffusivity of CO₂ in gaseous phase (m²/s), which is described as:

D_{G, C O 2} = θ_{g} τ_{g} D_{g}^{f},

(11)

where

D_{g}^{f}

is the diffusivity (m²/s) of the fth gas species in an ideal gaseous phase, which according to [27], is given by:

D_{g}^{f} = \frac{R T}{3 \sqrt{2} π P N_{a v o g} d^{2}} \sqrt{\frac{8 R T}{π M}},

(12)

where

P

is the gaseous phase pressure,

N_{a v o g}

is Avogadro’s number,

d

is the molecular diameter,

M

is the molecular weight of the fth gas species. The tortuosity in the gaseous phase (

τ_{g}

) is computed from [26]:

τ_{g} = \frac{θ_{g}^{7 / 3}}{n^{2}},

(13)

With the molecular diameter of the gaseous species of 10⁻¹⁰ (m),

D_{g}^{f}

is 1.17 × 10⁻⁴ m²/s.

Henry’s constant dependence on temperature was described according to Van ’t Hoff equation:

H (T) = H (T_{r e f}) e x p (C (\frac{1}{T} - \frac{1}{T_{r e f}})),

(14)

where

H (T_{r e f})

= 3.4∙10⁻⁴ (mol/(m³∙Pa)), C = 2400,

T_{r e f}

= 25 °C according to [28].

2.2.2. Two-Phase Flow of Miscible Fluids

In this study, the bentonite hydration was modelled assuming that CO₂ will be insignificant in the total gas phase and the gas phase will be mainly water vapour and air. Then, the mass and energy balance over unit volume equations were defined for each component (water, air, and enthalpy) of this type:

\frac{\partial m_{*}}{\partial t} + \nabla {\vec{l}}_{*} = Q, m_{*} = m_{w}, m_{A}, h

(15)

The mass of water, air, and enthalpy were defined the following [29,30]:

m_{w} = m_{L} + m_{V} = ρ_{w} \frac{{e S}_{w}}{1 + e} + ρ_{V} \frac{e (1 - S_{w})}{1 + e}

(16)

m_{A} = m_{A L} + m_{A G} = ρ_{A} n (1 - S_{w} (1 - H))

h = \frac{ρ_{s} h_{s}}{1 + e} + m_{L} h_{w} + m_{V} h_{V} + m_{A} h_{A}

where

m_{w}

is the total mass of water as a sum of liquid water and vapour (kg/m³);

m_{A}

is the total mass of air as a sum of dissolved air and air in the gas phase (kg/m³);

e

is the void ratio;

h_{w}, {h_{V}, h}_{A}

is the enthalpy of water, vapour, and air (J/m³);

ρ_{w}

is water density (kg/m³);

ρ_{V}

is vapour density (kg/m³);

ρ_{A}

is air density(kg/m³);

S_{w}

is water saturation (m³/m³); and

H

is dimensionless Henry’s constant for air.

Fluxes of water, air (kg/(m²∙s)) and enthalpy (W/m²) were defined considering advective diffusive transport in both phases:

{\vec{l}}_{w} = (ρ_{w} q_{L}) + (ρ_{v} q_{G}) + i_{v}

(17)

{\vec{l}}_{A} = (ρ_{A} {H q}_{L}) + (ρ_{A} q_{G}) + i_{A}

{\vec{l}}_{h} = ρ_{W} h_{L} q_{L} + h_{V} (ρ_{V} q_{G} + i_{V}) + h_{A} (ρ_{A} q_{G} + i_{A} + ρ_{A} {H q}_{L}) - λ \nabla T

where

λ

is material thermal conductivity.

λ = n S_{w} λ_{w} + X_{A L} n S_{w} λ_{A} + n (1 - S_{w}) (λ_{A} + λ_{v}) + (1 - n) λ_{s}

(18)

X_{A L} = \frac{m_{A L}}{m_{A L} + m_{L}}

(19)

where

λ_{s}

,

λ_{W}

,

λ_{v}

,

λ_{A}

are the thermal conductivities of solid, water, vapour, and air, respectively (W/(m·K)).

The advective flow rate for liquid and gas phases (m/s) was described by the extended Darcy law:

q_{L} = - \frac{k_{L} k_{r L}}{μ_{w}} (\nabla P_{L} + ρ_{w} g)

(20)

q_{G} = - \frac{k_{G} k_{r G}}{μ_{G}} (\nabla P_{G})

(21)

where

k_{L}

is the intrinsic permeability of a porous medium (m²);

k_{r L}

is relative permeability of the material for the liquid phase, respectively (unitless);

k_{G} {, k}_{r G}

are intrinsic permeability (m²) and relative permeability (unitless) of the material for the gas phase, respectively;

μ_{w}

,

μ_{g}

are the dynamic viscosity of liquid and gas phases (Pa∙s);

P_{L}

,

P_{G}

—liquid and gas pressures (Pa).

Diffusive fluxes (kg/(m²∙s)) were described [29]:

i_{V} = - ρ_{G} τ_{V} n (1 - S_{w}) D_{V} \nabla (\frac{ρ_{V}}{ρ_{G}})

(22)

i_{A} = - i_{V}

(23)

where

τ_{V}

is vapour tortuosity (unitless) and

D_{V}

is the vapour diffusion coefficient (m²/s).

Substituting all the terms (Equations (17)–(23)) in to the Equation (16) and differentiating left hand side of the Equation (16) with respect of selected primary variables (liquid pressure, gas pressures and temperature), the governing equations for miscible fluid flow under non isothermal conditions were derived, which then need to be solved over finite element mesh.

2.3. Numerical Model

As was presented above, the reactive transport modelling requires information on flow conditions, which are governed by materials’ hydraulic properties (intrinsic and relative permeability, initial porosity, and water retention curve), transport-related properties (tortuosity, solute diffusivity in water, and saturation), chemistry-related data (equilibrium constants of chemical reactions considered in the analysis, and activity coefficients) as well as porewater solution composition, amount of minerals in the porous medium, exchanged ion concentrations, data on surface complexation sites, and initial and boundary conditions. Data on hydraulic, thermal, and transport-related properties of liquid, gas, and bentonite are provided in Table 3.

For the analysis of the porewater chemical composition (concentrations of dissolved species and activities), the thermodynamic database ThermoChimie v11a [31] was used. The ThermoChimie database was first developed in 1995 by Andra (the French National Radioactive Waste Management Agency). Now, it is being developed jointly with Radioactive Waste Management (RWM) (Didcot, UK) and ONDRAF/NIRAS (Bruxelles, Belgium). As stated in [32], the thermodynamic database ThermoChimie provides an accurate and consistent set of data specifically chosen for use in modelling the behaviour of radionuclides in waste packages, engineered barriers, and both the near surface and deep geosphere.

In this study, the ion exchange reactions were modelled with the Gaines–Thomas convention [25]. Ion selectivity constants for montmorillonite and initial concentrations of exchanged ions were taken from [13]. The surface complexation was modelled using data for site capacities (Table 4) and protolysis constants for montmorillonite given in [23] (Table 5). The data on initial porewater composition, initial mineral volume fractions, initial concentrations of exchanged ions at 25 °C, selectivity constants, and total concentrations of surface complexation sites are summarized in Table 4.

The following initial and boundary conditions were assumed (see also Figure 3):

With the heater temperature of 100 °C and the initial temperature in granite of 12 °C, the temperature gradient establishes quickly across the bentonite layer [12], with ~45 °C on the bentonite–granite interface. Thus, fixed-temperature boundary conditions were applied on the bottom (z = 0.45 m) boundary (T = 100 °C), representing the heater–bentonite interface, and on the upper (z = 1.14 m) boundary (T = 45 °C), representing the bentonite–granite interface.
Before transport simulation, the equilibrated porewater compositions with mineral phases, cation exchange sites, and surface protonation at 45 °C were derived with Phreeqc.
Initial suction (P_G − P_L) of 117.8 MPa was prescribed in the bentonite. With the defined water retention curve, this suction corresponded to 60% saturation.
No flow boundary conditions were applied to the heater–bentonite interface for water and air, and Dirichlet-type boundary conditions were applied for the bentonite–granite interface (P_G = P_L = 0.1 MPa).
No flux boundary condition was applied on the bottom model boundary representing the heater-bentonite interface for the transport of dissolved species and CO₂.
A fixed concentration boundary condition was applied on the upper boundary, representing the bentonite–granite interface for the transport of dissolved species.
Initially, after the equilibration of predefined porewater composition, the initial ion exchange concentration and minerals resulted in log(pCO₂) = –2 with the ThermoChimie database.
A fixed partial CO₂ pressure (log(pCO₂) = −4) boundary condition was applied on the top boundary, representing the bentonite–granite interface for CO₂ transport following CO₂ concentration in granite porewater (Case C1).
A fixed partial CO₂ pressure (log(pCO₂) = −1) boundary condition was applied on the top boundary, representing the bentonite–granite interface for CO₂ transport (Case C2).

Figure 3. Boundary conditions: (a) TH; (b) for transport.

3. Results

As was introduced above, for the analysis of the THC state of the bentonite barrier of a geological repository in granite, several cases were defined. They differed in the geochemical model and the consideration of CO₂ gas transport. Firstly, the importance of chemical processes relevant to the bentonite and related porosity changes were analysed by investigating the impact of THC processes interactions on the transport of non-reactive Cl and other species. Then, modelling of the evolution of CO₂ across the bentonite barrier was presented and the THC state potential impact on long-term safety was discussed.

3.1. Mineral Dissolution/Precipitation, Cation Exchange and Surface Complexation Effect on THC State of Bentonite

The effect of chemical processes, such as mineral dissolution/precipitation, ion exchange and surface complexation, on the THC state of the bentonite barrier was assessed through the analysis of the results when chemical reactions were not considered in the model (Case A) and the comparison with the case when chemical reactions were taken into account (Case B). The temperature evolution and hydration progress over time across the barrier for Case A is presented in Figure 4. In this case, the porosity remains constant, and thus the water content is directly proportional to the saturation evolution.

As can be seen from Figure 4, the hydration took place under a constant temperature gradient. The overall saturation of the barrier increased in time as the groundwater is infiltrating in the bentonite despite the evaporation taking place near the heater. Based on the modelling results, the barrier did not reach the full water saturated state S_w < 0.95 (with the prescribed water retention curve) within the simulation time (18 years of hydration).

The diffusivity of the initially present Cl⁻ within the bentonite was dependent on inflowing groundwater (increasing with increasing water saturation) and temperature (increasing with increased temperature). The advective–dispersive transport of the dissolved Cl resulted in Cl⁻ concentration profiles presented in Figure 5a. The concentration profiles of another solute (Na⁺) are presented in Figure 5b.

The modelling results showed that both Cl⁻ and Na+ concentrations in the bentonite gradually decreased due to their transport from the bentonite towards the granite, where a much lower concentration was retained.

The consideration of the mineral dissolution/precipitation, cation exchange, and surface complexation reactions (Case B) resulted in slight changes in porosity (Figure 6). The most significant porosity changes (increases) were observed near the bentonite–granite interface (~3%), while some decrease in porosity was observed near the heater. These porosity changes were driven by changes in the volume fraction (Figure 6a) of minerals, such as calcite, anhydrite, and quartz. Calcite and quartz dissolution was responsible for the porosity increase near the granite, while the porosity decrease near the heater was governed mainly by the anhydrite precipitation. According to [13], anhydrite is more stable (lower solubility) than gypsum for temperatures above a threshold temperature ranging from 43 °C to 56 °C, and thus, no gypsum mineral was observed in the analysed system. Anhydrite transformation into gypsum would be expected during cooling (not simulated in the current study).

These porosity changes were considered in the governing equations describing the water mass balance in the bentonite; therefore, the water content was affected. As can be seen in the Figure 7, the porosity changes were reflected in the slightly higher water content near the outer boundary of the bentonite. The highest difference in the water content between Case A and Case B was less than 4%. In addition, the impact was limited in space (the affected region was up to few cm from the bentonite–granite boundary) and were not influenced the transport of non-reactive anion Cl⁻ (compare Figure 5a and Figure 8a)).

Figure 8b presents the modelling results of the solute (Na⁺) transport for Case B. The Na⁺ concentration in porewater was directly affected by ion exchange with K⁺, Ca²⁺, and Mg²⁺. However, it was also indirectly affected by the mineral dissolution/precipitation because these processes led to changes in the concentration of calcium and magnesium in the porewater. Comparison of the Na⁺ concentration between Case A (see Figure 5b) and Case B (see Figure 5b) indicates that in Case B, a decreased concentration of Na⁺ in the zone near the heater could be seen after 1 or 5 y of heating, but for the longer time, the concentration near the heater was predicted to be higher.

As demonstrated, the consideration of additional chemical processes taking place in the bentonite (Case B) had no significant influence on the non-reactive Cl⁻ transport, but it could be important for other radionuclides whose chemical interactions (sorption) with porous media depend on porewater pH (for example, Cs [3], Ni [2], etc.). For pH evolution, surface protonation and mineral dissolution/precipitation reactions are relevant. The simulated pH profile across the bentonite at different time steps is presented in Figure 9.

Based on the modelling results, the pH evolution in the bentonite was quite different between Cases A and B. In Case B, the pH increase was observed near the boundary with the granite and then equilibrated at ~pH = 7, while the peak was not observed in Case A, and a gradual decrease in the pH was recorded. This clearly indicates that omitting and not considering some relevant chemical processes in the bentonite could affect indirectly the transport of radionuclides which sorption is pH sensitive.

3.2. Influence of CO₂ Transport

Consideration of CO₂ transport (in liquid and gas phases) in the numerical analysis allowed the CO₂ distribution across the barrier to be obtained. CO₂ diffusion in the gas phase is a fast process and will govern the transfer of the total CO₂. In Case C1, the CO₂ initial partial pressure in the bentonite (pCO₂ = 10⁻² bar) was higher than that imposed by the granite, and therefore, the gas diffusion was oriented towards the bentonite–granite interface. Due to degassing near the heater, the CO₂ content increased to some extent and then quickly decreased to the boundary value (pCO₂ = 10⁻⁴ bar) (Figure 10a).

In Case C2, the granite imposed a much higher partial CO₂ pressure, and thus, the influx of CO₂ into the bentonite was observed. Over time, its dissolution decreased with an increasing temperature, and some buildup of CO₂ was observed (Figure 10b). This distribution equilibrated quickly across the barrier and the profile was showing the decreasing trend towards the bentonite–granite interface. CO₂ distribution is important for mineral calcite: as the pressure of CO₂ gas in a system increases then calcite dissolves, and the decrease of CO₂ gas pressure results in calcite precipitation. However, if water temperature increases, then the amount of CO₂ dissolved in water decreases and leads to calcite precipitation. Due to the decrease in water temperature, the amount of CO₂ dissolved in water increases and causes calcite to dissolve. In our modelled system, these processes were not taking place one by one but occurred at the same time. Thus, a combined effect could be expected. Additionally, they were accompanied by other processes, such as dissolution/precipitation of calcite, anhydrite/gypsum, and cation exchange, where Ca²⁺ ions were involved. The modelling results showed that an imposed higher pCO₂ at the bentonite–granite interface led to a more significant calcite precipitation and dissolution of quartz close to the interface but not in the region close to the heater. In the inner part of the bentonite layer (closer to the heater), the anhydrite precipitation took over the other minerals in terms of increased volume fraction (Figure 11) in both cases C1 and C2. Subsequently, the resulting changes in porosity are presented in Figure 12.

In the case of a low pCO₂ at the interface, the calcite dissolution at the bentonite–granite interface led to a slight porosity increase (Figure 12a) near the interface. A slight increase in porosity was also obtained in Case B without consideration of CO₂ gas transport. An imposed higher pCO₂ would lead to a decreased porosity, which was also more variable farther from the bentonite–granite interface. As the porosity variations were observed near the bentonite–granite interface, it seems that CO₂ de-gassing due to temperature–dependent solubility did not have a significant impact.

Nevertheless, the observed porosity changes did not impose a significant impact on non-reactive tracer (Cl⁻) transport based on the concentration profiles presented in Figure 13.

The concentration evolution of other dissolved species over the bentonite showed a different trend. As an example, Ca²⁺ concentration profiles are presented in Figure 14. These solutes were transported by advection–dispersion as well as Cl⁻, but they also participated in cation exchange reactions and mineral dissolution/precipitation. The overall impact of those chemical processes resulted in the solute distribution over the barrier different than that for Cl.

Based on the results presented in Figure 14, different boundary conditions in terms of pCO₂ at the interface affected Ca²⁺ concentration distribution in the bentonite barrier, especially at longer times. Following the simulated Ca²⁺ concentration, a faster decrease in the concentration over the bentonite under the impact of a higher CO₂ flux from the interface would be expected. Under these conditions, the modelled concentration of Ca²⁺ near the heater was more than twice as low by the end of the simulation. The transported Ca²⁺ towards the interface reacted with the dissolved CO₂ forming the calcite and lowering the calcium concentration in the porewater. Ca²⁺ cations are also incorporated in anhydrite, which tend to precipitate in the region closer to the heater.

The modelled concentration of sulphate SO₄²⁻ was also significantly affected by the interface boundary condition. The model predicted a higher concentration than the initial in a particular location in the bentonite layer (Figure 15) even after 1 year of heating. Such results demonstrate that the resulting concentration distribution of reacting species was not determined solely by the transport processes. Peak concentrations were not observed near either of the interfaces (bentonite–granite or bentonite–heater) but within the barrier. Based on the modelling results, in time, the location of the peak concentration of sulphate was closer and closer to the heater.

The interplay of a number of chemical processes was also reflected in the pH evolution as presented in Figure 16. The modelled pH evolution over the bentonite with a low pCO₂ at the interface was similar to that obtained in Case B without consideration of CO₂ gas transport. A pH increase to 9.5 was observed close to the bentonite–granite interface, and then it gradually reduced to the initial pH of the bentonite. This suggests that if the CO₂ influx to the bentonite is low, CO₂ gas diffusion representation could be omitted in the numerical model, and therefore, computational resources could be saved. Meanwhile, in the case with an imposed high pCO₂ (Case C2), a non-linear decrease in the pH was observed in the direction further away from that fixed on the interface. By the end of simulation, the pH in the outer part of the barrier was higher than the initial, but in the region close to the heater it was lower than the initial pH of the bentonite (~7). The higher imposed CO₂ influx resulted in an increased pH over a larger part of the bentonite in comparison to the fixed low CO₂ flux.

Observed differences in the concentration profiles of solutes (except Cl⁻) clearly indicate the importance of proper assessment of the conditions on the bentonite from the surrounding environment. While the solute concentration distribution was affected throughout the bentonite barrier, the mineral precipitation/dissolution and porosity were affected only close to the bentonite–granite interface. Therefore, it might be beneficial to extend the numerical model taking into account the surrounding host rock, as its properties and conditions prevailing in situ could be affected by the bentonite barrier and a transitional zone could be developed. The importance of the boundary conditions in terms of CO₂ should be kept in mind while interpreting the measurement results as the boundary conditions accepted in the numerical model should be representative of the system where the measurements came from. In the future, modelling of the FEBEX in situ test considering more detailed experiment geometry and heating history is foreseen before model application for long-term predictions.

Based on the modelling results, it seems that chemical processes taking place in the bentonite under elevated temperature did not affect significantly the barrier properties, such as porosity, relevant to radionuclides transport. From the long-term safety perspective, the increase in porosity is relevant as this could enhance radionuclide diffusivity from the repository. The transport of non-sorbing radionuclides, such as I-129, and Cl-36, would not be affected significantly by the porosity changes induced by cation exchange, surface complexation, and mineral dissolution/precipitation with and without considering CO₂ flux; however, porosity changes due to other processes (swelling, corrosion, etc.) should be assessed. For sorbing radionuclides—such as Cs-135, Sr-90, U-238, and others—considered chemical processes and accepted boundary conditions are important, as these will determine the chemical state of the bentonite porewater and radionuclide sorption conditions (pH, ionic strength, etc.).

4. Conclusions

A numerical model for THC processes in the bentonite barrier of the geological repository in granite was developed to analyse the behaviour of the THC state of the bentonite under repository conditions and its impact on long-term safety. The model was implemented in COMSOL Multiphysics, coupled with the geochemical solver Phreeqc via the iCP interface. The model formulation was based on the two-phase flow of miscible fluids (water and air) considering important phenomena such as heat transfer, gas dissolution, advective–diffusive transport in the gaseous phase, CO₂ transport, and chemical interactions (speciation, mineral dissolution/precipitation, ion exchange, and surface complexation). Three cases were analysed: a TH model with tracer transport (Case A), a THC model without consideration of CO₂ transport (Case B), and a THC model where CO₂ transport in a gas phase was taken into account (Case C). Two variants of Case C were distinguished depending on boundary conditions (low CO₂ pressure at the bentonite–granite interface (Case C1) and high CO₂ pressure at the bentonite–granite interface (Case C2). After the analysis of the modelling results, the following conclusions can be drawn:

Mineral dissolution/precipitation induced slight changes in the bentonite porosity. The most significant porosity increase was observed near the bentonite–granite interface (by ~3%) in Case B (THC without consideration of CO₂ gas transport);
Changes in porosity depend on the partial CO₂ pressure at the bentonite–granite boundary. In the case of a low pCO₂ at the bentonite–granite interface (Case C1), the calcite dissolution led to a slight porosity increase, while a higher pCO₂ (Case C2) led to decreased porosity. In the latter case, the changes were observed at the larger distance from the bentonite–granite interface;
Consideration of chemical processes taking place in the bentonite had no significant influence on non-reactive Cl transport for the considered period of time, but it would be important for other radionuclides whose sorption in porous media depends on porewater pH;
CO₂ de-gassing due to temperature-dependent solubility did not have a significant impact on porosity change;
Different boundary conditions in terms of pCO₂ at the bentonite–granite interface had a larger impact on the sulphate concentration compared to cation distribution in the bentonite barrier;
Modelling results presented in this study could be useful for further interpretation of the FEBEX test measurement results before model application for large scale predictions;
In the future, the modelling of the FEBEX in situ test considering more detailed experiment geometry and heating history is foreseen.

Author Contributions

Investigation, A.N. and D.G.; writing—original draft preparation, A.N.; visualisation, A.N.; writing—review and editing, D.G. and G.P. All authors have read and agreed to the published version of the manuscript.

Funding

The EURAD program has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 847593, as well as the Lithuanian State Budget.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

References

Bennett, D.G.; Gens, R. Overview of European concepts for high-level waste and spent fuel disposal with special reference waste container corrosion. J. Nucl. Mater. 2008, 379, 1–8. [Google Scholar] [CrossRef]
Zhao, X.; Qiang, S.; Wu, H.; Yang, Y.; Shao, D.; Fang, L.; Liang, J.; Li, P.; Fan, Q. Exploring the Sorption Mechanism of Ni(II) on Illite: Batch Sorption, Modelling, EXAFS and Extraction Investigations. Sci. Rep. 2017, 7, 8495. [Google Scholar] [CrossRef] [PubMed]
Belousov, P.; Semenkova, A.; Egorova, T.; Romanchiuk, A.; Zakusin, S.; Dorzhieva, O.; Tyupina, E.; Izosimova, Y.; Tolpeshta, I.; Chernov, M.; et al. Cesium Sorption and Desorption on Glauconite, Bentonite, Zeolite, and Diatomite. Minerals 2019, 9, 625. [Google Scholar] [CrossRef]
Thatcher, K.E.; Bond, A.E.; Norris, S. Pore pressure response to disposal of heat generating radioactive waste in a low permeability host rock. Int. J. Rock. Mech. Min. Sci. 2020, 135, 104456. [Google Scholar] [CrossRef]
Bumbieler, F.; Plua, C.; Tourchi, S.; Vu, M.-N.; Vaunat, J.; Gens, A.; Armand, G. Feasibility of constructing a full-scale radioactive high-level waste disposal cell and characterization of its thermo-hydro-mechanical behavior. Int. J. Rock. Mech. Min. Sci. 2021, 137, 104555. [Google Scholar] [CrossRef]
Müller, H.R.; Garitte, B.; Vogt, T.; Kohler, S.; Sakaki, T.; Weber, H.; Spillmann, T.; Herthrich, M.; Becker, J.K.; Giroud, N.; et al. Implementation of the full-scale emplacement (FE) experiment at the Mont Terri rock laboratory. Swiss J. Geosci. 2017, 110, 287–306. [Google Scholar] [CrossRef]
Li, X.L.; Bastiaens, W.; Van Marcke, P.; Verstricht, J.; Chen, G.J.; Weetjens, W.; Sillen, X. Design and development of large-scale insitu PRACLAY heater test and horizontal high-level radioactive waste disposal gallery seal test in Belgian HADES. J. Rock. Mech. Geotech. Eng. 2010, 2, 103–110. [Google Scholar] [CrossRef]
Grimsel Test Site. High Temperature Effects on Bentonite Buffers (HotBENT)—Aims & Objectives. Available online: https://www.grimsel.com/gts-projects/hotbent-high-temperature-effects-on-bentonite-buffers/hotbent-introduction (accessed on 8 March 2024).
Sánchez, M.; Pomaro, B.; Gens, A. Coupled THM analysis of a full-scale test for high-level nuclear waste and spent fuel disposal under actual repository conditions during 18 years of operation. Géotechnique 2023, 73, 418–438. [Google Scholar] [CrossRef]
Bosch, J.A.; Qiao, Y.; Ferrari, A.; Laloui, L. Thermo-hydro-mechanical analysis of the complete lifetime of the bentonite barrier in the FEBEX in-situ test. Geomech. Energy Environ. 2023, 34, 100472. [Google Scholar] [CrossRef]
Zhang, L.; Samper, J.; Montenegro, L. A coupled THC model of the FEBEX in situ test with bentonite swelling and chemical and thermal osmosis. J. Contam. Hydrol. 2011, 126, 45–60. [Google Scholar] [CrossRef]
Zheng, L.; Fernández, A.M. Prediction of Long-Term Geochemical Change in Bentonite Based on the Interpretative THMC Model of the FEBEX In Situ Test. Minerals 2023, 13, 1522. [Google Scholar] [CrossRef]
Samper, J.; Mon, A.; Montenegro, L. A revisited thermal, hydrodynamic, chemical and mechanical model of compacted bentonite for the entire duration of the FEBEX in situ test. Appl. Clay Sci. 2018, 160, 58–70. [Google Scholar] [CrossRef]
Zheng, L.; Xu, H.; Rutquist, J.; Reagan, M.; Birkholzer, J.; Villar, V.; Fernandez, A.M. The hydration of bentonite buffer material revealed by modeling analysis of a long-term in situ test. Appl. Clay Sci. 2020, 185, 105360. [Google Scholar] [CrossRef]
Sedighi, M.; Thomas, H.R.; Vardon, P.J. Reactive transport of chemicals in compacted bentonite under non-isothermal water infiltration. J. Geotech. Geoenviron. Eng. 2018, 144, 04018075. [Google Scholar] [CrossRef]
Samper, J.; Zheng, L.; Montenegro, L.; Fernandez, A.M.; Rivas, P. Coupled thermo-hydro-chemical models of compacted bentonite after FEBEX in situ test. Appl. Geochem. 2008, 23, 1186–1201. [Google Scholar] [CrossRef]
Samper, J.; Mon, A.; Ahusborde, E.; Yu, H.; Narkuniene, A.; Hokr, M.; Montenegro, L.; Amaziane, B.; Ossmani, M.E.; Xu, T.; et al. Multiphase flow and reactive transport benchmark for radioactive waste disposal. Environ. Earth Sci. (under review).
Nardi, A.; Idiart, A.; Trinchero, P.; de Vries, L.M.; Molinero, J. Interface COMSOL-PHREEQC (iCP), an efficient numerical framework for the solution of coupled Multiphysics and geochemistry. Comput. Geosci. 2014, 69, 10–21. [Google Scholar] [CrossRef]
Huertas, F.; Farina, P.; Farias, J.; Garcia-Sineriz, J.L.; Villar, M.V.; Fernandez, A.M.; Martin, P.L.; Elorza, F.J.; Gens, A.; Sanchez, M.; et al. Full-Scale Engineered Barriers Experiment Updated Final Report 1994–2004; Tech. Publ 05-0/2006; Enresa: Madrid, Spain, 2006; p. 590. [Google Scholar]
Villar, M.V.; Arnaud, D.; Besuelle, P.; Cernochova, K.; Collin, F.; Cuevas, J.; Cuss, R.; de Lesquen, C.; El Tabbal, G.; Graham, C.; et al. D7.2 HITEC. Updated State-of-the-Art on THM Behaviour of i) Buffer Clay Materials and of ii) Host Clay Materials; Deliverable D7.2 HITEC (2023); EURAD Project, Horizon: Bruxelles, Belgium, 2020; No 847593; p. 121. [Google Scholar]
Bradbury, M.H.; Baeyens, B. Porewater chemistry in compacted re-saturated MX-80 bentonite. J. Contam. Hydrol. 2003, 61, 329–338. [Google Scholar] [CrossRef]
Birgersson, M. A general framework for ion equilibrium calculations in compacted bentonite. Geochim. Cosmochim. Acta 2017, 200, 186–200. [Google Scholar] [CrossRef]
Bradbury, M.H.; Baeyens, B. A mechanistic description of Ni and Zn sorption on Na–montmorillonite: Part II. Modelling. J. Contam. Hydrol. 1997, 27, 223–248. [Google Scholar] [CrossRef]
Subsurface Flow Module. User’s Guide. Comsol Multiphysics. 2021. Available online: https://doc.comsol.com/6.1/doc/com.comsol.help.ssf/SubsurfaceFlowModuleUsersGuide.pdf (accessed on 20 March 2024).
Parkhurst, D.L.; Appelo, C.A.J. User’s Guide to PHREEQC (Version 2): A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations; USSG: Denver, CO, USA, 1999; 312p. [Google Scholar]
Millington, R.J.; Quirk, J.M. Transport in Porous Media. Transp. Porous Media 1960, 49, 377–378. [Google Scholar] [CrossRef]
Lasaga, A.C. Kinetic Theory in the Earth Sciences; Princeton Series in Geochemistry; Princeton University Press: Princeton, NJ, USA, 1998; 811p, ISBN 0-691-03748-5. [Google Scholar]
NIST Chemistry WebBook, SRD 69. Available online: https://webbook.nist.gov/cgi/cbook.cgi?ID=C124389&Mask=10 (accessed on 21 March 2024).
Navarro, V.; Asensio, L.; Gharbieh, H.; De la Morena, G.; Pulkkanen, V.M. Development of a multiphysics numerical solver for modeling the behavior of clay-based engineered barriers. Nucl. Eng. Technol. 2019, 51, 1047–1059. [Google Scholar] [CrossRef]
Pollock, D.W. Simulation of fluid flow and energy transport processes associated with high-level radioactive waste disposal in unsaturated alluvium. Water Resour. Res. 1986, 22, 765–775. [Google Scholar] [CrossRef]
Thermodynamic Database Andra—RWM—Ondraf Thermo_Chimie. Available online: www.thermochimie-tdb.com (accessed on 22 February 2023).
Rodríguez, J.; Colàs, E.; Duro, L.; Fuller, A.J.; Harvey, L.; Thermo Chimie. Track-Changes and Track-Error Document: From Version 10d to Version 11. Technical Report TCIII-2022-07E_vs1. April 2022.. Available online: https://www.thermochimie-tdb.com/pages/version.php (accessed on 22 February 2023).

Figure 1. FEBEX in situ test layout.

Figure 2. Model geometry.

Figure 4. Bentonite hydration modelling results: temperature (a) and water content (b) profiles at time steps 1, 5, and 18 years.

Figure 5. Concentration profiles at different times (1 y, 5 y, 18 y of heating) for Case A: (a) Cl⁻ A; (b) Na⁺.

Figure 6. Impact of consideration of mineral dissolution, cation exchange and surface complexation (Case B) on mineral volume fraction and porosity: (a) volume fraction profile of separate minerals; (b) porosity profile (after 1 year, 5 years, 18 years of heating).

Figure 7. Modelled profile of water content evolution in the bentonite due to hydration under nonisothermal conditions.

Figure 8. Concentration profiles at different times (1 y, 5 y, 18 y of heating) for Case B: (a) Cl⁻; (b) Na⁺.

Figure 9. pH profile in the bentonite at t = 1, 5 and 18 y: (a) without consideration of mineral dissolution, cation exchange, and surface complexation, Case A; (b) with consideration of mineral dissolution, cation exchange, and surface complexation, Case B.

Figure 10. Partial pressure of CO₂ across the bentonite at different times (after 0.001 years, 1 year, 5 years, 18 years of simulation): (a) at boundary pCO₂ = 10⁻⁴ bar (Case C1); (b) at boundary pCO₂ = 10⁻¹ bar (Case C2).

Figure 11. Simulated concentration of minerals by the end of simulation: quartz (1), calcite (2), anhydrite (3) for Case C1 (a) and Case C2 (b) (not to scale).

Figure 12. Impact of pCO₂ of incoming porewater (at bentonite–granite interface) on barrier porosity evolution after 1 y (t1), 5 y (t2), 18 y (t3) of heating in Case C1 (a) and Case C2 (b) (not to scale).

Figure 13. Simulated profiles of dissolved Cl⁻ concentration after 1 year, 5 years, 18 years of heating: (a) Case C1; (b) Case C2.

Figure 14. Simulated evolution of dissolved Ca²⁺ concentration after 1 year, 5 years, 18 years of heating: (a) Case C1; (b) Case C2.

Figure 15. Simulated evolution of dissolved SO₄²⁻ concentration in the bentonite: (a) Case C1; (b) Case C2.

Figure 16. Simulated profiles of pH: (a) Case C1; (b) Case C2.

Table 1. Discretisation.

Range	Mesh Size, m
0.45 m–0.46 m	0.000625
0.46 m–0.49 m	0.0025
0.49 m–0.59 m	0.005
0.59 m–1.00 m	0.01
1.00 m–1.10 m	0.005
1.10 m–1.13 m	0.0025
1.13 m–1.14 m	0.00625

Table 2. Summary of analysed cases.

Cases	Transported Species	Geochemical Model
Case A	Na, Cl	Aqueous complexation, acid–base reactions
Case B	Na, Cl, C, Ca, S, Si, Mg, K	Aqueous complexation, acid–base reactions, mineral dissolution/precipitation, ion exchange, surface complexation
Case C1	Na, Cl, C, Ca, S, Si, Mg, K, CO₂ gas diffusion	Aqueous complexation, acid–base reactions, mineral dissolution/precipitation, ion exchange, surface complexation, fixed pCO₂ at bentonite–granite boundary (log (pCO₂) = −4)
Case C2	Na, Cl, C, Ca, S, Si, Mg, K, CO₂ gas diffusion	Consideration of mineral dissolution/precipitation, ion exchange, fixed pCO₂ at bentonite–granite boundary (log (pCO₂) = −1)

Table 3. Data on hydraulic, thermal, transport related properties of liquid, gas, and bentonite used in the study.

Parameter [13,18,29]	Value or Dependence
Initial porosity of bentonite	0.41
Solid density (kg/m³) (T in Celsius)	$2780 e^{(- 2 \cdot 10^{- 5} (T - 12))}$
Water density (kg/m³) (T in Celsius)	$ρ_{w} = 998.2 \cdot e x p (5 \cdot 10^{- 7} (P_{L} - 100) - 2.1 \cdot 10^{- 4} (T - 12))$
Vapour and air density (kg/m³) (T in Celsius)	$ρ_{g} = ρ_{v} + ρ_{A}$ $ρ_{v} = \frac{e^{(0.06374 T - 0.1634 {\cdot 10^{- 3} T}^{2})}}{194.4} e^{- (2.16677 (P_{g} - P_{L})) / (ρ_{w} (T + 273.15))}$ $ρ_{A} = \frac{M_{a} P_{A}}{R T}$
Thermal properties [13,18,29]
Thermal conductivity of the liquid (W/m °C)	1.5
Thermal conductivity of the air (W/m °C)	2.6·10⁻²
Thermal conductivity of the vapour (W/m °C)	4.2·10⁻²
Thermal conductivity of the solid (W/m °C)	1.23
Specific heat of the water (J/kg °C)	4202
Specific heat of the air (J/kg °C)	1000
Specific heat of the vapour (J/kg °C)	1620
Specific heat of the solid (J/kg °C)	835.5
Vaporisation enthalpy (J/kg)	2.45·10⁶
Thermal compressibility of the water (°C⁻¹)	2.1·10⁻⁴
Thermal compressibility of the solid (°C⁻¹)	2·10⁻⁵
Hydraulic properties [13,18]
Intrinsic permeability for liquid flow (m²)	$k_{L} = k_{o} \frac{φ^{3}}{{(1 - φ)}^{2}} \frac{{(1 - φ_{o})}^{2}}{{φ_{o}}^{3}}$ $, k_{o}$ = 3.75 10⁻²¹ m²
Relative permeability to liquid	$k_{r L} = {S_{W}}^{3}$
Intrinsic permeability for gas flow (m²)	5·10⁻¹⁰
Relative permeability to gas (m²)	$k_{r G} = (1 - S_{W})^{3}$
Van Genuchten retention curve	$S_{w} = S_{r} + \frac{S_{s} - S_{r}}{{(1 + {(α (P_{g} - P_{l})}^{(\frac{1}{1 - b})})}^{b}}$ $S_{r} = 0.05$ $, S_{s} = 0.95$ $α = 5 \cdot 10^{- 5}$ ${kPa}^{- 1}, b = 0.21$
Liquid viscosity (kg/m∙s) (T in Celsius)	$661.2 \cdot 10^{- 3} {(T + 44)}^{- 1.562}$
Gas viscosity (kg/m∙s)	1.76·10⁻⁵
Vapour tortuosity factor	0.10
Henry’s constant for air (T in Celsius)	${H = ρ^{l} \cdot (T + 273.15) \cdot 7.1281183 \cdot 10}^{- 11} \cdot e^{\frac{1997.32}{(T + 273.15)}}$
Henry’s constant for CO₂ at 25 °C (mol/(m³∙Pa))	3.4·10⁻⁴
Transport-related properties [18,29]
Molecular diffusion in water $D_{0} (T)$ (m²/s)	$D_{o} (T) = D_{o} (T_{r e f}) \frac{T}{T_{r e f}} \frac{μ_{0}^{l}}{μ^{l}}$ $, D_{o} (T_{r e f}) = 2 \cdot 10^{- 11}$ $, T_{r e f}$ = 25 °C
Longitudinal dispersivity (m)	0.01
Molecular diameter of the gases species (m)	10⁻¹⁰
Vapour diffusivity (m²/s)	$D_{V} = 5.9 \cdot 10^{- 6} \frac{T^{2.3}}{P_{G}}$

Table 4. Initial porewater composition, initial mineral volume fractions, initial concentrations of exchanged ions, selectivity constants, and total concentrations of surface complexation sites [13].

	Bentonite	Granite
Initial porewater composition
pH	7.72	8.35
Na⁺ (mol/L)	1.3 × 10⁻¹	3.83 × 10⁻⁴
K⁺ (mol/L)	1.73 × 10⁻³	7.83 × 10⁻⁶
Ca²⁺ (mol/L)	2.23 × 10⁻²	1.83 × 10⁻⁴
Mg²⁺ (mol/L)	2.33 × 10⁻²	1.33 × 10⁻⁶
HCO₃⁻ (mol/L)	4.13 × 10⁻⁴	3.93 × 10⁻⁶
SO₄²⁻ (mol/L)	3.23 × 10⁻²	7.93 × 10⁻⁵
Cl⁻ (mol/L)	1.63 × 10⁻¹	1.33 × 10⁻⁵
SIO_2(aq) (mol/L)	1.13 × 10⁻⁶	1.43 × 10⁻⁴
Initial volume fraction of the minerals (%)
Calcite	1	5
Quartz	4.5	20
Anhydrite	0	0
Gypsum	0.016	0
Halite	0	0
Initial ion exchange concentration meq/100 g bentonite
Na⁺	21.10	-
K⁺	1.94	-
Ca²⁺	31.31	-
Mg²⁺	41.41	-
Total concentration of surface complexation sites (mol/kg)
$\equiv S^{s} O H$	2 × 10⁻³	-
$\equiv S^{w 1} O H$	4 × 10⁻³	-
$\equiv S^{w 2} O H$	4 × 10⁻³	-

Table 5. Data for cation exchange and surface complexation reactions.

	Selectivity Constant K_Na-cation [13]	Log K
Cation exchange reactions on montmorillonite
Na⁺ + X-K ⇔ K⁺ + X-Na	0.138
Na⁺ + 0.5 X₂-Ca ⇔ 0.5 Ca²⁺ + X-Na	0.294
Na⁺ + 0.5 X₂-Mg ⇔ 0.5 Mg²⁺ + X-Na	0.288
Or expressed with half reactions:
X⁻ = X⁻		0.0
K⁺ + X⁻ = KX		0.86
Ca²⁺ + 2X⁻ = CaX₂		1.064
Mg²⁺ + 2X⁻ = MgX₂		1.082
Surface complexation		[23]
$\equiv S^{s} O H_{2}^{+} \leftrightarrow \equiv S^{s} O H + H^{+}$		−4.5
$\equiv S^{s} O^{-} + H^{+} \leftrightarrow \equiv S^{s} O H$		7.9
$\equiv S^{w 1} O H_{2}^{+} \leftrightarrow \equiv S^{w 1} O H + H^{+}$		−4.5
$\equiv S^{w 1} O^{-} + H^{+} \leftrightarrow \equiv S^{w 1} O H$		7.9
$\equiv S^{w 2} O H_{2}^{+} \leftrightarrow \equiv S^{w 2} O H + H^{+}$		−6.0
$\equiv S^{w 2} O^{-} + H^{+} \leftrightarrow \equiv S^{w 2} O H$		−10.5

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Narkuniene, A.; Grigaliuniene, D.; Poskas, G. THC Modelling of Bentonite Barrier of Geological Repository in Granite and Its Impact on Long-Term Safety. Appl. Sci. 2024, 14, 7851. https://doi.org/10.3390/app14177851

AMA Style

Narkuniene A, Grigaliuniene D, Poskas G. THC Modelling of Bentonite Barrier of Geological Repository in Granite and Its Impact on Long-Term Safety. Applied Sciences. 2024; 14(17):7851. https://doi.org/10.3390/app14177851

Chicago/Turabian Style

Narkuniene, Asta, Dalia Grigaliuniene, and Gintautas Poskas. 2024. "THC Modelling of Bentonite Barrier of Geological Repository in Granite and Its Impact on Long-Term Safety" Applied Sciences 14, no. 17: 7851. https://doi.org/10.3390/app14177851

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

THC Modelling of Bentonite Barrier of Geological Repository in Granite and Its Impact on Long-Term Safety

Abstract

1. Introduction

2. Materials and Methods

2.1. Conceptual Model

2.2. Mathematical Model

2.2.1. Transport Processes

2.2.2. Two-Phase Flow of Miscible Fluids

2.3. Numerical Model

3. Results

3.1. Mineral Dissolution/Precipitation, Cation Exchange and Surface Complexation Effect on THC State of Bentonite

3.2. Influence of CO₂ Transport

4. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

Article Menu

THC Modelling of Bentonite Barrier of Geological Repository in Granite and Its Impact on Long-Term Safety

Abstract

1. Introduction

2. Materials and Methods

2.1. Conceptual Model

2.2. Mathematical Model

2.2.1. Transport Processes

2.2.2. Two-Phase Flow of Miscible Fluids

2.3. Numerical Model

3. Results

3.1. Mineral Dissolution/Precipitation, Cation Exchange and Surface Complexation Effect on THC State of Bentonite

3.2. Influence of CO2 Transport

4. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

3.2. Influence of CO₂ Transport