Numerical modelling of unsaturated MX-80 bentonite subjected to two different hydration paths and subsequent loading to high-pressures

. MX-80 bentonite has been considered as a suitable material for the construction of engineered barriers employed in deep geological radioactive waste repositories. These barriers are generally formed of compacted unsaturated bentonite, the latter experiencing a slow saturation due to its low permeability while interacting with the surrounding groundwater. In order to verify the long-term safety requirements of engineered barriers, their response to hydration has to be carefully assessed. As part of the recent European project BEACON (Bentonite Mechanical Evolution), the behaviour of MX-80 bentonite subjected to different hydration paths was investigated in a number of laboratory and field experiments and numerical studies. This paper is concerned with numerical simulations of two laboratory experiments performed during the project, with the objective of examining the predictive capabilities of the proposed numerical modelling approach. The experiments were selected due to the granular state of bentonite at its placement in the testing apparatus, which differed from the large number of previous experiments conducted on specimens of compacted bentonite blocks. The paper provides a brief introduction to the adopted modelling framework, a summary of calibrated parameters for the hydro-mechanical constitutive modelling and the results of numerical simulations, concluding that a satisfactory numerical simulation of the experiments was achieved.


Introduction
MX-80 bentonite has been examined over the years as a suitable material for the construction of engineered barriers employed in geological radioactive waste disposal applications [1 to 5]. These barriers are generally formed of unsaturated bentonite, the latter experiencing a slow saturation as a consequence of its interaction with the surrounding groundwater. In order to verify the long-term safety requirements of engineered barriers, their response to hydration is one of the key parameters to be carefully assessed.
As part of the European project BEACON (Bentonite Mechanical Evolution), the behaviour of MX-80 bentonite was investigated in a number of experimental laboratory settings, focused on characterising the dual-porosity structure of this expansive clay. One set of such tests, applying different hydration paths, was conducted at EPFL (École Polytechnique Fédérale de Lausanne) [6]. An unsaturated granular MX-80 bentonite was first poured into a high-pressure oedometer cell and subsequently subjected to either free-swelling or constant-volume hydration starting from hygroscopic conditions, before proceeding with the application of vertical stresses as high as 20 MPa. The results of these laboratory tests were subsequently used in the * Corresponding author: giuseppe.pedone@unitn.it BEACON's numerical modelling work package. The aim was to assess the predictive capabilities of different constitutive models by reproducing the hydromechanical processes investigated in the laboratory and comparing the simulation results with experimental observations. The present paper reports the results of twodimensional (2D) axisymmetric finite element (FE) analyses carried out to reproduce two oedometer tests performed at EPFL. The simulations were conducted with the Imperial College Finite Element Programme (ICFEP) [7,8], adopting the Imperial College dualporosity structure model IC-DSM [9,10] in combination with (i) a van Genuchten-type soil water retention (SWR) model that accounts for specific volume variations and (ii) a hydraulic conductivity (HC) model that assumes logarithmic variation of permeability with suction.
A description of the models used in the numerical simulations is first reported. The FE meshes generated for the analyses are subsequently introduced, together with model parameters and the boundary conditions adopted. A comparison between experimental and numerical results is then shown and conclusions are finally drawn.

Description of the models 2.1 Mechanical model (IC-DSM)
The Imperial College dual-porosity structure model (IC-DSM) [9,10] is an extension of the IC single structure model (IC-SSM) [11,12], that was developed based on the Barcelona Basic Modelling (BBM) framework [13]. The model is formulated for unsaturated clays and adopts two independent stress variables: (i) suction, ‫ݏ‬ = ‫ݑ‬ െ ‫ݑ‬ ௪ , and (ii) net stress, ߪ ത = ߪ ௧௧ െ ‫ݑ‬ , with ‫ݑ‬ and ‫ݑ‬ ௪ being the air and water pressures in the pores, respectively, and ߪ ௧௧ being the total stress.
To allow for a smooth transition from saturated to unsaturated states and vice versa, the model also introduces an equivalent suction, ‫ݏ‬ = ‫ݏ‬ െ ‫ݏ‬ , and an equivalent stress, ߪ = ߪ ത + ‫ݏ‬ , where ‫ݏ‬ is the air-entry value of suction. The model is further generalised in the ( ‫,ܬ‬ ‫,‬ ߠ, ‫ݏ‬ ) space, where ‫,ܬ‬ ‫,‬ ߠ are the invariants of the equivalent stress tensor, representing generalised deviatoric stress, mean equivalent stress and Lode's angle, respectively.
In order to enable the modelling of unsaturated expansive clays, a double-porosity structure formulation has been implemented in the model, in agreement with [14,15]. This formulation differentiates two levels of structure in the clay: the macro-structure, whose mechanical behaviour is ruled by the original IC-SSM formulation; and the micro-structure, assumed to be elastic, volumetric and fully saturated.

Behaviour of the micro-structure
Assuming the micro-structure to be fully saturated implies that it can be defined in terms of effective stresses, where the mean effective stress corresponds to ‫‬ ᇱ = ‫‬ + ‫ݏ‬ . The assumptions that the micro-structure is also volumetric and elastic imply that changes in ‫‬ ᇱ result in elastic volumetric micro-strains, οߝ ௩, : where the micro-structural bulk modulus, ‫ܭ‬ , is defined as: The bulk modulus ‫ܭ‬ is additional to the two bulk moduli associated with the macro-structure and defined in the IC-SSM formulation: ‫ܭ‬ ௦,ெ , associated with equivalent suction, and ‫ܭ‬ ,ெ , associated with mean equivalent stress. The three bulk moduli define the overall elastic soil behaviour in the double-structure formulation.
In Equation (2), ݁ is the micro-structural void ratio and ߢ is the micro-structural elastic compressibility parameter. For consistency, the following must be satisfied: where ݁ ெ is the macro-structural void ratio and ݁ is the overall void ratio of the material.
The IC-DSM model also introduces a void factor, ‫ܨܸ‬ = ݁ /݁, that expresses the degree of dominance of individual structural levels in the overall clay structure, hence enabling the quantification of the micro-structural evolution in the clay. ‫ܨܸ‬ is automatically set to zero when fully saturated conditions are attained, i.e. when ‫ݏ‬ = 0.

Interaction between micro-structure and macro-structure
Although the micro-structural volumetric deformation is elastic, it is assumed to contribute to the macrostructural volumetric plastic strains, οߝ ௩,ఉ , through an additional plastic mechanism: defined by the interaction function, ݂ ఉ , between the two levels of structure. The shape of this function is dependent on whether the micro-structure swells (S) or compresses (C) and is defined as: In Equation (5): ‫‬ ‫/‬ expresses the degree of openness of the structure in terms of the distance between the current stress state (represented by ‫‬ ) and the yield surface (represented by ‫‬ ), while ܿ ଵ , ܿ ଶ , ܿ ଷ and ܿ ௦ଵ , ܿ ௦ଶ , ܿ ௦ଷ are coefficients defining the shape of the interaction function.

Soil water retention (SWR) model
For the FE analyses herein presented, a non-hysteretic Van Genuchten-type [16] SWR model was adopted, formulated in terms of degree of saturation, ܵ , and equivalent suction [17]: where: ܵ is the residual degree of saturation; ߙ, ݉ and ݊ are fitting parameters controlling the shape of the water retention curve; ߰ is the parameter controlling the effect of the specific volume, ߥ.

Geometry and discretisation
Two oedometer tests conducted by Bosch et al. [6] on unsaturated MX-80 bentonite were analysed by means of FE simulations performed with the code ICFEP [7,8], with the objective of assessing predictive capabilities of the employed mechanical and hydraulic constitutive modelling introduced in Section 2. This modelling approach has been applied successfully in the simulations of experiments using unsaturated compacted bentonite blocks both in laboratory experiments [19] and in field prototypes [20], involving MX-80 and Febex bentonite, respectively. The experiments simulated here were instead conducted on unsaturated specimens prepared from granular MX-80 bentonite, which was poured into a high-pressure oedometer cell. After their preparation, the specimens were subjected to: -Stage 1: either free-swelling (test P1-3) or constantvolume hydration (test P2-1), starting from hygroscopic conditions; -Stage 2: gradual application of around 20 MPa vertical stress.
Due to the axisymmetric nature of the experiments under investigation, 2D axisymmetric FE simulations were undertaken. The domains analysed (diameter of 17.5 mm, height of 12.5 mm) were discretised using 8-noded quadrilateral displacement-based elements, with 4 pore pressure degrees of freedom at the corner nodes. The meshes generated are shown in Figures 1 and  2 for tests P1-3 and P2-1, respectively, together with the mechanical boundary conditions adopted in the two stages of each test.  The numerical simulations undertaken were hydromechanically fully coupled. Moreover, given that very large deformations (around 80% volumetric strains) were expected in one of the two tests analysed, a large displacement formulation was adopted.

Input parameters
The model parameters used in both analyses are summarised in Tables 1 and 2 (the double-structure formulation parameters are reported in italic in Table 1).  All parameters were derived from the laboratory data reported in [1][2][3][4], except for the values listed in bold in Table 1. The latter were directly based on the EPFL data reported by Bosch et al. [6], therefore selected to maximise the model predictive capabilities, as discussed in Section 6.
It should be noted that both hydro-mechanically coupled and fully drained simulations of the two experiments were performed, however, only the latter analyses are reported here. Consequently, the HC model parameters listed in Table 2 did not influence the results of the analyses discussed in the following.

Initial and boundary conditions
The initial conditions assumed in both analyses are reported in Table 3. They correspond to the as-poured and hygroscopic conditions provided in [6]. Following the testing procedure adopted at EPFL, the two FE analyses conducted were also divided in two stages (see Figures 1 and 2 for a summary of the mechanical boundary conditions applied): -Stage 1: during which hydration was simulated, either under a constant vertical stress of 21 kPa (Test P1-3) or under constant-volume conditions (Test P2-1). Hydration up to full saturation was simulated by imposing homogeneous and gradual suction reductions across the entire mesh. The suction reductions were assumed to take place under fully drained conditions, and, therefore, in these analyses the evolution of total stresses, void ratios, suctions and degrees of saturation with time was not investigated. -Stage 2: during which vertical total stresses were increased up to 20 MPa in the Test P1-3 analysis and 5 MPa in the Test P2-1 analysis (the latter stopped earlier because very small void ratios were attained).

Results
A comparison between laboratory data [6] and numerical results for Tests P1-3 and P2-1 is shown in Figure 3 in terms of variations of void ratio with total stress. Additional results obtained with the FE analyses are illustrated in Figures 4 and 5 in terms of suction and degree of saturation variations, respectively.

Test P1-3
The swelling strain predicted in Test P1-3 during Stage 1 is very similar to the one observed in the laboratory (see path A-B in Figure 3). Such a large swelling strain was obtained by adopting a micro-structural compressibility, ߢ , equal to 0.360. The value adopted here is about 3 times larger than the value that could be inferred from the laboratory data shown by Tang & Cui [3]. However, the micro-structural model parameters are generally difficult to experimentally evaluate with certainty, so they are often back calculated.
At the end of Stage 1, very small positive pore pressures (2 kPa) are reached (A-B in Figure 4), the material becomes fully saturated (A-B in Figure 5), and, therefore, the void factor, ‫,ܨܸ‬ automatically switches to zero [9,10]. As a consequence, during Stage 2, the double-porosity structure is not active and the hydromechanical behaviour follows the single structure framework [11].
The pre-yield response in Stage 2 is well reproduced (Figure 3), even though a behaviour stiffer than the one measured is predicted close to yielding, but this results from the use of a constant compressibility value, ߢ. After yielding, as expected, the numerical results show a linear compressibility response (in the logarithmic scale), while the laboratory data indicate a non-linear behaviour post-yield. The model could reproduce this behaviour if unsaturated conditions were retained in the specimen, as a double-porosity structure would remain active.
The slope of the virgin compression line at full saturation used in the analyses (Table 1) is different from the one suggested by [3] and was selected to match the average post-yield compressibility observed during Test P1-3. The selected virgin compression line allowed to obtain a final void ratio very close to the one measured at the end of the test (Figure 3).

Test P2-1
During Stage 1, hydration up to full saturation (path A-B' in Figures 4 and 5) takes place under confined conditions, so the void ratio remains constant and equal to the initial value of 0.85 (path A-B' in Figure 3). A significant increase in vertical total stresses is predicted during this stage, in accordance with the laboratory data provided by Bosch et al. [6].   As shown in Figure 3, the analysis tends to underpredict the swelling pressure at the end of Stage 1. The measured value is close to 3 MPa, while the predicted one is around 2.3 MPa. However, it is worth recalling that in Test P2-1, before hydration, Bosch et al. [6] applied an additional vertical stress in the range of 0.12 to 0.31 MPa, with the aim of ensuring good contact between sample and top piston. This additional loading was not modelled in the FE analysis in order to maintain constant-volume conditions, so a swelling pressure under-estimation in the range of 0.12 to 0.31 MPa was expected.
During Stage 1, the swelling pressure exceeds the yield stress, hence, at the beginning of Stage 2 (point B' in Figure 3), the stress state lies already on the virgin compression line of the material. As a consequence, the compressibility behaviour predicted in Stage 2 of Test P2-1 corresponds to the one predicted, post-yield, during Stage 2 of Test P1-3 ( Figure 3).

Conclusions
The paper describes the FE analyses carried out to reproduce the oedometer tests performed at EPFL as part of the BEACON project. A description of the analyses is provided, together with a summary of the constitutive modelling approach used in the simulations. A comparison between experimental and numerical results is also reported, showing that the adopted modelling approach can reproduce reasonably well the hydro-mechanical behaviour of the granular MX80 bentonite when subjected to different hydration processes (confined and unconfined) and subsequent loading.
The results suggest that, when the parameter values listed in Tables 1 and 2 are used, the swelling strain in oedometric conditions is well predicted (Stage 1, Test P1-3). However, when the same set of parameter values is adopted, the swelling pressure developing in confined conditions is under-predicted by around 0.5 MPa (Stage 1, Test P2-1). Swelling pressure predictions could be refined by enhancing the coupling effects between micro-structure and macro-structure (e.g. by modifying the coupling functions regulating the ߚ interaction mechanism summarised in Section 2.1.2).
After hydration, when loading takes place (Stage 2), both analyses show a linear compressibility response (in the logarithmic scale), as the double-porosity structure ceases to be active at the end of Stage 1. The predicted compressibility behaviour is representative of the material response, but it appears to be simpler than the one observed in the laboratory, where a non-linear variation of the void ratio with the logarithm of the vertical stress was observed. The predictive capabilities of the model, in this respect, could be further improved by maintaining a multi-porosity structure also in fully saturated conditions, as suggested by the Mercury Intrusion Porosimetry (MIP) data reported in Bosch et al. [6].