Impact statement
This study provides key contributions in processing real-world geomechanical monitoring data from deep geological cavities, applying GP-based global sensitivity analysis using time-dependent Sobol indices, and achieving efficient calibration of geomechanical models. The structured approach demonstrates that only a few critical material parameters need calibration to accurately reflect in-situ monitoring data within complex constitutive models of rock salt. This finding is particularly significant for safety-critical applications such as deep geological disposal, where precise modeling is essential for long-term safety and stability. The results emphasize the efficiency and accuracy of the GP-based surrogate model in simplifying the calibration process while maintaining high fidelity to real-world conditions.
1. Introduction
Deep repository material models are complex geological models that account for the mechanics of soil and rock, hydrological properties, thermal effects, and chemical interactions (Pitz et al., Reference Pitz, Grunwald, Graupner, Kurgyis, Radeisen, Maßmann, Ziefle, Thiedau and Nagel2023; Claret et al., Reference Claret, Prasianakis, Baksay, Lukin, Pepin, Ahusborde, Amaziane, Bátor, Becker and Bednár2024). They are highly parametrized and the associated numerical analyses are computationally expensive (Kurgyis et al., Reference Kurgyis, Achtziger-Zupani, Bjorge, Boxberg, Broggi, Buchwald, Ernst, Flügge, Ganopolski and Graf2024; Wojnarowicz et al., Reference Wojnarowicz, Madaschi and Laloui2024). Their complexity has prevented their widespread adoption for many-query tasks, that is, simulations to explore multiple scenarios and the effects of uncertainties on the deep repository prognosis (Kurgyis et al., Reference Kurgyis, Achtziger-Zupani, Bjorge, Boxberg, Broggi, Buchwald, Ernst, Flügge, Ganopolski and Graf2024). In this regard, we propose using Gaussian processes (GPs) as surrogates for high-fidelity geological models. GPs can approximate the outputs of complex simulations with much lower computational cost, enabling efficient calibration, design, and validation (Radaideh and Kozlowski, Reference Radaideh and Kozlowski2020; Myren and Lawrence, Reference Myren and Lawrence2021; Sung and Tuo, Reference Sung and Tuo2024). Speeding up the simulations with the use of a GP-based surrogate will not only allow for validated prognoses for the repository through the assimilation of data but also help in making informed decisions and enhancing the overall robustness and adaptability of the repository management process.
1.1. Constitutive models for deep geological repositories
Rock salt formations are one of the potential host rocks in Germany considered for secure long-term nuclear waste storage due to their distinctive mechanical and hydraulic properties (Bollingfehr et al., Reference Bollingfehr, Buhmann, Dörr, Filbert, Gehrke, Heemann, Keller, Krone, Lommerzheim, Mönig, Mrugalla, Müller-Hoeppe, Rübel, Weber and Wolf2017; StandAG, 2017). Ensuring the natural integrity of geological barriers is crucial for safety, therefore reliable numerical calculations that depend on advanced material models are essential (Bollingfehr et al., Reference Bollingfehr, Buhmann, Dörr, Filbert, Gehrke, Heemann, Keller, Krone, Lommerzheim, Mönig, Mrugalla, Müller-Hoeppe, Rübel, Weber and Wolf2017). The thermo-mechanical behavior of rock salt is typically assessed by means of laboratory tests, in particular, short-term triaxial compression tests and long-term creep tests (Langer, Reference Langer1985; Wittke, Reference Wittke2014; Fecker, Reference Fecker2018). To include all the characteristics of rock salt, various constitutive models have been developed to effectively simulate the hydraulic, thermal, and mechanical behavior of rock salt; refer, for example, to Schulze et al. (Reference Schulze, Heemann, Zetsche, Hampel, Pudewills, Günther, Minkley, Salzer, Hou, Wolters, Rokahr and Zapf2007) and Hampel et al. (Reference Hampel, Argüello, Hansen, Günther, Salzer, Minkley, Lux, Herchen, Düsterloh, Pudewills, Yildirim, Staudtmeister, Rokahr, Zapf, Gährken, Missal and Stahlmann2013, Reference Hampel, Lüdeling, Günther, Sun-Kurczinski, Wolters, Düsterloh, Lux, Yildirim, Zapf, Wacker, Epkenhans, Stahlmann and Reedlunn2022a). These constitutive models are able to capture for instance multiple creep phases, healing, and dependency on, for example, stress conditions, time, temperature, and humidity. Several different approaches to capturing the complex behavior of rock salt were compared and briefly described in past research projects (Hampel et al., Reference Hampel, Günther, Salzer, Minkley, Leuger, Zapf, Rokahr, Herchen, Wolters and Düsterloh2010, Reference Hampel, Herchen, Lux, Günther, Salzer, Winkley, Pudewills, Yildirim, Rokahr, Gährken, Missal and Stahlmann2016, Reference Hampel, Lüdeling, Günther, Salzer, Yildirim, Zapf, Epkenhans, Wacker, Gährken, Stahlmann, Sun-Kurczinski, Wolters, Herchen and Lux2022b). Among them is the constitutive model TUBSsalt, which was developed by the Institute for Geomechanics and Geotechnics (IGG, TU Braunschweig) and presented for the first time in Gährken et al. (Reference Gährken, Missal and Stahlmann2015). It has been shown in Hampel et al. (Reference Hampel, Lüdeling, Günther, Sun-Kurczinski, Wolters, Düsterloh, Lux, Yildirim, Zapf, Wacker, Epkenhans, Stahlmann and Reedlunn2022a), (Reference Hampel, Herchen, Lux, Günther, Salzer, Winkley, Pudewills, Yildirim, Rokahr, Gährken, Missal and Stahlmann2016), and (Reference Hampel, Lüdeling, Günther, Salzer, Yildirim, Zapf, Epkenhans, Wacker, Gährken, Stahlmann, Sun-Kurczinski, Wolters, Herchen and Lux2022b) that TUBSsalt, as well as other constitutive models, accurately capture the thermal and mechanical characteristics of rock salt. The application of such complex constitutive models requires expert knowledge and ideally data to infer the numerous model parameters. Therefore, this constitutive model is the focus of the presented calibration process.
1.2. Gaussian processes as efficient surrogates for computationally expensive models
GPs are non-parametric probabilistic models that use Bayesian inference to make predictions and learn from data. They are particularly useful for modeling complex input–output data relationships and for making predictions in situations where the data is noisy or incomplete (Kennedy and O’Hagan, Reference Kennedy and O’Hagan2000; Kennedy and O’Hagan, Reference Kennedy and O’Hagan2001; Gu and Wang, Reference Gu and Wang2018; Gramacy, Reference Gramacy2020; Teckentrup, Reference Teckentrup2020; Myren and Lawrence, Reference Myren and Lawrence2021). One of the key advantages of GPs over concurrent regression ansatzes, such as neural networks, is their ability to effectively handle uncertainty (Schulz et al., Reference Schulz, Speekenbrink and Krause2018) and their applicability in a low-data regime. GPs can generate probabilistic predictions of the model output, considering the uncertainty in both the input data and the model itself. This makes GPs a good alternative for calibrating complex and non-linear computational models (see, e.g., (Wu et al., Reference Wu, Kozlowski, Meidani and Shirvan2018; Mahdaviara et al., Reference Mahdaviara, Rostami, Keivanimehr and Shahbazi2021; Li et al., Reference Li, Hariri-Ardebili, Deng, Wei and Cao2023; Veasna et al., Reference Veasna, Feng, Zhang and Knezevic2023) and (Sung and Tuo, Reference Sung and Tuo2024) for a recent review on computer model calibration). Additionally, global sensitivity analysis for computationally expensive models, which assesses the influence of input parameters on model output, becomes computationally feasible with GP-based surrogates (Oakley and O’Hagan, Reference Oakley and O’Hagan2004; Marrel et al., Reference Marrel, Iooss, Laurent and Roustant2009; Srivastava et al., Reference Srivastava, Subramaniyan and Wang2017). GPs have been applied to sophisticated computer models in various fields, including nuclear physics (Kejzlar et al., Reference Kejzlar, Neufcourt, Nazarewicz and Reinhard2020), environmental science (Cheng et al., Reference Cheng, Konomi, Matthews, Karagiannis and Kang2021), and digital twining (Thelen et al., Reference Thelen, Zhang, Fink, Lu, Ghosh, Youn, Todd, Mahadevan, Hu and Hu2022, Reference Thelen, Zhang, Fink, Lu, Ghosh, Youn, Todd, Mahadevan, Hu and Hu2023).
In Sacks et al. (Reference Sacks, Schiller and Welch1989), the concept of using GPs as surrogates to predict model outputs at untried locations within the parameter space was introduced. The use of GPs for calibrating computer models was then pioneered by Kennedy and O’Hagan (Reference Kennedy and O’Hagan2000) and (Reference Kennedy and O’Hagan2001). Subsequent work by Bayarri et al. (Reference Bayarri, Berger, Paulo, Sacks, Cafeo, Cavendish, Lin and Tu2007a) and Higdon et al. (Reference Higdon, Kennedy, Cavendish, Cafeo and Ryne2004) refined this framework, with Higdon adopting a fully Bayesian approach. Extending these methodologies to multivariate data (Higdon et al., Reference Higdon, Gattiker, Williams and Rightley2008), and further advancements by others, such as Bayarri et al. (Reference Bayarri, Berger, Cafeo, Garcia-Donato, Liu, Palomo, Parthasarathy, Paulo, Sacks and Walsh2007b), have enhanced the calibration of models with multivariate outputs. (Santner et al., Reference Santner, Williams, Notz and Williams2003), and additional studies, including those by Fang et al. (Reference Fang, Li and Sudjianto2005), Loeppky et al. (Reference Loeppky, Sacks and Welch2009), and Baker et al. (Reference Baker, Barbillon, Fadikar, Gramacy, Herbei, Higdon, Huang, Johnson, Ma, Mondal, Pires, Sacks and Sokolov2022), provide in-depth discussions on the optimal design of computer experiments, focusing on the strategic layout of model evaluations in parameter space. Beyond calibration, GP emulators facilitate understanding variability in model outputs when parameters are uncertain, known as uncertainty analysis. The tutorial by O’Hagan (Reference O’Hagan2006) offers an accessible introduction to uncertainty analysis using GP emulators.
Unlike much of the existing literature on surrogate modeling for complex mechanical systems, which often relies on synthetic data, our study goes beyond the implementation of a surrogate modeling approach and exploits it to enable the calibration of the TUBSsalt constitutive model using
$ 14 $
years of real-world monitoring data collected in an open drift located in the northern main drift of Gorleben, Germany. The primary contribution of this article is to demonstrate the applicability and effectiveness of a GP-based data-driven methodology in addressing a real-world problem. Specifically:
-
• We develop a GP-based surrogate model that approximates the deformation behavior of drift in rock salt formations and verify its accuracy in closely replicating the high-fidelity model’s behavior.
-
• We perform a time-dependent sensitivity analysis using Sobol’ indices, utilizing the surrogate model to extract clear insights from the complex geomechanical constitutive model.
-
• We calibrate the model parameters using real in situ monitoring data from the Gorleben site, benefiting from the efficiency of the GP-based surrogate model.
The remainder of this article is structured as follows: Section 2 introduces the mechanical aspects of modeling a deep geological drift in a rock salt formation. The formulation of the calibration process from experimental and monitoring data as well as the global sensitivity analysis is presented in Section 3. GPs as surrogate models are described in Section 4. Section 5 contains the results of the training and validation as well as calibration with the surrogate model. The article concludes with a summary and outlook in Section 6.
2. Mechanical modeling and numerical solution of a drift in the deep geological formation of rock salt
This section presents the geomechanical model for a drift located on the north-western flank of the Gorleben salt dome in Germany. The latter has been explored with regard to its suitability as a location for a repository for high-level radioactive waste for decades. First, the drift location is described, and the assumptions for the computational model are outlined. Then, the kinematics, governing equations, and the constitutive model TUBSsalt for rock salt are presented. Finally, we outline the numerical solution of the boundary value problem.
2.1. Detailed site description and problem geometry
The mechanical model considered in this work is based on the cross-section of a drift located on the north-western flank of the Gorleben salt dome, a former salt exploration mine in Germany. Monitoring data were collected and provided by the German Federal Company for Radioactive Waste Disposal (BGE mbH).
The considered measurement cross-section is located in the northern main drift of Gorleben with a depth of
$ 840 $
m below the top edge of the ground. Excavation in the area of the measurement location was finished on October 19, 1999, without a recut being carried out afterwards. The monitoring data consists of a time series of convergence measurements, which indicate the change in distance between opposing fixed points inside the rock, from which the deformation rate of the drift contour is derived. At each measuring location, the horizontal and vertical distances are recorded periodically as a standard procedure. Consequently, the computational model adopted in this work represents a similar open emplacement drift of a deep geological repository based on the location where the monitoring data were obtained. Since this phase does not involve storing highly radioactive and heat-generating waste, the problem considered is purely mechanical.
Depending on crystallinity and the presence of secondary aggregates, the rock salt can be divided into homogeneous areas based on their viscous behavior, specifically the steady-state (secondary) creep rate (Hunsche et al., Reference Hunsche, Schulze, Walter and Plischken.d.). A vertical geological cross-section of the Gorleben salt dome and its homogeneous salt areas can be found in Bornemann et al. (Reference Bornemann, Behlau, Fischbeck, Hammer, Jaritz, Keller, Mingerzahn and Schramm2008). Since the measurement location is homogeneously surrounded by a salt formation known as Streifensalz z2HS2, interactions between different homogeneous areas are not taken into account. Therefore, the entire numerical model can be constructed under the assumption of a uniform rock salt material.
Figure 1a depicts the cross-section of the measurement location together with the vertical (2–4) and horizontal (1–3) measurement distances. Figure 1b corresponds to the computational model of the drift in FLAC3D with history locations for the evaluation of the displacements.
Figure 1. Cross-sectional representation of the drift area. a. Cross-section of the measurement locations from the Gorleben site, provided by the BGE mbH. b. A computational model of the drift in FLAC3D, showing the history locations used for evaluating displacements and the mesh used in the numerical solution.
2.2. Continuum mechanics
Modeling in continuum mechanics involves three ingredients: (i) balance equations, (ii) kinematics, and (iii) a constitutive model. This section introduces these three ingredients of continuum mechanics for the deep repository under the assumptions presented in the previous section. We follow the notation for general continuum mechanical equations introduced by Anand and Govindjee (Reference Anand and Govindjee2020), Itasca Consultants GmbH (2023). The constitutive equations for TUBSsalt are presented as they were introduced in Gährken et al. (Reference Gährken, Missal and Stahlmann2015) and Epkenhans et al. (Reference Epkenhans, Wacker and Stahlmann2022).
2.2.1. Balance equations and kinematics
In the current framework, we are focused only on the changes in mechanical quantities due to the excavation of the drift, while keeping the temperature constant, as there is no significant temperature development during the initial phase of the repository operation. Therefore, only the balance of linear momentum needs to be considered, which in its strong form and current configuration is given by
$$ \operatorname{div}\;\boldsymbol{\sigma} +\rho \hskip0.1em \mathbf{b}=\rho \hskip0.1em \frac{\hskip0.1em d\hskip0.1em \mathbf{v}}{\hskip0.1em \mathrm{d}\hskip0.1em t}\hskip0.44em \mathrm{in}\hskip0.24em \Omega, $$
where
$ \Omega $
denotes the computational domain. Further,
$ \boldsymbol{\sigma} $
denotes the Cauchy stress,
$ \rho $
the density,
$ \mathbf{b} $
acceleration caused by external body forces (here: gravitation),
$ \mathbf{v} $
the velocity vector, and
$ \mathrm{d}\hskip0.1em \mathbf{v}/\mathrm{d}\hskip0.1em t $
the acceleration. Dirichlet and Neumann-type boundary conditions are defined at boundaries
$ {\Gamma}_D $
and
$ {\Gamma}_N $
, respectively, as
$$ {\displaystyle \begin{array}{c}\mathbf{v}=\overline{\mathbf{v}}\hskip0.42em \mathrm{on}\hskip0.42em {\Gamma}_D\\ {}\boldsymbol{\sigma} \cdot \mathbf{n}=\mathbf{t}\hskip0.42em \mathrm{on}\hskip0.42em {\Gamma}_N,\end{array}} $$
where
$ \overline{\mathbf{v}} $
and
$ \mathbf{t} $
denote a prescribed velocity and traction and
$ \mathbf{n} $
the surface normal vector. Due to the distinct creep mechanisms of rock salt, a time-dependent problem needs to be considered. The mechanical initial conditions are defined as
$$ {\displaystyle \begin{array}{c}\mathbf{v}\left(t=0\right)=\hskip0.3em {\mathbf{v}}_0\\ {}\boldsymbol{\sigma} \left(t=0\right)=\hskip0.3em {\boldsymbol{\sigma}}_0.\end{array}} $$
The boundary value problem is complemented by the constitutive model TUBSsalt and will be solved using the commercial software FLAC3D (Itasca Consultants GmbH, 2023).
Only small deformations are considered in this work since no long-term analysis of the rock salt emplacement will be performed. Therefore, the kinematics are defined by
$$ \dot{\boldsymbol{\varepsilon}}=\frac{1}{2}\hskip0.1em \left(\operatorname{grad}\hskip0.32em \mathbf{v}+\operatorname{grad}\hskip0.32em {\mathbf{v}}^T\right), $$
where
$ \dot{\boldsymbol{\varepsilon}} $
is the strain rate tensor derived from the velocity vector
$ \mathbf{v} $
.
2.2.2. Constitutive model TUBSsalt
In this subsection, the constitutive model TUBSsalt is briefly introduced. This model is based on the rheological model shown in Figure 2 and is proficient in characterizing various thermo-mechanical aspects of rock salt such as primary, secondary, and tertiary creep, recovery creep, shear, creep and tension failure, healing, and the influence of temperature (Gährken et al., Reference Gährken, Missal and Stahlmann2015; Epkenhans et al., Reference Epkenhans, Wacker and Stahlmann2022).
Figure 2. Rheological model and corresponding strain components of the constitutive model TUBSsalt for rock salt. While
$ {\dot{\boldsymbol{\varepsilon}}}_{el} $
is represented by a spring,
$ {\dot{\boldsymbol{\varepsilon}}}_p $
,
$ {\dot{\boldsymbol{\varepsilon}}}_s $
,
$ {\dot{\boldsymbol{\varepsilon}}}_t $
,
$ {\dot{\boldsymbol{\varepsilon}}}_n $
as well as
$ {\dot{\boldsymbol{\varepsilon}}}_z $
are modelled by hardening or softening sliders and viscous dampers.
Based on the rheological model, the total strain rate
$ \dot{\boldsymbol{\varepsilon}} $
of the TUBSsalt constitutive model is split additively into six components as
$$ \dot{\boldsymbol{\varepsilon}}={\dot{\boldsymbol{\varepsilon}}}_{el}+\underset{{\dot{\boldsymbol{\varepsilon}}}_{vp}}{\underbrace{{\dot{\boldsymbol{\varepsilon}}}_p+{\dot{\boldsymbol{\varepsilon}}}_s+{\dot{\boldsymbol{\varepsilon}}}_t+{\dot{\boldsymbol{\varepsilon}}}_n+{\dot{\boldsymbol{\varepsilon}}}_z}}. $$
Here,
$ {\dot{\boldsymbol{\varepsilon}}}_{el} $
is the elastic strain rate and
$ {\dot{\boldsymbol{\varepsilon}}}_{vp} $
the inelastic, visco-plastic strain rate. The latter consists of the primary creep rate
$ {\dot{\boldsymbol{\varepsilon}}}_p $
, the secondary creep rate
$ {\dot{\boldsymbol{\varepsilon}}}_s $
, the tertiary creep and healing rate
$ {\dot{\boldsymbol{\varepsilon}}}_t $
, the strain rate from creep and shear failure
$ {\dot{\boldsymbol{\varepsilon}}}_n $
, and the strain rate from tension failure
$ {\dot{\boldsymbol{\varepsilon}}}_z $
. In Figure 2, the elastic strain rate
$ {\dot{\boldsymbol{\varepsilon}}}_{el} $
is represented by a spring, and the creep strain rates
$ {\dot{\boldsymbol{\varepsilon}}}_p $
,
$ {\dot{\boldsymbol{\varepsilon}}}_s $
, and
$ {\dot{\boldsymbol{\varepsilon}}}_t $
as well as the failure strain rates
$ {\dot{\boldsymbol{\varepsilon}}}_n $
and
$ {\dot{\boldsymbol{\varepsilon}}}_z $
are described by hardening or softening sliders with yield functions
$ F $
and viscous dampers, each characterised by an individual viscosity
$ \eta $
. For a visco-elastoplastic material model, the stress increment depends on the elastic strain component and can be written as
$$ \dot{\boldsymbol{\sigma}}=\mathbf{D}\left(\boldsymbol{\sigma}, \kappa \right)\cdot \left(\dot{\boldsymbol{\varepsilon}}-{\dot{\boldsymbol{\varepsilon}}}_{vp}\right), $$
where
$ \kappa $
is a loading-history parameter depending on visco-plastic strain rate
$ {\dot{\boldsymbol{\varepsilon}}}_{vp} $
. Rock salt is considered to be an isotropic material, and thus
$ \mathbf{D} $
is the isotropic stiffness tensor, depending on bulk and shear moduli
$ K $
and
$ G $
. Both of these properties decrease as the level of the damage-induced dilatancy
$ {\varepsilon}_{v,d} $
(2.14) increases. In the following, the specific forms of the strain components are introduced in more detail.
The elastic strain rate
$ {\dot{\boldsymbol{\varepsilon}}}_{el} $
is calculated according to Hooke’s law using the stiffness matrix
$ \mathbf{D} $
and the stress rate tensor
$ \dot{\boldsymbol{\sigma}} $
as
$$ {\dot{\boldsymbol{\varepsilon}}}_{el}={\mathbf{D}}^{-1}\cdot \hskip0.3em \dot{\boldsymbol{\sigma}}. $$
Primary and recovery creep only occur after a load change. An increase in the stress deviator causes high primary creep rates, which decrease as the deformation progresses until no more primary deformations occur. There holds
$$ {F}_p>0:{\dot{\boldsymbol{\varepsilon}}}_p=\frac{F_p}{\eta_p^{\ast }}\hskip0.3em \cdot \hskip0.3em \frac{\partial {\sigma}_{eq}}{\partial \boldsymbol{\sigma}} $$
where
$ {\dot{\boldsymbol{\varepsilon}}}_p $
represents the primary creep rate tensor,
$ {F}_p $
the yield function of primary creep,
$ {\eta}_p^{\ast } $
the current viscosity of primary creep and
$ \frac{\partial {\sigma}_{eq}}{\partial \boldsymbol{\sigma}} $
denotes the derivatives of the equivalent stress
$ {\sigma}_{eq} $
with respect to the components of the stress tensor. In case of a decrease of the stress deviator, the viscosity of primary creep
$ {\eta}_p^{\ast } $
is replaced by the viscosity for recovery creep
$ {\eta}_{p, rec} $
.
Secondary creep is always active as soon as a stress deviator is present and leads to a constant creep rate for a constant stress state. It is therefore also referred to as stationary creep and defined as
$$ {\dot{\boldsymbol{\varepsilon}}}_s=\frac{F_s\cdot {q}_s}{\eta_s}\hskip0.3em \cdot \hskip0.3em \frac{\partial {\sigma}_{eq}}{\partial \boldsymbol{\sigma}}, $$
where
$ {\dot{\boldsymbol{\varepsilon}}}_s $
represents the secondary creep rate tensor,
$ {F}_s $
the yield function of secondary creep,
$ {\eta}_s $
is the viscosity parameter for secondary creep, and
$ {q}_s $
is the temperature coefficient for secondary creep.
Tertiary creep initiates as soon as the stress deviator exceeds the dilatancy criteria represented by the yield function
$ {F}_t $
. Above this dilatancy strength, rock salt shows softening which is described by an accelerated reduction of the viscosity
$ {\eta}_t^{\ast } $
depending on the level of damage-induced dilatancy
$ {\varepsilon}_{v,d} $
(2.14):
$$ {F}_t>0:{\dot{\boldsymbol{\varepsilon}}}_t=\frac{F_t\cdot k}{\eta_t^{\ast}\cdot {q}_t}\hskip0.3em \cdot \hskip0.3em \frac{\partial {Q}_t}{\partial \boldsymbol{\sigma}}, $$
where
$ {\dot{\boldsymbol{\varepsilon}}}_t $
represents the tertiary creep rate tensor,
$ {\eta}_t^{\ast } $
is the current viscosity of tertiary creep,
$ k $
a coefficient for loading rate and stress state,
$ {q}_t $
a temperature coefficient, and
$ \frac{\partial {Q}_t}{\partial \boldsymbol{\sigma}} $
the directional derivatives of the potential function for tertiary creep
$ {Q}_t $
with respect to the components of the stress tensor. The increase in volume or dilatancy due to microcracks is taken into account in
$ \frac{\partial {Q}_t}{\partial \boldsymbol{\sigma}} $
.
The process of healing replaces the tertiary creep once a stress state falls below the dilatancy threshold and damage has already been determined. There holds
$$ {F}_t>-{\sigma}_z:{\dot{\boldsymbol{\varepsilon}}}_t=\frac{F_t\cdot {q}_v}{\eta_v^{\ast }}\hskip0.3em \cdot \hskip0.3em \frac{\partial {Q}_v}{\partial \boldsymbol{\sigma}}, $$
where
$ {\sigma}_z $
is the tensile strength,
$ {q}_v $
a temperature coefficient,
$ {\eta}_v^{\ast } $
the current viscosity of healing, and
$ \frac{\partial {Q}_v}{\partial \boldsymbol{\sigma}} $
the directional derivatives of the potential function for healing
$ {Q}_v $
with respect to the components of the stress tensor.
Creep and shear failure: Once the damage-induced dilatancy
$ {\varepsilon}_{v,d} $
(2.14) exceeds the failure volumetric strain
$ {\varepsilon}_{v,d,b,\ast } $
, time-dependent failure deformations occur in addition to the deformations from the creep components:
$$ {\varepsilon}_{v,d}\ge {\varepsilon}_{v,d,b\ast }:{\dot{\boldsymbol{\varepsilon}}}_n=\frac{F_n\cdot k}{\eta_n^{\ast}\cdot {q}_n}\hskip0.3em \cdot \hskip0.3em \frac{\partial {\sigma}_{eq}}{\partial \boldsymbol{\sigma}}, $$
where
$ {\dot{\boldsymbol{\varepsilon}}}_n $
is the strain rate tensor of the post-failure,
$ {\eta}_n^{\ast } $
the current viscosity parameter for healing, and
$ {q}_n $
a temperature coefficient.
Tension failure occurs when tensile material strength
$ {\sigma}_{z,0} $
is exceeded and is denoted by the strain rate tensor
$ {\dot{\boldsymbol{\varepsilon}}}_z $
.
An important characteristic of the constitutive model is the damage-induced dilatancy
$ {\varepsilon}_{v,d} $
, which is decisive for the evaluation of the geological integrity of the excavation damage zone of rock salt. As soon as the dilatancy strength of the material is exceeded, which is achieved when the yield function
$ {F}_t $
becomes greater than
$ 0 $
, softening takes place as a result of crack formation, which in turn creates pathways for radionuclides. In the first step, the tensor of the damage-induced strain increment
$ \Delta {\boldsymbol{\varepsilon}}_d $
can be determined using a timestep
$ \Delta t $
as
$$ \Delta {\boldsymbol{\varepsilon}}_d=\left({\dot{\boldsymbol{\varepsilon}}}_t+{\dot{\boldsymbol{\varepsilon}}}_n\right)\cdot \Delta t. $$
The increment of the damage-induced dilatancy
$ {\varepsilon}_{v,d} $
is equivalent to the first invariant of
$ \Delta {\boldsymbol{\varepsilon}}_d $
:
$$ \Delta {\varepsilon}_{v,d}={I}_1\left(\Delta {\boldsymbol{\varepsilon}}_d\right)=\Delta {\varepsilon}_{d, xx}+\Delta {\varepsilon}_{d, yy}+\Delta {\varepsilon}_{d, zz}. $$
For further details on quantities not explained in this section, such as the yield functions
$ {F}_p $
,
$ {F}_t $
, and
$ {F}_n $
, the potential functions for tertiary creep and healing
$ {Q}_t $
and
$ {Q}_v $
, as well as derivations of all equations of TUBSsalt provided here, the reader is referred to Epkenhans et al. (Reference Epkenhans, Wacker and Stahlmann2022). In summary, the constitutive model TUBSsalt comprises a total of
$ 25 $
material parameters, summarized in Appendix B, Table B1.
2.3. Numerical solution
The constitutive model TUBSsalt is implemented into the commercial software FLAC3D (Itasca Consultants GmbH, 2023), which is a program for three-dimensional engineering mechanics computations. FLAC3D relies on the finite difference method (FDM) and translates the continuum equations into ordinary differential equations for each element, which are then solved using an explicit central finite difference approach in time. The spatial discretization in FLAC3D is realized by meshing the continuum into hexahedral elements, with the option to use tetrahedral, wedge, and pyramid elements. It is recommended to use hexahedral elements, which are divided into tetrahedral elements to reduce the volumetric locking effect and thus achieve more accurate results. The balance of the linear momentum, Equation (2.1), is iterated to an equilibrium state, which is achieved when the unbalanced mechanical force for all the grid points in the model is negligibly low.
Time-dependent phenomena such as creep require a timestep
$ \Delta t $
to solve the equations of the TUBSsalt constitutive model. At the same time, for creep analysis, the state of equilibrium must be maintained, otherwise inertial effects may affect the solution. For this purpose, the unbalanced force is monitored in the model. The finite volume scheme can be summarised for each time step as follows: After new strain rates are determined from nodal velocities, new stresses are calculated from strain rates and previous stresses using constitutive equations. By applying the balance of the linear momentum, new velocities, and deformations are subsequently calculated from stresses and forces. More information on the FLAC3D solution algorithm can be found in Itasca Consultants GmbH (2023).
3. Model calibration and sensitivity analysis
In this section, we briefly state the calibration process of the mechanical model introduced in Section 2. We distinguish between two model calibration stages. The first stage assumes that stress–strain data are available (Section 3.1) and the constitutive model can be directly calibrated from common geomechanical laboratory tests. This approach is usually employed as the first calibration step to state some prior knowledge regarding the material parameters. The second stage uses in situ monitoring data of the drift convergence (Section 3.2). This approach is used for model parameter recalibration and considers the real conditions of the deep repository. As part of the calibration process, a method for global sensitivity analysis is introduced in Section 3.3.
3.1. Initial model calibration from mechanical testing data
Commonly, laboratory tests are used for the calibration of constitutive model parameters. In the case of rock salt and as shown in Table B1, short-term triaxial compression or extension tests, long-term creep tests, and healing tests as well as indirect tensile tests need to be conducted under variation of temperature and confining stress to be able to cover all strain parts described in Section 2.2.2. While strength tests use a constant axial strain rate to apply a load on the mostly cylindrical specimens with a fixed confining pressure, the specimens in creep and healing tests have a predefined stress state so that the time-dependent behavior can be observed.
Common to the aforementioned mechanical tests is that strain states are considered, which yield well-defined stress states. Hence, stress–strain data can be obtained from the mechanical tests, and the calibration of the constitutive model is a regression problem that can be cast as the optimization problem
$$ {\boldsymbol{\theta}}^{\ast }={\mathrm{arg}\ \mathrm{min}}_{\boldsymbol{\theta}}\sum \limits_{i=1}^N\parallel {\boldsymbol{\sigma}}^i-\boldsymbol{\sigma} \left({\boldsymbol{\varepsilon}}^i;\boldsymbol{\theta} \right){\parallel}^2, $$
with
$ N $
the number of stress–strain data pairs. Here, we consider the minimization of the least-squares error between measured stress
$ {\boldsymbol{\sigma}}^i $
and predicted stress for the associated strain value
$ {\boldsymbol{\varepsilon}}^i $
. The constitutive model is parametrized in the material parameters
$ \boldsymbol{\theta} $
, and the semicolon denotes parametrization. This approach was used to determine and calibrate all TUBSsalt parameters for Gorleben salt, resulting in a parameter set presented in Stahlmann et al. (Reference Stahlmann, Missal and Gährken2016).
However, the samples used in laboratory tests are no longer in their original state and the material tests only give insight into the material behavior in localized regions from where the sample has been extracted. As a consequence, model predictions typically do not match with real-world monitoring data, necessitating re-calibration of the constitutive model. Therefore, model calibration based on monitoring data is addressed in the following section.
3.2. Inverse problem formulation for the model re-calibration
We now address the task of model calibration from monitoring data. The objective is to identify the optimal set of material model parameters
$ \boldsymbol{\theta} $
that minimize the discrepancy (least-squares error) between model predictions and the in-situ monitoring data. This process is framed as an inverse problem, solved with an optimization method.
Let
$ {\mathbf{Y}}_{\mathrm{moni}}^T=\left[{\mathbf{y}}_{\mathrm{moni}}^T\left({t}_1\right),{\mathbf{y}}_{\mathrm{moni}}^T\left({t}_2\right),\dots, {\mathbf{y}}_{\mathrm{moni}}^T\left({t}_n\right)\right] $
denote the vector of experimental observations at different time instances for a specific location with
$ n $
the total number of time instances. The corresponding model prediction is denoted by
$ {\mathbf{s}}_i\left(\boldsymbol{\theta} \right)=\mathbf{s}\left({t}_i,\boldsymbol{\theta} \right) $
. Then, in analogy to Equation (3.1)), material parameters are identified from
$$ {\boldsymbol{\theta}}^{\ast }={\mathrm{arg}\ \mathrm{min}}_{\boldsymbol{\theta}}\;\sum \limits_{i=1}^n\parallel {\mathbf{y}}_{\mathrm{moni}}\left({t}_i\right)-{\mathbf{s}}_i\left(\boldsymbol{\theta} \right){\parallel}^2. $$
Note that in Equations (3.1)) and (3.2) we omit a weighting matrix, which is often introduced to account for the size of measurement errors. In our setting, the error sizes did not have a relevant influence. The optimization problem (3.2) is implemented using the SciPy Python library (Virtanen et al., Reference Virtanen, Gommers, Oliphant, Haberland, Reddy, Cournapeau, Burovski, Peterson, Weckesser and Bright2020), employing a global optimization algorithm. We opted for the differential evolution algorithm (Storn and Price, Reference Storn and Price1997), but also experimented with genetic and dual annealing algorithms. However, differential evolution proved to be faster and more reliable in this context.
3.3. Sensitivity analysis based on time-dependent Sobol’ indices
Mechanical models that are used in geomechanical contexts are often characterized by a large number of material parameters, and so is TUBSsalt. However, not every parameter in a data set may be sensitive to specific monitoring data. To focus only on the important ones and to reduce the number of model evaluations when solving (3.2), we perform a sensitivity analysis.
In this contribution, we compute sensitivities based on Sobol’ indices (Sobol, Reference Sobol2001; Saltelli, Reference Saltelli2002; Saltelli et al., Reference Saltelli, Annoni, Azzini, Campolongo, Ratto and Tarantola2010). The first-order Sobol’ index quantifies the contribution of each parameter
$ {\theta}_j $
to the variance in the predicted solution
$ {\mathbf{s}}_i\left(\boldsymbol{\theta} \right) $
, while averaging out the effects of other inputs. The first-order Sobol’ index for each time point
$ {t}_i $
is defined as:
$$ S{1}_i\left({\theta}_j\right)=\frac{V_{\theta_j}\left({E}_{\theta_{\sim j}}\left({\mathbf{s}}_i\left(\boldsymbol{\theta} \right)|{\theta}_j\right)\right)}{V\left({\mathbf{s}}_i\left(\boldsymbol{\theta} \right)\right)}, $$
where
$ V\left({\mathbf{s}}_i\left(\boldsymbol{\theta} \right)\right) $
is the variance of the solution at time
$ {t}_i $
, and
$ {E}_{\theta_{\sim j}}\left({\mathbf{s}}_i\left(\boldsymbol{\theta} \right)|{\theta}_j\right) $
is the conditional expectation given
$ {\theta}_j $
.
Similarly, the total-order Sobol’ index for each solution
$ {\mathbf{s}}_i\left(\boldsymbol{\theta} \right)=\mathbf{s}\left({t}_i,\boldsymbol{\theta} \right) $
and input parameter
$ {\theta}_j $
at time
$ {t}_i $
measures the total effect of
$ {\theta}_j $
on the variance of the solution, including interactions with other inputs. The total-order Sobol’ index is defined as:
$$ {ST}_i\left({\theta}_j\right)=1-\frac{V_{\theta_{\sim j}}\left({E}_{\theta_j}\left({\mathbf{s}}_i\left(\boldsymbol{\theta} \right)|{\theta}_{\sim j}\right)\right)}{V\left({\mathbf{s}}_i\left(\boldsymbol{\theta} \right)\right)}, $$
where
$ {E}_{\theta_j}\left({\mathbf{s}}_i\left(\boldsymbol{\theta} \right)|{\theta}_{\sim j}\right) $
is the conditional expectation given all parameters except
$ {\theta}_j $
.
Sobol’ indices
$ S{1}_i\left({\theta}_j\right) $
and
$ {ST}_i\left({\theta}_j\right) $
allow us to quantify the sensitivity of the solution
$ \mathbf{s}\left({t}_i,\boldsymbol{\theta} \right) $
at each time step
$ {t}_i $
to the different parameters
$ \boldsymbol{\theta} $
. This analysis helps us understand how each parameter affects the system’s behavior over time. To identify and potentially discard less influential parameters, we propose calculating global indicators: cumulated, time-averaged, and maximum Sobol’ indices. These indicators offer insights into the overall influence of the parameters throughout the time-dependent solution and are defined as follows:
-
• The Normalized cumulated Sobol’ indices are calculated by summing the influence of each input parameter over all time points and then normalizing by the maximum integrated value across all input parameters:
$$ \mathrm{Norm}.\int S{1}_i=\frac{\sum_{i=1}^nS{1}_i}{\max \left({\sum}_{i=1}^nS{1}_i,{\sum}_{i=1}^n{ST}_i\right)},\hskip1em \mathrm{Norm}.\int {ST}_i=\frac{\sum_{i=1}^n{ST}_i}{\max \left({\sum}_{i=1}^nS{1}_i,{\sum}_{i=1}^n{ST}_i\right)}. $$
-
• The time-averaged Sobol’ indices provide an overall measure of the importance of each input parameter across all time points. They are computed as the mean of the Sobol’ indices at each time step
$ {t}_i $
:
$$ \overline{S1}=\frac{1}{n}\sum \limits_{i=1}^n\;S{1}_i,\hskip1em \overline{ST}=\frac{1}{n}\sum \limits_{i=1}^n\;{ST}_i. $$
-
• The maximum Sobol’ indices identify the time steps at which each input parameter has the highest influence. These are particularly useful for pinpointing critical moments in the time-dependent behavior of the model:
$$ \max S{1}_i=\underset{t}{\max}\left(S{1}_i\right),\hskip1em \max {ST}_i=\underset{t}{\max}\left({ST}_i\right). $$
A general framework for time-dependent variance-based sensitivity analysis has been put forth in Alexanderian et al. (Reference Alexanderian, Gremaud and Smith2020). Therein, generalized Sobol’ indices for processes are defined via a decomposition of the covariance function of the process. If they are approximated with an unweighted quadrature on a uniform grid, a normalized version of the cumulated Sobol’ indices introduced above is recovered. An even more general framework, presented in Gamboa et al. (Reference Gamboa, Janon, Klein and Lagnoux2014), considers sensitivity analysis of functional outputs, covering time-dependent responses as a special case.
The convergence of the Sobol’ indices estimation can be assessed using the maximum relative change between successive iterations. The relative change for each Sobol’ index is calculated as the absolute difference between the current and previous indices, normalized by the previous indices plus a small constant to avoid division by zero. Specifically, the relative change for
$ S{1}_i $
is given by:
$$ \mathrm{S}{1}_i\;\mathrm{relative}\ \mathrm{change}=\frac{\left|\mathrm{S}{1}_{i\mathrm{curr}}-\mathrm{S}{1}_{i\mathrm{prev}}\right|}{\mathrm{S}{1}_{i\mathrm{prev}}+1\times {10}^{-10}}, $$
and the maximum value across all indices is monitored. The same applies to
$ {ST}_i $
. If the maximum change in both indices is below a predefined threshold (
$ 0.01 $
), the process is deemed converged, indicating that the sampling effort is sufficient for accurate sensitivity analysis.
A bottleneck is that a sensitivity analysis using Sobol’ indices with Saltelli’s algorithm (Saltelli et al., Reference Saltelli, Annoni, Azzini, Campolongo, Ratto and Tarantola2010), pick-and-freeze estimators (Gamboa et al., Reference Gamboa, Janon, Klein, Lagnoux and Prieur2016), or other methods requires a significant number of model evaluations. One way to bind the associated computational cost is by using surrogate models such as the Polynomial Chaos expansion (Sudret, Reference Sudret2008) or GPs. In this contribution, GPs are used and presented in the next section.
4. Gaussian processes as surrogate models for complex computer models
The model
$ {\mathbf{s}}_i\left(\boldsymbol{\theta} \right) $
represents the mechanical behavior of the deep geological repository introduced in Section 2. Evaluating
$ {\mathbf{s}}_i\left(\boldsymbol{\theta} \right) $
directly is computationally demanding, especially in many-query scenarios such as model calibration, optimization, or design. Surrogate modeling offers an efficient alternative to address the computational challenges associated with directly solving
$ {\mathbf{s}}_i\left(\boldsymbol{\theta} \right) $
. Among the available techniques, Proper Orthogonal Decomposition (POD) (Agarwal et al., Reference Agarwal, Urrea-Quintero, Wessels and Wick2024) and Physics-Informed Neural Networks (PINNs) (Anton et al., Reference Anton, Tröger, Wessels, Römer, Henkes and Hartmann2024) are particularly suitable for full-field simulations where capturing the entire spatial domain is crucial. However, GPs present a compelling option when the focus is on specific quantities of interest, such as an excavation convergence as in the present study, rather than the entire field. GPs are especially advantageous due to their ability to provide probabilistic predictions and handle smaller datasets effectively. In this study, we adopt GPs to develop a surrogate model, denoted as
$ {\hat{\mathbf{s}}}_i\left(\boldsymbol{\theta} \right) $
, to approximate the complex behavior of
$ {\mathbf{s}}_i\left(\boldsymbol{\theta} \right) $
.
In scenarios where
$ \mathbf{s} $
predicts a time-sequence output, we define a set of GPs,
$ {\left\{{\hat{\mathbf{s}}}_i\right\}}_{i=1}^n $
. Each
$ {\hat{\mathbf{s}}}_i $
corresponds to a distinct time instance
$ {t}_i $
, as:
$$ {\hat{\mathbf{s}}}_i\left(\boldsymbol{\theta} \right)\sim \mathcal{GP}\left({m}_i\left(\boldsymbol{\theta} \right),\hskip0.34em {k}_i\Big(\boldsymbol{\theta}, {\boldsymbol{\theta}}^{\prime}\Big)\right), $$
where
$ {m}_i\left(\boldsymbol{\theta} \right) $
is the mean function and
$ {k}_i\left(\boldsymbol{\theta}, {\boldsymbol{\theta}}^{\prime}\right) $
is the covariance (or kernel) function for the
$ i $
-th GP. This function quantifies the similarity between two sets of parameter configurations,
$ \boldsymbol{\theta} $
and
$ {\boldsymbol{\theta}}^{\prime } $
, for the specific time point. Each GP
$ {\hat{\mathbf{s}}}_i $
is independently trained on a subset of the data corresponding to its time point:
$$ \left[{\hat{\mathbf{s}}}_i\left({\boldsymbol{\theta}}_1\right),\hskip0.34em {\hat{\mathbf{s}}}_i\left({\boldsymbol{\theta}}_2\right),\dots, {\hat{\mathbf{s}}}_i\left({\boldsymbol{\theta}}_N\right)\right]\sim \mathcal{N}\left(\mathbf{0},{\mathbf{K}}_i\right), $$
where
$ {\mathbf{K}}_i $
is the covariance matrix for the
$ i $
th GP, computed using a kernel
$ {k}_i $
. Note that we could also introduce a priori correlation between GPs at different points in time; however, the independence assumption is very common. Here, for each GP, we use the Matérn kernel (with smoothness parameter
$ \nu =5/2 $
) defined as:
$$ {k}_i\left(\boldsymbol{\theta}, {\boldsymbol{\theta}}^{\prime}\right)={\sigma}_i^2\left(1+\sqrt{5}r+\frac{5}{3}{r}^2\right)\exp \left(-\sqrt{5}r\right), $$
where
$ r\left(\boldsymbol{\theta}, {\boldsymbol{\theta}}^{\prime}\right)=\sqrt{\sum_{d=1}^D{\left(\frac{\theta_d-{\theta}_d^{\prime }}{l_{d,i}}\right)}^2} $
, with
$ {l}_{d,i} $
and
$ {\sigma}_i^2 $
being the length-scales and variance parameter for the ith GP.
The next crucial step is the determination of the GP-based surrogate model hyperparameters. These hyperparameters include the length-scales
$ {l}_{d,i} $
, variance
$ {\sigma}_i^2 $
, and any other parameters specific to the chosen kernel function. The hyperparameters are typically optimized by maximizing the likelihood of the observed data under the GP model. This optimization can be formulated as:
$$ \underset{\sigma_i,{l}_{d,i}}{\max}\;\log\;\mathrm{\mathcal{L}}\left({\sigma}_i,{l}_{d,i}|{D}_i\right), $$
where
$ \mathrm{\mathcal{L}} $
is the likelihood of the training data
$ {D}_i $
given the hyperparameters. This process is often carried out using gradient-based optimization techniques.
Upon training each
$ {\hat{\mathbf{s}}}_i $
with a dataset
$ {D}_i=\left\{\left({\boldsymbol{\theta}}_1,{s}_{i1}\right),\dots, \left({\boldsymbol{\theta}}_N,{s}_{iN}\right)\right\} $
, where
$ {s}_{ij} $
is the output from
$ \mathbf{s} $
for input
$ {\boldsymbol{\theta}}_j $
at time
$ {t}_i $
, each GP can make predictions for unseen parameter sets
$ {\boldsymbol{\theta}}_{\ast } $
at its respective time instance. This ensemble of GPs allows for efficient evaluation of the time-sequenced outputs of
$ \mathbf{s}\left({t}_i,\boldsymbol{\theta} \right) $
without relying on computationally expensive numerical solutions.
The mean prediction
$ {\mu}_i\left({\boldsymbol{\theta}}_{\ast}\right) $
and the standard deviation
$ {\sigma}_i\left({\boldsymbol{\theta}}_{\ast}\right) $
at each time point
$ {t}_i $
are given by:
$$ {\mu}_i\left({\boldsymbol{\theta}}_{\ast}\right)={\mathbf{K}}_{\ast}^T{\mathbf{K}}_i^{-1}\mathbf{s},\hskip1em {\sigma}_i^2\left({\boldsymbol{\theta}}_{\ast}\right)={k}_i\left({\boldsymbol{\theta}}_{\ast },{\boldsymbol{\theta}}_{\ast}\right)-{\mathbf{K}}_{\ast}^T{\mathbf{K}}_i^{-1}{\mathbf{K}}_{\ast }, $$
where
$ {\mathbf{K}}_{\ast } $
is the covariance vector between the training inputs and
$ {\boldsymbol{\theta}}_{\ast } $
, and
$ {\mathbf{K}}_i $
is the covariance matrix of the training inputs.
These predictions not only provide the expected output of the model
$ \mathbf{s} $
for an unseen parameter set but also quantify the uncertainty associated with these predictions. This feature of GPs is particularly useful in decision-making processes where uncertainty plays a crucial role.
It is important to note that before training the GP-based surrogate model, we perform feature scaling on both input parameters and outputs to achieve zero mean and unit variance. This step is crucial for enhancing the surrogate model’s numerical stability, which is developed using the Scikit-learn library (Pedregosa et al., Reference Pedregosa, Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos, Cournapeau, Brucher, Perrot and Duchesnay2011). This library offers a ready-to-use module to build surrogate models based on GPs for multivariate nonlinear regression problems. GP hyperparameters are optimized using the Limited-memory Broyden–Fletcher–Goldfarb–Shanno with bound constraints (L-BFGS-B) algorithm; see Scikit-learn library documentation for more details. This optimization process is key to maximizing the log-marginal likelihood defined in (4.4), ensuring that the kernel parameters are finely tuned to best represent the underlying patterns of the data.
5. Numerical analysis
As described in Section 2, TUBSalt is a highly nonlinear time-dependent constitutive material model with
$ 25 $
material parameters. After a preliminary study, we train a GP-based surrogate model (Section 4) to capture the input–output relationship between selected parameters of the constitutive model TUBSsalt and the drift convergence at different time instances, conduct a sensitivity analysis and calibrate the constitutive model. This section is structured as follows: First, the high-fidelity model evaluation and dataset creation for training the GP-based surrogate model is presented in Section 5.1. Section 5.2 deals with the training of the GP-based surrogate model, the accuracy evaluation with testing data, and sensitivity analysis using Sobol’ indices. Finally, in Section 5.3, the model parameters are calibrated using in situ monitoring data.
5.1. High-fidelity model simulation
Monitoring data: The mechanical model described in Section 2 represents an open drift located in the northern main drift of Gorleben for which monitoring data is available. The first convergence measurement took place on October 20, 1999, 1 day after the excavation of the drift. This is important because it allows for the measurement of primary creep, which is most pronounced at the beginning of the excavation. On the day of the initial measurement, the horizontal distance (1–3) measured
$ 8.94 $
m, and the vertical distance (2–4) was
$ 6.00 $
m, as depicted in Figure 1a. The duration of the monitoring was approximately
$ 14.3 $
years, with a significant variation of time intervals between measurements. During the first 2 months, measurements were taken every 2–3 days to capture the high creep rates. Subsequently, the measurement interval was increased from 1 week to 1 month, and from the end of 2002 onwards, measurements were taken approximately every 6 months as a constant creep rate had been achieved. Since the measuring points of the convergence distances are anchored in the rock over an anchoring length, the deformations of the cavity are also measured at these fixed points in the rock salt. Those locations are defined as history points in the numerical model in FLAC3D, demonstrated in Figure 1b. The convergences
$ {\mathbf{y}}_{\mathrm{moni}}^T\left({t}_i\right)={\left[\Delta {u}_x\left({t}_i\right),\hskip0.34em \Delta {u}_z\left({t}_i\right)\right]}^T $
are measured as a vertical and horizontal distance between two points, therefore the displacements
$ {u}_x $
in horizontal and
$ {u}_z $
in vertical direction are obtained from the displacement vector
$ \mathbf{u} $
for the specific nodes. The uncertainty of the convergence measurements is around
$ \pm 0.5 $
mm, while observations are on the order
$ \mathcal{O}\left({10}^2\right) $
. Consequently, we consider the monitoring data to be noise-free.
Model setup: Since only one material is considered and the cross-section of the drift is almost symmetric, only the right symmetry half of the model needs to be discretized. Dirichlet boundary conditions are applied by fixing the left, right, and bottom sides of the system in their normal directions. We model the overburden with a load of
$ 16.774 $
MPa applied to the top of the model, which is a Neumann boundary condition. According to the geology of the site given in Bornemann et al. (Reference Bornemann, Behlau, Fischbeck, Hammer, Jaritz, Keller, Mingerzahn and Schramm2008), the load results from the different heights and densities of the rock salt, cap rock, tertiary and quaternary layer. The densities were obtained from Kock et al. (Reference Kock, Eickemeier, Frieling, Heusermann, Knauth, Minkley, Navarro, Nipp and Vogel2012). The initial stress state is derived from the overburden load and displacements for time
$ t=0 $
are set to zero. The initial stress state is assumed to be isotropic. The slight directional dependence of the monitoring data can be attributed to the differing dimensions of the drift cross-section in the
$ x $
and
$ z $
directions. The absence of nearby geotechnical structures and the sufficient distance between the measurement site and other side drifts, or homogeneous areas allow the model to be simplified to a 2D model. The dimensions in the
$ x $
and
$ z $
direction for the model are given by
$ {l}_x=50 $
m and
$ {l}_z=100 $
m, so that an influence of the drift on the system boundaries can be excluded. Temperature measurements conducted in the close area of the drift show an almost constant value of
$ 37{}^{\circ} $
C. Therefore, a constant temperature is assumed over the entire system as a simplification. The mesh density of the model was optimized by performing a convergence study, ensuring that increasing the number of elements results in no significant changes in the curves of simulated horizontal and vertical convergences. Before the time-dependent creep analysis, the model reaches an equilibrium state twice: first, after applying the initial and boundary conditions, and second, after setting the stresses in the zones inside the cross section to zero to simulate the excavation. The equilibrium state is reached as soon as the average ratio between out-of-balance and total forces falls below
$ 1\cdot {10}^{-6} $
.
Preliminary parameter selection: TUBSalt has a total of
$ 25 $
material parameters. Out of these,
$ 5 $
parameters can be read directly from laboratory stress–strain data, namely
$ {K}_0 $
,
$ {G}_0 $
,
$ {\varepsilon}_{v,d,b} $
,
$ {\sigma}_{z,0} $
and
$ \rho $
. In the following, these are considered as well-calibrated and independent of the monitoring data. The remaining
$ 20 $
parameters serve as potential parameters for calibration. We then conducted a preliminary, empirical study with only two simulations: several damage-associated strain components were deactivated, and the result was compared with a reference solution obtained with the full set of strain components (2.5). It has been found that the test case considered is dominated by primary and secondary creep mechanisms and damage-associated strain components are negligible. Details can be found in Appendix A. As a result, in the following we focus on a total of
$ 7 $
TUBSsalt parameters for primary and secondary creep:
$ {\eta}_p $
,
$ {\sigma}_{0, eq,p} $
,
$ {p}_p $
, and
$ {E}_p $
for primary creep and
$ {\eta}_s $
,
$ {\sigma}_{0, eq,s} $
, and
$ {p}_s $
for secondary creep, and thus
$ \boldsymbol{\theta} =\left[{\eta}_p,{\sigma}_{0, eq,p},{p}_p,{E}_p,{\eta}_s,{\sigma}_{0, eq,s},{p}_s\right] $
. The parameter ranges of those
$ 7 $
parameters, as well as the parameters that are considered as fixed, are summarized in Appendix B, Table B1.
Sampling of parameter values in FLAC3D: Following the approach in Section 3.1, a parameter set for the constitutive model TUBSsalt based on laboratory strength and creep tests has been determined in a project presented in Stahlmann et al. (Reference Stahlmann, Missal and Gährken2016) for the Gorleben site. These values serve as characteristic parameter values for Gorleben salt and also help define reasonable parameter ranges for those parameters that need calibration, see also Table B1 in Appendix B. To create training and testing data by means of a multitude of simulations, the parameter values of the selected TUBSsalt parameters are sampled in the given ranges in Table B1 assuming a uniform distribution. For this purpose, a SciPy Quasi-Monte Carlo generator is used to generate a Sobol sequence by creating low-discrepancy, quasi-random numbers (Sobol, Reference Sobol1967; Owen, Reference Owen2019). The loop for parameter sampling is implemented in the creep analysis of Flac3D, so that the saved equilibrium state of the model is recalled for every iteration. The creep analysis is performed for the same duration of the monitoring. To maintain the state of equilibrium, the average ratio between out-of-balance and total forces is limited to
$ 5\cdot {10}^{-6} $
. The timestep increases over time with a maximum of
$ \Delta t=5.56 $
h. With the given configuration, a single simulation run of the high-fidelity model requires approximately
$ 4 $
min.
Processing the full-field simulation results to convergences: After carrying out the simulations to create training and testing data, the displacements derived at the history locations are converted into convergences at the time instances for which monitoring data are available. To determine the vertical convergences, the amounts of the vertical displacements
$ {u}_z\left({t}_i\right) $
at the top and bottom are summed up together. To determine the horizontal convergences, the amount of the horizontal displacements
$ {u}_x\left({t}_i\right) $
is multiplied by two to account for the symmetry of the system. The displacements of an exemplary full-field simulation at
$ t=14.5 $
year is depicted in Figure 3. To compare the simulated convergences with the monitoring data, the convergences of the simulations must be set to zero at the time of the initial monitoring measurement. Therefore, the convergences from the first day after excavation of the drift can not be analysed and will be neglected. Convergences are determined at all points in time at which monitoring is available. As a result, Figure 4 shows the processed data of
$ 200 $
simulations. Figure 4a depicts the histograms obtained after sampling the input material parameters, while Figure 4b,c correspond to the vertical and horizontal convergences over time (blue lines) in comparison with the corresponding monitoring data (black squares) and the obtained probability density function (PDF) at selected time instances (red area). This figure illustrates how the uncertainty in material parameters propagates through the mechanical model implemented in FLAC3D, reflecting its nonlinearity. That is, uniformly sampled inputs result in heavily tailed outputs. The task of the GP-based surrogate is to represent this input–output relationship by learning the PDF of the convergences at different time instances as a function of the model parameters.
Figure 3. Full-field simulation at
$ t=14.5 $
yr corresponding to the geometry depicted in Figure 1b. a. Vertical
$ {u}_z $
and b. horizontal displacements
$ {u}_x $
. History locations are marked by gray dots. The simulation uses TUBSsalt parameters:
$ {\eta}_p=8.0\cdot {10}^4 $
MPa
$ \cdot $
d,
$ {E}_p=75 $
MPa,
$ {\sigma}_{0, eq,p}=30 $
MPa,
$ {p}_p=0.5 $
,
$ {\eta}_s=3.0\cdot {10}^7 $
MPa
$ \cdot $
d,
$ {\sigma}_{0, eq,s}=30 $
MPa, and
$ {p}_s=1.5 $
.
Figure 4. Input–output simulation data at the monitoring location. a. Histograms for the input material parameters
$ {\eta}_p $
,
$ {E}_p $
,
$ {\eta}_s $
,
$ {p}_s $
,
$ {\sigma}_{0, eq,p} $
,
$ {p}_p $
, and
$ {\sigma}_{0, eq,s} $
of TUBSsalt. b. Vertical and c. horizontal convergence trajectories over time. b. and c. show the monitoring data (black squares) and
$ 200 $
model realizations (blue lines) from the FLAC3D simulations along with the obtained PDF at selected time instances (red area).
To understand the individual and combined effect of the material parameters on the drift convergence, Sobol’ indices-based global sensitivity analysis can be performed as outlined in Section 3.3.
5.2. GP-based surrogate and sensitivity analysis
In this subsection, we focus on developing a surrogate model for predicting the temporal evolution of horizontal and vertical drift convergences based on the synthetic dataset described in the previous Section 5.1. The GP-based surrogate is constructed so that it predicts the convergences at those time instances
$ {t}_i $
for which monitoring data are available. As shown in Figure 4, the original density of the measurement data is significantly higher at the beginning of the monitoring. To prevent distortion of both the time-dependent Sobol’ indices and the subsequent calibration process, the available data points in the first
$ 3.6 $
years are filtered based on the time intervals in such a way that the measurement points are distributed as evenly as possible over time. As a result, our approach treats horizontal and vertical convergences at
$ 40 $
equidistant time instances as independent outputs, framing our problem a multivariate nonlinear regression with
$ 7 $
inputs (parameters of the TUBSsalt material model for rock salt) and
$ 80 $
outputs (the drift convergence at
$ 40 $
different selected time instances). This finally leads to an ensemble of
$ 80 $
GPs,
$ 40 $
each for both vertical and horizontal convergences. Each GP takes the
$ 7 $
material parameters as input and outputs the value of the convergence at its specific time instance as output.
GP-based surrogate model: We split the dataset into training (
$ 75\% $
) and testing (
$ 25\% $
) datasets, which results in
$ 150 $
model realizations for training and
$ 50 $
for testing the accuracy of the surrogate model. The selected amount of training data results from an investigation of the relation between the amount of training data and the accuracy of the surrogate model prediction. On the one hand, we experienced that fewer model realizations led to a significant decrease in the
$ {R}^2 $
coefficient for the testing dataset. On the other hand, if more model realizations were used, only a marginal increase in the
$ {R}^2 $
coefficient would be achieved compared with the computational cost required to train the surrogate model. Note that the number of model realizations can be further reduced and optimized using adaptive sampling strategies, see Fuhg et al. (Reference Fuhg, Fau and Nackenhorst2021) for a recent review. The generation of the GP-based surrogate model with
$ 7 $
parameters and
$ 150 $
training data samples takes
$ 12.3 $
s.
Figure 5 displays scatter plots comparing the true output versus predicted output for the GP-based surrogate model trained to predict vertical and horizontal convergence at different time instances. Each subplot represents a distinct time instance, in years, as indicated in their titles, respectively. The blue dots indicate training data, and the red dots indicate testing data. The black diagonal line represents perfect predictions. The high
$ {R}^2 $
scores for the training data across all components indicate that the model fits the training data extremely well. The testing
$ {R}^2 $
scores range from
$ 0.960 $
to
$ 0.983 $
for both vertical and horizontal convergence, indicating a very good generalization to unseen data. The consistently high
$ {R}^2 $
scores and tight clustering of data points confirm the model’s robustness and reliability in predicting convergence accurately at the selected time instances.
Figure 5. Accuracy evaluation of the GP-based surrogate model. The true versus predicted outputs for
$ 5 $
out of the
$ 40 $
GPs that constitute the surrogate model, for both a. vertical and b. horizontal convergence at different time instances. Each GP approximates the convergence at a different time instance as indicated in the title of each subplot. The blue dots represent training data, and the red dots represent testing data. The black diagonal line denotes the line of perfect prediction.
Sensitivity analysis: With the GP-based surrogate model at hand, we compute the first-order (S1) and the total-order (ST) Sobol’ indices. Figure 6 displays the results of the time-dependent sensitivity analysis. From the figure, it is clear that the effect of each parameter varies over time. The parameters for primary creep
$ {\eta}_p $
and
$ {E}_p $
have a significant impact at the beginning of the simulation, but this effect diminishes as time progresses. Conversely,
$ {\eta}_s $
and
$ {p}_s $
initially have little influence but become more impactful over time as they are parameters of secondary creep. In particular,
$ {\sigma}_{0, eq,p} $
,
$ {p}_p $
, and
$ {\sigma}_{0, eq,s} $
have little to no influence on the output variation. Therefore, these material parameters can be fixed by taking the reference value from Stahlmann et al. (Reference Stahlmann, Missal and Gährken2016) and removed from further analysis.
Figure 6. Time-dependent global sensitivity analysis. First and total order Sobol’ indices over time for a. vertical and b. horizontal convergences, showing the influence of the primary creep parameters (
$ {\eta}_p $
,
$ {\sigma}_{0, eq,p} $
,
$ {p}_p $
,
$ {E}_p $
) and secondary creep parameters (
$ {\eta}_s $
,
$ {\sigma}_{0, eq,s} $
,
$ {p}_s $
), presented from top to bottom, respectively.
The indicators defined in Section 4 are computed from the time series shown in Figure 6 to quantify the effect of each parameter on the model output variation more precisely. Figure 7 presents these indicators as bar plots for both the first-order (S1) and total-order (ST) Sobol’ indices. These bar plots clearly indicate which parameters most significantly affect the model outputs. Figure 7 complements our observation from Figure 6: the horizontal and vertical convergences were found to be sensitive only to the parameters
$ \boldsymbol{\theta} =\left[{\eta}_p,{E}_p,{\eta}_s,{p}_s\right] $
.
Figure 7. Aggregated global sensitivity analysis. Normalized first-order (a, b) and total-order (c, d) Sobol’ indices for vertical (a, c) and horizontal (b, d) convergences, showing the influence of the primary creep parameters (
$ {\eta}_p $
,
$ {\sigma}_{0, eq,p} $
,
$ {p}_p $
,
$ {E}_p $
) and secondary creep parameters (
$ {\eta}_s $
,
$ {\sigma}_{0, eq,s} $
,
$ {p}_s $
), from left to right, respectively. The cumulated first-order and total Sobol’ indices are normalized against their respective highest overall values for a fair comparison with the other indicators.
Convergence analysis of the Sobol’ indices was performed as specified in Section 3.3. The process begins with an initial set of
$ {2}^7=128 $
samples and iteratively increases the sample size by
$ {2}^7 $
in each step until either a maximum of
$ {2}^{16}=65536 $
samples is reached or convergence is achieved. From (3.8), it was concluded that
$ {2}^{14}=16384 $
samples are enough to assure Sobol’ indices analysis convergence according to criterion (3.8). In particular, the maximum changes for S1 and ST were
$ 0.0093 $
and
$ 0.0071 $
, respectively. The simulation time required for
$ {2}^{14} $
surrogate model evaluations was
$ 4.5 $
min.
5.3. Model calibration from in-situ monitoring data
The final task is to calibrate the four sensitive model parameters comprised in
$ \boldsymbol{\theta} $
using real drift convergence monitoring data. To achieve this, we embed the surrogate model trained in Section 5.2 into the optimization problem (3.2) using only the mean values of the GP-based surrogate prediction as defined in (4.5). Figure 8 displays the results obtained from the optimization problem. The results show that the simulation data generated from the high-fidelity model and the predictions from the GPs almost overlap, demonstrating the effectiveness of the surrogate model in calibration. It is crucial to highlight that model calibration is often a computationally demanding task, typically requiring multiple calls of the model. The global optimization of the remaining
$ 4 $
model parameters converged within
$ 29 $
iterations, which required
$ 1805 $
surrogate model evaluations. For comparison, training the GP-based surrogate model, used for both sensitivity analysis and calibration, only required
$ 150 $
model evaluations. The optimization took less than
$ 1 $
s to be completed on a standard laptop. This represents a dramatic reduction in the computational resources required to calibrate the deep-repository model, underscoring the efficiency of our surrogate modeling approach. The predicted parameter values agree very well with expert intuition so that the result can be considered reasonable.
Figure 8. Model calibration with monitoring data. Comparison of monitoring data (black squares), GP-based optimal surrogate model prediction (red dots), and corresponding FLAC3D simulation results using the optimal parameter values (blue line) for a. the vertical and b. horizontal convergences.
6. Conclusions
Deep geological disposal of hazardous materials requires robust numerical models to ensure long-term safety and stability. The calibration of such models with real-world monitoring data is essential for accurately reflecting in situ conditions and enhancing repository management. This study contributes to advancing digital twinning for deep geological disposals by automating the calibration process using GP-based surrogate models. The presented workflow combines sensitivity analysis, surrogate modeling, and optimization to enable efficient calibration of a mechanical model representing the behavior of an emplacement drift in rock salt formations located in the northern main drift of the Gorleben salt dome in Germany.
The results demonstrate the efficiency and accuracy of the proposed approach. Initially, training of the GP-based surrogate model with
$ 7 $
input parameters,
$ 150 $
training data samples and
$ 80 $
outputs took
$ 12.3 $
s. The subsequent accuracy evaluation yielded
$ {R}^2 $
scores between
$ 0.960 $
and
$ 0.983 $
. Afterwards, a sensitivity analysis using time-dependent Sobol’ indices was performed with
$ \mathrm{16,384} $
surrogate model calls within
$ 4.5 $
min to identify four relevant material parameters. Finally, the GP-based surrogate model was calibrated based on
$ 14 $
year of convergence measurements, including the convergence of the global optimization in
$ 29 $
iterations,
$ 1805 $
model evaluations for gradient construction, and a duration of less than
$ 1 $
s. The surrogate model prediction provided both very good agreement with the monitoring data and valid values for parameters of the constitutive model TUBSsalt.
This approach reduces the computational burden associated with traditional high-fidelity models and enables rapid, iterative updates to model parameters as new monitoring data becomes available. By enhancing the scalability and adaptability of numerical models, this work lays the foundation for integrating advanced surrogate modeling techniques into the management of deep geological repositories.
While the workflow demonstrated high efficiency and accuracy for the presented mechanical model, further developments are required to extend its applicability. Specifically, first, the GP-based surrogate model calibration and optimization method must be extended to account for higher-dimensional monitoring data, which could include extensometer and permeability measurements. Second, an alternative formulation would be developing a time-dependent surrogate model, which enables forecasting capabilities together with uncertainty propagation. Third, efforts have to be undertaken to account for the multi-physics nature of deep geological disposal. Apart from the purely mechanical model investigated here, a variety of models has been developed in the past that describe, for example, transport of radio-nuclides through fluid flow, heat generation from high-level radioactive waste, or hydration of a sealing structure.
Data availability statement
All source code and simulation data used for the presented benchmark studies are published at https://doi.org/10.5281/zenodo.15049774 (Paul et al., Reference Paul, Urrea-Quintero and Fiaz2025). Restrictions apply to the availability of the monitoring data, which were used under license for this study. Monitoring data are available from the authors with the permission of the Bundesgesellschaft für Endlagerung mbH (BGE).
Acknowledgments
We are grateful for the provision and preparation of monitoring data from the Bundesgesellschaft für Endlagerung mbH (BGE).
Author contribution
Conceptualization: L.P., JH.UQ., U.F., A.H., and H.W.; Methodology: L.P., JH.UQ., U.F., H.Y., and H.W.; Software: L.P., JH.UQ., U.F., and A.H.; Data curation: L.P., U.F., and A.H.; Data visualization: L.P., JH.UQ., and U.F.; Formal analysis: L.P., JH.UQ., and U.F.; Writing original draft: L.P., JH.UQ., U.F., A.H., H.Y., H.W., U.R., and J.S.; Funding acquisition: J.S., U.R., and H.W.; All authors approved the final submitted draft.
Funding statement
This research is funded by the German Federal Ministry for the Environment, Nature Conservation, Nuclear Safety, and Consumer Protection (BMUV) and managed by project management agency Karlsruhe (PTKA) under grant number 02E12102.
Competing interests
The authors declare none.
Ethical standard
The research meets all ethical guidelines, including adherence to the legal requirements of the study country.
A. Appendix—Preliminary parameter selection
A preliminary parameter selection was performed to identify relevant strain components of TUBSsalt given the considered monitoring data. For this purpose, several damage-associated strain components were deactivated, in particular tertiary creep
$ {\dot{\varepsilon}}_t $
, creep and shear failure
$ {\dot{\varepsilon}}_n $
and tension failure
$ {\dot{\varepsilon}}_z $
. A comparison between
$ 274 $
gridpoint displacements at
$ 40 $
time instances is given in Figure A1 below, in which the vertical (Figure A1a) and horizontal (Figure A1b) displacements were simulated once with and once without damage strains. From the comparison, it can be observed that no deviation from the diagonal line is visible, which is also confirmed by the
$ {R}^2 $
value close to
$ 1.0 $
. Thus, the considered test case is dominated by creep mechanisms and softening, and post-failure strains are considered negligible.
Figure A1. Comparison of displacements with and without damage-associated strains by evaluating
$ 274 $
gridpoints at
$ 40 $
time instances for a. vertical and b. horizontal displacements. The simulations use TUBSsalt parameters:
$ {\eta}_p=8.0\cdot {10}^4 $
MPa
$ \cdot $
d,
$ {E}_p=75 $
MPa,
$ {\sigma}_{0, eq,p}=30 $
MPa,
$ {p}_p=0.5 $
,
$ {\eta}_s=3.0\cdot {10}^7 $
MPa
$ \cdot $
d,
$ {\sigma}_{0, eq,s}=30 $
MPa and
$ {p}_s=1.5 $
.
B. Appendix—Material parameter
Table B1 contains parameter values and ranges for the constitutive model TUBSsalt. As reference serves a parameter set determined for Gorleben salt in Stahlmann et al. (Reference Stahlmann, Missal and Gährken2016). However, it should be noted that the parameter set has been identified using data from various salt formations of Gorleben, and not only from the homogeneous z2HS2 area. As no parameters for healing are given in Stahlmann et al. (Reference Stahlmann, Missal and Gährken2016),
$ {\eta}_v $
and
$ {m}_v $
are obtained from Epkenhans et al. (Reference Epkenhans, Wacker and Stahlmann2022).
Table B1. TUBSsalt material parameter under consideration of reference values from Stahlmann et al. (Reference Stahlmann, Missal and Gährken2016). The parameters highlighted in red are directly read from experimental data. For the parameters highlighted in blue, a range of values is indicated to perform the sensitivity analysis.
Open access
Comments
No Comments have been published for this article.