Numerical Benchmarking of GREAT Cell Experiments: Poly-Axial Stress Effects on Fluid Flow in Fractured Rock Using Smeared and Discrete Methods

1. Introduction

DECOVALEX is a long-standing international benchmarking project focusing on the systematic improvement of models for repository research [1,2,3,4,5]. Multi-physical/chemical processes, so-called thermo-hydro-mechanical-chemical (THMC) processes, which can be used to describe the temporal and spatial development of repository systems in terms of continuum mechanics, play a special role here. Because of the international participation, different host rocks have been addressed. Aside from clay and salt rocks, crystalline rocks have also been investigated [6,7,8] and are the focus of the present study. Benchmarking based on analytical solutions and/or code comparisons with academic, synthetic test examples has been broadly used to check the accuracy of numerical models and/or to test the correct implementation of numerical methods [9]. For DECOVALEX tasks, geotechnical experiments from underground laboratories are systematically selected to validate the numerical models using measured data. The transferability of the models from laboratory to in-situ scale plays a central role (experimental analysis) [10,11,12,13,14,15]. By validation, we mean the demonstration of physical adequacy of a model for a well-specified intent or context. This can be assessed via the ability of models to predict measurement results that were not used for model calibration (blind prediction, class A prediction). In recent DECOVALEX projects, tasks for the characterisation of repository systems in a geological context have also been defined and processed (performance assessment) [16]. The analysis of complete repository systems requires coping with the enormous computational effort. Therefore, methods of efficient computing (i.e. parallel computing) are becoming increasingly important for comprehensive safety analysis. Alternatively, methods are also being developed to find suitable surrogate models for the complex coupled THMC processes, using machine learning methods [17,18,19,20].

Despite the long tradition of DECOVALEX, the SAFENET Task brought fracture mechanics into the project for the first time in a systematic benchmarking during DECOVALEX 2023 [21]. Initially, model comparisons were carried out with defined analytical solutions [22], and mechanical (M), hydro-mechanical (HM) and thermo-mechanical (TM) experiments were carried out in the rock laboratories of the TU Bergakademie Freiberg [23], the University of Edinburgh [24,25,26] and KICT [27,28,29], respectively. In the present work, new test examples for fully coupled HM models will be developed on the basis of the GREAT (Geo-Reservoir Experimental Analogue Technology) cell laboratory experiments, which will be used for the experimental planning of corresponding experiments in the new DECOVALEX project 2027 (SAFENET 2).

The GREAT cell of the University of Edinburgh is a unique experimental facility for conducting triaxial test in a rotational stress field. This setup allows the measurement of fracture permeabilities under varying stress conditions. However, the complexity of this equipment makes it challenging to reproduce its results numerically. To overcome this, we developed a simplified benchmark version suited for numerical modeling. For the first time, numerical benchmarks will be proactively designed based on simplified examples derived from GREAT cell experiments (with laboratory replication planned for future work). This novel approach is addressing the ongoing challenge of accurately modeling fracture mechanics processes under varying pore pressures and mechanical loads (HM conditions).

Numerical modeling fracture mechanics under varying HM conditions is challenging, requiring precise numerical methods to replicate fracture behaviors observed in experiments like the GREAT cell. Numerical methods for fracture modeling are broadly classified into sharp and smeared fracture approaches [30]. Sharp fracture methods, such as the Extended Finite Element Method (XFEM) [31], Embedded Discontinuity Method [32,33,34], and Finite Element Method with universal meshes [35,36,37], explicitly represent fractures as discrete discontinuities. In contrast, smeared fracture methods, including Continuum Damage Mechanics [38], Gradient Damage Mechanics [39], and the phase-field method [40], model fractures as distributed regions integrated into a continuum framework using concepts of non-locality [41].

In the current research, two numerical methods are adopted to explore fracture behavior in porous media under conditions relevant to GREAT cell benchmarks: the Variational Phase-Field (VPF) method [40,42] and lower-dimensional interface elements (LIE) [43]. The VPF method provides a robust framework for simulating fracture nucleation, propagation, and branching under hydro-mechanical coupling. In contrast, the LIE approach explicitly represents fracture surfaces and their hydro-mechanical interactions using lower-dimensional elements. Together, these methods offer a comprehensive understanding of fracture processes in coupled hydro-mechanical systems within porous materials. Both approaches are implemented within the finite element framework of OpenGeoSys [44,45].

The paper is structured as follows. A brief description of the methodology, including the experimental concept of the GREAT cell experiments and the theoretical basis for the variational phase field (VPF-FEM) and lower-order interface element (LIE-FEM) methods, which represent two different approaches to fracture mechanics, is briefly summarised in Section 2.2. The benchmarking procedure is presented in Section 3, starting with the concept and parameter settings (Section 3.1). The results obtained by both the VPF-FEM and LIE-FEM methods for intact and fractured rock samples are presented in detail in the following subsections 3.2, 3.3. Fracture propagation is discussed in Section 3.4. In the discussion, Section 4, the results of the VPF-FEM and LIE-FEM methods are compared with respect to effective permeability, fracture width, and breakdown pressure. The conclusions, Section 5, summarise the main results of the present work and give an outlook on future work. Moreover, the present work contributes to digitisation and open science in THM modelling [46] by providing a set of Jupyter notebooks for documentation and reproducibility of the presented benchmarks (see Appendix and Supplementary Material).

2. Methodology

2.1. Experimental concept

The GREAT cell, depicted in Figure 2, is designed to reproduce subsurface conditions in the laboratory down to a depth of 3.5 km in terms of geomechanical stress and geothermal reservoir conditions on 200 mm diameter rock samples containing fracture networks, thereby enabling these predictions to be validated. The cell represents a major new development in experimental technology, uniquely creating a truly poly-axial rotating stress field, facilitating fluid flow through the samples, and using state-of-the-art fibre-optic strain sensing for high-resolution circumferential strain measurements combined with fluid pressure, capable of thousands of detailed measurements per hour. The GREAT cell will help to better understand and predict fracture propagation and subsequent fluid flow characteristics, which is critical for geo-energy technologies that exploit favourable geological conditions to store or extract fluids from the subsurface. It is well known that fracture permeability decreases non-linearly with increasing normal stress, but the relationship between shear displacement and fracture permeability is less well understood [24,26].

Figure 1. (a) Typical rock sample and (b) simplification for benchmarking of GREAT cell.

Figure 1. (a) Typical rock sample and (b) simplification for benchmarking of GREAT cell.

Figure 2. The GREAT cell facility.

Figure 2. The GREAT cell facility.

In [21], we compared numerical results with experimental data of GREAT cell. However, fracture mechanics under coupled conditions in a novel experimental facility pose challenges for accurate numerical reproduction. This study does not include experimental comparisons; instead, it focuses on developing a benchmark suite tailored for researchers, particularly within the DECOVLEX project. This suite can serve model comparison in a setting that closely resembles GREAT cell tests and can help design tests to be performed specifically for model calibration and validation.

2.2. Mathematical models

The theory of linear poro-elasticity proposed by Biot [47] is widely used to model the deformation of fluid-saturated porous media with a nearly-incompressible fluid (i.e., a fluid with constant density and a finite bulk modulus).Under quasi-static conditions, the governing equations for poro-elastic deformation and fluid flow in porous media and fractures are:

$$\begin{array}{cc}\hfill & {\rho }_{\mathrm{f}}\left({\displaystyle \frac{{\alpha }_{\mathrm{m}}-\varphi }{{\kappa }_{\mathrm{s}}}}+{\displaystyle \frac{\varphi }{{\kappa }_{\mathrm{f}}}}\right)pt+{\rho }_{\mathrm{f}}{\alpha }_{\mathrm{m}}(\nabla ·)t-\nabla ·\left({\rho }_{\mathrm{f}}{\displaystyle \frac{\mathit{K}}{\mu }}\left(\nabla p-{\rho }_{\mathrm{f}}\mathbf{g}\right)\right)\hfill \end{array}$$

$$\begin{array}{cccc}\hfill & -{\delta }_{\Gamma }(\mathbf{x}-{\mathbf{x}}_{\Gamma })q_l={\rho }_{\mathrm{f}}{q}_{\mathrm{m}}^{\mathrm{s}},\hfill & \hfill & \mathrm{in}\Omega \setminus \Gamma ,\hfill \end{array}$$

$$\begin{array}{cccc}\hfill & {\rho }_{\mathrm{f}}wt-{\nabla }_{\Gamma }·\left({\rho }_{\mathrm{f}}{\displaystyle \frac{w^3}{12\mu }}\left({\nabla }_{\Gamma }p-{\rho }_{\mathrm{f}}\mathbf{g}\right)\right)+{\rho }_{\mathrm{f}}{q}_{\mathrm{l}}={\rho }_{\mathrm{f}}{q}_{\mathrm{fr}}^{\mathrm{s}},\hfill & \hfill & \mathrm{in}\Gamma .\hfill \end{array}$$

where $\rho$ is the density of the porous medium, given by $\rho =\varphi {\rho }_{\mathrm{f}}+(1-\varphi ){\rho }_{\mathrm{s}}$, where $\varphi$ represents the porosity. ${\rho }_{\mathrm{f}}$ is the fluid density, and ${\rho }_{\mathrm{s}}$ is the solid density; p is the pore pressure, and ${\alpha }_{\mathrm{m}}$ is the Biot coefficient. ${\kappa }_{\mathrm{s}}$ and ${\kappa }_{\mathrm{f}}$ are the bulk moduli of the solid matrix and a nearly-incompressible fluid, respectively. $\mathit{K}$ is the permeability tensor, and $\mu$ is the dynamic viscosity; w is the fracture aperture. $q_{\mathrm{m}}^{\mathrm{s}}$ and $q_{\mathrm{fr}}^{\mathrm{s}}$ are source terms for the matrix and fracture, respectively, while $q_{\mathrm{l}}$ represents the leak-off term. ${\nabla }_{\Gamma }$ represents the gradient operator along the fracture plane, and ${\delta }_{\Gamma }(\mathbf{x}-{\mathbf{x}}_{\Gamma })$ is a Dirac delta function supported on the fracture surface.

The poro-elastic stress definition, based on Biot theory [47], is given by $\mathit{\sigma }={\mathit{\sigma }}^{\prime }-{\alpha }_{\mathrm{m}}p\mathit{I}$, where $\mathit{\sigma }$ is the total stress, ${\mathit{\sigma }}^{\prime }$ is the effective stress, and $\mathit{I}$ is the identity tensor. The constitutive relation for poro-elastic materials describes the incremental relationship between stress and strain as $\Delta {\mathit{\sigma }}^{\prime }=\mathcal{C}:\Delta \mathit{\epsilon }$, where $\Delta {\mathit{\sigma }}^{\prime }$ and $\Delta \mathit{\epsilon }$ represent incremental changes in effective stress and linearized strain, respectively, and $\mathcal{C}$ is the fourth-order elasticity tensor. The linearized strain tensor is defined as $\mathit{\epsilon }={\nabla }^s={\displaystyle \frac{1}{2}}(\nabla +{(\nabla )}^{\mathrm{T}})$, which expresses strain in terms of the displacement gradient while ensuring symmetry.

Section 2.2 combines mass conservation and Darcy’s law for fluid flow in the matrix (porous medium). Section 2.2 addresses mass conservation in fractures using the lubrication approximation. Boundary conditions for these equations involve specifying displacements or tractions, and pressures or fluxes on the relevant boundaries and fracture interfaces.

2.2.1. Variational Phase Field Approach

The variational phase-field method is widely used in computational fracture mechanics because it effectively models fracturing across different media. Its applications are brittle fracture [48,49,50,51,52,53], ductile fractures [54,55], fatigue failure [56,57], creep-induced fractures [58], dynamic fracture [59,60,61], desiccation cracks [62,63], CO₂/hydraulic fracturing [64,65,66], fractured reactive porous media [67,68], grain size effects [69], nuclear waste repositories [70,71], and fractography [72,73], among others. Its versatility and robustness make it a valuable tool for fracture propagation analysis.

Grounded in Griffith’s fracture criterion [74], the variational phase-field approach effectively addresses its inherent limitations, particularly the dependence on high stress singularities for crack initiation. Francfort and Marigo [42] redefined Griffith’s criterion by minimizing total energy, combining strain and fracture surface energy:

$$\mathcal{F}(\mathit{\epsilon }\left(\right),\Gamma ):={\int }_{\Omega \setminus \Gamma }{\Psi }_0\left(\mathit{\epsilon }\left(\right)\right)\mathrm{d}\Omega +{\int }_{\Gamma }{G}_{\mathrm{c}}\mathrm{d}\Gamma .$$

(1)

This model overcomes the original criterion’s limitations, enabling fracture propagation modeling without known paths and handling nucleation and branching. However, numerically implementing crack surface energy for complex geometries is challenging due to the evolving discrete crack set $\Gamma$.

To overcome this, Bourdin et al. [49] proposed a variational phase-field model with a scalar phase-field variable $v:\Omega \to [0,1]$ and a regularization length parameter $\mathsf{\ell }>0$. The regularized total energy is:

$${\mathcal{F}}_{\mathsf{\ell }}(\mathit{\epsilon },v):={\int }_{\Omega }\Psi (\mathit{\epsilon },v)\mathrm{d}\Omega +{\int }_{\Omega }{\displaystyle \frac{G_{\mathrm{c}}}{4c_{\mathrm{n}}}}\left({\displaystyle \frac{{(1-v)}^n}{\mathsf{\ell }}}+\mathsf{\ell }{|\nabla v|}^2\right)\mathrm{d}\Omega ,$$

(2)

with the ${\mathtt{AT}}_1$ model ($n=1$, $c_n=2/3$) and the ${\mathtt{AT}}_2$ model ($n=2$, $c_{\mathrm{n}}=1/2$) [50]. $\Psi (\mathit{\epsilon },v)$ is the degraded strain energy (see the Appendix H for various strain energy decompositions); For brevity, we denote the strain tensor $\mathit{\epsilon }\left(\right)$ simply as $\mathit{\epsilon }$.

To extend the formulation to hydro-mechanical fracturing, one can decompose the strain energy of a porous medium based on Biot’s work  [75] as follows:

$$\mathcal{F}(\mathit{\epsilon },v,\zeta ):=\underset{\mathrm{dry}\mathrm{strain}\mathrm{energy}}{\underbrace{{\int }_{\Omega }\Psi (\mathit{\epsilon },v)\mathrm{d}\Omega }}+\underset{\mathrm{wet}\mathrm{strain}\mathrm{energy}}{\underbrace{{\int }_{\Omega }{\displaystyle \frac{M}{2}}{\left(\zeta -{\alpha }_m\nabla ·\right)}^2\mathrm{d}\Omega }}+{\int }_{\Omega }{\displaystyle \frac{G_{\mathrm{c}}}{4c_{\mathrm{n}}}}\left({\displaystyle \frac{{(1-v)}^n}{\mathsf{\ell }}}+\mathsf{\ell }{|\nabla v|}^2\right)\mathrm{d}\Omega ,$$

(3)

the first term represents the strain energy of the dry, fractured medium (dry), modeled using phase-field methods. The second term accounts for poro-elastic coupling in a saturated, fractured medium (wet), with $\zeta ={\displaystyle \frac{p}{M}}+{\alpha }_{\mathrm{m}}\nabla ·$. Here, ${\displaystyle \frac{1}{M}}={\displaystyle \frac{{\alpha }_{\mathrm{m}}-\varphi }{{\kappa }_{\mathrm{s}}}}+{\displaystyle \frac{\varphi }{{\kappa }_{\mathrm{f}}}}$ (Biot modulus) measures fluid-saturated material compressibility.

Recent studies [76,77,78] show that ${\alpha }_{\mathrm{m}}$ varies in presence of fracture and in phase-field models, the drained bulk modulus degrades as $g\left(v\right)\kappa$, where $\kappa ={\displaystyle \frac{E}{3(1-2\nu )}}$. You and Yoshioka [65] proposed that Biot coefficients depend on the phase field and splitting strain energy approaches. Therefore, the Biot coefficients are expressed as functions of v and the strain tensor $\mathit{\epsilon }$:

$$\begin{array}{cc}\hfill \alpha (v,\mathit{\epsilon })& =1-\left[g\left(v\right)\mathcal{H}\left(\mathrm{tr}\right(\mathit{\epsilon }\left)\right)+\mathcal{H}\left(\mathrm{tr}\right(-\mathit{\epsilon }\left)\right)\right]{\displaystyle \frac{\kappa }{{\kappa }_{\mathrm{s}}}}.\hfill \end{array}$$

(4)

To model fluid flow in fractured porous media, considering both the matrix and fractures, the phase field variable plays a key role in adjusting the permeability tensor. Accordingly, we introduce an anisotropic permeability tensor $\mathit{K}$ defined as:

$$\mathit{K}\left(v\right)=K_{\mathrm{m}}\mathit{I}+{(1-v)}^{\xi }\left({\displaystyle \frac{w_{\mathrm{n}}^2}{12}}-K_{\mathrm{m}}\right)\left(\mathit{I}-{\mathit{n}}_{\Gamma }\otimes {\mathit{n}}_{\Gamma }\right),$$

(5)

where, $K_{\mathrm{m}}$ represents the intrinsic permeability of the matrix, $\xi \ge 1$ adjusts the impact of fractures on permeability, w denotes the fracture aperture following the cubic law, and ${\mathit{n}}_{\Gamma }$ denotes the normal vector to the fracture plane. By approximating the displacement jump vector $\mathbf{w}$ in terms of the maximum principal strain ${\epsilon }_1$ and a characteristic length scale $h_{\mathrm{e}}$, we derive:

$$w_{\mathrm{n}}\approx h_{\mathrm{e}}{\mathit{\epsilon }}_1,$$

(6)

in this approach, we adopt a concept similar to that proposed by You and Yoshioka [65], while the normal vector is defined based on the principal strain direction. This definition provides greater robustness near the fracture tip and avoids the challenges associated with using $\mathbf{n}=\nabla v/|\nabla v|$, resulting in a simpler and more effective solution for our model.

The ?? and Section 2.2 for hydro-mechanical fracturing using a variational phase field approach are summarized as follows:

$$\begin{array}{cc}\hfill \nabla ·\left({\mathit{\sigma }}^{\prime }(,v)-\alpha (v,\mathit{\epsilon })p\mathit{I}\right)+\rho \mathbf{g}=\mathbf{0}& \mathrm{in}\Omega ,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \left({\displaystyle \frac{\alpha (v,\mathit{\epsilon })-\varphi }{{\kappa }_{\mathrm{s}}}}+{\displaystyle \frac{\varphi }{{\kappa }_{\mathrm{f}}}}\right){\displaystyle \frac{\partial p}{\partial t}}+\alpha (v,\mathit{\epsilon })\nabla ·{\displaystyle \frac{\partial }{\partial t}}-\nabla ·\left({\displaystyle \frac{\mathit{K}\left(v\right)}{\mu }}\left(\nabla p-{\rho }_{\mathrm{f}}\mathbf{g}\right)\right)=q_{\mathrm{s}}& \mathrm{in}\Omega ,\hfill \end{array}$$

(8)

here, $\alpha (v,\mathit{\epsilon })$ represents the phase-dependent Biot coefficient; $\mathit{K}\left(v\right)$ denotes the phase field-dependent permeability tensor, Equation (5); p is the pressure, and the continuity of pressure across fractures ensures $p=p_{\mathrm{fr}}=p_{\mathrm{m}}$, where $p_{\mathrm{fr}}$ and $p_{\mathrm{m}}$ are fracture pressure and matrix pressure, respectively. Additionally, Effective stress ${\mathit{\sigma }}^{\prime }(,v)$ depends not only on displacement but also on the phase field variable v. Since the fluid is incompressible, Equation (8) is divided by ${\rho }_{\mathrm{f}}$ to obtain a formulation based on volumetric fluxes (compare Section 2.2).

2.2.2. Lower-dimensional interface elements

The lower-dimensional interface element approach is specifically designed for pre-existing fractures and fracture paths within the framework of the FEM [43]. This technique represents fractures using lower-order elements and incorporates a partial modification of the Extended Finite Element Method (XFEM) [31].

One of its key features is the local enrichment of the finite element approximation space of displacement, allowing for discontinuities across fractures. This is mathematically expressed as

$${}^h\left(\mathbf{x}\right)=\sum _{i\in I}{N}_i^m{\left(\mathbf{x}\right)}_i^{^}+\sum _{i\in I_{\mathrm{d}}}{N}_i^m\left(\mathbf{x}\right){\psi }_{\mathrm{d}}\left(\mathbf{x}\right){\widehat{\mathbf{w}}}_i,$$

(9)

where ${}^h\left(\mathbf{x}\right)$ denotes the FE displacement approximation at point $\left(\mathbf{x}\right)$, $\widehat{\mathbf{w}}$ is the displacement jump, N stands for the shape function, ${\psi }_{\mathrm{d}}\left(\mathbf{x}\right)$ the global enrichment function, I is the set of nodes in the domain, and $I_{\mathrm{d}}\in I$ is a subset of all nodes in the lower order interfaces.

In addition to displacement enrichment, the method accounts for stress within fractures by employing a localized constitutive law. For elastic deformation, this relationship is defined as

$$\mathrm{d}{\mathit{\sigma }}^{\mathrm{f}}=\mathcal{K}\mathrm{d}\mathbf{w},$$

(10)

where $\mathcal{K}$ is the stiffness tensor whose coefficient matrix in the fracture’s natural coordinate system reads

$$\mathcal{K}=\left[\begin{array}{cc}{k}_{\mathrm{tt}}& k_{\mathrm{tn}}\\ k_{\mathrm{nt}}& k_{\mathrm{nn}}\end{array}\right],$$

(11)

with $k_{\mathrm{tt}}$ and $k_{\mathrm{nn}}$ are the joint shear and normal stiffness, respectively. $k_{\mathrm{tn}}$ and $k_{\mathrm{nt}}$ govern the coupling effects between normal and shear displacements.

The governing equations for hydro-mechanical processes in fractured porous media involve a coupled system that accounts for both momentum balance and mass conservation within the porous matrix $\Omega$ and the the fracture interface $\Gamma$. The momentum balance equation, incorporating the enriched finite element approximation of displacement, is given by:

$$\nabla ·\left({\mathit{\sigma }}^{\prime }\left(\right)-\alpha p\mathit{I}\right)+\rho \mathbf{g}=\mathbf{0},\mathrm{in}\Omega .$$

(12)

Mass conservation in the porous matrix and the fracture follows:

$$\begin{array}{cc}\hfill \left({\displaystyle \frac{{\alpha }_{\mathrm{m}}-\varphi }{{\kappa }_{\mathrm{s}}}}+{\displaystyle \frac{\varphi }{{\kappa }_{\mathrm{f}}}}\right)pt& +{\alpha }_{\mathrm{m}}\nabla ·{\displaystyle \frac{\partial }{\partial t}}-\nabla ·\left({\displaystyle \frac{{\mathit{K}}_{\mathrm{m}}}{\mu }}(\nabla p-{\rho }_{\mathrm{f}}\mathbf{g})\right)\hfill \end{array}$$

$$\begin{array}{ccc}\hfill -{\delta }_{\Gamma }(\mathbf{x}-{\mathbf{x}}_{\Gamma })(q_{+}+q_{-})=q_{\mathrm{m}}^{\mathrm{s}},& \hfill & \hfill \mathrm{in}\Omega ,\end{array}$$

(13)

$$\begin{array}{cccc}\hfill wt+w{S}_{\mathrm{fr}}pt& -{\nabla }_{\Gamma }·\left(w{\displaystyle \frac{{\mathit{K}}_{\mathrm{fr}}}{\mu }}({\nabla }_{\Gamma }p-{\rho }_{\mathrm{f}}\mathbf{g})\right)+q_{+}+q_{-}=0,\hfill & \hfill & \mathrm{on}\Gamma .\hfill \end{array}$$

(14)

Here, $S_{\mathrm{fr}}={\varphi }_{\mathrm{f}}/{\kappa }_{\mathrm{f}}$ is the specific storage, where fluid compressibility, $1/{\kappa }_{\mathrm{f}}$, is considered in the storage term, while the fluid density ${\rho }_{\mathrm{f}}$ is assumed to be constant, following a nearly-incompressible formulation; ${\mathit{K}}_{\mathrm{fr}}$ and ${\mathit{K}}_{\mathrm{m}}$ represent the intrinsic permeabilities of the fracture and matrix, respectively. The terms $q_{+}$ and $q_{-}$ describe fluid exchange between the fracture and the surrounding matrix. Both mass balance equations are normalized by the fluid density ${\rho }_{\mathrm{f}}$.

3. Numerical Benchmarks

3.1. Concept

A comprehensive benchmarking exercise is performed to establish the models’ capabilities and required settings, such as the required discretization for comparable accuracy. These benchmarks are simplified versions of the GREAT cell experiments, focusing on general features like fracture patterns, rotating boundary conditions, and coupled processes. These benchmark exercises are conducted in 2D (plane strain), representing a horizontal section through the middle of the rock samples tested in the GREAT cell (Figure 1b). This 2D plane-strain approach allows the main characteristics of hydro-mechanical (HM) fracture mechanics to be studied efficiently. The benchmarking concept involves a stepwise increase in complexity, including both mechanical (M) and hydro-mechanical (HM) approaches.

We distinguish two main types of benchmarks: intact and fractured rock tests (Figure 3a). In the first type, an intact rock sample is subjected to rotating external stress conditions () and fluid is injected from a central borehole ():

• : Intact rock sample under rotating boundary conditions.
• : Intact rock sample with fluid injection.

For the fractured rock samples, two subtypes are also considered: static and propagating fractures. Static fractures are studied as full (, ) and half-fractured (, ) versions to have a symmetric and non-symmetric case. Again, mechanical and hydro-mechanical versions are studied, the latter with fluid injection into the rock fracture:

• : Fully fractured sample.
• : Half-fractured sample.
• : Fully fractured sample with inflow-outflow.
• : Half-fractured sample with inflow-outflow.

Finally, fracture propagation conditions are studied by injecting fluid into a pre-existing fracture () and fracture nucleation from borehole ():

• : Propagation of pre-existing fracture by fluid injection.
• : Fracture nucleation from borehole by fluid injection.

We repeated each of these benchmarks under different loading conditions and material properties. The material properties are provided in Table 1, and the fluid properties are listed in Table 2. The computational model incorporates two distinct elastic materials within its domain: a central circle ($r=0.097$ m) of rock surrounded by a rubber sheath in a 2D configuration.

The boundary conditions applied in the simulation include both Dirichlet and Neumann conditions. For the momentum balance equation in the sample domain, Dirichlet boundary conditions are used to fix displacement at specific points and prevent rigid body rotation by enforcing zero displacement at selected locations. Mathematically, these conditions are defined as:

$$\left\{\begin{array}{cc}{u}_x(0,y)=0,u_y(0,y)=0\hfill & \mathrm{for}y=-0.09894\mathrm{m},\hfill \\ u_y(x,0)=0,\hfill & \mathrm{for}x=-0.09894\mathrm{m}.\hfill \end{array}\right.$$

where $u_x$ and $u_y$ are the displacement components in the x and y directions, respectively. Neumann boundary conditions define the normal stress ${\sigma }_{\mathrm{rr}}$ on the Pressure-Exerting Elements (PEE’s) and Dynamic-Sealing Strips (DSS’s).

A triaxial loading system is used to generate a sinusoidal response in the sample. To study the effect of poly-axial rotating loading, this load is rotated (e.g., Load A is rotated 45° clockwise to generate Load B). The loading specifications for the PEE’s are detailed in Table 3. The stress on the DSS’s,

$${\sigma }_{\mathrm{DSS}}^i={\displaystyle \frac{{\sigma }_{\mathrm{PEE}}^i+{\sigma }_{\mathrm{PEE}}^{i+1}}{2}},$$

(15)

is determined as the average of the stresses applied to the adjacent PEE’s. The arrangement of PEEs and DSSs is provided in Figure 3b.

To illustrate this loading condition, Figure 4 shows the schematic of loading applied to the GREAT cell for loads A, B, and C (see Table 3). For simplicity in the figure, the DDS load is not included but can be approximated by the average of the stresses applied to the adjacent PEE’s.

3.2. Intact Rock Tests

3.2.1. Intact Rock Sample under Rotating Boundary Conditions–

The primary objective of the first benchmark exercise is to evaluate the mechanical deformation of different rock samples in the GREAT cell under various boundary conditions in a 2D plane strain setup. The benchmark considers a plane strain triaxial stress condition, where the in-plane principal stresses are defined as ${\sigma }_2=10$ MPa and ${\sigma }_3=1$ MPa. Meanwhile, the out-of-plane stress ${\sigma }_1$ is governed by the plane strain constraint, given by, ${\sigma }_1=\nu ({\sigma }_2+{\sigma }_3)$ for both Greywacke and Gneiss samples.

Figure 5 presents the trace of effective stress (${\sum }_{i=x,y,z}{\sigma }_{ii}^{\prime }$) and volumetric strain versus angle for intact samples: (a) Greywacke and (b) Gneiss, using the VPF-FEM method as an example. Note that in Figure 5, the stress distribution is illustrated across the entire domain, while the distribution of volumetric strain is shown within the region $\sqrt{x^2+y^2}<0.065$ m for the benchmarks. The strain in the PEE’s is relatively high due to their soft material properties, which is why the visualization of strain results is limited to the region near the core of the sample.

Figure 6a illustrates the volumetric strain versus angle for intact Greywacke and Gneiss samples subjected to triaxial loading conditions (Load A, see Table 1). The results indicate that VPF-FEM and LIE-FEM effectively replicate similar deformation responses, demonstrating strong agreement between the two methods for intact samples. Greywacke experiences greater negative volumetric strain compared to Gneiss under poly-axial loading because of its lower Young’s modulus and higher Poisson’s ratio, which contribute to increased bulk compressibility.

3.2.2. Intact Rock Sample with Fluid Injection–

This benchmark, which does not consider any fracture, is designed to verify the basic computational setting for hydro-mechanical simulations. In addition to the mechanical loads, a zero constant pore pressure is prescribed at the outer boundary. At the center of sample, fluid is injected at a rate of $Q_0^{\mathrm{v}}=2.085\times {10}^{-9}$ m³/s,

The hydro-mechanical simulations follows a two-stage process: a 3000 s equilibrium phase under mechanical loading to stabilize initial conditions, followed by a 500 s fluid injection phase to model fluid flow (see Section 4.1 for more details).

Figure 7 presents the pressure, trace of effective stress, and volumetric strain profiles of the Greywacke and Gneiss samples for the test under Load A at the final time $t_{\mathrm{f}}=3500\mathrm{s}$ (e.g. LIE-FEM). For the volumetric strain profile the results of $r>0.065$ m are again excluded. This benchmark highlights the hydro-mechanical contrast between Greywacke and Gneiss under fluid injection. Greywacke, with lower stiffness, higher permeability, and greater porosity, enables widespread strain, indicating higher bulk compressibility. In contrast, Gneiss, with higher stiffness, lower permeability, and lower porosity, experiences localized pressure buildup, higher stress concentration, and restricted strain, leading to lower compressibility.

Additionally, Figure 6b illustrates the volumetric strain at the probe-circle for different samples using both VPF-FEM and LIE-FEM. The results show that an increase in pressure at the center of the domain leads to a decrease in effective stress, resulting in a reduction in the magnitude of compressive volumetric strain compared to .

3.3. Fractured Rock Tests

3.3.1. Static Fracture Samples

In this section, numerical simulations are performed on static fracture samples using both VPF-FEM and LIE-FEM, comparing their results.

To simulate the effect of fracture orientation on strain without generating a new mesh for each case, the fracture is kept identical while the poly-axial loading is rotated. The material properties of the sample are listed in Table 1, and the PEE and DDS loads are detailed in Table 3.

Fully Fractured Sample–

In this benchmark, 2D plane strain numerical simulations are performed to establish a baseline model for assessing the impact of fracture orientation on strain distribution under poly-axial loading. A planar fracture within the specimen, $\Gamma =\left[-0.094,0.094\right]\times \left\{0\right\}$, is considered, under poly-axial loading applied on PEE’s and DSS’s (Table 3).

Figure 8 illustrates the profiles of stress and volumetric strain for the Greywacke sample in the test subjected to Load B at the final time $t=t_{\mathrm{f}}$ (LIE-FEM and VPF-FEM). As the Figure 8 indicates, the VPF-FEM results exhibit larger strain near fracture tip.

Figure 10a, Figure 10c, and Figure 10e illustrate the volumetric strain at the probe-circle for Greywacke and Gneiss samples under varying loading conditions. The spikes in the VPF-FEM and LIE-FEM curves indicate a jump in volumetric strain as the probe circle passes through the fracture locations.

The results demonstrate strong consistency away from the fracture. However, at the intersection of the probe-circle and the fracture, the VPE-FEM method shows a higher absolute value of volumetric strain. These differences arise from fundamental distinctions between the approaches. In the LIE-FEM model, both normal and shear fracture stiffness are explicitly accounted for, whereas in the VPF-FEM framework, fracture behavior under mixed loading is dictated by the chosen strain energy decomposition (see Appendix G). Here, the Volumetric–Deviatoric strain energy decomposition is adopted. A sensitivity analysis exploring the impact of different strain energy decompositions in VPF-FEM and fracture stiffness in LIE-FEM is provided in Section 4.2.

Half Fractured Sample–

In this section, all input data and boundary conditions provided in , except a single wing fracture, defined as $\Gamma =\left[0,0.04\right]\times \left\{0\right\}$, is considered.

Figure 9 illustrates the profiles for trace of stress and volumetric strain obtained from both the VPF-FEM and LIE-FEM methods under the test subjected to Load B at the final time $t=t_{\mathrm{f}}$. Both methods capture the stress singularity at the fracture tips.

Figure 10b, Figure 10d, and Figure 10f compare the volumetric strain versus angles at the probe-circle for Greywacke and Gneiss samples under various loading conditions. The results demonstrate strong overall agreement, with differences between the two methods diminishing as the probe-circle moves farther from the fracture tip.

Figure 9. : Trace of effective stress and volumetric strain profiles of the Greywacke sample for the test under Load B at the final time $t=t_{\mathrm{f}}$ (LIE-FEM & VPF-FEM).

Figure 9. : Trace of effective stress and volumetric strain profiles of the Greywacke sample for the test under Load B at the final time $t=t_{\mathrm{f}}$ (LIE-FEM & VPF-FEM).

Figure 10. /: Volumetric strain versus angle for fractured samples under different loading conditions as detailed in Table 3.

Figure 10. /: Volumetric strain versus angle for fractured samples under different loading conditions as detailed in Table 3.

Fully Fractured Sample with Fluid Injection–

In order to explore how fracture permeability varies across different materials and different load conditions, a series of hydro-mechanical simulations () is performed. Five distinct load conditions (Table 3) are considered, each designed to create different angles between the second principal stress direction (${\sigma }_2=10$ MPa) and the plane of PEE 5. A planar fracture within the specimen $\Gamma =\left[-0.094,0.094\right]\times \left\{0\right\}$, is analyzed.

For the fluid flow, the initial pressure throughout the domain is set as $p_0=0.1$ MPa. Dirichlet boundary conditions are imposed such that $p(-0.094,0)=p(0.094,0)=3.45$ MPa. After allowing the system to stabilize for 3000 s to reach equilibrium, fluid is injected at the center with an injection rate of $Q_0^{\mathrm{v}}(0,0)=4.167·{10}^{-7}{\mathrm{m}}^3/\mathrm{s}$ (25 ml/min) for 500 s. Furthermore, the initial fracture width is considered to be $w_{\mathrm{ini}}=1\times {10}^{-6}$ m.

Figure 11 presents the pressure, effective stress trace, and volumetric strain profiles of the Greywacke sample from the test under Load B at the final time step ($t=t_{\mathrm{f}}$), using LIE-FEM and VPF-FEM. As illustrated in Figure 11b, the stress within the recaptured region diminishes to nearly zero. This phenomenon occurs due to the Isotropic model used in VPF-FEM for this benchmark. When the phase field variable v reaches zero, it leads to complete stress attenuation at the fracture.

Figure 13a, Figure 13c, and Figure 13e illustrate volumetric strain versus angle for fractured samples under different loading conditions, as detailed in Table 3 for . A jump is observed where the probe-circle intersects the fracture. This jump is more significant in the VPF-FEM method for hydro-mechanical fractures due to the higher strain values within the fracture captured by VPF-FEM.

Figure 14a, Figure 14c, and Figure 14e show the fracture width along the fracture for fractured samples under different loading conditions. The fracture width in all cases is of a similar order; however, differences arise due to variations in the numerical models. In VPF-FEM, the width is calculated using a strain-based method, whereas in the LIE-FEM method, it is explicitly determined from displacement.

Figure 12. : Pressure, trace of effective stress, and volumetric strain profiles of the Greywacke sample for the test under Load B at the final time $t=t_{\mathrm{f}}$ (LIE-FEM & VPF-FEM).

Figure 12. : Pressure, trace of effective stress, and volumetric strain profiles of the Greywacke sample for the test under Load B at the final time $t=t_{\mathrm{f}}$ (LIE-FEM & VPF-FEM).

Half Fractured Sample with Fluid Injection–

In this section, all input data and boundary conditions provided in are used. A single wing fracture, defined as $\Gamma =\left[0,0.04\right]\times \left\{0\right\}$, is considered. For the boundary conditions related to fluid flow, a Dirichlet boundary condition with a pressure of $p(0.04,0)=3.45$ MPa is applied on the right side of the fracture. After allowing the system to stabilize for 3000 s to reach equilibrium, fluid is injected with an injection rate of $Q_0^{\mathrm{v}}(0,0)=4.167·{10}^{-7}{\mathrm{m}}^3/\mathrm{s}$ (25 ml/min) on the left side of the fracture.

Figure 12 presents the pressure, stress trace, and volumetric strain profiles of the Graywacke sample from the test under Load B at the final time step ($t=t_{\mathrm{f}}$), using LIE-FEM and VPF-FEM.

Figure 13b, Figure 13d, and Figure 13f illustrate volumetric strain versus angle for fractured samples under different loading conditions, as detailed in Table 3 for . As shown, there is no jump in this set of figures since the probe-circle does not pass through the fracture.

Figure 14b, Figure 14d, and Figure 14f illustrate the fracture width along the fracture for fractured samples under different loading conditions. Figure 14f (Load C, VPF-FEM) shows that the fracture width remains nearly unchanged. The higher stress perpendicular to the fracture plane under Load C prevents opening, likely due to the Isotropic model used in VPF-FEM for this benchmark. In contrast, LIE-FEM shows clear fracture opening for Load C.

Figure 13. /: Volumetric strain versus angle for fractured samples under different loading conditions as detailed in Table 3.

Figure 13. /: Volumetric strain versus angle for fractured samples under different loading conditions as detailed in Table 3.

Figure 14. /: Fracture width along the fracture for fractured samples under different loading conditions, as detailed in Table 3.

Figure 14. /: Fracture width along the fracture for fractured samples under different loading conditions, as detailed in Table 3.

3.4. Fracture propagation

3.4.1. Hydraulic Fracturing (pre-existing fracture)–

In this set of benchmarks, hydraulic fracturing is simulated under various boundary conditions. Fluid is injected at the center of an inclined fracture to model the fracturing process. To evaluate the impact of poly-axial stress boundary conditions on fracture propagation, three load scenarios (A, B, and C in Table 3) are defined. A planar, initially inclined fracture, $\Gamma =\left\{(x,y)\mid x=scos{45}^{\circ },y=ssin{45}^{\circ },s\in \left[-{\displaystyle \frac{0.02}{2}},{\displaystyle \frac{0.02}{2}}\right]\right\}$ is considered. The initial fracture width is set to $1\times {10}^{-7}$ m with an initial pressure of 0.1 MPa. The radial stress is applied at the boundary ( Table 3), and it is allowed to equilibrate over 3000 s. Subsequently, an injection rate of $Q_0=1\times {10}^{-5}\times h$ m²/s is applied at a single line element ($h=1\times {10}^{-3}$ m) at the center of the fracture. The rest of the material and fluid properties used in the simulations are presented in Table 1 and Table 2. Here, the Volumetric-Deviatoric strain energy decomposition is used.

Figure 15 presents the phase field and pressure profiles under different loading conditions for the Greywacke sample at $t=500$ s after the start of injection. The results demonstrate the influence of poly-axial loading on the direction of fracture propagation. In all cases, it is evident that the direction of hydraulic fracturing tends to align parallel to the maximum stress at the boundaries.

In Scenario Load A, where the minimum principal stress is applied in the y-direction, the fracture primarily propagates perpendicular to this direction (Figure 15a). However, a slight counterclockwise deviation from the horizontal direction observed due to the effects of poly-axial loading. Although the maximum principal stress acts horizontally, the overall stress distribution (see Figure 4a) influences the fracture path, causing it to deviate slightly from a purely horizontal trajectory. In Scenario Load B, where the minimum principal stress is aligned with the fracture plane, the fracture deviates into direction of maximum principal stress (Figure 15b). In Scenario Load C, with the minimum principal stress applied in the x-direction, the hydraulic fracture predominantly propagates vertically (Figure 15c). Therefore, the computational results align well with expected hydraulic fracture propagation mechanisms, demonstrating that the numerical model effectively captures key processes governing hydraulic fracture evolution under different boundary conditions.

3.4.2. Hydraulic fracturing (borehole)–

In the final benchmark, a benchmark set is introduced where fractures nucleate from a borehole under various poly-axial conditions. The injection rate is set to $Q_0=2\times {10}^{-7}\times 2\pi r$ m²/s at the borehole, where $r=0.005$ m. The material and fluid properties used in the simulations are detailed in Table 1 and Table 2. Note that no notch is considered in this case, and injection occurs directly inside the borehole.

Figure 16 displays the phase field and pressure profiles for the Greywacke sample under different loading conditions at $t=500$ s after injection began. The fracture nucleation and propagation tend to align parallel to the maximum overall stress at the boundaries.

4. Discussion

The VPF-FEM and LIE-FEM approaches are two different numerical methods of dealing with fracture problems. VPF-FEM is a smeared method; it represents fractures using a diffuse model represented by a phase-field variable. It is powerful because VPF-FEM can capture crack propagation without re-meshing, can handle crack branching and arbitrary crack paths, and is therefore suitable for modeling both fracture nucleation and propagation. However, to accurately capture the smeared fracture zone, a fine mesh resolution near the fracture is required. In contrast, the pre-existing fracture in LIE-FEM is explicitly embedded in the mesh and is therefore more discrete.

The primary distinction between these two approaches lies in their treatment of fracture mechanics. VPF-FEM does not explicitly define a fracture stiffness parameter; instead, it represents fractures using a degraded continuum approach, where stiffness gradually decreases based on a diffuse phase (damage) field. In contrast, LIE-FEM explicitly models fractures by incorporating stiffness parameters ($k_{\mathrm{nn}},k_{\mathrm{tt}}$), which dictate mechanical behavior under normal and shear loads. Regarding fracture propagation, both methods integrate fracture toughness but in different ways. In VPF-FEM, crack evolution is described through a phase-field variable, which captures the gradual transition from an intact material state to full fracture. This evolution is governed by an energy functional that incorporates fracture toughness, ensuring that crack propagation adheres to an energy-minimization principle, whereas LIE-FEM can employs cohesive zone models, where fracture toughness controls discrete fracture growth (though this aspect is not considered in this study). A notable advantage of VPF-FEM is its ability to model fracture nucleation, as demonstrated in benchmark , which examines hydraulic fracture initiation from a borehole. For shear behavior, VPF-FEM manages shear stiffness through energy decomposition, limiting independent control over shear strength. LIE-FEM, on the other hand, explicitly defines shear stiffness via $k_{\mathrm{tt}}$, enabling independent control over shear resistance, even in closed fractures.

4.1. , – Intact rock samples

In the benchmark cases without fractures, and , the results obtained with VPF-FEM and LIE-FEM are identical. This consistency is attributed to the identical application of rotating boundary conditions in both methods and equivalence of their underlying poro-elastic model formulations.

Effect of initial conditions on hydro-mechanical model behavior

Accurate hydro-mechanical simulations require careful attention to sequential boundary conditions, especially in coupled fluid flow modeling for fractured media. In this study, a two-step approach is used to ensure numerical stability and physical accuracy. The first step includes a 3000 s steady-state phase under mechanical loading to allow the system to reach equilibrium (Figure 17b). This phase minimizes unrealistic pore pressure variations (e.g. Figure 17a) and establishes stable initial conditions for subsequent simulations. In the second step, a 500 s fluid injection phase is performed to model the flow dynamics. This sequential approach ensures that the effects of mechanical deformation are fully accounted for before fluid flow is introduced, a consideration that is especially critical in low-permeability media where strong coupling between fluid flow and deformation can lead to numerical instabilities. For VPF-FEM models, careful initialization is critical to achieving realistic fracture nucleation. Without proper initialization, artificial fracture nucleation can occur during the steady-state phase, leading to nonphysical results.

4.2. , – Fracture mechanics

In this set of benchmarks, the influence of an existing fracture on mechanical deformation is evaluated. For smaller fractures (), the volumetric strain measured along the probe circle shows strong agreement between VPF-FEM and LIE-FEM methods. However, for larger fractures (), where the fracture intersects the probe circle, both methods exhibit jumps in volumetric strain. While their results remain consistent at locations further from the fracture. To address this, three studies are conducted: the effect of mesh size in VPF-FEM, strain energy decomposition in VPF-FEM, and the impact of fracture stiffness in LIE-FEM.

Effect of mesh size in VPF-FEM

A mesh refinement analysis is performed for VPF-FEM to evaluate its behavior relative to LIE-FEM in modeling both volumetric strain and stress. Here, the Volumetric-Deviatoric strain energy decomposition is used in VPF-FEM. Figure 18a shows the volumetric strain at the probe circle for different mesh sizes, compared to LIE-FEM results. Figure 18b and Figure 18c illustrate the von Mises stress and strain along a vertical line perpendicular to the fracture at $x=0$, highlighting the role of shear stress and strain in shaping these results. The findings reveal that mesh size does not significantly affect the magnitude of the volumetric strain jump. However, it can influence the maximum von Mises strain within the fracture ($y=0$). These differences arise from fundamental variations in how each method models fractures. LIE-FEM explicitly represents fractures, whereas VPF-FEM uses a smeared phase-field approach, leading to smoother transitions. The LIE-FEM results are particularly sensitive to the fracture stiffness parameters ($k_{\mathrm{nn}}$ for normal deformation and $k_{\mathrm{tt}}$ for shear deformation). In contrast, VPF-FEM results are strongly influenced by the strain energy decomposition choice.

Effect of strain energy decomposition approaches in VPF-FEM

To analyze the effect of strain energy decomposition on fracture behavior, the benchmark is repeated using four strain energy decomposition approaches for VPF-FEM (see Appendix H). Figure 19a presents the volumetric strain results, showing that the Isotropic and Volumetric-Deviatoric models exhibit volumetric strain jumps, similar with the LIE-FEM method, due to their ability to redistribute strain. In contrast, the Spectral and No-tension models display minimal volumetric strain jumps, as they confine fracture (propagation) strictly to tensile regions, preventing volumetric expansion near the crack zone. Figure 19b further illustrates the Von Mises stress distribution along the perpendicular line at (0,0) in the y-direction. It highlights that both Isotropic and Volumetric-Deviatoric approaches lead to nearly zero Von Mises stress and high Von Mises strain inside the fracture at $y=0$.

Effect of fracture stiffness in LIE-FEM

is evaluated using the LIE-FEM method across different fracture stiffness values. The Figure 21 illustrates that as fracture stiffness increases, the magnitude of jumps in volumetric straining within the LIE-FEM method decreases.

In LIE-FEM, all fractures incorporate a tension cut-off mechanism, ensuring that no traction is transferred across a fracture once the normal traction becomes non-negative. This enforces the condition ${\mathbf{t}}_{\Gamma }^{\prime }·{\mathbf{n}}_{\Gamma }\le 0$.

While the previous discussion highlights the performance differences between VPF-FEM and LIE-FEM, it is important to address an issue observed in the LIE-FEM method when using a quadrilateral mesh (see Figure A24a). This configuration introduces volumetric locking, which affects the accuracy of the strain and stress results. To address this, the B-bar method introduced by Hughes (1980) [79] has been used. However, this problem is not observed in the VPF-FEM method due to its formulation. More details on the B-bar method can be found in the Appendix G.

4.3. , – Effective permeability evolution under different poly-axial loading

This set of benchmarks evaluates the evolution of effective permeability under poly-axial loading, taking into account hydro-mechanical responses in the presence of fractures. Benchmarks and represent the most comprehensive scenarios, capturing the interactions between strain, fracture mechanics, and fluid flow. The volumetric strain plots for these cases indicate that both VPF-FEM and LIE-FEM can reproduce the hydro mechanical behavior of samples under different loading conditions.

The key concept explored here is the evolution of effective permeability under varying poly-axial loading conditions. In this setup, the fracture is oriented perpendicular to the main experimental element (PEE 5/5a, Figure 3b). As the angle $\alpha$ between the second principal stress ${\sigma }_2=10$ MPa and the PEE 5/5a plane increases, the normal stress acting on the fracture plane decreases monotonically. The normal stress ${\sigma }_{\mathrm{n}}$ for each load case is determined as follows:

$${\sigma }_{\mathrm{n}}={\displaystyle \frac{1}{2}}\left[\left({\sigma }_2+{\sigma }_3\right)+({\sigma }_2-{\sigma }_3)cos\left(2\alpha \right)\right],$$

where ${\sigma }_2=10\mathrm{MPa}$, ${\sigma }_3=1\mathrm{MPa}$, and $\alpha$ takes values of $0^{\circ },22.5^{\circ },{45}^{\circ },67.5^{\circ },\mathrm{and}{90}^{\circ }$ (Figure 3b). Effective permeability at the fracture center is calculated for each loading condition (Table 3) using the cubic law for ideal parallel plate fractures [80], $k={\displaystyle \frac{b^2}{12}}$, where b represents the effective hydraulic aperture.

The results shown in Figure 22 illustrate the evolution of permeability for Graywacke and Gneiss samples under full and half wing fractures ( and ). As expected, fracture permeability decreases with increasing normal stress in the fracture plane, in agreement with theoretical predictions. In , the VPF-FEM method predicts larger permeabilities than the LIE-FEM method. However, in , the LIE-FEM method predicts higher permeabilities. Despite these differences, both methods show good agreement—not only are they of the same order, but the differences remain small. This slight variation is attributed to the nature of width calculation in each method.

4.4. , – Hydraulic Fracturing

This section focuses exclusively on hydraulic fracturing results obtained with the VPF-FEM method, as fracture propagation is not modeled with the LIE-FEM method. Two benchmarks are considered in this context: , which examines fracture propagation from a pre-existing fracture, and , which simulates fracture nucleation from a central borehole (Figure 3a). The breakdown (peak or failure) pressure, which marks the onset of fracture propagation, is analyzed for Gneiss and Graywacke samples under different loading conditions (Load A, B, C, Figure 4). The results for these scenarios are presented in Figure 23a and Figure 23b. In the absence of pre-existing fractures, as in , the breakdown pressure reflects a more realistic scenario where it is primarily determined by the material properties of the specimens rather than the poly-axial loading conditions. Conversely, in the presence of pre-existing fractures (, Figure 23a), the breakdown pressure varies slightly depending on the orientation of the applied load relative to pre-existing fracture. Specifically, when the maximum applied stress aligns with the fracture plane, the breakdown pressure is lower (Load A in Figure 23a) due to reduced resistance to propagation. In contrast, when the loading direction is nearly perpendicular to the fractures, the breakdown pressure increases because the stress concentration at the fracture tips is lower, making fracture initiation more difficult and delaying propagation.

5. Conclusions

Based on the GREAT cell experiments [24,25,26,81], a systematic collection of benchmarks has been developed to test various aspects of fracture mechanics for crystalline rocks (Gneiss and Graywacke). The special feature of the GREAT cell compared to conventional triaxial tests is that polyaxial stress ratios (rotating stress conditions) can be applied to the specimens, so that in particular fracture properties can be investigated under variable stress conditions as they may occur in geotechnical applications such as deep geothermal energy, energy storage, and disposal.

The benchmark types include intact rocks (, ) as well as symmetrically (, ) and asymmetrically fractured samples (, ). The benchmark types were simulated using two different fracture mechanics approaches, the variational phase field (VPF) and the lower interface element (LIE) methods, both implemented in the finite element method (FEM). The differences between diffuse and explicit representations of fracture in the context of mechanical and hydro-mechanical processes have been worked out (see below). In addition, fracture propagation processes have been analyzed using VPF-FEM method (, ).

Key findings from the comprehensive benchmark study include

• For the mechanical benchmarks and , sinusoidal strain patterns are successfully reproduced by both models. Gneiss exhibits smaller strain amplitudes due to its higher stiffness, while Graywacke exhibits larger total compressive strains.
• In static fracture modeling (, ), stress singularities are captured by both methods, showing good agreement far from the fracture. However, the magnitude of volumetric strain jumps, when the probe circle passes through the fracture, is affected by the choice of fracture stiffness in LIE-FEM and strain energy decomposition in VPF-FEM. Consistent global deformation responses are produced by both approaches.
• For the hydro-mechanical cases and , stress redistribution is affected by permeability variations, with normal stress at the fracture plane decreasing as the principal stress angle $\alpha$ increases.
• For propagating fractures (and ), the direction of fracture propagation is accurately predicted by VPF-FEM, and the breakdown pressure is consistent with fracture mechanics expectations.

While LIE-FEM provides a more explicit representation of discontinuities, VPF-FEM excels in modeling fracture evolution. However, LIE-FEM is more computationally efficient because it increases the degrees of freedom only in the fracture regions, whereas VPF-FEM requires an additional phase field variable per node. In addition, LIE-FEM explicitly incorporates fracture stiffness, unlike VPF-FEM, which affects stress calculations at discontinuities. Despite previous analytical verification (see [22,30,43,65,68]), real-world inconsistencies persist due to differences in fracture width representation and fracture stiffness modeling, highlighting the need for further investigation.

This collection of benchmarks serves two purposes. First, it is a reference for model comparison with other fracture mechanics methods, such as discrete element methods, cellular automata, or, in the future, machine learning. On the other hand, it serves as a starting point for the evaluation of more complex GREAT cell experiments, e.g. considering 3D effects, fracture networks and anisotropic material behavior. The current benchmark collection initially covers mechanical and hydro-mechanical processes and will be extended to thermal processes in the next version as part of the DECOVALEX 2027 task SAFENET-2. This will allow to cover the whole range of thermo-hydro-mechanical (THM) fracture mechanics. The benchmark collection will also be successively published as a collection of Jupyter notebooks, thus facilitating the integration of further numerical methods and codes and contributing to open science.