Hydro-Mechanical Effects and Cracking in Opalinus Clay

Abstract

In this chapter, we investigate hydro-mechanical effects in the Opalinus Clay, especially those leading to cracking. We present a methodology comprising laboratory and field scale experiments, as well as the development and application of numerical approaches.

2.1 Objectives and Outline

In this chapter, we investigate hydro-mechanical effects in the Opalinus Clay, especially those leading to cracking. We present a methodology comprising laboratory and field scale experiments, as well as the development and application of numerical approaches. In Sect. 2.2, we present a short overview concerning the Opalinus Clay (OPA) formation, its behavior and applications. This potential host rock for radioactive waste disposal shows complex coupled behavior, especially related to shrinkage and shrinkage/desiccation-induced cracking.

Section 2.3 describes the Cyclic-Deformation (CD-A) experiment in the Mont Terri Rock laboratory, one of the focal points of the GeomInt2 project. The experiment allows the measurement of relative humidity, deformation and cracking at in-situ scale. Furthermore, a drilling campaign allowed to extract core material that has been used in different laboratory experiments described in this chapter. We monitor the fracture response of OPA and evaluate the influence of its anisotropy (bedding orientation) as well as the effect of different suctions on the load-displacement curves. These experiments are modeled within numerical approaches based on the methods described in Sect. 2.4. In the latter, a continuum based, unsaturated hydro-mechanical approach is extended with the phase-field method (Sect. 2.4.2) and with the discrete element method (Sect. 2.4.3).

After comparing the response of the material both at the laboratory and at the field scales, we propose an interactive visualization strategy in Sect. 2.6. Finally, we review the proposed methodology in Sect. 2.7 and discuss its extension, future applications and transfer.

2.2 Observations and Processes in Opalinus Clay

The Opalinus Clay formation (OPA) is considered as a potential host rock for radioactive waste disposal due to its low hydraulic permeability and ability to retain radionuclides. Opalinus Clay formations are located in Northern Switzerland (Nagra 2021) and South-Eastern Germany (Hoth et al. 2007; Bundesgesellschaft für Endlagerung (BGE) 2020). In comparison to other clay formations, the OPA is well characterized by numerous in-situ experiments, conducted in the framework of the international Mont Terri Rock Laboratory in Switzerland. Its applicability in various geotechnical applications including the disposal of radioactive waste have been investigated, e.g. Bossart et al. (2018) and references therein.

The formation was accumulated during the dogger period (Jura) sandy and sandy facies as well as related subfacies types (Lauper et al. 2018; Kneuker et al. 2021). OPA is an anisotropic and highly heterogeneous material and shows complex and coupled thermo-hydromechanical (THM) behavior, including shrinkage and swelling. Swelling is related to an expansion of the material. In the current study, we focus on shrinkage, which is related to bulk volume reduction induced by environmental changes such as variation of relative humidity (RH), that, in turn is related to changes in water content. This effect can be observed and measured both in laboratory and field scales (Wild et al. 2017; Bock 2000; Zhang et al. 2010). Furthermore, it can lead to shrinkage-induced or desiccation cracks. With focus on how such cracks develop, the behavior of OPA is investigated in the following via experimental and numerical approaches that are applied at different scales. In-situ monitoring and borehole extraction are carried out in the Mont Terri Rock Laboratory in order to characterize the THM-behaviour of the OPA (Bossart et al. 2018; Kneuker and Furche 2021). In particular, the CD (Cyclic Deformation) experiment was a valuable starting point to gather experimental data and experience at the in-situ scale and to investigate the hydraulic-mechanical coupling induced by swelling and shrinking of shaly facies of the OPA due to cyclic variations of air humidity (Matray and Möri 2012; Ziefle et al. 2017). In the follow-up project, the CD-A-experiment, these processes were investigated further in detail, this time in the sandy facies of the OPA. Consequently, the CD-A-experiment is one focal point of the GeomInt2 project and will be described in the following.

2.3 Monitoring and Data Integration

In the Cyclic-Deformation (CD-A) experiment (Ziefle et al. 2022c), coupled effects due to the excavation and seasonal changes in relative humidity (RH) are monitored regularly. The long-term measuring program focuses on a detailed characterization of the host rock as well as on the evolution of process variables like the pore pressure, the displacements and the water content amongst others. The evolution of the excavation damaged zone (EDZ), corresponding changes of the water content and monitoring of desiccation cracks within the first decimeters of the rock are the central points of the experiment. These measurements enable an increased understanding of coupled HM effects in the EDZ and provide a valuable basis for a multi-disciplinary investigation aiming at further development of numerical modelling approaches and constitutive equations to reliability capture coupled effects in the OPA.

The CD-A experiment is characterized by two horizontal niches (twins) that are excavated perpendicular to the strike of bedding and that are not covered with any kind of stabilization (like shotcrete). Consequently, a direct interaction between the climatic conditions and the host rock takes place. While the HM effects in the open twin are evoked by the excavation and the seasonal variation of temperature and relative humidity, the relative humidity in the closed twin remains constant due to a sealing at its entrance. With that, the desaturation effect is minimized as much as possible holding the relative humidity high. The different evolution of the measured relative humidity within the niches is pictured in Fig. 2.1. The sudden drops in the red curve (closed niche) are related to door openings and installations. The open niche shows shows a cyclic decrease and increase of RH correspondent to the winter and summer months, respectively.

Fig. 2.1

CD-A experiment: measured relative humidity [%] in the twin niches. For example, the period between September 2019 and February 2020 shows a decrease of the RH, i.e. desaturation. Note that the RH in the closed niche remains high and approximately constant

The changes in RH, especially those leading to an increase of suction, i.e., decrease in the capillary pressure, can lead to cracking. This phenomenon has been observed in both twins, but more intensely in the open niche. A detailed monitoring of the evolution of desiccation cracks along a horizontal line (scanline) in each niche is presented in Regard et al. (2021). The program includes monthly monitoring of the length and the aperture of the cracks initiated at the niche walls. The measurements have indicated significantly more desiccation cracks in the open than in the closed twin. The sum of crack apertures in the open niche is also larger than in its closed twin. Due to the increase of relative humidity in the summer months, the sum of the crack apertures decreases.

In addition to the RH and cracking measurements, the CD-A experimental campaign comprises measurements focusing on the characterization of the claystone as well as the observation of process variables like the convergence, pore pressure and permeability which are amongst others measured in the boreholes presented in Fig. 2.2. Furthermore, laboratory scale tests are part of the study and provide further understanding on the OPA behavior. These tests have been carried out at material from CD-A drill cores extracted in March 2021 in the context of the permeability measuring campaign as presented in Fig. 2.2. Multi-disciplinary investigations along the presented boreholes are discussed in Ziefle et al. (2022a) and underline the high heterogeneity of the host rock with different facies and sub-facies types and the impact of the excavation process and the ventilation conditions on the near field around the niches.

Fig. 2.2

Summary of boreholes from measuring campaign in 2022. Extensometers, mini-piezometers and permeability sensors are installed, among others

In the following, we describe the performed laboratory experiments that focus on anisotropy characterization and mechanical response of the material under different suctions.

2.3.1 Laboratory Experiments

2.3.1.1 Mechanical Characterization of Opalinus Clay

In the first phase of the GeomInt project, fracture propagation characteristics in Opalinus clay have been investigated both numerically and experimentally. In this regard, a set of three-point-bending experiments has been performed on rectangular samples of Opalinus clay. Given the status of bedding structure in obtained cores, only the configuration of loading being perpendicular to the bedding direction could be studied. Accordingly, the effect of structural anisotropy on sample strength as well as the fracture propagation behavior was investigated, albeit qualitatively, adopting different numerical methods. As an extension to the first phase of the project, a series of semi-circular bending experiments has been conducted in order to thoroughly investigate the direction-dependent strength characteristics of the Opalinus clay. The semi-circular configuration (Fig. 2.3) was chosen because of its minimal material requirements and relatively easier preparation procedure. The samples were taken from two different drill cores of sandy facies type, namely A24 and A26, near the CD-A twins (Kneuker and Furche 2021; Ziefle et al. 2022a, c).

Following ISRM suggested regulations (Kuruppu et al. 2014), the samples were prepared with radius of $36\times 10^{-3}$ m, thickness of $40\times 10^{-3}$ m and an artificial notch at center with length of $14\times 10^{-3}$ m. Given sensitivity of the clay material to water, all specimens were prepared under dry conditions using a diamond band saw with $2\times 10^{-3}$ m thickness. As summarized in Fig. 2.3, three type of samples, namely P-, Z- and S-samples, have been considered for the experiments. The naming convention is indicative of loading angle $\theta$, which represents the orientation between the isotropy plane and the loading direction. We have tested following loading directions $\theta = [0, 45, 90]^{\circ }$ that represent P-, Z- and S-samples, respectively.

Fig. 2.3

Schematic summary of SCB samples: (from left to right) P-samples, Z-samples and S-samples

In total, three samples, two from A26 cores and one from A24 core, for each loading angle were tested. The tests were performed under quasi-static loading conditions with strain rate of 0.03 mm/min. The typical laboratory results in terms of applied load versus the deflection of the samples obtained from A24 core drill are illustrated in Fig. 2.4a. The load-deflection curves in the three cases exhibited a similar pattern; a downward convex trend followed by an approximately linear increase in the load response until failure. Furthermore, the sharp drop of the load represents the brittle failure mechanism during the experiments. The mean value and the range between the maximum and minimum achieved failure loads are summarized in Fig. 2.4b. As shown, the maximum obtained load before failure, indicative of material toughness, increases as the angle between the loading direction and the isotropy plane increases. The considerable range among toughness of the samples with similar anisotropy structure could be attributed to various factors, e.g. micro defects generated in samples during the preparation process or the conditions at which cores were drilled.

Fig. 2.4

a Obtained Load-Deflection curve for Opalinus clay sample with different loading angles, b Mean value and range of the obtained peak loads for each loading angle. With increasing loading angle, an increase of the peak load is observed

Figure 2.5 shows the final fracture pattern of the three samples from A24 cores. Note that the crack does not propagate straight along the bedding direction for the P-sample due to the heterogeneity of the material. In the Z-sample, the crack propagation along bedding direction can be identified more clearly. Finally, the crack propagation in the S-sample also shows, initially, a strong dependency on the bedding direction, but continues its trajectory straight, as expected.

Fig. 2.5

Failure pattern and fracture path trajectory of a P-sample, b Z-sample and c S-sample. The cracks propagate approximately straight from the notch

2.3.1.2 OPA Characterization Under Different Humidity Conditions

As also observed at in-situ scale, in the CD-A experiment, changes in the RH lead to variations in the mechanical and fracture behaviors. We investigate this phenomenon in the laboratory by characterizing the mechanical response of the OPA at various RH. To do that, a customized desiccator cell has been configured and set up at the geomechanics and geotechnics laboratory of CAU. The idea originates from the first phase of the project, where the concept of a humidity controlled bending test setup was proposed (Sattari et al. 2021). As shown in Fig. 2.6, the customized desiccator cell contains four distinct parts, namely

• 1. Epoxy glass cell, equipped with RH and temperature sensors, to provide an isolated environment for the designed experiments

• 2. Salt bath

• 3. Interchangeable mechanical loading apparatus and

• 4. Adjustable camera tube fabricated for use of a micro-camera dedicated to optical displacement measurements.

Fig. 2.6

Customized Desiccator Cell equipped with relative humidity (RH) and temperature sensors, a salt bath, an interchangeable mechanical loading apparatus and an adjustable camera tube

The test procedure generally takes place in two steps; the equilibrium phase and the loading step. The sample is exposed to a certain RH during the equilibrium phase while being allowed to adjust with the environment. The humidity conditions inside the cell were achieved and controlled by applying different supersaturated salt solutions. For the Opalinus clay samples, the equilibrium normally occurs between 7 to 10 days. The equilibrium period varies due to different initial water content of the sample and the specified RH of the cell. The equilibrium conditions are monitored with consecutive weight measurements of each sample. The unique design of the cell provides an opportunity to assess mechanical characteristics of OPA under controlled humidity conditions with minimizing the environmental exposure of sample between equilibrium phase and loading step. In this section, anisotropic strength characteristics of Opalinus clay are studied in different humidity conditions using the semi-circular bending configuration. Table 2.1 summarizes the four different humidity conditions considered in this study. The obtained suction pressure for each salt solution using 3 samples each is estimated using Kelvin’s equation and based on the recorded RH and temperature during the experiments,

$$\begin{aligned} \psi = -\frac{R\,T}{V_{w0}\,M}\,ln\left( \frac{p}{p_0}\right) \end{aligned}$$

(2.1)

where the universal gas constant $R=8.314$ J/(mol K), T is the temperature in Kelvin, $V_{w0}$ is the specific volume of water and $M = 0.01801528$ kg/mol is the molar mass of the pore fluid. The term ${p}/{p_0}$ corresponds to the relative humidity (RH).

Table 2.1 Utilized salt solutions along with the achieved relative humidity and generated suction

The obtained experimental results are presented in terms of failure load as a function of imposed suction for each loading angle in Fig. 2.7. A general increasing trend for the material strength following the increase in the imposed suction and resultant decrease of water content can be clearly observed. This behavior is associated with the generation of capillary forces in a partially saturated porous medium due to the pressure difference between gas and fluid phases at the pores (Wild et al. 2015).

Fig. 2.7

Load-displacement curves under different RH. In general, the load increases with suction

Considering the aspects mentioned in this section, we aim at reproducing the shrinkage-behavior of OPA numerically in the following.

2.3.1.3 Additional Tests on Percolation

We utilized the drilled cores to carry out fluid percolation tests. These tests are important in the context of hydraulic fracturing. Knowing when the fluid will induce a crack aids in understanding critical conditions in the field. Fluid percolation experiments were carried out in the geomechanics and geotechnics laboratory of CAU Kiel. The main objective to these experiments is to investigate the characteristics of fluid driven fractures in Opalinus clay under effects of structural anisotropy as well as the anisotropy of the stress state. The utilized setup includes a true triaxial loading apparatus, resembling the in-situ stress conditions, and a syringe pump. More details related to experimental setup and the true-triaxial device can be found in Sattari et al. (2021). Cubic samples were prepared with length of $43\times 10^{-3}$ m and a scaled borehole drilled with diameter of $7\times 10^{-3}$ m and length of $20\times 10^{-3}$ m. Fluid is injected under constant flow rate into the sample through an adapter and a steel tube glued to the drilled borehole. A fluorescent oil was used as the injection fluid in order to trace even the tiniest fractures initiated during the process. Schematic representation of the applied principle stresses along with the bedding direction of the samples are shown in Fig. 2.8.

Fig. 2.8

Schematic representation of cubic S-, P- and Z-samples (left to right)

Given the position of the scaled borehole with respect to the bedding direction, three different loading configurations can be tested:

• 1. S-samples where the drilled borehole is perpendicular to the bedding direction

• 2. P-samples where the borehole is parallel to bedding and

• 3. Z-samples where bedding direction is inclined to the borehole.

Accordingly, Table 2.2 summarizes the number of tested samples for each bedding direction and the stress state in which the experiment was conducted.

Table 2.2 Summary of fluid percolation experiments on Opalinus clay

The experimental test starts after the in-situ stresses are applied. Once the mechanical equilibrium has been reached, the fluid injection then starts at flow rate of 10 ml/h. Upon observing the breakdown pressure, the experiment is stopped and the sample is retrieved. After careful investigation of failure patterns on the surface, the samples were sawed into number of cross section to have a better understanding of fracture propagation behavior inside it. Figure 2.9 shows the crack morphology on the surface, the outline of the cross sections, and the inner trajectories of fracture path indicated by florescent oil under UV light. It can be clearly seen that the pressure induced fractures inside Opalinus clay can either follow a straight path or propagate in branches and form a network. These observations are also reflected in the pressure response shown in Fig. 2.10. Generally, the drop in the pressure response is associated with increase in the occupied volume by the injection fluid followed by initiation and propagation of fractures inside the sample. In some cases, e.g., Fig. 2.10c, the initiated fracture reaches the surface of the sample rapidly and causes the depletion of the injection fluid out of the sample and a sharp breakdown in pressure response. In other cases, e.g., Fig. 2.10b, fractures propagate in tiny branches and form a network before reaching the sample surface which is reflected in the fluctuations of the pressure response.

Fig. 2.9

Fracture path trajectory, cross-section cut and inner fracture path trajectory indicated by florescent oil for a P-sample b S-sample and c Z-sample of Opalinus clay under isotropic in-situ stress of 12 MPa

Fig. 2.10

Comparison between pressure-time curve for different isotropy planes of Opalinus clay under mechanical stress configurations of a $S_\text {v}=12$ MPa, $S_\text {H}=12$ MPa, $S_\text {h}=12$ MPa, b $S_\text {v}=12$ MPa, $S_\text {H}=8$ MPa, $S_\text {h}=8$ MPa, and c $S_\text {v}=12$ MPa, $S_\text {H}=8$ MPa, $S_\text {h}=4$ MPa

2.4 Numerical Modeling Approach

To reproduce shrinkage effects observed in OPA in the laboratory and field scales (CD-A experiment), we utilize and extend the numerical methods evaluated in the first phase of the GeomInt project (Yoshioka et al. 2021; Kolditz et al. 2021a). The extensions are related to the anisotropic and cracking behaviors of the material.

Intuitively, cracks generate edge separations with creation of boundary surfaces over which the displacements are discontinuous. This makes numerical simulation of crack propagation challenging especially if we do not prescribe the location, topology, or time of their creation. For modeling crack propagation, several strategies are proposed in the literature. One way to deal with this sort of topology changes is by adaptively remeshing the domain (Bouchard et al. 2000; Meyer et al. 2004; Branco et al. 2015; Rangarajan et al. 2015). However, assigning boundary conditions to the newly created surfaces and including potential contact are not straightforward tasks. Another way to model crack propagation is to enrich elements that contain cracks to capture the displacement discontinuity, an approach known as eXtended/Generalized Finite Element Methods (XFEM, GFEM) (Moës et al. 1999; Belytschko et al. 2001; Gasser and Holzapfel 2005; Fries and Belytschko 2006; Belytschko et al. 2009; Khoei et al. 2012). Despite recent success (Shauer and Duarte 2022), the implementation efforts remain a daunting task especially when arbitrary fracture topologies are considered (Belytschko et al. 2001, 2003; Duarte et al. 2007).

Crack initiation needs to be prescribed in the methods discussed above. This can become a limitation in applications of geomechanical integrity, where we need to assess fracture initiation possibilities without any a-priori knowledge. The phase-field approach for brittle fracture represents the crack using a phase-field (damage) variable that results in a diffuse cracked interface. Furthermore, it does not require any special mesh treatment and prescription of the crack position. We apply and evaluate the response of both diffuse cracks using a continuum mechanics based approach combined with the phase-field and sharp cracks using the discrete element method. While the former is based on successive minimization of the total energy and cracks initiate as a result of the energy minimization, the latter finds fracture initiation based on the stress as the stress exceeds the strength defined between discrete elements.

2.4.1 Unsaturated Hydro-Mechanical Approach

The continuum formulation of the unsaturated HM-coupled problem outlined below is discussed in detail in Pitz et al. (2022) and is mainly based on a simplification of the approach from Grunwald et al. (2022). It is implemented in the current version of the open source code OpenGeoSys (Bilke et al. 2022).

The isothermal formulation is founded on the conservation of water volume in a deformable porous medium according to Biot’s consolidation theory (Biot 1941) considering unsaturated Darcy flow by the Richards equation (Richards 1931):

$$\begin{aligned} \begin{aligned} 0&= \rho _\text{ LR }S_\text{ L }\left( \phi \beta _{p\text{,LR }}+ S_\text{ L }(\alpha _\text{ B }-\phi ) \beta _{p\text{,SR }}\right) \frac{\text {d} {p_\text{ LR }}}{\text {d} t} - \,\textrm{div}\,\left( \rho _\text{ LR }\frac{{k^\text{ rel}_\text{ L }}{\textbf{k}_\text{ S }}}{\mu _\text{ LR }} \left( \,\textrm{grad}\,p_\text{ LR }- \rho _\text{ LR }\textbf{g} \right) \right) \\&+ \left[ \rho _\text{ LR }\phi + \rho _\text{ LR }p_\text{ LR }S_\text{ L }(\alpha _\text{ B }- \phi ) \beta _{p\text{,SR }}\right] \frac{\text {d} {S_\text{ L }}}{\text {d} t} + S_\text{ L }\rho _\text{ LR }\alpha _\text{ B }\,\textrm{div}\,\frac{\text {d} {\textbf{u}}}{\text {d} t}\,\\ \end{aligned} \end{aligned}$$

(2.2)

with the density of the liquid phase $\rho _\text{ LR }$, the liquid saturation $S_\text{ L }$, the porosity $\phi$, the liquid and solid compressibility $\beta _{p\text{,LR }}, \beta _{p\text{,SR }}$, respectively, the Biot coefficient $\alpha _\text{ B }$, the liquid pressure $p_\text{ LR }$, the intrinsic permeability tensor ${\textbf{k}_\text{ S }}$ and relative permeability ${k^\text{ rel}_\text{ L }}$, the liquid viscosity $\mu _\text{ LR }$ and the solid displacement vector $\textbf{u}$.

The capillary pressure $p_\text{ cap }$ is defined as the difference between the liquid pressure $p_\text{ LR }$ and the gas pressure $p_\text{ GR }$. Following the Richards equation, pressure gradients in the gas phase are neglected and by choosing the gas pressure to zero, we obtain following relationship

$$\begin{aligned} p_\text{ cap }= p_\text{ GR }- p_\text{ LR }= - p_\text{ LR }\, . \end{aligned}$$

(2.3)

The liquid saturation $S_\text{ L }$ is linked to the effective liquid saturation via $S_\text{ L }= S_\text {e}\left( S_\text {L,max}-S_\text {L,res} \right) + S_\text {L,res}$, which is in turn defined as a function of the capillary pressure according to van Genuchten (Van Genuchten 1980):

$$\begin{aligned} S_\text {e} = \left( 1+\left( \frac{p_\text{ cap }}{p_\text {e}} \right) ^n\right) ^{\frac{1}{n}-1}\,, \end{aligned}$$

(2.4)

with $p_\text {e}$ as the gas entry pressure and $n\in \left[ 0,1\right)$ as empirical material property. The relative permeability is defined according to van Genuchten/Mualem (Mualem 1976):

$$\begin{aligned} {k^\text{ rel}_\text{ L }}=\sqrt{S_\text {e}}\left( 1-\left( 1-S_\text {e}^{n/(n-1)}\right) \right) ^{(2n-2)/n}\,. \end{aligned}$$

(2.5)

The momentum balance equation is given by

$$\begin{aligned} \textbf{0} = \,\textrm{div}\,\left( \mathbf {\sigma }' - \alpha _\text{ B }S_\text{ L }p_\text{ LR }\mathbf {\textrm{I}}\right) + \left[ \left( 1- \phi \right) \rho _\text{ SR }+S_\text{ L }\phi \rho _\text{ LR }\right] \textbf{g}\,, \end{aligned}$$

(2.6)

where $\rho _\text{ SR }$ is the density of the solid phase, $\textbf{g}$ the earth’s gravitational acceleration vector and $\mathbf {\sigma }'$ represents the effective stress tensor, which can be calculated by the solid stiffness $\textbf{C}$ and the elastic strains: $\mathbf {\sigma }' = \mathcal {\textbf{C}} :\left( \mathbf {\epsilon }- \mathbf {\epsilon }_{\text {inelastic}} \right)$. In the current approach, inelastic strains are not taken into account, i.e. $\mathbf {\epsilon }_{\text {inelastic}} = 0$.

2.4.2 Extension with the Phase-Field Approach (Diffuse Cracks)

The coupled HM formulation is extended to account for cracking with the phase-field modeling approach for brittle fracture. With the phase-field approach (Bourdin et al. 2000, 2008), we introduce a scalar variable that represents the damage of the material $d:\Omega \rightarrow [0,1]$. A material, which is free from damage (intact) is represented by $d=0$, while a fully degraded material by $d=1$. Damage evolves in a diffuse manner, i.e. the crack phase-field variable evolves continuously, as shown on the right side of Fig. 2.11.

We couple the evolution of the crack phase-field d with the momentum balance equation by multiplying a degradation function g(d) based on the evolution of d with the effective stresses. Before doing so, we split the effective stresses into active $\mathbf {\sigma }'_+$ and non-active $\mathbf {\sigma }'_-$ parts computed from decomposing the elastic strain energy $\Psi$ (Amor et al. 2009; Miehe et al. 2010) as follows

$$\begin{aligned} \mathbf {\sigma }'= & {} g(d)\frac{\partial \Psi _+}{\partial \mathbf {\epsilon }} + \frac{\partial \Psi _-}{\partial \mathbf {\epsilon }} \nonumber \\= & {} g(d)\mathbf {\sigma }_+' + \mathbf {\sigma }_-' \end{aligned}$$

(2.7)

Note that only the active, crack-driving part is degraded by g(d). The momentum balance equation is rewritten as

$$\begin{aligned} \textbf{0} = \,\textrm{div}\,\left( g(d)\mathbf {\sigma }_+' + \mathbf {\sigma }_-' - \alpha _\text{ B }S_\text{ L }p_\text{ LR }\mathbf {\textrm{I}}\right) + \left[ \left( 1- \phi \right) \rho _\text{ SR }+S_\text{ L }\phi \rho _\text{ LR }\right] \textbf{g}\,. \end{aligned}$$

(2.8)

In this study, g(d) is given by

$$\begin{aligned} g(d) = (1-d)^2 + \eta , \end{aligned}$$

(2.9)

where $\eta$ is the residual stiffness (e.g., $\eta =1 \times 10^{-8}$ Pa) assigned in fully degraded part ($g(d)=0$) for numerical stability. An overview regarding the derivation and application of a similarly combined unsaturated HM and phase-field approaches is described in Cajuhi et al. (2018). The crack phase-field evolution is given by

Fig. 2.11

Discrete representation of fractures (left) and diffused phase-field representation of fractures (right)

$$\begin{aligned} g'(d)\Psi _+ ({\textbf{u}}) +\frac{\mathcal{G}_c}{4 c_w} \left( \frac{w'(d)}{\ell } + 2\ell \Delta d \right) = 0 , \end{aligned}$$

(2.10)

where

$$\begin{aligned} w(d) := (1-d)^n \quad \text {and} \quad c_w :=\int _0^1 (1-s)^{n/2} \, \textrm{d} s , \end{aligned}$$

(2.11)

and $n=1$ and 2 stand for so-called AT$_1$ and AT$_2$ models respectively (Pham et al. 2011; Bourdin et al. 2014), and $c_w = 2/3$ for AT$_1$ and 1/2 for AT$_2$.

In the phase-field literature, the split of the energy is sometimes referred to anisotropic response (Ambati et al. 2015). In this book, the term anisotropy refers to the bedding orientation and constitutive behavior. For referencing the energy decomposition (split), we use the term tension-compression asymmetry. There are a few popular split models in the literature, nevertheless, most of them can only be used for isotropic materials. As for anisotropic materials, the principal directions of the stress and strain no longer coincide. Hence, the orthogonality condition, a requirement for the split method to be variationally consistent, is not guaranteed. Thus, further treatment is required. In our recent work, we generalized two existing split models for isotropic materials in order to account for materials with anisotropic constitutive behavior. Borrowing some basic concepts from He and Shao (2019), we performed an energy preserving transformation so that the orthogonality condition is satisfied, and thus, a decomposition of an arbitrary anisotropic elastic model is feasible. We encourage those interested readers to follow our manuscript for a detailed clarification (Ziaei-Rad et al. 2023). Exemplarily, we show the verification of the framework in Sect. 2.5.1.2.

2.4.3 Extension with the Discrete Element Method (Discrete Cracks)

The coupled HM formulation is extended to account for cracking with the Finite Discrete Element Method (FDEM).

2.4.3.1 FDEM Mechanical Solver

Originally proposed by Munjiza (2004), the Finite Discrete Element Method (FDEM) utilizes principles of continuous numerical approaches in combination with those of discontinuous methods to simulate interaction of multiple deformable bodies (Lisjak et al. 2014). In other words, while the elastic deformation of discrete bodies and the initiation and propagation of fractures are captured using Finite Element analysis, the interaction between individual blocks of elements are realized by Discrete Element analysis. In FDEM, each block of the problem domain is discretized into a finite element mesh with zero-thickness joint elements being inserted on its intersection (common) edges of the adjacent elements as shown in Fig. 2.12. Then, an explicit finite difference time integration scheme is applied to solve the equation of motion for the discretized system which reads as

$$\begin{aligned} \textbf{M}\frac{\partial ^2 \textbf{x}}{\partial t^2} + \textbf{C}\frac{\partial \textbf{x}}{\partial t} + \textbf{F}_\text {int}(\textbf{x}) - \textbf{F}_\text {ext}(\textbf{x}) -\textbf{F}_c(\textbf{x}) = 0, \end{aligned}$$

(2.12)

where $\textbf{M}$ and $\textbf{C}$ denote the diagonal lumped mass and damping matrices, respectively. $\textbf{x}$ is the nodal displacement vector, and $F_\text {int}$, $F_\text {ext}$ and $F_c$ are the nodal internal, external and contact forces, respectively. Contact forces, $F_c$, are calculated between contacting discrete blocks either in the normal direction, using a penalty method, or in the tangential direction, employing Coulomb-type friction law. Accordingly, internal resisting forces, $F_\text {int}$, are calculated as summation of elastic forces and joint elements’ bonding forces.

Fig. 2.12

a FDEM discretization of the intact discrete blocks containing triangular finite elements and zero-thickness joint elements, b dual grid of FDEM flow solver

Following a cohesive-zone approach, the mixed-mode crack model consists of a typical stress-strain curve where the bonding stresses are proportional to the separation and the sliding of the fracture edges. More details of the crack model can be found in the literature (Mahabadi et al. 2012).

2.4.3.2 FDEM Flow Solver

In order to account for the saturated/unsaturated porous media flow within the framework of FDEM, a vertex-centered finite volume scheme is adopted. Based on this approach, the conservation of the fluid phase is enforced on a dual grid of polygonal control volumes. As illustrated in Fig. 2.12, each control volume is generated around vertices of the mechanical solver grid connecting centers of triangular elements through midpoints of primary grid. The adopted mass balance equation of fluid phase in the porous medium is a simplified version of Eq. 2.2 with incompressible solid phase. The hydro-mechanical coupling in FDEM is realized by a staggered approach alternating between solid and flow solvers. Within this approach, the valuation of the stress field is affected by the pore pressure according to the effective stress principle, while the flow calculations are affected by the volume change associated with each computational node.

2.5 Process Simulation Results

We show different simulation setups and their results using the methods described in Sect. 2.4. At first, we model the mechanical processes at laboratory scale considering the anisotropy of the Opalinus Clay and using the FDEM formulation (Sect. 2.4.3). The response of the samples is compared with the laboratory experiments discussed in Sect. 2.3.1. Cracks are modeled either discretely or diffusely and we compare both approaches. At the field scale, the deformation and pore pressure behaviors are described with the continuum approach from Sect. 2.4.1 and the hydro-mechanical effects such as the shrinkage due to the de-saturation are discussed. We model the formation of the excavation damaged zone (EDZ) with discrete cracks (formulation from Sect. 2.4.3). Finally, desiccation cracks in the field are modeled with the combined poromechanical phase-field approach (formulation from Sect. 2.4.2).

2.5.1 Laboratory Scale

2.5.1.1 Crack Paths Due to Deformation

As indicated by the achieved experimental results in Sect. 2.3.1.1, not only the deformation response but also the strength characteristics of Opalinus clay are strongly influenced by the bedded structure of the samples. In the 2D FDEM method, the anisotropic elastic behavior is captured by a transversely isotropic stress-strain law, while the strength directionality is modeled adopting a so-called smeared approach (Lisjak et al. 2014). Within this approach, the strength of each joint element is directly linked to its relative orientation with respect to the bedding layers. In this way, the strength parameters and fracture energies fluctuate linearly between the minimum limiting value, when the cohesive element is parallel to the bedding layer, and the maximum limiting value, when the cohesive element is perpendicular to the bedding layer.

Fig. 2.13

a Numerical model setup of SCB experiment and mesh details for b P-sample c Z-sample and d S-sample

Table 2.3 Material and input parameters for FDEM simulations of three point bending test of OPA

The FDEM models were calibrated to match the obtained experimental peak loads reported in Sect. 2.3.1.1. The laboratory-scale model along with the employed boundary conditions as well as the mesh details for P-, S- and Z-samples are illustrated in Fig. 2.13. The experimental loading conditions were realized by explicitly modeling the supports and the loading platen. In order to avoid the impact effect of loading platen, the loading velocity increases linearly from zero to 0.1 m/s in 1 ms. Following the procedure reported in the work of Tatone and Grasselli (2015), the calibrated mechanical material properties in the FDEM model are summarized in Table 2.3.

Figure 2.14 illustrates the numerical simulation results concerning the fracture path trajectories upon failure and the distribution of horizontal stress for the three cases. The crack propagation shows good agreement with the experimental observations depicted in Fig. 2.5.

Fig. 2.14

Fracture path trajectory (top) and corresponding distribution of horizontal stress [Pa] (bottom) at 0.5 mm of crack mouth opening displacement for Opalinus clay of P-, Z- and S-samples with a $\theta = 0^\circ$, b $\theta = 45^\circ$ and c $\theta = 90^\circ$, respectively

2.5.1.2 Pure Shear Test

In order to verify the anisotropic phase-field model detailed in Ziaei-Rad et al. (2023), we present a set of numerical examples of a single edge notched shear test. The geometry and boundary conditions for this examples are shown in Fig. 2.15. We consider a square plate ($0.001\,\text {m} \times 0.001\,\text {m}$) under plane strain conditions with a pre-notched crack at the middle height of the specimen. The length scale parameter was chosen as $\ell = 1\times 10^{-5}$ m which is small enough with respect to the specimen dimensions. The AT$_1$ model was used in the simulations. We used the parameters of Opalinus clay listed in Table 2.3. Note that the critical surface energy $G_c=100$ N/m was taken to be independent of the material orientation.

Fig. 2.15

Schematic of a cracked square plate under a displacement shear loading. An incremental displacement loading is applied on the top edge. The material principal axes ($\boldsymbol{e}_1,\boldsymbol{e}_2$) have an angle of $\alpha$ with the global coordinate system ($\boldsymbol{x},\boldsymbol{y}$)

As boundary condition, we apply fixed displacement (zero in all directions) at the bottom edge ($y=-L/2$). For simulating shear, we prescribe a displacement along the x-direction at the top edge ($y = L/2$). The computation was performed with a monotonic displacement controlled loading with a constant displacement increment $\Delta \boldsymbol{u} = 3\times 10^{-7}$ m.

Figure 2.16 depicts the crack paths for orthotropic volumetric-deviatoric (Amor et al. 2009) split under shear loading with three material orientation angles ($\alpha =-\pi /4$, 0, and $\pi /4$). Herein a significant difference in crack response is observed between the three cases.

Fig. 2.16

Cracked square plate under shear. Phase field contours of orthotropic volumetric-deviatoric for three material orientations. The same parameters were used for all cases. Here, we observe that the crack paths are strongly dependent on the material orientation. With $\alpha =\pi /4$ there exist straight cracks that grow towards the lower right corner, almost along the material principal axis with lower stiffness. The two other cases differ from an expected isotropic response

The volumetric-deviatoric decomposition assumes the volumetric compression does not contribute to the crack evolution. However, in case of anisotropy, the Poisson effect may pronounce more in a direction compared to the others, thus, even under a uniaxial compression, some crack driving forces can lead into damage development. Consequently, the crack response of a compression-dominated loading is depending on differences in directional stiffness of material. An implementation of the presented anisotropic phase-field formulation combined with the unsaturated HM approach is underway.

2.5.2 Field Scale

2.5.2.1 Hydro-Mechanical Effects and Model Response

In this section, we numerically model the hydro-mechanical coupled effects in the context of the CD-A experiment using the macroscopic approach presented in Sect. 2.4.1. The two-dimensional model consists of a cross section of the open and closed niches and surrounding domain, Fig. 2.17. We prepare an isotropic and an anisotropic model, in which the bedding direction is assumed to be approximately horizontal as monitored in the field (Ziefle et al. 2021). The OPA parameters and properties are extracted from several references in the literature (Bock 2000; Martin and Lanyon 2003; Garitte et al. 2017; Bossart et al. 2018; Zhang et al. 2019) and presented in Table 2.4. The transversal component of an anisotropic property is defined with the subscript 1, while its parallel components are defined with subscripts 2 and 3. These subscripts corresponds to a local coordinate system oriented according to the bedding layers, i.e. 1 is orthogonal to the bedding. We compute the isotropic equivalent properties as the average of the anisotropic properties, e.g. the isotropic equivalent elastic modulus corresponds to the average of its parallel and perpendicular components. The isotropic equivalent properties are denoted with the subscript $(.)_\text {iso}$ in Table 2.4. The initial stress values are based on in-situ stress measurements (Martin and Lanyon 2003). Their orientation correspond to the global coordinate system, i.e. in two-dimensional setting their horizontal and vertical components $\sigma _\text {h}$ and $\sigma _\text {v}$ lie on the global coordinate system (x, y) (Fig. 2.17), while in three-dimensional setting $\sigma _\text {h}, \sigma _\text {H}$ and $\sigma _\text {v}$ on (x, y, z) (Fig. 2.21). The initial pore pressure and in-situ stresses are applied on the domain. After instantaneous excavation, the seasonal boundary condition computed from the RH measurements (Fig. 2.1) with the Kelvin’s equation (Eq. 2.1) is applied at the niches. In the closed niche, the drops in the RH curve are related to installations and openings of the sealing door. We consider that the latter remains under high RH, i.e. neglect the drops, and assume a constant pore pressure of ${-}1$ MPa. The open niche follows a cyclic pressure evolution, interpolated from the first season of the RH curve, Fig. 2.1.

Table 2.4 Parameters and material properties used for simulating the CD-A experiment

Fig. 2.17

Setup of the two-dimensional model of the CD-A experiment consisting of open and closed niches. The niche radius is 1.15 m and distance between them is 10 m. The bedding direction is horizontal and indicated in the figure

Figure 2.18 shows the saturation results using the anisotropic model for March 2020 (first desaturation peak in winter) Fig. 2.18a and September 2020 (first re-saturation peak in summer) Fig. 2.18b. Note that the horizontal direction is more permeable and, consequently, the de-saturated zone is more extended in this direction when compared to its orthogonal counterpart. While the saturation in the closed niche remains high, the saturation near the wall of the open niche increases between March and September. At the same time, the overall desaturated zone extends along the domain, but still indicates a local effect. Except for direction dependency, the same conclusions are drawn in the isotropic model.

Fig. 2.18

Saturation changes in deformation results in the 2D model of the CD-A experiment. The region near the open niche dries, de-saturates, during winter and re-saturates during the summer months. Convergence is observed due to an increase of the deformation near the niches. Note that both effects are more pronounced in the open niche in comparison to the closed niche

The aforementioned tendencies, i.e. de-saturation and re-saturation during winter and summer, respectively, have been compared with the monitored water content and electrical resistivity results. A decrease in saturation has been associated to a decrease of water content and an increase of the resistivity, while the opposite occurs when the saturation increases (Amabile et al. 2020; de Jong et al. 2020; Cajuhi et al. 2021). These comparisons, including a geologic characterization of the niches have revealed that inhomogeneities like cracks, shear zones, carbonate-rich and sandy lenses subfacies in the claystone have a strong influence on the seasonal desaturation or re-saturation behavior in the near field. A detailed description of this study and its conclusions are reported in Ziefle et al. (2022a, b).

The numerical results show that the seasonal saturation changes affect the deformation behavior. We compare the strain response of the isotropic and anisotropic models, Figs. 2.19 and 2.20, respectively, in the first winter and summer peaks. On both niches, the contribution of the strains due to excavation is equal. The increase of deformation represents convergence of the niches. Note that the projected initial stress is the same in both isotropic and anisotropic models. In the isotropic model, the strain responses are similar, while the strains in vertical direction are larger than in the horizontal direction in the anisotropic model. Note that the strains decrease from winter to summer when comparing the responses of both models. This decrease is specially observed in the vertical and shear components, and particularly, in the open niche, where seasonal effects are more pronounced and accounted for in the model.

Fig. 2.19

Deformation response of isotropic model. Blue regions denote negative strains, while red ones relate to positive strains. The strain during the first desaturation peak (March 2020) are shown on the left, while the strains during the second peak are shown on the right (September 2020). An equal strain response is observed near the niches due to excavation. The change of saturation induces slightly larger deformation in the vicinity of the open niche

Fig. 2.20

Deformation response of anisotropic model. Blue regions denote negative strains, while red ones relate to positive strains. The strains due to the first desaturation peak (March 2020) are shown on the left, while the strains due to the second peak are shown on the right (September 2020). An equal strain response is observed near the niches due to excavation. The change of saturation induces larger deformation in the vicinity of the open niche. Due to anisotropic elastic properties, the strains in vertical direction are larger than in the horizontal direction

As the driving mechanism of the crack phase-field is the active elastic energy, which evolves with the strains (see Sect. 2.4.1), we expect to model cracks due to drying in the open niche using this approach. Desiccation cracks are expected to occur due to mode I loading, i.e. tension cracks. In the field, they were observed on the wall surfaces, i.e. where the rock is exposed and consequently drier in comparison to the undisturbed rock. Such cracks are expected to propagate from the niche walls into the depth of the rock. To further understand the stress distribution on the wall surfaces, we have prepared a 3D model of the experiment. For completeness, the model also considers the lock to the closed niche and the gallery that provides access to the niches. The seasonal boundary conditions are applied in the open niche, lock and gallery. The numerical model is shown in Fig. 2.21. The numerical model has been also utilized for interactive visualization purposes in the virtual laboratory prototype developed by Graebling et al. (2022a) and briefly discussed in Sect. 2.6. As aforementioned, the bedding direction follows the horizontal direction of the niches’ plane. With respect to the gallery, the bedding shows an average orientation of (147/57) degrees (Jaeggi et al. 2020; Ziefle et al. 2021) as pointed out in Fig. 2.21.

Fig. 2.21

Setup of the three-dimensional model of the CD-A experiment. In-situ pore pressure and stresses are applied as initial conditions. The closed niche is subjected to a constant pore pressure, while the latter varies seasonally in the open niche

For a comprehensive and comparative visualization, we plot the saturation behavior in the model, Fig. 2.22. The upper part of the model is not shown in the visualization. Similar to the two-dimensional model, the saturation decreased in March 2020. A saturation decrease is observed in the open niche, lock and gallery as a result of the pressure changes. Due to the application of a constant pore pressure (boundary condition) in the closed niche, the saturation remains very high, i.e approximately 1. Note that the extent of the de-saturated zone in the vicinity of the gallery, lock and open niche remains approximately constant in size. An excavation damaged zone is not set explicitly in the model.

Fig. 2.22

Saturation response in the three-dimensional model of the CD-A experiment. The saturation increases from March 2020 a to September 2020 b in the open niche and gallery

The saturation changes and a comparison between the modeled deformation and monitored convergence data is found in Ziefle et al. (2022c). As also observed in the 2D setting, the strain response of the model is affected by the excavation and by the seasonally change of air humidity. The strain response during the first de-saturation peak in March 2020 and re-saturation in September 2020 are shown in Fig. 2.23. Overall, all shown strain components increase up to winter and decrease in summer as observed, for example, in the strain component parallel to the axial direction of the gallery, Fig. 2.23a, b, and shear component, Fig. 2.23e, f. In the regions were the gallery and niches intersect, there is an increased influence of the saturation changes happening in the gallery. This can be visualized in all shown strain components, but specially in the direction parallel to the niches’ axis as depicted in Fig. 2.23c, d. Note that there is a development of positive strains (stretching) on the walls of the niches and gallery. These contribute to an increase of the energy required for crack propagation and contributes to the fact that cracks are prone to develop on the walls of the niches.

Fig. 2.23

Deformation response of the three-dimensional, anisotropic model. The results concerning the first desaturation peak (March 2020) are shown on the left, while the ones during the first de-saturation peak are shown on the right (September 2020). The strains are positive at the gallery and niches’ walls. Note that the strains are larger in the open niche

The described model behavior and mechanisms serve as a basis for model extension with respect to cracking. The proposed, extended model (Sect. 2.4.1) is first applied in a laboratory scale test related to the desiccation of a slab, which has been chosen due to the similarities with CD-A. Several tests concerning the energy splits, model formulation and understanding the influence of the numerical parameters involved in the coupled Richard mechanics phase-field approach have been carried out in Cajuhi et al. (2022). In the latter, the calibration and applicability of the coupled framework in the field scale are described in detail. The implemented framework is available in the open source software OpenGeoSys (Bilke et al. 2019, 2022). A general application of the framework in three-dimensional setting using a homogeneous geological window of the open niche and in context of desiccation cracking is discussed in Sect. 2.5.2.3. Before entering the fully coupled problem, we model and discuss the cracks initiated during the excavation using a FDEM approach in the following.

2.5.2.2 Excavation Damaged Zone

The Excavation Damaged Zone (EDZ) is considered as one of the main factors compromising the integrity in the near field of underground repositories. In the framework of FDEM, simulation of EDZ formation is captured with two separate model runs (Lisjak et al. 2014). In the first run, the in-situ stress state is applied by translating the specified stresses into nodal forces. After reaching the steady state, the displacement of outer boundaries are fixed allowing the model to maintain the in-situ stress conditions (Mahabadi et al. 2012). The excavation step would then be initiated employing a core replacement technique (Kavvadas 2005). Within this approach, the Young’s modulus of the core material, i.e., the excavation region, is reduced in a step-wise manner to mimic the gradual de-confinement of the rock mass resulting from the tunnel advancement.

Due to the heterogeneous geology of the Opalinus clay rock, we evaluate the influence of the bedding layer orientation in forming an EDZ. In the first case, as shown in Fig. 2.24, the numerical model setup includes a 15 m $\times$ 15 m square domain with the opening region at its center with internal radius of 1.15 m. The geological window is representative of a one-quarter model of CD-A twins in which the bedding direction is horizontal. Additionally, the EDZ formation is studied with bedding direction inclined to tunnel axis at 45 degrees. In the second case, due to the asymmetric nature of the problem, a full 30 m $\times$ 30 m square domain has been considered. The employed material properties for transversely isotropic Opalinus clay is assumed same as the ones reported in Table 2.3. For both numerical test cases the vertical and horizontal in-situ stress are equal to 6.5 MPa and 4.5 MPa, respectively.

Fig. 2.24

Two-dimensional FDEM geometry of CD-A experiment for modeling the excavation damaged zone

Figures 2.25 and 2.26 show the gradual increase of the failure pattern at different ratios of deconfinement, i.e., $\lambda$. Deconfinement is defined as the ratio between the decreased stress magnitude at the vicinity of the tunnel face to the initial magnitude of in-situ stress. Note that the crack propagation is parallel the bedding orientation in both cases. At the end of the simulation, we notice that the crack network formed with the horizontal orientation is more complex and connected in comparison to the case where the bedding is inclined with $45^\circ$. On the other hand, the fracture network extends deeper into the formation for the second case. The percentage of broken joints as a function of distance from tunnel face are plotted for both cases at $\lambda = 0.6$ in Fig. 2.27a. The steeper curve represents the $45^\circ$. Moreover, in Fig. 2.27b we show how the distribution of broken joints can be used in the creation of a model with permeability values that depend on the distance from the niche wall, i.e., the near the wall the more permeable the model, since more joints are broken.

Fig. 2.25

FDEM results of EDZ failure pattern near the open niche with horizontal bedding layer at a $\lambda = 0.25$, b $\lambda = 0.75$

Fig. 2.26

FDEM results of EDZ failure pattern near the open niche with $45^circ$ bedding layer at a $\lambda = 0.3$, b $\lambda = 0.6$

In other words, the complex crack network induced by the excavation can represent a zone of more instability and larger permeability in comparison to the intact rock. To further understand the latter, hydraulic models need to be taken into account and coupled with the mechanical formulation. Ongoing work is being performed to account for hydraulic effects in the FDEM approach.

Fig. 2.27

a Distribution of broken joints percentage around underground opening. The percentage of broken joints in the model with horizontal bedding direction is larger (y-axis) than that in the $45^\circ$ model. In the former, the cracks are less deep (x-axis) than those generated in the inclined model. Crack development motivates the EDZ abstraction concept and altered permeability profile shown in (b)

2.5.2.3 Desiccation Cracks

In this subsection, we use the fundamental equations from the continuum based methods (Sect. 2.4) to simulate desiccation cracks in three-dimensional setting. For the latter, no changes in the phase-field model are required, i.e. the same evolution equation is utilized and there is no need to adapt the algorithm to initiate the cracks or due to branching. The verification of the framework is based on the setup discussed in Cajuhi et al. (2018).

Since the open niche shows more pronounced hydro-mechanical effects, we propose a quarter three-dimensional model of a so called homogeneous geological window as shown in Fig. 2.28. The thickness of the model is the same as the radius of the niche (1.15 m). The total height/width is 5 m. A drying pore pressure related to the RH drop during the first winter is applied at the niche wall. At the remaining boundaries, no flow and symmetry boundary conditions are applied. Since the (free) movement of the niche wall, imposed by mechanical boundary conditions, can affect crack initiation and propagation and these conditions are not well defined in the field, we propose a free and a restrained setting. In the latter, the movement of the nodes at the front and back surfaces of the niche is hindered in the direction normal to the plane.

Fig. 2.28

Quarter model of the three-dimensional model of the CD-A experiment (geological window). In both free and restrained settings, symmetry boundary conditions are applied. In the restrained case, the movement of the niche wall is hindered by restraining the normal displacement of the nodes

In the following, the results concerning the degree of saturation and cracking are shown. For visualization purposes, a threshold filter is applied on the crack phase-field variable so that the elements where damage values equal or greater than $99\%$ are removed from the blended out from the mesh. When the wall is allowed to move freely, a single crack initiates and propagates along as shown in Fig. 2.29a. The crack initiates at a pore pressure corresponding to approximately $94\%$ RH and 80–81$\%$ saturation. As expected, restraining the front and back surfaces, and consequently the movement of the wall, induces earlier crack propagation and a more complex crack pattern. In the latter, the crack initiates approximately at $97\%$ RH and 90–92$\%$ saturation. Both cases lie in the expected drying RH range (Fig. 2.1), also reported in Regard et al. (2021), where crack development at the niche walls has been monitored and quantified.

Fig. 2.29

Joint visualization of saturation and crack phase-field responses of quarter model of the open niche. The model desaturates near the wall of the open niche (red color corresponds to drier conditions, while blue color to wetter conditions). For a direct visualization, the elements with damage equal or greater than $99\%$ are not plotted, i.e. the empty spaces at the wall correspond to cracks. A single desiccation crack initiates at the wall of the open niche in the free setting. In the restrained setting, a more complex crack pattern develops

The current phase-field evolution equation is not able to account for the mechanical and hydraulic anisotropies of the material. For this reason, the developed cracks do not necessarily propagate parallel to the expected bedding orientation as observed experimentally. Ongoing work is being carried out to extend the numerical approach with to anisotropic setting. Basis for the extension is the framework proposed in Ziaei-Rad et al. (2023) and discussed in Sect. 2.4. This extension should affect the crack orientation and induce cracks parallel to the bedding.

Additionally, the numerical simulation of desiccation cracking has been addressed in the framework of FDEM with coupled RM approach. The problem domain includes a 15 m $\times$ 15 m two-dimensional quarter model with 1.15 m niche radius. The assumed mechanical properties are the same as the ones considered in previous FDEM models dedicated to numerical simulation of fracturing in OPA. Fluid properties, porosity, values of permeability matrix and van Genuchten model parameters are all taken from Table 2.4. As mentioned in Sect. 2.4.3.2, the volume conservation equation and the mechanical equilibrium both are solved with finite difference explicit time integration scheme. This, consequently, forces smaller time steps for stability and convergence of the numerical solution. In the FDEM model, we impose a drying flux of $q=6.0 \times 10^{-9}$ m/s that acts on the niche wall for 24 h. The failure pattern as well as the saturation distribution after 24 h is shown in Fig. 2.30.

Fig. 2.30

FDEM results of desiccation fracturing around the open niche after a 12 h and, b 24 h simulation time

As can be seen, unsaturated areas, as a result of desiccation boundary condition, are only limited to the surrounding of the niche because of low permeability of Opalinus clay. Cracking starts from the top of the niche after approximately three hours of simulation and at approximately 3.5 MPa suction (RH between 97–98$\%$).

2.6 Data Analytics and Visualization

A prototype of the CD-A experiment using the 3D model presented in Sect. 2.5.2.1 has been prepared and discussed in detail in Graebling et al. (2022a). The model is placed in the Mont Terri rock laboratory and its results can be accessed via the proposed interactive tool. Furthermore, experimental data and its borehole positions can be combined and visualized together with the numerical model. In Graebling et al. (2022a), convergence measurements obtained from the boreholes via extensometers have been compared with the strains computed with the numerical model. The direct visualization and comparison of the data, combined with the experience to be able to understand where the data come from, i.e. by visualizing the borehole positions and the position of the experiment itself, increases the possibilities for exchange between different fields of expertise and with that, a multidisciplinary approach. Further details on the approach and examples of the visualization are discussed in Chap. 5.

2.7 Summary, Conclusions and Research Transfer

2.7.1 Summary

This chapter presented a workflow for understanding the Opalinus Clay characteristics, especially with respect to its shrinkage behavior. We carried out several laboratory experiments concerning the mechanical characterization of the rock. The results have shown and confirmed the strong influence of the anisotropy on the fracture behavior. In the SCB samples where the initial crack is parallel to the bedding, the crack propagates more easily, at a lower peak load, when compared to the experiments performed in the samples with initial crack perpendicular to the bedding. The tests have also shown the combined interaction of bedding and heterogeneities in the OPA since the cracks do not propagate uniformly in all samples. We also computed these experiments using FDEM and phase-field models coupled with the mechanics. Experience has also been gathered with respect to sample preparation and cutting procedures. For increasing the understanding on the effect of humidity variations, tests under different humidity conditions have been carried out. The tests were performed in a customized desiccator cell which has been equipped with different sensors to measure the controlled RH and temperature. In summary, the experimental tests show that an increase of suction (decrease of RH) leads to an increase of the strength of the material—as long as no fractures are induced that lead to an overall weakening of the material. The obtained data can be used as input in the simulations of drying/wetting and will affect the fracture behavior of the material in coupled hydro-mechanical simulations.

The laboratory tests have a direct connection with the field, since the samples were obtained from drilled boreholes in the vicinity of the in-situ CD-A experiment. In this region, the main lithology is characterized by the sandy facies of the OPA. We presented and used material properties obtained from the literature and from the experiments performed in laboratory. These data have been used as input in the various presented numerical models. The monitored data at the CD-A twin niches have been used for comparison with water content and convergence data. These results have been described briefly and are object of further investigations (Ziefle et al. 2022b, c). The twin niches at the CD-A experiment provide a good platform for comparing the hydro-mechanical effects under seasonal (open niche) and controlled high RH conditions (closed niche). Furthermore, the niches are equipped with several monitoring sensors, e.g. extensometers, permeability, water content, crackmeter, among others. The open niche shows more pronounced shrinkage effects in comparison to the closed niche, for example with respect to desiccation crack development.

The development of excavated damaged zones and further hydro-mechanical effects near the excavations, for example the initiation and propagation of desiccation cracks during ventilation, are motivating facts for developing and extending numerical approaches with respect to cracking and evaluating their applicability in the in-situ scale. With this objective in mind, we proposed continuum and discontinuum based models. A comparative study between these approaches for the laboratory scale is shown in Sattari et al. (2022). We described a continuum approach for modeling the unsaturated HM-coupled problem that has been implemented in Bilke et al. (2019, 2022). This approach has been applied in two- and three-dimensional settings. Aim of the simulations was to understand the effects of de-saturation and re-saturation in isotropic and anisotropic conditions for the hydraulic (permeability) and the mechanic (elasticity properties). The isotropic properties have been computed as an average of the anisotropic ones. We have shown the results for the first de-saturation peak, i.e. largest suction, and for the first re-saturation peak, i.e. largest pore pressure. The de-saturated zone increases up to the smallest pressure and part of it remains up to the re-saturation peak. We have observed the influence of the saturation changes on the deformation behavior and have stated that these are more pronounced in the open niche, where the seasonal behavior plays an important role and desaturation is stronger. These conclusions are valid in both isotropic and anisotropic models. However, in the latter, the desaturated zone extends further due to larger permeability along the horizontal direction.

The numerical results concerning the hydro-mechanical effects in 2D setting have shown that the contour and vicinity of the open niche are drier in comparison to the domain as a whole. This represents a driving mechanism for desiccation cracks. To obtain further information on the drying behavior of the walls, the model has been extended to 3D. This model is composed of open niche, closed niche, lock to the closed niche and gallery. We have observed an increase of the strains on the walls of the niches and gallery, with more pronounced effects in the latter and in the open niche.

Note that the models discussed up to now have not considered the development and extent of an EDZ. For modeling the crack networks formed during the excavation, we have used an FDEM approach. The approach combines the finite element analysis with a discontinuous representation of the crack, i.e. through zero-thickness joint elements that are inserted in the mesh. We observed a more complex crack network when modeling an horizontal bedding layer on the cross section of the niches that might represent a zone of weaker stability and larger permeability.

The unsaturated HM approach has been extended to account for cracks due to drying. The model has been applied into an homogeneous geological window, consisting of a quarter of the open niche. The free and restrained niche wall movement have been simulated to account for the full range of boundary conditions in the field. In the former, a single crack has been formed, while the latter induced a more complex crack pattern. The crack propagation is not direction dependency The effect of the hydraulic and mechanical anisotropies have not been taken into account in this model. The obtained crack patterns do not follow the bedding direction since the model is isotropic. Nevertheless, the quarter model has been able to predict desiccation cracks in the expected and monitored RH range (94–97$\%$).

Finally, an interactive visualization of data in the field proportions has been proposed. The 3D model presented in this chapter has been inserted in this virtual prototype of the Mont Terri Laboratory, where the CD-A experiment takes place. The numerical results of saturation and deformation can be rendered and directly compared with borehole data, its position and monitored data provided by the sensors.

2.7.2 Conclusions and Research Transfer

We have proposed a workflow consisting of laboratory, field scale experiment and visualization of data. The laboratory experiments have provided new data related to the mechanical and fracture behaviors of the sandy facies of the OPA. Furthermore, a relationship between peak load and suction has been established for this material. The data obtained experimentally can be used as input for coupled hydro-mechanical simulations. For that, e.g. the strength needs to be determined for each tested suction range. An important model improvement would be to use strength-suction dependent material properties for simulating the material in laboratory and field scales.

The FDEM method showed good agreement with the measured laboratory experiments and has been also tested in the field scale for dealing with the EDZ. The numerically obtained EDZ shows complex crack networks, which might have an important effect on the permeability changes near the excavations. For further understanding this behavior, the hydraulic effects need to be coupled with the fracture-mechanical approach. This is object of ongoing work. A comparison between the response of the coupled hydro-mechanical methods (continuum and discontinuum based) is underway (Sattari et al. 2022). With that, we will evaluate the computational costs related to the inclusion of the discrete elements (joints) and of the addition of the phase-field evolution equation in each formulation.

The extension of the continuum framework based on the unsaturated HM approach using the phase-field model for brittle fracturing has the advantage and ability of simulating desiccation cracks. For a more resolved crack, the meshes need to be refined and this fact is directly connected to increased computational effort. Methods for automatic mesh refinement need to be further developed and explored in the context of HM coupled simulations. In the context of the unsaturated HM approach, the use of the presented anisotropic phase-field formulation is being evaluated and implemented in the open source framework. Furthermore, ongoing research is planned to account for swelling effects and crack closing as well as cracks due to swelling.

The long-term planning and execution of the CD-A experiment is an important aspect for testing the developed/extended modeling approaches and comparing its results to field data. Furthermore, the interactive data visualization of in-situ data is an important tool for communication and data exploration. The prototype presented in Graebling et al. (2022a) can be extended to account for further boreholes and sensors.

An interesting aspect for further research is the application of the presented methodology on other types of clays and rocks. This methodology can aid on finding the links between the different spatial and temporal scales. The combination of laboratory and field experiments, fundamental and applied research as well as interactive visualization represents an added value for different applications such as those in the geo-energy and radioactive waste disposal.