Development of a THMC code for bentonites in COMSOL Multiphysics

The proposed use of active clays for the isolation of radioactive wastes in deep geological repositories has been followed by a deeper understanding of this type of soils. This increased knowledge has led to the need for both conceptual and numerical models capable of capturing the main trends in behaviour and the different couplings between different physical-chemical phenomena. In addition, the model must have a high degree of flexibility that enables it to accommodate future developments or new relevant phenomena. This work presents a numerical THMC code developed entirely on the COMSOL Multiphysics numerical implementation platform, which provides the required adaptability. This model includes, for the first time in this environment, a reactive transport model in unsaturated porous media for a relevant geochemical system (consistent with the MX-80 bentonite) together with a THM model based on a double porosity approach. The chemical potentials of water and solutes are used for the definition of thermodynamic equilibria between both porosity levels. Trends in the behaviour of a bentonite sample under oedometric conditions are satisfactorily simulated in response to a process of saturation and change in salinity conditions. Variations in swelling pressure, porosity distribution or dissolution/precipitation of the main accessory minerals are analysed and explained by means of the proposed conceptual model.


Introduction
Active clays, especially compacted bentonite, are considered one of the most important engineering barrier elements in many deep geological repository concepts for spent nuclear fuel due to their characteristics such as low hydraulic conductivity, high swelling potential and high retention capacity. [1][2][3][4]. This material will be used in several repository concepts in crystalline and argillaceous host rock, usually both as buffer and/or backfill elements. MX-80 has been one of the most studied bentonites [5].
Hydro-mechanical models have been used to study the macroscopic behaviour of this bentonite focused on changes on deformability and water flow [6]. However, several studies [7,8] show that the hydro-mechanical response (hydraulic conductivity and swelling pressure) of the bentonites depends on the level of salinity and chemical composition of the surrounding groundwaters in deep geological repositories. The composition of these waters can be very different [9,10]. This aspect together with the complex geochemical system of the MX-80 [11] makes it necessary to develop models that consider the effect of chemical conditions [12].
In this context, the main objective of this work is the development of a THMC code that allows predicting the behaviour of the bentonite barrier from an unsaturated state during an extensive timeframe. For this purpose, the model of Navarro et al (2017) [13] has been improved by the implementation of an extended geochemical model typical of MX-80 bentonite [14] and the surrounding groundwaters. The code has been implemented in COMSOL Multiphysics. A synthetic case has been simulated to check the scope of the numerical tool: a saturation test with low salinity water of a compacted MX-80 sample. The results obtained in the numerical analysis show the main expected trends of behaviour of the bentonite.

Conceptual model
A double porosity model [15] was selected to conceptualize the structure of the compacted bentonite: (i) macrostructure as inter-aggregate void space and microstructure as the intra-aggregate void space. The BExM [16] model was used as stress-strain constitutive framework. BBM [17] was used to simulate the macrostructural strains. Both models were adapted following [13].
In the proposed model, macrostructural water (free water) can flow under conventional hydrodynamic gradients. The chemical potential of macrostructural water is defined as where μ VO is the chemical potential of pure water at temperature T (absolute value), WMM is the molar mass of water, ρ w is the density of free water, s M and s MO are de matric and osmotic suctions in the macrostructure respectively. The microstructural space is saturated with water that can be considered immobile under conventional hydrodynamic gradients. In the same way, the chemical potential of the microstructural water can be defined as where ρ m is the microstructural water density, p is the net mean stress, Δs mNCC that defines the increase in microstructural suction due to the contribution to the chemical potential by the ions in excess of the CEC (defined in [18]) and s mS is the "structural" suction of the microstructure and it is defined by a state function ( Figure 1) characterised by microstructural void ratio [19]. constitutive law (π = π(e m )). Adapted from [20].
Adopting equilibrium between structural levels, the chemical potentials are equal in macrostructural, Eq. 1, and microstructural water, Eq. 2, (assuming ρ w = ρ m ) Organizing the Eq. 3 in a different arrangement, and identifying s mS with π The right-hand side term can be defined as a "boundary swelling pressure", π Β . When the equilibrium between macrostructural and microstructural water is not fulfilled, π and π B are not equal. The water mass transfer between the macro and microstructural level is then defined by a kinetic law defined in [21] given by the expression where H and C are constitutive parameters. H determines the transfer coefficient when the water mass exchange process is finalized (π = π B ) and C determines how the mass transfer term changes as π reaches π B .
The mechanical equilibrium equation (with its associated boundary conditions) can be solved if the changes in the total strain (dε ε ε ε) are known. This magnitude can be computed using an additive formulation given by the expression where dİ M the increase of macrostructural strain, dİ m is the increase of microstructural strain, and dİ Mm defines the rearrangement of the macrostructure because of dİ m . It is accepted that the rearrangement of the macrostructural level does not induce strains in the microstructural level.
The macrostructural mechanical constitutive model, defined in [22], is based on the associated plasticity approach of the Barcelona Basic Model (BBM [23,24]). Given the von Mises stress q, the slope of the critical state line M (a material parameter), the net mean stress p, and the net mean yield stress at the current suction p O , the yield function F is defined as where the increase in tensile strength with suction p S was calculated using the expression S M p k s = ⋅ with k being a material parameter and s M being the macrostructural suction. The value of the net mean yield stress at the current suction p O was calculated using the expression where the reference stress p C is a material parameter, and the evolution of the saturated pre-consolidation stress p O * is calculated from the macrostructural plastic strain rate using the hardening law where λ(0), r, β, κ io and α i are material parameters.
The macrostructural elastic volumetric strain dε MV e was obtained using the expression e M MV p s ds dp where the "basic" value of the "stress" bulk modulus K p was obtained using the expression and the "suction" bulk modulus K s was obtained with the expression Here, κ S ( s M , p ) is the elastic stiffness for changes in suction, which was calculated as follows where κ So , α Sp , p REF and α SS are material parameters.
The microstructural constitutive model is given by the state surface depicted in Figure 1, where the parametrization of the function π = π (e m ) can be found in [19,25].
The mass balance of the macrostructural water and the chemical master species are defined using the equations described in [26].
The micro-macro coupling term dε ε ε ε Mm includes both the macrostructural plastic volumetric strain, dε ε ε ε Mm1 , and the free swelling, dε ε ε ε Mm2 . The first term is caused by the macrostructural packing when the microstructural void ratio changes [27] and can be computed as proposed in the Barcelona Expansive Model (BExM [28]) from the expression where the function f β is the interaction function and β is the subscript that denotes the interaction function to be chosen: I for the "suction increase" conditions and D for the "suction decrease" conditions [29]. The second term, dε ε ε ε Mm2 , is caused by the destructuration of the microstructure at low confinement conditions. The increase of e M due to free swelling processes (identified as "Δe Mm2 " to differentiate it from the increase of e M due to other processes) occurs only when e m is high; that is, when π and the confinement π B are low. The model proposed assumes that, for Δe Mm2 to occur, π B must be less than a reference value π B,REF , which is a material parameter of the model. Under these conditions, in a monotonic swelling path, the following applies where stiffness κ π (dimensionless) is another material parameter for this swelling caused by microstructural destructuration. Multicomponent diffusion transport has been adapted following the formulation from [30]. The chemical speciation in macrostructural water has been calculated following the mathematical procedure described in [22]. The activity coefficients in macrostructural water have been estimated using the LLNL (Lawrence Livermore National Laboratory) model for charged species [31] and Drummond approach for uncharged species [32]. The chemical composition in microstructural water is calculated assuming Donnan equilibrium between both structural levels [18]. The activity coefficients in microstructural water were assumed equal to the macrostructural one.
The dissolution/precipitation mineral rates are defined by the general equation provided by Lasaga [33].

Numerical test description
A synthetic case consisting on a saturation test with low salinity water of a compacted MX-80 sample has been simulated in order to check the scope of the numerical tool. The test comprised a double symmetrical infiltration at the top and the bottom of the sample (length of 5 cm and diameter of 5 cm).

Geochemical model
The geochemical model is composed on 11 master species and 32 secondary species. Parameters of these species (the geochemical reactions with their equilibrium constants, hard core diameter and diffusion coefficients in free water) can be found in [14].
The accessory mineral contemplated were calcite, gypsum quartz and K-feldspar (Table 1).

Hydromechanical model
The values of the hydromechanical parameters of the constitutive equations were the same as those used for [13].

Initial and boundary conditions
The infiltration test consisted on a double symmetrical infiltration at the top and the bottom of the sample (length of 5 cm and diameter of 5 cm). Table 2 show the initial properties of the MX-80 sample. The initial conditions for the mechanical problem are shown in Table 3. Table 4 presents the chemical composition of the initial porewater and the inflow solution (top/bottom) A restriction in the displacement in all directions (sample totally confined) and isothermal conditions (25 ºC) were assumed during the test.

Results and discussion
In order to check the scope of the numerical tool, the synthetic saturation test was simulated for one year. Figure 2 presents the spatial distributions at different times of a set of characteristic variables of hydromechanical behaviour: microstructural and macrostructural void ratio, degree of saturation, liquid pressure, net mean stress and radial stress. Since the configuration of the problem can be considered symmetrical, only half of the sample is represented. Increases in both void ratios are observed at the initial stages of the saturation process at the boundaries and a progressive equilibration over time. Saturation is reached approximately after 75 days. The evolution of net mean stress and radial stress are almost parallel, pointing to a uniform distribution of axial (vertical) stress since no friction is added to the model. Figure 3 shows the spatial distribution of different variables related with the chemical behavior of the sample such as: total concentration of calcium and sulphate in micro and macrostructural water, gypsum concentration and pH. Dissolution of gypsum is produced as can be observed in Figure 3c. The dissolution of the total content of gypsum is reached after approximately 320 days. During this process the concentration of calcium increases (more significatively in microstructural water). The sulphate concentration decreases due to the low sulphate concentration in the boundaries. The pH is stabilized around in 7.7, a higher value than the saturating water. This behavior can be explained by the buffer effect produced by the carbonates/calcite system. The Figure 4 presents the time evolution several variables. The axial (vertical) stress, increases in two phases. The first phase occurs when the specimen is saturated (after approximately 75 days) and then experiences an additional increment when all the gypsum is dissolved and all the sulphate resulting from the dissolution leaves the specimen (after 320 days). Average void ratios present a symmetrical evolution, as expected due to the confined conditions and the average ionic strength evolves similarly to swelling pressure. An initial drop is reached after saturation and a further decrease can be observed after the total dissolution of gypsum. A further increase would be observed after the dissolution of the calcite, which presents a slower kinetics.

Conclusions
The present work has shown an application of a THMC model of the behavior of active clays. Specifically, a test of saturation of a sample of bentonite MX-80 has been simulated with a solution of lower ionic strength than the initial porewater. The model has been able to reproduce all the expected hydromechanical and geochemical behaviors, as well as all the couplings between the two phenomena correctly. It is of special interest to see how the axial stress curve, swelling pressure, presents a behavior in two phases, due to the coupling of the hydromechanical phenomena of saturation and geochemical dissolution of the minerals that compose the bentonite.