Analyzing the correlation between stochastic fracture networks geometrical properties and stress variability: a rock and fracture parameters study

In this study, the correlation between geometric properties of the fracture network and stress variability in a fractured rock was studied. Initially, discrete fracture networks were generated using a stochastic approach, then, considering the tensorial nature of stress, the stress field under various tectonic stress conditions was determined using finite-difference method. Ultimately, stress data were analyzed using tensor-based mathematical relations. Subsequently, the effects of four parameters including rock tensile strength, rock cohesion, fracture normal stiffness and fracture dilation angle on the stress perturbation distribution were evaluated. The obtained results indicated that stress perturbation and dispersion are directly related to fracture density, which is expressed as the number of fractures per unit area utilizing the window sampling approach. It was also demonstrated that they are inversely related to power-law length exponent which represents the length of fracture. It was observed that stress distribution, among the evaluated parameters, is more sensitive to the fracture normal stiffness and the effects of rock parameters on stress distribution are negligible. It was concluded that the highest stress distribution is created when the fracture network is dense with fractures having high length and low normal stiffness value.


Introduction
One of the most crucial issues in rock mechanics studies and geomechanical topics in hydrocarbon reservoir engineering is the evaluation of in-situ stresses and factors affecting stress perturbation. These geomechanical topics include hydraulic fracturing, wellbore stability, determination of rock mass permeability and optimum mud weight, prevention of sand production, selection of appropriate strategies for well completion, reservoir production rate, as well as the evaluation of earthquake potential (Zoback 2007;Hudson et al. 1988;Latham et al. 2013;Matsumoto et al. 2015;Zang et al. 2009;Dehghan et al. 2017;Farsimadan et al. 2020). Moreover, it is important to consider the stress conditions for accurate prediction of the mechanical behavior of jointed rock mass ). Many studies have been conducted to measure in-situ stresses and research shows that fractures play an important role in the distribution and perturbation of tectonic stresses (Bruno et al. 1994;Day-Lewis et al. 2010;Rajabi et al. 2017;Schoenball et al. 2017;Hickman et al. 2004; Barton et al. 1994;Sahara et al. 2014;Khodaei et al. 2020). Bruno et al. (1994) showed significant azimuth changes of maximum horizontal stress in a reservoir with the depth and location of a subsurface structure using field data analytically and finite element modeling. Day-Lewis et al. (2010) investigated the direction of maximum horizontal compressive stress as a function of depth in two research wells near the San Andreas Fault in central and southern California. They found that the stress orientation shows the scale-invariant fluctuations at distances of 10 cm to several kilometers of the fault. The similarity between the scale of stress orientation fluctuations and the magnitude of earthquake frequency with the size of faults showed that these fluctuations are controlled by the stress perturbation that is due to the slip of faults in the crust under critical stress in the proximity of faults. Sahara et al. (2014) performed an analysis on the occurrence of breakout, breakout orientation and fracture characteristics. They observed that breakout orientation in the form of anomalies, centrally occurs adjacent to the fault cores and decreases with distance from the fault core. The pattern of breakout orientation in the proximity of natural fractures shows that the rotation of breakout relative to the direction of mean horizontal stress (S hmin ) is strongly influenced by the fracture orientation. Breakouts are also commonly found asymmetrically in zones with high fracture densities. In addition to the principal stresses heterogeneities, breakout heterogeneities are affected by mechanical heterogeneities, such as weak zones with different elastic modulus, rock strength and fracture patterns. Rajabi et al. (2017) performed the first analysis of tectonic stress in Clarence-Moreton Basin in New South Wales area of Australia. Their observations suggested that structures can play an important role in controlling stresses so that stress perturbation as a result of faults and fractures highly influences wellbore stability and permeability of reservoir rock, especially for safe and sustainable extraction of methane gas from gas reservoirs in this region.
Extensive research has been previously done on the generation of discrete fractures networks using stochastic approaches to determine the properties of jointed rock and their geomechanical and hydromechanical behavior (Min et al. 2004;Baghbanan et al. 2007Baghbanan et al. , 2008Li et al. 2019), as well as geomechanical modeling using finite difference method (Oda et al. 1993;Kobayashi et al. 2001;Rutqvist et al. 2013). Among the most important researches on generating a discrete fracture network using stochastic approaches is the generation of a two-dimensional stochastic fracture pattern using field data in the Sellafield area of the United Kingdom to determine the equivalent permeability tensor of the fractured rock mass by Min et al. (2004). Also, Baghbanan et al. (2007Baghbanan et al. ( , 2008, generated a discrete fracture network with the assumption that the geometry of the fractures are stochastically distributed in the network; therefore, they could investigate the hydraulic and hydromechanical behaviors that are dependent on the rock mass size and estimate the representative elementary volume. Concerning geomechanical continuum modeling, Oda et al. (1993) developed an elastic stress/strain relationship using crack tensor method to determine the influence of joints on the elastic behavior of rock mass in a   three-dimensional finite element model. Kobayashi et al. (2001) investigated the coupling of mechanical and hydraulic behavior of fractured rock mass during a hypothetical shaft sinking in the Sellafield area of the United Kingdom using a continuum approach (finite element method) and combining the crack tensor theory and the Barton and Bandis model. Rutqvist et al. (2013) investigated the coupling of geomechanical model and fluid flow in the fractured rock using a continuum model and crack tensor approach.
In these studies, the investigation of stress perturbation is based on an approach that utilizes scalar/vector formulations which analyzes the magnitude and orientation of principal stresses, individually. However, in nature, stress is a tensor and its magnitude and orientation must be taken into consideration simultaneously. As far as we are aware, the first use of a tensor-based approach in this matter can be credited to Lei et al. (2018). In this study, by combining numerical and mathematical analysis, local stress variability in fractured rocks was investigated under static stress loading conditions. Moreover, in previous studies, limited research has been done on the effect of different rock and fracture parameters on the stress field perturbation (regarding the geometrical properties of fracture). In this study, four different parameters including rock tensile strength, rock cohesion, fracture normal stiffness and fracture dilation angle were evaluated.

Research methodology
In this section, in order to determine the correlation between the geometric properties of fracture network (including fracture density and fracture length) and stress variability, as well as to investigate the effects of rock and fracture parameters on the distribution of local stress perturbation, a similar approach to the one that was proposed by Lei et al. (2018) was followed. They generated different DFN realizations based on the fracture intensity. Then, they used FEM-DEM approach as the numerical method and   (Figueiredo et al. 2015) Rock Elastic modulus tensor-based mathematical formulation to analyze the stress data. Their study aimed to investigate how the intensity and the connectivity of the fracture network would affect the stress variability. Of all the rock and fracture parameters, they studied only the effect of friction coefficient parameter on the stress distribution. In this work, the discrete fracture network was first generated by stochastic realization based on the fracture density. Then, the finite difference method in the continuum approach was utilized as the geomechanical modeling and numerical method. Subsequently, different realizations of the discrete fracture network were subjected to orthogonal tectonic stress loading. Finally, in order to analyze the stress, the correlations between the geometric properties of fracture and stress distribution and the effects of different parameters were evaluated using tensor-based equations. Figure 1 shows the methodologies and analyses used to achieve the goals of this research.

Stochastic approach in the generation of fracture network
In the stochastic approach, fractures are considered as straight lines in two-dimensional model and as planar discs/ polygons in three-dimensional model and the geometrical properties, such as location, frequency, length, orientation and fracture aperture are considered as dependent variables on the probability distribution of the outcrop mapping. Orientational data can be processed using Fisher, normal, or even uniform distribution functions (Einstein et al. 1983) and fracture lengths may display negative exponential, lognormal, gamma, or power-law distributions (Davy 1993;Bonnet et al. 2001).
In the present study, the fracture network length was determined using the power-law according to the following equations (Bonnet et al. 2001;Lei et al. 2016): whereas in this equation, n(l) is the number of fractures, a is the power-law length exponent, α is the density term, l is the fracture length and l min and l max are the smallest and largest fracture lengths.
In theory, in two-dimensional model, a is limited to the interval of [1, ∞), however extensive measurements based on navigation maps show that in natural fracture systems, a generally varies between 1.3 and 3.5 (Bonnet et al. 2001). The density term (i.e., α) is dependent on the total number of fractures in the system and is a function of different fracture orientations (Davy et al. 2010).
The fracture frequency can be described from two aspects of fracture density and fracture intensity and can be determined using the P ij system, where i is the dimension of the sample and j is the dimension of the measurement. In this paper, the fracture density (P 20 ) was determined using a window sampling approach that is defined as the number of fractures per unit area (Zeeb et al. 2013): In this equation, N is the number of fractures, A is the surface area and L is the size of the model.
In this study, using open-source software of Alghalandis Discrete Fracture Network Engineering (ADFNE) (Alghalandis 2014, 2016, 2017), a fracture network in the size of 1 m × 1 m was stochastically generated. Determination of the location and orientation of fractures were done using a uniform distribution function and Fisher distribution function completely randomly. The largest and smallest fracture lengths were calculated using the power-law model to be 50 m and 0.02 m, respectively. Considering five different values for the power-law length exponent from 1.5 to 3.5 and 3 fracture densities of 80 m −2 , 2 × 80 m −2 and 4 × 80 m −2 , a total of 15 different fracture network realizations were generated. Figure 2 illustrates the different realizations of the fracture networks randomly generated at different densities.

Numerical method and geomechanical parameters
In this study, FLAC 2D software (Itasca Consulting Group Inc 2017) was employed to determine the stress distribution in response to tectonic stress loading that included three types of orthogonal effective principal tectonic stress loading at stress ratios of SR = 1, SR = 2 and SR = 3 (Fig. 3). The most   important advantage of this software is the complete perseverance of not only the fracture dead-ends after meshing but also the backbone of the fracture. Figure 4 shows a comparison of the shape of the fracture network after meshing using finite difference method (FLAC 2D software) and solely the backbone.
For intact rock and fracture shear stress/strain behavior, the Mohr-Coulomb model was used. Characteristics of a limestone sample were considered as rock mechanical properties. Table 1  Fractures with material fillings can be considered as an equivalent model of solids in which the elastic modulus of fracture is calculated using the following equation: whereas in this equation, d is the size of the mesh element. Figueiredo et al. (2015) developed a simple model with a vertical fracture to validate the use of mesh element size and quantity and whether it has the ability to estimate the stress correctly or not. It was found that the difference between the model made with the analytical solution was less than 5% and this model is quite suitable for the calculation of stress distribution and stress concentration around the fracture.
In this paper, the model size is considered to be 1 m × 1 m with 40,000 meshes per square meter (200 × 200 mesh) and with the mesh element length of d = 0.5 cm in accordance with the geometry and model in the study of Figueiredo et al. (2015).

Determination of perturbation and dispersion of stress field using mathematical equations
Stress data analysis was performed using the recent developed tensor-based mathematical formulas (Gao 2017;Gao (3) Harrison 2016Harrison , 2018. Euclidean distance was used for the distribution of local stress perturbation, which represents the distance between the local stress tensor, S i and the mean stress tensor, S̅ .
In two-dimensional models of the stress tensor field S, which is comprised of n part of the stress size, the stress tensor of part i is formulated as follows: The mean stress field is also determined using the following equation: Figure 5 shows that the mean stress tensor is equivalent to the far-field stress tensor (different ratios of tectonic stress loading). Figure 6 illustrates the geometry and meshing of the finite difference model used in the present study in a discrete fracture network sample (a = 1.5 and P 20 = 80 m −2 ) and the way of determination of the local stress tensor at two different positions. Table 2 presents the way of calculation of local stress tensor and d(S i , S̅ ) for the same fracture sample at the stress ratio SR = 3, (in the main parameters).
For the S i stress tensor, the form of the distinct components of the s i stress vector is as follows:

Results and data analysis
The analysis of results obtained regarding the effect of different rock and fracture parameters on 2D fractured rock's mechanical behavior under several different tectonic stress loading conditions is presented. In the next step, following the stochastic generation of the fracture network, which comprised different realizations that were differentiated considering the various distinct lengths and densities of fractures that were considered, the numerical modeling was performed. In this process finite-difference method was employed and three boundary loading conditions of 1) S xx = 10 MPa, S yy = 10 MPa, 2) S xx = 20 MPa, S yy = 10 MPa and 3) S xx = 30 MPa, S yy = 10 MPa were applied orthogonally.
Sequentially, the analysis of stress data was carried out using mathematical formulations to examine the effect of the fracture network's geometric properties and rock and T fracture parameters on the distribution of stress. Figures 7,8,9 and 10 show the relationship between the fracture density and the power-law length exponent a (fracture length indicator) with the mean local stress perturbation md(S, S̅ ) and the effective variance V e (stress distribution scattering indicator) at different values of k n , ψ F , C R and σ tR , respectively. It is conspicuous that two rock mass parameters, C R and σ tR , have minimal effect on the variations of md (S, S̅ ) and V e therefore, their influence can be ignored. The fracture normal stiffness, k n , can be defined as the normal spring stiffness of each joint element. It mainly acts as the link between the normal stress on each element to the normal displacement. The data explicitly indicate that the fracture normal stiffness has a significant influence on the stress perturbation distribution, md (S, S̅ ) and V e . In general, as k n increases, the values of md (S, S̅ ) and V e decrease (the difference in the decrease in V e values is more significant from k n = 500 GPa/m to k n = 1000 GPa/m than from the k n = 1000 GPa/m to k n = 2000 GPa/m). Also, the difference in V e values for different ψ F is higher at high fracture density (P 20 = 320 m −2 ). In this study, no change in md (S, S̅ ) and V e was observed in all fracture network realizations by varying the tensile strength (σ tF ) from zero to 5 MPa.
Moreover, from these figures and Table 3, it was found that the mean local stress perturbation and effective variance, in the order of highest to lowest, are sensitive to fracture density, fracture length and fracture normal stiffness. Figures 11, 12 and 13 show the distribution of local stress perturbation in fractured rocks at different k n values (as the most important parameter affecting the perturbation and distribution of the local stress) at tectonic stress ratios of 1, 2 and 3, respectively. The stress distribution at the stress ratio of 1 (SR = 1), which is isotropic, is quite uniform, while fluctuations are observed in anisotropic stresses (SR = 2 and SR = 3). With the increase in tectonic stress, local stress perturbations are detected at the tip and intersection of fractures. Generally, in the generated fracture networks, local Fig. 11 Distribution of local stress perturbation at stress ratio of 1 at different normal stiffness values for (a) fracture density of P 20 = 80 m −2 , (b) fracture density of P 20 = 160 m −2 and (c) fracture density of P 20 = 320 m −2 ◂ stress perturbations decreased with the increase in powerlaw length exponent (a) (decrease in fracture length) and normal stiffness (k n ) and it increased with an increase in fracture density (P 20 ).

Conclusion
The geometrical properties of fracture network and different rock and fracture parameters have a significant influence on the behavior of fractured rocks. In this study, these effects on stress distribution in a geological environment were investigated.
1. The fracture network with fracture density of 320 m −2 , power length exponent of 1.5, at 30-10 tectonic stress loading and normal stiffness of 500 GPa/m (as the parameter that most critically affects the stress distribution), has the highest local stress perturbation compared to other fracture networks presented in this study.
2. As the stress ratio S ∞ xx ∕S ∞ yy increases, the stress perturbation becomes more noticeable at the tip and intersection of fractures. 3. The mean distribution of local stress perturbation (md(S, S̅ )) and effective variance (V e ) have a direct relationship with the increase in stress ratio and fracture density and an inverse relationship with increasing power-law length exponent (fracture length reduction). 4. Investigation of the effect of four different parameters of rock and fracture revealed that k n has a significant effect on stress distribution. The rock parameters, including C R and σ tR , have minute effect on the variations of md (S, S̅ ) and V e . Generally, with increasing k n , the values of md(S, S̅ ) and V e are decreased. 5. To give an instance, in one of the studied cases (P 20 = 320 m −2 , a = 1.5 and SR = 3), with increasing k n value from 500 to 2000 GPa/m, md(S, S̅ ) and V e were decreased from 11.5 to 8.2 MPa (28.69%) and from 30 to 12.1 MPa 2 (59.66%), respectively. It can be stated that the difference in the decrease of V e for k n = 500 GPa/m to k n = 1000 GPa/m is higher than k n = 1000 GPa/m to k n = 2000 GPa/m and the difference in values of V e for ψ F = 0° and ψ F = 15° is observed to a certain extent at the highest fracture density studied (P 20 = 320 m −2 ).

Compliance with ethical standards
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://creat iveco mmons .org/licen ses/by/4.0/.

Fig. 13
Distribution of local stress perturbation at stress ratio of 3 at different normal stiffness values for (a) fracture density of P 20 = 80 m −2 , (b) fracture density of P 20 = 160 m −2 and (c) fracture density of P 20 = 320 m −2 ◂