Geomechanical Integrity of Host and Barrier Rocks -
Experiment, Modeling and Analysis of Discontinuities

COMSOL created mesh with DFN-mesh, matrix-mesh, and tunnel geometry [TUBAF] In addition to the open-source DFN generators (dfnWorks and frackit), a commercial software “Comsol Multiphysics” was utilized for complex mesh generation. The DFN approach integrated with Comsol was used to generate fractures, connect them to the matrix and account for underground structures. The generated mesh files can be easily exported for the further use in OpenGeoSys framework. With the help of "FEconv" program [1], the generated networks can be converted from the COMSOL format (.mphtxt) into a compatible OpenGeoSys format (.vtu). The custom workflow can be employed to generate a variety of fractures having different orientations. This approach complements the already existed open-source workflows, making it feasible to create more complex geometries.

Right scale: velocity in the fractures, Left scale: velocity in the matrix (OGS simulation)
In the following picture the flow to the tunnel inside a fractured rock is simulated (slice in the centre of the geometry). The inflow area are the fractures on the right side of the geometry. The major flow occurs via the highly permeable fractures while drainage into the tunnel occurs through the matrix.

In linking GeomInt2 to other projects, members of our department (Geotechnical Institute Technische Universität Bergakademie Freiberg) bridge the gap between mostly quasistatic THM models and methods for wave propagation. These results may lead to acceleration spectra, supporting the further design of the DGR and its components according to KTA or Eurocode [2].
Mesh of modelled domain with two exemplary fault planes (a), heterogeneous permeability distribution (b), maximum change in Coulomb failure stress (ΔCFS < 0 indicates failure) during a glacial cycle on first fault plane (c) and on second fault plane (d) [2]

[1] Iban Constenla, Victor Sande, Francisco Pena (2012). FEconv – Finite Element mesh Conversor. http://victorsndvg.github.io/FEconv/index.html

[2] Kern, D., Deng, T., Magri, F., Malkovsky, V. I., Nagel, T., (2022). Advancing transient simulation of hydro-mechanically coupled systems in geological disposal applications. DAEF 2022 (Deutsche Arbeitsgemeinschaft Endlagerforschung), 4th to 6th July 2022, Köln.

[TUBAF] Within the framework of another project at TUBAF, several different DFN software packages were tested with respect to their ability to create complex fracture networks and to generate input compatible with OpenGeoSys. There, the DFNs will be used to simulate flow and transport in fracture networks. In GeomInt-2, these methods form the basis for extension to HM-coupled problems using the methods developed for single fractures in GeomInt.

Left: velocity field of a 3D geometry with DFN-network Right: pressure field of a DFN-network for only fracture flow

Two articles were published recently using methods developed in GeomInt and GeomInt2. Li et al. [1] use OpenGeoSys to model large-scale cavern structures in rock salt by making use of high-performance computing capabilities of the code and respective clusters at German and Chinese research institutions. Kafle et al. [2] simulated evolving localization patterns in fine-grained soils due to saturation changes in order to assess slope stability in light of precipitation as well as water level fluctuations in a hydro-reservoir. A fully coupled hydro-mechanical model for partial saturation was linked to the OGS#MFront [3] interface and an algorithm for strength reduction.

[1] Li, J., Zhang, N., Xu, W., Naumov, D., Fischer, T., Chen, Y., Zhuang, D., & Nagel, T. (2022). The influence of cavern length on deformation and barrier integrity around horizontal energy storage salt caverns. Energy, 123148. https://doi.org/10.1016/j.energy.2022.123148
[2] Kafle, L., Xu, W., Zeng, S., & Nagel, T. (2022). A numerical investigation of slope stability influenced by the combined effects of reservoir water level fluctuations and precipitation: A case study of the Bianjiazhai landslide in China. Engineering Geology, 297(May 2021), 106508. https://doi.org/10.1016/j.enggeo.2021.106508
[3] Helfer, T., Bleyer, J., Frondelius, T., Yashchuk, I., Nagel, T., & Naumov, D. (2020). The MFrontGenericInterfaceSupport project. Journal of Open Source Software, 5(48), 2003. https://doi.org/10.21105/joss.02003

[UoS] Hydro-mechanical coupling of fractured porous media is a challenging task. Properties of the fracture domain vary distinctly from the properties characterizing the intact poro-elastic media and lead to an overall numerically “stiff” system of PDEs. To guarantee numerical stability special treatment of the resulting system of equations is required. Due to the distinct numerical properties in the fracture and intact region, monolithic schemes result in ill-conditioned systems that require direct solution strategies. Direct numerical solvers lack scalability in terms of parallel computing and are therefore not ideal for computationally large problems. Hence, we proposed a new partitioned coupling scheme that incorporates advanced quasi-newton stabilization for the coupling of both numerical domains that can then be solved by individual solvers and standard preconditioning using highly parallelized standard solvers. Our work has been summarized in a publication that has recently been accepted [1].

Flow field obtained for the fracture domain using the proposed partitioned coupling scheme to hydro-mechanically couple the poro-elastic and fracture domain applying a quasi-Newton based approach

[1] Schmidt, P., Jaust, A., Steeb, H. and Schulte, M., Simulation of flow in deformable fractures using a quasi-Newton based partitioned coupling approach, Computational Geosciences, (2022)

[UoS] Two articles were recently published within the GeomInt project. In Schmidt et al. (2021) [1], the earlier proposed hybrid-dimensional flow model was extended by an in-situ stress state relation to capture normal contact mechanisms and was applied to pressure and flow transients from experiments conducted by the collaborative project STIMTEC at the Reiche Zeche underground laboratory. The study could relate the transient change in fracture permeability and non-linear aperture changes to the length of a fracture and the contribution of its geometrical and normal stiffness. In Schmidt et al. (2021) [2], we had the chance to apply the hybrid-dimensional model to a unique set of pressure, flow and fracture deformation data obtained from pumping operations at Grimsel test site. Based on numerical fits of the measurement data, we could show that fracture deformation contributes to the compensation of injected fluid volume in large amounts, whereas the contribution of the fluid’s compressibility and therefore pressure diffusion is negligible. In this context we discussed the limits of pressure diffusion models for the characterization of fractured reservoirs and highlighted the importance of hydro-mechanical modelling even under small to moderate pressure amplitudes.

[1] Schmidt, P., Steeb, H. & Renner, J. Investigations into the opening of fractures during hydraulic testing using a hybrid-dimensional flow formulation. Environ Earth Sci 80, 497 (2021), https://doi.org/10.1007/s12665-021-09767-4
[2] P Schmidt, N Dutler, H Steeb, Importance of fracture deformation throughout hydraulic testing under in-situ conditions, Geophysical Journal International, 2021, ggab354,

[TUBAF] One of the main problems in multiple fracture propagation simulations is the computation time. To overcome the problem a verified and efficient method like LIE method (lower-dimensional interface elements) is implemented in OpenGeoSys. The LIE method is verified in a two-dimensional setup with the Sneddon analytical solution, and also a three-dimensional setup is validated with the Mont Terri underground laboratory experiment. The 3D-single fracture propagation model requires several days simulation run time. The simulation of a multiple fracture setup could take around several weeks or even months due to the strong coupling between hydraulics and mechanics and the large number of elements. Thus, it is necessary to find an efficient solution to reduce the computation time and have a more usable setup. We build an efficient workflow to create a multiple fracture mesh for numerical simulations.

[UoS] Conventional hydraulic field tests on single fractures and/or fractured reservoirs provide data sets consisting of pressure and flow transients. Experimentally measured data of the transient fracture opening are barely available and limited to information obtained in the injection borehole. Recent experiments conducted at the Grimsel Test Site provide a consistent set of pressure, flow and fracture opening data. The transient fracture aperture changes were measured in a distance of seven meters from the injection borehole. Based on the unique experimental data set we could perform an inverse numerical analysis of the hydro-mechanical properties of the tested fracture by applying the proposed hybrid-dimensional model for flow processes in deformable fractures. The numerical data successfully reproduce the whole set of measured pressure, flow and fracture opening transients. Further, we could show the importance of fracture deformation throughout hydraulic testing and the limitations/failure of purely pressure diffusion-based models. The comprehensive study on hydro-mechanical fracture flow processes under in-situ conditions led to a recently submitted manuscript.

Comparison of experimental pressure, flux, fracture deformation and total injected fluid volume protocols with numerical results

The random surface generator provided by Candela et al. 2009 was used to analyze the characteristics of results for different kinds of surfaces.

There are two possible curve shapes: curves with a clear peak and curves which reach a plateau.

Shear tests have been conducted using basalt. The sample was multiply sheared using different boundary
conditions. A photography of the rock sample was taken after the last shear test. Comparing it with
the simulation results shows a good congruency, see Fig. 1. One major difference between experiment
and simulation is that the left upper corner was entirely destroyed in the experiment whereas in the
simulation just a damaged zone was detected.

Figure 1: Comparison of the basalt surface after the last shear test with the simulation result of the
destructed zones in red.

 Characterization of fractures by inverse analysis using numerical calculations gain complexity once injected flow rates lead to mechanical induced changes of the fracture volume. Throughout step rate tests at the Reichen Zeche test site conducted by collaborating researchers of the STIMTEC project, non-linear pressure-flow relationships could be found in the measurement protocols indicating changes in the fracture’s permeability induced by the injected fluid volume. In order to characterize the tested fractures and estimate the stiffness and initial aperture, calculations of a fully coupled hybrid-dimensional model has been applied. The numerical simulations could successfully reproduce the measurement data and non-linear pressure-flow relationships be explained with changes of the fracture volume and permeability. Throughout this process the characteristic fracture-stiffness is a crucial aspect since it controls the volumetric response to the applied transient pressure perturbations. Visualization of the fracture-stiffness emphasize its non-constant behaviour in space and in time related to the deformation state of a fracture. The cross-link of numerical simulations and experimental field data provided a unique possibility to study hydro-mechanical effects throughout testing procedures on the field scale.

Figure 1 Comparison of experimental pressure, flow and injected fluid volume protocols with numerical results

Figure 2 Stiffness evolution in time over fracture length

TUBAF

The random surface generator provided by Candela et al. 2009¹ was used to analyze the characteristics of results for different kinds of surfaces. Two parameters describe the roughness of the surfaces: the Hurst exponent and the standard deviation of the height data. The Hurst exponent, H, which is a dimensionless measure of the fractal characteristics of a surface, varies between 0 and 1. A small value of H corresponds to a rough surface whereas a higher value characterizes a smooth surface. The standard deviation of the height data, σ(h), indicates the overall height distribution of the generated surface. A series of calculations has been carried out, and averaged results are displayed as surfaces, as can be seen in the figures.

There are two possible curve shapes: curves with a clear peak and curves which reach a plateau. It is interesting to see which parameters result in which curve type. The coloured edges indicate identical edges for comparison.

¹Thibault Candela et al. Characterization of Fault Roughness at Various Scales: Implications of Three-Dimensional High Resolution Topography Measurements. In: Pure and Applied Geophysics 166:10 (Oct. 2009), pp. 1817-1851. ISSN: 1420-9136. DOI: 10.1007/s00024-009-0521-2. URL: https://doi.org/10.1007/s00024-009-0521-2.

The research group at TU Bergakademie Freiberg is currently working on shear tests in the rock mechanical lab. Lab tests formerly conducted within the GeomInt-project were done on coarse-grained granite (Saxony, Germany). In order to supplement the knowledge on HM-coupled behavior on crystalline rock with (very) fine grains, basalt rock from Thuringia (Germany) was recently chosen as sample material (Fig. 1).

Fig. 1: Close-up of fine-grained basalt rock (long side of image is 20 cm)
We adopt a multi-stage testing technique and systematically vary several test boundary conditions. In a primary quasi-static stage we apply constant normal load on the shear plane while stimulating with different shear velocities. The investigated range of conditions for stress and speed is 1 to 2.3 [MPa] and 1 to 100 [mm/min], respectively. Subsequently, shear plane behaviour under dynamic loading conditions is examined. Thus, a sinusoidal dynamic component is superimposed to the static normal stress on shear plane. Again, we sweep through frequencies and amplitudes to examine a variety of dynamic boundary conditions. Mechanical characteristics (behavior) and deformation (wear) of the shear plane is monitored during and after each test stage. By visual inspection wear of the shear surface, asperity degradation and macroscopic damage can be investigated and categorized (Fig. 2).

Fig. 2: Visual inspection of shear plane before testing (left) and after dynamic loading (right)
Quantification of damage pattern is done with the aid of a 3D-surface scanner. The cloud of points in 3D-space is triangulated using different meshing tools and statistically analyzed afterwards (Fig. 3).

Fig. 3: Triangulated point clouds gained by 3D-scanning of the shear plane before testing (left) and after dynamic loading (right)

After completion of these experiments the next step will be to conduct shear tests on rock samples using constant normal stiffness boundary conditions. Normal Stiffness of Shear-Plane will be varied in range of approx. 1.0 to 5.0 [MPa/mm]. We aim to document the dilatancy behaviour of a shear plane in crystalline rock under varying confining conditions which may be seen as restraint due to the surrounding rock mass. Thus we can draw conclusions regarding the size and width of possible flow paths.

Fracture flow processes in deformable fractures requires a tight hydro-mechanical coupling. Non-local fluid pressure induced fracture volume changes lead to inverse water level fluctuations throughout hydraulic well stimulations. Besides volumetric phenomena local fracture permeability changes due to its deformation state may influence the pressure diffusion. Consistent coupling of fluid flow and fracture deformation is numerically challenging and results in a highly non-linear flow formulation of the governing equations shown in Fig. 1.

Figure 1: Governing equation for fracture flow in deformable fractures.

From a numerical point of view different schemes have been evaluated in terms of stability and efficiency. An iterative scheme allowing for non-conformal meshes for fluid and solid domains is tested against a strongly coupled scheme that solves the problem monolithically. Due to the number of iterations needed to reach convergence like shown in Fig. 2 the strongly coupled scheme is preferred for complex problems especially in three dimensions.

Figure 2: Comparison of iterative vs strongly coupled scheme.

The influence of permeability changes throughout fracture stimulation on the pressure diffusion evolution is finally studied on a three-dimensional fracture geometry showing that a deformation dependent permeability has a high impact on the pressure evolution. The experimental set-up is shown in Fig. 3 along with graphs of the pressure distribution and permeability for different time steps.

UFZ

The following three methods have been implemented in OGS for fracture propagation in crystalline formation (Fig. 3):
(a) Lower-dimensional Interface Element (LIE)
(b) Variational Phase-Field (PF)
(c) Non-Local Damage (NLD)
They have been verified against the closed form solution, which is shown for the LIE case in Fig. 4. The study outcome is summarised in a manuscript submitted to International Journal of Geomathematics, which is currently under review.

Figure 3: In OGS implemented methods to simulate fracture propagation.

Figure 4: Verification of the LIE method.

TUBAF

Numerical simulations were done using the surface model from surface scans and basic rock parameters as only inputs. Ideas from Casagrande et al. 2017¹ were used and further developed. Fig. 5 and Fig. 6 show first results.

Figure 5: Peak and residual shear stress against shear displacement. The peak stress represents a typical curve for a first shear cycle and the residual stress a curve for a repeated shear cycle.

Figure 6: Left: Lab data of a direct shear test. Right: Numerical simulation. It is obvious there is need for further developments to capture the behaviour of the lab experiments in the simulation.

[1] Davide Casagrande et al. "A New Stochastic Approach to Predict Peak and Residual Shear Strength of Natural Rock Discontinuities". In: Rock Mechanics and Rock Engineering (Aug. 2017), pp. 69-99.

Fluid flow through fracture networks may cause phenomena such as precipitation and crystallization, modifying the effective properties of the porous media. This may be advantageous in applications like CO2 sequestration, but in the case of clogging injection and extraction wells, also may cause problems in
applications like geothermal energy generation.
Fig. 1 shows the geometry of a micro-fluidic cell obtained with XRCT. On the left side, the initial state of the cell is shown. The right side shows the cell after calcite has formed.

Figure 1: Fractures with asperities. Left: Initial state micro-fluidic cell. Right: Micro-fluidic cell with calcite precipitates.
The calcite introduces asperities to the previously smooth fracture walls (Fig. 1 - bottom). To study the influence that the asperities play in the fluid flow through the porous media, a 2-dim domain is generated applying a random-walk approach (Fig. 2).

Figure 2: Fluid flow in fluid-filled fracture with asperities.
Fig. 2 (left) shows the simulation of a periodic flow through an artificially generated fracture with asperities. Fig. 2 (right) shows a close-up of the solid-fluid interface where the effect that asperities cause in the fluid flow is more clear.

An accurate numerical simulation of fluid-filled fractures in porous media (like reservoir rocks) requires to model the coupling between the viscous fluid in the fractures and the poroelastic skele-ton [4]. In the proposed numerical approach, the porous medium is described with the quasi-static poroelastic equations taking into account rigid, i.e. undeformable solid grains and a negligible effective pore fluid density. The present model concentrates on inherent hydro-mechanical coupling effect while leak-off phenomena are currently not taken into account.

Figure 1: Hydro-mechanical coupling of fluid-filled fractures embedded in a poroelastic matrix.

Present investigations of (deformable) fractures includes large (solid) deformations and failure modelling. The use of classical mesh-based methods (like e.g. FEM) implies computationally expensive geometrical meshing and remeshing approaches. Furthermore, employing mesh-based methods may entail other error sources, such as mesh distortions and the alignment of elements along the crack orientation. To overcome these difficulties, we use the (Lagrangian) meshless Smoothed Particle Hydrodynamics (SPH) method.

In the past, SPH has been widely applied for various types of problems in Computational Fluid Dynamics (CFD) especially if free surfaces and/or interfaces play an important role [1, 3]. We have implemented the SPH method into the high performance software HOOMD-Blue. This software for particle simulations is parallelized for GPU and CPU computations. The implementation performance is tested with scalability tests for both parallelizations schemes. Fig. 2 presents the result of strong scalability analysis carried out for single phase fluid flow through a sphere array with high porosity.

Figure 2: Strong scalability test. Left: CPU parallelization. Right: GPU parallelization.

In our present investigation, the fluid-filled fractures are modelled with the Navier-Stokes equation and solved with a quasi-incompressible single phase (SPH) scheme. For the surrounding fluid-saturated porous matrix material, a (hybrid) quasi-static poroelastic approach is proposed. The surface-coupled model has been validated with Terzaghi's 1-dim consolidation benchmark problem [2] (Fig. 3) as well as with more sophisticated 2-dim problems.

Figure 3: Left: Terzaghi's 1-dim consolidation benchmark problem. Right: Pressure diffusion.

First, an analysis is carried out for finding the minimum number of SPH particles for the discretization of the fracture. For creeping flow conditions, the chosen SPH discretization has to recover the parabolic velocity profile of fluid-flow in fractures and thus the physically-sound momentum interaction.

Figure 4: Particles discretization resolution for simulation of Poiseuille flow.

Second, we present (Fig. 5) initial results for the velocity field v(x, t) inside a fracture with low aspect ratio (width/length → 0:5).

Figure 5: Left: Spatial discretization of matrix material and fracture. Center: Velocity component
in e_x (horizontal) direction at t → 0. Right: Evolving velocity component in e_x direction at t » 0.

In the next step, we carry out fractures and fracture networks including complex fracture morphology (e.g. asperities) of single cracks and/or fracture growth. Furthermore, fractures with asperities will also be studied (Fig. 6).

Figure 6: Fracture with asperities geometry.

References

[1] J J Monaghan. Smoothed particle hydrodynamics. Reports on progress in physics, 68(8):1703, 2005.
[2] M Osorno and H Steeb. Smoothed particle hydrodynamics modelling of poroelastic media. PAMM, 16(1):469{470, 2016.
[3] R. Sivanesapillai, N. Falkner, A. Hartmeier, and H. Steeb. Pore-scale SPH study on the evolution of interfacial areas in porous media during primary drainage and imbibition. Advances in Water Resources, 95:212{234, 2016.
[4] C Vinci, J Renner, and H Steeb. A hybrid-dimensional approach for an ecient numerical modeling of the hydro-mechanics of fractures. Water Resources Research, 50(2):1616{1635, 2014.

The triangulated surface scans of a rock joint are analyzed with Matlab to characterize the roughness. The surface consists of  ≈ 300 000 triangles. Roughness parameters include Z_2s, R_s and the parameters C, θ*_max, A₀ (according to Grasselli et al. 2002).

Figure 1: Curve t to get the parameters C, θ*_max and A₀. The shear direction is the positive x-direction.

Figure 2: Aperture map at initial positions. The two surfaces are not matching perfectly resulting in non-zero apertures. At some points the surfaces are simply not identical. At others it is the positioning which has to be improved.

Figure 3: Displacements of 2mm in x-direction (left) and y-direction (right) with additional dilation of 0:3mm have a signicant influence on the spatial aperture distribution. Colors indicating aperture state (see Fathi et al. 2016): green=open, red=overlap and blue=just in contact (here an absolute value of the aperture smaller than threshold of 0:05mm is considered as contact).

References

[1] Ali Fathi et al. \Geometric Eect of Asperities on Shear Mechanism of Rock Joints". In: Rock Mechanics and Rock Engineering 49.3 (Mar. 2016), pp. 801{820. issn: 1434-453X. doi: 10.1007/s00603-015-0799-6. url:
https://doi.org/10.1007/s00603-015-0799-6.
[2] G. Grasselli, J. Wirth, and P. Egger. \Quantitative three-dimensional description of a rough surface and parameter evolution with shearing". In: International Journal of Rock Mechanics and Mining Sciences 39.6 (2002), pp. 789{800. issn: 1365-1609. doi: https://doi.org/10.1016/S1365- 1609(02)00070- 9. url: http://www.sciencedirect.com/science/article/pii/S1365160902000709.
[3] Diego Mas Ivars and Katarina Noren. \Geometric Eect of Asperities on Shear Mechanism of Rock Joints". In: Svensk Karnbranslehantering AB SKB R-13-29 (May 2017). issn: 1402-3091. url: http://www.skb.com/
publications/ (visited on 12/18/2017).
[4] Van-Manh Nguyen, Heinz Konietzky, and Thomas Fruhwirt. \New Methodology to Characterize Shear Behavior of Joints by Combination of Direct Shear Box Testing and Numerical Simulations". In: Geotechnical and Geological Engineering 32.4 (Aug. 2014), pp. 829{846. issn: 1573-1529. doi: 10.1007/s10706-014-9761-8. url: https://doi.org/10.1007/s10706-014-9761-8.