The Modeling of Laboratory Experiments with COMSOL Multiphysics Using Simplified Hydromechanical Model

Coupled physical processes will take place in a multibarrier disposal system for spent nuclear fuel and high-level radioactive waste. The knowledge of these processes (thermal, hydraulic, mechanical, chemical, microbiological, etc.) as well as the scope and scale of their interactions is fundamental for the safety assessment of a disposal facility. Numerical modeling is an important component in the process of acquiring and deepening the knowledge of coupled processes, while experimental evidence isimportant for model validation. This article will present a hydro-mechanical model developed by the Lithuanian Energy Institute (LEI) in the framework of H2020 project BEACON (Bentonite Mechanical Evolution). The non-linear elastic model developed in COMSOL Multiphysics (Burlington, MA 01803, USA) was applied to predict the swelling behavior of largescale oedometer experiments (MGR) performed by Research Centre for Energy, Environment and Technology (CIEMAT, Spain). In these experiments on bentonite hydration at isochoric conditions, a sample was made of two layers of calcium bentonite (FEBEX type) having initially different hydromechanical characteristics: one layer made of pellets and the other of a compacted block. Satisfactory agreement between the modeling results and the experimental data were obtained, especially for water intake and sample saturation.


Introduction
Every member state of the European Union (EU) is responsible for the implementation of a safe and sustainable solution for spent nuclear fuel and radioactive waste management and disposal, as is stated in the EU legislation (EU directive 2011/70/EURATOM). An international consensus exists on geological disposal as a mature solution for the final management step of generated and highly radioactive waste. The concept of geological disposal is based on a multibarrier system which comprises the engineered barrier system (EBS) consisting of the waste form, waste canisters, backfill, buffer, plugs, seals, and the natural barrier provided by the repository host rock. Each of the engineered barriers will serve safety functions (such as isolation, containment) due to their appropriate physical, chemical and/or mechanical properties. Thus, the knowledge and understanding of such properties and their evolution over time is crucial in the context of geological disposal.
A bentonite-based material is often considered as a possible engineered barrier in various geological repository concepts. This consideration is driven by bentonite's suitable properties such as its swelling capacity, low hydraulic conductivity under saturation conditions, cation exchange capacity for radionuclide sorption, and self-sealing ability [1]. To a large extent these properties are determined by bentonite's mineralogical composition. Despite the existence of various mineralogical compositions of bentonite, the most commonly investigated candidates for engineered barriers are calcium (Ca) and sodium (Na) bentonites. The representative examples of Ca-bentonite are FEBEX from Almeria deposits in Spain [2][3][4][5] and FoCa from France [6,7]. The examples of Na-bentonite are MX-80 from Wyoming (USA) [8][9][10] and Gaomiaozi (GMZ) from China [11,12].
Over the past few decades, a number of national and international research projects and studies have been carried out to develop the understanding of the complex behavior of bentonite-based materials. A broad range of aspects related to thermo-hydro-mechanicalchemical processes was investigated in numerous studies focusing on high temperature effect [13,14], liquid/gas flow properties [15][16][17], mechanical [18,19] or chemical [20,21] aspects. While much research has been conducted on compacted bentonite material, comparativelyfewer investigations have been performed on other forms of bentonite (granular [22], pellets of different sizes and forms [23][24][25][26], pellet-powder mixtures [6,[27][28][29]).
The research activities range from experimental investigations at different scales to numerical modeling activities. Experimental evidence at various scales gives insights about the behavior of engineered materials, but for a relatively short timeframe. Nonetheless, these pieces of evidence are very important for the development and validation of numerical models which could provide answers from a long-term perspective. Numerical models of bentonite behavior under various conditions frequently couple the Richards or multiphase flow equations with mechanical models, from simplified ones such as linear elastic [30,31] or non-linear elastic models [31][32][33] to more complex elasto-plastic models [34,35] or models based on hypoplasticity [36]. In the last few decades, elasto-plastic models, such as the Barcelona Basic Model [37] or the Barcelona Expansive Model [38][39][40], have become leaders in the simulation of bentonite behavior. However, the increasing complexity of the models has led to an increasing number of parameters, which is computationally challenging. At the same time the calibration of the models may require more extensive experimental programs compared to relatively simpler models. Thus, the demand for the development of relatively simpler models for the reasonable prediction of bentonite behavior still remains.
In 2017 the H2020 project BEACON (Bentonite Mechanical Evolution) was initiated [41]. The objective of the project is to develop the understanding of fundamental processes that lead to bentonite homogenization and to improve the capabilities of the numerical models. In the earlier studies of the EBS, the interaction between separate components was neglected and the "ideal" final state was assumed. However, the initial mechanical properties of the installed unsaturated EBS, consisting of pellets, blocks and voids, will be different after full saturation. Thus, from the repository safety perspective, it is very important to demonstrate that the performance of the present designs of backfill, buffer, plugs and seals is appropriate. As noted in [41], the analysis of the existing knowledge base in the field revealed current limitations in the predictive capability of numerical models, as the issue of swelling and homogenization is challenging from a conceptual and a numerical point of view. New experiments presenting real bentonite EBS have been performed in the project to support the development of numerical models. It is important to note that within this project several bentonite materials (MX-80, FEBEX) of different forms (block, pellets, pellets-powder mixture) have been investigated experimentally and numerically.
The Lithuanian Energy Institute (LEI) has been participating in the modeling activities of the BEACON project. A non-linear elastic hydro-mechanical model capable of representing hydro-mechanical behavior and the interaction between different forms of bentonites was developed in COMSOL Multiphysics (USA) and validated against experiments performed by project partners. Some modeling results of MX-80 type bentonite behavior under hydration are analyzed and compared with experimental results in [42]. In this article, modeling results of the hydro-mechanical behavior and the interaction of FEBEX bentonite forms (compacted block and pellets) are presented and validated against large-scale oedometer experiments (MGR) [43,44].

Methodology
This study was devoted to the numerical analysis of hydro-mechanical behavior of FEBEX type bentonite blocks and pellets at laboratory scale. Five large-scale oedometer  [43,44]. The aim of the hydration experiments was to follow evolution over time of bentonite swelling and homogenization and factors affecting them: hydration rate and sample layout. Approximately half of cylindrical cell was completed with a compacted bentonite block having an average dry density of 1.60 g/cm 3 and the other part with bentonite pellets having an average dry density close to 1.30 g/cm 3 . The sample was hydrated with deionized water through the bottom of the cell by keeping constant water pressure or constant water flux. Two sample layouts were considered: (a) pellets in the upper part and the bentonite block in the lower part (MGR21-MGR24) and (b) the bentonite block in the upper part and pellets in the lower part of the oedometer cell (MGR27). The axial pressure was always measured on the top of the cell. More details about MGR experiments can be found in [43,44].
The Richards equation was considered for the water flow modeling in the MGR experiments. It was assumed that bentonite mechanical responses in terms of deformation or/and developed swelling pressure are mainly governed by bentonite saturation. Wetting-induced swelling was modelled as elastic strain and its impact on porosity change was assessed. Reduction of stiffness during the hydration of bentonite with water occurs, thus Young's modulus dependency on saturation was considered in the model. The hydro-mechanical (HM) model included couplings to consider the impact of mechanical deformations on porosity change and subsequently on change of water balance, specific moisture capacity, storage coefficient and permeability. Plastic deformations of bentonite were not considered within this model. The non-linear elastic HM model was implemented in the numerical tool COMSOL Multiphysics. The comparison of the modeling and experimental results as well as the discussion are presented in Section 3 of this article.

Mathematical Model
During hydration the top outlet remained open to atmosphere, as is indicated in [44], air was not pressurized, thus, it was assumed that the pore space not filled with liquid contains an immobile fluid (gas) at atmospheric pressure. Subsequently, the unsaturated state could be reasonably represented by Richards equation instead of full two-phase flow model. Partial differential equations for the mass balance of liquid water were solved: where ρ is the density of fluid, C m is specific moisture capacity, g is the acceleration of gravity, S e is effective saturation, S is the storage coefficient, p is the dependent variable (pressure of liquid phase), k s is absolute permeability at saturated conditions, µ is the dynamic viscosity of fluid, k r is relative permeability, D represents elevation and Q m is fluid source (positive) or sink (negative). Within an unsaturated material, hydraulic properties change as fluids move through the medium, filling some pores and draining others. It was assumed that the pore space not filled with liquid contains an immobile fluid (gas) at atmospheric pressure. The difference between the pressure potentials of gas (p gas ) and liquid (p liq ) in unsaturated conditions is called capillary pressure (p c ) (suction): The relation between capillary pressure (suction) and saturation was represented by the van Genuchten formulation [45]: where p 0 is air entry pressure, m and n are van Genuchten fitting parameters. Different capillary pressure curves were defined for the compacted block and pellets in the model as shown in Figure 1. Parameter values for van Genuchten water retention curves fitting the experimentally measured data were reported in [46,47] and were taken into the model considering initial dry densities of the bentonite block and pellets.
where p0 is air entry pressure, m and n are van Genuchten fitting parameters. Different capillary pressure curves were defined for the compacted block and pellets in the model as shown in Figure 1. Parameter values for van Genuchten water retention curves fitting the experimentally measured data were reported in [46,47] and were taken into the model considering initial dry densities of the bentonite block and pellets. Bentonite hydration highly depends on the sample's hydraulic properties (permeability and saturation with water). Unsaturated hydraulic permeability is dependent on effective saturation Se. While measurement of saturated permeability is straightforward, direct estimation of unsaturated permeability is complicated. Usually it is estimated from back-calculation expressing unsaturated permeability as the product of relative permeability kr and saturated permeability ks and fitting the relative permeability function with measured flow data.
Saturated permeability ks for the bentonite block and pellets was derived from saturated hydraulic conductivity Ksat, which was defined according to empirical relationship as a function of dry density ρdas reported in [46]: The dependency of relative permeability on effective saturation was expressed as a power law: In the model, the exponent n=3 [46] was selected for the FEBEX bentonite block, and the exponent n=1.9 [47] for the pellets. Bentonite hydration highly depends on the sample's hydraulic properties (permeability and saturation with water). Unsaturated hydraulic permeability is dependent on effective saturation S e . While measurement of saturated permeability is straightforward, direct estimation of unsaturated permeability is complicated. Usually it is estimated from back-calculation expressing unsaturated permeability as the product of relative permeability k r and saturated permeability k s and fitting the relative permeability function with measured flow data.
Saturated permeability k s for the bentonite block and pellets was derived from saturated hydraulic conductivity K sat , which was defined according to empirical relationship as a function of dry density ρ d as reported in [46]: The dependency of relative permeability on effective saturation was expressed as a power law: In the model, the exponent n = 3 [46] was selected for the FEBEX bentonite block, and the exponent n = 1.9 [47] for the pellets.
The mechanical behavior of bentonite was represented using a pure elastic constitutive model based on generalized Hooke's law complemented with hydration (swelling)-induced strain in analogy with thermally induced strain in [23]: where σ ij is the effective stress tensor, C e ijkl are the components of the fourth-order stiffness tensor of material properties, ε kl is the total strain tensor and ε sw kl is the swelling strain tensor. This is a rather conservative representation of strains without distinguishing whether they are related to changes in microstructure or macrostructure; even existence of double structure is confirmed in compacted bentonite. The focus is the overall behavior of a system made of two bentonite forms under wetting and the ability of its representation with averaged estimates such as degree of saturation. With an assumption of reversible deformations, the material stiffness can be determined by Young's modulus E and Poisson's ratio ν. Young's modulus in both bentonite layers was assumed to be dependent on water saturation as a result of model calibration.
The bentonite-based material swells after absorption of the water. Under free swell conditions, the swelling results in the volume change, while under confined conditions, the swelling pressure builds up against the retaining cell walls. The elastic swelling model [6] was adopted to model MGR experiments, where swelling strain was modelled as proportional to the change of water saturation: where β sw is a moisture swelling coefficient for isotropic material, ∆S w is the change of water saturation, δ ij -Kronecker-delta. Taking into account the initial dry density of the pellet layer and experimentally measured swelling pressure reported in [47], the swelling pressure expected from pellet layer is relatively low (0.5-0.7 MPa) compared to that from bentonite block of 1.6 g/cm 3 (5 MPa). This suggests that the overall swelling of the sample will be mainly governed by the block. Thus, it was assumed that the pellets will swell into the void space around them during hydration, but there will be no swelling-induced stress of the pellet layer as a whole, reducing the number of model parameters to be calibrated/fitted. The swelling coefficient in the bentonite block was assumed to be dependent on water saturation as a result of model calibration.
With the assumption of the slow deformation of a solid [48], the porosity change was evaluated as a function of volumetric strain ε v (ε v = ε 11 + ε 22 + ε 33 ): The following couplings were considered in this hydro-mechanical model: • Feedback on water mass balance through water source/sink due to change of porosity [48]: • The porosity change impact on flow equation parameters: specific moisture capacity C m and storage coefficient S: where χ f is water compressibility, χ p is matrix compressibility, h is pressure head (h = p ρg ).

Numerical Model
The modeling of CIEMAT experiments (MGR21-MGR24 and MGR27 [43,44])were performed with a numerical model developed in COMSOL Multiphysics (Burlington, MA, USA). COMSOL Multiphysics is a general-purpose platform for modeling engineering applications. It allows conventional physics-based user interfaces and coupled systems of partial differential equations for simulation with the finite element method.

Geometry and Discretization
MGR experiments were performed in a constant-volume oedometer with an internal cell diameter of 100 mm and an internal height of approximately 100 mm [43,44]. The initial heights of the bentonite block layer varied slightly between 49.4 mm and 50.1 mm and the initial heights of the pellet layer between 49.7 mm and 50.4 mm, as indicated in Table 1. The numerical model was realized in 2-D axisymmetric conditions. Analyzed domains were discretized into 1502 triangular mesh elements, as presented in Figure 2. The interface between the bentonite pellets and the block was discretized more to reduce numerical errors and to obtain more accurate modeling results. Differences in the sample layout were taken into account:

•
The material representing the layer of bentonite block was assigned to the upper domain of the computational model and the bentonite pellets was assigned to the lower domain for the modeling of the MGR21-MGR24 experiments (layout Up: block/Lo: pellets).

•
For the modeling of the MGR27 experiment, the material properties were assigned in the opposite way (layout Up: pellets/Lo: block).

Numerical Model
The modeling of CIEMAT experiments (MGR21-MGR24 and MGR27 [43,44])were performed with a numerical model developed in COMSOL Multiphysics (Burlington, MA, USA). COMSOL Multiphysics is a general-purpose platform for modeling engineering applications. It allows conventional physics-based user interfaces and coupled systems of partial differential equations for simulation with the finite element method.

Geometry and Discretization
MGR experiments were performed in a constant-volume oedometer with an internal cell diameter of 100 mm and an internal height of approximately 100 mm [43,44]. The initial heights of the bentonite block layer varied slightly between 49.4 mm and 50.1 mm and the initial heights of the pellet layer between 49.7 mm and 50.4 mm, as indicated in Table 1. The numerical model was realized in 2-D axisymmetric conditions. Analyzed domains were discretized into 1502 triangular mesh elements, as presented in Figure 2. The interface between the bentonite pellets and the block was discretized more to reduce numerical errors and to obtain more accurate modeling results. Differences in the sample layout were taken into account:


The material representing the layer of bentonite block was assigned to the upper domain of the computational model and the bentonite pellets was assigned to the lower domain for the modeling of the MGR21-MGR24 experiments (layout Up: block/Lo: pellets).  For the modeling of the MGR27 experiment, the material properties were assigned in the opposite way (layout Up: pellets/Lo: block).  In all numerical models representing MGR21-MGR24, MGR27 experiments it was represented that hydration by water takes place through the bottom boundary of the computation model.

Input Parameters
The parameters necessary to model the HM behavior of bentonite pellets and block are summarized in Table 2.

Initial and Boundary Conditions
The initial conditions and duration of each MGR experiment were different, and it was reflected in the modeled cases accordingly, see Table 3. Water supply to the samples during the MGR experiments was represented by setting appropriate boundary conditions for the numerical model. As already mentioned above, the hydration took place through the bottom part of the oedometer cell. Therefore, no flow hydraulic conditions were set in the model for the top and side boundaries. For the bottom boundary, the Dirichlet boundary condition was applied (a constant water pressure of 15 kPa representing unlimited water flow rate) for the MGR21, MGR23, MGR24 and MGR27 experiments. For the model representing the MGR22 experiment, the bottom boundary was assigned the Neumann boundary condition (a limited/constant water inflow rate of 0.05 cm 3 /h). For all model boundaries, the zero-displacement condition in normal direction was assigned, i.e., friction was not taken into account and possible interaction between bentonite and the oedometer walls was not captured. As a consequence of such assumption, the same distribution of axial swelling pressure will be expected in the block and pellets.  The results of the modeling and the measurements show that a different mode of water supply through the bottom part of the sample (unlimited/predefined water flow rate) had influence on different trends of water intake evolution. The results of the models representing experiments with a constant water pressure on the hydration boundary (MGR21, MGR23, MGR24 and MGR27) showed a much faster water intake than the model with limited water supply (constant water inflow rate, experiment MGR22). According to the modeling results, half of the water volume necessary for full saturation of the sample (approximately 130 cm 3 ) was taken in around 120 days in the case of MGR22, while it took a much shorter period of time (around 4 days) for the models representing MGR21, MGR23 and MGR24 experiments and around 30 days for the model representing the MGR27 experiment.

Evolution of Water Intake over Time
The comparison between the modeling results of the MGR22 and MGR23 experiments with the longest duration (266 day and 210 days, respectively) and of the same sample layout showed that besides different evolution of the water intake over time,adifferent volume of water was also necessary to fully saturate the sample: 250 and 270 cm 3 , respectively. The experimentally measured data show the same trend. As there is a relation of the bentonite's mechanical behavior to the water content in the material, the water supply mode could be considered as an important aspect to be considered for predicting bentonite barrier behavior.
The models representing the MGR21, MGR23, MGR24 experiments ((same layout Up: block/Lo: pellets, same hydration mode (constant water pressure at sample bottom), but different duration (31 days, 210 days, 14 days)) showed a very similar behavior of the water intake. However, the experimentally measured water intake varied to some extent, especially between the results of MGR21 and MGR23. The samples of these tests were not identical; differences in the initial saturation of pellets were larger than that for the bentonite block. Thus, the experiments with identical samples would be interesting to study, as they could explain the cause of the measured differences in water intake (initial conditions of pellets and/or measurement uncertainty). In general, a satisfactory agreement between the modeling results of the MGR22 experiment and the experimentally measured data was achieved (layout Up: block/Lo: pellets, constant water inflow rate condition). The results of the modeled experiments with layout Up: block/Lo: pellets and unlimited The comparison between the modeling results of the MGR22 and MGR23 experiments with the longest duration (266 day and 210 days, respectively) and of the same sample layout showed that besides different evolution of the water intake over time, adifferent volume of water was also necessary to fully saturate the sample: 250 and 270 cm 3 , respectively. The experimentally measured data show the same trend. As there is a relation of the bentonite's mechanical behavior to the water content in the material, the water supply mode could be considered as an important aspect to be considered for predicting bentonite barrier behavior.
The models representing the MGR21, MGR23, MGR24 experiments ((same layout Up: block/Lo: pellets, same hydration mode (constant water pressure at sample bottom), but different duration (31 days, 210 days, 14 days)) showed a very similar behavior of the water intake. However, the experimentally measured water intake varied to some extent, especially between the results of MGR21 and MGR23. The samples of these tests were not identical; differences in the initial saturation of pellets were larger than that for the bentonite block. Thus, the experiments with identical samples would be interesting to study, as they could explain the cause of the measured differences in water intake (initial conditions of pellets and/or measurement uncertainty). In general, a satisfactory agreement between the modeling results of the MGR22 experiment and the experimentally measured data was achieved (layout Up: block/Lo: pellets, constant water inflow rate condition). The results of the modeled experiments with layout Up: block/Lo: pellets and unlimited water supply (constant water pressure condition) were also in line with the ex-perimental data. In these cases, a slightly faster water intake to the sample was observed in the modeling results; however, the total absorbed volume of the water in the saturated samples agreed well.
The modeling results presented in Figure 3 also allow us to compare the water intake by samples from different layouts. The models for MGR21, MGR23, MGR24 experiments (layout Up: block/Lo: pellets) provided results indicating a faster overall sample saturation than for the MGR27 experiment (layout Up: pellets/Lo: block). In the MGR27 experiment, water supply through a compacted block was modeled. The bentonite block with a density higher than that of the whole layer of pellets had a much lower permeability and that was the main reason for the slower overall water intake compared to the experiments MGR21, MGR23 and MGR24 (layout Up: block/Lo: pellets). These observations indicate the potential importance of the sample layout with regard to the hydration boundary.

Evolution of Average Saturation over Time
The evolution of average saturation over time is strongly related to the water intake curves presented in Figure 3, while the higher volume of water supplied to the sample lead to the faster average saturation and vice versa. The modeled and experimentally measured evolution of average saturation over time in the whole sample is presented in Figure 4. water supply (constant water pressure condition) were also in line with the experimental data. In these cases, a slightly faster water intake to the sample was observed in the modeling results; however, the total absorbed volume of the water in the saturated samples agreed well.
The modeling results presented in Figure 3 also allow us to compare the water intake by samples from different layouts. The models for MGR21, MGR23, MGR24 experiments (layout Up: block/Lo: pellets) provided results indicating a faster overall sample saturation than for the MGR27 experiment (layout Up: pellets/Lo: block). In the MGR27 experiment, water supply through a compacted block was modeled. The bentonite block with a density higher than that of the whole layer of pellets had a much lower permeability and that was the main reason for the slower overall water intake compared to the experiments MGR21, MGR23 and MGR24 (layout Up: block/Lo: pellets). These observations indicate the potential importance of the sample layout with regard to the hydration boundary.

Evolution of Average Saturation over Time
The evolution of average saturation over time is strongly related to the water intake curves presented in Figure 3, while the higher volume of water supplied to the sample lead to the faster average saturation and vice versa. The modeled and experimentally measured evolution of average saturation over time in the whole sample is presented in Figure 4. The modeled and experimentally measured trends correlate quite well. However, a slightly faster average saturation was predicted by a numerical model compared to the experimental data of MGR21, MGR23, MGR24 and MGR27 experiments (unlimited water supply, constant water pressure). The modeled average saturation was 0.91 by the end of the MGR21 experiment (after 34 days), while the experimental data showed a value of 0.88. The modeled value was ~3% higher than the measured one. By the end of the MGR24 experiment (after 17 days), the modeled and estimated average saturation values were 0.8 and 0.72, respectively; the difference was ~11%. The modeled time of full average saturation in the MGR23 experiment was at about t = 150 days from the beginning of hydration. Meanwhile, the experimental data showed full sample saturation after ~175 days. Within The modeled and experimentally measured trends correlate quite well. However, a slightly faster average saturation was predicted by a numerical model compared to the experimental data of MGR21, MGR23, MGR24 and MGR27 experiments (unlimited water supply, constant water pressure). The modeled average saturation was 0.91 by the end of the MGR21 experiment (after 34 days), while the experimental data showed a value of 0.88. The modeled value was~3% higher than the measured one. By the end of the MGR24 experiment (after 17 days), the modeled and estimated average saturation values were 0.8 and 0.72, respectively; the difference was~11%. The modeled time of full average saturation in the MGR23 experiment was at about t = 150 days from the beginning of hydration. Meanwhile, the experimental data showed full sample saturation after~175 days. Within the numerical model, the saturation is prescribed as a ratio of actual volumetric water content and maximum volumetric water content (it could reach the maximum value of 1.0). Meanwhile, the experimental measure-ments reported it as 1.04. A similar discrepancy was also obtained for the MGR22 experiment (model with limited water supply, constant inflow boundary condition). The time of full average saturation in the whole sample was represented well in the model for the MGR22 experiment, i.e., it was observed after 216 days of simulation compared to approximately 226 days according to the experimental data.
Despite the MGR27 experiment (layout Up: pellets/Lo: block) having the longest duration, the modeling results and the experimental data showed that the whole sample was still not fully saturated after even 278 days. The modeling results show a value 0.93 of average saturation at the end of the MGR27 experiment, while the estimated average saturation based on measurements was 0.96. The difference between the modeled and the measured average saturation was rather small (~4%). The final state of the sample of such layout (Up: pellets/Lo: block, saturation through block) was not fully saturated, in contrast to the samples with the upside-down layout (Up: block/Lo: pellets, saturation through pellets). This again indicates the importance of the sample layout with regard to the hydration boundary.

Evolution of Axial Swelling Pressure over Time
This subsection is devoted to the discussion of axial swelling pressure build-up as a result of bentonite hydration under volume confined conditions. It is important to note that the sensor for axial swelling measurement had been installed on the top of the MGR oedometer. For MGR21-MGR24 experiments, the sensor was in the contact with the compacted bentonite block, while for the MGR27 experiment, the sensor was in the contact with bentonite pellets. The modeled and measured axial pressure for MGR21-MGR24 experiments are presented in Figure 5. the numerical model, the saturation is prescribed as a ratio of actual volumetric water content and maximum volumetric water content (it could reach the maximum value of 1.0). Meanwhile, the experimental measurements reported it as 1.04. A similar discrepancy was also obtained for the MGR22 experiment (model with limited water supply, constant inflow boundary condition). The time of full average saturation in the whole sample was represented well in the model for the MGR22 experiment, i.e., it was observed after 216 days of simulation compared to approximately 226 days according to the experimental data. Despite the MGR27 experiment (layout Up: pellets/Lo: block) having the longest duration, the modeling results and the experimental data showed that the whole sample was still not fully saturated after even 278 days. The modeling results show a value 0.93 of average saturation at the end of the MGR27 experiment, while the estimated average saturation based on measurements was 0.96. The difference between the modeled and the measured average saturation was rather small (~4%). The final state of the sample of such layout (Up: pellets/Lo: block, saturation through block) was not fully saturated, in contrast to the samples with the upside-down layout (Up: block/Lo: pellets, saturation through pellets). This again indicates the importance of the sample layout with regard to the hydration boundary.

Evolution of Axial Swelling Pressure over Time
This subsection is devoted to the discussion of axial swelling pressure build-up as a result of bentonite hydration under volume confined conditions. It is important to note that the sensor for axial swelling measurement had been installed on the top of the MGR oedometer. For MGR21-MGR24 experiments, the sensor was in the contact with the compacted bentonite block, while for the MGR27 experiment, the sensor was in the contact with bentonite pellets. The modeled and measured axial pressure for MGR21-MGR24 experiments are presented in Figure 5. Many research projects have shown that the dry density and water saturation of bentonite are among the factors influencing its swelling capacity (swelling pressure). Higher density and higher saturation result in higher swelling pressure that will develop if the Many research projects have shown that the dry density and water saturation of bentonite are among the factors influencing its swelling capacity (swelling pressure). Higher density and higher saturation result in higher swelling pressure that will develop if the sample is hydrated under volume confined conditions. The modeling results show that the axial swelling pressure in the two shortest experiments (MGR21 and MGR24) reached around 2.3 MPa (after 34 days) and around 1.1 MPa (after 14 days), respectively. These model-derived values overestimated the measured values by~18% and~25%, respectively.
The corresponding final average saturation of the bentonite block for both experiments (MGR21 and MGR24) obtained by the numerical model was 0.91 and 0.78, while the average dry density was 1.50 g/cm 3 and 1.57 g/cm 3 , respectively. These results indicate that the degree of saturation had higher importance than dry density on swelling pressure in unsaturated conditions of the bentonite block. The experimentally measured axial swelling pressures were in line with the modeling results for the MGR21 and MGR24 experiments.
The modeling results of the two longest experiments, MGR22 (266 days) and MGR23 (210 days), with the same sample layout (Up: block/Lo: pellets) show that the development of axial swelling pressure on the top of the sample (in contact with the bentonite block) depends on water supply to the sample (hydraulic boundary condition for water intake). The sample with unlimited water supply (constant water pressure condition on bottom boundary, MGR23) became saturated much faster compared to sample saturation with limited water supply (constant water flow rate condition on bottom boundary, MGR22). The development of axial swelling pressure in the MGR23 experiment was faster and the peak pressure value was reached earlier. However, the modeled peak values of axial swelling pressure (around 3 MPa) were similar for both experiments.
The final modeled average dry density of the compacted block in fully saturated samples (MGR22 and MGR23) was around 1.49 g/cm 3 . Taking into account the correlation between swelling pressure and dry density obtained in small standard oedometers for FEBEX bentonite [46], the maximum swelling pressure of around 2.8 MPa is expected for such dry density. The modeled and experimentally measured peak axial pressure observed for MGR22 and MGR23 experiments correlated well with this expected maximum swelling pressure based on empirical relationship. The peak value of axial swelling pressure on the top of the sample was represented well in the model; however, it should be noted that the overall profile of the swelling pressure evolution differed from the experimentally measured one. The experimental pattern of swelling pressure development in tests MGR21-MGR24 seems to prevail that was observed for compacted FEBEX bentonite block [19] with initial peak and higher values at saturation. This increase/decrease/increase could potentially be attributed to the double structure of bentonite, as the strains induced by microstructure cause the changes in macrostructure and the overall change depends on the suction level. As the current model did not consider micro and macro structures separately, this trend was not captured.
The modeled and measured axial swelling pressure on the top of the sample (in contact with bentonite pellets) in the MGR27 experiment (Up: pellets/Lo: block) is presented in Figure 6.
For the test MGR27, where the swelling pressure was measured at pellet zone, the smooth increase was observed. Such trend could be expected for the pellets if the initial water content is higher as reported in [6]. Swelling pressure development increase/decrease/ increase trend or more smooth development for pellet powder mixture was found to be dependent on dry density, initial water content. The modeled peak axial swelling pressure in the top of the MGR27 sample was around 2.5 MPa, while the experimentally measured peak pressure was only around 1.4 MPa. It was more than two times lower than peak axial pressure on the top of the sample of the experiments with different layout (Up: block/Lo: pellets): MGR22 and MGR23. The axial swelling pressure sensor was not installed in the bottom part of MGR cells, and therefore it was not possible to compare the developed swelling pressure on the top and the bottom part of the samples. The theoretically expected swelling pressure in the bottom part (bentonite block) of the MGR27 experiment would be around 2.2 MPa, taking into account the final average density of the bentonite block (1.454 g/cm 3 ) [46] or even higher taking into account that the bentonite block was not fully saturated at the end of the MGR27 experiment. A similar large-scale oedometer experiment for MX-80 bentonite with similar sample configuration (a compacted block in the upper part and pellets in the lower part) was performed by POSIVA (Finland) within the BEACON project. The reported results [49] show that that the peak axial swelling pressure in the pellet layer was much lower compared to the bentonite block under fully saturated conditions. It is highly possible that different swelling pressure also developed in pellets and block zones in the analyzed MGR experiments. A possible explanation of the differences in swelling pressure measured at different sample locations could be friction between the sample and the oedometer cell. For the test MGR27, where the swelling pressure was measured at pellet zone, the smooth increase was observed. Such trend could be expected for the pellets if the initial water content is higher as reported in [6]. Swelling pressure development increase/decrease/increase trend or more smooth development for pellet powder mixture was found to be dependent on dry density, initial water content. The modeled peak axial swelling pressure in the top of the MGR27 sample was around 2.5 MPa, while the experimentally measured peak pressure was only around 1.4 MPa. It was more than two times lower than peak axial pressure on the top of the sample of the experiments with different layout (Up: block/Lo: pellets): MGR22 and MGR23. The axial swelling pressure sensor was not installed in the bottom part of MGR cells, and therefore it was not possible to compare the developed swelling pressure on the top and the bottom part of the samples. The theoretically expected swelling pressure in the bottom part (bentonite block) of the MGR27 experiment would be around 2.2 MPa, taking into account the final average density of the bentonite block (1.454 g/cm 3 ) [46] or even higher taking into account that the bentonite block was not fully saturated at the end of the MGR27 experiment. A similar large-scale oedometer experiment for MX-80 bentonite with similar sample configuration (a compacted block in the upper part and pellets in the lower part) was performed by POSIVA (Finland) within the BEACON project. The reported results [49] show that that the peak axial swelling pressure in the pellet layer was much lower compared to the bentonite block under fully saturated conditions. It is highly possible that different swelling pressure also developed in pellets and block zones in the analyzed MGR experiments. A possible explanation of the differences in swelling pressure measured at different sample locations could be friction between the sample and the oedometer cell.
The numerical model did not take friction into account, and therefore the modeled axial swelling pressure was uniform along the height of the sample. In such a way the predicted peak pressure in MGR27 (see Figure 6) overestimated the experimentally measured pressure on the top of the sample (in contact with pellets). Experiments with pressure sensors on top and at the bottom of the cell would be of interest to confirm such a trend for layered samples with different forms of bentonite and one side hydration. If the impact The numerical model did not take friction into account, and therefore the modeled axial swelling pressure was uniform along the height of the sample. In such a way the predicted peak pressure in MGR27 (see Figure 6) overestimated the experimentally measured pressure on the top of the sample (in contact with pellets). Experiments with pressure sensors on top and at the bottom of the cell would be of interest to confirm such a trend for layered samples with different forms of bentonite and one side hydration. If the impact of friction is significant, then further model development should consider this effect. Nevertheless, estimating the peak pressure in the whole sample is crucial for repository safety and the model was capable of representing it.

Distribution of Water Content along the Height of the Sample
The vertical line at r = 2.5 cm along the height of the sample was selected for the comparison of the modeling results and the measured data obtained after dismantling of MGR experiments. As the experimental set-up was dismantled after different durations of experiments, the results are presented in separate graphs.
The modeling results and measured data of the distribution of water content in the bentonite sample along the distance from the hydration surface are presented in Figure 7. The solid horizontal profiles in the figure indicate the initial distribution of water content.
The vertical line at r = 2.5 cm along the height of the sample was selected for the comparison of the modeling results and the measured data obtained after dismantling of MGR experiments. As the experimental set-up was dismantled after different durations of experiments, the results are presented in separate graphs.
The modeling results and measured data of the distribution of water content in the bentonite sample along the distance from the hydration surface are presented in Figure 7. The solid horizontal profiles in the figure indicate the initial distribution of water content. As has been noted, the MGR samples were hydrated through the bottom part during all the experiments. The saturation along the height of the samples after the two shortest experiments (MGR21 and MGR24) was quite different, and therefore the modeled values of water content vary in a wide range (between 0.17 and 0.38). For the experiments where the sample was fully saturated by the end (MGR22 and MGR23) or almost saturated (MGR27), the modeled distribution of the water content was observed to be distributed much more homogeneously (varying between 0.3 and 0.35). Although the sample was finally saturated (MGR22 and MGR23 experiments), the water content along it was not fully homogenized, i.e., it slightly decreased upwards from the hydration surface (bottom of sample). The steep change at the interface of the layers was observed mainly for tests As has been noted, the MGR samples were hydrated through the bottom part during all the experiments. The saturation along the height of the samples after the two shortest experiments (MGR21 and MGR24) was quite different, and therefore the modeled values of water content vary in a wide range (between 0.17 and 0.38). For the experiments where the sample was fully saturated by the end (MGR22 and MGR23) or almost saturated (MGR27), the modeled distribution of the water content was observed to be distributed much more homogeneously (varying between 0.3 and 0.35). Although the sample was finally saturated (MGR22 and MGR23 experiments), the water content along it was not fully homogenized, i.e., it slightly decreased upwards from the hydration surface (bottom of sample). The steep change at the interface of the layers was observed mainly for tests MGR22 and MGR23. It seems to the authors that this could be attributed to the assumption of elastic behavior of both materials. Further model development should proceed to represent irreversible strains, primarily in the pellet zone. The modeling results agree quite well with the experimental data in saturated or almost saturated samples (MGR22, MGR23 and MGR27 experiments), as demonstrated in Figure 7. However, the numerical model overestimated the water content (mainly in the bentonite block) in the case of unsaturated samples (MGR21 and MGR24).
The modeling results presented in Figure 7b,c can be used to analyze the impact of the water supply mode on the water content distribution. As the figure shows, the measured water content in the block layer was quite similar, while it differed to some extent for the pellet layer (MGR22 and MGR23 experiments). The numerical model results show a more similar distribution for both experiments: with unlimited water supply (MGR23) and limited water supply (MGR22). However, the initial water content of the pellets was not identical for these experiments, and thus it was not possible to reasonably distinguish whether the hydration mode had an influence on the final water content distribution or whether it was somewhat hidden by the differences in the initial conditions.
A comparison of the water content distribution within samples with different layouts (MGR23, layout Up: block/Lo: pellets and MGR27, layout Up: pellets/Lo: block) could be produced from the results presented in Figure 7c,d. As the results show, the hydration through the bottom of the sample with layout Up: pellets/Lo: block led to a slightly more even distribution of water content along the sample height in comparison to the sample with hydration through the pellet layer. The initial water content in the samples was identical, which means that the impact of the sample layout had a limited effect on the final water content vertical distribution along the sample.

Distribution of Dry Density along the Height of the Sample
The modeling results and measured data of the distribution of dry density in the bentonite sample along the distance from the hydration surface are presented in Figure 8. The solid horizontal profiles in the figure represent the initial distribution of dry density. saturated samples (MGR21 and MGR24).
The modeling results presented in Figure 7b,c can be used to analyze the impact of the water supply mode on the water content distribution. As the figure shows, the measured water content in the block layer was quite similar, while it differed to some extent for the pellet layer (MGR22 and MGR23 experiments). The numerical model results show a more similar distribution for both experiments: with unlimited water supply (MGR23) and limited water supply (MGR22). However, the initial water content of the pellets was not identical for these experiments, and thus it was not possible to reasonably distinguish whether the hydration mode had an influence on the final water content distribution or whether it was somewhat hidden by the differences in the initial conditions.
A comparison of the water content distribution within samples with different layouts (MGR23, layout Up: block/Lo: pellets and MGR27, layout Up: pellets/Lo: block) could be produced from the results presented in Figure 7c,d. As the results show, the hydration through the bottom of the sample with layout Up: pellets/Lo: block led to a slightly more even distribution of water content along the sample height in comparison to the sample with hydration through the pellet layer. The initial water content in the samples was identical, which means that the impact of the sample layout had a limited effect on the final water content vertical distribution along the sample.

Distribution of Dry Density along the Height of the Sample
The modeling results and measured data of the distribution of dry density in the bentonite sample along the distance from the hydration surface are presented in Figure 8. The solid horizontal profiles in the figure represent the initial distribution of dry density. The modeled dry density distribution along the height of the samples from the shortest experiments (MGR21 and MGR24) is presented in Figure 8a. As the figure shows, the modeled dry density distribution along the height varied in a wide range. Besides, the modeled values were significantly changed compared to the initial ones, especially in the compacted bentonite block layer. In the saturated samples (MGR22 and MGR23), the modeled dry density was more homogenous (varying between 1.37 and 1.49 g/cm 3 ). The modeled dry density distribution in the MGR27 experiment was the most homogeneous (varying between 1.4 and 1.47 g/cm 3 ) among all analyzed cases.
The modeling results varied to some extent compared to the experimentally measured values. The best correlation of dry density distribution was obtained in the shortest experiment MGR24 (14 days) (see Figure 8a). In an experiment of longer duration (MGR21, 34 days), the agreement was not so good, especially for the bentonite pellet layer (see Figure 8a). Potentially, it is related to assumed hydro-mechanical behavior and parameter values of the bentonite pellets. For experiments with finally fully saturated samples (MGR22 and MGR23), two uniform different distributions of dry densities were modeled in the block layer and the pellet layer, while the experimental data showed a gradual increase in dry density from the hydration surface upwards. Nevertheless, the modeled final values of dry density were similar to the experimental data, especially for the bentonite pellets. For the MGR27 experiment, almost uniform distributions of dry density were predicted by the numerical model in the block and pellets, which means underestimating the dry density in the pellet layer and overestimating it in the bentonite block layer.
The modeling results presented in Figure 8b,c can be used to analyze the impact of the water supply mode on the water content distribution. As the figure demonstrates, the measured dry density in the block layer was quite similar, while it differed to some extent for the pellet layer (MGR22 and MGR23 experiments). The numerical model results show a more similar distribution for both experiments: with unlimited water supply (MGR23) and limited water supply (MGR22). However, the initial water content of the pellets was not identical for these experiments, and thus it was not possible to distinguish whether the hydration mode had an effect on the final dry density distribution or whether it was somewhat hidden by the differences of initial conditions.
A comparison of the modeled water content distribution within samples with different layouts (MGR23, layout Up: block/Lo: pellets and MGR27, layout Up: pellets/Lo: block) could be produced from the results presented in Figure 8c,d. The modeling of the hydration through the bottom of the sample in both layouts provided quite similar distributions of dry density along the sample height. The initial water content in the samples was identical, and thus the numerical model did not predict significant differences in dry density distribution for both analyzed sample layouts.
The first attempts were made by the LEI team to model the system of two different forms of bentonite with a simplified model. Taking into account the obtained results, it can be seen that model's predictive capabilities are limited. The model output could be treated more as indicatory of trends (e.g., full saturation time under same hydration conditions, tendency of occurrence/absent of homogenization) but not the absolute values. Subsequently, further model developments are needed with the main focus of the consideration of friction (for laboratory scale experiments), the representation of irreversible strains and the role of the microstructure on the overall behavior of different forms of bentonite: capability to represent water retention dependency on dry density should be added too. The new experimental findings about pellets and pellet-powder mixtures from the BEACON project are expected to support further model development. Different material layouts, hydration conditions should be explored further too.

Conclusions
A non-linear elastic hydro-mechanical model was developed in COMSOL Multiphysics (Burlington, MA, USA) in the framework of H2020 project BEACON (Bentonite Mechanical Evolution). This model was applied to predict the swelling behavior of five large-scale oedometer experiments (MGR) performed by Research Centre for Energy, Environment and Technology (CIEMAT, Madrid, Spain), where two layers of FEBEX bentonite (pellets and a compacted block) having initially different hydro-mechanical characteristics were hydrated at isochoric conditions. The following conclusions can be drawn from the comparison of the obtained modeling results and the experimental data:

•
In general, a satisfactory agreement was obtained between the modeling results and the experimental data, especially for water intake to the samples and average water saturation in all analyzed MGR experiments; • Although the modeled evolution of axial swelling pressure on the top of the bentonite block differed from that experimentally measured in MGR22 and MGR23, the model did a good representation of the peak value of around 3 MPa; • The development of axial swelling pressure was faster and the peak pressure value was reached earlier for the experiments with unlimited water supply to the samples. The modeling results of the MGR22 and MGR23 experiments revealed that the hydration strategy in terms of limited/unlimited water supply had a significant impact on the axial pressure evolution but not on the final value; • The numerical model predicted the same axial pressure profiles at the bottom and at the top of the samples as wall friction, which was not taken into account during modeling. However, the peak axial swelling pressure in the whole sample was represented well, which is crucial for repository safety; • The modeled distributions of water content and dry density along the bentonite sample from the hydration surface varied to some extent compared to the experimentally measured values. However, the trend of the distributions and the final values of these parameters were in line with the experimental data.

•
Taking into account the obtained results, it could be concluded that the model could be used for similar analysis in the future to some extent. However, further more precise investigations into the behavior of bentonite material in different forms, layouts and hydration conditions and subsequent model development are needed.