# EURAD state-of-the-art report: development and improvement of numerical methods and tools for modeling coupled processes in the field of nuclear waste disposal

Francis Claret, N. Prasianakis, Attila Baksay, D. Lukin, Guillaume Pépin, Etienne Ahusborde, Brahim Amaziane, G. Bátor, D. Becker, A. Bednár, et al.

► To cite this version:

Francis Claret, N. Prasianakis, Attila Baksay, D. Lukin, Guillaume Pépin, et al., EURAD state-of-the-art report: development and improvement of numerical methods and tools for modeling coupled processes in the field of nuclear waste disposal. Frontiers in Nuclear Engineering, 2024, 3, pp.1437714. 10.3389/fnuen.2024.1437714. [hal-04679659v2](https://hal.science/hal-04679659v2)

HAL Id: hal-04679659

https://hal.science/hal-04679659v2

Submitted on 18 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN ACCESS

EDITED BY
Bernd Grambow,
UMR6457 Laboratoire de Physique
Subatomique et des Technologies Associées
(SUBATECH), France

REVIEWED BY
Pierre Toulhoat,
Académie des Technologies, France
Allan Moreira Mulin Leal,
ETH Zürich, Switzerland

*CORRESPONDENCE

F. Claret,

<figure><img src="image_3.png" /></figure>

f.claret@brgm.fr

RECEIVED 24 May 2024

ACCEPTED 08 August 2024

PUBLISHED 09 October 2024

CITATION

Claret F, Prasianakis NI, Baksay A, Lukin D, Pepin G, Ahusborde E, Amaziane B, Bátor G, Becker D, Bednár A, Béreš M, Bérešová S, Böthi Z, Brendler V, Brenner K, Březina J, Chave F, Churakov SV, Hokr M, Horák D, Jacques D, Jankovský F, Kazymyrenko C, Koudelka T, Kovács T, Krejčí T, Kruis J, Laloy E, Landa J, Ligurský T, Lipping T, López-Vázquez C, Masson R, Meeussen JCL, Mollaali M, Mon A, Montenegro L, Pisani B, Poonoosamy J, Pospiech SI, Sači Z, Samper J, Samper-Pilar A-C, Scaringi G, Sysala S, Yoshioka K, Yang Y, Zuna M and Kolditz O (2024) EURAD state-of-the-art report: development and improvement of numerical methods and tools for modeling coupled processes in the field of nuclear waste disposal.

Front. Nucl. Eng. 3:1437714.

doi: 10.3389/fnuen.2024.1437714

COPYRIGHT

© 2024 Claret, Prasianakis, Baksay, Lukin, Pepin, Ahusborde, Amaziane, Bátor, Becker, Bednár, Béreš, Bérešová, Böthi, Brendler, Brenner, Březina, Chave, Churakov, Hokr, Horák, Jacques, Jankovský, Kazymyrenko, Koudelka, Kovács, Krejčí, Kruis, Laloy, Landa, Ligurský, Lipping, López-Vázquez, Masson, Meeussen, Mollaali, Mon, Montenegro, Pisani, Poonoosamy, Pospiech, Sači, Samper, Samper-Pilar, Scaringi, Sysala, Yoshioka, Yang, Zuna and Kolditz. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

EURAD state-of-the-art report:
development and improvement
of numerical methods and tools
for modeling coupled processes
in the field of nuclear waste
disposal

F. Claret¹*, N. I. Prasianakis², A. Baksay³, D. Lukin⁴, G. Pepin⁵,
E. Ahusborde⁶, B. Amaziane⁶, G. Bátor⁷, D. Becker⁸, A. Bednár⁷,
M. Béreš⁹, S. Bérešová⁹, Z. Böthi³, V. Brendler¹⁰, K. Brenner¹¹,
J. Březina¹², F. Chave¹¹, S. V. Churakov², M. Hokr¹², D. Horák⁹,
D. Jacques¹³, F. Jankovský¹⁴, C. Kazymyrenko¹⁵, T. Koudelka¹⁶,
T. Kovács⁷, T. Krejčí¹⁷, J. Kruis¹⁷, E. Laloy¹³, J. Landa¹², T. Ligurský⁹,
T. Lipping¹⁸, C. López-Vázquez¹⁹, R. Masson¹¹,
J. C. L. Meeussen²⁰, M. Mollaali²¹, A. Mon²², L. Montenegro²²,
B. Pisani²², J. Poonoosamy²³, S. I. Pospiech¹⁰, Z. Saâdi²⁴,
J. Samper²², A.-C. Samper-Pilar²², G. Scaringi¹⁶, S. Sysala⁹,
K. Yoshioka²¹, Y. Yang²³, M. Zuna¹⁴ and O. Kolditz²¹

¹BRGM, Orléans, France, ²Laboratory for Waste Management, Paul Scherrer Institute, Villigen, Switzerland, ³TS Enercon, Budapest, Hungary, ⁴SURAO, Praha, Czechia, ⁵Andra, Chatenay Malabry, France, ⁶Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, Pau, France, ⁷Social Organization for Radioecological Cleanliness, Veszprem, Hungary, ⁸GRS gGmbH, Braunschweig, Germany, ⁹Institute of Geonics, Czech Academy of Sciences, Ostrava, Czechia, ¹⁰Helmholtz Zentrum Dresden-Rossendorf, Dresden, Germany, ¹¹CNRS, Inria, LJAD, University Côte d'Azur, Nice, France, ¹²Technical University of Liberec, Liberec, Czechia, ¹³Belgian Nuclear Research Centre, Mol, Belgium, ¹⁴Waste Disposal Processes and Safety Department, Radioactive Waste and Decommissioning Division, ÚJV Řež, a.s., Husinec, Czechia, ¹⁵EDF R&D, IMSIA - UMR CNRS 9219, Institut des Sciences de la mécanique et Applications industrielles, Palaiseau, France, ¹⁶Institute of Hydrogeology, Engineering Geology and Applied Geophysics, Faculty of Science, Charles University, Prague, Czechia, ¹⁷Faculty of Engineering, Czech Technical University in Prague, Prague, Czechia, ¹⁸Faculty of Information Technology and Communication Sciences, Tampere University, Pori, Finland, ¹⁹LatinGEO Lab IGM+ORT, Universidad ORT Uruguay, Montevideo, Uruguay, ²⁰Nuclear Research and Consultancy Group, Petten, Netherlands, ²¹Helmholtz Centre for Environmental Research UFZ, Leipzig, Germany, ²²Civil Engineering School, Interdisciplinary Center for Biology and Chemistry (CICA), Universidade de A Coruña, Campus de Elviña, A Coruña, Spain, ²³Institute of Energy and Climate Research, Nuclear Waste Management (IEK-6), Forschungszentrum Juelich, Germany, ²⁴Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSE-ENV/SPDR/UEMIS, Fontenay-aux-Roses, France

The Strategic Research Agenda (SRA; https://www.ejp-eurad.eu/publications/eurad-sra) of the European Joint Programme on Radioactive Waste Management (EURAD; https://www.ejp-eurad.eu/) describes the scientific and technical domains and sub-domains and knowledge management needs of common interest between EURAD participant organizations. Theme number 7 is entitled “Performance assessment, safety case development and safety analyses.” A list of research and development priorities and activities of common interest to be addressed within EURAD for theme 7 have been
established. Amongst others, the Understanding and modelling of multi-physical Thermo-Hydro-Mechanical-Chemical coupled processes (THMC) occurring in radioactive waste disposal is a major and permanent issue to support optimization of design and safety case abstraction. To tackle this challenge a research work package entitled "DONUT: Development and improvement of numerical methods and tools for modelling coupled processes" has been conducted within the EURAD join programming initiative. The purpose of this work package is to improve/develop methods or numerical tools in order to go a step further in development of (i) relevant, performant and cutting-edge numerical methods that can easily be implemented in existing or new tools, in order to carry out high-performance computing to facilitate the study of highly coupled processes in large systems, (ii) numerical scale transition schemes for coupled processes, (iii) innovative numerical methods to carry out uncertainty and sensitivity analyses. In this paper the work carried out within the DONUT work package is put in perspective regarding the existing concept and literature on the field. It does not pretend to be exhaustive but rather to put emphasis on particular issues tackled during the project.

KEYWORDS

numerical tools, numerical methods, upscaling, uncertainties treatment, machine
learning, coupled processes, reactive transport

1 Introduction

Demonstrating safety of a nuclear waste disposal over time comparable to geological timescales relies on a rigorous, complex and iterative scientific approach referred to as “long-term behavior science” that is based upon three pillars: experiments, modeling, and

natural/archaeological analogs (see Bildstein et al. (2019) for a review).
In that context, improvement of numerical methods and tools for
modelling multi-physical Thermo-Hydro-Mechanical-Chemical
coupled processes (THMC) occurring in radioactive waste disposal
is a major and permanent issue to support its optimization of design.
Indeed, numerical simulations are necessary either to make predictive

<figure><img src="image_4.png" /><figcaption>FIGURE 1
Thermal, hydraulic, mechanical, chemical processes concomitance affecting repository and geological environment during operational and post-closure period.</figcaption></figure>
multi-physical assessments, at time frames and space scales larger than experiments can cover, or to model the experiment itself to verify, calibrate and improve numerical models. These numerical simulations require integrating, in a consistent framework, an increasing scientific knowledge acquired for each of the individual components of a system for radioactive waste disposal. As described and detailed in Bildstein et al. (2019) and Claret et al. (2022) components of the so-called “multiple-barrier system” from the waste form to the biosphere include a combination of waste containers, engineered barriers, and natural barriers that will interact during the lifetime of the repository. The thermal, hydraulic-gas, mechanical and chemical processes and their timeline interactions are illustrated on Figure 1. This implies to consider couplings of different and non-linear processes, applied to a wide range of materials with contrasting properties as a function of time and space in ever-larger systems.

In this context, the development of cutting-edge and efficient numerical methods is thus necessary, in the scope of development and use of powerful and relevant numerical tools for assessments. In parallel, it is also necessary to manage the uncertainties associated to the input data that feed the models, and to the representation of the processes, in order to assess the range of variability of the results and to identify the main parameters and processes that drive the behaviour of the systems of interest. Managing uncertainties in these complex systems requires the improvement and the development of innovative, appropriate and efficient methods.

Within the EURAD joint programming initiative, the work package DONUT (Development and Improvement Of NNumerical methods and Tools for modelling coupled processes) ambitions to go a step further in the development:

* of relevant, performant and cutting-edge numerical methods that can easily be implemented in existing or new tools, in order to carry out high-performance computing to facilitate the study of highly coupled processes in large systems. These methods and their implementation in tools will be mainly applied to reactive transport, 2-phase flow, and THM modelling in porous and fractured media;

* of numerical scale transition schemes for coupled processes (meso to macro scale), supporting the study of specific multi-scale couplings such as chemo-mechanics;

* of innovative numerical methods to carry out uncertainty and sensitivity analyses;

In this paper the work carried out within the DONUT work package is put in perspective regarding the existing concept and literature existing on the field. It does not pretend to be exhaustive but rather to put emphasis on particular issues tackled during the project. The latter has involved 31 beneficiaries, consisting of Research Institutes, Universities, Technical Support Organization and Waste Management Organizations from 8 EU members countries and Switzerland.

## 2 Numerical methods for high performance computing of coupled processes

Development of new methods and new tools to efficiently solve numerically highly coupled problems arising in the context of waste

repository management depends on the coupled process to be tackled. After, we briefly comment on the coupled processes that have been considered in the project, we present the foundations of the numerical methods that have been developed as well as the existing tools (or codes) that have served as a basis for further developments.

### 2.1 Coupled processes

This section is devoted to the description of the physical models to be studied.

#### 2.1.1 Spatial coupling of gas free flows and two-phase porous media flows

Flow and transport processes within domains consisting of a porous medium and an adjacent free-flow regions find widespread applications in various industrial and environmental contexts. This is particularly evident in the context of radioactive waste disposal in deep geological repositories, where models of this nature become imperative for predicting mass and energy exchange at the interface between the repository and the excavated ventilation tunnels particularly relevant in the presence of microcracks. Typically, the porous medium, initially saturated with the liquid phase, undergoes drying due to suction (involved by tunnel ventilation) or crack flow in the vicinity of the interface. The microcracking network acting then as a vector for the propagation of diffusive flows. To effectively model such a complex physical process, it is essential to consider the flow of both liquid and gas phases within the porous medium, incorporating phenomena such as vaporization of the water component in the gas phase and dissolution of gaseous components in the liquid phase. Moreover, the model should potentially account for conversion of diffusive into free flow and vice versa throughout the microcracked layers and tunnel. Within the tunnel, a free single-phase gas flow can be considered, assuming instantaneous vaporization of the liquid phase at the interface. This single-phase gas flow needs to be compositional and non-isothermal to accurately reflect changes in relative humidity and temperature within the tunnel, factors that strongly influence the liquid flow rate at the interface.

If many works have been performed to model and discretize the coupling of single phase Darcy and free flows, there is very little work on the coupling of a two phase gas liquid compositional Darcy flow with a single phase compositional free flow. Such a coupled model has been recently proposed in Mosthaf et al. (2011) and Vanderborght et al. (2017) (see also Birgle et al., 2018; Beaude et al., 2019) using proper matching conditions at the interface between the porous medium and the free flow regions. This model is the starting point of our work and has been simplified to take into account the physical characteristics of our problem focusing on the drying processes at the interface between the nuclear waste repository and the ventilation excavated tunnels.

To obtain our simplified model, we assume that the longitudinal dimensions of the tunnels are large compared with their diameters allowing to reduce the model in the tunnel to a 1D free flow as in Brenner et al. (2016). In this work we extend this approach to account for a fracture network embedded in the matrix domain. The model couples the liquid-gas Darcy flow in the 3D matrix and 2D
fractures with the 1D free gas flow along the tunnel taking into account the liquid vaporisation and the convective boundary layer at the porous medium–tunnel interface. The discretization is based on a resolved interface with a non-matching 1D mesh along the tunnel centerline. It is applied to simulate the drying process between an operating tunnel and the radioactive waste storage rock with explicit representation of the fractures in the excavated damaged zone (see Brenner et al., 2022).

## 2.1.2 Capillary gas entry pressure and hysteresis

The study of the performance of the host rock and engineered barriers materials (bentonite plugs, backfill...) is crucial for radioactive waste confinement in a deep geological disposal (e.g., safety case of CIGEO project developed by Andra, France). Reliable models are needed for water and gas transport in a clay host rock formation (e.g., COx: Callovo-Oxfordian), and in any clay-based engineered materials used as barriers in the disposal, as well as their impact on radionuclide transport within the disposal. There is a substantial experimental evidence showing the existence of non-zero capillary gas entry pressure ($P_{c,ge}$) and hysteresis in hydraulic properties of these materials Andra (2012) but their modelling by exact mathematical formulations, and their impact on two-phase (water-gas) flow have never been studied. The impact of non-zero capillary gas entry pressure and hysteresis on maximum gas pressure and on desaturation range around the repository during the gas transitory period must be assessed.

In this work, we studied first the effect of capillary gas-entry pressure by implementing a new van Genuchten-Mualem (VGM) model for hydraulic properties in iTOUGH2 Finsterle (2007) and TOUGH2-MP (Zhang et al., 2008), called VGM-PE (Amri, 2021; Amri et al., 2022). This model has been first proposed for unsaturated single liquid water phase flow by Vogel et al. (2001) and Ippisch et al. (2006) and extended here to two-phase liquid-gas flow by incorporating a new function of the relative permeability to gas of any porous material. In addition, two methods for the simultaneous optimization of the parameters of the three hydraulic functions (water retention curve and relative permeability to water and gas curves) have been proposed (single objective and multi-objective minimization functions) and verified through fitting of COx-hydraulic properties data (Amri, 2021; Amri et al., 2022). Implementation of the VGM-PE model in iTOUGH2 and TOUGH2-MP needed some code modifications with numerical verification and experimental validation. Therefore, two new numerical approaches (PC-method and SG-method; Amri (2021); Amri et al. (2022)) have been developed for solving the two-phase flow problem incorporating VGM-PE. Both numerical methods were verified through a numeric experiment of 1D capillary rise of water in a vertical clay column proposed by Vogel et al. (2001). Results of the unsaturated water profiles by using either iTOUGH2/EOS3 or TOUGH2-MP/EOS3 code showed that the uncertainty on waterfront advance in the column is highly underestimated when $P_{c,ge}$ is neglected (relative error of more than 65%). Moreover, the experimental validation of VGM-PE at the scale of a borehole (dozen of meters), excavated in the Callovo-Oxfordian (COx) clay rock (ANDRA PGZ1 experiment for nitrogen gas injection in the COx - Bure Underground Research Laboratory; de La Vaissière (2011)), showed that both numerical methods based on either iTOUGH2/EOS3 or on TOUGH2-MP/EOS3 simulate

fairly well the gas ($N_2$) pressure measured within the chamber, and that the accumulation of gas within the chamber becomes higher when $P_{c,ge} \neq 0$. This leads to the increase of gas pressure and high desaturation of the EDZ and COx around the borehole at later times. After numerical verification and experimental validation of VGM-PE, a study of the contribution of a 3D HLW disposal cell (scale of some hundreds of meters) to the hydrogen transport within the backfilled access drift has been carried out to demonstrate the impact of $P_{c,ge}$ on the maximum gas pressure simulated within the backfilled drift. The 3D HLW disposal cell model geometry and parametrization are those proposed in EURAD-benchmark WP GAS, Task 4.2 (Levasseur et al., 2021). Simulation of this model by the SG-method fails to reach the maximum simulation time of 100,000 years. The numerical solution diverges earlier at 817 years accompanied by a gas pressure build up to 20 MPa, and a high desaturation of the drift backfill, the EDZ and host rock around the drift. However, simulations with $P_{c,ge} \neq 0$ assuming a linear approximation of the water retention curve near full liquid water saturation, as specified in the benchmark (Levasseur et al., 2021), were possible until 100,000 years and showed that the maximum gas pressure within the backfilled drift increases by about 3,5 MPa, i.e., from 9,6 MPa for parametrization $P_{c,ge} = 0$ to 13,1 MPa for parametrization $P_{c,ge} \neq 0$. This high increase of gas pressure is accompanied by a substantial drainage of water in the backfilled drift leading to a high gas saturation of about 0,623. The same observations have been made by Saàdi (2022) for his simulation of FORGE benchmark Wendling et al. (2019) model of hydrogen migration in a small repository “module” of 100 HLW disposal cells. Notice, however, that this gas pressure value can be higher than the lithostatic pressure of the host rock layer. Notice also that gas “piston effect” on water displacement in the drift network Saàdi et al. (2020) could be significant due to the high gas velocities simulated.

The second part of the study dealt with the problem of the two-phase flow with hysteresis and with $P_{c,ge} = 0$ (Amri et al., 2023). It consisted of a comparative study of five different conceptual hysteretic models (with and without gas entrapment) on data of water retention curves of five different materials (3 sands, bentonite, cement) collected from the literature, and on the other, of an application of these models to an intermediate laboratory scale experiment of Soltrol-220 (non-wetting fluid) injection (at a rate of 0,7 mL/s in the time interval $0 < t < 5,5$ h) in a sand 40/50 aquifer reservoir saturated by Glycerol-water (wetting fluid). This experiment has been carried out by Trevisan et al. (2014) to mimic supercritical $CO_2$ ($scCO_2$) trapping in an aquifer. Simulations of this experiment performed using iTOUGH2/EOS10 module (equivalent to T2VOC code developed by Falta et al. (1995)) with the new five implemented hysteretic models, showed that the conceptual hysteretic model with gas entrapment, Beriozkin and Mualem (2018), enhances better the simulation of scanning curves, and simulates well the distribution of Soltrol-220 plume in the aquifer during the post-injection times with less CPU time. This hysteretic model is promising for simulations of hydrogen migration at the large spatial scale of the repository for radioactive waste.

## 2.1.3 Thermo-hydro-mechanics for porous media flows

We shall start with the analysis of hydro-mechanical (HM) processes, which occur in rocks as well as in buffer and sealing materials in potential deep underground nuclear waste repositories.
The basic model is Biot's type poroelasticity, which couples the linear elasticity and Darcy's flow in saturated media. Even this simplest HM model shows that there are important questions about the robustness of the discretization (in space and time) and the efficient solution of the arising linear systems. Those points will be discussed in further detail in Section 2.2. Deep geological repositories of spent nuclear fuels (DGR) in many countries will be constructed in rock masses of the granite type. Such rocks are naturally fractured owing to their brittleness. The fracture network in rocks can create preferential paths for the transport of radionuclides and can be influenced by underground openings and related stress changes in rocks. This motivates us to consider HM models for porous media with fractures where fractures are considered as lower-dimensional objects. The simplest approach arises from the linear Biot model of poroelasticity where a weaker and more permeable material is considered for fracture representation. The corresponding mixed-dimensional model consisting of the rock matrix and the fracture network was derived in Březina and Stebel (2024) by an averaging over the fracture thickness and by prescribing appropriate interface conditions. Such a treatment extends previous work related to porous flow without mechanics, see e.g., Martin et al. (2005), Berre et al. (2019). Advanced HM models with fracture networks can account for nonlinear effects such as the dependence of fracture conductivity on the hydraulic aperture or the non-penetration condition prescribed on fractures. The first example is a model developed in Blaheta et al. (2020). Its flow part arises from Martin et al. (2005) and its mechanical part is based on linear elasticity in the rock matrix and frictionless contact on fractures. The second example is a model from Stebel et al. (2024) which extends the linear approach from Březina and Stebel (2024). Unlike Blaheta et al. (2020), displacements are also defined on fractures and thus a simplified linear friction can be considered there together with non-penetration conditions. A similar model, which incorporates contacts on fractures and a fracture growth model, relies on an alternative approach with degenerated finite elements, commonly known as cohesive zone modeling, as seen in Carrier and Granet (2012) and in R7.01.25 (2019). Other examples of nonlinear HM models on fractures can be found in Ganis et al. (2014), Garipov et al. (2016), Khoei et al. (2018), Settgast et al. (2017).

While cited above theoretical frameworks exist, in our opinion, there is still a lack of extended works dealing comprehensively with the continuum-fracture poroelastic model, particularly incorporating a consistent elasticity model on fractures Fontana (2022). This becomes especially crucial in the context of waste repositories, where the fracture domain is employed to describe larger fracture zones. These zones are nonempty and can be considered as a continuum with diminished strength.

Changes in temperature may also induce important effects that, often, cannot be neglected, requiring the use of thermo-hydro-mechanical (THM) coupled models. The research teams at Charles University and the Czech Technical University have developed and numerically implemented a coupled multiscale THM model suitable for describing the behavior of expansive soils under non-isothermal conditions Mašín (2017). In particular, the model is intended to be applied to bentonite buffers in deep geological nuclear waste repositories. The model uses the theory of soil hypoplasticity and a double-structure

approach Alonso et al. (1999) – backed by various micro-mechanical studies–in which soil macrostructure (silt-size aggregates of clay particles) and microstructure (internal structure of the aggregates) are treated separately. Four coupled components are thus described: a mechanical model for the microstructure, a hydraulic model for the microstructure, a mechanical model for the macrostructure, and a hydraulic model for the microstructure. These models are coupled at each structural level, and the behavior of the two structural levels is linked through the double-structure coupling function Mašín (2013). In particular, the mechanical behaviour of the macrostructure accounts for partial saturation, and its hydraulic response is based on a void ratio-dependent water retention model. The microstructure is considered fully saturated, with its mechanical behavior being reversible volumetric, governed by the Terzaghi effective stress principle. However, this simplification has proven to limit the model performance in some cases, such as during oedometric swelling, and necessitates improvements.

Each of the structural levels responds differently to temperature change, suction change, and mechanical action. The developed model thus considers a temperature dependence for all four model components. As for the mechanical behavior of the macrostructure, temperature influences the position of the normal compression line of a soil. The normal compression line represents the maximum void ratio of the soil at the given mean stress and temperature; therefore, a heating-induced shift in the position of the normal compression line implies heating-induced soil compaction of loosely compacted soil Mašín (2017). This effect is considered in the model by assuming that the asymptotic state boundary surface depends on temperature. Temperature also affects the mechanical response of the microstructure. It influences the basal spacing and, thus, both the amount of crystalline water and the diffuse double-layer thickness. Depending on conditions and clay type, the combinations of these effects can result in either swelling or contraction. This is accounted for in the model by a single parameter, calibrated through experiments.

In the water retention behavior, it is difficult to distinguish experimentally the response of the macrostructure and the microstructure Sun et al. (2022). Because the microstructure may be considered saturated up to large values of suction, temperature-induced microstructural swelling also implies an increase in water retention capacity, and microstructural shrinkage implies a decrease in water retention capacity. Therefore, a simplified approach is chosen, in which the effect of temperature on water retention capacity of the microstructure is considered to depend solely on the change in surface tension of water with temperature, without considering the effect of a variable grain-water interface contact angle Mašín (2017). Global water retention capacity is thus implied both by the retention of the macrostructure and by the specific volume of the (saturated) microstructure, which are linked through the double-structure coupling framework.

The model is included into the balance equations, where a variable number of components can be taken into account, i.e., grains, liquid water, vapour, dry air or their combinations. The grains and liquid water can be assumed compressible or incompressible.

In addition to the expansive behaviour, one of the other specific features of bentonite is the so-called oversaturation effect. It denotes
the phenomenon that bentonite can store more water than the amount corresponding to its pore space if the conventional value of the water density is considered (1,000 kg/m3). In Michalec et al. (2021), a THM model capturing this phenomenon has been proposed. It describes the water content within bentonite by the gravimetric water content, not by the usual degree of saturation related to filling the pore volume. This model has been validated on Full-scale Engineered Barriers Experiment (FEBEX) and successfully compared with other approaches, see Gens et al. (2021).

Another problem related to coupled models for bentonite is the thermodynamic consistency. Even though thermodynamic principles form one of the cornerstones of reliable mathematical modelling, they have often been disregarded in modeling of swelling clays. The thermodynamic consistency of HM coupling for the widely used Barcelona Basic Model has been analysed in Ligurský and Michalec (2023), and it has been found very problematic. A more promising conceptual framework for swelling together with its coupling to water retention processes and even chemical effects can be found in Pulkkanen (2019).

The study of geomaterial behavior have always been a significant concern for EDF as a major player in the energy sector. The critical importance of understanding how geomaterials respond in various industrial contexts, was initially driven by their impact on the integrity of infrastructure such as dams and nuclear facilities. More recently with emergence of various projects of underground waste management and consequent concerns about its environmental implications, the research domain was enriched by studies under extreme loading conditions. In this context, since their initial introduction in the finite element software code_aster (www.code-aster.org) approximately 25 years ago, the model for thermo-hydro-mechanical (THM) coupling has undergone continuous improvement in tandem with code evolution. The current open-source version of the code, that became HPC compatible, offers users a range of constitutive relations governing the thermo-hydro-mechanical behavior in both saturated and unsaturated media. The chosen modeling approach is grounded in the presentation of porous environments, particularly as elaborated by Coussy and Ulm (1995). Constitutive equations for the behavioural laws are derived from thermodynamic considerations, where each component of the porous environment—mechanical, hydraulic, and thermal—is treated as a state variable. Subsequently, relations describing hydraulic and thermal quantities are introduced, taking into account the robust couplings between these phenomena along with mechanical deformations through the hypothesis of Bishop effective stress Bishop (1959). The adopted scheme R7.01.10 (2019), R7.01.11 (2019) offers a significant advantage in that the resulting relationships can be coupled with any mechanical behavior law, providing a pathway for describing highly nonlinear phenomena. We highlight further two prominent classes of mechanical constitutive relations available in code_aster designed to represent rock formations and soils: Drucker-Prager R7.01.16 (2019) and Hoek-Brown R7.01.18 (2019) type. The introduction of both has some common background. First, geomaterials exhibit an asymmetric traction-compression curve and deform plastically under lower traction, displaying greater “elasticity” threshold under shear stress induced by hydrostatic compression. Consequently, the yield function is typically formulated with the assistance of the

hydrostatic stress invariant, which is asymmetric under stress reversal. Drucker-Prager-type constitutive relation R7.01.16 (2019) features a straightforward linear yield function and are commonly employed to represent rock/soil behavior under significant compression. However, these laws have limitations in situations where stress component values approach the traction loading domain, as observed in tunnel excavation. The LKR constitutive relation R7.01.24 (2019), Kleine (2007), Raude (2015) is a viscoplastic model developed by EDF to characterize the behavior of rocks, particularly callovo-oxfordian for the nuclear waste storage project Cigéo in East France. In the opposition to the Drucker-Prager constitutive relation, here the plasticity criterion is based on the quadratic Hoek-Brown yield function (R7.01.18, 2019; Hoek and Brown, 1980), allowing for a more accurate representation of the traction loading domain. Recent studies Di Pietro et al. (2022) have revealed a direct link between the Drucker-Prager and Hoek-Brown-like yield surfaces, making the introduction of the latter less artificial and more connected to the nonlinear nature of constitutive elastic media. In a similar vein, the reasoning can be extended to derive the cohesive zone constitutive relations for cracked media, preparing the ground for a unified framework describing both cohesive cracks and bulk geomaterials Fontana (2022).

The primary industrial application in the context of nuclear waste disposal management is associated with simulating tunnel excavations. The unique feature of this simulation lies in the complexity of the loading path of vault top of the tunnel which undergo the evolution from high tri-compression to bi-traction during the excavation process. Nuclear waste disposal aims to limit the diffusion of radionuclides, the priority should be given then to considering the physical phenomena of transport in cracked or highly plastified environments. Although the THM models and constitutive relations currently implemented in the code_aster are perfectly able to capture the main phenomena individually or at small scale, the large-scale simulation of fully coupled processes remains outside of scope due to its numerical complexity. These models concurrently face challenges of poor robustness and low numerical efficiency, particularly when dealing with active multi-scaling mechanisms of instationary diffusive front propagation. Addressing this research area should take precedence in future work.

### 2.1.4 Thermodynamically consistent modeling of corrosion

A mechanistic model of electro-chemical corrosion in deep geological storage of nuclear waste has been proposed in Bataillon et al. (2010). This model named Diffusion Poisson Coupled Model (DPCM) assumes that the carbon steel canister is covered by a dense spinel iron oxyde layer. The oxyde layer is in direct contact with the oxydizing medium or with a layer of corrosion products obtained by a precipitation mechanism, which can be modelled as a saturated porous medium, and with the undamaged steel. A fine description of the anaerobic corrosion process arising at the surface of the steel canister is necessary to predict both the evolution of the canister along time, as well as to predict the evolution of the chemical composition of the fluid in the porous medium (e.g., dihydrogen concentration) to be coupled with reactive transport models. While the corrosion process impacts the chemical composition of the fluid, the evolution of the fluid via reactive transport models impacts the boundary conditions of the
corrosion models which depend on the fluid composition, for instance through its pH.

From a mathematical point of view, it consists in a drift-diffusion-Poisson system of equations on a moving domain with strongly coupled boundary conditions. The main difficulty for the mathematical analysis comes from the moving boundaries whose velocity depends on the concentrations of the chemical species and on the electrical potential. Only results on simplified models have been obtained so far. Existence of solutions for some simplified models which do not take into account the moving boundary interfaces have been obtained in Chainais-Hillairet and Lacroix-Violet (2012), Chainais-Hillairet and Lacroix-Violet (2015) The existence of pseudo-stationary states is established in Chainais-Hillairet and Gallouët (2016) under some simplifying assumptions on the model.

In opposition to the simpler case of concrete carbonation (Aiki and Muntean, 2013, Zurek, 2019), no Lyapunov functional has been found so far for the DPCM model despite it is expected that the model fulfils Onsager's reciprocal relations which ensures that the physical energy decays along time (and is thus a Lyapunov functional). Besides, there has been recent advances in the mathematical description of models with generalized gradient flow structure (roughly speaking, model satisfying Onsager's reciprocal relation). Otto's seminal works (Jordan et al., 1998; Otto, 2001) has been a key step to understand the mathematical structure of generalized gradient flows where the evolution of the physical configuration relies on transport of chemical species, as it is the case for corrosion models. The Poisson-Nernst-Planck system has been shown in Kinderlehrer et al. (2017) to have the generalized gradient flow structure highlighted in Jordan et al. (1998). A first extension of this approach on moving domains was performed in Portegies and Peletier (2010) and extended to a toy model for corrosion in Merlet et al. (2022). Finally, a framework was recently proposed in Mielke (2011) and Mielke et al. (2016) allowing to go beyond Onsager's framework, i.e., when the relation between the forces and the evolution are no longer linear.

In the framework of the EURAD project the DPCM model has been updated to make it consistent with thermodynamics. The dissipation of the electrochemical Helmholtz free energy is guaranteed, first for a simplified two-species model set on a fixed domain (Cancès et al., 2023a), then for the full model on a moving domain (Cancès et al., 2023b). Moreover, the control of the free energy and of its dissipation rate allowed (Cancès et al., 2023a) to rigorously establish the global in time existence of a weak solution to the reduced two-species model. Such a mathematical analysis result for the full model presented in Cancès et al. (2023b) is not yet available so far. Let us further stress that the numerical simulations presented in Cancès et al. (2023a) for the reduced two-species model are better behaved than the methods based on the original DPCM approach as the new model preserves physical bounds on the charge carrier densities that could be violated in some specific configurations by the classical DPCM model.

The second axis where significant contributions have been proposed in on the design of numerical methods for drift-diffusion systems with vacancy diffusion. In such models, the convection due to the electric field is of nonlinear nature. Such a vacancy diffusion process has been adopted for iron and oxygen atoms in the magnetite in the last generation models (Cancès et al.,

2023a; Cancès et al., 2023b). The design of fast, accurate and thermodynamically consistent numerical methods for such vacancy diffusion as therefore needed. An extension to this context of the so-called *square-root approximation* (SQRA) finite volume scheme (Lie et al., 2013) has been proposed in *Cancès and Venel (2023)* for a single nonlinear Fokker-Planck equation. Extension to the more challenging context of multi-species diffusion accounting for cross-diffusion has been carried out in *Cancès et al. (2023c)*, where a very robust hybridization with the Scharfetter-Gummel scheme was further proposed.

## 2.2 Numerical methods

This section briefly described the state of the art on the numerical methods for the suitable approximation of the physical processes mentioned above.

### 2.2.1 Numerical coupling strategies for multi-physics processes

When one aims at simulating coupled multiphysics processes, different strategies can be relevant. One can either solve the full system all at once in what one would call a fully-coupled or fully implicit strategy, or one can solve each physical process apart and sequentially. This second approach should provide similar results as the fully coupled one if one iterates the sequential procedure until convergence, while allowing for a modular approach where each module is devoted to one physics.

Reactive transport modeling (RTM) plays a crucial role in the understanding of coupled thermal, hydraulic and chemical (THC) processes for the performance assessment of geological disposal of radioactive waste (Claret et al., 2018; Bildstein et al., 2019). In some situations of interest this context, such processes are governed by non-isothermal multiphase multicomponent flows with reactive transport in porous media. At the continuum scale, reactive transport models for porous media rely on constitutive equations that make it possible to describe the porous medium with respect to its macroscopically measurable properties such as permeability, dispersivity, diffusibility, and so on Steefel et al. (2015a). The significance of reactive transport constitutive equations and their treatment in saturated or non-saturated conditions have been recently reviewed by Claret et al. (2018). From the numerical point of view, modeling THC phenomena leads to a highly nonlinear coupled system of degenerate partial differential equations (PDEs) governing the compositional multiphase flow to algebraic or ordinary differential equations (ODEs) modeling respectively equilibrium and kinetic reactions. Several approaches can be found in the literature to deal with coupling between flow, transport, chemistry. The first approach is a fully coupled procedure and consists in solving a nonlinear system gathering all equations at each time step (Brunner and Knabner, 2019; Knodel et al., 2022). The alternative method follows a sequential paradigm and is commonly referred to as the Operator Splitting approach. The prevalence of the operator splitting method is evident in numerous numerical codes cited in the literature (Claret et al., 2018; Bildstein et al., 2019). Notably, several sequential codes, including HYTEC (van der Lee et al., 2003), ORCHESTRA (Meeussen, 2003), TOUGHREACT (Xu et al., 2012), CORE2D
V4 (Samper et al., 2012), Crunchflow (Steefel et al., 2015a), and OpenGeoSys (Kolditz et al., 2012), have found extensive application in reactive transport modeling for radioactive waste management. In these models, the operator splitting method enables the sequential solution of flow, reactive transport (or possibly, flow, transport, and chemistry) at each time step, potentially within an iterative loop. This approach allows for a modular programming structure in reactive transport modeling, offering significant convenience in integrating various physical and chemical processes. The modular structure facilitates the implementation of dedicated modules for hydrological and chemical aspects, optimizing the use of codes specifically designed for each side of reactive transport. A growing trend in the literature highlights the external coupling of such codes and the development of new reactive transport packages. For example, HYTEC, with its modular structure, incorporates the geochemical module CHESS for chemical processes (van der Lee et al., 2003). Many works employ the external geochemical package PHREEQC (Parkhurst and Appelo, 2013) for coupling with flow and transport codes. Instances include the coupling of PHREEQC with HYDRUS1 by Šimůnek et al. (2006) and its coupling with COMSOL MultiPhysics in a variable saturated and non-isothermal reactive transport problem by Nardi et al. (2014). The flexibility of PHREEQC in recognizing various instantaneous and kinetic reactions through basic scripts is a motivating factor for coupling it with specialized transport codes. Notably, PHREEQC has been successfully coupled with the open-source software OpenFOAM in continuum/pore scale simulations (Soulaine et al., 2022). Developers' efforts have led to the creation of PhreeqcRM, a C++ class with methods that integrate the full capacity of PHREEQC, making the coupling with object-oriented codes convenient (Parkhurst and Wissmeier, 2015). This class, identified in the code, enables a seamless integration of PHREEQC with transport simulators that use an operator-splitting approach. Within the DONUT project, a C++ version of the ORCHESTRA chemical solver, which has the full functionality of the original Java version, has been made available on GITHUB. No comparative study exists yet to quantify the accuracy loss for sequential approach in comparison with fully implicit approach, but its gain in implementation and saving in computing time are fairly obvious. However, sequential approaches introduce operator splitting errors and restrictions on the time step are mandatory to ensure mass conservation for instance. In the context of the simulator DuMuX (Koch et al., 2020), a sequential approach was developed and implemented for coupling two-phase flows and reactive transport, in a high performance computing context. More precisely, in Ahusborde et al. (2017), a sequential approach was developed to study the migration of hydrogen produced by the corrosion reaction in deep geological radioactive waste repository. This approach splits the global problem into two sub-problems. The first sub-problem computes a two-phase compositional flow where only species present in both phases are treated implicitly. The second sub-problem calculates a reactive transport problem where flow properties are given by the first step. Both subproblems were solved using a fully implicit manner (Ahusborde et al., 2018) but in a sequential procedure. However, some questions about the coupling between all the physical phenomena described above have also to be studied to assess if sequential approaches are suited to model very coupled phenomena or if fully implicit strategies would be more

relevant. For example, when simulating fast kinetics or strong chemical interaction, higher levels of coupling, implicitness may be useful as this can eliminate operator-splitting error and expedite nonlinear convergence. Obtaining solutions to these large systems of equations calls for a massively parallel implementation to reduce both computation time and memory requirements. In this context, in Ahusborde et al. (2019b), a fully implicit approach for a single phase multicomponent flow with reactive transport has been developed and validated via numerous 2D and 3D test cases including high performance computing. In Ahusborde et al. (2019a), the global implicit approach developed in Ahusborde et al. (2019b) was extended to deal with reactive two-phase flow for a one-dimensional test case proposed in Seigneur et al. (2018). This test case deals with porosity and permeability changes due to dissolution/precipitation dissolution/precipitation of minerals that can be linked to a process such as the atmospheric carbonation of partially saturated concrete. In Ahusborde et al. (2021), the fully implicit method has been used in a High Performance Computing context, with simulation using up to 768 processors. In Ahusborde et al. (2023), the sequential and implicit strategies developed in Ahusborde et al. (2018) and Ahusborde et al. (2021), respectively have been extended to non-isothermal flow to solve THC processes, still in a HPC context.

Besides, THM problems are solved in SIFEL (Koudelka et al., 2017; Koudelka et al., 2018) by two approaches: (a) a staggered approach, in which transport and mechanical parts are solved independently and in sequence (an iterative procedure is used to eliminate unbalanced forces); (b) a fully coupled approach, in which a complete thermo-hydro-mechanical stiffness matrix is assembled and solved using a Newton-Raphson iterative scheme in a single step. For both approaches, research is currently devoted to eliminating unwanted pore water pressure oscillations that arise in the model when full saturation is approaching.

Concerning the spatial coupling between free gas flows in the galleries and coupled liquid gas Darcy flows in the surrounding matrix and fracture network, domain decomposition appears to be a natural approach compatible with a modular implementation. A nonlinear domain decomposition algorithm based on the previous works (Masson et al., 2016; Bingle et al., 2018) has been designed in the spirit of Masson et al. (2020).

### 2.2.2 Linear and nonlinear algebraic solvers

Space and time discretization of HM and THM problems provides time-stepping numerical scheme with linear or nonlinear systems to be solved in each time step. These systems have a block structure corresponding to the selection of basis state variables. In the case of large problems, it can be more efficient to solve such systems iteratively by some Krylov space method with suitable preconditioner (Saad, 2003). A general way for construction of preconditioners begins with exploiting the block structure, keeping some diagonal blocks but replacing others by the corresponding Schur complements (Ferronato et al., 2010; Haga et al., 2010; Axelsson et al., 2012; White et al., 2016). Such preconditioners have to be subsequently modified by simplification of the Schur complements (Luber, 2022; Luber and Sysala, 2024). For implementation, subsystems arising from the diagonal blocks are usually solved by suitable inner iterations. Specific problems are the robustness of such techniques with
respect to oscillations and contrast in problem coefficients (Hong and Kraus, 2018; Hong et al., 2023; Luber and Sysala, 2024) when materials are heterogeneous. Similar preconditioners can be used for solving nonlinear system arising e.g., for flow in partly saturated media, when the problems are linearized by various Newton-type techniques (Bergamaschi and Putti, 1999). The stochastic character of the problems arising in many geotechnical applications brings another challenge to be solved, see Blaheta et al. (2020) for more details.

The solvers mentioned in the paragraph above can be interpreted as the so-called monolithic solvers of systems with block structure. As their alternative, one can also use block-splitting procedures separating hydro, mechanical and eventually other processes, see e.g., Mikelic and Wheeler (2013), Both et al. (2017), Březina and Stebel (2024), Almani et al. (2023), Hong et al. (2023). The block-splitting procedures are mostly combined with Uzawa's algorithms to guarantee the convergence. The splitting techniques are convenient if we solve HM problems with fractures and consider non-penetration contact conditions on the fractures (Stebel et al., 2024). In such a case, the problems do not lead to a block system of equations but to a system containing variational inequalities in mechanics. Therefore, the splitting techniques enable us to use state-of-the-art solvers on contact subproblems. For example, one can use the penalty method (Garipov et al., 2016) or methods based on Lagrange's multipliers and the corresponding duality analysis, see e.g., Franceschini et al. (2016), Franceschini et al. (2020), Franceschini et al. (2022) or Stebel et al. (2024) where the MPGP (Modified Proportioning with Gradient Projections) algorithm proposed in Dostál (2009) has been used.

In the numerical approximation of the THM problem provided by SIFEL (Koudelka et al., 2017; Koudelka et al., 2018), the construction of the stiffness matrix-to be consistent with the global Newton-Raphson iterative scheme for the adopted Runge-Kutta adaptive model integration-is carried out by two alternative approaches: (a) approximation using the linear part of the thermo-hydro-mechanical hypoplastic model; (b) numerical estimation of the stiffness matrix using perturbation. As neither of the approaches offers quadratic convergence because of various numerical problems, the convergence can be rather slow. Hence efforts are being made to find a reliable way to construct the stiffness matrix to obtain robust and fast converging simulations.

The Newton-Raphson method is used in the staggered as well as in fully coupled formulation. There are two implementations of the method. The classical method, where the system matrix is updated at every iteration, is usually more time consuming because of permanent matrix factorisation. On the other hand, the modified method uses the same system matrix for several steps, and the matrix factorisation is needed only when the matrix is updated. The modified Newton-Raphson method usually requires a shorter computing time. Despite the fact that modifications of the Newton-Raphson method often lose local superlinear convergence, they can be faster and simpler in some cases than Newton-Raphson method if they are combined with appropriate iterative solvers of linearized systems. It was demonstrated in Karátson et al. (2024) where the quasi-Newton method has been successfully combined with the deflated conjugate gradient method and algebraic multigrid.

Reactive transport modeling leads to a highly nonlinear coupled system of degenerate partial differential equations to algebraic or

ordinary differential equations. Due to the strong coupling of the flow and reactive transport equations, a specific attention was paid by many author to improve the convergence of the nonlinear and linear solvers, even more when considering fully implicit approaches. For instance, in Fan et al. (2012), after the Jacobian matrix assembly, a series of algebraic reductions (Schur complements), are performed to reduce the number of discrete equations to be solved. This strategy is commonly used in compositional reservoir simulations. In Brunner and Knabner (2019) and Knodel et al. (2022), the reduction technique described in Hoffmann et al. (2012) for reactive single phase flow was extended to the case of two-phase reactive flow. Unlike the classical strategy, their general transformation method does not only eliminate unknown equilibrium reaction rates. It also potentially reduces the nonlinear coupled part of the problem, allowing the use of large time steps and avoiding the potential drawbacks of sequential approaches.

### 2.2.3 Discretization methods

As concerns the discretization of HM problems, it was even very early observed that two-field (displacement-pressure) P1-P1 finite element (FE) discretization (i.e., both the displacement and the pressure are approximated by piecewise linear fields on simplex meshes) is not robust and can produce spurious oscillations due to a locking (e.g., Lewis and Schrefler, 1998). P2-P1 elements (i.e., the approximate displacement field is now piecewise quadratic) help but they are more expensive in computations. The reason for these effects is nowadays explained by inheriting saddle point structure. Later, three-field (displacement-velocity-pressure) formulation was used with discretization e.g., by P1-RT0-P0 elements (RT0 denotes the lowest order Raviart-Thomas element, and P0 denotes piecewise constant fields). This discretization provides conservativeness for flow, but it was recently shown that a fully stable discretization with the lowest order FE for elasticity requires abandoning conformity and using the stabilized Crouzeix-Raviart elements or the discontinuous Galerkin formulations, see Yi (2013), Frigo et al. (2021), Luber (2022), Luber and Sysala (2024). It is also possible to discretize elasticity by the mixed FE and use four-field formulation, see Bean et al. (2017). For higher-order time discretization of the Radau type and its stability, we refer to Axelsson et al. (2015a, 2015b).

Finite Volume methods are natural to discretize physical balance laws in a locally conservative way. Mathematical analysis tools have been developed in the last decades to assess their robustness and accuracy (see for instance Eymard et al. (2000)). It is now well understood that the Scharfetter-Gummel numerical fluxes (Scharfetter and Gummel, 1969), dedicated to the numerical approximation of convection-diffusion fluxes, allow to capture in an accurate way the long-time behaviour of drift-diffusion systems (Chatard, 2011; Bessemoulin-Chatard and Chainais-Hillairet, 2017). This motivated the choice of Scharfetter-Gummel based fluxes in Bataillon et al. (2012) for the simulation of the DPCM corrosion model (Bataillon et al., 2010). These particular fluxes were also used in the numerical investigations on reduced models (Chainais-Hillairet and Lacroix-Violet, 2015; Chainais-Hillairet and Gallouët, 2016; Zurek, 2019; Cancès et al., 2023a).

In the framework of the DONUT workpackage, the design of numerical methods for drift-diffusion systems with vacancy
diffusion has been carried out. In such models, the convection (for instance due to some electric field) is of nonlinear nature. Such a vacancy diffusion process has been adopted for iron and oxygen atoms in the magnetite in the last generation models (Cancès et al., 2023a; Cancès et al., 2023b).

The design of fast, accurate and thermodynamically consistent numerical methods for such vacancy diffusion as therefore needed. An extension to this context of the so-called *square-root approximation* (SQRA) finite volume scheme (Lie et al., 2013) has been proposed in Cancès and Venel (2023) for a single nonlinear Fokker-Planck equation. Extension to the more challenging context of multi-species diffusion accounting for cross-diffusion has been carried out in Cancès et al. (2023c), where a very robust hybridization with the Scharfetter-Gummel scheme was further proposed.

When a system of partial differential equations is cast under a gradient-flow structure, this gradient-flow structure suggests canonical time-discrete approximations and provides a context for the design of numerical methods which preserve structural properties Minimizing movement schemes (Jordan et al., 1998) provide naturally stable time-discretizations of generalized gradient flows. They can be seen as generalizations of the Backward Euler scheme and are now common for studying generalized gradient flows (Ambrosio et al., 2008). The exploitation of such fully discrete minimizing movement schemes for numerical purpose is quite recent and goes back to Benamou et al. (2016). Then several improvements have been proposed to improve the efficiency of such approaches (Cancès et al., 2020; Li et al., 2020; Carrillo et al., 2022).

Most of the methods that have been designed in order to maintain their thermodynamical consistency require some strong regularity assumptions on the mesh that are not suitable in the context of porous media flows in complex geometries. Such frameworks motivated the development of advanced finite volume methods, among which the nodal Vertex Approximate Gradient (VAG) scheme introduced by Eymard et al. (2012). The discretization of the 3D non-isothermal compositional liquid gas Darcy flow is based on VAG scheme (Xing et al., 2017). The coupling with the 1D free gas flow is based on non-matching meshes between the 3D-2D Darcy flow domain and the 1D domain in order to account for complex geometries of the tunnels (Brenner et al., 2022).

Another type of discretization method, for reactive transport simulations is the lattice Boltzmann method (Higuera et al., 1989; Qian et al., 1992). It is a special discretization of the kinetic Boltzmann equation with applications of LBM range from diffusion dominated reactive transport simulations in porous media to multiphase fluid dynamics of turbulent flow (Succi, 2001; Kang et al., 2002; Chen et al., 2003). Recently, it has been implemented to study the flow properties of a compacted illite digital rock (Pazdniakou and Dymitrowska et al., 2018), as well as for the cross scale modelling of dissolution and precipitation (Fazeli et al., 2020; Lönartz et al., 2023). Implementation of the classical nucleation theory allows to trace the formation of first critical nuclei up to the micrometer scale as was shown for the case of barite and Celestine precipitation (Prasianakis et al., 2017). Models that account for nucleation and growth of particles are very sensitive to microscopic parameters like the surface tension and the size of first critical nuclei (Noguera et al., 2010). An application to

radioactive waste disposal and more specific the Radium-bearing barite precipitation has been recently presented (Curti et al., 2019). Upon upscaling the pore-level algorithm can be used to investigate the coupled advection-diffusion-reaction-precipitation processes (Churakov and Prasianakis, 2018). The precipitation and the dissolution of minerals has a direct effect to the mass transport and vice versa (Prasianakis et al., 2018; Poonoosamy et al., 2019). A decrease in the porosity will result to increased resistance to flow, which under a hydraulic gradient could lead to increase in stresses. A direct coupling of lattice Boltzmann simulators for fluid flow with mechanics simulators still remains a challenging task.

Despite the fact that modifications of the Newton-Raphson method often lose local superlinear convergence, they can be faster and simpler in some cases than Newton-Raphson method if they are combined with appropriate iterative solvers of linearized systems. It was demonstrated in Karátson et al. (2024) where the quasi-Newton method has been successfully combined with the deflated conjugate gradient method and algebraic multigrid.

Capillary effects play a crucial role for the migration of multiphase fluids in heterogeneous media (van Duijn et al., 1995; Cancès, 2010; Bourgeat et al., 2013). Discontinuities of the conservative variables caused by spatial heterogeneities make it necessary to design appropriate strategies concerning the discretization (Hoteit and Firoozabadi, 2008; Cancès, 2009; Bastian, 2014; Eymard et al., 2014) and the nonlinear solvers (Wang and Tchelepi, 2013; Brenner et al., 2017). As proposed in Hoteit and Firoozabadi (2008) and Bastian (2014), a natural choice for the primary variable is the capillary pressure rather than the saturation, requiring adaptation of the existing tools like TOUGH2-MP (Pruess et al., 1999). These adaptation are already at the basis of the recent contribution (Amri et al., 2020). Future work will be dedicated to the validation of this new version of TOUGH2-MP on PGZ-experiment by adding an adequate hysteresis model in the compositional two-phase flow and transport equations.

### 2.2.4 Error control

In scientific investigations where precision is paramount, such as in the study of complex physical processes or the design of engineering structures, error control plays a pivotal role. It provides a means to quantify and manage uncertainties, reinforcing the confidence in the predictive capabilities of numerical simulations. Additionally, the incorporation of error control measures facilitates the optimization of computational efficiency, guaranteeing the computational resources are allocated judiciously to areas where refinement is most warranted. Ensuring the convergence of nonlinear numerical simulations for large-scale structures is typically a challenging task. In practice, when faced with difficulties, one might consider relaxing the convergence criteria while refining the loading steps. This fundamental engineering approach becomes more intricate when dealing with HM (or THM) simulations, cited above. The complexity arises from the multi-physical coupling phenomenon, where relaxing the convergence criteria inevitably results in alterations to the physical outcomes. While it is relatively straightforward to make reasonable guesses about the validity of results in pure mechanical or hydraulic simulations, coupling introduces a nearly impossible mission due to the intricate interplay of all phenomena involved. One might achieve physically relevant results from a mechanical
standpoint, such as structural deformations, while simultaneously obtaining inaccurate outputs for pressure or hydraulic flux distribution or vice versa. This discrepancy is primarily attributed to the fact that, in standard finite elements code, sources of convergence errors—be they hydraulic, mechanical, or temporal refinements—are combined to generate a single scalar stopping criterion. The proposed solution involves disentangling these errors to facilitate a more nuanced analysis. Furthermore, a heuristic algorithm could be implemented, automatically halting the code execution once the optimal convergence for the given mesh is achieved, thereby providing a more software agnostic experience for finite elements code users here are two predominant methodologies available for tackling the error estimation problem, presented above. These distinct approaches offer varied strategies and perspectives, each with its own set of advantages and considerations. The first, the classical “a priori” one has a drawback of relying on the problem regularity, which is not known a priori for most of the industrial applications. On a contrary “a posteriori” error analysis (see for example, Ainsworth and Oden, 2000) seems to be more relevant as the solution regularity is “reconstructed” based on the current finite elements solution. In most recent developments, one reconstructs the associated fluxes (stress for displacement, hydraulic flux for pressure) which dispose some important mathematical properties (smooth elements transition, symmetry, see Riedlbeck et al. (2017); Di Pietro et al. (2022)). The principal advantage of the following technique is the easy comparisons of all error’s sources. The last property leading to the possible straight implementation of different stopping algorithms. In the last decade a posteriori error estimation, based on flux reconstruction on patches, has been implemented in the various open-source finite elements software, both in academic, FreeFem++<fnref n="1" /> and industrial, code_aster<fnref n="2" /> contexts. In the latter case the ultimate aim is to assess the efficiency of the current technique for industrial applications. Encouraging results have been obtained, particularly in a tunnel excavation-like test case, where the error estimation closely approaches perfect efficiency compared to a reference study conducted on a highly refined mesh (Riedlbeck et al., 2017; Botti and Riedlbeck, 2018). A more recent results were obtained through the extension of a posteriori error estimation techniques to account for contact non-linearity, as described in Di Pietro et al. (2022) and Fontana (2022). During a tunnel excavation process, multiple contact interfaces can emerge, either due to construction geometry or through time evolution caused by crack propagation, progressive material softening, and plastification. It is essential to consider contact and potentially friction between two bodies to provide a more realistic description of the interface zones in civil engineering construction. The combination of both errors estimations approaches cited above should logically lead to the relevant criterion for multi-physical radioactive waste management applications. However, the current obstacle hindering the extension of a posteriori error estimation strategy to large-scale calculations is the notably high computational cost. Two primary reasons contribute to this perceived “poor” numerical

performance. The first one is typical software/hardware entanglement bottleneck, which is linked to the initial sequential programming scheme that generated the data-processing architecture incompatibility, specifically with sub-optimal memory usage. The second reason is more profonde and is related to the presence of spurious oscillations (Escoffier and Kazymyrenko, 2023; Kazymyrenko et al., 2023) in the discretized solution within the cohesive zone region. In the case of surface breaking cracks, the geometric singularity on its extremity reduces the regularity of the continuous solution affecting then discretized model quality. The oscillations contribute to the observed issues with the robustness of highly non-linear interface modelling under mixed compressive loading conditions, typical in underground excavation simulations, resulting in erroneous flux reconstruction close to the crack boundaries. While the issue was noted across most common industrially deployed discretization schemes, it could potentially be mitigated by selecting the appropriate finite elements to represent crack interfaces. However, as of now, a definitive solution remains elusive. It is imperative to identify and tackle these challenges in order to surmount current limitations and enhance the effectiveness of error estimation in large-scale industrial simulations.

Finally, we conclude that the convergence of the nonlinear and linear solvers can be improved using a *posteriori* estimators to achieve computational savings by stopping timely the linear and nonlinear solvers (Vohralík and Wheeler, 2013; Riedlbeck et al., 2017; Ben Gharbia et al., 2019), but some technical locks need to be addresses before these procedures being parallelised in order to meet industrial needs in numerical simulations.

### 2.2.5 Going parallel

The relevance of parallelization for numerical simulations lies in its capacity to dramatically enhance computational efficiency by distributing tasks across multiple processors units, thereby accelerating the solution of complex physical models and facilitating the analysis of data set variations. The development of more powerful and energy-efficient computer architectures is widening the use of large scale high performance computing techniques. Notable increases in computational power are being achieved by many core systems with powerful accelerators such as general-purpose graphics processing units (GP-GPUs). These systems also feature increasingly intricate hierarchical cache and memory designs, enabling the provision of the necessary data flow essential for computing units. These advanced hardware configurations present challenges for implementing and scaling current simulation algorithms while also presenting opportunities for the creation of new ones. In 2022, the world's first public exascale computer was announced, it has achieved the long anticipated peak performance of one exaflops $10^{18}$, representing a 100-fold enhancement compared to supercomputers. This advancement stems not from faster processors, but from a significantly larger number of processors/cores. The current breakthrough is poised to significantly impact large-scale computational science and engineering. Parallel algorithms will need to demonstrate significantly higher levels of concurrency while maintaining robust scalability. When addressing transient problems, where information progresses forward in time, sequential exploitation has traditionally been favored. However, the vast potential for

<footnote n="1">1 https://freefem.org/</footnote>

<footnote n="2">2 www.code-aster.org</footnote>
parallelism undoubtedly motivates a shift in this paradigm. One of the key drivers behind harnessing greater levels of parallelism will be to reduce time-to-solution.

All alternative approaches to parallelizing numerical simulations involve task decomposition, wherein computational tasks are partitioned into smaller subtasks capable of concurrent execution across multiple processors or cores. Several iterative methods have been investigated extensively to date. Lions et al. (2001) employed time discretization of a partial differential equation (PDE) that lends itself well to parallel implementations. The method, based on an Euler scheme, combines coarse resolutions and independent fine resolutions in time in the same spirit as standard spatial approximations. Emmett and Minion (2012) also applied a parallelization of numerical methods for PDEs in the temporal direction. The method was iterative with each iteration consisting of deferred correction sweeps performed alternately on fine and coarse space-time discretizations. With regard to direct methods, time-parallel methods can be found in Gander (2015). In general these methods can exploit low levels of concurrency. Moreover, Badia and Olm (2018) proposed a parallel-in-time solver for ordinary differential equations that relies on the well-known Schur complement method in linear algebra. Chen et al. (2019) used multipoint flux approximation (MPFA) method to discretize the spatial variables on staggered grid in order to accurately describe the flow behaviour in the porous media whose permeability is anisotropic, in which case the permeability tensor is full. For staggered grid in 2-D space, pressure values were defined at the centre of the cells and velocity components were located on edges. They were considered the permeability fields which contain different anisotropy angles in the simulation and fractures are represented explicitly by volumetric grid cells and the numerical algorithm is parallelized in order to reduce the simulation time. The algorithm is implemented in parallel using FORTRAN and the linear system arising from the pressure equation is solved by the aggregation-based algebraic multigrid solver (Notay, 2010; Napov and Notay, 2012).

As previously discussed in section 2.1.2, the computation of HM or THM models typically employs a time-stepping scheme involving the solution of linear or nonlinear systems at each time step. While the time-stepping process itself is sequential, opportunities for parallelization emerge during the solution of these systems. The first possibility is to use preconditioners based on approximations of the Schur complement systems (Luber, 2022; Luber and Sysala, 2024). It enables us to solve multi-physical processes in parallel. The second possibility is to apply domain decomposition methods and related parallel solution on subproblems arising from the splitting of the multi-physical system. For example, one can use the Schwarz domain decomposition methods (Blaheta et al., 2010; Blaheta et al., 2018b) or the FETI-based domain decomposition methods (Hapla et al., 2016; Kružík et al., 2020a; Dostál et al., 2022). Another domain decomposition method reflecting fracture propagation was suggested in Settgast et al. (2017). In addition preconditioners may incorporate inner iterative solvers, for which efficient parallel solution methods are readily accessible (Ordonez et al., 2023; Ordonez Egas, 2022).

Utilizing the finite element method for solving coupled thermo-hydro-mechanical problems results in a significant quantity of unknowns due to the definition of displacements, pore pressures,

and temperature at each node. Consequently, the total number of unknowns in the algebraic equation systems surpasses that of systems corresponding to purely mechanical problems, which only discretize displacement. Additionally, the high number of unknowns in each node result in a very broad bandwidth of the system "rigidity" matrix. These factors necessitate parallel processing, particularly for large engineering structures. Parallelisation of a finite element code can be done efficiently with the help of either domain decomposition methods (DDM) or alternatively splitting method. In both cases the resolved domain is split into smaller subdomains which are treated independently on particular processor. In DDM method each subproblem is solved independently and the continuity of data on the borders of subdomains is enforced during supplementary operation. In splitting method, the decomposition is used mainly in order to reduce memory consumption due to parallelisation of elementary calculation (element by element in finite elements method calculation), the total stiffness matrix is then assembled as a whole and finally the global algebraic equation is then solved using second level parallelisation relied on the other parallelisation schemes.

With respect to highly nonlinear and nonsymmetric systems of algebraic equations, the most suitable DDM is the Schur complement method, which is a variant of the Gaussian elimination. Splitting method may be less efficient then the DDM scheme but it allows more versatile usage while the second level parallelisation of linear algebraic system is completely independent. For the latter task one could rely on external expertise such as know-how of worldwide community of software developers proposing a highly stable and scalable open source solvers as MUMPS MUltifrontal Massively Parallel Sparse direct Solver (Amestoy et al., 2001; Amestoy et al., 2006) or PETSc Portable, Extensible Toolkit for Scientific by Computation (Balay et al., 2019).

Since 2019 all new versions of Code_Aster are distributed as fully parallelizable. The adopted parallelization is the splitting method described above. While all elementary (in the sense of finite elements method) calculations are naturally parallelizable, the parallelization of a *posteriori* errors necessitates in deep comprehension of Code_Aster architecture. Error estimator that seems to be suitable for hydro-mechanical coupling problems described in the recent works of Riedlbeck et al. (2017) and Botti and Riedlbeck (2018) is based on patch reconstruction scheme. It means that each elementary calculation is not a simple mesh element integral as in the standard finite element procedure, but a patch calculation, patch being an agglomeration of mesh elements connected to the given node. The parallelization scheme adopted in Code_Aster must then be modified in order to allow new kind of elementary calculation be more efficient. Given the usually rigid architecture of industrial finite elements code the full parallelization of a *posteriori* error estimation is still a rather challenging task. The ORCHESTRA chemical solver is basically a module that operates on cells with a given total chemical composition. Internally this solver module does not make use of parallel processing, and it has no spatial or discretization issues as it operates on single spatial points in a transport domain. However, it was designed to be used in the form of multiple parallel instances that can operate in parallel on large groups of cells/nodes. The ORCHESTRA chemical solver is well suited for this type of parallelization, as it has a very small
memory footprint, and is prepared for the existence of multiple
instances in a single memory space. Furthermore, the module is
designed to make efficient use of the chemical composition
calculated in each cell during a previous time step by using this
condition as start estimation in the iteration procedure. This greatly
reduces the number iterations, and thus calculation time, in
transport systems where chemical changes are concentrated in
specific spatial zones.

The coupled liquid gas Darcy and free gas flows is implemented
in the open source parallel code ComPASS (Xing et al., 2017; Lopez
et al., 2018; Beaude et al., 2018). The COMPASS parallelization is
based on mesh partitioning and MPI communications. It uses the
Petsc and Hypre libraries for the solution of the linear systems. Since
the size of the 1D free gas flow model is negligible compared with the
size of the 3D Darcy flow model, the domain decomposition
approach is adapted to a simple parallel implementation of the
coupled model.

Another numerical method which exhibits a significant
parallelization advantage is the lattice Boltzmann method. The
method operates typically in a regular cartesian grid, and the
operations are local at each grid point, requiring information
only from the neighboring computational grid points. This allows
to decompose the computational domain and split the
computational costs in an efficient way among several computing
units. Moreover, the algorithm shows significant acceleration when
deployed in GP-GPU architectures with a speed up of one to two
orders of magnitude compared to a single CPU core. The
parallelization efficiency is maintained also in hybrid CPU/GP-
GPU supercomputer architectures with the weak scaling being
practically linear (Latt et al., 2021; Krause et al., 2021; Safi et al.,
2017; Bauer et al., 2021).

## 2.3 Available tools

Let us describe the available tools that will be at the basis of the
developments carried out in the context of the project.

CALIPSO is a Fortran90 code for the simulation of corrosion of
iron based alloys relying on the discretization presented in Bataillon
et al. (2012) of the DPCM model (Bataillon et al., 2010). Its
development was initiated by ANDRA in 2008. A stabilized
version has been registered in 2014 at the French Agency for
Program Protection<fnref n="3" /> under reference IDDN. FR.OO 1. I
60009.000.s.P.20 I 5.000.30000. It is currently maintained by
Christian Bataillon (CEA) and Laurent Trenty (ANDRA).

COMPASS<fnref n="4" /> is an open source parallel code co-developped by
the partners of the ANR CHARMS project BGRM, UCA-CNRS-
Inria-LJAD, MdS, Storengy and LJLL (Xing et al., 2017; Beaude et al.,
2018; Lopez et al., 2018). It has the following features:

* Physics

* Multiphase compositional thermal Darcy flow model

* 2D discrete fracture or fault network coupled with the
  surrounding 3D matrix

* Multi-branch well model

* Numerical methods

* Parallel Vertex Approximate Gradient (VAG) scheme on polyhedral meshes accounting for the flow in the 3D matrix and the 2D fault network.

* Fully implicit Coat's type formulation of multiphase compositional thermal models.

* GMRES or BiCGStab linear solver with CPR-AMG preconditioner

* Computer science

* SPMD Paradigm

* Fortran 2003 + C/C++ + MPI

* Libraries Metis, Petsc + Hypre, vtk, HDF5

* High level python API

ORCHESTRA was originally developed as a Java code that was
specifically designed to allow the development new chemical
adsorption model types and to make it possible to combine these
chemical models with reactive transport processes. It is available
from orchestra. meeussen.nl. Its unique property, in contrast wil all
other available chemical solver codes, is that all chemical and
physical model code is defined in object oriented text files
outside the executable code. These text files are read as input
during runtime and can be extended or adapted by end users.
This implies that users can change the chemical model library
without the need to change the solver code, and the other way
around that the solver code can be changed without affecting the
chemical model capabilities. This approach/structure, that has been
described in Meeussen (2003), has resulted in a solver calculation
kernel that is relatively small and simple, but also extremely fast. The
simplicity of the kernel allowed it to be translated within this project
from Java to a C++, which makes the solver much easier to use from
other, e.g., transport codes, or from python. Because the C++ version
can read the same input files as the Java version, the complete
chemical model library (text file) is automatically available for both
the Java and the C++ version. So Java and C++ version have the
same chemical model capabilities. The (open) source code is
available on github: hansmeeussen/orchestraCPP (github.com).

The chemical model library contains a large set of chemical
models, including a large number of state of the art adsorption
models, and has proven to be very calculation efficient in reactive
transport benchmarks (Steefel et al., 2015b). Modelling of multi-
component transport is also implemented in the ORCHESTRA code
which is able to calculate multicomponent transport of individual
dissolved species taking into account the effect of developing electric
potentials. The ORCHESTRA approach additionally allows to
combine diffusion and advection, and can take into account
multiple pore domains with different diffusion properties (e.g.,
Donnan and surface (Meeussen, 2023; Jenni et al., 2021). As the
framework is written in Java which allows ORCHESTRA to run on
different operating systems, and the currently developed C++
version makes it relatively easy to combine the solver module
with C++/FORTRAN code. The module is fully prepared for
parallel processing, and can be used in a form with multiple
independent solver instances. The ORCHESTRA solver module,
both in Java and C++ form, has been extensively tested within the

<footnote n="3">³ https://www.app.asso.fr/</footnote>

<footnote n="4">⁴ http://www.anr-charms.org/page/compass-code</footnote>
Machine Learning Benchmark exercise in which it demonstrated to give results consistent with PHREEQC and GEMS calculations.

DuMu<sup>X</sup>: For THC numerical modeling, several sequential (Ahusborde et al., 2017; Ahusborde et al., 2018) and fully implicit (Ahusborde et al., 2019a; Ahusborde et al., 2019b; Ahusborde et al., 2021; Ahusborde et al., 2023) modules have been developed and integrated in DuMu<sup>X</sup> framework. DuMu<sup>X</sup> (DUNE for Multi-{Phase, Component, Scale, Physics, ...} Flow and Transport in Porous Media) is a free and open-source object-oriented software written in C++ (Koch et al., 2020). The environment is equipped with efficient solvers and massively parallel computation capability provided by DUNE (Bastian et al., 2008a; Bastian et al., 2008b).

TRIAX–The coupled thermo-hydro-mechanical model for expansive clays is implemented in an in-house element test driver developed since 2001 at Charles University and elsewhere by David Mašín and co-workers. The driver covers a wide range of constitutive models and testing conditions. Virtually, any element test can be simulated using the procedure generated by Janda and Mašín (2016); the majority of standard and less-standard testing conditions are pre-defined (including complex thermo-hydro-mechanical tests in unsaturated conditions, rate-dependent loading, plotting response envelopes, hypoplastic asymptotic state boundary surfaces, and elasto-plastic yield surfaces). Loading stages can be combined to form complex loading paths. The driver can also run any constitutive model implemented in umat format. The driver can be freely downloaded at https://soilmodels.com/triax/. The download package includes the source code, pre-compiled Windows and Linux executables, and a wide range of input files associated with various publications.

ExCalibre is an online tool (freely accessible at https://soilmodels.com/excalibre/) which enables users to automatically calibrate advanced soil constitutive models (Kadlíček et al., 2018). It adopts a unique model-specific calibration algorithm which remains consistent with the physical meaning of material parameters rather than using a blind parameter optimisation. It thus reproduces the thinking of an experienced engineer calibrating the model, leading to a reliable parameter set to be used in geotechnical simulations. After calibration, the user can manually change the parameters to fine-tune the calibration or to investigate the effect of individual parameters on model predictions. A simple Excel spreadsheet based on the provided template is used as input. In addition to automatic calibration, ExCalibre also acts as an element test driver simulating standard laboratory tests.

SIFEL (Simple Finite Elements) is an open-source finite-element code that has been developed since 2001 at the Department of Mechanics of the Faculty of Civil Engineering of the Czech Technical University. SIFEL comprises several components which can be used as independent computer codes, or can be linked to obtain additional functionalities. All parts are released under a GNU General Public Licence (GPL) and can be downloaded at http://mech.fsv.cvut.cz/sifel/. Within the SIFEL package, in particular, MEFEL is the component which deals with mechanical problems, such as linear and non-linear statics, eigen-vibration analysis, forced dynamics, time dependent problems of rheology. The coupled transport of heat, moisture and other species can be modelled in the part TRFEL. Namely, it contains heat transport, moisture transport, coupled heat and moisture transfer,

moisture and salt transfer. Coupled thermo-hydro-mechanical models, including models for expansive clays, are implemented in the part METR which interacts with the MEFEL and TRFEL. All parts of the code are equipped with one, two and three-dimensional finite elements and material models. The material models can be isotropic, orthotropic or generally anisotropic.

Code_Aster: Today, all the poro-mechanical models used in EDF are capitalized in Code_Aster (EDF code in Open-Source: www.code_aster.org). Code_Aster is integrated in multifunctional finite elements platform Salome_Meca. Code_Aster offers a full range of multiphysical analysis and modelling methods that go well beyond the standard functions of a thermo-mechanical calculation code: from seismic analysis to porous media via acoustics, fatigue, stochastic dynamics, etc. Its modelling, algorithms and solvers are constantly under construction to improve and complete them (1. 200.000 lines of code, 200 operators). Resolutely open, it is linked, coupled and encapsulated innumerable ways. Code_Aster provides full thermo-hydro-mechanical coupling framework for stationary and transient calculations. It relies on highly efficient external linear algebra solvers. It worth to mention that the software is under quality assurance QA (independent validations, reference of 2. 000 test cases, 13.000 pages of documentation, source management, qualification of version, etc.).

GEM-Selektor: (http://gems.web.psi.ch/GEMS3/) is a graphical user interface of GEMS3K solver of chemical equilibria (Wagner et al., 2012; Kulík et al., 2013). It greatly facilitates configuring the chemical system, the input of composition recipes, the computation of single equilibria or process simulations, viewing/exporting the results as tabulated data or plots, and managing the thermodynamic data. Chemical system definitions can be exported per mouse click into GEMS3K i/o files for use in coupled codes such as GEMSFITS and reactive-transport simulators Comsol-GEM, OpenGeoSys-GEM and CSMP++GEM. Simplified 1-D reactive transport models can be performed directly in GEMS-Selektor using its GEM2MT module. GEMS is open-source software, available at https://bitbucket.org, with binaries distributed from http://gems.web.psi.ch as freemium together with the built-in PSI-Nagra and SUPCRT98 chemical thermodynamic databases. In addition, third-party thermodynamic databases (e.g., Cemdata'18, HERACLES, Mines'19), maintained by their owners (https://www.empa.ch/cemdata; http://tdb.mines.edu), are available as plugins for GEM-Selektor. The GEMS Development Team currently involves 11 members from 6 institutions, and strives to implement innovative concepts, modern algorithmic frameworks, and tools to improve thermodynamic data, all in order to ensure the state-of-the-art functionality of GEMS for the next decade.

PERMON PERMON (http://permon.vsb.cz/) is a set of solvers combining quadratic programming and domain decomposition methods. It makes use of and extends the PETSc framework (https://www.mcs.anl.gov/petsc/) for numerical computations. The core of PERMON is formed by the PermonQP package which is able to solve large scale quadratic programming (QP) problems. Among the main applications are contact problems of mechanics (Kružík et al., 2020a; Kružík et al., 2020b) and Support Vector Machines (SVM) machine learning (Pecha et al., 2023). Contact problems can be decomposed by means of the TFETI (Total Finite Element Tearing and Interconnecting) non-overlapping domain decomposition method. Within the EURAD project,
PERMON was combined with the Flow123d library (see below) in order to solve large-scale HM problems for porous media with fractures where non-penetration contact conditions were considered (Stebel et al., 2024).

Flow123d (https://flow123d.github.io/) is a versatile simulator of transport and mechanical processes in fractured porous media. The novelty of this software stems from its comprehensive support for computations on mixed meshes comprising simplicial elements of different dimensions, enabling the integration of continuum models with discrete fracture networks. It handles a spectrum of physical processes, including (un)saturated Darcy flow via mixed-hybrid finite elements, as well as solute transport complemented by sorption, dual porosity, decays, and simple reactions and discretized either by finite volume methods or the discontinuous Galerkin (DG) method. The DG method is similarly employed in the heat transport model. The mechanical model with contacts on the fractures is addressed by finite elements. The software features a preconditioned iterative nonlinear hydro-mechanical coupling with sequentially coupled solute and heat transport processes. On the technical front, Flow123d now features full MPI-aware containerization, supporting both Docker and Singularity, which enhances its portability and deployment efficiency across different computing environments. It is scalable up to hundreds of processors using MPI and PETSc, ensuring high performance in complex simulations. The software also boasts a flexible input interface based on the YAML file format and supports output in GMSH and VTK formats. Recent enhancements include an optimized, vectorized assembly over patches that fit into cache memory, significantly boosting computational efficiency. All these functionalities are implemented in C++ and Python, and the software is distributed under the GNU Public License v3, emphasizing its commitment to open-source development and collaboration.

OpenGeoSys (OGS)<fnref n="5" /> is a scientific open-source initiative for the numerical simulation of thermo-hydro-mechanical/chemical (THMC) processes in porous and fractured media (Kolditz and Bauer, 2004; Kolditz et al., 2012). The basic concept of OGS consist on providing a flexible numerical framework, using primarily the Finite Element Method (FEM) for solving multi-field coupled processes with application in different scientific and technical disciplines (Rink et al., 2012; Bilke et al., 2014; Rink et al., 2014). For example, OGS has been successfully applied in the fields of regional, contaminant and coastal hydrology, technical and geothermal energy systems, geotechnical engineering, energy storage, CO₂ sequestration/storage and nuclear waste management and disposal (Kalbacher et al., 2007; Kalbacher et al., 2012; Kolditz et al., 2021). The actual version OGS-6 is providing complete workflows starting from data integration, HPC for coupled process simulation and using virtual reality (VR) concepts for data analytics. OGS-6 is developed and maintained platform-independently using professional software engineering tools such as version management (GitHub) and containerization (e.g., Docker, Singularity). A strict code review is conducted for quality assurance completed by unit testing and

comprehensive benchmarking. OGS provides open interfaces for combining with other simulators (e.g., GEMS, iPhreeqc for geochemical processes) including Python bindings (Li et al., 2014; He et al., 2015; Kosakowski and Watanabe, 2014). A recent overview of OGS software engineering can be found in Bilke et al. (2019) and Lehmann et al. (2024).

CORE²D V5 (Fernández, 2017) is a code for transient saturated and unsaturated water flow, heat transport and multicomponent reactive solute transport under both local chemical equilibrium and kinetic conditions in heterogeneous and anisotropic media. CORE²D V5 is the 5th version of the code and is available upon request. The flow and transport equations are solved with Galerkin finite elements and an Euler scheme for time discretization. The solute transport equation accounts for advection, molecular diffusion and mechanical dispersion. The chemical formulation is based on the ion association theory and uses an extended version of the Debye-Hückel equation (B-dot) for the activity coefficients of aqueous species. The following chemical reactions are considered: aqueous complexation, acid-base, redox, mineral dissolution/precipitation, cation exchange, surface complexation and gas dissolution/exsolution. CORE²D V5 relies on the “com” thermodynamic database of EQ3/6 (Wolery, 1992). The code also allows the use of other thermodynamic databases. CORE²D V5 is based on the sequential iteration approach to solve for chemical reactive solute transport. Iterations are repeated until prescribed convergence criteria are attained (Samper et al., 2009). The Gaines-Thomas convention is used for cation exchange. Surface complexation is modelled by using three types of protonation/deprotonation sites, SS-OH, SW1-OH and SW2-OH, as proposed by Bradbury and Baeyens (1997). CORE²D V5 takes into account the changes in porosity due to mineral dissolution/precipitation reactions and their feedback effect on the flow and transport parameters under isothermal and nonisothermal conditions (Águila et al., 2020; Samper et al., 2021). CORE²D V5 has been used to model laboratory tests (Yang et al., 2008; Mon, 2017), in situ experiments (Molinero-Huguet et al., 2004; Zhang et al., 2008) and perform long-term simulations for geological high-level radioactive waste repositories (Lu and Zhou, 2011). It has been extensively benchmarked against other codes (Poonoosamy et al., 2021). CORE²D V5 results can be postprocessed with public-domain post-processing tools such as ParaView.

INVERSE-FADES-CORE V2 (Mon, 2017) is a finite element code for modelling non-isothermal multiphase flow, heat transport and multicomponent reactive solute transport under both chemical equilibrium and kinetic conditions in deformable media. The code takes into account the mass balance of water, air, solid and enthalpy; the transport of solids and the mechanical equilibrium. The solute transport equation accounts for advection, molecular diffusion and mechanical dispersion. INVERSE-FADES-CORE V2 solves both forward and inverse multiphase flow and multicomponent reactive transport problems in 1-, 2- and 3-D axisymmetric porous and fractured media. The state variables of the forward model include the liquid and the gas pressures and temperature, which are solved with a Newton-Raphson method. The inverse problem is solved by minimizing a generalized least-squares criterion with a Gauss-Newton-Levenberg-Marquardt method

<footnote n="5">5 www.opengeosys.org</footnote>
<figure><img src="image_5.png" /><figcaption>FIGURE 2
Vision for scale transition schemes: from process to management scales (figure sources: https://developer.nvidia.com, BGR (OGS simulation), Google Earth).</figcaption></figure>

(Dai and Samper, 2004). The forward routines of INVERSE-FADES-CORE have been widely verified with analytical solutions and other reactive transport codes. Mon (2017) implemented reactive gas transport in INVERSE-FADES-CORE V2 by including additional mass balance equations for the reactive gaseous species in the gaseous phase. INVERSE-FADES-CORE V2 has been benchmarked with other codes for benchmarking problems dealing with high-level radioactive waste management (Samper et al., 2024).

## 3 Scale transition schemes for coupled processes

Scale transitions with related, subsequent models are the only way to convey multi-physics, chemical as well as microbiological process information between different scales of interest (Bilke et al., 2019). Due to coupling effects and non-linear phenomena, upscaling methods are a particular challenge for coupled THM/C processes. In general we can distinguish between different directions in multi-scale modeling, e.g., (1) development of inherent upscaling methods (e.g., coarse graining), (2) building connected model chains (see Figure 2), (3) using computational power (i.e., High-Performance-Computing) to increase model domains (Link to previous chapter 2). As there is no perfect way for scale transition, different directions have been followed.

### 3.1 Multiscale modelling

Figure 2 serves as an explanation for (2) building connected model chains. It starts from modeling at the process scale using

methods down to the atomistic, molecular scales such as Lattice-Boltzmann, Lattice-Element approaches, Molecular Dynamics or Particle Methods. These modelling techniques are able to represent fundamental processes at very small (atomistic) scales, but computationally very expensive and difficult to upscale to measurement scales. The opposite end of the scales is represented by so-called complexity-reduced models. At this scale, questions of long-term environmental impacts of radioactive waste repositories have to be addressed. To this purpose data from various sources have to be combined including geological and geographical information as well as socio-economic aspects. For large scale analysis, stochastic approaches and network simulation are being used. Only a small portion from small scale simulations can be used at this scale. However, fundamental process understanding will be a key for more reliable long-term predictions.

A possible link between the small and large scales, in particular for multi-physics problems such as THM/C, is continuum mechanics as it is based on first-order principles (conservation laws of mass, momentum and energy). On the other side continuum mechanics invoke thermodynamic consistent, constitutive laws, which may describe complex material behavior at multiple scales (e.g., micro-mechanics). In addition to fundamental continuum mechanical approaches, numerical methods using significant computational power (HPC) are available for multi-scale purposes. However, it has to be clearly stated that computational mechanics for sure can handle a certain range of process scales but not the whole required scale chain as shown in Figure 2. Therefore, the development of reliable and computational efficient model chains is important to further improve performance and safety assessment tools of radioactive waste management.
<figure><img src="image_6.png" /><figcaption>FIGURE 3
Conceptual approach for fractured porous media (Berre et al., 2019).</figcaption></figure>

In addition to building model chains over scales, as mentioned above, the second challenge for THM/C modeling is the structural complexity of potential host rocks for nuclear waste repositories. A key structural element for geological formations is the presence of discontinuities (e.g., faults, fractures, cracks). Characterizing the state of the stress (i.e., spatial stress directions and magnitudes) in a field is, if not impossible, extremely difficult. Even if the general stress field of a potential site is reasonably estimated, it is nearly impossible to map a precise description of the existing fractures and faults and to predict their potential evolution when temperatures, pore pressures, mechanical stresses are changing. For these reasons, characterization of the earth model is typically combined with stochastic approaches involving multiple realizations in order to bracket the reality. For upscaling, so-called fractured-porous media, different methods have been developed in past (Berre et al., 2019, Figure 3): (1) Discrete-Fracture-Network (DFN) models explicitly take into account the presence of discontinuities, (2) Single and Multi-Continua approaches use effective, statistically representative parameters for describing fracture behavior (REV concept). The main idea of Multi-Continua models is to couple multiple continua via exchange functions, where every continuum is representing a typical part, i.e., fracture systems, rock matrix. (3) Discrete-Fracture-Matrix (DFM) models are combining both discrete and continual behaviors. Discrete approaches (1&3) require very precise numerical methods for accurate representation of exchange processes between fracture and matrix compartments.

## 3.2 Multi-scale experiments

Dissolution and precipitation of minerals and its impact on the transport of fluids and solutes in porous media is a key factor for the

long-term behaviour of natural and engineered systems in the subsurface. The implementation of such coupled processes into numerical reactive transport codes requires a mechanistic process understanding and model validation with quantitative experiments (Tartakovsky et al., 2008; Katz et al., 2011; Poonoosamy et al., 2015; Poonoosamy et al., 2019). These experiments are usually based on a simple chemical setup with well-defined kinetic and initial transport parameters. The discrepancies between experimental observations highlighted the limitations of continuum scale models when it comes to predicting changes in the transport properties in porous media due to mineral precipitation and consequent porosity changes (Chagneau et al., 2015; Poonoosamy et al., 2015; Poonoosamy et al., 2019). The primary cause for these discrepancies is that reactive transport operates at the continuum scale and considers upscaled parameters. Consequently, the interdependencies between porosity and transport parameters for instance, are accounted for in numerical models by using empirical relationships developed on the continuum scale (Hommel et al., 2018). These relationships were either fitted to experimental data or based on theoretical considerations (Poonoosamy et al., 2019; Seigneur et al., 2019). However, solute transport depends on more complicated pore space changes (pore geometry, pore connectivity) than the upscaled parameter “porosity” as accounted for in continuum-scale reactive transport approaches. The averaging process also includes geochemical parameters such as reactive surface area and reaction rates. The heterogeneity in porous media likely controls the mass fluxes and residence time in the reactive media and therefore are known to control the effective rates of mineral dissolution by reducing effective mineral surface areas physically (e.g., in fractures, or at partial saturation) and chemically (Berkowitz et al., 2016; Deng et al., 2018b; Deng et al., 2018a; Steefel et al., 2015a; Jung and Navarre-Sitchler, 2018b; Jung and Navarre-
Sitchler, 2018a). Moreover, the reactive surface area of minerals constituting the porous media is not directly experimentally accessible and is often used as a fitting parameter in reactive transport simulations (Noiriel et al., 2009). The development of processes-based relationships rationalizing the effective transport properties of porous media and the reactivity of minerals remains one of the challenges of reactive transport models (Prasianakis et al., 2017; Prasianakis et al., 2018; Noiriel and Soulaine, 2021). The combination of high-quality imaging of evolving porous media and pore scale modelling is one of the key tools to develop these process-based relationships. For instance, Steinwinder and Beckingham used pore scale modeling on X-ray tomography images of an Alberta Basin Sandstone to study the impact of pore and pore throat size distribution variations on porosity and permeability evolution during coupled mineral dissolution and precipitation reactions, deriving a constitutive equation linking porosity and permeability. Various imaging methodologies have also been developed for the visualization of reactive flows in porous media, e.g., the positron emission tomography (PET, Kulenkampff et al., 2018), magnetic resonance imaging (Fridjonsson et al., 2011; Poonoosamy et al., 2020), and neutron tomography (Shafizadeh et al., 2015). However, X-ray tomography is the most widely applied tool (Spanne et al., 1994; Coles et al., 1996; Coles et al., 1998; Lindquist et al., 2000; Noiriel et al., 2004; Blunt et al., 2013; Menke et al., 2021). The principal reason is that besides conferring advantages such as non-invasive, non-destructive, and dynamic assessment of flow paths, it can also assess mineralogical changes. Recently, microfluidics or Lab-on-chip experiments (Harrison et al., 2017; Soulaine et al., 2017; Soulaine et al., 2022; Poonoosamy et al., 2019) have gained a lot of interest to unravel coupled hydrogeochemical processes at the pore scale. Lab on a chip can mimic the pore structure of rocks and enable the monitoring of mineral nucleation, precipitation, and dissolution using optical microscopy and spectroscopic techniques e.g., Raman spectroscopy (Poonoosamy et al., 2020) or synchrotron-based tomography (Deng et al., 2020). Lab-on a chip paired with pore-scale modelling has proven to be a powerful tool for rationalizing coupled processes e.g., evaluate the impact of hydraulic heterogeneity on precipitation processes (Poonoosamy et al., 2020), decode oscillatory zoning exhibited by solid solutions crystallizing in porous media under diffusive transport regime (Poonoosamy et al., 2021), or parameterize porosity-diffusivity relationships with respect to coupled mineral dissolution-precipitation reactions (Poonoosamy et al., 2022; Lönartz et al., 2023). Such experiments can be used to calibrate pore scale models prior to simulations over a wide range of Peclet and Damköhler numbers. The numerical datasets generated by these pore scale models can be further used to train convolutional neural network models, which can then be further applied to derive process-based relationships (Prasianakis et al., 2020).

## 3.3 Scale transitions schemes for hydro-mechanical processes

In the following section we briefly describe numerical methods which are actively developed by the project partners in order to tackle scale transition problems.

### 3.3.1 Two common approaches for numerical simulation of fractures generalized/extended finite-element (G/XFEM) and Variational Phase Field (VPF) methods

#### 3.3.1.1 Generalized/extended finite-element-method

In recent years, generalized (Strouboulis et al., 2000) or extended (Belytschko and Black, 1999) finite element methods (G/XFEM) and phase-field models (Bourdin et al., 2000) for the description of existing and developing discontinuities and singularities within continuum mechanical approaches have established (Belytschko et al., 2009) themselves ahead of all others in recent years. Both methods differ fundamentally and have their own strengths and weaknesses. G/XFEM locally extends the approach and test function space by formulations that can map the discontinuous course of the solution and introduces corresponding additional local degrees of freedom. Usually, this approach is combined with so-called level set methods, which help to localize the discontinuity and thus ultimately determine in which elements the solution space has to be extended. This approach allows the approximation of discontinuous solutions on comparatively coarse grids, but requires extensive programming infrastructures for the treatment of flexible additional degrees of freedom, level sets and other aspects, which require a not inconsiderable implementation effort, especially in branched crack systems. In contrast, the phase field method represents a sharp interface such as discrete crack with a smoothly varying scalar field which is regularized by the internal length. Thus the discontinuity is not described as a low-dimensional unit, but by the means of an additional variable defined on the entire solution area, which originates in certain characteristics of an energetically motivated variational formulation in Francfort and Marigo (1998). This approach can be implemented into existing continuum codes without additional algorithmic effort where it is immediately available for solving arbitrary-dimensional problems. Crack formation, propagation, branching and unification also occur without additional effort.

G/XFEM (Belytschko and Black, 1999) was originally developed for crack propagation problems and was also applied in the geotechnical context, e.g., for multiphase flows (Chessa and Belytschko, 2003; Mohammadnejad and Khoei, 2013) and heat transport (Khoei et al., 2012; Shao et al., 2019). Current developments of generalized and extended finite element methods in the context of hydraulic stimulation are mainly concerned with the efficient coupling of solid-state and flow-mechanical problems (Yazid et al., 2009; Watanabe et al., 2012; Meschke and Leonhart, 2015).

#### 3.3.1.2 Variational phase-field models for fractures

The theoretical basis of the so-called phase-field model for fracture was incepted in the work by Francfort and Marigo (1998) that revisited Griffith's fracture criteria as an energy minimization problem (variational approach), and later on its numerical solution was proposed by regularizing the energy functional with a smoothly varying scalar (phase-field variable) and a regularization length parameter (internal length) by Bourdin et al. (2000). Although their work had gone relatively unnoticed in the mechanics community for a decade or so, its popularity has been surging since a ceremonial publication by Miehe et al. (2010) who coined the method "phase-field model". Its brief
<figure><img src="image_7.png" /><figcaption>FIGURE 4
Concept of the Variational Phase Field (VPF method) for handling discontinuities in the context of coupled processes.</figcaption></figure>

narrative history can be found in a summary article in Bourdin et al. (2000). In recent years, the phase-field model enjoys its wide applicability ranging from brittle and cohesive fracture (Bourdin et al., 2008; Amor et al., 2009; Hakim and Karma, 2009; Kuhn and Müller, 2010; Pham et al., 2011; Verhoosel and de Borst, 2013; Vignollet et al., 2014; Ambati et al., 2015a; Marigo et al., 2016; Wu, 2017; Sargado et al., 2018; Tanné et al., 2018). To ductile fracturing (Ambati et al., 2015b; Miehe et al., 2015; Kuhn et al., 2016; Alessi et al., 2017) fatigue (Alessi et al., 2018), desiccation fracture (Maurini et al., 2013; Cajuhi et al., 2018), dynamic fracturing (Bourdin et al., 2008; Borden et al., 2012; Schlüter et al., 2014; Li et al., 2006), and hydraulic fracturing (Bourdin et al., 2012; Miehe et al., 2015; Yoshioka and Bourdin, 2016; Ehlers and Luo, 2017; Lee et al., 2017; Santillán et al., 2017; Chukwudozie et al., 2019).

A set of crack (Γ) in a domain (Ω) are discontinuities where a jump in quantities such as displacement is observed (Figure 4). Therefore, when represented in a discretized form (e.g., finite element), the discretization must honor the discontinuous geometry of Γ, and it can get easily untractable when the crack set evolves with time (Figure 4, left). In phase-field models for fracture, the discontinuity (Γ) is represented by a smoothly varying scalar (v as in Figure 4, right) so that the discretization does not need to track the geometry change but the change in the scalar field instead. Naturally a special care must be taken in setting up such a scalar field in order to preserve the energy equivalence in the two systems and the convergence of the phase-field representation to the discrete fracture as the regularization approaches zero (Chambolle, 2004; Bourdin et al., 2008).

First extension of the phase-field model to hydraulic fracturing was proposed in Bourdin et al. (2012) by adding the fluid driven pressure work on the crack lips and Wheeler et al. (2014) further adapted it by adding the poro-elastic effects. In these models, while the pressures in the non-fractured medium (Ω/Γ) and the fracture (Γ) are distinguished, their profiles are considered spatially uniform within the domains. Phase-field fracture models coupled with Darcy-Reynolds flow system were developed in Mikelić et al. (2015a), Miehe et al. (2015), Miehe et al. (2016), Santillán et al. (2017), Xia et al. (2017), Lee et al. (2017) and Darcy-Stokes type in Ehlers and Luo (2017), Heider and Markert (2017). Other “homogenization” approaches include a model that uses the phase-field variable as an indicator in Mikelić et al. (2015b) or as a weight function for permeability averaging in Miehe et al. (2016). Heider and Markert (2017) and Ehlers and Luo (2017) developed the

phase-field fracture model within the framework of the Theory of Porous Media (TPM) and considered the phase-field variable as a weight between the Darcy flow in porous media and the Stokes flow in fracture. Alternatively, coupling with fluid flow has been simplified by linking to an external standalone fluid flow simulator in Wick et al. (2016), Yoshioka and Bourdin (2016) where a phase-field variable was used to differentiate the fracture from the non-fracture domain. Instead of using the phase-field variable as an indicator or a weight, Chukwudozie et al. (2019) applied the phase-field calculus to homogenize the system while conserving the mass.

Another mass conserving approach was proposed by Santillán et al. (2017) where an independent 1D domain is generated dynamically out of the main computational domain based on certain the phase-field threshold and the fluid flows are solved separately and iteratively.

Figure 5 depicts a graphical representation of various numerical methods for the description of evolution of discontinuities, i.e., (a) lower-interface-method (LIE), (b) variational-phase field method (VPF), and non-local-deformation (NLD) methods, which will be also used in the context of various EURAD work packages (e.g., WP-GAS).

### 3.3.2 Multi-continua methods for hydro-mechanical processes

#### 3.3.2.1 Modelling approach

So-called multi-continuum models are often used to better determine the heterogeneity of fractured or karstied aquifers. For example, double-continuum models take advantage of the different hydraulic properties of the rock matrix and the fissure system. While fissure systems are characterized by high permeabilities and low storage capacities, the rock matrix has low permeabilities and high storage capacities. Thus, the fissure system and the rock matrix are attributed to separate “overlapping” continua, which are coupled with each other via source-sink terms in the entire model area (see Figure 3).

Similar to the multiphase concept of a porous medium, several state variables are defined at each point in the multi-continuum model. In addition, the individual continua also have different material parameters. If one wanted to use the category of phase for a double-continuum model, it would be an allocation from a macroscopic point of view: thus not in the microscopic sense fluid and solid, but fissure system and matrix. The key point in multi-
<figure><img src="image_8.png" /><figcaption>FIGURE 5
Graphical comparison of numerical methods for description of discontinuities (Yoshioka and Bourdin, 2016).</figcaption></figure>

continuum modelling is the description of the interaction between the individual continua. In a sense, the exchange relationships represent non-equilibrium conditions. The recording of the local dynamics of exchange processes (gradient approaches) is problematic. Compared to the classical single-continuum approaches, multi-continuum models are better able to map heterogeneities. The differentiated continua have to process smaller structural variabilities, which has an advantageous effect on the averaging technique (e.g., size of the REV). Average state variables can be determined separately for the fissure and matrix continuum rather than an average state variable for the entire fissure-matrix system.

The numerical implementation of multi-continuum models can be done on the one hand by coupling separate equations for fissure and matrix media via exchange terms (Bibby, 1981; Huyakorn et al., 1983; Birkhölzer et al., 1993a; Birkhölzer et al., 1993b) or on the other hand by simultaneously solving the equation systems using iterative techniques (Pruess and Narasimhan, 1982; Neretnieks and Rasmuson, 1984).

Multi-continuum models are widely used, for example, in the interpretation of pump tests in solid rock or the delayed hydraulic response to precipitation events in karst aquifers (Sauter et al., 2006), the hazard assessment of groundwater deposits in near-surface solid rock areas and the exploration of oil and gas deposits (Gilman et al., 1983). An advantage for the model calibration is the justifiable exploration effort to adapt the limited number of model parameters. Typical geological formations for which classical single-continuum models may fail and require multi-continuum approaches include porous sandstones, karst aquifers, loose rocks with clay

intercalations or silty portions, and fissured bedrock. The application potential of the multi-continuum concept to process simulation in hierarchical fissure or pore systems (e.g., intergranular pores of the grain skeleton and micropores in the rock grains) seems to be interesting. In the process, fracture systems of different scale sizes are to be assigned individual continua (Clemens et al., 1996). More recent efforts in model development also concern the coupling of discrete fissure models with multi-continuum models.

### 3.3.2.2 Hydromechanical problems

Dual and multiple continua for hydromechanical problems follows from original models for this kind of conceptual representation of individual processes of solute transport with dead-end pores or mechanics of composites. More details of development of the porous/fracture flow and solute transport models is given in the (section 3.3.2.3). Interaction of stress and deformation with fluid flow, the general used theory is by Biot (1941) and the generalisation into the dual porosity concept of poroelasticity is by Aifantis (1977).

The recent work aims to various applications, with shares of reservoir engineering, hydraulic fracturing, and deep geological repository. More frequent are hierarchical porous media than upscaled fractured rock. Ashworth and Doster (2019) compares various formulas of previous authors, especially 3 works in details, including a table form of coefficient formulas list, and evaluates the importance of coupling. Example calculation is with Mandel problem (Nguyen and Abousleiman, 2010, Figure 6).

Liu et al. (2018) derive theoretical formulation and analytical solutions to the inclined wellbore problem and axisymmetric
<figure><img src="image_9.png" /><figcaption>FIGURE 6
Problem used for benchmarking dual-continuum poroelastic
models (Ashworth and Doster, 2019).</figcaption></figure>

Mandel-Type problem of dual-porosity, dual-permeability poro-
chemo-electro-elasticity. The model is able to explain phenomena
which could seem anomalous without consideration of the coupling
effects. Kim et al. (2018) simulates hydraulic fracturing process by
means of 3D thermo-poro-mechanical model including the dual
porosity approach. The movement of water is retarded with respect
to the fracture propagation. It also theoretically demonstrates which
of the effects can be detected by electromagnetic geophysical
methods. Vasilyeva et al. (2019) fits the method to perforated
domains, which establish quite active sub-topic within the dual-
continuum upscaling. Non-local multi-continuum method (Chung
et al., 2018) is used which is specific in construction of local basis
functions capturing also the non-local effects. Special attention is to
generalisation of the boundary conditions on perforations. Choo
et al. (2016) and Sanchez et al. (2016a) deal with unsaturated
conditions of the dual-porosity media and model the hydro-
mechanical and thermos-hydro-mechanical processes,
respectively. Kim et al. (2015) suggest an algorithm for coupling
of chemo-thermo-poro-mechanics, for fractured rock. Wang et al.
(2014) simulates the fractured rock, with effect of stress on flow and
transport, as an issue of the deep geological repository safety. The
transport is represented by multiple continua (MINC) and the
hydro-mechanics by Oda (1986) crack tensor theory, denoting
the method as Ex-MINC (extended) - Figure 7. The benchmark
problem was also used to successfully compare this upscaled
homogenization concept with several solution on full fracture
network (small-scale view).

3.3.2.3 Lumped Parameter Modeling

The term Lumped Parameter Modeling (LPM) is not commonly
used in context of repository models, but many of them are built on

these principles. The meaning is that a system is represented by a single
value of certain quantity (either a parameter or a state variable) for
whole or a component, in an average or integral meaning, contrary to
point values representing infinitesimally small volume. The classical use
of LPMs is to create the model based on system observation, capturing
its behavior as a link between “input” and “output”, or control and
response quantities—understanding the system as “black-box”.
Alternatively, LPMs can be seen as an upscaling method, using
simplified and numerically efficient model with less degrees of
freedom whereas the lumped behavior is based on solution of
detailed equations controlling the system (physically-based model).
The principal ingredient of such approach is then a procedure of
integral postprocessing of quantity fields resulting from the
controlling equations.

Besides LPM, also a term “compartment model” is used in
similar meaning, but including flexibility such that the meaning of a
compartment can be either a lumped subsystem, or a discretization
element (Figure 8). Such case is the concept in Goldsim software
(Golder Associates, 2014), which is also considered as one of the
options in Czech performance assessment of spent nuclear fuel
repository (Vopálka et al., 2006).

There are not many recent studies using lumped models for flow
and solute transport, which are meant also as basis for the
performance assessment models. Therefore the references cover a
longer period and contain work which was not further developed in
the recent 5–10 years. The state of the art report is limited to the
partial topics of how the continuous field quantities of 2D/3D
models are integrated into the lumped parameter models and
how the models with various features related to LPMs are used
in performance assessment.

The principle of “integrating” the physics into the compartment
block in/out behavior in Goldsim software (Golder Associates, 2014)
is the use of analytical solution, which calculates a “breakthrough” of
radionuclide controlled by homogeneous fracture and rock
parameters for the respective compartment. Similar concept of
analytical solution is introduced by Chopra et al. (2015) which
improves the referred existing (mostly widely used) solutions with
adding the proper hydrodynamic dispersion solution.

The older work of Gylling et al. (1997) represents a reference
case of performance assessment model which has features of both
distributed and lumped model. The channel network model
(CHAN3D) for solute transport in fractured granite uses direct
input of data (not as upscaling), but the channels are defined with
strong abstraction of physical meaning (regularized geometry with
transforming of parameters to the same effective properties).

A lumped model denoted as Tank model for regional
groundwater flows is proposed for a case study of the
groundwater flow system (Kazumba et al., 2008). The aquifer is
divided into Tanks within which the average values of the
groundwater levels are assumed to be representative in the
Tanks. For each Tank the mass-balance equations are written.
The Quasi-Newton optimization technique together with
Akaike's Information Criterion, AIC, are employed in order to
balance the model/observation fit and reliability of the parameter
estimation. The results suggest disconnection of some regional
aquifers as well as inflow of seawater to the deeper layers.

The model representing the groundwater pathway in risk
assessment or performance assessment developed by Robinson
<figure><img src="image_10.png" /><figcaption>FIGURE 7
Scheme of the extended multiple interacting continua (Ex-MINC) method (Wang et al., 2014).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 8
Scheme of a compartment capturing a solute transport process in the lumped-parameter model sense (redrawn after Golder Associates, 2014).</figcaption></figure>

and Chu (2013) is called the Residence Time Distribution (RTD) Mixing Model (RTDMM), allows for an arbitrary distribution of fluid travel times to be represented, to capture the effects on the breakthrough curve of flow processes such as channelized flow and fast pathways and complex three-dimensional dispersion. Mathematical methods are derived directly from the theory of residence time distributions. A simple mixing model is presented, along with the basic equations required to enable an arbitrary RTD to be reproduced using the model. The advantages include easy incorporation into a multi-realization probabilistic simulation; computational burden similar to 1D model with the same number of grid cells; and straightforward implementation into existing flow and transport modelling codes. Examples include

cases of a radioactive decay chain, dual porosity transport and sorption.

The case study (Hagedorn et al., 2018) examines the nitrate pollution vulnerability of groundwater in sedimentary aquifers of California's South Coast Range using stepwise logistic regression modelling. The study highlights the utility of two parameters, Dissolved Oxygen and Modern Water Proportion, for predicting nitrate contamination vulnerability. Authors conclude a possibility to use the results to guide field investigations in areas with little data but similar hydrogeology and land use.

Overview paper (Rechard et al., 2014) summarizes modelling of radionuclide transport in the unsaturated and saturated zone conducted between 1984 and 2008 to evaluate the performance
of Yucca Mountain disposal system. The listed models start at 1-D transport for a single porosity media without lateral dispersion in both the saturated and unsaturated zone for the first assessment in 1984. Gradually the models progressed until 2008 to the 3-D transport model the saturated zone, combined with 1-D transport results, which evaluated decay of radionuclides, in order to evaluate compliance with groundwater protection requirements. Uncertainty in flow within the unsaturated and saturated zone was generally important to explaining the spread in the individual dose performance measure.

For Nuclear fuel cycle and nuclear waste disposal decisions, Huff (2017) developed the Cyder disposal environment and repository module, which implements models of medium likelihood of hydrological radionuclide to support assessment suitable for fuel cycle simulation simulator “Cyclus”. By means of four modules with close-form solution of physics it is possible to reduce the computational need.

The lumped Rapid Saline Groundwater Exfiltration Model (RSGEM) was developed to support operational water management of freshwater resources in coastal lowlands, a need exists for a rapid, well-identifiable model to simulate salinity dynamics of exfiltrating groundwater (Delsman et al., 2017). The RSGEM was applied to a field site in the coastal region of the Netherlands, parameter estimation and uncertainty analysis were performed using generalized likelihood uncertainty estimation. The model showed good correspondence to measured groundwater levels, exfiltration rates and salinity response. Moreover, RSGEM results were very similar to a detailed, complex groundwater flow and transport model previously applied to this field site.

Kawasaki and Ahn (2008) also employs the residence time distribution and includes stochastic methods like random walk. The evaluation can be modularised and so it convenient e.g., for efficient representation of multiple waste packages. It also avoids artificial numerical dispersion for coarse discretisation.

Murakami and Ahn (2011) develops a model with array of compartments, among which a transition probability matrix describes radionuclide transport (Markov-chain process). Transition probabilities are upscaled from finite-element flow and transport solutions. Authors suggest efficiency of the method for modelling effects of parameter changes. Paper of Cadini et al. (2010) brings another example of compartment model with stochastic features, where detailed local-scale modelling results “feed” a global-scale analysis of the repository, at reasonable computational expenses.

The topic was recently revisited in the EURAD project (work package DONUT), where a postprocessing method was developed, which derives the lumped-parameter model inputs from a 3D flow and transport model containing detailed hydrogeological data of a site (Landa and Hokr, 2022). The method was demonstrated on nine repository candidate sites in the Czech Republic, calculating the transport path characteristics for Goldsim software. The acquired data could be also used directly to compare the sites or to narrow their selection.

## 3.4 Scale transitions schemes for hydro-chemical processes

### 3.4.1 Lattice Boltzmann method (LBM)

The lattice Boltzmann method (LBM) is a special discretization of the Boltzmann equation operating at the level of velocity

probability distribution functions $f_i$, with focus on fluid dynamics and mass and heat transport phenomena. At every distribution function corresponds a discrete velocity vector (Frisch et al., 1986; Succi et al., 1989). Different discretization schemes lead to different number of discrete velocities, there exist different lattice models (Higuera et al., 1989; Qian et al., 1992). For the recovery of Navier-Stokes fluid motion equations, for two-dimensional simulations the standard D2Q9 square lattice with 9 discrete velocities is usually implemented, and for three dimensional simulations the D3Q27 cubic lattice with 27 discrete velocities, due to their simplicity and robustness in complex geometry domains (see Figure 9). The diffusion equation can also be recovered on the same or smaller velocity sets.

For reactive transport simulations in geological relevant flows, the mass is transported via diffusion and advection. For the modelling of the advection-diffusion equation the passive scalar approach is usually implemented. Such a model is composed of a basis fluid medium (a population set which has the dynamics as described by the Navier-Stokes equations, plus a passive scalar coupled population set that simulates the diffusion of ions. The basic lattice Boltzmann discretization approximates the Navier-Stokes equations to a certain accuracy, which for very high flow velocities can become significant and can manifest itself as the loss of Galilean invariance. There are several ways to overcome this issue. For example, the isothermal guided equilibrium nine-velocity model (D2Q9 lattice), can be selected as the basis model for the advective flow field as well as for the field that advect and diffuse (Prasianakis et al., 2009; 2017). This model provides enhanced Galilean invariance and rotational isotropy on the standard lattices.

The discrete velocities of populations $f_i$ for $i = 0–8$, are $c_i = (0, 0)$ for $i = 0, c_i = (\pm 1, 0)$ and $(0, \pm 1)$ for $i = 1–4$, and $c_i = (\pm 1, \pm 1)$ for $i = 5–8$ (Qian et al., 1992). The following population-moments correspond to the density of the solution $\rho$ and the momentum $j_a$ in the direction $a = x, y$ as shown in Equation 1:

$$ \sum_{i=0}^{8} f_i = \rho, \sum_{i=0}^{8} c_{ia} f_i = j_a \qquad (1) $$

The guided equilibrium populations $f_i^{eq}$ are given in a closed form, where $T_0 = 1/3$ as shown in Equation 2:

$$ f_i^{eq} = \rho \prod_{a=x,y} \frac{(2c_{ia}^2 - 1)}{2c_{ia}^2} (c_{ia}^2 - 1 + c_{ia}u_a + u_a^2 + T_0) \qquad (2) $$

The Boltzmann BGK equation is solved as shown in Equation 3:

$$ \partial_t f_i + c_{ia} \partial_a f_i = -\frac{1}{\tau} (f_i - f_i^{eq}) \qquad (3) $$

where $\tau$ is the relaxation parameter and $\mu = \tau\rho T_0$ is the resulting macroscopic dynamic viscosity. BGK stands for the Bhatnagar-Gross-Krook collision model as depicted in the right hand side of the aforementioned equation and describes the relaxation of populations $f_i$ to their equilibrium state $f_i^{eq}$ with relaxation time $\tau$. For the advection-diffusion reaction equations of the reactive species a D2Q9 model can also be implemented. The solutes (primary and secondary species) will advect and diffuse along the streamlines of the flow field. The equilibrium populations $\xi_i^{eq}$ for the reactive species $\xi_i$ are given as shown in Equation 4:
<figure><img src="image_12.png" /><figcaption>FIGURE 9 Standard lattices (left) two dimensional 9-velocity lattice (D2Q9), (right) the three dimensional 27-velocity lattice (D3Q27).</figcaption></figure>

<figure><img src="image_13.png" /></figure>

$$ \xi_i^{eq} = C_i \prod_{\alpha=x,y} \frac{(1 - 2c_{i\alpha}^2)}{2c_{i\alpha}^2} (c_{i\alpha}^2 - 1 + c_{i\alpha}u_{\alpha} + T_0), \quad (4) $$

where $C_i$ is the concentration of the considered ions and $u_\alpha$ is obtained from the basis model. The relevant population moment that corresponds to the concentration is as shown in Equation 5:

$$ \sum_{i=0}^{8} \xi_{i} = C_{i} \qquad (5) $$

Advantages of the LBM for reactive transport geochemical applications are the ability to simulate evolving geometries of high complexity (due to dissolution or precipitation), the local mass and momentum conservation, as well as the parallelization efficiency. The high scalability of LBM codes is related to the locality of the mesh based particle-particle interactions, which requires only next neighbor communication between the grid points. To that, it should be added that LB codes are usually simple in terms of coding and they do not depend on external libraries. This allows the incorporation of additional physics (electrochemical, multiphase etc.) in a relatively fast pace, as soon as the theoretical model is expressed in the lattice Boltzmann framework.

### 3.4.1.1 Classical nucleation theory embedded in a cross scale LB model for precipitation

Mineral dissolution and precipitation is a multiscale process that covers several length and time scales. According to the classical nucleation theory, at the molecular level, when a supersaturation level is exceeded the ions start forming critical nuclei which serve as the bases of crystal formation. The energetic barrier that has to be overcome is different depending on the mechanism of precipitation. Heterogeneous precipitation (HET) proceeds on an existing template, and the formation of critical occurs at lower supersaturations, compared to the homogeneous nucleation (HON). Homogenous nucleation takes place in the bulk solution and requires no substrate, but higher level of supersaturation. In Figure 10 (left) a schematic representation of this mechanism is plotted. The blue curve shows the HET nucleation path and the red curve the HON. In Prasianakis et al. (2017), the classical nucleation

theory was implemented within the LBM framework. Once the pore structure is known (e.g., x-ray tomogram), the reactive surface areas for heterogeneous nucleation are resolved at the micrometer level. For the homogeneous nucleation, a sub-lattice model has been implemented. The formation of critical nuclei depends on saturation index. The number of molecules that form the critical nuclei and induction time are modelled by CNT, and the precipitation proceeds on spherical particles which grow until they fill the space. The surface area for precipitation is the sum of the surface of the spheres. A schematic of this process is shown in Figure 10 (right). At the moment only 2D simulations have been conducted, and only for barite and celestite minerals whose precipitation follows the classical nucleation paths. While challenging, this model can be potentially applied for cases where minerals follow non-classical nucleation paths. The parametrization of HET and HON nucleation mechanisms takes into account also the contact angles of different phases and the surface tensions.

### 3.4.1.2 Electrochemical couplings at lattice Boltzmann framework

The dominating pores for ion transport in tight porous media (such as high-performance concrete and compacted clays) are normally at several nanometers size. At this scale, the interface effect, especially the electrical double layer effect, has an important role on ion transport. For example, the relative effective diffusivities of different species in cement paste follow the order $CL^{-} \approx HTO/O_{2} > Na^{+} > Cs^{+}$ (Yang et al., 2019a). However, this is completely different from the order: $CL^{-} < HTO/O_{2} < Na^{+}$ obtained by the classical Nernst-Planck equation as shown in Equations 6, 7:

$$ \frac{\partial C_i}{\partial t} + \nabla \cdot J_i = 0 \qquad (6) $$

$$ J_i = -D_{i,0} \nabla C_i - D_{i,0} \frac{z_i e C_i}{kT} \nabla \psi \qquad (7) $$

where $C_i$, $D_{i,0}$, $J_i$, $\psi$ and $z_i$ denote the aqueous concentration (mol/m³), the diffusion coefficient in free water (m²/s), flux (mol/m²/s), local electrical potential (V) and valence of the i'th species, respectively; t is time (s), e the absolute charge of the electron
<figure><img src="image_14.png" /><figcaption>FIGURE 10</figcaption></figure>

Cross scale precipitation modelling. (left) Classical nucleation Theory and energetic barriers to produce the critical nuclei for homogeneous (HON) and heterogeneous (HET) nucleation mechanisms. The two dimensional 9-velocity lattice (D2Q9), (right) Sub-lattice model for homogeneous nucleation which allows to model the evolution of reactive surface areas.

(C), k the Boltzmann constant (J/K) and T the absolute temperature (K). This classical model cannot accurately capture the ionic steric effect and the ion-ion correlation effect (Gillespie, 2014). Therefore in order to consider the steric effect and the ion-ion correlation effect, an excess chemical potential $\mu_i^{ex}$ term is usually added to the classical chemical potential as shown in Equation 8:

$$J_i = -D_{i,0} \nabla C_i - D_{i,0} \frac{z_i e C_i}{kT} \nabla \psi - D_{i,0} \frac{C_i}{kT} \nabla \mu_i^{ex} \quad (8)$$

This equation is the modified Nernst-Planck equation. However, the estimation of excess chemical potential in the Stern layer is still developing. The model was set up by Bikerman (1942) to correct the ion concentrations in the Stern layer. After Bikerman's model, numerous continuous-scale theoretical models were developed to calculate the excess chemical potential in the Stern layer, but still each theory has its limitations (Lu and Zhou, 2011; Gillespie, 2014; Qiao et al., 2016). The exact ion distribution in the EDL at an arbitrarily charged surface or ionic strength of the electrolyte solution can be obtained by molecular-scale simulation. Therefore, the statistical information obtained by the molecular-scale simulation can be used to calculate the dimensionless excess chemical potential $\mu_i^{ex}/kT$ as shown in Equation 9:

$$\frac{1}{kT} \mu_i^{ex}(x) = -\frac{ez_i}{kT} \psi^{GCMC}(x) - \ln\left(\frac{C_i^{GCMC}(x)}{C_{i,\infty}}\right) \quad (9)$$

where the equilibrium ion concentration distribution $C_i^{GCMC}(x)$ and the electrical potential distribution $\psi^{GCMC}(x)$ away from the solid surface are obtained by the grand canonical Monte Carlo scheme at atomic scale. This dimensionless excess chemical potential $\mu_i^{ex}/kT$ from the grand canonical Monte Carlo simulation is then applied to fitting with the nonlinear functions $F(x, C_{i,\infty})$ to get the relationship of $\mu_i^{ex}/kT$ with respect to the distance $x$ and bulk concentration $C_{i,\infty}$. This upscaling strategy combines pore-scale modeling with atomic-scale modeling by adding the correction term for the excess chemical potential into the classical Nernst-Planck equation to account for steric and ion-ion correlation effects. After this, the modified Nernst-Planck equation is solved by the Lattice-Boltzmann method to simulate ion transport in tight porous media at pore scale and to calculate the effective diffusion coefficient of each species.

### 3.4.1.3 Partial saturation

In unsaturated nano porous media, electrical double layer (EDL) effects caused by the charged surface of clays can significantly enhance cation diffusion but reduce anion diffusion (Wu, 2017; Yang et al., 2019b). On the other hand, since the polarizability of anions is normally larger than that of cations, anions are closer to the air/water interface than cations in partially saturated porous media (Olivieri et al., 2018; Jungwirth and Tobias, 2022). As the thickness of the water layer reaches the nanoscale, the air-water and water-solid interfaces may both have a significant influence on ion distribution and further affect ion transport in partially saturated porous media. A novel pore-scale model to solve coupled Poisson-Nernst-Planck equations has been established to simulate solute diffusion in unsaturated compacted clays. Results indicate the critical influence of the electrical double layer effect on ionic radionuclide transport in nanopores of clay-based materials (Yang et al., 2024). By passing we note that the mobility of gases in partially saturated clays has been also evaluated using molecular dynamics modelling in Owusu et al. (2022), Owusu et al. (2023), and the results can be upscaled directly to pore-level methods.

### 3.4.2 Multi-continua methods for hydro-chemical processes

Similar conceptual models of solute transport are covered by slightly variable terminology-double/dual porosity, double/dual continuum, two-region model. Contrary, the terms describing the configuration can be the same used as related to various physical processes. An exception is often such that dual-permeability and dual-porosity distinguish a relation to either the fluid flow or the solute transport.

Examples of early work are Coats and Smith (1964) with explanation of stagnant water in dead-end pores and, from another application field, Bedford and Stern (1972) developing theory for mechanics of composite material. The common idea for conceptual models is that a heterogeneous material is represented macroscopically by two or more interacting continua, each with separate degrees of freedom for a state physical quantity in a geometric point of space, understood in sense of representative elementary volume. Various terminologies are for denoting the meaning of continua, for solute transport it is typically the mobile and immobile zone.
<figure><img src="image_15.png" /><figcaption>FIGURE 11
Concept of the multiple interacting continua (MINC) method (Jiang and Younis, 2015) redesigned form of original (Pruess and Narasimhan, 1982).</figcaption></figure>

In the context of upscaling, the multi-continuum models are used as homogenized counterpart of fractured rock, capturing the different behaviour of fracture volume and matrix pores, as introduced e.g., by Barenblatt and Zheltov (1960) and Warren and Root (1963).

One of the further direction of concept development was the “multiple interacting continua” (MINC) by Pruess and Narasimhan (1982). It assumes a nested sequence of continua which represent a virtual distance from the fracture into the matrix, i.e., one degree of freedom for the immobile zone is multiplied to represent the micro-scale heterogeneity more accurately. The models are in quite continuous interest in time and it is a live topic both for application of methods for specific case studies and for further developments and extensions of theory or calculation efficiency. The number of relevant publications in the last 5 year creates a significant share of the total number over the history (up to 30%).

The recent studies are oriented to various types of improvements. A lot of studies is related to the reservoir engineering especially gas/oil in shales, CO₂ enhancing or CO₂ storage, with effort to suit the methods and data for specific type of rocks.

Examples of generic topics (Figure 11):

* Procedures to obtain dual- or multi-continuum parameters from local

* Relevant pore-scale properties which are (controlling) for determination of macro-scale parameters

* Selection of number of continuum levels and their association to the material

* Extension of methods from standard flow and transport equations to mode complex processes (non-linear, coupled problems)

* Numerical efficiency improvement

### 3.4.2.1 Phenomena capturing, conceptual improvement

Vogler et al. (2018) compares a hierarchy of upscaled dual porosity models and by comparison with tracer experiment they show that the dominant error arises from neglecting of local convection and that the properly upscaled model can explain the experiment without need to consider non-Fickian diffusion (see section 3.3.6). Wang et al. (2016) introduce a transient transfer function which is based on discrete fractures. Wang et al. (2020)

derives a new scheme based on the Generalized Multiscale Finite Element Method (Efendiev et al., 2013) and establish a triple-continuum model for fractured-vuggy carbonate reservoirs. The method improves efficiency. This is accomplished by systematically generating an approximation space through solving a series of local snapshot and spectral problems. The resulting eigenfunctions can pass the local features to the global level when acting as basis functions in the global coarse problems.

### Outline of simplified unit cell

<figure><img src="image_16.png" /><figcaption>FIGURE 12
Example of small-scale cell used for calculation of macro-scale dual porosity parameters (Ceriotti et al., 2019).</figcaption></figure>

### 3.4.2.2 Small scale models/upscaling parameters

Chen et al. (2019) develop a method to measure the dual porosity parameters, in relation to hydrocarbon sources in shales. It is improved with respect to the previous work through different upstream and downstream. Steiner and Mittermeir (2017) evaluate the matrix-fracture transfer rate, from a numerical calculation on a fine scale single matrix block model. They conclude that a general formula it not yet established but the suggested type is perspective for further generalisation. Ceriotti et al. (2019) introduce a technique to obtain two effective parameters of the upscaled dual porosity model through numerical solution of a simplified unit cell (Figure 12). Porta et al. (2015) evaluate experiments in “synthetic” porous media (glass beads) and uses a pore-scale 3D image. They resolves the former limitation that relying solely on
information about the geometrical attributes of the pore space obtained through NMR and electric resistivity measurements does not lead to a reliable characterization of the full set of model parameters.

### 3.4.2.3 Numerical methods

Ghahfarokhi (2017) shows that a structured mesh orientation has effect on the Oda's method accuracy even when the known correction is used. Wang et al. (2017) applied the “proper orthogonal decomposition”, which is useful for model reduction and introduced a new version with mass conservation. Next, Wang et al. (2018) improves the solution algorithm to ensure the global mass conservation. Karimi-Fard and Durlofsky (2016) introduced a special technique of fine-grid discretization, upscaling, and coarse-grid discretization, which can be used with “arbitrary” simulator for accurate solution with little computing requirement. Chung et al. (2018) provide rigorous and accurate non-local solution (in the oversampled region).

### 3.4.3 Hybrid-scale approach

Reactive transport modeling has established itself as a powerful, versatile and essential tool for tackling the complexity of the multi-space and temporal-scale issues engineers and scientists face. Active topics of research dealing with reactive transport modelling include the development of pore scale and hybrid, or multiple continua, models to capture the scale dependence of coupled reactive transport processes. At the continuum scale, the reactive transport constitutive equations rely on the description of the porous medium with respect to its macroscopically measurable properties (e.g., permeability, diffusability). At the pore scale each point of the space is occupied by either a fluid or a solid phase and macroscopically measurable properties do not apply. Coupling the pore scale and the continuum scale is a challenge and recently porous Media4Foam plateform (Soulaine et al., 2021) has been developed in order to do it in a smart manner using the Darcy-Brinkman-Stokes equation. The developed platform was validated against existing benchmark. At the continuum scale the reference benchmark have been published in Xie et al. (2015) and Poonoosamy et al. (2021). The obtained results have been published in Pavuluri et al. (2022). At the pore scale, the reference benchmark was the one published by Molins et al. (2021) and the result discussed in Soulaine (2021). At the hybrid-scale the ability of porous Media4Foam modelling framework was demonstrated simulating dissolution and precipitation processes in fractured porous media. Here, the reactive medium consisted of celestite grains that reacted with a barium chloride solution injected into the system, leading to the dissolution of celestite and the growth of barite. Differences in mineral precipitation has been observed in dissolution patterns by varying the injection rates.

### 3.4.4 Coarse graining

Coarse-grained models aim at simulating the behavior of complex systems using a simplified (coarse-grained) representation. Coarse-grained approaches are widely used in atomistic/molecular modeling. Simplification is essential, as the large numbers of atoms/molecules cannot be represented. Coarse-grained models are also often used for multiscale modeling. Coarse-grained approaches are attempting to bridge neighbor scales in a consistent continuum-mechanical way.

We provide a few recent examples to elucidate the *status quo* and potential of the multiscale methodology (Tahmasebi and Kamrava, 2018). Recently, Churakov and Prasianakis (2018) presented a review on a holistic process-based approach for reactive mass transport with an emphasis on mineral reactions in porous media. The paper reports the current *status quo* of mineral-fluid interface chemistry at molecular level and develops a road map for the general challenge: How to translate atomistic phenomena into macroscopic levels for practical model applications?

Gurikov et al. (2016) present a cellular automaton model for large-scale simulations at molecular level. They attempt to describe mass transport processes in porous media, i.e., diffusion, adsorption and directed flow in porous media, with large representative elementary volumes (REV). The framework could be successfully used to model diffusion processes in disordered porous structures with phase separation phenomena. Essential features of the molecular processes could be represented for physically meaningful time scales.

Recently, Schmidt and Steeb (2019) introduced a multiscale approach for hydro-mechanical coupling of fluid-filled fractures using a hybrid-dimensional element formulations with non-conformal meshes. In their hybrid approach they use parabolic velocity profile with the fracture to capture the fluid mechanics behavior at small scale and avoiding an explicit discretization. The fracture dimension is reduced by one order and implemented as interface elements to the modeling domain. The mass and momentum balances are described by modified balance equations, respectively. Different coupling schemes (staggered vs. monolithic) have been implemented and tested for testing the accuracy of the hybrid scheme. A number of test cases have been used for validation of the newly developed approach. Alyaev et al. (2019) using a similar approach for simulating single-phase flows in porous media with inertial effects based on a control volume heterogeneous multiscale method.

Coarse-grained methods are also used in related areas such as sediment transport, when modeling flow over mobile porous media. Biegert et al. (2018) couple particle and fluid phase flow with a coarse-grained approach. They conducted a grain-resolved direct numerical simulation of multiphase flow for the coupled particle-fluid systems. The particle phase is introduced by a boundary method where particle interactions can be taken into account via a discrete element method. The governing equations for fluid flow and stress budgets are derived by double averaging procedures.

### 3.4.5 Porosity permeability alterations in reacting environments and upscaling approaches

Precipitation and dissolution reactions in porous media dominate and control a large number of geochemical processes and industrial applications. The precipitation and dissolution of minerals from aqueous solutions alters the pore space and its connectivity. This has a strong effect on the mass convection and diffusion through the porous medium. For example, when a mineral precipitates/dissolves on the reactive porous surface, the overall porosity decreases/increases leading to a subsequent decrease/increase in permeability and effective diffusivity. At the same time the connectivity of the pores can also change in a way to block or to facilitate the mass diffusion processes.
<figure><img src="image_17.png" /><figcaption>FIGURE 13</figcaption></figure>

Evolution of permeability for different reactive transport dissolution scenario in a porous medium. Porosity permeability correlations can be extracted via fitting for use with macroscopic codes. Blue points represent the random set of geometries permeabilities that can serve as starting points (star). A geometry with 20% porosity is selected. Temporal evolution of the structure for three different cases were Peclet and Damkohler number vary.

It is common that the field scale simulations (continuum scale) use a simplified description of the processes that actually take place at the pore-level scale (microscale). For example, the change in permeability and diffusivity of a porous medium due to reactions is correlated to a global change of the bulk porosity in the medium, using empirical laws. It has been demonstrated, that smooth simplified relations are frequently inadequate to predict accurately the temporal evolution of the system of interest (Steefel et al., 2015a).

The predictability of the models can be greatly improved after resolving experimentally and numerically the physical phenomena and process mechanisms occurring at the microscale, followed by an upscaling procedure to the continuum scale (Noiriel and Daval, 2017; Deng et al., 2018a; Prasianakis et al., 2018; Molins et al., 2019). For example, depending on the conditions, a calcite rock can be dissolved following a face dissolution mechanism or via wormhole mechanism (Hoefner and Fogler, 1988; Portier et al., 2009; Menke et al., 2016). The resulting porosity-permeability correlation for the same rock can be greatly different (Prasianakis et al., 2018). The flow conditions can be described by the Peclet and Damkohler numbers. In Figure 13 such an example is shown, where the LBM framework has been used. Different flow conditions on the same rock have been implemented, which resulted in distinct permeability evolution paths.

Capturing the dynamic processes of porosity clogging has been the subject of several studies where pore-level simulators have been used in conjunction with computer generated pore and fracture geometries (Fazeli et al., 2019; Fazeli et al., 2020; Wang et al., 2023), or with specialized experiments (Lönartz et al., 2023).

At the moment, it is computationally very challenging to extract such correlations, since three dimensional reactive transport calculations are needed.

Upscaling of mineral precipitation and dissolution processes in porous media affecting the pore structure and subsequently the flow properties in the continuum scale is a challenging task to be considered in mathematical models. To model such processes, one needs conservation laws for mass, momentum and solute in time-dependent domains where the evolution of the sharp interface separating two domains (mobile and immobile) is not known a priori. Therefore, a scheme accounting for moving/free boundary conditions, where the boundaries and hence the domains evolve, needs to be considered.

In modelling such coupled processes, two spatial flow scales are normally considered: pore (Stokes/Reynolds) and macro (Darcy) scale (van Noorden, 2009). Precipitation/dissolution reactions are considered using reactions rates where multiple chemical components could play a role. One of the most promising numerical schemes of recent is to use the phase-field model to represent moving sharp interfaces (Xu and Meakin, 2011) and it has been extended to porous media (Bringedal et al., 2019).

# 4 Tools and methods to quantify/derive uncertainties induced by coupled processes

## 4.1 Understanding of uncertainties in disposal systems and quantification of uncertainties

### 4.1.1 Introduction

Treatment of uncertainties in the performance assessment of geological nuclear waste disposal has been recognized as an
important issue since the development of the concept in the late 80's (Bonano and Cranwell, 1988). Different kind of uncertainties are described such as the one related to the datasets, the modelling approach itself, the conceptual uncertainty, etc. This is why in EURAD a strategic study entitled "Understanding of uncertainty, risk and safety (UMAN)" has been conducted. It has allowed (i) to define and categorize uncertainties (Brendler and Pospiech, 2023), (ii) compile, review, compare and refine strategies, approaches and tools for the management of uncertainties in the safety analysis and the safety case (Plischke and Rohlig, 2023), (iii) develop a common agreed understanding among different actors of national disposal programmes on strategies and approaches for uncertainties management by sharing knowledge and experience (Kaempfer et al., 2023; Haverkate et al., 2024). If uncertainties can be epistemic or aleatory, sources of uncertainty cannot always be clearly assigned to one of the two types of uncertainty (Brendler and Pospiech, 2023) and glossary inside for a clear definition.

The subject of couplings between the thermal, hydrological and mechanical processes in rocks has become an important subject mainly due to the modeling requirements for the design and performance assessment of radioactive waste repositories and other engineering fields such as gas/oil projects, geothermal energy exploitation, and environment impact evaluation in general. Coupled THM (thermo-hydro-mechanical) and coupled THMC (thermo-hydro-mechanical and (geo)chemical) processes are the ones most often required in the numerical models (Hudson et al., 2005; Hudson et al., 2017).

The requirement of coupling the processes in a geological system can be derived from processes that affect each other, which occurs mainly in two ways. The first is direct coupling in which one process induces directly the development of another process, represented by a cross term involving both processes in the governing equations. The other way is indirect coupling in which one process changes the parameters controlling another process.

A key process in the early stage of repository construction is the creation of an excavation damaged zone (EDZ) resulting in increased permeability around the tunnels. Performance assessment requires the knowledge of spatial-temporal evolution of EDZ properties. During the Exploitation Stage, decay heat from radioactive waste may induce large temperature gradients in the buffer zone and in the near-field rock. At the same time, resaturation in the near-field occurs which is a slow process controlled by low rock permeability and buffer characteristics. Heat generation may also lead to water vapor convection near the heat sources. This process induces transient pore pressure build-up resulting from differential thermal expansion. The pressure build-up affects effective stresses in the engineered and natural barrier system. Heating may impact rock strength and may induce chemical reactions with impact on rock properties as well. The Long-Term Post-Closure Stage is the subject of long-term performance and safety assessment. The key process is the likely failure of the engineered barriers that may lead to the release and subsequent transport of radionuclides and any nonradioactive contaminants. The transport characteristics in the tunnel backfill and near-field rock during this period depend on the coupled processes and effects that occurred during the previous stages.

Main sources of uncertainties related to radwaste disposal and their assessment that numerical modelers must consider:

* The extremely long (geological) timescale over which system evolution must be evaluated and predicted. Thus, the long-term evolution of safety functions of the repository cannot be well supported by comparably short-term laboratory, in situ and/or URL investigations.

* The inherent variability of geological features makes the characterization of host rock formations challenging at the required level of detail based on small-scale experiments and laboratory tests.

* One must consider the coupled behavior of thermo-hydraulic-mechanical as well as (geo)chemical processes to be investigated with respect to spatial scale and time frame of investigation.

* The perturbation effect of relatively rapid (geological) changes, e.g., climate change, earthquake as well as heat generation on the long-term evolution of the repository system. From a numerical modelling point of view, this requires either hybrid modelling approaches or the availability to dynamically change the time-stepping algorithms of the models.

* The applicability of the models used for assessing the system evolution. Thus, the risk related to uncertainty regarding the conceptual model used to describe interconnected physical processes must be evaluated.

In this state-of-the-art report, we provide a brief overview of coupled processes and the associated uncertainties with respect to underground nuclear waste repositories. First, we summarize the relevant THM and THMC processes and their general method of resolution. This is followed by characterizing the uncertainties of these phenomena. The report ends with conclusions regarding the status of handling these uncertainties.

## 4.1.2 Descripion of the coupled processes in the radioactive disposal system

### 4.1.2.1 Thermal-hydro-mechanically and chemically coupled processes and their resolution

One of the main difficulties in the prediction of the spatial-temporal evolution of a radioactive waste disposal is to consider all relevant processes and their couplings (Millard, 2017). These include mechanical, hydraulic, thermal as well as geochemical processes which may be extended to gaseous and biological phenomena. The THMC-processes and their connections are illustrated in Figure 14.

The basis of mechanical processes is that rock mechanical properties are changed due to repository excavation; and the rock mass may contain water. Thus, under loading, the stresses developing in the rock medium are linked to both solid and fluid components, leading to the concept of effective stresses.

The hydraulic part depends on rock saturation. In the simplified case of a saturated rock, water flow is not only controlled by pore pressure gradient (Darcy flow), but also by pore volume change induced by rock deformation. Water expansion due to temperature change may be considered as well.

The evolution of temperature is governed by heat release from the waste. This is described by heat diffusion due to developing temperature gradients and advection associated with fluid flow. Chemical processes are dominated by dissolution and precipitation enhanced by thermal gradients.
<figure><img src="image_18.png" /><figcaption>**FIGURE 14**
Schematic illustration of the coupled thermal-hydro-mechanical-chemical processes relevant to safety assessment of nuclear waste repositories (modified from Gaucher et al. (2015)).</figcaption></figure>

Changes in water pressure and chemical reactions affect water and gas flow paths and contaminant transport characteristics.

Gaseous processes include gas generation, pressure build-up as well as gas lift. Furthermore, extensive pressure build-up may lead to the failure of some parts of the Engineered Barrier System (EBS) further enhancing the release of radionuclides from the repository. Gas (mainly methane) could be generated as a result of biodegradation of organic material, hydrogen gas generation is foreseeable from the corrosion of metallic components and some gas could be generated from hydrolysis.

Consequently, the assessment of radioactive waste disposal evolution requires the knowledge of physical variables controlling thermo-hydro-mechanical-chemical processes at any time. Thus, the resolution of a coupled system of partial differential equations is required. There are two types of approaches to resolve such problems (Ismail-Zadeh, 2010).

In the first type, all variables are solved simultaneously. In simultaneous solution methods, all the equations are considered as part of a single system. However, direct solution of these equations is time consuming, especially when these equations are

non-linear, and the problem is three-dimensional, which is the case generally for radioactive waste disposal systems. When the equations are linear and fully coupled, the simultaneous approach is the best.

However, the second method, the iterative approach may be computationally much more efficient. The equations can be so complex and non-linear that coupled methods are too challenging to apply. It may then be preferable to treat each equation as if it has only a single unknown explicitly, temporarily treating the other variables with the currently available values for them. The equations are then solved in turns, repeating the cycle until all equations are satisfied. Since some terms, e.g., the coefficients and source terms that depend on other variables change as the simulation proceeds, iterative solvers are much more efficient and so these are preferred. Iterations performed on each equation are called inner iterations (e.g., explicit integration scheme for solving the equation of motion). In order to obtain a solution that satisfies all the equations, the coefficient matrices and source vector must be updated after each cycle and this process is then repeated. These cycles are called outer iterations.
### 4.1.2.2 Sources of uncertainty in coupled processes associated with radioactive waste repositories

For the performance assessment of a radioactive waste repository, analyses of the coupled THM(C) phenomena must be considered. However, many of these processes are non-linear and the constitutive relationships between them are not well known. Thus, various constitutive relationships are used that are applicable for different geological settings. What coupled THM or THMC processes should be included in a simulation model depend on site-specific safety criteria and scenarios. Therefore, one should characterize uncertainty arising from conceptual model approaches.

Birkholzer et al. (2019) point out that the various sources of uncertainty associated with numerical modelling of evolution of radioactive waste repositories lead to the realization of validating numerical modelling results against laboratory and field data. For this purpose, several laboratory experiments have been carried out using brendler samples of different sizes ranging at cm and m scale. Field measurements are also conducted at experimental sites, either at locations where repositories are proposed to be situated or at areas with representative geologic conditions. Several subsurface facilities, referred to as underground research laboratories (URL) are designed to conduct such tests. The goal of these tests is to study coupled processes and check modeling results often at the spatial scale envisioned for the waste emplacement. Numerical simulations with different modeling approaches are then performed to test the capability of the model against measured data. The goal of these simulations is to improve modelling approaches, to increase model confidence as well as to estimate uncertainty of various parameters governing coupled processes. However, it must be pointed out that geological settings and site characteristics cannot be known completely due to limitations associated with site observations and field tests. Another source of uncertainty lies in different key processes governing the evolution of the disposal system at different stages of the repository and spatial scales (Figure 14). Birkholzer et al. (2019) and Tsang et al. (2012) summarize key processes and issues for different phases of repository lifetime, such as early repository implementation phases (Construction and Open Drift Stage), the period of thermal disturbance after waste emplacement (Exploitation Stage) and a long-term period where much of the early perturbations are no longer valid and the system is approaching a late time equilibrium (Long-Term Post-Closure Stage).

### 4.1.3 Uncertainties of coupled THMC models for the EBS and needs for further research

THMC models of the EBS have conceptual model uncertainties, regarding the identification of the relevant processes, parameter uncertainties and uncertainties in the initial and boundary conditions. Usually, these uncertainties are analysed by means of sensitivity analyses. The relevance of physical and chemical processes are analysed by performing simulation runs with and without the consideration of a given physico-chemical process. Similarly, the uncertainty in THMC parameters is commonly ascertained by performing systematic sets of sensitivity runs in which the values of the parameters are increased and decreased one-at-a-time.

Amongst others example of the remaining uncertainties of the coupled THMC models of the EBS that need further research include:

1. Multiphase flow:

a. The complex porosity structure of the bentonite which requires the use of dual and triple continuum models (Sanchez et al., 2005; Sanchez et al., 2012; Sanchez et al., 2016b)

b. The retention curve of the bentonite: hysteresis and temperature dependence

2. Reactive solute transport

a. Relevance of reactive gases such as CO₂(g) and O₂(g)

b. Solute diffusion from the EBS into the host rock

c. Alternative conceptual models of solute diffusion in the EBS:

i. The single interlayer-porosity Donnan equilibrium model (Birgersson, 2017)

ii. The electrochemical model based on the Nernst-Planck equation (Idiart and Pekala, 2016)

iii. The multiple porosity (interlayer, double layer and free) model (Samper et al., 2009; Zheng and Samper, 2015).

iv. The use of different effective diffusion coefficients for each chemical species while preserving charge balance

3. Geochemistry

a. Precipitation of the secondary mineral phases at the canister-bentonite interphase such as Fe-phyllosilicates (chlorite, cronstedtite and berthierine); Fe-rich smectites and zeolites such as phillipsite, chabazite and merlionite;

b. Precipitation of secondary mineral phases at the concrete-bentonite interphase such as CSH, CASH, zeolites such as phillipsite, chabazite and merlionite, and iron silicates

c. Fe reduction in the bentonite

d. The uncertainties in the kinetic rate laws, rate constants, reactive surfaces and catalytic effects of kinetic mineral dissolution/precipitation. The dynamic update of the reactive surface areas of the minerals;

4. Coupled effects

a. Changes in porosity and other transport and chemical parameters caused by mineral precipitation/dissolution reactions

b. Simulating canister corrosion in a more realistic manner by adopting a dynamic corrosion front;

c. The coupling of the evolution of the bentonite microstructure with the changes in the composition of the exchange complex

To date Global sensitivity and uncertainty quantification methods have been widely used in performance assessment studies to quantify the prediction uncertainties caused by uncertainties in flow and transport parameters (Bonano and Cranwell, 1988; Helton et al., 2012). A few studies have analyzed the uncertainties in using reactive transport models. Matsumoto and Neyama (2008) developed a preliminary uncertainty analysis tool for performance assessment calculations considering uncertainties in the groundwater concentrations of Na, Ca, HCO₃ and Cl. Criscenti et al. (1996) and Denison and Garnier-Laplace (2005) applied Monte Carlo methods for geochemical calculations. Their focus were: (a) the significance of uncertainties in predicted pH values for a simple buffer solution and for calcite saturation indices, and (b) the effects of database parameter uncertainty on uranium (IV) equilibrium calculations (i.e., Uranium (IV))
speciation), respectively. Ceriotti et al. (2017) presented a model to estimate the uncertainty associated with the amounts of CO₂ generated by carbonate-clays reactions in subsurface systems. Recently, Ayoub et al. (2020) presented a comprehensive uncertainty and sensitivity analysis of Cs sorption reactive transport modeling with 3 cation exchange sites in the near field of a deep geological repository for nuclear waste based on the Morris and Sobol method. Ayoub et al. (2020) analyzed the uncertainties in: 1) Cation exchange selectivities and 2) Clay pore water cation concentrations. They were able to identify the most important uncertain parameters affecting the transport of Cesium and the combination of parameters values leading to the maximum Cesium concentrations at a specified location using classification trees. Finally, they proposed a Cesium isotherm and a K<sub>d</sub> uncertainty range based on a large number of numerical simulations.

### 4.1.4 Approximation of the probability density function

The classical Monte Carlo is an attractive method for uncertainty propagation due to its noninvasive nature. However this flexibility is offset by slow convergence, leading to high computational costs. The Multi-Level Monte Carlo (MLMC) method introduced by Giles (2008) provides an efficient way to estimate the expectation of a variable resulting from an approximation. Basic idea is to do high volume of samples using a less accurate but computationally cheaper approximation, complemented by far fewer samples of the difference between the observable and its approximation. If this difference has a significantly smaller variance, we obtain an estimate of the same accuracy but at the fraction of the computational cost. The MLMC method has been already studied and analyzed in connection with variety of PDE problems: elliptic equations (Barth et al., 2013; Cliffe et al., 2011; Abdulle et al., 2013), parabolic equations (Barth et al., 2013), multiphase flow (Muller et al., 2013; Lu et al., 2016). Recently, an MLMC approach consisting of a neural network metamodel and several PDE levels has been investigated for Darcy flow (Špetlík and Březina, 2021) and solute transport (Špetlík and Březina, 2022) problems, demonstrating decreased cost in certain cases. The MLMC method offers potential as an efficient estimator for variance-based sensitivity analysis methods such as Sobol, discussed in detail in Mycek and De Lozzo (2019). Unlike performance assessment (PA) based on metamodels, PA utilizing multiresolution or multifidelity models is advantageous for its comprehensive characterization of total estimator error, which implicitly includes the approximation error. Multiresolution description of fractured media, is feasible via Discrete Fracture-Matrix approach, as demonstrated in simulators like Flow123d. However, the effective construction of coarse approximation levels relies on fast and robust numerical upscaling methods. The upscaling of Darcy flow using neural networks, as developed by Špetlík et al. (2024), could facilitate the application of MLMC to practical 3D performance assessment problems.

A methodology for addressing uncertainties in a broader context as required in PA will be confronted with (as outlined above) fragmentary information, with heterogeneous spatial (and sometimes also temporal) geochemical composition (minerals, solution, hydrodynamic characteristics), and even with occasional scientific dissent about governing processes and other boundary

conditions. This applies (with increasing magnitude) for the parameterization of thermodynamics of solids (and solid solutions), ion-ion interactions, and properties of accessible/ reactive surfaces. Strategies to overcome such problems are very diverse and depend to a high degree on the efforts necessary to perform additional lab experiments or field sampling/ characterization. Ideally it builds on previous knowledge of phenomenology and processes in the respective geochemical systems. The range of strategies given below is sorted (approximately) from very “cheap” to “expensive” - in the sense of fiscal, temporal and personal resources requested. Starting at the 2nd alternative, a variety of pdfs is usually considered. For (namely, aqueous) systems with sufficient independent parameter determinations Gaussian (or log Gaussian) distributions or as. more robust alternatives trigonal or uniform pdfs are used, the latter often being the only sensible approach when dealing with bimodal parameter distributions. Moreover, in the not so exotic case of numerical data themselves only being estimations (or even expert guesses), only the first three tiers are applicable at all.

1. “Do nothing” option: ignore the respective uncertainty as the underlying process may hopefully not be essential (should of course be verified) and the harm of ignoring may be smaller than using some estimates not really defensible.

2. “Expert guess” which can be anything between pure assumptions and the application of well-established informal rules based on general geochemical principles.

3. “Estimates” stemming from pure correlations in the sense of empirical fits, over chemical analogs (an often exploited analogy covers Ln (III) vs. An (III)), the usage of structure-based incremental systems, linear free energy relationships, or quantum chemical approximations (to name just the most prominent examples).

4. “Re-computation” of parameters from other (hopefully similar) models (e.g., going from the Pitzer-approach to SIT parameters and vice versa–this direction being the more challenging one).

5. “Re-fitting” of a pool of already existing raw data with uncertainties as by-product.

6. Setting up an experimental program to derive the missing parameters in the lab.

### 4.1.5 Examination of uncertainties of sorption of radionuclides in rocks

Examination of uncertainties of sorption of radionuclides in rocks.

In the following the international status with respect to the assessment of uncertainties within the field of sorption of radionuclides in rocks is highlighted, sketching also some possible development lines. The first three sub-chapters describe the types of uncertainty usually encountered in conventional geochemistry. Thereafter, consequences with respect to predictive modelling applications are outlined.

#### 4.1.5.1 Aqueous speciation

It is common to encounter uncertainties in thermodynamic data (e.g., in a recent europium data compilation Jordan et al. (2004) illustrate the effect of experimental errors on thermodynamic
parameters and activity coefficients), with the majority of studies focusing on the reproducibility of experimental data. However, systematic errors, which are more challenging to quantify, are often overlooked. Furthermore, the formulations of chemical reactions, which describe the reactions from educt to product, are often based on mere postulations. It is possible that a chemical system can be described by two mutually incompatible sets of species that represent experimental findings equally well. This results in an overlap of numerical and conceptual uncertainty. Furthermore, in natural systems, some reactions are strongly kinetically hindered, especially if they involve the restructuring of molecules or ions. This is regularly observed in redox systems, for example, involving Pu, Se, U or Tc species. Finally, the uncertainty related to thermodynamic data at temperatures higher than 25°C is also of major concern. In this context, there is a paucity of data, particularly in the form of calorimetric studies or speciation experiments. Furthermore, the temperature dependence is often described by purely empirical temperature-dependent-functions.

With regard to the parameters required in activity models to describe ion-ion interactions in highly saline solutions, the situation differs between the two major approaches. While in the SIT (Specific Ion Interaction Theory–Ciavatta (1980)) the most comprehensive compilations contain uncertainties Rand et al., the situation in the Pitzer model (Pitzer, 1991) is more complicated. Due to the high connectivity and interdependence of the various types of binary and ternary Pitzer coefficients, it is uncommon for uncertainties to be provided. This is because there is currently no established procedure for determining sensible uncertainties for them. This may be a concern in highly saline environments such as salt rock or certain clay formations, whereas in all other host rock scenarios, the SIT approach is the most appropriate.

### 4.1.5.2 Solids (minerals)

The situation here is analogous to that described for aqueous species. A further source of uncertainty is the exact stoichiometry of the solid (minerals in nature are rarely “pure”) and their physical properties, such as specific surface area (SSA). This is connected to the problem that dissolution experiments to yield solubility constants usually start from well-crystalline, aged phases. In nature, however, many of these phases do not form at all under hydrothermal conditions or undergo a transition from amorphous intermediates (with much higher solubilities) to eventually well-crystalline ones. Thus, it is helpful to distinguish between those processes that involve the dissolution of “matured” well-crystalline phases and the precipitation (and dissolution) of “fresh” amorphous phases. With regard to the solubilities of a variety of phases, values have been determined for both variants. These values can be used for the purpose of establishing correlations to close respective data gaps.

With regard to temperature-dependent-function, the situation is more favourable than for aqueous solutions (and considerably more favourable than for surface species), given the availability of a wealth of enthalpic and entropic data for solid phases.

In addition to the aforementioned uncertainties regarding mineral phases, spatial heterogeneities (including accessible surfaces) within the rocks may also contribute considerably to the total uncertainty retrieved from the solid phase. While the majority of uncertainties relevant for geochemical processes, especially on the mineral surfaces, are epistemic, a

few-specifically those associated with mineralogical heterogeneity-are aleatoric. Other uncertainties where additional experiments may not reduce uncertainties include pore size distributions, accessibility of pore volumes or reactive surface areas. In conclusion, the spatially heterogeneous distribution of minerals, their grain sizes, exact chemical composition, their exposure to voids or fluid pathways within the rocks, and the probability of fluids actually passing along their surfaces, may present significant challenges in even estimating the uncertainty resulting from them.

### 4.1.5.3 Surface speciation

Regarding surface speciation, the publication of uncertainties is not yet well established in this field. If it were, the same restrictions would apply as already outlined above for aqueous and solid species. Surface reactions are further limited by the existence of a variety of competing models describing sorption via surface complexation models and ion exchange. These include non-electrostatic, diffuse double layer, constant capacitance, triple layer or basic stern models (Stumm et al., 1992). The respective parameters are usually not directly convertible, thus the parameter pool to derive probability density functions (pdfs) is diminished. Another drawback is the relation to mineral characteristics being themselves not well established, e.g., surface protonation or again specific surface area, or surface binding site density. An overview of the data situation can be found in Ochs and Payne (2012) as well as in the mineral-specific sorption database REST (Brendler et al., 2003). Data scatter evidenced as multi-modal distribution patterns is often related to alternative experimental methods. A classical case is the determination of SSA of clay minerals by either BET measurements or via the EGME method.

### 4.1.5.4 Established uncertainty assessment methodologies

Despite the difficulties encountered in geochemistry, there is at least a well-established procedure for treating parameter uncertainties. The most recent edition of the Thermochemical Database of the Nuclear Energy Agency (NEA TDB), published in 2013 and 2020 (Lemire et al., 2013; 2020), contains in Annex C, entitled “Assigned Uncertainties”, comprehensive procedures for computing uncertainties when combining or averaging several independent sources of parameter values. Finally, the derivation of uncertainties when re-formulating chemical reaction equations (to accommodate, for instance, switches in the master species set required in different speciation codes) must follow the common procedures of error propagation. The uncertainties provided for the thermodynamic values recommended by the NEA TDB are provided in a consistent and traceable manner, thus setting international standards. The same can be said of the values provided by the most recent PSI/Nagra TDB (Thoenen et al., 2014), which to a large extent actually relies on the NEA TDB but follows a different approach with respect to SIT coefficients. Other TDBs that are frequently employed in the context of nuclear waste disposal include ThermoChimie (Giffaut et al., 2014), the Yucca Mountain Project TDB (Gdowski et al., 2002), and the HKL TDB (Helgeson et al., 1981). These TDBs either adhere to the general guidelines established by the NEA TDB or do not report any uncertainties.
### 4.1.5.5 Predictive modelling

It is also worth noting that, at present, no geochemical speciation or reactive transport code is utilising the available information about uncertainties. In order to do so, separate routines must be written, "enclosing" the appropriate code and adding such functionalities (Stockmann et al., 2017). The transfer of geochemically induced uncertainties to PA is typically accomplished through the use of lookup tables (e.g., for solubilities or distribution coefficients). More recent developments have sought to address the issue of multidimensionality by employing uncertainty and sensitivity analysis, the introduction of surrogate functions, neural networks, or other machine learning approaches. For further details, please refer to section 4.2 and section 4.3.

A particular challenge arises in the treatment of parameters that exhibit a high degree of dependence. In geochemistry, the correlation matrices derived from multiple parameter fits are typically not considered. This may result in the aggregation of uncertainties, which can lead to significant overestimation of combined uncertainties (Meinrath, 1997).

### 4.1.6 Derivation of procedures for pdf of corresponding parameters for chemical models

Wide set of data for radionuclide transport processes is available for central European crystalline massive conditions. Specific conditions of transport processes specific for Bohemia crystalline rock massive conditions were described in Havlová (2013), defining that Central European crystalline rock massifs are specific by groundwater of lower salinity and rock structures and textures less distorted by orogenic long-term processes in comparison with Swedish and Finish crystalline rock or with Alpine granites as in Grimsel Test Site. Studies of radionuclide transport properties on crystalline rock has been performed in the Czech Republic since late 90's (Passas et al., 1996; Buňatová, 1998; Havlová, 1999) being focused namely, on radionuclide sorption ($Cs^i$, $Sr^{II}$, $Np^{III}$, $I^{-1}$) on generic granites under synthetic groundwater conditions. Those studies, mostly using batch sorption experiments on crushed materials, were complemented by studies focused on aerobic/ anaerobic behaviour of redox sensitive species (Se, Tc), using more comprehensive methods, like surface scanning and spectroscopic studies for element (Se) specific sorption at the rock surface (Videnská and Havlová, 2012; Videnská et al., 2014b; Videnská et al., 2014a).

Migration of radionuclides in dynamic system is described with the transport parameters such as the retardation and distribution coefficients (R and $K_d$), Peclet number (Pe) or hydrodynamic dispersion coefficient (D). 2D- or 3D-models deal with longitudinal hydrodynamic dispersion coefficient. These parameters can be determined using dynamic column technique, which has many advantages in comparison with static batch experiments. Namely, it enables: (i) direct determination of contaminant retardation coefficients dissolved in aqueous phase, (ii) closer approximation of the real conditions in the natural system, (iii) study of less disturbed samples, (iv) study of desorption process following immediately after sorption, and so on. Retardation of both non-sorbing species ($HTO$, $Cl^{-1}$, $Se^{IV/VI}$) and sorbing species ($Sr^{II}$) in dynamic granitic column system has been presented (Palagyi et al., 2013; Videnska et al., 2013; Videnská et al., 2015; Štamberg et al., 2014).

Systematic approach toward potential site selection brought into the focus some of the potential DGR sites since 2015. The most comprehensive approach to the rock condition characterisation under conditions close to the DGR emplacement in the depth of 550 m brought Bukov URF survey, including a wide set of sorption and diffusion experiments for several radionuclides ($Cs^i$, $Sr^{II}$, $Se^{VI}$, $U^{VI}$, $HTO$, $Cl^{-1}$, $I^{-1}$ - Bukovská and Verner (2017). The above mentioned study has been complemented with more generic diffusion studies, being performed either in generic Czech massive granite samples (Havlová et al., 2012) or on samples from Grimsel Test Site (Soler et al., 2015).

Direct tool for potential assessment of fracture flow quantification and sensitivity analyses is available in ÚJV laboratory. The physical model of granite block (dimensions: 80×50×40 cm) was assembled and the generated artificial fracture was characterised by means of 3D laser scanner. The goal of the tracer tests, using different types of tracers ($Cl^{-1}$, $I^{-1}$, $HTO$, $Pb^{II}$) has been to develop tools for evaluation of migration and retention of potential contaminants in the rock environment as the key input values for the safety assessment of anthropogenic activities and also to enable test mathematical codes and assess the uncertainty of reaction transport processes determination. Subsequently the block was instrumented with on-line measuring devices (conductivity and pressure sensors), allowing performance of migration experiments. Detailed description of the fracture parameters (Figure 15) serves as a key input for modelling purposes (Jankovec et al., 2019; Jankovský et al., 2019a) using both conventional software (MODFLOW, FEFLOW) and in-house code (Flow123d). However, natural fracture flow represents more complex environment, including surface coatings and miscellaneous topographic domains within the surface (Jankovský et al., 2019b; Missana et al., 2008; Huber et al., 2014).

However, laboratory experiments are not able to reflect the real conditions in the rock massif all the time. Therefore, the tracer tests both under static or dynamic conditions *in situ* are to look into upscaling the transport processes. Such an experiment can be focused both on studying radionuclide diffusion into the rock matrix (Soler et al., 2015) or on advective flow experiments in fracture, including modelling in MODFLOW/MT3D code (Zuna et al., 2017).

## 4.2 Methods to evaluate uncertainties

### 4.2.1 Adjoint state methods

Adjoint state (AS) methods have been widely used in optimization and inverse problems in: 1) Groundwater flow models (Townley and Wilson, 1985; Sun, 2013); 2) Conservative solute transport models (Samper, 1986; Ackerer et al., 2014) and 3) Geostatistics (Samper and Neuman, 1989). The adjoint state variables play the role of Lagrange multipliers. The advantages of the AS method include: 1) The calculation of the adjoint state requires a single calculation of the forward problem, regardless of the number or parameters; 2) The method is well suited for computing derivatives or local sensitivities for a large number of parameters; and 3) The AS equations can be derived for the continuous formulation or for the discretized numerical equations. Delay et al. (2017) presented a comparison of discrete and continuous adjoint states for groundwater flow in