Effects of Glacial Isostatic Adjustment on Fault Reactivation and Its Consequences on Radionuclide Migration in Crystalline Host Rocks

Kern, Dominik; Magri, Fabien; Malkovsky, Victor; Steffen, Holger; Nagel, Thomas Article — Published Version Effects of Glacial Isostatic Adjustment on Fault Reactivation and Its Consequences on Radionuclide Migration in Crystalline Host Rocks Environmental Modeling & Assessment Provided in Cooperation with: Springer Nature Suggested Citation: Kern, Dominik; Magri, Fabien; Malkovsky, Victor; Steffen, Holger; Nagel, Thomas (2024) : Effects of Glacial Isostatic Adjustment on Fault Reactivation and Its Consequences on Radionuclide Migration in Crystalline Host Rocks, Environmental Modeling & Assessment, ISSN 1573-2967, Springer International Publishing, Cham, Vol. 30, Iss. 1, pp. 177-192, https://doi.org/10.1007/s10666-024-09997-3 This Version is available at: https://hdl.handle.net/10419/323311 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by/4.0/ Environmental Modeling & Assessment (2025) 30:177–192 https://doi.org/10.1007/s10666-024-09997-3 RESEARCH Effects of Glacial Isostatic Adjustment on Fault Reactivation and Its Consequences on Radionuclide Migration in Crystalline Host Rocks Dominik Kern1·Fabien Magri2,3 ·Victor Malkovsky4·Holger Steffen5·Thomas Nagel1 Received: 9 November 2023 / Accepted: 6 September 2024 / Published online: 24 September 2024 © The Author(s) 2024 Abstract To assess the robustness of a safety case for a deep geological repository (DGR), it is necessary to analyze a range of scenarios covering likely, less likely, and hypothetical future developments. Crystalline rock can, under ideal conditions, provide a suitable hydrogeologic barrier due to its extremely low matrix permeability. However, this host rock is often fractured, which can compromise its hydro-mechanical (HM) barrier function. We quantify how faults that are prone to reactivation during glacial events can affect radionuclide migration around a DGR in a crystalline host rock. We extend a previously developed finite element model of coupled fluid flow and radionuclide transport to numerically solve the component transport problem before and after fault reactivation. Assuming that fault reactivation is triggered by changes in mechanical boundary conditions, we derive heterogeneous permeability distributions in the reactivated faults by evaluating the Coulomb failure stress criterion of finite element solutions of a complementary hydro-mechanical problem. Specifically, we evaluate the consequences of glacial isostatic adjustment (GIA) during a glacial cycle. We find that the increased permeability in the reactivated faults accelerates the migration of radionuclides along the fault by channeling the flow, while it is reduced in the direction perpendicular to the fault. The channeling observed is also a result of heterogeneous permeability enhancement, and the flow fields differ from those of the previous model which postulated a homogeneous permeability enhancement. Although the proposed numerical workflow has been applied to the case of GIA, it is adaptable to study hydro-mechanical processes induced by seismic events or by hydrofracking in enhanced geothermal systems. Keywords Deep geological repository ·Radionuclide migration ·Advection–diffusion transport ·Glacial isostatic adjustment · Fault reactivation ·Coulomb failure stress ·Poroelasticity 1 Introduction To date, the safest solution for the disposal of nuclear waste is considered to be storage in deep geological repositories (DGRs). Regardless of differences in national regulations, the time scales of radioactive decay require that nuclear waste be isolated from the biosphere for long periods of time, on the order of millions of years. A safety assessment must therefore consider not only the current state of the site, but also possible future states. Geoscientific studies provide evidence of alternating cold and warm periods on geological time scales. During such periods, changes in precipitation rates, continental glaciation, and marine transgressions can alter hydraulic, geochemical, Extended author information available on the last page of the article and geomechanical states and properties [1], including in situ stresses [2]. Seismic events, ground rupture due to faulting, uplift, or subsidence due to glaciation can nucleate fractures. Glaciation, expansion of permafrost beneath and beyond the ice sheets, as well as thermo-elastic effects, can also generate large stress changes and hydraulic gradients with the potential to cause rock failure [3]. These couplings between glacial loading and rock mechanical response affect hydraulic permeability in terms of opening and closing existing pathways. In addition, new fractures may develop, providing additional pathways for flow and radionuclide transport [4–8]. Based on palaeodata, this conceptual paper presents numerical simulations investigating the effects of glacial cycles (i.e., advance and retreat of a glacier) on the hydrogeological regime of a fictitious repository in a crystalline host rock. The significance of fault activation has been 123 178 D. Kern et al. emphasized in recent research, along with its implications for thermo-hydro-mechanical (THM) modeling of glacierinduced subsidence/uplift, subsurface particle transport, and tracer breakthrough curves [9,10]. It is essential to comprehend these processes in order to evaluate the long-term security of deep geological repositories. Fault activation, for instance, can change radionuclide and tracer breakthrough curves as well as drastically impact pollutant movement paths. Furthermore, lithological heterogeneity and hydrological gradients, among other factors, have a significant impact on subsurface particle transport, which is crucial in determining the dispersion of contaminants over time. The stability of nuclear waste sites is impacted by the mechanical response of rock formations to past and future glacial cycles, which can be qualitatively estimated by THM modeling of glacial subsidence and uplift. Although the THM models presented here have significant limitations due to the assumptions and unknowns involved, they provide insights into the processes of HM coupling and fault reactivation. Consequently, they play an important role in the overall safety assurance and design of DGRs, but should not be considered independent performance assessments. In contrast to previous studies using either one-dimensional ordinary differential equations (ODE) or two-dimensional finite element method (FEM) models [11–15], or DEM models [16,17], we employ three-dimensional FEM models, being competitive in their corresponding problem classes [18,19]. The only limitation to detecting fault reactivation in postprocessing is that the proposed approach does not account for other fracturing than in the existing faults and subsequent permeability and stiffness changes caused by reactivation. In a previous study [20], it was shown that the orientation of faults can affect their influence on the migration of a contaminant plume and that ecological hazard could arise depending on the distance between the repository and the fault. The simulations presented here focus on the effects of a variable stress state on the hydraulic permeability of the faults and consequent radionuclide migration, which were previously neglected. The main novelty is the combination of a hydro-mechanical and component transport simulation in one specific model instead of generalized relations between them. The assumptions about the repository, geomechanical parameters of the rock mass, and the boundary conditions are detailed in Sect.2, followed by the description of the physical model and its discretization into finite elements and time steps in Sect.3. The results of exemplary scenarios for the dispersion of radionuclides in a worst-case scenario of radionuclide leakage are presented in Sect.4and discussed in Sect.5. 2 Data and Assumptions Geological processes potentially affecting host rock properties vary largely on spatial and temporal scales [21,22]. Because DGRs are typically located in seismically inactive regions, we place emphasis on a scenario in which the geodynamic process of glacial isostatic adjustment (GIA) [23]isthe driving force behind stress-induced permeability changes. GIA affects continental-scale areas, even outside of the glaciated region, over many thousands of years. A putative geological structural model representative of a DGR in crystalline rocks is shown in Fig. 1. The weight of the ice sheet pushes the lithosphere downwards and bulges the lithosphere around the ice sheet. As the ice sheet retreats, the crust rebounds. Due to the viscoelastic nature and therefore time-dependent behavior of the mantle underneath, this rebound is time-delayed and can therefore last many thousand years after the ice vanished [24]. Unlike local effects below the ice sheet [25,26], these movements affect areas far from the ice sheet. Understanding these processes and their effects therefore requires modeling on a global scale with a dedicated GIA model [27,28]. We Fig. 1 DGR scale and location in the host rock (bottom left zoom) with respect to a putative repository site affected by glaciation. The depth of the DGR (green) is about a hundred meters and some hundred meters away from potential faults 123 Effects of Glacial Isostatic Adjustment on Fault Reactivation... 179 Table 1 Geometry of the site, faults, and DGR [20]Feature Domain Meridional fault 6.5km<x<6.6km, −∞ <y<∞,−∞ <z<∞ Latitudinal fault −∞ <x<∞,5.5km<y<5.6km, −∞ <z<∞ Deep geological repository 7.0km<x<7.3km, 3.5km<y<5.0km, −125 m <z<−75 m Site 0.0km<x<12.0 km, 0.0km<y<18.0km, −1000 m <z<120 m to 430 m evaluate results from the well-established ICE-6G_C(VM5a) model [29,30] that describes geodynamically constrained Earth structures, adjusted to the last glaciation. A submodel technique is used to transfer data from the GIA scale to the repository-scale model. Specifically, the results of the GIA simulation are interpolated to a kilometer-scale site model as prescribed time-dependent displacement boundary conditions. We assume that the terrain at the site is not glaciated, but that the ice sheet extends at least a few hundred kilometers and is thus of continental scale. The surface topography is inferred from a digital elevation model causing pressure-driven groundwater flow (i.e., regional flow). Based on previous work by Malkovsky et al. [20], we consider two representative faults in the model, a latitudinal (east– west) and a meridional (north–south) fault. Either one of these prototypical faults or both are reactivated. In addition, we consider heterogeneous fault reactivation with a spatially variable permeability in the fault depending on local stress changes and compare it to the previously studied case of homogeneous fault reactivation. As indicated by experimental results [31,32], we run HM simulations to detect regions of fault reactivation, in contrast to the theoretical limit cases where the homogeneous fault permeability value is taken as the maximum value that can occur in the heterogeneous faults. The DGR and the faults are geometrically discretized as cuboids with the faults being thin structures cutting the entire domain from top to bottom and side to side. Their locations and dimensions are listed in Table 1and are in line with the previous study [20]. Note that z=0 is an arbitrary reference level in the modelling domain not necessarily coinciding with the sea level. Figure 2shows the hydraulic permeability within the domain, which is representative of all the heterogeneously distributed parameters, in both intact rock and fault. After fault reactivation, the heterogeneity gets more pronounced in the faults. The 3D permeability field is based on field measurements [33] and computed from a stochastic model [34]. For some of the remaining parameters, minima and maxima are known from tests. Their distributions are derived from the permeability distribution by linear interpolation. Porosity is assumed to increase with permeability, while density and Young’s modulus are assumed to decrease. From the spectrum of radionuclides in high active waste, we have chosen the long-lived Americium-241 (half-life t1/2=432.6a) as representative. We considered the conservative case, when the radionuclide is carried by groundwater in the most mobile form without any retardation, since experimental work [35] suggests that they move through gneiss in highly mobile colloidal form. In addition, we neglect pore diffusion, since Fig. 2 Mesh for the component transport simulation, exaggerated scaling in z-direction (factor 3). The distances between DGR (green) and faults are 400m and 500m. The two vertical slices show the initial heterogeneous permeability distribution [34]. The magnified image shows the mesh refinement towards the faults 123 180 D. Kern et al. dispersivity dominates the component transport in rock. All physical parameters of the numerical example are listed in Table 2. 3 Models and Methods The problem is split into two loosely coupled parts for computational efficiency because of disparate time scales and mesh resolution requirements. At first, we perform coupled hydro-mechanical simulations assuming GIA-driven deformation as the primary cause of stress changes to estimate the resulting permeability changes. The time scale of these simulations is given by glacial cycles. Here, we refer to the last cycle (110 ka). Then, we perform component transport simulations of a leakage scenario with and without fault reactivation. For both types of simulations, we use the finite element code OpenGeoSys [37,38], which is fully open-source. In our implementation, we applied the freely available MKL library by Intel, particularly utilizing the PARDISO direct solver for linear equation systems due to its robustness and efficiency. It should be noted that Intel MKL is freeware, not open-source, as detailed in its license agreement. However, OGS and the models presented here can be solved with linear solvers from open-source libraries. Figure 2shows the mesh used for the component transport simulation and highlights the subdomains of the DGR and the assumed faults. The faults are modeled as full-dimensional domains (fault thickness h>0) of porous media. 3.1 Hydromechanical Simulation We use the following equation set (u-pform) for quasi-static, fully saturated, poroelastic media [39] with solid displacement vector uand pore pressure pas primary variables: α∂ε ∂t+S∂p ∂t+∇·q=0 with ε=∇·u and q=−k μ∇p−ρfg,(1) ∇·σ+ρg=0with σ=C:ε−αp1 and ε=∇ s⊗u.(2) The fluid compressibility Cfand the solid compressibility Csare included in the storativity: S=φCf+(α −φ)Cs,(3) where φdenotes porosity and αthe Biot coefficient. Further, qdenotes specific discharge vector, ggravity vector, Cfourth-order elasticity tensor, εstrain tensor, εvolume strain, and 1second-order unit tensor. The remaining parameters are listed in Table 2. We impose atmospheric pressure at the top and impermeable boundaries elsewhere (watersheds, impermeable bedrock) for the hydraulic (Eq.1). Further, we impose atmospheric pressure as vertical stress on the top boundary, vanishing displacements in normal directions at Table 2 Physical properties of host rock, faults, and DGR [20, 36] Name Symbol Value and units Fluid density ρf1·103kg m−3 Fluid viscosity μ1·10−3Pa s Fluid compressibility Cf4·10−10 Pa−1 Solid density ρsρs=ρs,max −k−kmin kmax−kmin ρs,max −ρs,min ρs,min =2.7·103kg m−3,ρs,max =3.5·103kg m−3 Young’s modulus (bulk, drained) EE=Emax −k−kmin kmax−kmin Emax −Emin Emin =5.63 ·1010 Pa, Emax =7.87 ·1010 Pa Poisson’s ratio (bulk, drained) ν0.27 Porosity φφ=φmin +k−kmin kmax−kmin φmax −φmin φmin =0.15%, φmax =0.5% Biot coefficient α1.0 Permeability (rock) kInterpolation of field-data [34] kmin =8·10−18 m2,kmax =1·10−13 m2 Permeability (homogeneous fault) khom 1·10−12 m2 Permeability (heterogeneous fault) khet See Eq.10 Retardation factor R1 (neutral tracer) Pore diffusion D00m 2s−1(no pore diffusion) Longitudinal dispersivity αL100 m Transversal dispersivity αT100 m Decay constant (241Am) κ16 ·10−4a−1 123 Effects of Glacial Isostatic Adjustment on Fault Reactivation... 181 all sides (north, east, south, west) and a prescribed, timedependent vertical displacement uGIA z(t)at the bottom for the mechanical (Eq. 2). The time-dependent vertical displacement at the lower boundary represents the forcing condition due to GIA. This boundary condition is derived from the ICE model, which provides the displacement field resulting from the varying glacier loads. This approach is justified because it directly incorporates the effects of glacial loading and unloading on the Earth’s surface, which are more accurately represented by changes in displacement than by changes in stress at the surface. The imposed vertical displacement uGIA z(t)captures the viscoelastic response of the Earth’s crust to glacial cycles. To further illustrate this boundary condition, we provide a plot of the imposed time-varying displacement at the lower boundary in Fig. 3. We notice a dominating rigid body mode (domain moves as a whole), but also differences between vertex displacements leading to shear deformation. In order to establish equilibrium, we first perform a presimulation in which we move from an unstressed initial state to a resting hydrostatic state. Before any glacial loading is applied, this pre-simulation ensures that the initial stress level represents a balanced state. To ensure stability at the beginning of the simulation, the initial conditions for faults are set so that they are not close to failure. To check for fault reactivation, we compute the load on the faults in terms of normal stress and shear stress from the simulations σn=(n⊗n):σ,(4) τabs =σn−σnn,(5) with normal vector nof the fault plane. Additionally, we assume a Coulomb model for the strength of the faults τfail =−σntan ϕ+τcoh,(6) with angle of internal friction ϕand cohesion τcoh. Tensile stresses are defined positive, σn>0. We evaluate the Couloumb failure stress (CFS) criterion [40], either in an absolute or relative formulation τCFS =τfail −τabs,(7) τrCFS =1−τabs τfail ,(8) assuming zero cohesion τcoh =0. That is a conservative assumption on the uncertain cohesion of existing faults on long time scales which is lower than the cohesion in undamaged rock [22]. If the CFS drops below the critical value τCFS <0(9) at any moment, these points are considered irreversibly activated. The updated permeability of each finite element inside the fault region depends linearly on how many of its nodes got activated kel =nactivated ntotal kfault,(10) where nactivated and ntotal denote numbers of activated nodes and total number of nodes per element, respectively. In our simulations, Eq.10 governs the postprocessing stage, where Fig. 3 Vertical displacements that were interpolated from the global GIA model to our modeling domain 123 182 D. Kern et al. the permeability of the elements within the fault zone is updated based on nodes activating owing to fault reactivation. Given the ratio of activated nodes to total nodes, this formulation implies an increase in fault permeability when fault failure is triggered. It is crucial to remember that Eq.10 and spatial discretization in our simulation architecture are responsible for capturing the localized increase in permeability that occurs within the fault zone following reactivation. The permeability of elements within the fault zone increases dramatically in comparison to the intact rock or the pre-reactivation state of the fault, even if only a small percentage of nodes “break” and cause fault reactivation. The models were discretized with Taylor-Hood elements (quadratic approximations of displacements and linear approximation of pressure). An Euler-backward time-discretization advances the transient equations. The HM-equations were set up monolithically and solved by Newton iterations. The specific values for the numerical solution procedure are listed in Table 3. 3.2 Component Transport Simulation We use the same hydraulic Eq.1without mechanical coupling ε=0 and an advection-dispersion model for the component transport [41] S∂p ∂t=∇·k μ∇p−ρfg,(11) ∂ ∂tRc+q·∇c φ=∇·D∇c−κRc,(12) with retardation factor R=1+ρsKd(1−φ) φ(13) and dispersion tensor D=D0I+(αL−αT)q⊗q q+αTqI,(14) which simplifies on our assumption of a neutral tracer and vanishing pore diffusion. Table 3 Algorithmic and discretization parameters of the hydromechanical simulation (HM) Name Value Time step 1000 a Finite elements 219,813 tetrahedral elements Newton iteration termination εp2<1·10−4Pa and εui2<1·10−8m Table 4 Algorithmic and discretization parameters of the component transport simulation (HC) Name Value Time step 1 a Finite elements 1,361,603 linear tetrahedral elements Picard iteration termination εp2<1·10−6Pa and εc/c02<1·10−4 Coupling iterations termination εp2<1·10−4Pa and εc/c02<1·10−2 In a transient state, this storage modeling without mechanical coupling may underestimate fluid flows; however, since we are in a quasi-static regime, we assume that the grain constituents of the formation are incompressible. All boundaries are impermeable to radionuclide transport and fluid flow (bedrock, water divides). We assume a general waste model [20] whose concentration decays following a time-dependent law cBC leak =c0e−κt,(15) with initial concentration c0and decay constant κ. The models were discretized with linear finite elements (linear approximation of pressure and concentration). An Euler-backward time-discretization with stabilization by the flux-corrected transport scheme (FCT) [42] advances the transient equations. The FCT scheme was implemented as a staggered scheme solving the separate processes with Picard iterations. The algorithm is modified to avoid overdiffusive results. This numerical procedure has been validated on benchmark1. The specific values for the numerical solution procedure are listed in Table 4. The mesh has a variable resolution in order to be able to resolve the flow in the narrow faults while keeping computational times reasonable. We considered two element layers in the faults, as apparent from the zoomed region in Fig. 2. The meshing, for both hydromechanics and component transport simulations, was a tradeoff between accuracy and computational costs. We have performed a convergence analysis on selected subproblems (see AREHS report2). Time discretization by backward Euler is sufficient for the types of partial differential equations we have in our quasi-static modeling. 4 Results Hydromechanical simulations (HM) are used to evaluate fault reactivation. Afterwards, component transport simu1OGS benchmark on Flux-corrected Transport. 2report https://zenodo.org/records/11367280 Coming online soon. 123 Effects of Glacial Isostatic Adjustment on Fault Reactivation... 183 lations (HC) reveal the effects of mechanically induced permeability changes on radionuclide contamination. As outlined in Sect.2, different fault activation scenarios are considered. 4.1 Hydromechanical Simulation of Fault Reactivation As the models are designed to cover one glacial cycle, each simulation starts from an equilibrium state reached 110,000 years ago and runs in time steps of 1000 years to the present. The models calculate the relative Coulomb failure stress and the corresponding slip direction. If the relative Coulomb failure stress Eq.9drops below zero once, then the corresponding node is irreversibly considered reactivated (no healing). A postprocessing tool finds all reactivated nodes and performs element-wise permeability updates according to Eq. 10. By assumption, the reactivation affects only the predefined fault zones (localization). Figure4shows the updated permeability field. Reactivation occurs in inclined, mostly continuous bands stretching from 400 m depth up to the surface. It can be seen that on average, the permeability of the original fault (not reactivated) is k=1·10−16 m2with a maximum value of k=1·10−14 m2. GIA-induced stress field changes induce an increase in permeability of reactivated faults up to four orders of magnitude (maximal value kf=1·10−12 m2) in our example. An important result of our models is that fault reactivation primarily affects the upper part of faults. This finding is plausible, because the displacement boundary condition leads to deviatoric stress changes of similar magnitude across the domain, whereas isotropic stress increases linearly with depth. Other influences include initial stress levels, mechanical properties of fault zones, heterogeneity, and the relationship between fault permeability variations and stress redistribution. When analyzing the initial stress state, it is important to consider the circumstances that existed before the fault was reactivated. Our simulations start from a state of equilibrium reached around 110,000 years ago, which gives us a “fictive” distribution of initial stress within the geological formation. An important component in defining favorable conditions for fault reactivation is how the stress changes over time as a result of glacial loading and unloading. 4.2 Component Transport Simulations of Leakage Similarly to Malkovsky et al. [20], we assume that RN leakage occurs 1000 a prior to fault activation and run the component transport simulation for a further 2000 a. Consequently, all simulations share the same initial phase without reactivated faults starting at t=0auntil t=1000 a, when one of the faults or both are reactivated. We then continue with simulations from t=1000 auntil t=3000 afor different fault reactivation scenarios. The component transport simulations start from an uncontaminated domain, i.e., zero RN concentration, and assume leakage according to the model Eq. 15 with reference concentration c(t=0a)=c0. These simulations start from steady-state hydraulic conditions and use a constant time step of 1 a. It should be noted that these time steps, and any time periods mentioned in the description of the results, do not refer to an actual geological time scale, but rather should be seen as a time interval that allows transient fluid flow patterns to develop and RN plumes to propagate in the system. RN concentrations are therefore normalized to illustrate Fig. 4 Updated permeability along the fault planes derived from Coulomb failure stress calculations (HM models). The maximal values (yellow patterns) correspond to fault reactivation. The continuous high permeability zones (yellow/orange stripes) provide preferential pathways for channelized hydraulic flow. The dashed lines indicate orthogonal slices through the center of the DGR (green projections) and a fourth slice called middle plane to be used for 2D plots 123 184 D. Kern et al. relative concentration changes over time and should not be interpreted as the actual RN concentration in the biosphere. The calculated RN distribution at t=1000 a,asshownin Fig. 5, is used as initial radionuclide distribution in all fault reactivation simulations. At this time, the concentration at the DGR sides decayed to c=0.202c0. We observe a westward and southwestward migration of radionuclides following the regional hydraulic flow. The concentration isoline c/c0= 10−5traveled about 1 km in the fastest direction during these 1000 years. The following scenarios, in which the heterogeneous permeability distribution in potentially reactivated faults is determined from HM simulations according to the CFS criterion Eq.9, are considered: Scenario 1: Meridional fault (N-S) is reactivated. Scenario 2: Latitudinal fault (W-E) is reactivated. Scenario 3: Both faults are reactivated. Each of these scenarios is also run for the case of homogeneous faults, i.e., the faults’ permeability is set equally in the entire fault domain (Table 2). In addition, a reference case Fig. 5 Radionuclide distribution at t=1000 aon orthogonal cross sections through the center of the DGR (xDGR =7150 m, yDGR =4250 m, zDGR =−100 m), indicated in redontheleft.ThisRN distribution is used to initiate the simulations with fault reactivation (scenarios 1–3) and without fault reactivation (reference case). Faults are indicated in gray 123 Effects of Glacial Isostatic Adjustment on Fault Reactivation... 191 Preand postprocessing scripts, as well as input and result files, are available on request from the corresponding author. Declarations Competing interests The authors declare no competing interests. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. References 1. Neuzil, C. (2012). Hydromechanical effects of continental glaciation on groundwater systems. Geofluids, 12(1), 22–37. https://doi. org/10.1111/j.1468-8123.2011.00347.x 2. Martin, C., & Lanyon, G. (2003). Measurement of in-situ stress in weak rocks at Mont Terri Rock Laboratory, Switzerland. International Journal of Rock Mechanics and Mining Sciences, 40(7–8), 1077–1088. https://doi.org/10.1016/S1365-1609(03)00113-8 3. Boulton, G., Chan, T., Christiansson, R., Ericsson, L.O., Hartikainen, J., Jensen, M.R., Stanchell, F.W., Wallroth, T.: Thermohydro-mechanical (THM) impacts of glaciation and implications for deep geologic disposal of nuclear waste. In: Stephanson, O. (ed.) Elsevier Geo-Engineering Book Series vol. 2, pp. 299–304. Elsevier, (2004). 10.1016/S1571-9960(04)80057-0 4. Zhang, X., Ma, F., Dai, Z., Wang, J., Chen, L., Ling, H., & Soltanian, M. R. (2022). Radionuclide transport in multi-scale fractured rocks: A review. Journal of Hazardous Materials, 424, 127550. https://doi.org/10.1016/j.jhazmat.2021.127550 5. Bense, V., Gleeson, T., Loveless, S., Bour, O., & Scibek, J. (2013). Fault zone hydrogeology. Earth-Science Reviews, 127, 171–192. https://doi.org/10.1016/j.earscirev.2013.09.008 6. Person, M., Bense, V., Cohen, D., & Banerjee, A. (2012). Models of ice-sheet hydrogeologic interactions: A review. Geofluids, 12(1), 58–78. https://doi.org/10.1111/j.1468-8123.2011.00360.x 7. Jeong, W. C., Kim, J. Y., & Song, J. W. (2004). A numerical study on the influence of fault zone heterogeneity in fractured rock media. KSCE Journal of Civil Engineering, 8, 575–588. https://doi.org/10. 1007/BF02899582 8. Hudson, D.: Unsaturated flow and transport through a fault embedded in fractured welded tuff. Water Resources Research 40(4) (2003). https://doi.org/10.1029/2003WR002571 9. O’Brien, G., Bean, C., McDermott, F.: A numerical study of passive transport through fault zones.Earth and Planetary Science Letters 214(3-4), 633–643 (2003). https://doi.org/10.1016/S0012821X(03)00398-4 10. Crawford, J., Löfgren, M.: Modelling of radionuclide retention by matrix diffusion in a layered rock model. Technical Report R-1722, SKB (2019) 11. Jobmann, M.: Site-specific evaluation of safety issues for high-level waste disposal in crystalline rocks. final report of ursel project. Technical Report TEC–28-2015-AB, DBE Technology GmbH, Peine (2016) 12. Niemeyer, M., Hugi, M., Smith, P., Zuidema, P.: Kristallin-i performance assessment: First results from sensitivity studies. In: Geological disposal of spent fuel and high level and alpha bearing wastes, (1993) 13. Keesmann, S., Noseck, U., Buhmann, D., Fein, E., Schneider, A.: Modellrechnungen zur Langzeitsicherheit von Endlagern in Salzund Granitformationen. Technical Report GRS-206, Gesellschaft für Anlagenund Reaktorsicherheit, Germany (2005) 14. Malkovsky, V., Liebscher, A., Nagel, T., & Magri, F. (2022). Influence of tectonic perturbations on the migration of long-lived radionuclides from an underground repository of radioactive waste. Environmental Earth Sciences, 81(23), 537. https://doi.org/10. 1007/s12665-022-10635-y 15. Nasir, O., Fall, M., Nguyen, S. T., & Evgin, E. (2013). Modeling of the thermo-hydro-mechanical-chemical response of sedimentary rocks to past glaciations. International Journal of Rock Mechanics and Mining Sciences, 64, 160–174. https://doi.org/10.1016/j. ijrmms.2013.08.002 16. Fälth, B., Hökmark, H.: Modelling end-glacial earthquakes at Olkiluoto. expansion of the 2010 study. Technical Report POSIVAWR-12-08, Posiva, Finland (2012) 17. Kwon, S., & Min, K.-B. (2021). Fracture transmissivity evolution around the geological repository of nuclear waste caused by the excavation damage zone, thermoshearing and glaciation. International Journal of Rock Mechanics and Mining Sciences, 137, 104554. https://doi.org/10.1016/j.ijrmms.2020.104554 18. Matiatos, I., Varouchakis, E. A., & Papadopoulou, M. P. (2019). Performance evaluation of multiple groundwater flow and nitrate mass transport numerical models. Environmental Modeling & Assessment, 24, 659–675. https://doi.org/10.1007/s10666-0199653-7 19. Kadeethum, T., Lee, S., & Nick, H. (2020). Finite element solvers for Biot’s poroelasticity equations in porous media. Mathematical Geosciences, 52, 977–1015. https://doi.org/10.1007/s11004-02009893-y 20. Malkovsky, V., Nagel, T., Kern, D., Magri, F.: Radionuclide migration from an underground radioactive waste repository under the influence of tectonic fault emergence: The Nizhnekanskiy Massif (Siberia, Russia) example. Environmental Modeling & Assessment, 1–12 (2023). https://doi.org/10.1007/s10666-023-09893-2 21. Gudmundsson, A. (2011). Rock fractures in geological processes. Cambridge University Press.https://doi.org/10.1017/ CBO9780511975684 22. Wyllie, D.C.: Foundations on rock: Engineering practice, Second Edition. CRC Press, (2003). https://books.google.de/books? id=4RE4DwAAQBAJ 23. Steffen, H., & Wu, P. (2011). Glacial isostatic adjustment in Fennoscandia–A review of data and modeling. Journal of Geodynamics, 52(3–4), 169–204. https://doi.org/10.1016/j.jog.2011.03. 002 24. Kierulf, H. P., Steffen, H., Barletta, V. R., Lidberg, M., Johansson, J., Kristiansen, O., & Tarasov, L. (2021). A GNSS velocity field for geophysical applications in Fennoscandia. Journal of Geodynamics, 146, 101845. https://doi.org/10.1016/j.jog.2021.101845 25. Fischer, U.H., Bebiolka, A., Brandefelt, J., Cohen, D., Harper, J., Hirschorn, S., Jensen, M., Kennell, L., Liakka, J., Näslund, J.-O., et al.: Radioactive waste under conditions of future ice ages. In: Haeberli, W., Whiteman, C. (eds.) Snow and Ice-Related Hazards, Risks, and Disasters, pp. 323–375. Elsevier, (2021). https://doi.org/ 10.1016/B978-0-12-817129-5.00005-6 26. Egholm, D. L., Pedersen, V. K., Knudsen, M. F., & Larsen, N. K. (2012). Coupling the flow of ice, water, and sediment in a glacial 123 192 D. Kern et al. landscape evolution model. Geomorphology, 141, 47–66. https:// doi.org/10.1016/j.geomorph.2011.12.019 27. Wu, P., Steffen, R., Steffen, H., Lund, B.: Glacial isostatic adjustment models for earthquake triggering, pp. 383–401. Cambridge University Press, (2021). https://doi.org/10.1017/9781108779906. 029 28. Mitrovica, J. X., Milne, G. A., & Davis, J. L. (2001). Glacial isostatic adjustment on a rotating earth. Geophysical Journal International, 147(3), 562–578. https://doi.org/10.1046/j.1365-246x. 2001.01550.x 29. Argus, D. F., Peltier, W. R., Drummond, R., & Moore, A. W. (2014). The Antarctica component of postglacial rebound model ICE6G_C (VM5a) based on GPS positioning, exposure age dating of ice thicknesses, and relative sea level histories. Geophysical Journal International, 198(1), 537–563. https://academic.oup.com/gji/ article-pdf/198/1/537/17366620/ggu140.pdf 30. Peltier, W.R., Argus, D.F., Drummond, R.: Space geodesy constrains ice age terminal deglaciation: The global model. Journal of Geophysical Research: Solid Earth 120(1), 450–487 (2015). https://doi.org/10.1002/2014JB011176 https://arxiv.org/ abs/https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/ 2014JB011176 31. Rutqvist, J., Graupner, B., Guglielmi, Y., Kim, T., Maßmann, J., Nguyen, T. S., Park, J.-W., Shiu, W., Urpi, L., Yoon, J. S., et al. (2020). An international model comparison study of controlled fault activation experiments in argillaceous claystone at the Mont Terri Laboratory. International Journal of Rock Mechanics and Mining Sciences, 136, 104505. https://doi.org/10.1016/j.ijrmms. 2020.104505 32. Cappa, F., Guglielmi, Y., & De Barros, L. (2022). Transient evolution of permeability and friction in a slowly slipping fault activated by fluid pressurization. Nature communications, 13(1), 3039. https://doi.org/10.1038/s41467-022-30798-3 33. Ozerskiy, A.: Hydrogeology of the Archean crystalline rock massif in the southern part of the Yenisseyskiy Ridge (Siberian Craton). Universal Journal of Geoscience 5, 151–155 (2017). https://doi. org/10.13189/ujg.2017.050505 34. Malkovsky, V.I., Ozerskiy, A.Y.: Stochastic model of filtration properties distribution for enclosing rocks of an underground repository of radioactive waste on the basis of pumping tests. In: Zharikov, A.V., Chizhova, I.A. (eds.) Proceedings of the 15th international conference physico-chemical and petrophysical researches in earth sciences, pp. 159–162. IGEM RAS, (2014) 35. Malkovsky, V., Yudintsev, S., & Aleksandrova, E. (2018). Leaching of radioactive waste surrogates from a glassy matrix and migration of the leaching products in gneisses. Radiochemistry, 60, 648–656. https://doi.org/10.1134/S1066362218060140 36. Kochkin, B.T., Malkovsky, V.I., Yudinzev, S.W.: Scientific basis for assessing the safety of geological isolation of long-lived radioactive waste (in Russian). IGEM RAS, (2017) 37. Bilke, L., Fischer, T., Naumov, D., Lehmann, C., Wang, W., Lu, R., Meng, B., Rink, K., Grunwald, N., Buchwald, J., Silbermann, C., Habel, R., Günther, L., Mollaali, M., Meisel, T., Randow, J., Einspänner, S., Shao, H., Kurgyis, K., … Garibay, J.: OpenGeoSys. Zenodo (2022). https://doi.org/10.5281/zenodo.7092676 38. Bilke, L., Flemisch, B., Kalbacher, T., Kolditz, O., Helmig, R., & Nagel, T. (2019). Development of open-source porous media simulators: Principles and experiences. Transport in Porous Media, 130(1), 337–361. https://doi.org/10.1007/s11242-019-01310-1 39. Zienkiewicz, O.C.: Computational geomechanics with special reference to earthquake engineering. Wiley, (1999). https://books. google.de/books?id=JPRDAQAAIAAJ 40. Zoback, M. D. (2007). Reservoir geomechanics. Cambridge University Press.https://doi.org/10.1017/CBO9780511586477 41. de Marsily, G.: Quantitative hydrogeology: Groundwater hydrology for engineers. Academic Press, (1986). https://books.google. de/books?id=0fhOAAAAMAAJ 42. Boris, J.P., Book, D.L.: Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works. Journal of Computational Physics 11(1), 38–69 (1973). https://doi.org/10.1016/00219991(73)90147-2 43. Kurylyk, B. L., MacQuarrie, K. T., & McKenzie, J. M. (2014). Climate change impacts on groundwater and soil temperatures in cold and temperate regions: Implications, mathematical theory, and emerging simulation tools. Earth-Science Reviews, 138, 313–334. https://doi.org/10.1016/j.earscirev.2014.06.006 44. Bense, V., Ferguson, G., Kooi, H.: Evolution of shallow groundwater flow systems in areas of degrading permafrost. Geophysical Research Letters 36(22) (2009). https://doi.org/10.1029/ 2009GL039225 45. Cacace, M., & Blöcher, G. (2015). Meshit–A software for three dimensional volumetric meshing of complex faulted reservoirs. Environmental Earth Sciences, 74, 5191–5209. https://doi.org/10. 1007/s12665-015-4537-x Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Authors and Affiliations Dominik Kern1·Fabien Magri2,3 ·Victor Malkovsky4·Holger Steffen5·Thomas Nagel1 BDominik Kern [email protected]g.de Fabien Magri [email protected]und.de Victor Malkovsky malko[email protected] Holger Steffen holger[email protected] Thomas Nagel [email protected] 1Institute for Geotechnics, TU Bergakademie Freiberg, Freiberg, Germany 2Division Research/International, The Federal Office for the Safety of Nuclear Waste Management (BASE), Berlin, Germany 3Institute of Geological Sciences, Freie Universität Berlin, Berlin, Germany 4Institute of Geology of Ore Deposits, Petrography, Mineralogy, and Geochemistry, Russian Academy of Sciences, Moscow, Russia 5Geodetic Infrastructure, Lantmäteriet, Gävle, Sweden 123