# Assessing Stress Variability in Fractured Rock Masses with Frictional Properties and Power Law Fracture Size Distributions

Etienne Lavoine, Philippe Davy, Caroline Darcel, Diego Mas Ivars, Hossein A Kasani

► To cite this version:

Etienne Lavoine, Philippe Davy, Caroline Darcel, Diego Mas Ivars, Hossein A Kasani. Assessing Stress Variability in Fractured Rock Masses with Frictional Properties and Power Law Fracture Size Distributions. Rock Mechanics and Rock Engineering, 2024, 57, pp.2407-2420. 10.1007/s00603-023-03683-8. insu-04426484

HAL Id: insu-04426484

https://insu.hal.science/insu-04426484

Submitted on 30 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Assessing Stress Variability in Fractured Rock Masses with Frictional Properties and Power Law Fracture Size Distributions

Etienne Lavoine<fnref n="1" /> Ⓒ · Philippe Davy<fnref n="2" /> · Caroline Darcel<fnref n="1" /> · Diego Mas Ivars<fnref n="3" />,<fnref n="4" /> · Hossein A. Kasani<fnref n="5" />

Received: 30 June 2023 / Accepted: 14 November 2023

© The Author(s) 2023

## Abstract

The presence of fractures in rock masses plays a major role in its stress state and its variability. Each fracture potentially induces a stress perturbation, which is correlated to its geometrical and mechanical properties. This work aims to understand and quantitatively predict the relationship between fractured systems and the associated stress fluctuations distribution, considering any regional stress conditions. The approach considers the rock mass as an elastic rock matrix into which a population of discrete fractures is embedded—known as a Discrete Fracture Network (DFN) modeling approach. We develop relevant indicators and analytical solutions to quantify stress perturbations at the fracture network scale, supported by 3D numerical simulations, using various fracture size distributions. We show that stress fluctuations increase with fracture density and decrease as a function of the so-called stiffness length, a characteristic length that can be defined as the ratio between Young's modulus of the matrix and fracture stiffness. Based on these considerations we discuss, depending on DFN parameters, which range of fractures should be modeled explicitly to account for major stress perturbations in fractured rock masses.

## Highlights

* Stress fluctuations in fractured rocks are predicted quantitatively from a tensorial approach.

* At the fracture scale, stress fluctuations depend on fracture size, orientation with respect to the applied remote stress field, and mechanical properties.

* At the network scale, stress fluctuations depend on fracture density, as well as size, orientation and mechanical property distributions.

**Keywords** Stress fluctuations · Fractured rock mass · Fractures · DFN

## 1 Introduction

Evaluating the in-situ stress state is a key point of geomechanical site modeling for many industries such as nuclear waste disposal (Figueiredo et al. 2023; Martin 2007). The stress state at any location results from a combination of regional stress conditions (McGarr and Gay 1978; Zoback 1992), and stress perturbations potentially induced by any significant heterogeneity such as topography (McTigue and Mei 1981; Savage and Swolfs 1986), material heterogeneity (Lei and Gao 2019), or geological structures (Martin and Chandler 1993).

Extensive field data (Barton and Zoback 1994; Yale 2003) have shown the role played by faults and joints in the perturbation of regional stress field. Indeed, the normal and shear displacements on fracture planes cause a deformation

<figure><img src="image_1.png" /></figure>

<footnote>Etienne Lavoine
e.lavoine@itasca.fr</footnote>

<footnote n="1">¹ Itasca Consultants S.A.S, Fractory, University of Rennes, Rennes, France</footnote>

<footnote n="2">² University of Rennes, CNRS, Géosciences Rennes, UMR 6118, 35000 Rennes, France</footnote>

<footnote n="3">³ Swedish Nuclear Fuel and Waste Management Company (SKB), Solna, Sweden</footnote>

<footnote n="4">⁴ Division of Soil and Rock Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden</footnote>

<footnote n="5">⁵ Nuclear Waste Management Organization (NWMO), Toronto, ON, Canada</footnote>
# Assessing Stress Variability in Fractured Rock Masses with Frictional Properties and Power Law Fracture Size Distributions

Etienne Lavoine<fnref n="1" /> Ⓒ · Philippe Davy<fnref n="2" /> · Caroline Darcel<fnref n="1" /> · Diego Mas Ivars<fnref n="3" />,<fnref n="4" /> · Hossein A. Kasani<fnref n="5" />

Received: 30 June 2023 / Accepted: 14 November 2023

© The Author(s) 2023

## Abstract

The presence of fractures in rock masses plays a major role in its stress state and its variability. Each fracture potentially induces a stress perturbation, which is correlated to its geometrical and mechanical properties. This work aims to understand and quantitatively predict the relationship between fractured systems and the associated stress fluctuations distribution, considering any regional stress conditions. The approach considers the rock mass as an elastic rock matrix into which a population of discrete fractures is embedded—known as a Discrete Fracture Network (DFN) modeling approach. We develop relevant indicators and analytical solutions to quantify stress perturbations at the fracture network scale, supported by 3D numerical simulations, using various fracture size distributions. We show that stress fluctuations increase with fracture density and decrease as a function of the so-called stiffness length, a characteristic length that can be defined as the ratio between Young's modulus of the matrix and fracture stiffness. Based on these considerations we discuss, depending on DFN parameters, which range of fractures should be modeled explicitly to account for major stress perturbations in fractured rock masses.

## Highlights

* Stress fluctuations in fractured rocks are predicted quantitatively from a tensorial approach.

* At the fracture scale, stress fluctuations depend on fracture size, orientation with respect to the applied remote stress field, and mechanical properties.

* At the network scale, stress fluctuations depend on fracture density, as well as size, orientation and mechanical property distributions.

**Keywords** Stress fluctuations · Fractured rock mass · Fractures · DFN

## 1 Introduction

Evaluating the in-situ stress state is a key point of geomechanical site modeling for many industries such as nuclear waste disposal (Figueiredo et al. 2023; Martin 2007). The stress state at any location results from a combination of regional stress conditions (McGarr and Gay 1978; Zoback 1992), and stress perturbations potentially induced by any significant heterogeneity such as topography (McTigue and Mei 1981; Savage and Swolfs 1986), material heterogeneity (Lei and Gao 2019), or geological structures (Martin and Chandler 1993).

Extensive field data (Barton and Zoback 1994; Yale 2003) have shown the role played by faults and joints in the perturbation of regional stress field. Indeed, the normal and shear displacements on fracture planes cause a deformation

<figure><img src="image_2.png" /></figure>

<footnote>Etienne Lavoine
e.lavoine@itasca.fr</footnote>

<footnote n="1">¹ Itasca Consultants S.A.S, Fractory, University of Rennes, Rennes, France</footnote>

<footnote n="2">² University of Rennes, CNRS, Géosciences Rennes, UMR 6118, 35000 Rennes, France</footnote>

<footnote n="3">³ Swedish Nuclear Fuel and Waste Management Company (SKB), Solna, Sweden</footnote>

<footnote n="4">⁴ Division of Soil and Rock Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden</footnote>

<footnote n="5">⁵ Nuclear Waste Management Organization (NWMO), Toronto, ON, Canada</footnote>
of the surrounding matrix resulting in a stress concentration at the fracture tips and a stress shadow above and below its central part (Jaeger et al. 2009). All fractures produce stress fluctuations, but their impact is related to their geometrical and mechanical properties, as well as their stress conditions (Homberg et al. 1997). At the network scale, stress interactions make the spatial distribution of the stress field even more complex, as any fracture may affect the loading conditions of the surrounding fractures (Kachanov 1989; Thomas et al. 2017). Understanding the controlling factors for stress fluctuations and interactions represents a key step for the prediction of fracture development (Healy et al. 2006; Kachanov 2003; Olson 1993). The inherent uncertainty in stress measurements and the scarcity of available data compared to the large rock volume of interest motivates the use of numerical simulation to understand the role played by fractures at the rock mass scale (Hakami et al. 2022).

In this paper, we quantitatively analyze and identify controlling factors for stress fluctuations at the network scale, supported by three-dimensional (3D) numerical simulations. We first perform a simple analysis from a single fracture, before addressing the network scale, using various fracture size distributions. We show that the intensity of the stress perturbation at the network scale can be predicted, knowing fracture geometrical and mechanical parameters, and applied stress conditions. Based on these considerations, we discuss the range of fractures that should be represented explicitly when performing geomechanical simulations, to correctly account for stress fluctuations while keeping manageable simulation time.

# 2 Methodology

## 2.1 Discrete Fracture Network (DFN) Model

For geological environments, fracture networks are often characterized by a wide distribution of fracture sizes (Bonnet et al. 2001). We denote the fracture density distribution $n(l, \theta)$, as the statistical characterisation of a fractured system defined from fractures orientation $\theta$ and size $l$ (Selroos et al. 2022), so that $n(l, \theta)dld\theta$ is the number of fractures of size and orientation in the range $[l, l + dl]$ and $[\theta, \theta + d\theta]$ per unit volume of rock. Natural fracture systems, especially in crystalline rocks, can often be described by a power law size distribution model (Bour 2002; Davy 1993; Davy et al. 1990):

$$n(l, \theta) = \alpha(\theta)l^{-a}, \quad (1)$$

with $\alpha(\theta)$ a density term that only depends on fracture orientation and $a$ the power law exponent of the fracture size distribution. This exponent is usually between 3 and 4 in

crystalline rocks (Bonnet et al. 2001; Darcel et al. 2006; Davy et al. 2010). Fracture density can be measured in different ways depending on the dimension of the measurement region and the dimension associated to fractures (Dershowitz and Herda 1992). For fracture networks made of circular fractures of minimum and maximum diameters $l_{min}$ and $l_{max}$, the total fracture intensity $P_{32}$ is defined as the total fracture surface per unit volume:

$$P_{32} = \frac{\pi}{4} \iint_{l,\theta} n(l, \theta)l^2 dld\theta. \quad (2)$$

Also, fracture connectivity may be assessed by the percolation parameter $p$ (Balberg et al. 1984; Bour and Davy 1998), expressed as the total excluded volume around fractures per unit volume:

$$p = \frac{\pi^2}{8} \iint_{l,\theta} n(l, \theta)l^3 dld\theta. \quad (3)$$

The DFN is statistically connected if $p$ is larger than the percolation threshold $p_c$. Considering a 3D DFN made of disk-shaped fractures, the percolation threshold lies within the range 0.7–2.8 (Balberg 1985; Bour and Davy 1998; De Dreuzy et al. 2000).

In this study, we define a set of DFN models following a uniform orientation distribution (all fracture orientations are equally represented) and a power law size distributions of exponent $a = [3, 4]$ with $l_{min} = 0.4$ m and $l_{max} = 4$ m, embedded in a cubic system of size $L = 5$ m, from low to high percolation parameter. These values are chosen to cover a wide range of fracture sizes while keeping the simulation time manageable (the larger the range of fracture sizes, the more fractures and blocks, and the longer the simulation time and memory requirements), and to avoid fractures that completely cut through the domain. Figure 1a, b shows the fracture size distribution of the generated DFNs. We also define “equivalent” DFN models with constant fracture size $l_{a=3}$ and $l_{a=4}$ so that corresponding DFNs following power law size distribution of exponent $a = 3$ and $a = 4$ have the same $P_{32}$ and $p(l_{a=3} = 1.56$ m and $l_{a=4} = 1.02$ m). For each set of parameters, only one realization is generated. In total, 44 DFN models are generated.

In the following, we refer to simulations with the following conventions. A simulation based on a DFN with power law size distribution is noted with the letter “a” and constant size models with a “l”. The density of the network is noted with “p” and the stiffness length with “ls” (see Sect. 2.2). For example, the simulation based on DFN with power law exponent $a = 3$, percolation parameter $p = 4$ and stiffness length $l_s = \infty$ is referred by “a3p4ls∞”. Figure 1 shows examples of generated DFN models.
of the surrounding matrix resulting in a stress concentration at the fracture tips and a stress shadow above and below its central part (Jaeger et al. 2009). All fractures produce stress fluctuations, but their impact is related to their geometrical and mechanical properties, as well as their stress conditions (Homberg et al. 1997). At the network scale, stress interactions make the spatial distribution of the stress field even more complex, as any fracture may affect the loading conditions of the surrounding fractures (Kachanov 1989; Thomas et al. 2017). Understanding the controlling factors for stress fluctuations and interactions represents a key step for the prediction of fracture development (Healy et al. 2006; Kachanov 2003; Olson 1993). The inherent uncertainty in stress measurements and the scarcity of available data compared to the large rock volume of interest motivates the use of numerical simulation to understand the role played by fractures at the rock mass scale (Hakami et al. 2022).

In this paper, we quantitatively analyze and identify controlling factors for stress fluctuations at the network scale, supported by three-dimensional (3D) numerical simulations. We first perform a simple analysis from a single fracture, before addressing the network scale, using various fracture size distributions. We show that the intensity of the stress perturbation at the network scale can be predicted, knowing fracture geometrical and mechanical parameters, and applied stress conditions. Based on these considerations, we discuss the range of fractures that should be represented explicitly when performing geomechanical simulations, to correctly account for stress fluctuations while keeping manageable simulation time.

# 2 Methodology

## 2.1 Discrete Fracture Network (DFN) Model

For geological environments, fracture networks are often characterized by a wide distribution of fracture sizes (Bonnet et al. 2001). We denote the fracture density distribution $n(l, \theta)$, as the statistical characterisation of a fractured system defined from fractures orientation $\theta$ and size $l$ (Selroos et al. 2022), so that $n(l, \theta)dld\theta$ is the number of fractures of size and orientation in the range $[l, l + dl]$ and $[\theta, \theta + d\theta]$ per unit volume of rock. Natural fracture systems, especially in crystalline rocks, can often be described by a power law size distribution model (Bour 2002; Davy 1993; Davy et al. 1990):

$$n(l, \theta) = \alpha(\theta)l^{-a}, \quad (1)$$

with $\alpha(\theta)$ a density term that only depends on fracture orientation and $a$ the power law exponent of the fracture size distribution. This exponent is usually between 3 and 4 in

crystalline rocks (Bonnet et al. 2001; Darcel et al. 2006; Davy et al. 2010). Fracture density can be measured in different ways depending on the dimension of the measurement region and the dimension associated to fractures (Dershowitz and Herda 1992). For fracture networks made of circular fractures of minimum and maximum diameters $l_{min}$ and $l_{max}$, the total fracture intensity $P_{32}$ is defined as the total fracture surface per unit volume:

$$P_{32} = \frac{\pi}{4} \iint_{l,\theta} n(l, \theta)l^2 dld\theta. \quad (2)$$

Also, fracture connectivity may be assessed by the percolation parameter $p$ (Balberg et al. 1984; Bour and Davy 1998), expressed as the total excluded volume around fractures per unit volume:

$$p = \frac{\pi^2}{8} \iint_{l,\theta} n(l, \theta)l^3 dld\theta. \quad (3)$$

The DFN is statistically connected if $p$ is larger than the percolation threshold $p_c$. Considering a 3D DFN made of disk-shaped fractures, the percolation threshold lies within the range 0.7–2.8 (Balberg 1985; Bour and Davy 1998; De Dreuzy et al. 2000).

In this study, we define a set of DFN models following a uniform orientation distribution (all fracture orientations are equally represented) and a power law size distributions of exponent $a = [3, 4]$ with $l_{min} = 0.4$ m and $l_{max} = 4$ m, embedded in a cubic system of size $L = 5$ m, from low to high percolation parameter. These values are chosen to cover a wide range of fracture sizes while keeping the simulation time manageable (the larger the range of fracture sizes, the more fractures and blocks, and the longer the simulation time and memory requirements), and to avoid fractures that completely cut through the domain. Figure 1a, b shows the fracture size distribution of the generated DFNs. We also define “equivalent” DFN models with constant fracture size $l_{a=3}$ and $l_{a=4}$ so that corresponding DFNs following power law size distribution of exponent $a = 3$ and $a = 4$ have the same $P_{32}$ and $p(l_{a=3} = 1.56$ m and $l_{a=4} = 1.02$ m). For each set of parameters, only one realization is generated. In total, 44 DFN models are generated.

In the following, we refer to simulations with the following conventions. A simulation based on a DFN with power law size distribution is noted with the letter “a” and constant size models with a “l”. The density of the network is noted with “p” and the stiffness length with “ls” (see Sect. 2.2). For example, the simulation based on DFN with power law exponent $a = 3$, percolation parameter $p = 4$ and stiffness length $l_s = \infty$ is referred by “a3p4ls∞”. Figure 1 shows examples of generated DFN models.
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Fracture size distributions and examples of generated DFN realizations. **a**, **b** Power-law fracture size distributions with exponent $a = 3$ and $a = 4$ (dotted lines represent the DFN realizations and</figcaption></figure>

dashed lines represent the density model) in a log-log plot. **c–f** DFN realizations of the a3p4, a4p4, I3p4, I4p4 models

## 2.2 The Synthetic Rock Mass Approach

We use the Synthetic Rock Mass (SRM) approach (Mas Ivars et al. 2011) for simulating the mechanical behaviour of the fractured rock mass. In the SRM approach, the rock is modeled as an assembly of deformable blocks delimited by fractures using a DFN representation. The mechanical behaviour of the rock mass is governed by the

interaction of the deformable rock blocks and fractures. Each block contact is divided into sub-contacts, where interaction forces between the blocks are applied. Simulations are performed using 3DEC (Itasca Consulting Group 2020), a three-dimensional numerical software dedicated to discontinuum modeling, based on the distinct element method (DEM). The cubic system of size *L* containing the
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Fracture size distributions and examples of generated DFN realizations. **a**, **b** Power-law fracture size distributions with exponent $a = 3$ and $a = 4$ (dotted lines represent the DFN realizations and</figcaption></figure>

dashed lines represent the density model) in a log-log plot. **c–f** DFN realizations of the a3p4, a4p4, I3p4, I4p4 models

## 2.2 The Synthetic Rock Mass Approach

We use the Synthetic Rock Mass (SRM) approach (Mas Ivars et al. 2011) for simulating the mechanical behaviour of the fractured rock mass. In the SRM approach, the rock is modeled as an assembly of deformable blocks delimited by fractures using a DFN representation. The mechanical behaviour of the rock mass is governed by the

interaction of the deformable rock blocks and fractures. Each block contact is divided into sub-contacts, where interaction forces between the blocks are applied. Simulations are performed using 3DEC (Itasca Consulting Group 2020), a three-dimensional numerical software dedicated to discontinuum modeling, based on the distinct element method (DEM). The cubic system of size *L* containing the
DFN is progressively cut into smaller and smaller blocks when fractures are added to the sample (Fig. 2).

The block assembly is then meshed block by block with a target mesh size *h* that should be at least two times smaller than the smallest fracture contained in the system to correctly describe displacements on the fracture plane. The matrix assumed to be elastic, with Young's modulus *E* and Poisson's ratio *v*. Fracture frictional properties are assigned to the area fraction corresponding to the fractures at the interfaces between the blocks. In the following, we consider simple elastic behaviour of the fracture plane, where the shear resistance stress *τ*<sub>f</sub> in the fracture plane is expressed as a linear function of shear displacement *t* on the fracture walls, modeled by a constant shear stiffness *k*<sub>s</sub>. If the applied stress is compressive, the fracture walls cannot penetrate as normal stiffness *k*<sub>n</sub> is generally much larger than *k*<sub>s</sub> and normal displacement is negligible in comparison to shear displacement. In the following, the surrounding matrix deformation is entirely due to fracture shearing, as we only apply compressive stress, considering *k*<sub>n</sub> ≫ *k*<sub>s</sub>. To focus on the importance of fracture length, we define the fracture stiffness length *l*<sub>s</sub> = *E*/*k*<sub>s</sub> (Davy et al. 2018). Decreasing the stiffness length *l*<sub>s</sub> is equivalent to increasing fracture plane resistance (increasing *k*<sub>s</sub>). The SRM specimens are loaded with conditions such that the maximum principal stress *σ*<sub>a</sub> is vertical and the confinement stress *σ*<sub>l</sub> is isotropic in the lateral directions (*x* and *y* directions). Numerical parameters used to perform the numerical simulations are summarized in Table 1. The elastic properties of the rock matrix are selected from the granite properties of the Forsmark site in Sweden (Hakami et al. 2022; SKB 2008). DFN density and shear stiffness are selected over a wide range of values to perform sensitivity analysis.

<table><caption>**Table 1** Parameters of the SRM specimens</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Domain size</td><td>L</td><td>5 m</td></tr><tr><td>Mesh size</td><td>h</td><td>0.2 m</td></tr><tr><td>DFN percolation parameter</td><td>p</td><td>0.5–8</td></tr><tr><td>Young's modulus</td><td>E</td><td>76.9 GPa</td></tr><tr><td>Poisson's ratio</td><td>ν</td><td>0.23</td></tr><tr><td>Fracture stiffness length</td><td>l<sub>s</sub> = E/k<sub>s</sub></td><td>∞–0.2 m</td></tr><tr><td>Maximum principal stress</td><td>σ<sub>a</sub></td><td>3 kPa</td></tr><tr><td>Confinement stress</td><td>σ<sub>l</sub></td><td>1 kPa</td></tr></tbody></table>

## 2.3 Stress Fluctuations Quantification

The stress field in the fractured rock mass in response to applied tectonic stress is obtained by solving the linear elastic geomechanical problem. At any point, i.e., mesh element, the stress field is defined by a tensor and corresponding invariants. For example, Von Mises stress σ<sub>e</sub> quantifies the intensity of the deviatoric component of a stress tensor. It is defined from stress tensor σ̄ as:

$$ \sigma_e(\bar{\sigma}) = \frac{1}{\sqrt{2}} \sqrt{(\sigma_1 - \sigma_2)^2 + (\sigma_1 - \sigma_3)^2 + (\sigma_2 - \sigma_3)^2}, \quad (4) $$

with σ₁, σ₂ and σ₃ the principal stress components of the stress tensor σ̄. At any point *x* in the 3D cartesian space, σ<sub>e</sub>(*x*) refers to the Von Mises stress at position *x*.

It is possible to analyze distributions of local values, as well as deviation from the mean stress tensor (σ̄) in the volume *V*:

$$ \langle \bar{\sigma} \rangle = \frac{1}{V} \int_V \bar{\sigma}(x) dV, \quad (5) $$

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Numerical setup used to perform simulation from DFN generation (left), block domain discretization (center), and stress field resolution (right). The initial cubic sample is cut by the DFN into smaller</figcaption></figure>

blocks, meshed with a size *h* that must be much smaller than the average spacing between fractures
DFN is progressively cut into smaller and smaller blocks when fractures are added to the sample (Fig. 2).

The block assembly is then meshed block by block with a target mesh size *h* that should be at least two times smaller than the smallest fracture contained in the system to correctly describe displacements on the fracture plane. The matrix assumed to be elastic, with Young's modulus *E* and Poisson's ratio *v*. Fracture frictional properties are assigned to the area fraction corresponding to the fractures at the interfaces between the blocks. In the following, we consider simple elastic behaviour of the fracture plane, where the shear resistance stress *τ*<sub>f</sub> in the fracture plane is expressed as a linear function of shear displacement *t* on the fracture walls, modeled by a constant shear stiffness *k*<sub>s</sub>. If the applied stress is compressive, the fracture walls cannot penetrate as normal stiffness *k*<sub>n</sub> is generally much larger than *k*<sub>s</sub> and normal displacement is negligible in comparison to shear displacement. In the following, the surrounding matrix deformation is entirely due to fracture shearing, as we only apply compressive stress, considering *k*<sub>n</sub> ≫ *k*<sub>s</sub>. To focus on the importance of fracture length, we define the fracture stiffness length *l*<sub>s</sub> = *E*/*k*<sub>s</sub> (Davy et al. 2018). Decreasing the stiffness length *l*<sub>s</sub> is equivalent to increasing fracture plane resistance (increasing *k*<sub>s</sub>). The SRM specimens are loaded with conditions such that the maximum principal stress *σ*<sub>a</sub> is vertical and the confinement stress *σ*<sub>l</sub> is isotropic in the lateral directions (*x* and *y* directions). Numerical parameters used to perform the numerical simulations are summarized in Table 1. The elastic properties of the rock matrix are selected from the granite properties of the Forsmark site in Sweden (Hakami et al. 2022; SKB 2008). DFN density and shear stiffness are selected over a wide range of values to perform sensitivity analysis.

<table><caption>**Table 1** Parameters of the SRM specimens</caption><thead><tr><td>Parameter</td><td>Symbol</td><td>Value</td></tr></thead><tbody><tr><td>Domain size</td><td>L</td><td>5 m</td></tr><tr><td>Mesh size</td><td>h</td><td>0.2 m</td></tr><tr><td>DFN percolation parameter</td><td>p</td><td>0.5–8</td></tr><tr><td>Young's modulus</td><td>E</td><td>76.9 GPa</td></tr><tr><td>Poisson's ratio</td><td>ν</td><td>0.23</td></tr><tr><td>Fracture stiffness length</td><td>l<sub>s</sub> = E/k<sub>s</sub></td><td>∞–0.2 m</td></tr><tr><td>Maximum principal stress</td><td>σ<sub>a</sub></td><td>3 kPa</td></tr><tr><td>Confinement stress</td><td>σ<sub>l</sub></td><td>1 kPa</td></tr></tbody></table>

## 2.3 Stress Fluctuations Quantification

The stress field in the fractured rock mass in response to applied tectonic stress is obtained by solving the linear elastic geomechanical problem. At any point, i.e., mesh element, the stress field is defined by a tensor and corresponding invariants. For example, Von Mises stress σ<sub>e</sub> quantifies the intensity of the deviatoric component of a stress tensor. It is defined from stress tensor σ̄ as:

$$ \sigma_e(\bar{\sigma}) = \frac{1}{\sqrt{2}} \sqrt{(\sigma_1 - \sigma_2)^2 + (\sigma_1 - \sigma_3)^2 + (\sigma_2 - \sigma_3)^2}, \quad (4) $$

with σ₁, σ₂ and σ₃ the principal stress components of the stress tensor σ̄. At any point *x* in the 3D cartesian space, σ<sub>e</sub>(*x*) refers to the Von Mises stress at position *x*.

It is possible to analyze distributions of local values, as well as deviation from the mean stress tensor (σ̄) in the volume *V*:

$$ \langle \bar{\sigma} \rangle = \frac{1}{V} \int_V \bar{\sigma}(x) dV, \quad (5) $$

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Numerical setup used to perform simulation from DFN generation (left), block domain discretization (center), and stress field resolution (right). The initial cubic sample is cut by the DFN into smaller</figcaption></figure>

blocks, meshed with a size *h* that must be much smaller than the average spacing between fractures
which is approximately equal to the applied remote stress (Gao et al. 2017).

Most studies characterize stress perturbation by analyzing separately the principal stress magnitudes and orientation variations (Hakala et al. 2019; Hakami 2006; Valli et al. 2011). However, given the tensorial nature of the stress, both analyses are not independent. Gao and Harrison (2016) proposed to have a single descriptor of the stress fluctuations by calculating the difference between the local and mean tensors. This is the approach we have implemented by calculating the local stress dispersion $D_e(x)$ as defined by Gao and Harrison (2018):

$$D_e(x) = \left\| \overline{\overline{\sigma}}(x) - \langle \overline{\overline{\sigma}} \rangle \right\|_F, \qquad (6)$$

where $\| \cdot \|_F$ refers to the Frobenius norm (also called Euclidean norm). Global stress dispersion $D_E$ is then defined as the quadratic mean of local stress dispersion over the whole volume:

$$D_E = \sqrt{\frac{1}{V} \int_V D_e(x)^2 dV}. \qquad (7)$$

To our knowledge, most numerical studies of fracture-induced stress perturbation based on this tensorial approach are made in 2D (Khodaei et al. 2020, 2021a; b; Lei and Gao 2018), which reduces the complexity of spatial organization and anisotropy inherent to fracture networks. It has been shown that stress variability depends on applied stress ratio, fracture and matrix mechanical parameters, fracture spatial organization and connectivity. Lei and Gao (2018) show that stress variability is more dominated by matrix resistance if fractures are disconnected, but more dependent on frictional sliding of fractures if the system is well connected.

If the global stress dispersion quantifies the intensity of stress fluctuations at the network scale, we need another indicator to quantify its localization. Depending on fracture network properties, the percentage of rock volume affected by this stress fluctuations may vary. To measure stress fluctuations localization, we use the participation ratio, inspired from (Davy et al. 1995; Edwards and Thouless 1972; Maillot et al. 2016), and defined as:

$$dQ = \frac{1}{\Omega} \cdot \frac{(\int_{\Omega} Q d\Omega)^2}{(\int_{\Omega} Q^2 d\Omega)}, \qquad (8)$$

where $Q$ is a measure of local stress perturbation (either stress dispersion or Von Mises stress), and $\Omega$ the integration volume. We first define $dD_e$ the stress dispersion participation ratio, with $Q = D_e$ and $\Omega = V$, the SRM volume. As local stress dispersion $D_e$ is positive, it is impossible to know if the deviation from applied remote stress corresponds to a

stress enhancement or diminution. We then introduce two Von Mises stress participation ratio $d\sigma_e^+$ and $d\sigma_e^-$, where the quantity $Q = \sigma_e(x) - \sigma_e(\langle \overline{\overline{\sigma}} \rangle)$ is integrated over $V^+$ and $V^-$ respectively, which corresponds to the volumes where the Von Mises stress is respectively larger and lower than the Von Mises stress of the mean stress tensor $\sigma_e(\langle \overline{\overline{\sigma}} \rangle)$. The participation ratio is a measure of the percentage of volume that is affected by the stress perturbation. If the deviation from applied remote stress is perfectly homogeneous in the volume, the ratio is equal to 1. On the other hand, if the stress perturbation is localized in a small volume, the participation ratio will tend to zero as the contribution of the perturbed volume to the total volume.

# 3 Results

We first perform a sensitivity analysis on fracture parameters at the single fracture scale, before going to the network scale by performing numerical simulations on the 44 DFN realizations described in Sect. 2.1.

## 3.1 Stress Fluctuations Induced by an Isolated Fracture

We consider an isolated disk-shaped fracture of diameter $l$, positioned at the center of the cubic domain of size $L$, with a strike direction aligned with $\bar{y}$ axis and tilted with a dip angle $\theta$ (Fig. 3a). Maximum principal and confining stresses $\sigma_a$ and $\sigma_l$ are applied vertically and horizontally, respectively, to the specimen boundaries. The stress field is solved in the matrix, and Von Mises stress $\sigma_e$ and local stress dispersion $D_e$ are computed in each mesh element of the model. The relative displacement between blocks on opposite sides of the fracture induces a stress perturbation around the fracture. Figure 3b shows that the Von Mises stress is maximum near the fracture tip while it decreases around the central part of the fracture, known as shadow effect. The Von Mises stress measured far from the fracture is $\sigma_a - \sigma_l = 2$ kPa (green color in Fig. 3b). Figure 3c shows the local stress dispersion field that is null far from the fracture, and positive around the fracture. It is not possible to know if the deviation from the applied remote stress corresponds to a stress enhancement or diminution when only looking at the local stress dispersion.

Figure 4a, c show the probability density functions of normalized Von Mises stress $\sigma_e^*$ (normalized by the Von Mises value of the applied remote stress) for various fracture size $l$ and stiffness length $l_s$ for the isolated fracture case with $\theta = 45^\circ$, in a log-linear representation. The distribution is a two-tails distribution, dominated by its central part around $\sigma_e^* = 1$, which reflects the limited volume of
which is approximately equal to the applied remote stress (Gao et al. 2017).

Most studies characterize stress perturbation by analyzing separately the principal stress magnitudes and orientation variations (Hakala et al. 2019; Hakami 2006; Valli et al. 2011). However, given the tensorial nature of the stress, both analyses are not independent. Gao and Harrison (2016) proposed to have a single descriptor of the stress fluctuations by calculating the difference between the local and mean tensors. This is the approach we have implemented by calculating the local stress dispersion $D_e(x)$ as defined by Gao and Harrison (2018):

$$D_e(x) = \left\| \overline{\overline{\sigma}}(x) - \langle \overline{\overline{\sigma}} \rangle \right\|_F, \qquad (6)$$

where $\| \cdot \|_F$ refers to the Frobenius norm (also called Euclidean norm). Global stress dispersion $D_E$ is then defined as the quadratic mean of local stress dispersion over the whole volume:

$$D_E = \sqrt{\frac{1}{V} \int_V D_e(x)^2 dV}. \qquad (7)$$

To our knowledge, most numerical studies of fracture-induced stress perturbation based on this tensorial approach are made in 2D (Khodaei et al. 2020, 2021a; b; Lei and Gao 2018), which reduces the complexity of spatial organization and anisotropy inherent to fracture networks. It has been shown that stress variability depends on applied stress ratio, fracture and matrix mechanical parameters, fracture spatial organization and connectivity. Lei and Gao (2018) show that stress variability is more dominated by matrix resistance if fractures are disconnected, but more dependent on frictional sliding of fractures if the system is well connected.

If the global stress dispersion quantifies the intensity of stress fluctuations at the network scale, we need another indicator to quantify its localization. Depending on fracture network properties, the percentage of rock volume affected by this stress fluctuations may vary. To measure stress fluctuations localization, we use the participation ratio, inspired from (Davy et al. 1995; Edwards and Thouless 1972; Maillot et al. 2016), and defined as:

$$dQ = \frac{1}{\Omega} \cdot \frac{(\int_{\Omega} Q d\Omega)^2}{(\int_{\Omega} Q^2 d\Omega)}, \qquad (8)$$

where $Q$ is a measure of local stress perturbation (either stress dispersion or Von Mises stress), and $\Omega$ the integration volume. We first define $dD_e$ the stress dispersion participation ratio, with $Q = D_e$ and $\Omega = V$, the SRM volume. As local stress dispersion $D_e$ is positive, it is impossible to know if the deviation from applied remote stress corresponds to a

stress enhancement or diminution. We then introduce two Von Mises stress participation ratio $d\sigma_e^+$ and $d\sigma_e^-$, where the quantity $Q = \sigma_e(x) - \sigma_e(\langle \overline{\overline{\sigma}} \rangle)$ is integrated over $V^+$ and $V^-$ respectively, which corresponds to the volumes where the Von Mises stress is respectively larger and lower than the Von Mises stress of the mean stress tensor $\sigma_e(\langle \overline{\overline{\sigma}} \rangle)$. The participation ratio is a measure of the percentage of volume that is affected by the stress perturbation. If the deviation from applied remote stress is perfectly homogeneous in the volume, the ratio is equal to 1. On the other hand, if the stress perturbation is localized in a small volume, the participation ratio will tend to zero as the contribution of the perturbed volume to the total volume.

# 3 Results

We first perform a sensitivity analysis on fracture parameters at the single fracture scale, before going to the network scale by performing numerical simulations on the 44 DFN realizations described in Sect. 2.1.

## 3.1 Stress Fluctuations Induced by an Isolated Fracture

We consider an isolated disk-shaped fracture of diameter $l$, positioned at the center of the cubic domain of size $L$, with a strike direction aligned with $\bar{y}$ axis and tilted with a dip angle $\theta$ (Fig. 3a). Maximum principal and confining stresses $\sigma_a$ and $\sigma_l$ are applied vertically and horizontally, respectively, to the specimen boundaries. The stress field is solved in the matrix, and Von Mises stress $\sigma_e$ and local stress dispersion $D_e$ are computed in each mesh element of the model. The relative displacement between blocks on opposite sides of the fracture induces a stress perturbation around the fracture. Figure 3b shows that the Von Mises stress is maximum near the fracture tip while it decreases around the central part of the fracture, known as shadow effect. The Von Mises stress measured far from the fracture is $\sigma_a - \sigma_l = 2$ kPa (green color in Fig. 3b). Figure 3c shows the local stress dispersion field that is null far from the fracture, and positive around the fracture. It is not possible to know if the deviation from the applied remote stress corresponds to a stress enhancement or diminution when only looking at the local stress dispersion.

Figure 4a, c show the probability density functions of normalized Von Mises stress $\sigma_e^*$ (normalized by the Von Mises value of the applied remote stress) for various fracture size $l$ and stiffness length $l_s$ for the isolated fracture case with $\theta = 45^\circ$, in a log-linear representation. The distribution is a two-tails distribution, dominated by its central part around $\sigma_e^* = 1$, which reflects the limited volume of
<figure><img src="image_4.png" /><figcaption>**Fig. 3** a Numerical setup for the isolated fracture case, and b 2D slice of Von Mises stress (rainbow color scale is indicated) and c stress dispersion (cool-warm color scale is indicated) fields. 2D slices cut view</figcaption></figure>

is defined from cutting plane perpendicular to the y axis at the center of the model

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Probability density functions of a, c normalized Von Mises stress and b, d normalized local stress dispersion for different a, b percolation parameters and c, d stiffness lengths, in the case of an isolated fracture</figcaption></figure>

<figure><img src="image_6.png" /></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 3** a Numerical setup for the isolated fracture case, and b 2D slice of Von Mises stress (rainbow color scale is indicated) and c stress dispersion (cool-warm color scale is indicated) fields. 2D slices cut view</figcaption></figure>

is defined from cutting plane perpendicular to the y axis at the center
of the model

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Probability density functions of a, c normalized Von Mises stress and b, d normalized local stress dispersion for different a, b percolation parameters and c, d stiffness lengths, in the case of an isolated fracture</figcaption></figure>

<figure><img src="image_7.png" /></figure>
influence around the fractures. Values where $\sigma_e^* < 1$ refer to the shadow zone around the fracture, while, values where $\sigma_e^* > 1$ refer to the stress increase at the fracture tips. Figure 4a shows that the larger the fracture size $l$, the higher the two tails of the distribution, the larger the stress perturbation. Note that stress shadowing cannot reduce the Von Mises stress by more than half of the applied Von Mises remote stress, while stress increase seems to have no limit (up to $3\sigma_e^*$ for $l = 2.4$ m). Figure 4c shows that the smaller $l_s$, the smaller the two tails of the distribution (distribution is dominated by its central part). Indeed, fracture shear stiffness limits the fracture plane displacement, which limits the stress perturbation induced by the fracture.

Figure 4b, d shows the normalized local stress dispersion $D_e^*$ (normalized by the Frobenius norm of remote

stress) for various fracture size $l$ and stiffness length $l_s$ for the isolated fracture case, in a log-linear representation. The local stress dispersion measured far from the fracture (blue color in Fig. 3c) is null. The distribution is dominated by these unperturbed zones (distribution peak is at $D_e^* = 0$) and decreases for large $D_e^*$ (the larger $D_e^*$, the smaller the probability of occurrence). The larger the fracture size $l$, the larger the stress perturbation, and the more large $D_e^*$ values are represented in the distribution. On the other hand, the smaller $l_s$, the smaller the stress perturbation, and the less large $D_e^*$ values are represented in the distribution.

We perform a sensitivity analysis to identify the dependency of global stress dispersion $D_E$ on fracture size $l$, dip angle $\theta$, stiffness length $l_s$ and applied stress ratio $\sigma_a/\sigma_l$ (Fig. 5). Global stress dispersion increases with fracture size

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Evolution of global stress dispersion $D_E$ with a fracture size $l$, b fracture length ratio $l/l_s$, c fracture dip angle $\theta$, and d applied stress ratio $\sigma_a/\sigma_l$. Numerical results (dots) are correctly matched by Equ. 9 (dashed lines)</figcaption></figure>
influence around the fractures. Values where $\sigma_e^* < 1$ refer to the shadow zone around the fracture, while, values where $\sigma_e^* > 1$ refer to the stress increase at the fracture tips. Figure 4a shows that the larger the fracture size $l$, the higher the two tails of the distribution, the larger the stress perturbation. Note that stress shadowing cannot reduce the Von Mises stress by more than half of the applied Von Mises remote stress, while stress increase seems to have no limit (up to $3\sigma_e^*$ for $l = 2.4$ m). Figure 4c shows that the smaller $l_s$, the smaller the two tails of the distribution (distribution is dominated by its central part). Indeed, fracture shear stiffness limits the fracture plane displacement, which limits the stress perturbation induced by the fracture.

Figure 4b, d shows the normalized local stress dispersion $D_e^*$ (normalized by the Frobenius norm of remote

stress) for various fracture size $l$ and stiffness length $l_s$ for the isolated fracture case, in a log-linear representation. The local stress dispersion measured far from the fracture (blue color in Fig. 3c) is null. The distribution is dominated by these unperturbed zones (distribution peak is at $D_e^* = 0$) and decreases for large $D_e^*$ (the larger $D_e^*$, the smaller the probability of occurrence). The larger the fracture size $l$, the larger the stress perturbation, and the more large $D_e^*$ values are represented in the distribution. On the other hand, the smaller $l_s$, the smaller the stress perturbation, and the less large $D_e^*$ values are represented in the distribution.

We perform a sensitivity analysis to identify the dependency of global stress dispersion $D_E$ on fracture size $l$, dip angle $\theta$, stiffness length $l_s$ and applied stress ratio $\sigma_a/\sigma_l$ (Fig. 5). Global stress dispersion increases with fracture size

<figure><img src="image_8.png" /><figcaption>**Fig. 5** Evolution of global stress dispersion $D_E$ with a fracture size $l$, b fracture length ratio $l/l_s$, c fracture dip angle $\theta$, and d applied stress ratio $\sigma_a/\sigma_l$. Numerical results (dots) are correctly matched by Equ. 9 (dashed lines)</figcaption></figure>
$l$, while it decreases with the fracture size ratio $l/l_s$. Also, it depends on the fracture orientation and applied stress ratio as the applied force on the fracture plane depends on the relative orientation of the fracture with respect to the applied stress tensor. The fixed values when tuning a parameter for the sensitivity analysis are: fracture size $l = 2$ m, stiffness length $l_s = \infty$, fracture dip angle $\theta = 45^\circ$ and applied stress ratio $\sigma_a/\sigma_l = 3$. Based on this sensitivity analysis, we propose an analytical solution for the global stress dispersion $D_{E,f}$ induced by a uniformly loaded disk-shaped fracture:

$$D_{E,f}(l, \theta, l_s) = \frac{3\pi (\sigma_a - \sigma_l) \cos(\theta) \sin(\theta)}{8 (1 + \frac{l}{l_s}) \sqrt{\frac{4}{3} \pi \left(\frac{l}{2}\right)^3 / L^3}} \quad (9)$$

## 3.2 Correlations Between Fracture Network Properties and Stress Fluctuations

Considering the DFN models defined in Sect. 2.1, numerical simulations are performed according to parameters summarized in Table 1. Maximum principal and confining stresses $\sigma_a$ and $\sigma_l$ are applied vertically and horizontally, respectively, to the specimen boundaries. The stress field is solved in the matrix, and Von Mises stress $\sigma_e$ and local stress dispersion $D_e$ are computed in each mesh element of the model. Figure 6 shows the probability density functions of normalized Von Mises and normalized local stress dispersion in a linear plot, considering percolation parameters from low to high for the a4ls$\infty$ models, and different stiffness lengths $l_s$ for the a4p4 models.

Figure 6a shows that the larger the fracture density (the larger the DFN percolation parameter), the wider the

<figure><img src="image_8.png" /><figcaption>Fig. 6 Probability density functions of a, b normalized Von Mises stress and c, d normalized local stress dispersion, as a function of a, c percolation parameter $p$ for the a4ls$\infty$ models and b, d stiffness length $l_s$ for the a4p4 models</figcaption></figure>
$l$, while it decreases with the fracture size ratio $l/l_s$. Also, it depends on the fracture orientation and applied stress ratio as the applied force on the fracture plane depends on the relative orientation of the fracture with respect to the applied stress tensor. The fixed values when tuning a parameter for the sensitivity analysis are: fracture size $l = 2$ m, stiffness length $l_s = \infty$, fracture dip angle $\theta = 45^\circ$ and applied stress ratio $\sigma_a/\sigma_l = 3$. Based on this sensitivity analysis, we propose an analytical solution for the global stress dispersion $D_{E,f}$ induced by a uniformly loaded disk-shaped fracture:

$$D_{E,f}(l, \theta, l_s) = \frac{3\pi (\sigma_a - \sigma_l) \cos(\theta) \sin(\theta)}{8 \left(1 + \frac{l}{l_s}\right) \sqrt{\frac{4}{3} \pi \left(\frac{l}{2}\right)^3 L^3} \quad (9)$$

## 3.2 Correlations Between Fracture Network Properties and Stress Fluctuations

Considering the DFN models defined in Sect. 2.1, numerical simulations are performed according to parameters summarized in Table 1. Maximum principal and confining stresses $\sigma_a$ and $\sigma_l$ are applied vertically and horizontally, respectively, to the specimen boundaries. The stress field is solved in the matrix, and Von Mises stress $\sigma_e$ and local stress dispersion $D_e$ are computed in each mesh element of the model. Figure 6 shows the probability density functions of normalized Von Mises and normalized local stress dispersion in a linear plot, considering percolation parameters from low to high for the a4ls$\infty$ models, and different stiffness lengths $l_s$ for the a4p4 models.

Figure 6a shows that the larger the fracture density (the larger the DFN percolation parameter), the wider the

<figure><img src="image_9.png" /><figcaption>Fig. 6 Probability density functions of a, b normalized Von Mises stress and c, d normalized local stress dispersion, as a function of a, c percolation parameter $p$ for the a4ls$\infty$ models and b, d stiffness length $l_s$ for the a4p4 models</figcaption></figure>
distribution of Von Mises stresses. Also, for high percolation parameters, the distribution mode (i.e. the most probable value) is larger than the Von Mises stress of the applied remote stress; the distribution is also asymmetric with respect to its mode. The same effect is observed on Fig. 6b, with increasing $l_s$ at constant fracture density.

We next compare the relative variations of the distributions with increasing fracture density. To do this, we normalize the distribution by the sum of the influence volume of each fracture, which is approximately a sphere of diameter similar to the fracture diameter. This normalization volume is proportional to a percolation parameter, as it is defined as the total excluded volume around fractures per unit volume. If interactions between fractures can be neglected, the differences between the distributions should be explained by a change of the affected volume only, and normalized

distribution should overlap. Normalized distributions are plotted in Fig. 7a, c. The distribution tails for low percolation parameters ($p < 2$) overlap and are approximately equal to the sum of the perturbations induced by the different fractures taken independently of each other, meaning that stress interactions are negligible over this fracture density range. Figure 7b, d shows the normalized probability density function of Von Mises stress and stress dispersion for all $p4ls∞$ models. The distributions overlap, meaning that regardless of the fracture size distribution, only the excluded volume occupied by fractures determines the distribution of stresses in the surrounding rock matrix.

We now calculate the global dispersion considering that fracture stress interactions are negligible, meaning that fractures have independent contributions to the stress fluctuations. Under these assumptions, the global stress dispersion

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Normalized probability density functions of normalized Von Mises stress (a, b) and normalized local stress dispersion (c, d), as a function of percolation parameter $p$ for $p4ls∞$ models (a, c) and DFN</figcaption></figure>

type for $p4ls∞$ models (b, d). Distributions are plotted in linear-log axis to facilitate curves comparison
distribution of Von Mises stresses. Also, for high percolation parameters, the distribution mode (i.e. the most probable value) is larger than the Von Mises stress of the applied remote stress; the distribution is also asymmetric with respect to its mode. The same effect is observed on Fig. 6b, with increasing $l_s$ at constant fracture density.

We next compare the relative variations of the distributions with increasing fracture density. To do this, we normalize the distribution by the sum of the influence volume of each fracture, which is approximately a sphere of diameter similar to the fracture diameter. This normalization volume is proportional to a percolation parameter, as it is defined as the total excluded volume around fractures per unit volume. If interactions between fractures can be neglected, the differences between the distributions should be explained by a change of the affected volume only, and normalized

distribution should overlap. Normalized distributions are plotted in Fig. 7a, c. The distribution tails for low percolation parameters ($p < 2$) overlap and are approximately equal to the sum of the perturbations induced by the different fractures taken independently of each other, meaning that stress interactions are negligible over this fracture density range. Figure 7b, d shows the normalized probability density function of Von Mises stress and stress dispersion for all $p4ls∞$ models. The distributions overlap, meaning that regardless of the fracture size distribution, only the excluded volume occupied by fractures determines the distribution of stresses in the surrounding rock matrix.

We now calculate the global dispersion considering that fracture stress interactions are negligible, meaning that fractures have independent contributions to the stress fluctuations. Under these assumptions, the global stress dispersion

<figure><img src="image_10.png" /><figcaption>**Fig. 7** Normalized probability density functions of normalized Von Mises stress (a, b) and normalized local stress dispersion (c, d), as a function of percolation parameter $p$ for $p4ls∞$ models (a, c) and DFN</figcaption></figure>

type for $p4ls∞$ models (b, d). Distributions are plotted in linear-log axis to facilitate curves comparison
to the square is additive, i.e., it is equal to the sum of the contribution of each individual fracture loaded by the applied remote stress.

We derive an equation for stress dispersion, noted $\tilde{D}_E$, from an analytical reasoning based on fracture density distribution $n(l, \theta)$ where $l$ and $\theta$ refer to the fracture size and orientation, respectively:

$$ \tilde{D}_E^2 = \iint_{l,\theta} n(l, \theta) D_{E,f}^2(l, \theta, l_s) dld\theta. \quad (10) $$

If fracture size and orientation distributions $n(l)$ and $n(\theta)$, respectively, are not correlated (i.e., $n(l, \theta) = n(l) \cdot n(\theta)$), hence Equ. 10 becomes:

$$ \tilde{D}_E^2 = F_\theta \int_{l_{\min}}^{l_{\max}} \frac{n(l)l^3}{(1+l/l_s)^2} dl, \quad (11) $$

where $F_\theta = \int_\theta \tau^2(\theta)n(\theta)d\theta$, is an orientation factor. If fracture size is much smaller than stiffness length ($l \ll l_s$), the network stress dispersion is proportional to the square root of the percolation parameter:

$$ \tilde{D}_E \propto \sqrt{F_\theta p}. \quad (12) $$

Otherwise, it depends on the fracture size distribution. For constant size network, Equ. 11 becomes:

$$ \tilde{D}_E^2 = F_\theta \cdot \frac{N_f}{V} \cdot \frac{l^3}{(1+l/l_s)^2}, \quad (13) $$

with $N_f$ the number of fractures contained in the volume $V$. For power law fracture size distribution, the network global stress dispersion is proportional to the power law exponent of the fracture size distribution. If $a = 3$ or $a = 4$, then Equ. 11 becomes respectively:

$$ \tilde{D}_E^2 = \alpha F_\theta \left( \frac{l_s}{1+l_{\min}/l_s} - \frac{l_s}{1+l_{\max}/l_s} \right). \quad (14) $$

$$ \tilde{D}_E^2 = \alpha F_\theta \left( \ln \left( \frac{l_{\max}}{l_{\min}} \right) + \ln \left( \frac{l_s + l_{\min}}{l_s + l_{\max}} \right) + \frac{1}{1+l_{\max}/l_s} - \frac{1}{1+l_{\min}/l_s} \right). \quad (15) $$

We also compute the global stress dispersion $D_E^*$ for the set of DFN realizations (Sect. 2.1) by applying the principle of superposition of the fracture contribution:

$$ D_E^{*2} = \sum_f D_{E,f}^2, \quad (16) $$

where $D_{E,f}$ is the global dispersion generated by the fracture $f$. To estimate the deviation of the actual stress dispersion

$D_E$ from $D_E^*$, we define $\Delta D_E = (D_E - D_E^*)/D_E$, as an indicator of the stress interactions between fractures.

Figure 8a shows that, for $l_s = \infty$, the global stress dispersion $D_E$ evolves as the square root of the percolation parameter, in agreement with Equ. 12 (dashed line in Fig. 8a). The principle of superposition correctly approximates the global stress dispersion for low percolation parameter ($\Delta D_E \approx 0$) but underestimates it when $p > 2.5$ (Fig. 8b). This emphasizes the role of fracture connectivity in increasing the stress fluctuations, showing that fracture clusters can be considered as “meta-fractures”, and produce more stress fluctuations than the sum of the fractures that compose them. Figure 8c shows that global stress dispersion decreases with shear stiffness $k_s$ in a way that depends on fracture size distribution, as shown by Equs. 13–15 that are reported as dashed lines.

Figure 9a shows the evolution of stress dispersion and Von Mises stress participation ratio with percolation parameter in the case where $l_s = \infty$. According to Equ. 8, the participation ratio measures the extent of stress perturbations, by computing the percentage of affected volume. For low percolation parameter, the participation ratios tend to 0, meaning that the stress perturbation is very localized in a small volume around the fracture themselves. The larger the percolation, the larger the participation ratio, until a plateau is reached. For percolation parameter much larger than the percolation threshold ($p \gg p_c$), the whole volume is statistically filled with intersecting fractures. Hence, adding new fractures in the system increases the stress perturbation intensity (measured by $D_E$), but not its localization. Figure 9b shows the volume ratios that is concerned by a stress increase or decrease, respectively noted $V^+$ and $V^-$. For low fracture density, the volume $V^+$ (respectively $V^-$), which corresponds to the volumes where the Von Mises stress is larger (respectively lower) than the Von Mises stress of the mean stress tensor, slightly decreases (respectively increases) with percolation parameter until a plateau is reached around 70% (respectively 30%) of the total volume, reached at $p \sim 2.5$. Also, both participation ratio and volume ratio are independent of the fracture type, meaning that stress perturbation localization is independent of the fracture size distribution.

## 4 Discussion

Fractured rock masses typically consist of many discontinuities over a wide range of fracture sizes (Bonnet et al. 2001), which poses numerical problem when building geomechanical models based on an explicit representation of the fractures. Identifying fractures that are responsible for major stress fluctuations is of great interest to simplify geomechanical numerical models by modeling explicitly only a part of the fractured system (Wang and Lei 2021). Several
to the square is additive, i.e., it is equal to the sum of the contribution of each individual fracture loaded by the applied remote stress.

We derive an equation for stress dispersion, noted $\tilde{D}_E$, from an analytical reasoning based on fracture density distribution $n(l, \theta)$ where $l$ and $\theta$ refer to the fracture size and orientation, respectively:

$$\tilde{D}_E^2 = \iint_{l,\theta} n(l, \theta) D_{E,f}^2(l, \theta, l_s) dld\theta. \quad (10)$$

If fracture size and orientation distributions $n(l)$ and $n(\theta)$, respectively, are not correlated (i.e., $n(l, \theta) = n(l) \cdot n(\theta)$), hence Equ. 10 becomes:

$$\tilde{D}_E^2 = F_\theta \int_{l_{\min}}^{l_{\max}} \frac{n(l)l^3}{(1+l/l_s)^2} dl, \quad (11)$$

where $F_\theta = \int_\theta \tau^2(\theta)n(\theta)d\theta$, is an orientation factor. If fracture size is much smaller than stiffness length ($l \ll l_s$), the network stress dispersion is proportional to the square root of the percolation parameter:

$$\tilde{D}_E \propto \sqrt{F_\theta p}. \quad (12)$$

Otherwise, it depends on the fracture size distribution. For constant size network, Equ. 11 becomes:

$$\tilde{D}_E^2 = F_\theta \cdot \frac{N_f}{V} \cdot \frac{l^3}{(1 + l/l_s)^2}, \quad (13)$$

with $N_f$ the number of fractures contained in the volume $V$. For power law fracture size distribution, the network global stress dispersion is proportional to the power law exponent of the fracture size distribution. If $a = 3$ or $a = 4$, then Equ. 11 becomes respectively:

$$\tilde{D}_E^2 = \alpha F_\theta \left( \frac{l_s}{1 + l_{\min}/l_s} - \frac{l_s}{1 + l_{\max}/l_s} \right). \quad (14)$$

$$\tilde{D}_E^2 = \alpha F_\theta \left( \ln \left( \frac{l_{\max}}{l_{\min}} \right) + \ln \left( \frac{l_s + l_{\min}}{l_s + l_{\max}} \right) + \frac{1}{1 + l_{\max}/l_s} - \frac{1}{1 + l_{\min}/l_s} \right). \quad (15)$$

We also compute the global stress dispersion $D_E^*$ for the set of DFN realizations (Sect. 2.1) by applying the principle of superposition of the fracture contribution:

$$D_E^{*2} = \sum_f D_{E,f}^2, \quad (16)$$

where $D_{E,f}$ is the global dispersion generated by the fracture $f$. To estimate the deviation of the actual stress dispersion

$D_E$ from $D_E^*$, we define $\Delta D_E = (D_E - D_E^*)/D_E$, as an indicator of the stress interactions between fractures.

Figure 8a shows that, for $l_s = \infty$, the global stress dispersion $D_E$ evolves as the square root of the percolation parameter, in agreement with Equ. 12 (dashed line in Fig. 8a). The principle of superposition correctly approximates the global stress dispersion for low percolation parameter ($\Delta D_E \approx 0$) but underestimates it when $p > 2.5$ (Fig. 8b). This emphasizes the role of fracture connectivity in increasing the stress fluctuations, showing that fracture clusters can be considered as “meta-fractures”, and produce more stress fluctuations than the sum of the fractures that compose them. Figure 8c shows that global stress dispersion decreases with shear stiffness $k_s$ in a way that depends on fracture size distribution, as shown by Equs. 13–15 that are reported as dashed lines.

Figure 9a shows the evolution of stress dispersion and Von Mises stress participation ratio with percolation parameter in the case where $l_s = \infty$. According to Equ. 8, the participation ratio measures the extent of stress perturbations, by computing the percentage of affected volume. For low percolation parameter, the participation ratios tend to 0, meaning that the stress perturbation is very localized in a small volume around the fracture themselves. The larger the percolation, the larger the participation ratio, until a plateau is reached. For percolation parameter much larger than the percolation threshold ($p \gg p_c$), the whole volume is statistically filled with intersecting fractures. Hence, adding new fractures in the system increases the stress perturbation intensity (measured by $D_E$), but not its localization. Figure 9b shows the volume ratios that is concerned by a stress increase or decrease, respectively noted $V^+$ and $V^-$. For low fracture density, the volume $V^+$ (respectively $V^-$), which corresponds to the volumes where the Von Mises stress is larger (respectively lower) than the Von Mises stress of the mean stress tensor, slightly decreases (respectively increases) with percolation parameter until a plateau is reached around 70% (respectively 30%) of the total volume, reached at $p \sim 2.5$. Also, both participation ratio and volume ratio are independent of the fracture type, meaning that stress perturbation localization is independent of the fracture size distribution.

## 4 Discussion

Fractured rock masses typically consist of many discontinuities over a wide range of fracture sizes (Bonnet et al. 2001), which poses numerical problem when building geomechanical models based on an explicit representation of the fractures. Identifying fractures that are responsible for major stress fluctuations is of great interest to simplify geomechanical numerical models by modeling explicitly only a part of the fractured system (Wang and Lei 2021). Several
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Evolution of **a**, **c** global stress dispersion $D_E$ and **b**, **d** deviation from superposition principle $\Delta D_E$, with **a**, **c** percolation parameter for all $l_s \infty$ models and with **b**, $d = 1/l_s$ for all p4 models. Filled dots refer to global stress dispersion $D_E$ computed by numerical simulations,</figcaption></figure>

while empty dots refer to approximated global stress dispersion $D_E^*$ (neglecting interactions) for the corresponding DFN realizations. The dots colors refer to the DFN type. Dashed lines refer to analytical solutions of $\bar{D}_E$ (Eqs. 12–15).

studies have focused on understanding the link between DFN parameters and stress fluctuations (Khodaei et al. 2021b; Lei and Gao 2018), but none of them propose a direct quantification of this process. In this study, we propose analytical solutions of global stress dispersion derived from the principle of superposition in Eq. 10. We show that for percolation parameter larger than the percolation threshold, the principle of superposition tends to underestimate the global stress dispersion, as connected fractures induce larger fracture plane displacements. This is coherent with findings of Lei and Gao (2018) who state that stress variability is more dominated by matrix resistance if fractures are disconnected. Still, these analytical solutions can be used to identify the range of fracture sizes that are responsible for most of the

stress perturbation in the system, depending on the fracture density distribution. Figure 10 shows the global stress dispersion $D_E$ and the corresponding cumulative contribution, for fractures of sizes ranging from $l_{\min} = 1$ to an upper limit $x$, considering power law fracture size distribution exponent $a = 3$ and $a = 4$, and a stiffness length $l_s = \infty$ and $l_s = 10$. This shows the high impact of the fracture size distribution on the global stress dispersion. Considering $l_s = \infty$, i.e., frictionless fractures, one can see that all fractures scales have important contribution for $a = 4$. This exponent is characteristic of self-similar (i.e., statistically similar at all scales) three-dimensional fracture networks (Bour 2002). For $a = 3$, this contribution to stress dispersion increases with fracture size, because the proportion of 'large' fractures inducing
<figure><img src="image_11.png" /><figcaption>**Fig. 8** Evolution of **a**, **c** global stress dispersion $D_E$ and **b**, **d** deviation from superposition principle $\Delta D_E$, with **a**, **c** percolation parameter for all $l_s \infty$ models and with **b**, $d = 1/l_s$ for all p4 models. Filled dots refer to global stress dispersion $D_E$ computed by numerical simulations,</figcaption></figure>

while empty dots refer to approximated global stress dispersion $D_E^*$ (neglecting interactions) for the corresponding DFN realizations. The dots colors refer to the DFN type. Dashed lines refer to analytical solutions of $\bar{D}_E$ (Eqs. 12–15).

studies have focused on understanding the link between DFN parameters and stress fluctuations (Khodaei et al. 2021b; Lei and Gao 2018), but none of them propose a direct quantification of this process. In this study, we propose analytical solutions of global stress dispersion derived from the principle of superposition in Eq. 10. We show that for percolation parameter larger than the percolation threshold, the principle of superposition tends to underestimate the global stress dispersion, as connected fractures induce larger fracture plane displacements. This is coherent with findings of Lei and Gao (2018) who state that stress variability is more dominated by matrix resistance if fractures are disconnected. Still, these analytical solutions can be used to identify the range of fracture sizes that are responsible for most of the

stress perturbation in the system, depending on the fracture density distribution. Figure 10 shows the global stress dispersion $D_E$ and the corresponding cumulative contribution, for fractures of sizes ranging from $l_{\min} = 1$ to an upper limit $x$, considering power law fracture size distribution exponent $a = 3$ and $a = 4$, and a stiffness length $l_s = \infty$ and $l_s = 10$. This shows the high impact of the fracture size distribution on the global stress dispersion. Considering $l_s = \infty$, i.e., frictionless fractures, one can see that all fractures scales have important contribution for $a = 4$. This exponent is characteristic of self-similar (i.e., statistically similar at all scales) three-dimensional fracture networks (Bour 2002). For $a = 3$, this contribution to stress dispersion increases with fracture size, because the proportion of ‘large’ fractures inducing
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Evolution of **a** stress dispersion and Von Mises stress participation ratios and **b** affected volume ratio $V^+/V$ and $V^-/V$. The dots colors refer to the DFN type</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10** **a** Global stress dispersion $D_E$ and **b** corresponding cumulative contribution, for fractures of sizes ranging from $l_{min} = 1$ to an upper limit $x$, considering power law fracture size distribution exponent $a = 3$ and $a = 4$, and stiffness length $l_s = \infty$ and $l_s = 10$</figcaption></figure>

important stress perturbation also increases. Moreover, it shows that, for frictional fractures with constant stiffness length, most of the contribution is made by fractures smaller than stiffness length $l_s$ (considered as frictionless). The stiffness length thus represents a characteristic length above which the associated stress perturbation becomes limited. This shows that special attention should be paid to the parameters controlling the contribution of fractures to stress fluctuations to simplify geomechanical models, and that a criterion based only on fracture size may be too simplistic (Wang and Lei 2021).

# 5 Conclusion

In this paper, we study the way stress is redistributed around fractures and assess stress variability in fractured rock masses, considering various fracture size distributions, densities, and frictional properties, supported by 3D numerical simulations.

Looking at one isolated fracture embedded in an elastic isotropic matrix, a detailed sensitivity analysis on different parameters showed the dependence of stress fluctuations
<figure><img src="image_12.png" /><figcaption>**Fig. 9** Evolution of **a** stress dispersion and Von Mises stress participation ratios and **b** affected volume ratio $V^+/V$ and $V^-/V$. The dots colors refer to the DFN type</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10** **a** Global stress dispersion $D_E$ and **b** corresponding cumulative contribution, for fractures of sizes ranging from $l_{min} = 1$ to an upper limit $x$, considering power law fracture size distribution exponent $a = 3$ and $a = 4$, and stiffness length $l_s = \infty$ and $l_s = 10$</figcaption></figure>

important stress perturbation also increases. Moreover, it shows that, for frictional fractures with constant stiffness length, most of the contribution is made by fractures smaller than stiffness length $l_s$ (considered as frictionless). The stiffness length thus represents a characteristic length above which the associated stress perturbation becomes limited. This shows that special attention should be paid to the parameters controlling the contribution of fractures to stress fluctuations to simplify geomechanical models, and that a criterion based only on fracture size may be too simplistic (Wang and Lei 2021).

# 5 Conclusion

In this paper, we study the way stress is redistributed around fractures and assess stress variability in fractured rock masses, considering various fracture size distributions, densities, and frictional properties, supported by 3D numerical simulations.

Looking at one isolated fracture embedded in an elastic isotropic matrix, a detailed sensitivity analysis on different parameters showed the dependence of stress fluctuations
on fracture geometrical and mechanical parameters (size, orientation, friction, etc.). Based on this sensitivity analysis, an analytical solution for global stress dispersion induced by a uniformly loaded fracture is proposed. For simplicity and to develop a clear analytical framework, we consider the fracture resistance to be elastic. Extension to more complex frictional laws is straightforward, as we show that the stress perturbation induced by a fracture depends on the force applied to its plane, which can be easily recovered from the surrounding stress field and the fracture orientation. However, extension to more complex elastic-plastic constitutive models of the rock matrix is not as simple, since the superposition principle used in Eq. 10 is only applicable to elastic conditions (Timoshenko 1951).

When looking at the network scale, because the contribution of individual fractures is nearly a sphere surrounding the fractures, a good proxy for fracture stress perturbation is the percolation parameter. We show that the distributions of Von Mises stress and local stress dispersion are driven by the DFN percolation parameter, regardless of the fracture size distribution. Moreover, for low fracture density (below the percolation threshold) the global stress dispersion can be reasonably approximated by summing the contribution of each individual fracture. For percolation parameter larger than the percolation threshold, the principle of superposition tends to underestimate the global stress dispersion, as connected fractures induce larger fracture plane displacements. Still, approximation from the principle of superposition can be used as a lower bound of global stress dispersion. This represents a key step forward the understanding of stress fluctuations in fractured rock mass.

Also, stress dispersion decreases when decreasing the stiffness length $l_s$ in a way that depends on fracture size distribution. The stiffness length is identified here as a characteristic length for stress fluctuations. Identifying characteristic length scales is critical to predict stress fluctuations, as a fracture may redistribute stress differently depending on whether it is reactivated or not, which depends on its geometrical and mechanical parameters. When dealing with dense multiscale fracture networks, identifying which fractures are responsible for major stress perturbations is of great interest for numerical models, as modeling every fracture is far beyond the reach of current numerical capabilities. The use of dual-scale models, where only a part of the fracture system is modeled explicitly (the other one being modeled implicitly using effective properties) is of great interest to tackle this issue. The developed analytical solutions at the network scale thus represent a key step towards geomechanical model simplifications.

**Acknowledgements** This work was funded by the Swedish Nuclear Fuel and Waste Management Company, Svensk Kärnbränslehantering

AB (SKB), and the Nuclear Waste Management Organization (NWMO) in Toronto, Canada. The authors wish to thank Rima Ghazal for her help with the numerical models in 3DEC.

**Author Contributions** EL, PD, CD, HK and DMI conceived the idea of this study. EL developed scripts for 3DEC simulations. EL, PD, and CD performed analysis of numerical results. EL took the lead in writing manuscript, all authors providing critical feed-back.

**Funding** This work was funded by the Swedish Nuclear Fuel and Waste Management Company, Svensk Kärnbränslehantering AB (SKB), and the Nuclear Waste Management Organization (NWMO) in Toronto, Canada.

**Data availability** Not applicable.

## Declarations

**Conflict of interest** The authors declare that they have no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Balberg I (1985) Universal percolation-threshold limits in the continuum. Phys Rev B 31(6):4053

Balberg I, Anderson C, Alexander S, Wagner N (1984) Excluded volume and its relation to the onset of percolation. Phys Rev B 30(7):3933

Barton CA, Zoback MD (1994) Stress perturbations associated with active faults penetrated by boreholes: possible evidence for near-complete stress drop and a new technique for stress magnitude measurement. J Geophys Res Solid Earth 99(B5):9373-9390

Bonnet E, Bour O, Odling NE, Davy P, Main I, Cowie P, Berkowitz B (2001) Scaling of fracture systems in geological media. Rev Geophys 39(3):347-383

Bour O (2002) A statistical scaling model for fracture network geometry, with validation on a multiscale mapping of a joint network (Hornelen Basin, Norway). J Geophys Re 107(B6):1

Bour O, Davy P (1998) On the connectivity of three-dimensional fault networks. Water Resour Res 34(10):2611-2622

Darcel C, Davy P, Bour O, De Dreuzy J (2006) Discrete fracture network for the Forsmark site. Swedish Nuclear Fuel and Waste Management Co., London

Davy P (1993) On the frequency-length distribution of the San Andreas fault system. J Geophys Res Solid Earth 98(B7):12141-12151

Davy P, Sornette A, Sornette D (1990) Some consequences of a proposed fractal nature of continental faulting. Nature 348(6296):56-58

Davy P, Hansen A, Bonnet E, Zhang SZ (1995) Localization and fault growth in layered brittle-ductile systems: Implications for
on fracture geometrical and mechanical parameters (size, orientation, friction, etc.). Based on this sensitivity analysis, an analytical solution for global stress dispersion induced by a uniformly loaded fracture is proposed. For simplicity and to develop a clear analytical framework, we consider the fracture resistance to be elastic. Extension to more complex frictional laws is straightforward, as we show that the stress perturbation induced by a fracture depends on the force applied to its plane, which can be easily recovered from the surrounding stress field and the fracture orientation. However, extension to more complex elastic-plastic constitutive models of the rock matrix is not as simple, since the superposition principle used in Eq. 10 is only applicable to elastic conditions (Timoshenko 1951).

When looking at the network scale, because the contribution of individual fractures is nearly a sphere surrounding the fractures, a good proxy for fracture stress perturbation is the percolation parameter. We show that the distributions of Von Mises stress and local stress dispersion are driven by the DFN percolation parameter, regardless of the fracture size distribution. Moreover, for low fracture density (below the percolation threshold) the global stress dispersion can be reasonably approximated by summing the contribution of each individual fracture. For percolation parameter larger than the percolation threshold, the principle of superposition tends to underestimate the global stress dispersion, as connected fractures induce larger fracture plane displacements. Still, approximation from the principle of superposition can be used as a lower bound of global stress dispersion. This represents a key step forward the understanding of stress fluctuations in fractured rock mass.

Also, stress dispersion decreases when decreasing the stiffness length $l_s$ in a way that depends on fracture size distribution. The stiffness length is identified here as a characteristic length for stress fluctuations. Identifying characteristic length scales is critical to predict stress fluctuations, as a fracture may redistribute stress differently depending on whether it is reactivated or not, which depends on its geometrical and mechanical parameters. When dealing with dense multiscale fracture networks, identifying which fractures are responsible for major stress perturbations is of great interest for numerical models, as modeling every fracture is far beyond the reach of current numerical capabilities. The use of dual-scale models, where only a part of the fracture system is modeled explicitly (the other one being modeled implicitly using effective properties) is of great interest to tackle this issue. The developed analytical solutions at the network scale thus represent a key step towards geomechanical model simplifications.

**Acknowledgements** This work was funded by the Swedish Nuclear Fuel and Waste Management Company, Svensk Kärnbränslehantering

AB (SKB), and the Nuclear Waste Management Organization (NWMO) in Toronto, Canada. The authors wish to thank Rima Ghazal for her help with the numerical models in 3DEC.

**Author Contributions** EL, PD, CD, HK and DMI conceived the idea of this study. EL developed scripts for 3DEC simulations. EL, PD, and CD performed analysis of numerical results. EL took the lead in writing manuscript, all authors providing critical feed-back.

**Funding** This work was funded by the Swedish Nuclear Fuel and Waste Management Company, Svensk Kärnbränslehantering AB (SKB), and the Nuclear Waste Management Organization (NWMO) in Toronto, Canada.

**Data availability** Not applicable.

## Declarations

**Conflict of interest** The authors declare that they have no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Balberg I (1985) Universal percolation-threshold limits in the continuum. Phys Rev B 31(6):4053

Balberg I, Anderson C, Alexander S, Wagner N (1984) Excluded volume and its relation to the onset of percolation. Phys Rev B 30(7):3933

Barton CA, Zoback MD (1994) Stress perturbations associated with active faults penetrated by boreholes: possible evidence for near-complete stress drop and a new technique for stress magnitude measurement. J Geophys Res Solid Earth 99(B5):9373–9390

Bonnet E, Bour O, Odling NE, Davy P, Main I, Cowie P, Berkowitz B (2001) Scaling of fracture systems in geological media. Rev Geophys 39(3):347–383

Bour O (2002) A statistical scaling model for fracture network geometry, with validation on a multiscale mapping of a joint network (Hornelen Basin, Norway). J Geophys Re 107(B6):1

Bour O, Davy P (1998) On the connectivity of three-dimensional fault networks. Water Resour Res 34(10):2611–2622

Darcel C, Davy P, Bour O, De Dreuzy J (2006) Discrete fracture network for the Forsmark site. Swedish Nuclear Fuel and Waste Management Co., London

Davy P (1993) On the frequency-length distribution of the San Andreas fault system. J Geophys Res Solid Earth 98(B7):12141–12151

Davy P, Sornette A, Sornette D (1990) Some consequences of a proposed fractal nature of continental faulting. Nature 348(6296):56–58

Davy P, Hansen A, Bonnet E, Zhang SZ (1995) Localization and fault growth in layered brittle-ductile systems: Implications for
deformations of the continental lithosphere. J Geophys Res Solid Earth 100(B4):6281–6294

Davy P, Le Goc R, Darcel C, Bour O, de Dreuzy JR, Munier R (2010) A likely universal model of fracture scaling and its consequence for crustal hydromechanics. J Geophys Res Solid Earth 115(B10):1

Davy P, Darcel C, Le Goc R, Mas Ivars D (2018) Elastic properties of fractured rock masses with frictional properties and power law fracture size distributions. J Geophys Res Solid Earth 2018:1

De Dreuzy J-R, Davy P, Bour O (2000) Percolation parameter and percolation-threshold estimates for three-dimensional random ellipses with widely scattered distributions of eccentricity and size. Phys Rev E 62(5):5948–5952

Dershowitz WS, Herda HS (1992) Interpretation of fracture spacing and intensity. In: Proceedings the 33th US symposium on rock mechanics (USRMS), Santa Fe, New Mexico, American Rock Mechanics Association

Edwards J, Thouless D (1972) Numerical studies of localization in disordered systems. J Phys C Solid State Phys 5(8):807

Figueiredo B, Sjöberg J, Mattila J, Hakala M, Suikkanen J (2023) Analysis and determination of the stress field at the Olkiluoto site. In: Proceedings IOP conference series: earth and environmental science, vol 1124. IOP Publishing, London, p 012002

Gao K, Harrison JP (2016) Mean and dispersion of stress tensors using Euclidean and Riemannian approaches. Int J Rock Mecha Min Sci 85:165–173

Gao K, Harrison JP (2018) Scalar-valued measures of stress dispersion. Int J Rock Mech Min Sci 106:234–242

Gao K, Harrison JP, Lei Q, Latham J-P (2017) Investigating the relationship between far-field stress and local values of the stress tensor. Proc Eng 191:536–542

Hakala M, Ström J, Valli J, Juvani J (2019) Structural control on stress variability at Forsmark: SKB R-19-23. Swedish Nuclear Fuel Waste Management Company, Stockholm, Sweden

Hakami H (2006) Numerical studies on spatial variation of the in situ stress field at Forsmark—a further step. Site Descript Model Forsmark Stage 2:1

Hakami E, Mas Ivars D, Darcel C (2022) Methodology for rock mechanics modelling of the Forsmark site: Swedish Nuclear Fuel and Waste Management Co.

Healy D, Jones RR, Holdsworth RE (2006) Three-dimensional brittle shear fracturing by tensile crack interaction. Nature 439(7072):64–67

Homberg C, Hu J, Angelier J, Bergerat F, Lacombe O (1997) Characterization of stress perturbations near major fault zones: insights from 2-D distinct-element numerical modelling and field studies (Jura mountains). J Struct Geolo 19(5):703–718

Itasca Consulting Group I (2020) 3DEC—three-dimensional distinct element code, Ver.7.0: Minneapolis, Itasca

Jaeger JC, Cook NG, Zimmerman R (2009) Fundamentals of rock mechanics. Wiley, London

Kachanov M (1989) Three-dimensional problems of strongly interacting arbitrarily located penny-shaped cracks. Int J Fract 199:1

Kachanov M (2003) On the problems of crack interactions and crack coalescence. Int J Fract 120(3):537–543

Khodaei M, Delijani EB, Hajipour M, Karroubi K, Dehghan AN (2020) On dispersion of stress state and variability in fractured media: effect of fracture scattering. Model Earth Syst Environ 2020:1–8

Khodaei M, Delijani EB, Dehghan AN, Hajipour M, Karroubi K (2021a) Stress/strain variability in fractured media: a fracture geometric study. Geotech Geol Eng 2021:1–20

Khodaei M, Delijani EB, Hajipour M, Karroubi K, Dehghan AN (2021b) Analyzing the correlation between stochastic fracture

networks geometrical properties and stress variability: a rock and fracture parameters study. J Pet Explor Prod 11(2):685–702

Lei Q, Gao K (2018) Correlation between fracture network properties and stress variability in geological media. Geophys Res Lett 2018:1

Lei Q, Gao K (2019) A numerical study of stress variability in heterogeneous fractured rocks. Int J Rock Mech Min Sci 113:121–133

Maillot J, Davy P, Le Goc R, Darcel C, De Dreuzy J-R (2016) Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models. Water Resour Res 52(11):8526–8545

Martin CD (2007) Quantifying in situ stress magnitudes and orientations for Forsmark. Forsmark stage 2.2: Swedish Nuclear Fuel and Waste Management Co.

Martin CD, Chandler NA (1993) Stress heterogeneity and geological structures. Proc Int J Rock Mech Min Sci Geomech Abstracts 30:993–999

Mas Ivars D, Pierce ME, Darcel C, Reyes-Montes J, Potyondy DO, Young RP, Cundall PA (2011) The synthetic rock mass approach for jointed rock mass modelling. Int J Rock Mech Min Sci 48(2):219–244

McGarr A, Gay N (1978) State of stress in the earth's crust. Annu Rev Earth Planet Sci 6(1):405–436

McTigue DF, Mei CC (1981) Gravity-induced stresses near topography of small slope. J Geophys Res Solid Earth 86(B10):9268–9278

Olson JE (1993) Joint pattern development: effects of subcritical crack growth and mechanical crack interaction. J Geophys Res Solid Earth 98(B7):12251–12265

Savage WZ, Swolfs HS (1986) Tectonic and gravitational stress in long symmetric ridges and valleys. J Geophys Res Solid Earth 91(B3):3677–3685

Selroos I-O, Mas Ivars D, Munier R, Hartley L, Darcel C, Davy P, Trinchero P, Cottrell M, Follin S, Hermanson J, Pointe PL, Libby S (2022) Discrete fracture network modelling handbook volume I: concepts, data and interpretation methods. Swedish Nuclear Fuel and Waste Management Co., London

SKB (2008) Site description of Forsmark at completion of the site investigation phase SDM-Site Forsmark. TR-08-05, Swedish Nuclear Fuel and Waste Management Co., London

Thomas RN, Paluszny A, Zimmerman RW (2017) Quantification of fracture interaction using stress intensity factor variation maps. J Geophys Res Solid Earth 2017:1

Timoshenko S (1951) Theory of elasticity. Oxford University Press, Oxford

Valli J, Kuula H, Hakala M (2011) Modelling of the in situ stress state at Olkiluoto Site, Western Finland. Posiva Oy

Wang L, Lei Q (2021) A dual-scale fracture network model for computing hydro-mechanical properties of fractured rock. Comput Geotech 138:104357

Yale DP (2003) Fault and stress magnitude controls on variations in the orientation of in situ stress. Geol Soc Lond Spec Ations 209(1):55–64

Zoback ML (1992) First-and second-order patterns of stress in the lithosphere: the world stress map project. J Geophys Res Solid Earth 97(B8):11703–11728

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
deformations of the continental lithosphere. J Geophys Res Solid Earth 100(B4):6281–6294

Davy P, Le Goc R, Darcel C, Bour O, de Dreuzy JR, Munier R (2010) A likely universal model of fracture scaling and its consequence for crustal hydromechanics. J Geophys Res Solid Earth 115(B10):1

Davy P, Darcel C, Le Goc R, Mas Ivars D (2018) Elastic properties of fractured rock masses with frictional properties and power law fracture size distributions. J Geophys Res Solid Earth 2018:1

De Dreuzy J-R, Davy P, Bour O (2000) Percolation parameter and percolation-threshold estimates for three-dimensional random ellipses with widely scattered distributions of eccentricity and size. Phys Rev E 62(5):5948–5952

Dershowitz WS, Herda HS (1992) Interpretation of fracture spacing and intensity. In: Proceedings the 33th US symposium on rock mechanics (USRMS), Santa Fe, New Mexico, American Rock Mechanics Association

Edwards J, Thouless D (1972) Numerical studies of localization in disordered systems. J Phys C Solid State Phys 5(8):807

Figueiredo B, Sjöberg J, Mattila J, Hakala M, Suikkanen J (2023) Analysis and determination of the stress field at the Olkiluoto site. In: Proceedings IOP conference series: earth and environmental science, vol 1124. IOP Publishing, London, p 012002

Gao K, Harrison JP (2016) Mean and dispersion of stress tensors using Euclidean and Riemannian approaches. Int J Rock Mecha Min Sci 85:165–173

Gao K, Harrison JP (2018) Scalar-valued measures of stress dispersion. Int J Rock Mech Min Sci 106:234–242

Gao K, Harrison JP, Lei Q, Latham J-P (2017) Investigating the relationship between far-field stress and local values of the stress tensor. Proc Eng 191:536–542

Hakala M, Ström J, Valli J, Juvani J (2019) Structural control on stress variability at Forsmark: SKB R-19-23. Swedish Nuclear Fuel Waste Management Company, Stockholm, Sweden

Hakami H (2006) Numerical studies on spatial variation of the in situ stress field at Forsmark—a further step. Site Descript Model Forsmark Stage 2:1

Hakami E, Mas Ivars D, Darcel C (2022) Methodology for rock mechanics modelling of the Forsmark site: Swedish Nuclear Fuel and Waste Management Co.

Healy D, Jones RR, Holdsworth RE (2006) Three-dimensional brittle shear fracturing by tensile crack interaction. Nature 439(7072):64–67

Homberg C, Hu J, Angelier J, Bergerat F, Lacombe O (1997) Characterization of stress perturbations near major fault zones: insights from 2-D distinct-element numerical modelling and field studies (Jura mountains). J Struct Geolo 19(5):703–718

Itasca Consulting Group I (2020) 3DEC—three-dimensional distinct element code, Ver.7.0: Minneapolis, Itasca

Jaeger JC, Cook NG, Zimmerman R (2009) Fundamentals of rock mechanics. Wiley, London

Kachanov M (1989) Three-dimensional problems of strongly interacting arbitrarily located penny-shaped cracks. Int J Fract 199:1

Kachanov M (2003) On the problems of crack interactions and crack coalescence. Int J Fract 120(3):537–543

Khodaei M, Delijani EB, Hajipour M, Karroubi K, Dehghan AN (2020) On dispersion of stress state and variability in fractured media: effect of fracture scattering. Model Earth Syst Environ 2020:1–8

Khodaei M, Delijani EB, Dehghan AN, Hajipour M, Karroubi K (2021a) Stress/strain variability in fractured media: a fracture geometric study. Geotech Geol Eng 2021:1–20

Khodaei M, Delijani EB, Hajipour M, Karroubi K, Dehghan AN (2021b) Analyzing the correlation between stochastic fracture

networks geometrical properties and stress variability: a rock and fracture parameters study. J Pet Explor Prod 11(2):685–702

Lei Q, Gao K (2018) Correlation between fracture network properties and stress variability in geological media. Geophys Res Lett 2018:1

Lei Q, Gao K (2019) A numerical study of stress variability in heterogeneous fractured rocks. Int J Rock Mech Min Sci 113:121–133

Maillot J, Davy P, Le Goc R, Darcel C, De Dreuzy J-R (2016) Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models. Water Resour Res 52(11):8526–8545

Martin CD (2007) Quantifying in situ stress magnitudes and orientations for Forsmark. Forsmark stage 2.2: Swedish Nuclear Fuel and Waste Management Co.

Martin CD, Chandler NA (1993) Stress heterogeneity and geological structures. Proc Int J Rock Mech Min Sci Geomech Abstracts 30:993–999

Mas Ivars D, Pierce ME, Darcel C, Reyes-Montes J, Potyondy DO, Young RP, Cundall PA (2011) The synthetic rock mass approach for jointed rock mass modelling. Int J Rock Mech Min Sci 48(2):219–244

McGarr A, Gay N (1978) State of stress in the earth's crust. Annu Rev Earth Planet Sci 6(1):405–436

McTigue DF, Mei CC (1981) Gravity-induced stresses near topography of small slope. J Geophys Res Solid Earth 86(B10):9268–9278

Olson JE (1993) Joint pattern development: effects of subcritical crack growth and mechanical crack interaction. J Geophys Res Solid Earth 98(B7):12251–12265

Savage WZ, Swolfs HS (1986) Tectonic and gravitational stress in long symmetric ridges and valleys. J Geophys Res Solid Earth 91(B3):3677–3685

Selroos I-O, Mas Ivars D, Munier R, Hartley L, Darcel C, Davy P, Trinchero P, Cottrell M, Follin S, Hermanson J, Pointe PL, Libby S (2022) Discrete fracture network modelling handbook volume I: concepts, data and interpretation methods. Swedish Nuclear Fuel and Waste Management Co., London

SKB (2008) Site description of Forsmark at completion of the site investigation phase SDM-Site Forsmark. TR-08-05, Swedish Nuclear Fuel and Waste Management Co., London

Thomas RN, Paluszny A, Zimmerman RW (2017) Quantification of fracture interaction using stress intensity factor variation maps. J Geophys Res Solid Earth 2017:1

Timoshenko S (1951) Theory of elasticity. Oxford University Press, Oxford

Valli J, Kuula H, Hakala M (2011) Modelling of the in situ stress state at Olkiluoto Site, Western Finland. Posiva Oy

Wang L, Lei Q (2021) A dual-scale fracture network model for computing hydro-mechanical properties of fractured rock. Comput Geotech 138:104357

Yale DP (2003) Fault and stress magnitude controls on variations in the orientation of in situ stress. Geol Soc Lond Spec Ations 209(1):55–64

Zoback ML (1992) First-and second-order patterns of stress in the lithosphere: the world stress map project. J Geophys Res Solid Earth 97(B8):11703–11728

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.