WP 3: Pathways through stress transitions


Overleaf logo 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.


The random surface generator provided by Candela et al. 20091 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.





1Thibault 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).

Fig1 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).

Fig2 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).

Fig3 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.

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.

Fig 2 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.


Figure 3: Analysis of three connected fractures embedded in a three dimensional porous medium. Comparison of pressure evolution with a fixed and deformation dependent fracture permeability.


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.

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

LIE Figure 4: Verifi cation of the LIE method.


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

shear stress - disp curve 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.

Scherkurve & shear normal  stress 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.

fig1 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).

fig2 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 eff ect while leak-off phenomena are currently not taken into account.

figure 1 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 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 Figure 3: Left: Terzaghi's 1-dim consolidation benchmark problem. Right: Pressure di ffusion.

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 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 Figure 5: Left: Spatial discretization of matrix material and fracture. Center: Velocity component
in ex (horizontal) direction at t → 0. Right: Evolving velocity component in ex 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 Figure 6: Fracture with asperities geometry.


[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 Z2s, Rs and the parameters C, θ*max, A0 (according to Grasselli et al. 2002).

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

By comparison of the lower and the upper surface of a rock joint the mechanical aperture is calculated and zones where the fracture is open or closed can be identifi ed. Interesting is the case when one surface is translated like it will be done in the shear experiments. Though this first tests are merely geometrical the typical small contact areas after some displacement are observed. Another observation is that  flow paths (i.e. apertures greater zero) seem to form preferably perpendicular to the shear direction. This behavior is in accordance to results reported in the literature (e.g. Ivars et al. 2017, Nguyen et al. 2014) and should be taken in mind when interpreting shear tests because usually the water flow is parallel to the shear direction in the experimental setups.

Figure 2 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 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).


[1] Ali Fathi et al. \Geometric E ect 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:
[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 E ect 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.