THM-modelling benchmark initiative on the effects of temperature on the disposal of heat-generating radioactive waste in clay formations

RESEARCH PAPER THM-modelling benchmark initiative on the effects of temperature on the disposal of heat-generating radioactive waste in clay formations Eric Simo 8,1 •Christophe de Lesquen 2 •Rocio Paola Leon-Vargas 1 •Minh-ngoc Vu 2 •Simon Raude 3 • Ginger El Tabbal 3 •Arnaud Dizier 4 •Suresh Seetharam 4 •Asta Narkuniene 5 •Fre ´de ´ric Collin 6 • Hangbiao Song 6 •Antonio Gens 7 •Fei Song 7 •Alexandru-Bogdan Tatomir 10 •Thomas Nagel 8,9 • Jo ¨rg Buchwald 8,9 Received: 22 December 2023 / Accepted: 1 December 2024 The Author(s) 2025 Abstract Understanding the thermo-hydro-mechanical (THM) behaviour of clay formations, which are being considered as potential hosts for the disposal of radioactive waste in Europe, is important for the feasibility and the safe disposal of heat-generating radioactive waste in deep geological clayey formations. To ensure the reliability of numerical models and tools used to predict the THM evolution of such systems over long periods of time, it is necessary to verify and validate these tools through benchmarking initiatives. As part of the joint European Programme on Radioactive Waste Management (EURAD), the Influence of Temperature on Clay-based Material Behaviour (HITEC) work package invited seven teams from across Europe to participate in a benchmark initiative to assess the expertise and capabilities of these teams and their numerical tools to predict THM evolution in different clay materials. The results of this benchmark showed that all teams involved were able to adequately model the THM behaviour of heat-generating repository systems in clay formations, and the tools they used can be considered verified for the solution of coupled THM equations under variable boundary conditions relevant for nuclear waste disposal. Keywords Boom clay Callovo–Oxfordian clay Geological disposal High-level radioactive waste Opalinus Clay  Safety assessment Thermo-hydro-mechanical modelling 1 Introduction Argillaceous formations are considered across Europe as a potential host for the geological disposal of radioactive waste. This is due to their self-sealing properties, their low permeability, their small molecular diffusion and their retention capacity for radionuclides. Thermo-hydro-mechanical processes (THM processes) play an important role in the evolution of a repository in clay formations. The excavation of the repository mine leads to a damaged zone in the surrounding host geological formation. This excavation-damaged zone is characterized by a higher hydraulic conductivity that enhances hydraulic flow within the repository. Moreover, the heat generated by the waste leads to thermally induced stresses in the rock and potentially further increases the excavation-damaged zone. The heat propagation also gives rise to excess pore pressures in the host clayey formation due to the dilation of water in the porous host medium that may lead to tensile damage even in the far field of the repository system. A deep understanding of the THM behaviour of potential host clayey formations is therefore important for the design and safety assessment of repository systems in these formations. The assessment of the long-term safety of repository systems requires the use of numerical models and tools capable of reproducing the THM behaviour of the clay host Eric Simo, Christophe de Lesquen, Rocio Paola LeonVargas, Minh-ngoc Vu, Simon Raude, Ginger El Tabbal, Arnaud Dizier, Suresh Seetharam, Asta Narkuniene, Fre ´de ´ric Collin, Hangbiao Song, Antonio Gens, Fei Song, AlexandruBogdan Tatomir, Thomas Nagel, Jo ¨rg Buchwald have contributed equally to this work. Extended author information available on the last page of the article 123 Acta Geotechnica https://doi.org/10.1007/s11440-024-02502-w(0123456789().,-volV)(0123456789().,-volV) rock that has been observed through experimental investigations and predicting the THM evolution of the repository systems over the assessment period of up to one million years. It follows that the reliability of such predictive tools should be verified and validated [14]. Verification is the initial step towards validation, which will subsequently involve the comparison with in situ data from field experiments-endeavours. By verification, one means the mathematical correctness of a given numerical model. To verify such models, one can, for example, compare the output of numerical simulations against their analytic counterparts for a given problem. Given the increased complexity of numerical models employed in the realm of THM coupled safety assessment of deep geological repositories in clay, there are little analytic solutions against which computer codes can be verified. Ensuring the validity of simulation results and the subsequent reliability of decision-making become therefore increasingly challenging. This predicament is further compounded when multiple numerical codes are employed, adhering to the two-person rule during numerical analyses. Regrettably, no comprehensive guidelines or standards have been published to address the quality of research and commercial codes currently utilized in the safety assessment of repositories worldwide. Several benchmarking initiatives have been carried out internationally in various research projects to address verification by code comparison on specific topics relevant to the safety assessment of repository, see DECOVALEX [1–3], BENVASIM [13,16], Reactive Transport Model Benchmarking [10,12]. In the scope of the project HITEC as part of the European joint programme on radioactive waste (EURAD), seven modelling teams from across Europe were invited to participate in a benchmark initiative to assess the expertise and capabilities of those teams and their numerical tools to predict the THM evolution in different clay host rocks in the near field and the far field of the disposal zone. They present realistic numerical benchmarks of typical safety assessment problems aiming at demonstrating that the different codes commonly employed in Europe yield consistent results. Consequently, based on the findings, these codes can be deemed verified within the framework of ‘‘Validation & Verification’’ (V&V). The results of this benchmarking initiative are the subject of the present paper. 2 Presentation of the work package HITEC of EURAD HITEC is the acronym for ‘‘Influence of Temperature on Clay-based Material Behaviour’’ that is the seventh work package (WP) of the EURAD Project. The EURAD project is a European joint programme on radioactive waste management that helps to develop a robust and sustained science, technology and knowledge management programme that supports the timely implementation of radioactive waste management (RWM) activities of the EU member states and serves to foster mutual understanding and trust between Joint Programme participants [8]. The work package HITEC of EURAD aims to develop and document an improved thermo-hydro-mechanical (THM) understanding of clay-based materials (host rocks and buffers) exposed at high temperatures ([100 C) or having experienced high-temperature transients for extended durations. The WP’s raison d’e ˆtre is to evaluate whether or not elevated temperature limits (of 100 150 C) are feasible for a variety of geological disposal concepts for heatgenerating radioactive waste. HITEC studies clay host rock formations exposed to temperatures of up to 120 C, documents and establishes the possible extent of elevated temperature damage in the near or far field and also indicates the likely consequences of any such damage. The WP also looks at bentonite buffers and determines the temperature influence on buffer swelling pressure, hydraulic conductivity, erosion or transport properties to see when the buffer safety functions start to be unacceptably impaired [8]. For the disposal of heat-generating radioactive waste, it is important to understand the consequences of the heat produced on the properties of the natural and engineered clay barriers and on their long-term performance. Most safety cases for disposal concepts that involve clay currently consider a temperature limit of 90 100 C. Being able to tolerate higher temperatures, whilst still ensuring an appropriate performance, would have significant advantages (e.g. shorter cooling times at surface intermediate storage facilities, more efficient packaging, fewer disposal containers, fewer transport operations, smaller facility footprints, etc.). This WP interrogates the validity of the currently applied thermal limits and also the importance of the accuracy of the assumed radiological waste properties and consequently marks a first step towards optimization of the architecture of deep geological disposal facilities [8]. Laboratory experiments carried out in the scope of WP HITEC are studying the effects of increased temperature in the near field, such as fracturing and self-sealing in the excavated-damaged zone (subtask 2.1), while other dedicated experiments are looking at the far-field effects, where thermal loading and the generation of overpressures may result in the loss of integrity of the geological barrier (subtask 2.2). Subtask 2.3 is focusing on the modelling of these effects. Two sets of 2D models were created to study the impact of temperature on the behaviour of the three clayey host formations (Boom Clay, Callovo–Oxfordian (COx) claystone and Opalinus Clay (OPA)) in the near Acta Geotechnica 123 field and in the far field. Some of the laboratory experiments were also simulated to improve the understanding of their THM behaviour. Three in situ heating tests (PRACLAY in the Boom Clay, FE in the Opalinus Clay and ALC1605 in the COx) were finally modelled, taking into account the results from the first two steps. The present work describes the results of 2D benchmark studies that have been realized to test the numerical capabilities available across Europe to model the THM behaviour of the three clay host rocks [8]. 3 Benchmark description The benchmarking exercise selected for this work consists of modelling the near-field effects of excavation and heating around the disposal gallery in a repository for radioactive waste in a clay formation. The disposal concept chosen here is inspired by the Swiss and French repository concepts. In the Swiss design, the gallery will be backfilled with bentonite, whereas the French one leaves some room to facilitate a possible retrieval of the waste packages. In addition, the highly exothermic packages are separated by spacing buffers to reduce the thermal load. In order to focus on the host rocks and to be able to compare the models, the boundary condition around the gallery is set at the interface between the gallery and the host rock [8]. Three subcases are proposed [8]: •Isotropic stress conditions with isotropic thermoelasticity •Anisotropic stress conditions with cross-anisotropic (i.e. transversely isotropic) thermo-elasticity •Anisotropic stress conditions with elastoplastic/damage models (the choice of the model is left to the modelling teams) The present work will focus on the first two subcases whereas the third case will be the subject of further publications. The goal with the first two subcases is first to compare the numerical codes on a fixed exercise by tightly imposing the boundary conditions and the mechanical constitutive laws and properties. The first subcase refers to a simple elastic isotropic case, which considers isotropic hydraulic, mechanical and thermal host rock properties as well as an isotropic far-field stress field. The fully coupled thermal-hydro-mechanical (THM) elastic behaviour of Callovo–Oxfordian claystone, Boom Clay and Opalinus Clay is analysed. The isotropic host rock properties are collected in Table 1. The second subcase takes into account the anisotropy of the host rocks’ properties and the in situ stress field. The principal directions are oriented parallel and perpendicular to the bedding plane that is assumed, for all three host rocks, coincident with the horizontal and vertical directions, respectively. Along the vertical principal direction, all host rocks exhibit lower intrinsic permeability, elastic modulus and thermal properties. The anisotropic host rock properties are collected in Table 2. Solid and water phase properties are the same as in the isotropic elastic case. As already mentioned, the model consists of a cross section of a heating gallery and host rock perpendicular to the gallery axis. Due to symmetry, only a quarter of the gallery cross section is modelled. The simulated region measures 100 m in both xand ydirections and a plane strain condition (zfx;y;zg¼0) is adopted. A representation of the domain is shown in Fig. 1together with the specified Table 1 Isotropic thermo-hydro-mechanical parameters for the Callovo–Oxfordian claystone, Opalinus and Boom Clay [8] Parameters Units COx OPA Boom Clay Solid phase density qs kg m32690 2340 2639 Saturated bulk density q kg m32386 2030 2000 Porosity n– 0.18 0.13 0.39 Intrinsic perm. Km22:30 1020 3:01020 2:83 1019 Young’s modulus E MPa 7000 6000 300 Poisson’s ratio m– 0.3 0.3 0.125 Biot coefficient b– 0.8 0.6 1 Thermal conductivity k Wm1 K1 1.67 1.85 1.47 Linear thermal expansion coefficient as K11:25 1051:71051:0105 Solid specific heat capacity cp Jkg1 K1 790 995 769 Table 2 Anisotropic thermo-hydro-mechanical parameters for the Callovo–Oxfordian claystone, Opalinus and Boom Clay [8] Parameters Units COx OPA Boom Clay Kkto bedding m23:91020 51020 41019 K?to bedding m21:31020 11020 21019 Ekto bedding MPa 8000 8000 400 E?to bedding MPa 5000 4000 200 mkto bedding ðmkkÞ – 0.21 0.35 0.125 m?to bedding ðm?kÞ – 0.35 0.25 0.25 G?to bedding MPa 2500 2300 80 kkto bedding Wm1K11.88 2.4 1.65 k?to bedding Wm1K11.25 1.3 1.31 Acta Geotechnica 123 observation points. The corresponding coordinates of these points are summarized in Table 5. The initial and boundary conditions are divided into three phases to reproduce the evolution of the THM behaviour in the near field of the gallery from the excavation of the gallery up to the heating of the rock (compare Fig. 2). Excavation phase: The total stress at the gallery wall AE is linearly decreased in 24 h from the initial stress to 5% for COx and OPA and to 50% for Boom Clay, while the total in situ horizontal stress remains constant along CB and the total in situ vertical stress is prescribed on CD. Dirichlet boundary conditions are assigned to the symmetry boundaries AB and DE. The pore water pressure at the gallery wall is linearly reduced in 24 h to 0.1 MPa, while on boundaries BC and CD a constant pore pressure of 4.7 MPa for COx and OPA and 2.25 MPa for Boom Clay is prescribed. Impervious boundaries are assigned on the symmetry boundaries AB and DE. Waiting phase: The boundary conditions remain unchanged with respect to the end of the excavation phase and are held for six months, allowing water to drain towards the tunnel. No thermal flux is applied at this stage. Heating phase: The mechanical boundary conditions are similar to the previous phases. Impervious hydraulic boundary conditions are now assigned to the gallery wall. In this phase, a heat flux of 200 Wm1is applied to the gallery wall over a period of ten years. The initial and boundary conditions are summarized in Tables 3and 4. The boundary conditions around the gallery are set at the interface between the gallery and the host rock. The water properties to be considered in this benchmarking exercise are a density of 1000 kg m3, a specific heat capacity of 4180 Jkg1K1and a compressibility at 40 Cof4:5104MPa1. The evolution of the volumetric thermal expansion coefficient of water is given as a polynomial fit of temperature under atmospheric pressure following the Eq. [11]: aw104½C1¼4106½C4T30:001½C3T2þ 0:1404½C2T0:3795½C1 with Tbeing the temperature in C. The evolution of the water viscosity as a function of the temperature under atmospheric conditions is approximated using Vogel’s formula: l¼exp AþB CþT  with Tbeing the temperature in K and lin mPas1. The fitting parameters have the following values: A¼3:719½,B¼578:919 K and C¼137:546 K [9] Seven modelling teams from across Europe were requested to build a near-field 2D generic model to simulate the three subcases presented above. The teams were organized in such a way that each host rock was studied by several teams and thus allowing a comparison of results and increasing confidence in the modelling work. The teams were ANDRA, BGE, EDF, EIG EURIDICE, LEI, University of Lie `ge and UPC Barcelona. Table 6gives an overview of the different teams and their corresponding codes. In total six numerical codes have been used by the different teams for the modelling of the different benchmarks. 4 Implementation of thermo-hydromechanical (THM) models in various numerical codes In the implementation of thermo-hydro-mechanical (THM) governing equations for porous media, the various numerical codes involved in this benchmark adopt different methodologies based on their core numerical strategies and intended applications. Below, we discuss the key features of the THM equations implemented in the codes FLAC3D, OpenGeoSys, COMSOL, LAGAMINE, Code_Aster, and Code_Bright. We use capital letters for the divergence (Div) and gradient (Grad) operators to clearly distinguish Fig. 1 Observation points in the analysis domain, see Table 5for coordinates [8] Acta Geotechnica 123 between the Lagrangian formulation (following a moving material point) and the Eulerian formulation (focused on a fixed spatial point). This convention helps to avoid confusion and emphasizes the specific frame of reference being used in the different codes. 4.1 Energy balance In CODE_ASTER the energy balance derived by [7]as presented in the previous section is expressed in the Lagrangian framework by: X p;c hmp cmp c |fflfflfflfflfflffl{zfflfflfflfflfflffl} Energy transfer due to mass þ _ Q0 |{z} External heat source þX p;c Div hmp cMp c  |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} Divergence of enthalpy flux þDiv q |fflffl{zfflffl} Divergence of heat flux X p;c Mp cFm |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} Work done by mass forces ¼H |{z} Energy dissipation term where hmp crepresents the specific enthalpy of component c in phase p;mp cdenotes the mass of component cin phase p; _ Q0is the rate of heat addition from an external source; Div hmp cMp c  indicates the divergence of the enthalpy flux of component cin phase p; Div qis the divergence of the Fig. 2 Boundary conditions at different stages of the analyses [8] Table 3 Summary of the boundary conditions applied on the gallery wall (compare Fig. 2)[8] Phases Mechanical conditions Hydraulic conditions Thermal conditions T0T0þ24 hðÞ: excavation rrelease up to 5 % of the initial in situ stress for COx and OPA, and up to 50 % for Boom Clay P0!Patm (0.1 MPa) in 24 h No flow at the borehole wall T0þ24 hðÞT0þ6 monthsðÞ: waiting rat 5 % of the initial in situ stress P0¼Patm As above T0þ6 monthsðÞT0þ10 yearsðÞ: heating As above No flow Const. power (200 Wm1) Table 4 Assumed initial conditions in the three host rocks [8] Parameters Components COX OPA Boom Clay Isotropic case Total stress MPa r012.5 12.5 4.5 Pore pressure MPa P04.7 4.7 2.25 Temperature CT022 22 16.5 Anisotropic case Total stress MPa rxx 12.4 2.2 3.825 ryy 12.7 4.0 4.5 rzz 16.4 6.5 3.825 Pore pressure MPa P04.7 2.1 2.25 Temperature CT022 22 16.5 Acta Geotechnica 123 heat flux vector q;Mp crepresents the hydraulic flux of component cin phase p, corresponding to win the Eulerian configuration; Fmdenotes the force per unit mass acting on the components due to external fields; and Hrepresents the energy dissipation term, accounting for heat generation and internal dissipation. In LAGAMINE, the energy balance equation is formulated in a weak form for an arbitrary virtual temperature field Twithin the current deformed configuration Xt, whose boundary is denoted by Ct r.: Z Xt _ St T |{z} Enthalpy evolution Tft T;i |{z} Heat flow oT oxt idXt ¼Z Xt Qt T |{z} Heat sink term TdXtZ Ct qT  qt T |{z} heat flux TdCt In Code_Bright, the balance of internal energy for the medium is defined by: o otESqSð1/Þ |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} Energy of the solid phase þElqlSl/ |fflfflfflffl{zfflfflfflffl} Energy of the liquid phase þEgqgSg/ |fflfflfflfflffl{zfflfflfflfflffl} Energy of the gas phase 2 6 6 6 6 6 4 3 7 7 7 7 7 5 þdiv ic iþjEs jþjEl jþjEg j  |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Energy fluxes : conduction and advection ¼fQ |{z} Internal or external energy supply where ic idenotes the conductive energy flux within the porous medium; fQrepresents sources of energy, both internal and external; jE icorresponds to the advective energy flux resulting from mass movement; and Esignifies the specific internal energy. In COMSOL, we have: o otqcpT |ffl{zffl} Energy storage term ðheat capacityÞ 2 6 6 6 6 4 3 7 7 7 7 5þdiv kgrad TðÞ |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} Divergence of heat flux ðconductionÞ þdiv qcpuT  |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} Divergence of convective heat transfer ¼Q |{z} Heat source or sink in FLAC3D the thermal energy balance is expressed as: Table 5 Coordinates of observation points in meter Points Coordinates Points Coordinates P1 (1.25, 0.0) P8 (0.0, 1.9) P2 (1.9, 0.0) P9 (0.0, 2.5) P3 (2.5, 0.0) P10 (0.0, 6.25) P4 (6.25, 0.0) P11 (0.0, 50.0) P5 (50.0, 0.0) P12 (0.0, 100.0) P6 (100.0, 0.0) P13 (0.88, 0.88) P7 (0.0, 1.25) P14 (1.33, 1.33) P15 (1.77, 1.77) Table 6 Involved modelling teams and their respective codes Organization Country Numerical code Boom Clay COx OPA ANDRA France COMSOL, Code_Aster U BGE Germany OpenGeoSys, FLAC3D UUU EDF France Code_Aster UU EIG EURIDICE Belgium COMSOL U LEI Lithuania COMSOL U ULiege Belgium Lagamine UU UPC Spain CODE_BRIGHT UU Table 7 Correspondence of energy balance terms across software Concept LAGAMINE Code_Aster Code_Bright COMSOL FLAC3D OpenGeoSys Energy storage _ St TTDiv hmp cMp c  HESqSð1/ÞþElqlSl/þEgqgSg/qcpTnTqcpT Heat flux (conductive) ft T;iqic ikgrad TqT ikgrad T Divergence of heat flux  qt TDiv qdivðic iþjEs jþjEl jþjEg jÞdivðkgradTÞDivðqT iÞdivðkgradTÞ Heat flux (advective) – Pp;chmp cmp cjEs jþjEl jþjEg jDivðqcpuTÞ– divðqcpuTÞ Heat source or sink Qt T _ Q0fQQqT vQ Acta Geotechnica 123 Div qT i |fflfflfflfflffl{zfflfflfflfflffl} Heat flux divergence þqT v |{z} Volumetric heat source intensity ¼onT ot |{z} Rate of heat storage per unit volume Two process models in OpenGeoSys were employed in this benchmarking: a non-isothermal Richards flow model coupled with mechanics, and a thermo-hydraulic (TH) model using thermo-mechanical storage coefficients [5]. The first model follows a monolithic approach, deriving hydraulic-mechanical couplings from three-dimensional effective stress and omitting vapour diffusion based on benchmark assumptions. The second TH model incorporates mechanical effects within the mass balance of the TH process via thermo-mechanical storage coefficients, assuming uniaxial strain (xx ¼yy ¼0) and constant vertical stress (rzz), while also considering thermal stresses from constrained transverse thermal expansion. The energy balance for the non-isothermal Richards flow model is given below: o otqcpT  |fflfflfflfflffl{zfflfflfflfflffl} Energy storage ðheat capacityÞ þdiv kgradTðÞ |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} Divergence of heat fluxðconductionÞ þdiv qcpuT  |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} Convective heat transfer ¼Q |{z} Heat source or sink Table 7provides a concise comparison of energy balance components across six THM simulation software: LAGAMINE, Code_Aster, Code_Bright, COMSOL, FLAC3D, and OpenGeoSys, each employing unique formulations. Energy storage is represented differently, with LAGAMINE and Code_Bright considering phase-specific energies, while COMSOL, OpenGeoSys, and FLAC3D use simpler forms such as qcpT. Code_Aster employs a divergence term related to enthalpy flux. Heat flow (conduction) is consistently expressed via Fourier’s law across software, commonly as kgrad T, with specific forms like ft T;iin LAGAMINE and ic iin Code_Bright. Heat flux divergence describes the distribution of conductive heat, often using divergence operators; Code_Bright uniquely integrates conduction and advection to account for multiphase effects. Convective heat transfer, noted in Code_Aster, Code_Bright, COMSOL, and OpenGeoSys, includes heat transported by fluid flow, while FLAC3D omits this term, focusing on solid mechanics. Each software includes a heat source or sink term, with labels as Qt T(LAGAMINE), _ Q0(Code_Aster), fQ(Code_Bright), Q(COMSOL, OpenGeoSys), and qT v(FLAC3D), allowing for a range of energy interactions. This overview emphasizes each software’s adaptability, from simplified single-phase models to multi-phase approaches for comprehensive THM simulations. 4.2 Fluid mass balance In LAGAMINE, The fluid mass balance equation is expressed in a weak form for any admissible virtual pore pressure field p w: Z Xt _ Mt wp w |fflffl{zfflffl} Time derivative of fluid mass ft w;i op w oxt i |fflfflffl{zfflfflffl} Mass flow gradient 0 B B B B B B B @ 1 C C C C C C C A dXt ¼Z Xt Qt wp w |fflffl{zfflffl} Fluid sink term dXtZ Ct qw  qt wp w |ffl{zffl} Boundary mass flux dCt Here, ft w;irepresents the mass flow rate, Qt wdenotes a term accounting for fluid sinks, and Ct qwis the boundary segment where the incoming fluid mass per unit area, given by  qt w,is specified. Additionally, _ Mt wdenotes the time rate of change of fluid mass. The fluid mass balance in FLAC3D is defined by: Div qi |fflfflfflffl{zfflfflfflffl} Divergence of fluid flux þqv |{z} Volumetric fluid source intensity ¼of ot |{z} Time derivative of fluid content In COMSOL, the water mass (mw) balance equation includes additional water source terms to take into account the different coupling processes (HM, TH and TM) and is expressed as follows: omw ot |{z} Time derivative of water mass þdiv qwqw ðÞ |fflfflfflfflfflffl{zfflfflfflfflfflffl} Divergence of mass flux due to water flow ¼QH |{z} Water source In Code_Bright, The total mass balance of water in the liquid phase is given by: Acta Geotechnica 123 o otxw lqlSl/þxw gqgSg/ hi |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Rate of change of water mass in liquid and gas phases þjl;w i;iþjg;w i;i  |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} Divergence of water mass flux in liquid and gas phases ¼fw |{z} External water supply Here, xw land xw grepresent the mass fraction of the component (water) relative to the total mass of the liquid and gas phases, respectively. Sland Sgare the saturation degrees of the liquid and gas phases, indicating the fraction of pore volume each phase occupies. qland qgdenote the densities of the liquid and gas phases, respectively. jl;w i represents the mass flux of water within the liquid phase, while jg;w irepresents the mass flux of water within the gas phase. Lastly, fwindicates an external water source. The liquid water mass balance equation for hydraulics in Code_Aster is defined by: omw ot |{z} Rate of change of liquid water mass content þoMw;i oxi |fflffl{zfflffl} Divergence of liquid water mass flow ¼qv |{z} Volumetric source term In OpenGeoSys, the fluid mass balance equation is given by: omw ot |{z} Rate of change of water mass content þdiv jw |fflffl{zfflffl} Divergence of water mass flux ¼qv |{z} Volumetric source or sink term where:mw: represents the mass of water per unit volume in the porous medium. jw: water mass flux vector, defined by Darcy’s law. qv: volumetric source or sink term, representing water sources or sinks per unit volume. Table 8summarizes the correspondence of fluid mass balance terms across different THM codes. In LAGAMINE, the weak form of the fluid mass balance includes explicit terms for the time derivative of fluid mass, mass flow gradient, fluid sink term, and boundary mass flux. In Code_Aster and OpenGeoSys, the formulations emphasize the rate of change of water mass content and divergence of water mass flow, aligning with terms in FLAC3D and COMSOL where fluid storage and transport mechanisms are central. Code_Bright introduces specific terms for the liquid and gas phase fluxes and the external water supply, necessary for multiphase behaviour. These distinctions reflect differences in handling boundary conditions and explicit or implicit source and sink terms across the codes. 4.3 Mechanical balance equation for momentum conservation The mechanical balance equation, characterizing the conservation of momentum in porous media, considers both the effective stresses within the solid matrix and the body forces acting upon it. This formulation is generally valid across multiple THM modelling codes and is expressed as follows: rij;jþbi¼0 where rij is the component of stress, birepresents body force. Table 8 Correspondence of fluid mass balance terms across THM codes Concept LAGAMINE Code_Aster Code_Bright COMSOL FLAC3D OpenGeoSys Time derivative of mass _ Mt womw=oto otxw lqlSl/þxw gqgSg/ hi omw=otof=otomw=ot Mass flow / flux ft w;i=oxt ioMw;i=oxijl;w i;iþjg;w i;idiv ðqwuÞDiv qidiv jw Volumetric source / sink term Qt wqvfwQHqvqv Boundary mass flux  qt w–– ––– Acta Geotechnica 123 4.4 Challenges in comparing implementations Comparing the implementations of THM processes in various numerical codes is challenging due to differences in frameworks (Eulerian vs. Lagrangian), numerical schemes, and handling of multi-physics couplings. Each code has unique methodologies in formulating balance equations, applying weak or strong forms, and managing coupling effects, especially in multiphase flow, thermal conduction, and mechanical deformation. This variability stems from the distinctive use cases each software targets, whether in geotechnical applications, environmental engineering, or multi-phase fluid transport in porous media. Therefore, benchmarking initiatives are essential for verifying and comparing the performance and accuracy of the THM models across different codes. 5 Constitutive equations for the skeleton In elasticity theory, the Hooke’s stress-strain relation is formulated as follows: eij ¼Cijkl rkl where eand rare respectively the strain and stress tensor. Cijkl are the coordinates of the fourth-order compliance tensor. (C) is the corresponding compliance matrix in Kelvin notation which we will specify below for the materials used in this benchmark.. In the following, we are presenting the compliance matrices of the materials used in this benchmark. The Isotropic Elastic Compliance Matrix For linear isotropic elastic materials, the compliance matrix CC depends on two independent material properties, such as the Young’s modulus Eand the Poisson’s ratio m. In the most general isotropic form, the compliance matrix is expressed as: C¼ 1 Em Em E000 m E 1 Em E000 m Em E 1 E000 000 1 2G00 0000 1 2G0 00000 1 2G 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 In this expression, the shear modulus G is given by: G¼E 2ð1þmÞ Transverse Isotropic Elastic Compliance Matrix For cross-anisotropic or transverse isotropic elasticity, five independent material parameters are required to describe the elastic behavior: Ejj;E?;mjjjj;mjj?;and Gjj? Here, the subscripts jjand ? refer to directions parallel and perpendicular to the isotropic planes, respectively. In this case, the elastic compliance matrix is defined as [17]: C¼ 1 Ekm?k E?mkk Ek 000 mk? Ek 1 E?mk? Ek 000 mkk Ekm?k E? 1 Ek 000 0001 2Gk? 00 0000 1 2Gkk 0 00000 1 2G?k 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 The symmetry of the compliance tensor imposes the following relationship between the Poisson’s ratios and Young’s moduli: mk? Ek¼m?k E?The shear modulus in the isotropic planes ðGjjjjÞcan be expressed as: The symmetry of the stress and strain tensors implies that the other shear moduli are equal: Gk? ¼Gkk 6 Analytical benchmark The first step towards model verification is to compare the accuracy of the THM models against well-defined analytical solutions. This section presents a benchmark analysis for the analytical solution of the coupled thermo-hydromechanical (THM) consolidation problem, originally derived by [4] and later corrected for the effective stress term by [6]. The aim of this benchmark is to validate the theoretical formulation of the involved THM models through coordinated efforts by several teams. The benchmark analysis will focus on three primary variables: temperature, pore pressure, and displacement of the solid skeleton. Problem overview The problem describes a heat source embedded in a fully fluid-saturated porous medium. Due to spherical Acta Geotechnica 123 expansivity of OPA clay, along with a lower Young’s modulus, which allows for pore pressure dissipation and reduces the thermal expansivity contrast between water and the clay matrix. This overpressure causes stress relief, leading to tensile stress at the gallery wall with a maximum value of 4.5 MPa at the end of the simulation. At P2, located 0.65 m behind the wall in the horizontal direction, a compressive stress of more than 4 MPa can be observed at the end of the simulation. 7.4 Boom Clay: elastic anisotropic case The results of the benchmark simulation for Boom Clay, taking into account anisotropic material properties, are presented in Fig. 9. EUR, ULG and BGE with two process modell participated in this benchmark. In exception of BGE results computed with the hydrothermal model, the results coming from THM process models are consistent for all variables. Pore pressure evolution coming from the hydrothermal model in OpenGeoSys with incorporation of thermo-mechanical effects overestimate the pore overpressure in the heating phase when compared to the results obtained with THM formulations. Further investigations are necessary to adjust this model. Overall, the evolution of pore pressure is similar to the isotropic case, with all pore pressure evolution curves converging towards a value of 2.8 MPa at the end of the simulation. Regarding the mechanical variables, good consistency has been obtained for the displacement and stress results. Multiple iterations were required to achieve a good consistency in this benchmark case. Discrepancies arose in the early stages from ambiguities in the benchmark specifications regarding the anisotropic Poisson ratios. It’s important to note that nuij 6¼ nuji for transversely isotropic materials like argillaceous materials. To properly calculate anisotropic components of the Poisson ratios, the orientation of both the numerical model and isotropic plane must be taken into account. Effective communication between teams was key to resolving this issue. As a result of thermal pressurization, tensile stresses occur at the gallery wall in the horizontal direction. The effective tensile stresses reach a value of nearly 1 MPa at the end of the simulation, which may cause tensile damage Fig. 9 Benchmark results for the elastic anisotropic case in Boom Clay [8] Acta Geotechnica 123 in the rock if the tensile strength of the rock is exceeded. However, this issue is limited to the boundary because compressive stresses are observed at P2. 7.5 Callovo-Oxfordian claystone: elastic anisotropic case The results of the subcase 2 for COx clay, assuming anisotropic material properties, are shown in Fig. 10. Seven teams participated in this benchmark (AND, BGE, EDF, EUR, LEI, ULG and UPC), and overall, good agreement has been achieved by all teams. The temperature, displacement and deviatoric stress evolution curves are coincident at all points. During the heating phase, five teams predicted the same pore pressure evolution trend. Results of BGE and EDF deviate from this trend. The reason for this discrepancy has not yet been found, but differences in the anisotropic Poisson ratios used by these two teams are suspected. The evolution of the effective vertical and horizontal stresses shows small discrepancies at each observation points that can be related to the difference in mesh size and element formulation. This will be further elaborated in the discussion section. The anisotropy of the COx clay appears to amplify the thermally induced pore pressure, as a maximum value of about 11 MPa was observed in the anisotropic case, compared to about 10 MPa in the isotropic case. The effective horizontal stresses reach approximately 7.5 MPa at the gallery wall, decreasing to less than 1 MPa at P2, which is located 0.65 m behind the gallery wall. At P4, which is located 0.5 m behind the gallery wall, the effect of thermal pressurization is lowest. Fig. 10 Benchmark results for the elastic anisotropic case in COx claystone [8] Acta Geotechnica 123 7.6 Opalinus claystone: elastic anisotropic case The results of the final benchmark for Opalinus Clay, taking into account anisotropic material properties, are presented in Fig. 11. The same teams that worked on the isotropic case for Opalinus Clay also participated in this benchmark. The overall results from all teams are in good agreement, with consistent predictions for temperature evolution. UPC predicted the highest pore pressure evolution in this benchmark. EDF and BGE obtained almost identical pore pressure evolution results. In contrast to COx, the results for the anisotropic case in Opalinus Clay show that less overpressure is generated compared to the isotropic case. A maximum overpressure value of about 6 MPa was obtained in this benchmark, while the isotropic case had a maximum value of 8.5 MPa. This may be explained by the higher permeability and higher Poisson ratio in the horizontal direction of the anisotropic Opalinus Clay. For Opalinus Clay under anisotropic conditions, the displacement results show that the surrounding region moves away from the gallery wall due to thermal expansion in the heating phase. This can be observed at P2, P3 and P4, with P4, located 6.5 m deep in the rock, moving 3 mm away from the gallery wall at the end of the simulation. The effective stress evolution resulting from the pore overpressure shows tensile stresses at the boundary of the gallery wall at P1, which vanish at P2. This effect has been observed in all benchmarks and has been discussed already in previous sections. The comparative assessment of all stress results is in good agreement for all teams except for observation point P1 where UPC predict less effective vertical and deviatoric stresses during the waiting phase. 8 Discussion Overall, the results of the benchmark studies conducted by the various teams involved showed a high level of consistency in the thermo-hydro-mechanical response of the clay-based materials being studied. These results were more consistent in the isotropic cases, wheareas some discrepancies were observed in the anisotropic benchmarks. In the anisotropic cases, the assumptions taken for Fig. 11 Benchmark results for the elastic anisotropic case in OPA claystone [8] Acta Geotechnica 123 computing the anisotropic components of the Poisson ratios in respect to the coordinate systems used in the different codes have a huge effect on the outcome of the benchmark. There are responsible for the discrepancies in the pore pressure evolution. This can be avoided if the Poisson’s ratios are consistently defined in all codes. Other differences in the results obtained by the different teams may be due to a variety of factors, including differences in the governing equations and different solution schemes and codes used, and differences in the numerical meshes employed. A more mathematical investigation is needed to quantify the effect of these factors. This is tedious to perform as much assumptions in this regard are not published and some numerical codes used in this study are not open source. To better understand the effect of the mesh size and element formulation on the obtained results, a mesh sensitivity study has been performed. Eight different meshes were developed in this study that can be divided into two groups. The first group consists of linear meshes whereas the second group is made of quadratic counterparts derived from the linear meshes. In each group, three meshes were made of quadrilateral elements. The second mesh is derived from the first one with 3040 elements by dividing each element by four to obtain 12160 elements. The third mesh is generated in the same manner by dividing each element of the second mesh by four to obtain 48460 elements. The fourth mesh has been generated using triangular elements. All these meshes were employed to carry out the isotropic benchmark for Boom Clay using OpenGeoSys by BGE. The results of this study are shown in Fig. 12. For all variables, a mesh dependency on the results is observed especially at the observation point P2 where the gradients are highest. The discrepancy in the pore pressure results is observed during the excavation (until 1d) and waiting phase (until 6 months) whereas in the heating phase results Fig. 12 Comparison of the THM response in OPA at the same distances from the gallery wall parallel (P2, P3) and perpendicular to bedding (P7, P8) based on UPC results Acta Geotechnica 123 from all meshes are identical. In the excavation phase, some oscillations are observed for all quadrilateral meshes with linear shape functions. This is also the case of the quadratic mesh with 3040 quadrilateral cells. It is important to notice that results with triangular meshes do not show this artefact. Thus, triangular element formulation seems to be more robust for THM simulation at least in OpenGeoSys. This needs to be tested for other codes. The stress evolution curves shows also some scattering of results at observation point P1 for all linear meshes. This is probably the result of the extrapolation problem due to the high gradients near the gallery wall that has been noticed by ULG for the anisotropic cases for Boom Clay and COx in combination with the low approximation order of the derived quantity stress. ULG results for these anisotropic cases at P1 presented in this work have been extrapolated from the integration point to the gallery wall to avoid this issue. The fact that the meshes of the second group do not show this problem further supports this hypothesis. The additional integration points coming with the quadratic elements and the higher approximation order of the stresses help to increase the quality of the extrapolation of the results to the boundary wall. Based on this study, one can conclude that the discrepancy observed in the benchmark results may be partly explained by the different meshes and element formulations employed by the teams. This is especially the case for results scattering observed near the gallery at P2 and to a smaller extent at P3 as one can be observed for instance in the results of the anisotropic benchmark cases, see for example Fig. 10. Figure 12 shows also that the results of the pore pressure evolution are consistent for all the meshes in the heating phase. Thus the discrepancy observed in the pore pressure results cannot be explained by its dependency on the spatial discretization. The replication of this study based on an anisotropic benchmark case may be necessary to confirm these conclusions. Further analyses may be needed to understand the other underlying mechanisms leading to the discrepancies. However, it can be concluded at this stage that all of the teams involved in this study were able to accurately model the THM evolution of heat-generating repository systems in clay formations and that the tools and techniques they used can be considered verified in this context based on their ability to produce similar results. The anisotropic effects observed in the benchmark studies were more pronounced in the temperature evolution and mechanical response of the clay-based materials but were not evident in the pore pressure evolution. Instead, the pore pressure appeared to homogenize during the heating phase, reaching similar asymptotic values at the various observation points. This suggests that the anisotropic effects may have a greater influence on the temperature and mechanical behaviour of the materials, but not on the pore pressure, see compare P2 vs. P8 and P3 vs. P9 in Fig. 13. For the isotropic case, the results at the points situated at the same distance from the gallery wall are identical as shown in Fig. 13. According to the benchmark studies, an increase in pressure was observed during the heating phase in all cases. This phenomenon, known as thermal pressurization, occurs due to the difference in the thermal expansion of water and the clay matrix. The water is constrained by the less expansive clay matrix, leading to an accumulation of pressure. The resulting tensile stresses that may occur in the clay formation as a result of thermal pressurization were observed to be local and to vanish quickly with depth in the rock in the studied benchmarks. However, it is Fig. 13 Comparison of the THM response in OPA at the same distances from the gallery wall parallel (P2, P3) and perpendicular to bedding (P7, P8) based on UPC results Acta Geotechnica 123 important to note that these results will not be the same in a typical repository configuration with backfilled galleries and multiple canisters disposed of adjacent to one another. Despite this, the results suggest that these tensile stresses may remain local even in this benchmark case and may disappear over time as the thermal power of the radioactive waste decays and the pore pressure decreases. Depending on the tensile strength of the clay, tensile stresses may cause local damage to the clay. Thermo-hydro-mechanical (THM) simulations can be time-consuming and computationally demanding, particularly when used to assess the safety of repository systems that require the consideration of large geological formations in the numerical analysis. In these cases, the use of numerically efficient methods may be essential to handle the computational demands of such simulations. The TH model with thermo-mechanical storage coefficients proposed by [5] and implemented in the OpenGeoSys has been shown to be able to adequately reproduce the pore overpressure for the isotropic benchmark in Boom Clay. The model slightly overpredicted the thermal induced pore pressure response in the anisotropic case. Further work is necessary to improve the model predictions in such situation. Nevertheless, this model has been found to result in a significant speed-up, up to two orders of magnitude, compared to traditional THM simulations. The actual speed-up depends on the size of the model and the most significant benefits are typically seen in models with a large number of degrees of freedom, where complexity reduction is particularly important. However, for a detailed evaluation of the stress and displacement evolution, fully coupled THM simulations remain necessary [5]. 9 Conclusions The EURAD work package HITEC seeks to enhance the understanding of thermo-hydro-mechanical processes in clay-based materials subjected to elevated temperatures. To this end, a benchmark initiative was conducted to assess the current state of modelling THM phenomena in clay materials. The results of this initiative, which involved teams from across Europe, showed that these teams and their numerical tools are capable of accurately predicting the THM behaviour of clay-based materials such as Boom Clay, Opalinus Clay, and Callovo-Oxfordian clay. In general, the results of the isotropic benchmarks were highly consistent, while some discrepancies were observed in the anisotropic benchmarks, potentially due to differences in the governing equations and codes used, as well as mesh size, element formulations and interpolation algorithms. In conclusion, the results of this benchmark demonstrate the expertise and capabilities of the participating teams in modelling THM phenomena for the safety assessment of repository systems in clay formations. The involved codes can be seen as verified within the framework of ‘‘Validation & Verification’’ (V&V). Further studies that take into account effects of plasticity and permeability changes due to damage are necessary to better reproduce benchmarks at repository conditions. This will be the subject of further publication of the involved teams. Acknowledgements The project leading to this work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 847593. This paper describes objective technical results and analysis. The authors alone are responsible for the contents of this study. Author contribution First author wrote the main manuscript. Other authors contributed by providing part related to their numerical analyses. Funding Open Access funding enabled and organized by Projekt DEAL. Data availability No datasets were generated or analysed during the current study. Declarations Conflict of interest The authors declare no conflict of interest. References 1. Birkholzer Jens T, Bond Alexander E (2022) DECOVALEX2019: an international collaboration for advancing the understanding and modeling of coupled thermo-hydro-mechanicalchemical (THMC) processes in geological systems. Int J Rock Mech Min Sci 154(154):105097. https://doi.org/10.1016/j. ijrmms.2022.105097 2. Birkholzer Jens T, Bond Alexander E, Hudson John A, Lanru Jing, Chin-Fu Tsang, Hua Shao, Olaf Kolditz (2018) Decovalex2015: an international collaboration for advancing the understanding and modeling of coupled thermo-hydro-mechanicalchemical (THMC) processes in geological systems. Environ Earth Sci 77(14):56. https://doi.org/10.1007/s12665-018-7697-7 3. Birkholzer Jens T, Chin-Fu Tsang, Bond Alexander E, Hudson John A, Lanru Jing, Ove Stephansson (2018) 25 years of DECOVALEX-Scientific advances and lessons learned from an international research collaboration in coupled subsurface processes. Int J Rock Mech Min Sci 122:103995. https://doi.org/10. 1016/j.ijrmms.2019.03.015 4. Booker JR, Savvidou C (1985) Consolidation around a point heat source. Int J Numer Anal Methods Geomech 9(2):173–184 5. Buchwald J, Kaiser S, Kolditz O, Nagel T (2021) Improved predictions of thermal fluid pressurization in hydro-thermal models based on consistent incorporation of thermo-mechanical effects in anisotropic porous media. Int J Heat Mass Transf 172:121127. https://doi.org/10.1016/j.ijheatmasstransfer.2021. 121127 6. Chaudhry AA, Buchwald J, Kolditz O, Nagel T (2019) Consolidation around a point heat source (correction and verification). Int J Numer Anal Methods Geomech 43(18):2743–2751. https:// doi.org/10.1002/nag.2998 Acta Geotechnica 123 7. Coussy O (2004) Poromechanics. John Wiley and Sons, Chichester, UK 8. Christophe (Andra) de Lesquen, Minh-ngoc (Andra) Vu, Eric (BGE) Simo, Alexandru (BGE) Tatomir, Paola Le ´on (BGE) Vargas, Pierre (CNRS-UGrenoble) Be ´suelle, Stefano (CNRSUGrenoble) Dal Pont, Alice (CNRS-UGrenoble) di Donna, Nicola ´s (CNRS-UGrenoble) Zalamea, Simon (EDF) Raude, Ginger (EDF) El tabbal, Arnaud (EURIDICE) Dizier, Suresh (SCK CEN) Seetharam, Asta (LEI) Narkuniene, Fre ´de ´ric (ULie `ge) Collin, Hangbiao (ULie `ge) Song, Abhishek (ULie `ge) Rawat, Antonio (UPC) Gens, Fei (UPC) Song HITEC-Deliverable-7.6, Andra, BGE, CNRS-UGrenoble, EDF, EURIDICE, SCK CEN, LEI, ULie `ge, UPC 9. DDBST (2023) Liquid dynamic viscosity calculation by vogel equation. Dortmund Data Bank Software & Separation Technology (DDBST). URL http://ddbonline.ddbst.de/VogelCalcula tion/VogelCalculationCGI.exe?component=Water. Accessed 25 Dec 2020 10. Dwivedi D, Arora B, Molins S, Steefel CI (2016) Benchmarking reactive transport codes for subsurface environmental problems. Groundw Assess Model Manag. https://doi.org/10.1201/ 9781315369044-19 11. Kell George S (1975) Density, thermal expansivity, and compressibility of liquid water from 0. deg. to 150. deg. correlations and tables for atmospheric pressure and saturation reviewed and expressed on 1968 temperature scale. J Chem Eng Data 20(1):97–105. https://doi.org/10.1021/je60064a005 12. Lehmann C, Bilke L, Buchwald J, Graebling N, Grunwald N, Heinze J, Meisel T, Renchao L, Naumov D, Rink K, Sen O, Selzer P, Shao H, Wang W, Zill F, Nagel T, Kolditz O (2024) Openworkflow-development of an open-source synthesis-platform for safety investigations in the site selection process; [openworkflow - entwicklung einer open-source-synthese-plattform fu ¨r sicherheitsuntersuchungen im standortauswahlverfahren]. Grundwasser 29(1):31–47. https://doi.org/10.1007/ s00767-024-00566-9 13. Lux KH, Rutenberg M, Feierabend J, Czaikowski O, Friedenberg L, Maßmann J, Pitz M, Sentis ML, Graupner BJ, Hansmann J, Hotzel S, Kock I, Rutqvist J, Hu M, Rinaldi AP (2023) BenVaSim–International Benchmarking for Verification and Validation of TH2M Simulators with Special Report December 2021, Chair for Waste Disposal Technologies and Geomechanics Clausthal University of Technology 14. Oberkampf William L, Trucano Timothy G, Charles H (2004) Verification, validation, and predictive capability in computational engineering and physics. Appl Mech Rev 57(5):345–384. https://doi.org/10.1115/1.1767847 15. OGS-Community (2024) Saturated point heat source benchmark. 1115https://www.opengeosys.org/docs/benchmarks/th2m/satur atedpointheatsource/. Accessed 11 Sep 2024 16. Pitz M, Kaiser S, Grunwald N, Kumar V, Buchwald J, Wang W, Naumov D, Chaudhry AA, Maßmann J, Thiedau J, Kolditz O, Nagel T (2023) Non-isothermal consolidation: a systematic evaluation of two implementations based on multiphase and Richards equations. Int J Rock Mech Min Sci. https://doi.org/10. 1016/j.ijrmms.2023.105534 17. Pardoen B (2015) Hydro-mechanical analysis of the fracturing induced by the excavation of nuclear waste repository galleriesusing shear banding. ULie `ge - Universite ´de Lie `ge http:// reflexions.ulg.ac.be/en/NuclearWasteStorage Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Authors and Affiliations Eric Simo 8,1 •Christophe de Lesquen 2 •Rocio Paola Leon-Vargas 1 •Minh-ngoc Vu 2 •Simon Raude 3 • Ginger El Tabbal 3 •Arnaud Dizier 4 •Suresh Seetharam 4 •Asta Narkuniene 5 •Fre ´de ´ric Collin 6 • Hangbiao Song 6 •Antonio Gens 7 •Fei Song 7 •Alexandru-Bogdan Tatomir 10 •Thomas Nagel 8,9 • Jo ¨rg Buchwald 8,9 &Eric Simo [email protected] 1 BGE TECHNOLOGY GmbH, Eschenstrasse 55, 31224 Peine, Germany 2 ANDRA, 1/7 rue Jean Monnet, 92290 Cha ˆtenay-Malabry, France 3 EDF R&D, 7 Bvd Gaspard Monge, 91120 Palaiseau, France 4 Belgian Nuclear Research Centre (SCK CEN), Boerentang 200, B-2400 Mol, Belgium 5 Lithuanian Energy Institute, Breslaujos str., 44403 Kaunas, Lithuania 6 Urban and Environmental Engineering Research Unit, Universite ´de Lie `ge, Alle ´edelaDe ´couverte 9, 4000 Lie `ge, Belgium 7 Department of Civil and Environmental Engineering, Universitat Politecnica de Catalunya, Jordi Girona 1-3, 08034 Barcelona, Spain 8 Geotechnical Institute, TU Bergakademie Freiberg, GustavZeuner-Straße 1, 09599 Freiberg, Germany 9 Helmholtz Centre for Environmental Research-UFZ, Permoserstraße 15, 04318 Leipzig, Germany 10 BGE mbH, Eschenstrasse 55, 31224 Peine, Germany Acta Geotechnica 123