Applied Science, Faculty of
Non UBC
Mining Engineering, Keevil Institute of
DSpace
Ojeda, P.; Elmo, D.; Rogers, S.; Brzovic, A. Discrete Fracture Network (DFN) Analysis to Quantify the Reliability of Borehole-Derived Volumetric Fracture Intensity. Geosciences 2023, 13, 187.
Ojeda, Pedro
Elmo, Davide
Rogers, Steve
Brzovic, Andres
Volumetric fracture intensity (P32) is a parameter that plays a major role in the mechanical and hydraulic behaviour of rock masses. While methods such as Ground Penetrating Radar (GPR) are available to map the 3D geometrical characteristics of the fractures, the direct measurement of P32 at a resolution compatible with geotechnical applications is not yet possible. As a result, P32 can be estimated from the borehole and surface data using either simulation or analytical solutions. In this paper, we use Discrete Fracture Network (DFN) models to address the problem of estimating P32 using information from boreholes (1D data). When calculating P32 based on Terzaghi Weighting, it is common practice to use drill run lengths and limit the minimum angle between the borehole and the intersected fractures. The analysis presented in this paper indicated that limiting the minimum angle of intersection would result in an underestimation of the calculated P32. Additionally, the size of the interval has a significant impact on the variability of the calculated P32. We propose a methodology to calculate P32 using variable lengths, depending on the angle between the fractures and the borehole. This methodology allows the capture of the spatial variation in intensity and simultaneously avoids artificially increasing or decreasing the intensity sampled along borehole intervals. Additionally, this work has addressed the impact of boundary effects in DFN models and proposes a methodology to mitigate them.
Citation: Ojeda, P.; Elmo, D.; Rogers,S.; Brzovic, A. Discrete FractureNetwork (DFN) Analysis to Quantifythe Reliability of Borehole-DerivedVolumetric Fracture Intensity.Geosciences 2023, 13, 187. https://doi.org/10.3390/geosciences13060187Academic Editors: Jesús Martínez-Frías and Hongyuan LiuReceived: 26 April 2023Revised: 13 June 2023Accepted: 15 June 2023Published: 18 June 2023Copyright: © 2023 by the authors.Licensee MDPI, Basel, Switzerland.This article is an open access articledistributed under the terms andconditions of the Creative CommonsAttribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).geosciencesArticleDiscrete Fracture Network (DFN) Analysis to Quantify theReliability of Borehole-Derived Volumetric Fracture IntensityPedro Ojeda 1, Davide Elmo 2,*, Steve Rogers 1 and Andres Brzovic 11 WSP, Vancouver, BC V6Z 2M1, Canada; pedro.ojeda@wsp.com (P.O.); steve.rogers2@wsp.com (S.R.);andres.brzovic@wsp.com (A.B.)2 NBK Institute of Mining Engineering, University of British Columbia, Vancouver, BC V6T 1Z4, Canada* Correspondence: davide.elmo@ubc.caAbstract: Volumetric fracture intensity (P32) is a parameter that plays a major role in the mechanicaland hydraulic behaviour of rock masses. While methods such as Ground Penetrating Radar (GPR)are available to map the 3D geometrical characteristics of the fractures, the direct measurement of P32at a resolution compatible with geotechnical applications is not yet possible. As a result, P32 can beestimated from the borehole and surface data using either simulation or analytical solutions. In thispaper, we use Discrete Fracture Network (DFN) models to address the problem of estimating P32using information from boreholes (1D data). When calculating P32 based on Terzaghi Weighting, it iscommon practice to use drill run lengths and limit the minimum angle between the borehole and theintersected fractures. The analysis presented in this paper indicated that limiting the minimum angleof intersection would result in an underestimation of the calculated P32. Additionally, the size of theinterval has a significant impact on the variability of the calculated P32. We propose a methodology tocalculate P32 using variable lengths, depending on the angle between the fractures and the borehole.This methodology allows the capture of the spatial variation in intensity and simultaneously avoidsartificially increasing or decreasing the intensity sampled along borehole intervals. Additionally, thiswork has addressed the impact of boundary effects in DFN models and proposes a methodology tomitigate them.Keywords: volumetric fracture intensity; discrete fracture network; borehole data; boundary effects1. IntroductionA key component of understanding the behaviour of a rock mass is the characterizationof discontinuities (i.e., structural features such as joints, bedding planes, faults, shear zones,veins). Discontinuities in a rock mass are usually described as an assemblage rather thanindividually [1], using specific properties, namely discontinuity orientation, size, shape,spatial location, and intensity.The necessity of modelling discontinuities for different engineering and geologicalapplications has led to the increasing use of the Discrete Fracture Network (DFN) approach,both as a stand-alone tool or integrated within more complex geomechanical models, asdiscussed in [2]. Note that the term fracture is herein applied in a broad context to referto a variety of discontinuities at different scales (faults, joints, veins, etc.). DFN modellinghas become increasingly popular in recent years amongst geotechnical practitioners andengineers, including the generation of synthetic rock mass properties [3,4], geomechanicalsimulation of open pits [5,6], quantification of rock mass pre-conditioning [7], estimation ofrock bridge percentage for stability analysis [8] and rock mass fragmentation and calculationof fragmentation distribution at cave mine scale [9].One of the most important inputs for DFN analysis is fracture intensity, which can beexpressed in various ways, including (i) the total sum of fracture surface area per unit ofvolume, P32; (ii) the total sum of fracture length per unit area, P21; and (iii) the number offractures per unit of length, P10.Geosciences 2023, 13, 187. https://doi.org/10.3390/geosciences13060187 https://www.mdpi.com/journal/geosciencesGeosciences 2023, 13, 187 2 of 22P32 is a nondirectional intrinsic measure of fracture intensity, and these characteristicsmake P32 the preferred measure of fracture intensity in DFN modelling [9]. While P32cannot be directly measured, it can be obtained from P10 or P21 using analytical solutionsor simulations [10].The simulation method has a number of key limitations, namely:• It is necessary to know beforehand the fracture size distribution and orientation togenerate the simulated DFN models.• The proportionality between 1D/2D and 3D fracture intensity depends on the relativeorientation of the fractures to the orientation of the survey line or planes.Zhang and Einstein [11] proposed a methodology to estimate the volumetric intensityof discontinuities by combining the data sampled on boreholes with the data sampled onexposed rock faces. The limitations of this methodology are that it requires assumptionsabout the fractures’ location and size distribution and the probability that a discontinuity,with its centroid in an objective volume, will intersect a borehole. Most importantly, itrequires exposed rock faces, which are not always available, particularly in the initial stagesof a project. In addition, even if surface data are available, this may not represent the under-ground conditions. Wang [12] proposed a numerical approximation between P10 and P32;the limitation of this methodology is that it assumes that the fracture population in a singleset follows a Univariate Fisher hemispherical probability distribution, with the methodnot being suitable for use with other types of probability distributions. To overcome thoselimitations Chilès et al. [13] developed a methodology based on Terzaghi Weighting [14] toestimate P32 directly from borehole intensity. Chilès’ method [13] has been widely adoptedin the DFN community. This paper specifically addresses the question remaining of usingself-assumed practical industry standards to excuse the lack of geological information (e.g.,fracture size).2. Chilès’ Methodology to Calculate P32 from Borehole IntensityWhen estimating fracture intensity from a linear or planar survey, fractures that areperpendicular to the survey are easily observed, while oblique fractures are harder toobserve. The relative difference between the orientation of the sampling borehole and theorientation of the intersecting fractures introduces a bias in favour of those fractures thatare perpendicular to the borehole. This orientation bias is illustrated in Figure 1, in whichα represents the acute angle between the scanline and the fracture. The apparent spacing(D’) ratio to the true spacing (D) is 1/sinα. In other words, fractures at an acute angle α areunderrepresented by a factor of “sinα”; this factor is known as the Terzaghi Weighting [14].Geosciences 2023, 13, x FOR PEER REVIEW 2 of 23 of volume, P32; (ii) the total sum of fracture length per unit area, P21; and (iii) the number of fractures per unit of length, P10. P32 is a nondirectional intrinsic measure of fracture intensity, and these characteristics make P32 the preferred measure of fracture intensity in DFN modelling [9]. While P32 can-not be directly measured, it can be obtained from P10 or P21 using analytical solutions or simulations [10]. The simulation method has a number of key limitations, namely: • It is necessary to know beforehand the fracture size distribution and orientation to generate the simulated DFN models. • The proportionality between 1D/2D and 3D fracture intensity depends on the relative orientation of the fractures to the orientation of the survey line or planes. Zha g and Einstein [11] propos d a meth d logy to estimat the volumetric inten-sity of discontinuities by combining the data sampled on boreholes with the data sampled on exposed rock faces. The limitations of this methodology are that it requires assump-tions about the fractures’ location and size distribution and the probability that a discon-tinuity, with its centroid in an objective volume, will intersect a borehole. Most im-portantly, it requires exposed rock faces, hich are not always available, particularl in the initial stag s of a project. In addition, even if surf ce data are vailable, this may not represent the un erground conditions. W ng [12] proposed a numerical approximation between P10 and P32; the limitation of this methodology is that it assumes that the fracture population i a single set follows a Univariate Fisher hemispherical prob bility distribu-tion, with the method not being suitable for use with other types of probability distribu-ti ns. To overcome tho limitations Chilès et al. [13] developed a method logy based on Terz ghi Weighting [14] to estimate P32 directly from boreh le intensity. Chilès’ method [13] has been widely adopted in the DFN community. This paper sp cifically ad ress s the question remai ing of using self-assumed practical industry sta dards to excuse the lack of geologic l information (e.g., fracture size). 2. Chilès’ Methodology to Calculate P32 from Borehole Intensity When estimating fracture intensity from a linear or planar survey, fractures that are perpendicular to the survey are easily observed, while oblique fractures are harder to ob-serve. The relative difference between the orientation of the sampling borehole and the orientation of the intersecting fractures introduces a bias in favour of those fractures that are perpendicular to the borehole. This orientation bias is illustrated in Figure 1, in which α represents the acute angle between the scanline and the fracture. The apparent spacing (D’) ratio to the true spacing (D) is 1/sinα. In other words, fractures at an acute angle α are un-derrepresented by a factor of “sinα”; this factor is known as the Terzaghi Weighting [14]. Figure 1. Apparent spacing D’ associated with true spacing D when the fractures (direction θ) form an angle α with the scanline (direction β). Adapted from Chilès et al. [13]. Terzaghi [14] proposed a correction for spacing (Equation (1)), in which d represents the spacing, N is the number of fractures intersected by the survey, and L is the length of the survey. Note that the Terzaghi Correction assumes a zero-thickness sampling survey. βi . t i ’ i t it t i t f t ( i ti ) fan angle α ith the scanline (direction β). Adapted fro Chilès et al. [13].ti r -t ic ess sa ling survey.d =NL sin α(1)Geosciences 2023, 13, 187 3 of 22According to Chilès et al. [13], Equation (1) can be rewritten in terms of the predictedvolumetric intensity P32.Pˆ32 =1L∑n11|cos(ωi − β)| (2)Chilès et al. [13] presented a generalization of Equation (2) for the case in whichfractures are sampled along a line with varying orientation (βi) or several lines with totallength L:Pˆ32 =1L∑n11|cos(ωi − βi)| (3)where L is the length of the scanline, βi is the local orientation of the scanline at the locationof fracture i, N is the number of fractures, and ωi the pole of fracture i. Note that in 3D,βi andωi are unit vectors and the factor equal to the cosine of the angle formed by the unitvectors, βi and ωi corresponds to the absolute value of the inner product <ωi, βi>.Equation (3) can be written in terms of the acute angle between the scanline and thefracture i (αi) asPˆ32 =1L∑n11sin(αi)(4)Since the weighting factor tends to infinity as α approaches zero, applying a maxi-mum limit to α (minimum bias angle) is expected. Priest [15] suggests using a maximumweighting of 10; this weight corresponds to an α value of approximately 6◦, while otherauthors recommend an α angle between 5◦ and 20◦ [16]. In practice, a limit of 15◦ is oftenused to improve the robustness of the estimator. However, this has the inconvenience ofintroducing some bias by discarding valid data [13].3. DFN Models and Boundary EffectsThe engineering community often disregards boundary effects when generating DFNmodels, even though their impact can be significant when the volume used to create theDFN model is not large enough compared to the rock mass volume under consideration.Priest [15] recommended generating fractures in a volume region much larger than thevolume of interest to minimize these boundary effects. Samaniego and Priest [17] indicatedthat a generation volume of four times the volume of interest would be sufficient to reduceboundary effects. Note that depending on the model size, this increase in the size ofthe generation volume may significantly increase the computation time. On the otherhand, Junkin et al. [18] indicate that it could be possible to mitigate boundary effects byincluding mapped traces or fracture seeds on the exterior boundaries of a DFN region.They proposed to sample the DFN model using planes to estimate P21 intensity in areaswhere the boundary effect does not exist and then use this intensity to stochastically seedtraces on the boundary of the DFN.This Section is dedicated to studying the boundary effects that may influence the inten-sity of DFN models. DFN models were generated using FracMan Software Version 8 [19].A new methodology is proposed to minimize and correct the effect of boundary effects onDFN models’ intensity.3.1. DFN Modelling to Investigate Boundary Effects on IntensityTo determine and quantify the impact of boundary effects on volumetric fractureintensity (P32), a set of DFN models was generated within a box of 150 m per side.P32 values were then calculated using progressively smaller boxes (Figure 2). The inputP32 values were compared with the P32 obtained for the different volumes. This modellingexercise aimed to quantify the generation box’s effect on the models’ target intensity.Geosciences 2023, 13, 187 4 of 22Geosciences 2023, 13, x FOR PEER REVIEW 4 of 23 Figure 2. Generation box (outer box) and sampling boxes (inner boxes) are used to calculate P32. Only some of the sampling boxes are shown for illustrative purposes. Fracture locations were generated using both centre and surface points. For centres, the generation point is the centre of the fracture, while in surface points, the generation point of the fracture is a random point within the fracture. According to the FracMan Manual [19], this option can reduce boundary effects when generating P32 values. Table 1 summarizes the properties used in the models. Table 1. Input parameter used in the models to Investigate boundary effects. Property Note Orientation Bootstrapped using a concentration parameter of 80 (Figure 3) Spatial model Enhanced Baecher with generation locations at fracture centres and surface point Intensity P32 (m−1) 2 Size: Lognormal Distribution Radius Mean (Xmean): 2-5-10-15-20-25 (m) and standard deviation (SD) of 40 % of the mean radius size Fracture Shape Hexagon with a constant aspect ratio of 1 Number of realizations 100 per scenario (1200 realizations) Figure 2. Generation box (outer box) and sampling boxes (inner boxes) are used to calculate P32.Only some of the sampling boxes are shown for illustrative purposes.Fracture locations were generated using both centre and surface points. For centres,the generation point is the centre of the fracture, while in surface points, the generationpoint of the fracture is a random point within the fracture. According to the FracManManual [19], this option can reduce boundary effects when generating P32 values. Table 1summarizes the properties used in the models.Table 1. Input parameter used in the models to Investigate boundary effects.Property NoteOrien ation Bootstrapped usi g a concentration parameter of 80 (Figure 3)Spatial model Enhanced Baecher with generation locations at fracturecentres and surface pointIntensity P32 (m−1) 2Size: Lognormal Distribution Radius Mean (Xmean): 2-5-10-15-20-25 (m) and standarddeviation (SD) of 40% of the mean radius sizeFracture Shape Hexagon with a constant aspect ratio of 1Number of realizations 100 per scenario (1200 realizations)Geosciences 2023, 13, x FOR PEER REVIEW 5 of 23 Figure 3. Stereonet with orientation used for bootstrapping, 563 Poles including Terzaghi Weighting. 3.2. Results Figure 4 summarizes the results for a mean fracture radius size of 20 m and a P32 input of 2 m−1. The results show a difference of up to 30% between the input and target intensity, thus demonstrating that boundary effects can affect the intensity. P32 values sam-pled at a scale corresponding to the generation box yield the input P32. Still, as the volume of the sampling region is progressively reduced (assuming the coordinates of the centre of every generation volume remain the same), the average P32 value increases, showing an asymptotic behaviour. Contrary to what may be expected, the models with fractures gen-erated using their centres as reference present a better agreement with the input P32, and the average curve stabilizes faster than the one for the models in which fractures are gen-erated using the surface points option. Figure 4. Boundary effect and effect of the sampling box in sampled P32. (a) Models generated using fractures by centres approach. (b) Models generated using fractures generated by the surface points approach. NSEWSymbol FeaturePole VectorsColor Density Concentrations0.00 - 0.550.55 - 1.101.10 - 1.651.65 - 2.202.20 - 2.752.75 - 3.303.30 - 3.853.85 - 4.404.40 - 4.954.95 - 5.50Contour Data Pole VectorsMaximum Density 5.07%Contour Distribution FisherCounting Circle Size 1.0%Plot Mode Pole VectorsVector Count 1102 (563 Entries)Hemisphere LowerProjection Equal Angle-20%-10%0%10%20%30%40%50%60%1.61.82.02.22.42.62.83.03.220 40 60 80 100 120 140 160VariationSampled P32(m-1)Sampling Box Side Length (m)a) Fractures Generated in Centre - Size LogNorm (Mean:20m, SD: 8m)Min. Max. Average-20%-10%0%10%20%30%40%50%60%1.61.82.02.22.42.62.83.03.220 40 60 80 100 120 140 160VariationSampled P32(m-1)Sampling Box Side Length (m)b) Fractures Generated in Surface Points - Size LogNorm (Mean:20m, SD: 8m)Min. Max. AverageFigure 3. Stereonet with orientation used for bootstrapping, 563 Poles including Terzaghi Weighting.Geosciences 2023, 13, 187 5 of 223.2. ResultsFigure 4 summarizes the results for a mean fracture radius size of 20 m and a P32 inputof 2 m−1. The results show a difference of up to 30% between the input and target intensity,thus demonstrating that boundary effects can affect the intensity. P32 values sampledat a scale corresponding to the generation box yield the input P32. Still, as the volumeof the sampling region is progressively reduced (assuming the coordinates of the centreof every generation volume remain the same), the average P32 value increases, showingan asymptotic behaviour. Contrary to what may be expected, the models with fracturesgenerated using their centres as reference present a better agreement with the input P32,and the average curve stabilizes faster than the one for the models in which fractures aregenerated using the surface points option.Geosciences 2023, 13, x FOR PEER REVIEW 5 of 23 Figure 3. Stereonet with orientation used for bootstrapping, 563 Poles including Terzaghi Weighting. 3.2. Results Figure 4 summarizes the results for a mean fracture radius size of 20 m and a P32 input of 2 m−1. The results show a difference of up to 30% between the input and target intensity, thus demonstrating that boundary effects can affect the intensity. P32 values sam-pled at a scale corresponding to the generation box yield the input P32. Still, as the volume of the sampling region is progressively reduced (assuming the coordinates of the centre of every generation volume remain the same), the average P32 value increases, showing an asymptotic behavi ur. Contrary to w at may be exp cted, the models with fractures gen-erated using their centres as reference present a etter agreem nt with the input P32, and the average curve stabilizes faster than the one for the models in which fractures are gen-erated using the surface points option. Figure 4. Boundary effect and effect of the sampling box in sampled P32. (a) Models generated using fractures by centres approach. (b) Models generated using fractures generated by the surface points approach. NSEWSymbol FeaturePole VectorsColor Density Concentrations0.00 - 0.550.55 - 1.101.10 - 1.651.65 - 2.202.20 - 2.752.75 - 3.303.30 - 3.853.85 - 4.404.40 - 4.954.95 - 5.50Contour Data Pole VectorsMaximum Density 5.07%Contour Distribution FisherCounting Circle Size 1.0%Plot Mode Pole VectorsVector Count 1102 (563 Entries)Hemisphere LowerProjection Equal Angle-20%-10%0%10%20%30%40%50%60%1.61.82.02.22.42.62.83.03.220 40 60 80 100 120 140 160VariationSampled P32(m-1)Sampling Box Side Length (m)a) Fractures Generated in Centre - Size LogNorm (Mean:20m, SD: 8m)Min. Max. Average-20%-10%0%10%20%30%40%50%60%1.61.82.02.22.42.62.83.03.220 40 60 80 100 120 140 160VariationSampled P32(m-1)Sampling Box Side Length (m)b) Fractures Generated in Surface Points - Size LogNorm (Mean:20m, SD: 8m)Min. Max. AverageFigure 4. Boundary effect and effect of the sampling box in sampled P32. (a) Models generatedusing fractures by centres approach. (b) Models generated using fractures generated by the surfacepoints approach.The variability increases as the sampling box size decreases, meaning that the ratiobetween sampling box size and fracture size influences the variability. For the case in whichfractures are generated at the centre point, there is a marked change in the slope of thecurves for the maximum and minimum values when the edge length of the generationvolume is twice the fracture radius.The difference between input P32 and sampled P32 can be explained if one considersthe underlying Baecher spatial model [20], which is used by default by many engineersand practitioners and it assumes fractures to be located uniformly in space, meaning thatthe fractures with generation points located outside the generation box and at a distancefrom the boundary smaller than a fracture diameter are not generated.In a hypothetical infinite generation space, those fractures would extend back intothe original boundary box, but in a constrained model, this results in a lower intensitynear the boundary. This effect was already described by Priest [15], who recommendedusing a smaller volume of interest located at the centre of the generation volume. Notethat the size of the generation box will depend on the fracture size distribution. Samaniegoand Priest [17] indicated that using a sampling area located at the central quarter of thegeneration volume minimizes the boundary effects. Depending on the model size, thismay significantly increase the computation time. Furthermore, since the boundary effect isrelated to the assumed fracture size, this may somewhat reduce the boundary effect.Geosciences 2023, 13, 187 6 of 223.2.1. Minimizing Boundary EffectsThe asymptotic behaviour observed (Figure 4) indicates that it is possible to define acorrection factor for the input P32. This correction factor can be defined as the ratio betweenthe input intensity and the average P32 value, represented by the asymptotic section of thecurve (Figure 5a). This factor will depend on the size of the volume of interest, the size ofthe generation box, the fracture size distribution and the fracture orientation.Geosciences 2023, 13, x FOR PEER REVIEW 6 of 23 The variability increases as the sampling box size decreases, meaning that the ratio between sampling box size and fracture size influences the variability. For the case in which fractures are generated at the centre point, there is a marked change in the slope of the curves for the maximum and minimum values when the edge length of the generation volume is twice the fracture radius. The difference between input P32 and sampled P32 can be explained if one considers the underlying Baecher spatial model [20], which is used by default by many engineers and practitioners and it assumes fractures to be located uniformly in space, meaning that the fractures with generation points located outside the generation box and at a distance from the boundary smaller than a fracture diameter are not generated. In a hypothetical infinite generation space, those fractures would extend back into the original boundary box, but in a constrained model, this results in a lower intensity near the boundary. This effect was already described by Priest [15], who recommended using a smaller volume of interest located at the centre of the generation volume. Note that the size of the generation box will depend on the fracture size distribution. Samaniego and Priest [17] indicated that using a sampling area located at the central quarter of the generation volume minimizes the boundary effects. Depending on the model size, this may significantly increase the computation time. Furthermore, since the boundary effect is related to the assumed fracture size, this may somewhat reduce the boundary effect. 3.2.1. Minimizing Boundary Effects The asymptotic behaviour observed (Figure 4) indicates that it is possible to define a correction factor for the input P32. This correction factor can be defined as the ratio between the input intensity and the average P32 value, represented by the asymptotic section of the curve (Figure 5a). This factor will depend on the size of the volume of interest, the size of the generation box, the fracture size distribution and the fracture orientation. Figure 5. Effect of correction factor in target intensity. (a) Correction factor calculation for a volume of interest (VOI) of 100 m per side, using an initial input P32 of 2 m−1. (b) Target intensity P32 of 2 m−1 obtained for a volume of interest of 100 m per side, using a correction factor of 0.84 applied to the initial input intensity. For example, let us consider the model in Figure 4a, with a mean fracture size of 20 m and fractures generated at the centre point. For this model, initial DFNs were developed considering an initial input intensity P32 of 2 m−1. Then average P32 values were calculated using progressively smaller boxes. To define a target intensity P32 of 2 m−1 in a volume of interest of 100 m per side, one would need to adjust the input intensity of the model by a correction factor, which for the specific models used in the analysis was calculated as 84% = = Figure 5. Effect of correction factor in target intensity. (a) Correction factor calculation for a volumeof interest (VOI) of 100 m per side, using an initial input P32 of 2 m−1. (b) Target intensity P32 of2 m−1 obtained for a volume of interest of 100 m per side, using a correction factor of 0.84 applied tothe initial input intensity.For example, let us consider the model in Figure 4a, with a mean fracture size of 20 mand fractures generated at the centre point. For this model, initial DFNs were developedconsidering an initial input intensity P32 of 2 m−1. Then average P32 values were calculatedusing progressively smaller boxes. To define a target intensity P32 of 2 m−1 in a volume ofinterest of 100 m per side, one would need to adjust the input intensity of the model by acorrection factor, which for the specific models used in the analysis was calculated as 84%of he target intensity within th volu e of inter st. This reduction of 84% correspondsto the factor in which the average curve of P32 stabilizes (2.00/2.38 = 0.84 of the averagesampled P32 presented in Figure 5a). When using this methodology, obtaining the targetP32 of 2 m−1 on the volume of interest is possible. The curve stabilizes approximately 50 mfrom the border of the generation box (Figure 5b). This distance is slightly lower than threetimes the mean fracture radius (60 m). Based on these results, models were generated fora range of fracture sizes (Table 1), using a generation box size of 100 m plus three timesthe mean fracture radius, assuming that the volume of interest corresponds to a box of100 m per side. The results of this exercise are presented in Figure 6 and show that for theassumed fracture orientation and size, it is possible to minimize the boundary effects whenusing a generation box size of 100 m plus three times the mean fracture size.Geosciences 2023, 13, 187 7 of 22Geosciences 2023, 13, x FOR PEER REVIEW 7 of 23 of the target intensity within the volume of interest. This reduction of 84% corresponds to the factor in which the average curve of P32 stabilizes (2.00/2.38 = 0.84 of the average sam-pled P32 presented in Figure 5a). When using this methodology, obtaining the target P32 of 2 m−1 on the volume of interest is possible. The curve stabilizes approximately 50 m from the border of the generation box (Figure 5b). This distance is slightly lower than three times the mean fracture radius (60 m). Based on these results, models were generated for a range of fracture sizes (Table 1), using a generation box size of 100 m plus three times the mean fracture radius, assuming that the volume of interest corresponds to a box of 100 m per side. The results of this exercise are presented in Figure 6 and show that for the assumed fracture orientation and size, it is possible to minimize the boundary effects when using a generation box size of 100 m plus three times the mean fracture size. Figure 6. Effect of box size on intensity, using generation box sizes depending on fracture mean radius. 3.2.2. Required Number of Realizations To investigate the number of realizations required to obtain a meaningful sampled average P32 in the volume of interest, the cumulative average of the sampled P32 was plot-ted against the number of realizations (Figure 7). This provides a graphical representation that determines whether the number of realizations converges to a constant sampled P32. The number of realizations required will depend on the size of the fractures. When fractures are small, the plot converges quickly, meaning that just a small number of reali-zations are needed to obtain a meaningful average P32. Conversely, more realizations are required for the results to converge when the fracture size increases. The analysis shows that for a mean fracture radius of 25 m, approximately 30 to 50 realizations are required. In comparison, for a mean fracture radius of 2 m, as few as 10 realizations may be sufficient to obtain a representative mean P32. Note that the number of realizations should increase if the plots do not converge. 1.81.922.12.22.32.42.520 40 60 80 100 120 140 160 180Sampled P32(m-1)Sampling Box Side Length (m)Generation Box Size Depending on Fracture Mean RadiusXmean: 2 m Xmean: 5 m Xmean: 10 m Xmean: 15 m Xmean: 20 m Xmean: 25 mFigure 6. Effect of box size on intensity, using generation box sizes depending on fracture mean radius.3.2.2. Required Number of RealizationsTo investigate the number of realizations required to obtain a meaningful sampledaverage P32 in the volume of interest, the cumulative average of the sampled P32 was plottedagainst the number of realizations (Figure 7). This provides a graphical representation thatdetermines whether the number of realizations converges to a constant sampled P32.Geosciences 2023, 13, x FOR PEER REVIEW 8 of 23 Figure 7. Cumulative average of sampled P32 for a box of 100 m per side (volume of interest). 3.2.3. Effect of Input Intensity on Target Intensity The effect of the input intensity was investigated by generating models with different intensities and comparing the variation in percentage between the input P32 and the P32 sampled in region volumes of different sizes. For this exercise, fractures were generated in a generation volume with a side length of 125 m, using the properties presented in Table 1. Table 2 shows the results of this exercise for a mean fracture radius of 25 m. Table 2. Percentage variation between input P32 and sampled P32 for fractures with a mean radius of 25 m. Side Length (m) Note 1 (m−1) 2 (m−1) 3 (m−1) 4 (m−1) 5 (m−1) 25 27.6% 24.2% 23.9% 23.9% 23.7% 50 25.0% 24.0% 23.9% 24.0% 24.1% 75 23.8% 23.4% 23.3% 23.3% 23.4% 100 18.4% 18.3% 18.3% 18.3% 18.4% 125 0.0% 0.0% 0.0% 0.0% 0.0% It can be observed in Table 2 that for a defined box size, the variation in percentage, with respect to the input intensity, is relatively constant, suggesting that if the other prop-erties are maintained constant, the correction factor is independent of the fracture inten-sity; therefore, the same correction factor can be applied to models with different intensity values. This will be useful in Section 4, where the linear and volumetric intensity relation-ship will be investigated for various intensities. 3.3. Methodology Proposed to Mitigate Boundary Effects on Intensity A methodology is proposed to mitigate the impact of boundary effects on intensity when using DFN models generated in regions with different sizes. The process is de-scribed in the step-by-step guide below and illustrated in Figure 8. 1. Define a volume of interest, which corresponds to the volume in which we want to define the target P32. 2. Generate fractures in a region bigger than the volume of interest. The size of the gen-eration volume depends on the mean fracture size and its standard deviation. A sen-sitivity analysis could be performed to estimate the size of the generation volume. 2.002.052.102.152.202.252.302.352.402.452.500 10 20 30 40 50 60 70 80 90 100Cumulative Average of Sampled P32(m-1)RealizationSampl d P32 Converge cy - Box 100 m per side Xmean: 2 m Xmean: 5 m Xmean: 10 m Xmean: 15 m Xmean: 20 m Xmean: 25 mFigure 7. Cumulative average of sampled P32 for a box of 100 m per side (volume of interest).The number of realizations required will depend on the size of the fractures. Whenfractures are small, the plot converges quickly, meaning that just a small number of real-izations are needed to obtain a meaningful average P32. Conversely, more realizations arerequired for the results to converge when the fracture size increases. The analysis showsthat for a mean fracture radius of 25 m, approximately 30 to 50 realizations are required. Incomparison, for a mean fracture radius of 2 m, as few as 10 realizations may be sufficient toobtain a represent tive mean P32. Note that t e number of e lizations should increase ifthe plots do not converge.Geosciences 2023, 13, 187 8 of 223.2.3. Effect of Input Intensity on Target IntensityThe effect of the input intensity was investigated by generating models with differentintensities and comparing the variation in percentage between the input P32 and the P32sampled in region volumes of different sizes. For this exercise, fractures were generated ina generation volume with a side length of 125 m, using the properties presented in Table 1.Table 2 shows the results of this exercise for a mean fracture radius of 25 m.Table 2. Percentage variation between input P32 and sampled P32 for fractures with a mean radiusof 25 m.Side Length(m)Note1 (m−1) 2 (m−1) 3 (m−1) 4 (m−1) 5 (m−1)25 27.6% 24.2% 23.9% 23.9% 23.7%50 25.0% 24.0% 23.9% 24.0% 24.1%75 23.8% 23.4% 23.3% 23.3% 23.4%100 18.4% 18.3% 18.3% 18.3% 18.4%125 0.0% 0.0% 0.0% 0.0% 0.0%It can be observed in Table 2 that for a defined box size, the variation in percent-age, with respect to the input intensity, is relatively constant, suggesting that if the otherproperties are maintained constant, the correction factor is independent of the fractureintensity; therefore, the same correction factor can be applied to models with differentintensity values. This will be useful in Section 4, where the linear and volumetric intensityrelationship will be investigated for various intensities.3.3. Methodology Proposed to Mitigate Boundary Effects on IntensityA methodology is proposed to mitigate the impact of boundary effects on intensitywhen using DFN models generated in regions with different sizes. The process is describedin the step-by-step guide below and illustrated in Figure 8.1. Define a volume of interest, which corresponds to the volume in which we want todefine the target P32.2. Generate fractures in a region bigger than the volume of interest. The size of thegeneration volume depends on the mean fracture size and its standard deviation. Asensitivity analysis could be performed to estimate the size of the generation volume.3. Run an adequate number of realizations with an initial P32 close to the target P32.The number of realizations to perform will depend on the fracture size distribution.Cumulative average plots can be used to estimate the appropriate number of real-izations. If the plots do not converge, then the number of realizations should beincreased. Conversely, if the curves converge quickly to the average sampled P32,then it is possible to decrease the number of realizations. This will reduce the totalcomputation time, obtaining similar results.4. Calculate P32 in the volume of interest and other smaller volumes to check that thecurve of average sampled P32 has stabilized for the volume of interest. If the curvehas not stabilized, it is necessary to increase the size of the generation region (go backto step 2).5. If for the volume of interest, the curve is in the asymptotic section, calculate theaverage P32 in the volume of interest and the ratio between input P32 and the averageP32 in the volume of interest. This radio corresponds to the correction factor.6. Rerun the model using the initial P32 multiplied by the correction factor.7. Check that the P32 in the volume of interest corresponds to the target P32 within anacceptable margin of error. If the error exceeds the acceptable error, go back to step 6.Generating several DFN models and keeping only the ones that comply with theacceptable margin of error is recommended. For this paper, a difference of ±1% of thetarget intensity has been adopted.Geosciences 2023, 13, 187 9 of 22Geosciences 2023, 13, x FOR PEER REVIEW 10 of 23 Figure 8. Flowchart with the methodology proposed to minimize boundary effect on intensity. 4. Variability of Linear Fracture Intensity and Its Relationship with Volumetric Intensity The relationship between the volumetric intensity (P32) and the linear fracture inten-sity (P10) was investigated using a series of DFN models of known intensity. The DFN models were then sampled well with synthetics to obtain P10. In total, 48 DFN models were generated using as a volume of interest of 100 m × 100 m × 100 m box region (centred in the origin) and the properties presented in Table 3. To minimize boundary effects, fractures were generated using the methodology proposed in Section 3. Fractures were generated within a region volume with a variable side length equal to 100 m plus three times the mean fracture radius. Once the DFN models were generated, the intensity P32 was sampled in the volume of interest and only models with a difference smaller than 1% of the target intensity were considered. 1. Define a volume of interest (VOI)2. Generate fractures in a region bigger than the VOI3. Run adequate number of realizations and generate cumulative average plots of sampled P32 in the VOI Does cumulative average plots converge?No: Increase the number of realizationsYes: Number of realizations is adequate4. Calculate P32 in the VOI and other smaller volumes to check that the average sampled P32 curve has stabilized in the VOI. Has the average sampled P32 curve stabilized in the VOI?No: Increase the size of the generation regionYes: The size of the generation region is adequate5. Calculate the correction factor as the ratio between input P32 and the average P32sampled in the VOI 6. Rerun the model using the initial P32 multiplied by the correction factor 7. Check that the P32 in the VOI corresponds to the target P32 within an acceptable margin of error.Error of sampled P32 in the VOI is greater than acceptable margin of error: Model is considered inadequateError of sampled P32 in the VOI is within acceptable margin of error: Model is considered adequateFigure 8. Flowchart with the methodology proposed to minimize boundary effect on intensity.. riability of Linear Fracture Intensity and Its Relationship with Volumetric IntensityThe relationship between the volumetric intensity (P32) and the linear fracture intensity(P10) was investigated using a series of DFN models of known intensity. The DFN modelswere then sampled well with synthetics to obtain P10.In total, 48 DFN odels were generated using as a volume of interest of 100 m× 100 m× 100 m box region (centred in the origin) and the properties presented in Table 3. Toinimize boundary effects, fractures were generated using the methodology proposed inSection 3. Fractures were generated within a region volume with a variable side lengthequal to 100 m plus three times the mean fracture radius. Once the DFN models weregenerated, the intensity P32 was sampled in the volume of interest and only models with adifference smaller than 1% of the target intensity were considered.Geosciences 2023, 13, 187 10 of 22Table 3. Input parameter used in the models to investigate the reliability of borehole-derived intensity.Property CommentsOrientation Bootstrapped using a concentration parameter of 80 and theorientation presented in Figure 3Spatial model Enhanced Baecher with generation locations at fracture centreIntensity P32 (m−1) 1-2-3-4-5-6-7-8Size: Lognormal Distribution Radius Mean (Xmean): 2-5-10-15-20-25 (m) and standarddeviation (SD) of 40% of the mean radius sizeFracture Shape Hexagon with a constant aspect ratio of 1The DFN models were sampled using 10 synthetic wells (Figure 9), and the orientationsof the wells (boreholes) were defined, trying to cover the whole stereonet area. Figure 9shows the fracture orientations used as input and blue lines representing the blind zones forα = 0◦ of each well. Terzaghi [14] introduced the term blind zone of a drill hole to describethe locus of the poles of joints parallel to the drill hole and less likely to be observed. Notethat for illustrative purposes, the blind zones of each well are presented only as a line(α = 0◦), but all poles near that line (normally ± 20◦ from the line defined by α = 0◦) areless likely to be observed in that particular well.Geosciences 2023, 13, x FOR PEER REVIEW 11 of 23 Table 3. Input parameter used in the models to investigate the reliability of borehole-derived inten-sity. Property Comments Orientation Bootstrapped using a concentration parameter of 80 and the orientation presented in Figure 3 Spatial model Enhanced Baecher with generation locations at fracture centre Intensity P32 (m−1) 1-2-3-4-5-6-7-8 Size: Lognormal Distribution Radius Mean (Xmea ): 2-5-10-15-20-25 (m) and standard deviation (SD) of 40 % of the mean radius size Fracture Shape Hexagon with a const t aspect ratio of 1 The DFN models were sampled using 10 synthetic wells (Figure 9), and the orienta-tions of the wells (boreholes) were defined, trying to cover the whole stereonet area. Figure 9 shows the fracture orientations used as input and blue lines representing the blind zones for α = 0° of each well. Terzaghi [14] introduced the term blind zone of a drill hole to describe th locus of the poles of joints parall l to the drill hole and l ss likely to be ob-served. Note that f r illustrative purp ses, the blind zones of each well a e presented only as a line (α = 0°), but all poles near that line (normally ± 20° from the line defined by α = 0°) are less likely to be observed in that particular well. Figure 9. (Left) Location of synthetic wells, oblique view from the southwest. (Right) Stereonets with orientation used for bootstrapping. Blind zones for α = 0° of each well are presented as blue lines while the ±20° from the line defined by α = 0°are presented as black cones. When sampling the DFN models with the synthetic wells, a linear relationship be-tween input P32 and sampled P10 was observed (Figure 10). As expected, greater sampled P10 dispersion is observed for greater values of input P32. The linear relationship between P32 and P10 was already demonstrated using DFN modelling by Dershowitz et al. [10]. Figure 9. (Left) Location of synthetic wells, oblique view from the southwest. (Right) Stereonetswith orientation used for bootstrapping. Blind zones for α = 0◦ of each well are presented as bluelines while the ±20◦ from the line defined by α = 0◦are presented as black cones.When sampling the DFN models with the synthetic wells, a linear relationship betweeninput P32 and sampled P10 was observed (Figure 10). As xpected, greater s mpled P10dispersion is observed for greater values of input P32. The linear relationship between P32and P10 was already demonstrated using DFN modelling by Dershowitz et al. [10].Geosciences 2023, 13, 187 11 of 22Geosciences 2023, 13, x FOR PEER REVIEW 12 of 23 Figure 10. Sampled P10 in wells for a given input P32. Cases for a Mean Radius of 2 m and 25 m. When the models are generated using the methodology proposed to mitigate the boundary effects, the relationship between P32 and P10 is independent of the fracture size used (Figure 11). While this statement suggests that P32 can be estimated directly from borehole intensity, the generation of DFN models requires fracture size as an input and the output of DFN models must be validated against rock exposures. Figure 11. (Left) Sampled P10 probability distribution for an input P32 of 2 m−1, black dashed lines correspond to normal distributions using the mean and SD of each data set. (Right) Mean (μ) and standard deviation (σ) of sampled P10 follow a linear relationship. The results show that the mean and standard deviation of the sampled P10 follows a linear relationship. Therefore, estimating the P10 variability for any given P32 using a nor-mal probability function for a given fracture orientation is possible. Since P10 and P32 fol-low a linear relationship, using probability functions and interpolation, it is then possible to estimate the variability of P32 for any given P10 (Figure 12). . l 10 in ells for a given input P32. Cases for a Mean Radius of 2 m and 25 m.t l t i t t l t iti t tr ff ts, t r l ti s i t 3 is i t f t fr t r sise ( i re 11). ile t is state e t s ests t at ca e esti ate irectl frborehole intensity, the generation of DFN models requires fracture size as an input and theoutput of DFN models must be validated against rock exposures.eosciences 2023, 13, x F P I 12 of 23 i re 10. Sa le 10 i ells f r a i e i t 32. ases f r a ea a i s f 2 a 25 . e t e els re e er te si t e et l r se t iti e t e r effects, t e rel ti s i et ee 32 10 e e e e c e i e ( i ). t st t t t 32 sti i tl fr re le i te sit , t e e er ti f els re ires fr ct re si e s i t t e t t f els st e li te i st r c e res. i re 11. ( eft) Sa le 10 r a ilit istri ti f r a i t 32 f 2 −1, lac as e li es c rres t r al istri ti s si t e ea a S f eac ata set. ( i t) ea ( ) a sta ar e iati ( ) f sa le 10 f ll a li ear relati s i . e res lts s t t t e e st r e i ti f t e s le 10 f ll s li e r rel ti s i . eref re, esti ti t e 10 ri ilit f r i e 32 si r-l r ilit f cti f r i e fr ct re rie t ti is ssi le. i ce 10 32 f l-l li e r rel ti s i , si r ilit f cti s i ter l ti , it is t e ssi le t esti te t e ri ilit f 32 f r i e 10 ( i re ). i re 11. ( eft) Sa le P10 r a ilit istri ti for an input P32 of 2 m− , lac as e li escorrespond to normal distributions using the mean and SD of each data set. (Right) Mean (µ) andstandard deviation (σ) of sampled P10 follow a linear relationship.The results sho that the ean and standard deviation of the sampled P10 followsa linear relationship. Therefore, estimating the P10 variability for any given P32 usi g anormal probability function for a given fracture orientation is possible. Since P10 and P32f llow a linear relationship, using probability functions and interpolation, it is t ssi lt ti t t i ilit 32 f r i e 10 (Figure 12).Geosciences 2023, 13, 187 12 of 22G osciences 2023, 13, x FOR PEER REVIEW 13 of 23 Figure 12. (Left) Example of estimated P32 variability for a P10 value of 2.5 m−1. (Right) estimated P32 probability distribution. 5. Volumetric Intensity Derived from Borehole Intensity It is common practice in the DFN community to calculate linear intensity using drill run lengths and restrict the minimum bias angle to 15° when calculating volumetric in-tensity [13]. This section presents the impact of changing interval length and the minimum bias angle for the calculated P32. A methodology to calculate P32 from borehole intensity is then presented. The main purpose of this methodology is to capture the spatial variation in intensity while avoiding artificially increasing or decreasing the intensity of the inter-vals. 5.1. Effect of the Minimum Bias Angle Considered to Estimate P32 from P10 Using the DFN model presented in Section 4, the volumetric intensity (P32) was cal-culated from the linear intensity (P10), using Equation (4). Note that P32 was calculated using the whole well length, and the analyses focused on evaluating the effect of the min-imum bias angle α in the estimation of P32. For this, three cases were considered: minimum α of 15°, 5° and 1°. When limiting the alpha angle to a minimum of 15°, the calculated P32 is underesti-mated (Table 4), especially in boreholes parallel to the fractures (up to 24% for Well 5). On the other hand, when the minimum bias angle decreases, the ratio starts to become closer to one. It is also worth noticing that the standard deviation (SD) increases when the min-imum angle decreases, suggesting that the dispersion of the data increases for low angles. Table 4 shows a good agreement between the input P32 and the calculated P32 is obtained using a minimum angle of 1°. These results align with the recommendations by Chilès et al. [13] of not discarding data by introducing a minimum bias angle. The problem with this approach is that it is limited to calculating average P32 when several boreholes or scanlines with different orientations are available, but it does not al-low for calculating the spatial variability of intensity along the wells. Table 4. Ratio between input P32 and calculated P32 from borehole intensity. Well Minimum α of 15° Minimum α of 5° Minimum α of 1° Average Ratio SD of Ra-tio Average Ratio SD of Ra-tio Average Ratio SD of Ra-tio Well 1 1.12 0.10 1.04 0.11 1.02 0.12 Well 2 1.15 0.11 1.05 0.12 1.02 0.13 Well 3 1.19 0.11 1.07 0.12 1.03 0.14 Well 4 1.14 0.11 1.05 0.12 1.03 0.13 Well 5 1.24 0.13 1.08 0.14 1.03 0.16 0.00.51.01.52.02.53.03.54.04.55.02 3 4 5 6 7 8P10(m-1)P32 (m-1)P32 Variability for a P10 of 2.5m⁻¹ - All Fracture SizesVariability of 100 Realizations Average of 100 RealizationsP32 Range P32 Variability0.000.010.020.030.040.050.060.070.080.092.0 3.0 4.0 5.0 6.0 7.0 8.0Probability DensityP32 (m-1) P32 Variability for a P10 of 2.5m⁻¹Figure 12. (Left) Example of esti ated P32 variability for a P10 value of 2.5 m−1. (Right) estimatedP32 probability distribution.5. Volumetric Intensity Derived from Borehole IntensityIt is common practice in the DFN community to calculate linear intensity using drillrun lengths and restrict the minimum bias angle to 15◦ when calculating volumetric inten-sity [13]. This section presents the impact of changing interval length and the minimumbias angle for the calculated P32. A methodology to calculate P32 from borehole intensity isthen presented. The main purpose of this methodology is to capture the spatial variation inintensity while avoiding artificially increasing or decreasing the intensity of the intervals.5.1. Effect of the Minimum ias ngle Considered to Estimate P32 from P10Using the DF el presented in Section 4, the volumetric intensity (P32) wascal fr the linear intensity (P10), using Equation (4). Note that P32 was calculatedusing the whole ll length, and the analyses focus d on evaluating the effect of theminimum bias angle α in the estimation of P32. For this, three cases were considered:minimum α of 15◦, 5◦ and 1◦.When limiting the alpha angle to a minimum of 15◦, the calculated P32 is underesti-mated (Table 4), especially in boreholes parallel to the fractures (up to 24% for Well 5). Onthe other hand, when the minimum bias angle decreases, the ratio starts to become closer toone. It is also worth noticing that the standard deviation (SD) increases when the minimumangle decreases, suggesting that the dispersion of the data increases for low angles. Table 4shows a good agreement between the input P32 and the calculated P32 is obtained using aminimum angle of 1◦. These results align with the recommendations by Chilès et al. [13] ofnot discarding data by introducing a minimum bias angle.Table 4. Ratio between input P32 and calculated P32 from borehole intensity.WellMinimum α of 15◦ Mini um α of 5◦ Minimum α of 1◦AverageRatioSD ofRatioAverageRatioSD ofRatioAverageRatioSD ofRatioWell 1 1.12 0.10 1.04 0.11 1.02 0.12Well 2 .15 0.1 1.05 0.12 1.02 0.13Well 3 1.19 0.11 1.07 0.12 1.03 0.14Well 4 1.14 0.11 1.05 0.12 1.03 0.13Well 5 1.24 0.13 1.08 0.14 1.03 0.16Well 6 1.11 0.10 1.04 0.10 1.01 0.11Well 7 1.11 0.10 1.04 0.10 1.02 0.12Well 8 1.16 0.12 1.05 0.12 1.02 0.13Well 9 1.15 0.11 1.05 0.12 1.02 0.13Well 10 1.13 0.10 1.05 0.11 1.03 0.12Total 1.15 0.12 1.05 0.12 1.02 0.13Geosciences 2023, 13, 187 13 of 22The problem with this approach is that it is limited to calculating average P32 whenseveral boreholes or scanlines with different orientations are available, but it does not allowfor calculating the spatial variability of intensity along the wells.5.2. Effect of Interval Length on the Calculated P32 from P10Because of the tendency to calculate linear intensity using run lengths, smaller angleswill result in intervals with high artificial variability in the calculated intensity. Note thatthe practice of using run lengths as reference intervals has no scientific basis and is largelydictated by empirical field practices, which are difficult to change. It is suggested to uselinear intensity values determined for geotechnical domains, which may include severaladjacent cores runs.The effects of varying the interval length on the calculated P32 were investigated usingsimple DFN models with constant fracture orientations. When the fracture orientation andthe sampling length are constant, Equation (4) can be expressed asPˆ32 =1L∑n11sin α=NL sin α=P10sin α(5)A volume of interest was defined within a box of 100 m per side, and DFN modelswere generated using a box size according to the recommendations presented in Section 3,and using the properties presented in Table 5. DFN models were sampled with a 100 m longvertical well. Under these conditions, the α angle is equal to the fracture’s input plunge,and the results are independent of the fracture trend.Table 5. Input parameters used in the models to investigate the effects of interval length for thecalculated P32.Property CommentsOrientation Constant Plunge of 90◦-45◦-30◦-15◦-5◦-2-1◦Spatial model Enhanced Baecher with generation locations at fracture centreIntensity P32 (m−1) 1-2-4-8Size: Lognormal Distribution Radius Mean (Xmean): 2-10-20 (m) and standard deviation(SD) of 40% of the mean radius sizeFracture Shape Hexagon with a constant aspect ratio of 1Number of realizations 10 realizations per model (840 realizations in total from thecombinations of intensity, size and α)The well was discretized in regular intervals of 2 m, 5 m, 10 m, 20 m, 50 m and 100 m.At the same time, boxes centred in each interval, and with the same side length as theintervals were used to calculate the actual P32 per interval (sampled P32). The actual P32 perinterval was then compared with the P32 calculated using Equation (5), based on the P10values sampled per Well interval. Note that for this exercise, estimating the P32 correctionfactor was unnecessary, since the actual P32 in each box was compared with P32 calculatedin each interval. For the same reason, the actual and sample P32 values are higher than theinput P32 values presented in Table 5.Figure 13 presents the results of comparing actual P32 and calculated P32 for differentinterval sizes. For small α angles, the dispersion increases significantly for small sizeintervals, with many values out of the graph scale. On the other hand, for 50 m and 100 mintervals, the values tend to follow the line that defines the linear fit of the data. Notethat the slope of the line is close to one for all cases, even for those with great dispersion.This effect can be explained, since for small α angles the probability of intersecting asmall interval is low. Therefore, many of the small intervals present P32 values of zero.Contrastingly, the intensity for intervals intersecting a fracture is extremely high. Thoseartificial low and high intensities diverge from the actual intensity, but when they areaveraged, they result in P32 values close to the actual intensity (Table 6). This effect wasalready discussed by Hekmatnejad et al. [21], who found that on a composite scale (intervalGeosciences 2023, 13, 187 14 of 22scale), fluctuations in calculated P32 are large and there may be a significant deviationbetween the calculated and the actual value of P32. However, on average, the values tendtoward the actual P32.Geosciences 2023, 13, x FOR PEER REVIEW 15 of 23 Contrastingly, the intensity for intervals intersecting a fracture is extremely high. Those artificial low and high intensities diverge from the actual intensity, but when they are averaged, they result in P32 values close to the actual intensity (Table 6). This effect was already discussed by Hekmatnejad et al. [21], who found that on a composite scale (inter-val scale), fluctuations in calculated P32 are large and there may be a significant deviation between the calculated and the actual value of P32. However, on average, the values tend toward the actual P32. Figure 13. Actual P32 compared to calculated P32 using different interval sizes and α angles. The black line corresponds to the linear fit of the data. Table 6. All cases combined (α = 1, 2, 5, 15, 30, 45 and 90°), using a minimum α of 1°, equivalent to a maximum factor of 57.3. P32 (m−1) Average Actual P32 (m−1) Calculated P32 (m−1) Box Size per Side Well Interval 2m 5m 10m 20m 50m 100m 2m 5m 10m 20m 50m 100m 1 1.2 1.1 1.1 1.1 1.1 1.1 1.2 1.2 1.2 1.2 1.2 1.2 2 2.3 2.3 2.2 2.2 2.2 2.2 2.3 2.3 2.3 2.3 2.3 2.3 4 4.5 4.5 4.4 4.4 4.4 4.4 4.5 4.5 4.5 4.5 4.5 4.5 8 8.9 8.9 8.9 8.8 8.8 8.8 8.9 8.9 8.9 8.9 8.9 8.9 1 1.2 1.1 1.1 1.1 1.1 1.1 1.2 1.2 1.2 1.2 1.2 1.2 2 2.3 2.3 2.2 2.2 2.2 2.2 2.3 2.3 2.3 2.3 2.3 2.3 It can be observed that even when the P32 calculated in small intervals present a high dispersion, their mean tent to represent the actual P32 (Table 6). This interesting result sug-gests that it would be possible to reduce the dispersion by using longer intervals for small α angles and smaller intervals for greater α angles. To investigate the appropriate length to use for each α angle, the maximum calcu-lated P32 was plotted for each interval length. Figure 14 shows the maximum calculated Figure 13. Actual P32 compare t calculated P32 using different i terval sizes and α angles. Theblack line corresponds to t e li ear fit of the data.Table 6. All cases co i e (α = 1, 2, 5, 15, 30, 45 and 90◦), using a mini um α of 1◦, equivalent to amaxi factor of 57.3.P32 (m−1)Average Actual P32 (m−1) Calculated P32 (m−1)Box Size per Side Well Interval2 m 5 m 10 m 20 m 50 m 100 m 2 m 5 m 10 m 20 m 50 m 100 m1 1.2 1.1 1.1 1.1 1.1 1.1 1.2 1.2 1.2 1.2 1.2 1.22 2.3 2.3 2.2 2.2 2.2 2.2 2.3 2.3 2.3 2.3 2.3 2.34 4.5 4.5 4.4 4.4 4.4 4.4 4.5 4.5 4.5 4.5 4.5 4.58 8.9 8.9 8.9 8.8 8.8 8.8 8.9 8.9 8.9 8.9 8.9 8.91 1.2 1.1 1.1 1.1 1.1 1.1 1.2 1.2 1.2 1.2 1.2 1.22 2.3 2.3 2.2 2.2 2.2 2.2 2.3 2.3 2.3 2.3 2.3 2.3It can be observed that even when the P32 calculated in small intervals present a highdispersion, their mean tent to represent the actual P32 (Table 6). This interesting resultsuggests that it would be possible to reduce the dispersion by using longer intervals forsmall α angles and smaller intervals for greater α angles.To investigate the appropriate length to use for each α angle, the maximum calculatedP32 was plotted for each interval length. Figure 14 shows the maximum calculated P32 perinterval for inpu intensities of 1 m−1 and 8 m−1. The data suggest that for α angles greaterthan 15 ◦, intervals with lengths ranging from 10 m to 15 m can be used. While for α anglesGeosciences 2023, 13, 187 15 of 22between 2◦ and 15 ◦, intervals with a length between 20 and50 m need to be used; finally,for α angles between 1◦ and 2◦, an interval length greater than 50 m is required.Geosciences 2023, 13, x FOR PEER REVIEW 16 of 23 P32 per interval for input intensities of 1 m−1 and 8 m−1. The data suggest that for α angles greater than 15 °, intervals with lengths ranging from 10 m to 15 m can be used. While for α angles between 2° and 15 °, intervals with a length between 20 and50 m need to be used; finally, for α angles between 1° and 2°, an interval length greater than 50 m is required. Figure 14. Maximum calculated P32 per interval length for input intensities of 1 m−1 and 8 m−1. 5.3. Proposed Methodology to Calculate P32 from Borehole Intensity Based on the analyses performed, a methodology to calculate P32 from borehole in-tensity is proposed. This methodology can be easily implemented using code, and the idea behind this is to allow the calculation of P32 in intervals small enough to capture the spatial variation in intensity, but at the same time without artificially increasing or decreasing the intensity of the interval. This is very useful when the interval intensity is used as input in block models, in which artificial changes in intensity may affect the result of the interpo-lation values in the block model. The following methodology is proposed to calculate P32 from borehole intensity: • Since P32 is an additive variable, it is possible to discretize a borehole in different overlapping intervals and calculate P32 as the addition of the P32 values calculated in each interval, as long as the fractures belonging to each interval are not double counted. • Based on the foregoing, it is possible to calculate P32 using the correction proposed by Chilès et al. [13] (Equation (4)), using regular intervals whose length depends on 0204060801001201400 20 40 60 80 100 120Max. Calcualted P32(m-1)Interval Length (m)Maximum Calculated P32 per Interval Length (Input P32= 8 m⁻¹)α =1°α =2°α =5°α =15°α =30°α =45°α =90°0102030405060700 20 40 60 80 100 120Max. Calcualted P32(m-1)Interval Length (m)Maximum Calculated P32 per Interval Length (Input P32= 1 m⁻¹)α =1°α =2°α =5°α =15°α =30°α =45°α =90°Figure 14. Maximum calculated P32 per interval length for input intensities of 1 m−1 and 8 m−1.5.3. Proposed Methodology to Calculate P32 from Borehole IntensityBased on the analyses performed, a methodology to calculate P32 from boreholeintensity is proposed. This methodology can be easily implemented using code, and theidea behind this is to allow the calculation of P32 in intervals small enough to capturethe spatial variation in intensity, but at the same time without artificially increasing ordecreasing the intensity of the interval. This is very useful when the interval intensity isused as input in block models, in which artificial changes in intensity may affect the resultof the interpolation values in the block model.The following methodology is proposed to calculate P32 from borehole intensity:• Since P32 is an additive variable, it is possible to discretize a borehole in differentoverlapping intervals and c lculate P32 a the addition of the P32 values calculated ineach interval, as long s the fractures belonging to each interval are not double counted.• Based on the foregoing, i is possibl to calculate P32 using the correction propos dby Chilès et al. [13] (Equation (4)), using regular intervals whose length depends onth magnitude of the acute angle (α) between the scanline and the fracture. Table 7present the recommended length intervals for each angle.• Add the P32 values calculated using the smallest intervals; this addition correspondsto the total P32.Geosciences 2023, 13, 187 16 of 22Table 7. Recommended interval length to calculate P32.α Angle (◦) Interval Length to Calculate P32 Comments15≤ α 10 m to 20 m Use a length similar to the size of the block model cells2≤ α <15 20 m to 50 m P32 must be calculated within the same Structural Domain1≤ α <2 50 m to total borehole length It can be extended at the full well length; P32 must be calculatedwithin the same Structural DomainNote that if only the average intensity is required, it is recommended to calculateP32 using the whole length of the borehole or the length within the structural domain ofinterest. It is still recommended to limit the minimum α angle to 1◦, to avoid extremehigh-weighting factors or division by zero.The proposed methodology was tested using the DFN model and synthetic wellspresented in Section 4. The interval lengths used, depending on the acute angle betweenthe well and the fracture, are presented in Table 8.Table 8. Intervals used to test the methodology proposed to calculate P32 depending on α angle.α Angle (◦) Interval Length to Calculate P32 Maximum Weighing Factor per Fracture15≤ α 15 m 3.92≤ α <15 30 m 28.71≤ α <2 90 m (well total length) 57.3Independent of the borehole direction, a good agreement was observed between theinput P32 and the median (50% percentile) of the calculated P32. Figure 15 presents acomparison between the proposed methodology and the P32 values calculated in Well 5using different interval sizes and minimum α angles, while Table 9 summarizes the resultobtained for each one of the 10 wells analyzed.Geosciences 2023, 13, x FOR PEER REVIEW 18 of 23 Figure 15. P32 values were calculated in Well 5 using different interval sizes and minimum α angles. Table 9. Comparison between the proposed methodology and the P32 values calculated in each one of the 10 wells using different interval sizes and minimum α angles. Min. α (°) Int. Length (m) Well 1 2 3 4 5 6 7 8 9 10 1° 3 Linear Fit Factor (1) 1.00 1.00 1.01 1.01 1.01 1.00 1.00 1.01 1.00 1.01 R2 0.47 0.44 0.41 0.44 0.37 0.48 0.48 0.43 0.43 0.47 % Varia-tion Me-dian (2) Min. 6% 7% 9% 7% 10% 6% 6% 8% 7% 7% Average 11% 13% 15% 13% 18% 11% 10% 14% 13% 12% Max. 24% 25% 28% 28% 38% 22% 20% 29% 25% 22% 15° 3 Linear Fit Factor (1) 1.11 1.14 1.17 1.13 1.22 1.11 1.11 1.15 1.14 1.12 R2 0.61 0.61 0.59 0.59 0.55 0.61 0.62 0.59 0.58 0.63 % Varia-tion Me-dian (2) Min. 12% 14% 17% 14% 20% 11% 11% 15% 14% 12% Average 15% 17% 20% 17% 24% 14% 14% 18% 17% 15% Max. 24% 25% 28% 28% 38% 22% 20% 29% 25% 22% 1° 15 Linear Fit Factor (1) 1.00 1.00 1.01 1.01 1.01 1.00 1.00 1.01 1.00 1.01 R2 0.81 0.80 0.77 0.80 0.75 0.82 0.82 0.79 0.79 0.82 % Varia-tion Me-dian (2) Min. 2% 2% 3% 3% 3% 2% 2% 3% 2% 3% Average 4% 5% 6% 5% 6% 4% 4% 5% 4% 5% Max. 10% 10% 13% 11% 15% 9% 8% 11% 9% 9% 15° 15 Linear Fit Factor (1) 1.11 1.14 1.17 1.13 1.22 1.11 1.11 1.15 1.14 1.12 R2 0.89 0.88 0.87 0.88 0.86 0.89 0.89 0.87 0.88 0.89 % Varia-tion Me-dian (2) Min. 10% 13% 15% 12% 19% 10% 10% 13% 13% 11% Average 11% 13% 16% 13% 20% 11% 10% 14% 13% 11% Max. 14% 16% 19% 16% 23% 13% 11% 16% 15% 13% Figure 15. P32 values were calculated in Well 5 using different interval sizes and minimum α angles.Geosciences 2023, 13, 187 17 of 22Table 9. Comparison between the proposed methodology and the P32 values calculated in each oneof the 10 wells using different interval sizes and minimum α angles.Min.α (◦)Int.Length(m)Well 1 2 3 4 5 6 7 8 9 101◦ 3Linear Fit Factor (1) 1.00 1.00 1.01 1.01 1.01 1.00 1.00 1.01 1.00 1.01R2 0.47 0.44 0.41 0.44 0.37 0.48 0.48 0.43 0.43 0.47% VariationMedian (2)Min. 6% 7% 9% 7% 10% 6% 6% 8% 7% 7%Average 11% 13% 15% 13% 18% 11% 10% 14% 13% 12%Max. 24% 25% 28% 28% 38% 22% 20% 29% 25% 22%15◦ 3Linear Fit Factor (1) 1.11 1.14 1.17 1.13 1.22 1.11 1.11 1.15 1.14 1.12R2 0.61 0.61 0.59 0.59 0.55 0.61 0.62 0.59 0.58 0.63% VariationMedian (2)Min. 12% 14% 17% 14% 20% 11% 11% 15% 14% 12%Average 15% 17% 20% 17% 24% 14% 14% 18% 17% 15%Max. 24% 25% 28% 28% 38% 22% 20% 29% 25% 22%1◦ 15Linear Fit Factor (1) 1.00 1.00 1.01 1.01 1.01 1.00 1.00 1.01 1.00 1.01R2 0.81 0.80 0.77 0.80 0.75 0.82 0.82 0.79 0.79 0.82% VariationMedian (2)Min. 2% 2% 3% 3% 3% 2% 2% 3% 2% 3%Average 4% 5% 6% 5% 6% 4% 4% 5% 4% 5%Max. 10% 10% 13% 11% 15% 9% 8% 11% 9% 9%15◦ 15Linear Fit Factor (1) 1.11 1.14 1.17 1.13 1.22 1.11 1.11 1.15 1.14 1.12R2 0.89 0.88 0.87 0.88 0.86 0.89 0.89 0.87 0.88 0.89% VariationMedian (2)Min. 10% 13% 15% 12% 19% 10% 10% 13% 13% 11%Average 11% 13% 16% 13% 20% 11% 10% 14% 13% 11%Max. 14% 16% 19% 16% 23% 13% 11% 16% 15% 13%ProposedMethodologyLinear Fit Factor (1) 1.00 1.00 1.01 1.01 1.01 1.00 1.00 1.01 1.00 1.01R2 0.88 0.87 0.87 0.87 0.85 0.89 0.88 0.87 0.87 0.88% VariationMedian (2)Min. 1% 1% 1% 1% 1% 0% 1% 1% 1% 2%Average 2% 2% 3% 3% 2% 1% 2% 2% 2% 3%Max. 5% 5% 6% 6% 6% 4% 4% 6% 2% 5%(1) Linear fit factor between P32 and calculated P32 (Input P32 = factor x Calculated P32), note that this factor isequivalent to 1/slope of the linear fit lines presented on Figure 15; (2) Percentage variation between input P32 andmedian (Percentile 50%) of the calculated P32.The methodology herein proposed is the one that provides the best result in termsof a good agreement between input P32 and calculated P32, and at the same time main-taining a relatively low variability. Note that when a minimum angle of 15◦ is consideredthe mean P32 values calculated are underestimated. This is especially relevant consider-ing it is common practice to use a minimum angle of 15◦ [13] to restrict the maximumweighting factor.An interval of 3 m corresponds to the typical length of a drill run. It is commonpractice (but not necessarily best practice) to calculate linear intensity on a run basis andthen estimate P32 using the same run length. The analyses show that using a length of3 m, even when restricting the minimum angle to 15 ◦, produces a high variability on thecalculated P32, meaning that the calculated P32 will often be under- or overestimated. Notethat, in practice, depending on the drilling operation and rock quality, the actual drillingrun lengths are variable and may be significantly smaller than 3 m. This means that if P32Geosciences 2023, 13, 187 18 of 22is calculated using those shorter intervals, the P32 variability may be much higher than theone calculated using 3 m intervals. These results agree with the work by Elmo and Steadand Yang et al. [22,23] on the problems of using run lengths to calculate RQD (rock qualitydesignation [24], which is an empirical rock mass quality parameter linked to fracturefrequency (i.e., linear intensity).As another validation method, P32 values were calculated using a grid with a totalvolume and location equal to the volume of interest, using cells of 10 m per side. In thecells, P32 is calculated as the area of fractures within the grid cell divided by the grid cellvolume. Figure 16 presents an example of the P32 values calculated in the grid for onerealization using an input P32 of 6 m−1 and a cross-section with a comparison between thegrid values and the values calculated in the wells.Geosciences 2023, 13, x FOR PEER REVIEW 19 of 23 Proposed Meth-odology Linear Fit Factor (1) 1.00 1.00 1.01 1.01 1.01 1.00 1.00 1.01 1.00 1.01 R2 0.88 0.87 0.87 0.87 0.85 0.89 0.88 0.87 0.87 0.88 % Varia-tion Me-dian (2) Min. 1% 1% 1% 1% 1% 0% 1% 1% 1% 2% Average 2% 2% 3% 3% 2% 1% 2% 2% 2% 3% Max. 5% 5% 6% 6% 6% 4% 4% 6% 2% 5% (1) Linear fit factor between P32 and calculated P32 (Input P32= factor x Calculated P32), note that this factor is equivalent to 1/slope of the linear fit lines presented on Figure 15 (2) Percentage variation between input P32 and median (Percentile 50%) of the calculated P32. As another validation method, P32 values were calculated using a grid with a total volume and location equal to the volume of interest, using cells of 10 m per side. In the cells, P32 is calculated as the area of fractures within the grid cell divided by the grid cell volume. Figure 16 presents an example of the P32 values calculated in the grid for one realization using an input P32 of 6 m−1 and a cross-section with a comparison between the grid values and the values calculated in the wells. When the proposed methodology is applied, it is possible t capture the spatial var-i tion of P32 nd simultaneously obtain a good agreement between the actual P32 (P32 val-ues calculated in the grid) and the P32 calculated along the wells. On the other hand, when P32 is calculated in 3 m intervals, high variability is observed, with artificially low or high values that do not agree with the actual values calculated on the grid cells. Figure 17 pre-sents the result of 100 realizations as cumulative frequency curves of actual P32 in the grid cells and the calculated P32 in wells. The proposed methodology is the one that produces the best results, especially when compared with the common practice of limiting the min-imum angle to 15° and calculating P32 using the drill run length (3 m). Figure 16. P32 values calculated in the grid for one realization, cross-sections, and comparison with calculated P32 values in wells. Actual P32 Calculated in GridCross-Section View from SW – Proposed MethodologyCross-Section View from SW – Min. α to 15°in 3 m IntervalsActual P32 m-1 In CellsCalculated P32 m-1 in WellsPlan View – Cros-Section Location and WellsCros Section Location Presented in Colour RedFigure 16. P32 values calculated in the grid for one realization, cross-sections, and comparison withcalculated P32 values in wells.When the proposed methodology is applied, it is possible to capture the spatialvariation of P32 and simultaneously obtain a good agreement between the actual P32 (P32values calculated in the grid) and the P32 calculated along the wells. On the other hand,when P32 is calculated in 3 m intervals, high variability is observed, with artificially low orhigh values that do not agree with the actual values calculated on the grid cells. Figure 17presents the result of 100 realizations as cumulative frequency curves of actual P32 inthe grid cells and the calculated P32 in wells. The proposed methodology is the one thatproduces the best results, especially when compared with the common practice of limitingthe minimum angle to 15◦ and calculating P32 using the drill run length (3 m).Geosciences 2023, 13, 187 19 of 22Geosciences 2023, 13, x FOR PEER REVIEW 20 of 23 Figure 17. Cumulative frequency curves for actual P32 in grid cells and calculated P32 in wells for an input P32 of 6 m−1 and a mean fracture radius of 20 m. The values presented correspond to the results of 100 realizations. 6. Summary and Conclusions A main component of rock mass characterization is the characterization of disconti-nuities, which play a major role in rock mass’s mechanical and hydraulic properties. This paper focused on volumetric fracture intensity (P32), which remains one of the key prop-erties of fracture characterization for DFN modelling, albeit it is not possible to measure it directly. The problem of estimating P32 from information observed along boreholes was ad-dressed by generating a series of DFN models and then comparing the volumetric inten-sity of those models with the volumetric intensity calculated from borehole intensity, us-ing a methodology based on the methodology proposed by Chilès et al. [13]. This paper first investigated the boundary effects on DFN models, and the authors have proposed a methodology to mitigate boundary effects for P32. Subsequently, the pa-per has focused on presenting a new method for calculating P32 using borehole fracture intensity. This methodology allows the capturing of the spatial variation in intensity while at the same time avoiding increasing or decreasing artificially P32 calculated in borehole intervals. This is very useful when the interval intensity is used as input for block model-ling of P32 values. The main findings of this research can be summarized as follows: • The analyses confirmed that DFN models present boundary effects on intensity, something that was already discussed by Priest [15], but that it is often disregarded when building DFN models. The boundary effect on the fracture intensity and other derived analyses may be considerable. • It is possible to minimize the boundary effect for DFN models when the generation box is defined as a function of the fracture size and the dimensions of the volume of interest. A correction factor for the input P32 can be defined to correct the intensity generated in the volume of interest. This factor will depend on the size of the volume of interest, the size of the generation box, the fracture size distribution and the frac-ture orientations. • The analysis showed that when the boundary effects are mitigated, the relationship between P10 and P32 is linear and independent of the fracture size. This implies that it is possible to estimate volumetric intensity directly from borehole data. Note that this 0%10%20%30%40%50%60%70%80%90%100%0 2 4 6 8 10 12Cumulative FrequencyP32 (m-1) Grid Values for an Input P32 of 6 m-1 and Mean Radius of 20 mGrid Proposed Methodology Min. α: 1°in 3m IntervalsMin. α: 15°in 3m Intervals Min. α: 1°in 15m Intervals Min. α: 15°in 15m IntervalsFigure 17. Cumulative frequency curves for actual P32 in grid cells and calculated P32 in wells foran input P32 of 6 m−1 and a mean fr ct re radius of 20 m. The values presented correspond to theresults of 100 realizations.6. Summary and ConclusionsA main component of rock mass characterization is the characterization of disconti-nuities, which play a major role in rock mass’s mechanical and hydraulic properties. Thispaper focused on volumetric fracture intensity (P32), which remains one of the key prop-erties of fracture characterization for DFN modelling, albeit it is not possible to measureit directly.The problem of estimating P32 from information observed along boreholes was ad-dressed by generating a series of DFN models and then comparing the volumetric intensityof those models with the volumetric intensity calculated from borehole intensity, using amethodology based on the methodology proposed by Chilès et al. [13].This paper first investigated the boundary effects on DFN models, and the authorshave proposed a methodology to mitigate boundary effects for P32. Subsequently, thepaper has focused on presenting a new method for calculating P32 using borehole fractureintensity. This methodology allows the capturing of the spatial variation in intensity whileat the same time avoiding increasing or decreasing artificially P32 calculated in boreholeintervals. This is very useful when the interval intensity is used as input for block modellingof P32 values.The main findings of this research can be summarized as follows:• The analyses confirmed that DFN models present boundary effects on intensity, some-thing that was already discussed by Priest [15], but that it is often disregarded whenbuilding DFN models. The boundary effect on the fracture intensity and other derivedanalyses may be considerable.• It is possible to minimize the boundary effect for DFN models when the generation boxis defined as a function of the fracture size and the dimensions of the volume of interest.A correction factor for the input P32 can be defined to correct the intensity generated inthe volume of interest. This factor will depend on the size of the volume of interest, thesize of the generation box, the fracture size distribution and the fracture orientations.• The analysis showed that when the boundary effects are mitigated, the relationshipbetween P10 and P32 is linear and independent of the fracture size. This implies thatit is possible to estimate volumetric intensity directly from borehole data. Note thatthis result should not be co strued to suggest that fracture size is not an importantarameter; indeed, fract re size remains a required input for generating a DFN model.Geosciences 2023, 13, 187 20 of 22• The analysis showed that simulations can quantify the variability and distribution ofthe volumetric intensity for a given linear intensity. Note that more than a unique P32value can be obtained for a given P10 value, but with this methodology, it is possibleto quantify the probability of obtaining a certain P32 for a given P10.• The general practice of constraining the minimum bias angle to 15◦ and calculat-ing P32 using small intervals (drill run lengths) results in an underestimation anda high variability in the calculated P32. When using the methodology proposed byChilès et al. [13] to calculate P32, the interval’s size greatly impacts the variability ofthe calculated P32. Using small intervals will artificially increase the variability of thecalculated P32, even when the average P32 tends to be close to the actual P32. On theother hand, when limiting the minimum angle between the fracture and the boreholeto a minimum of 15◦, the mean calculate P32 values are underestimated.• The proposed methodology to calculate P32 using variable lengths, depending onthe angle between the fractures and the borehole, allows the capture of the spatialvariation in intensity and at the same time avoid increasing or decreasing artificiallythe intensity of the intervals. This is very useful when the interval intensity is used asinput for interpolating P32 values in block models.• When only average volumetric intensity is required, it is recommended to calculateP32 using the total length of the borehole or the length within the structural domain ofinterest. In this case, limiting the minimum α angle to 1◦ is recommended to avoidextremely high-weighting factors or division by zero.Limitations and AssumptionsThe research described previously is subject to the following limitations and assumptions:• The main limitation of this study is that the analyses presented are purely simulatedand have not yet been applied to real data. The concept of validating modelling resultswith field data does not apply to P32 values since it is impossible to directly comparecalculated P32 with the actual P32 of the rock mass (which remains unknown).• All fractures have been to be assumed planar and extremely thin. Although fracturesin nature can be curved and present thickness, it is common practice to assume thatdiscontinuities are planar [25,26].• It was assumed that the fractures are hexagonal with a constant aspect ratio of 1.Although the assumption that fractures are equidimensional is commonly used, someinvestigators indicate that fractures are in reality not equidimensional [26].• It was assumed in the models that fractures are randomly and independently dis-tributed in space. In reality, fracture intensity may be dependent on geological features;for example, a common geological observation is that the distance to major faultscontrols fracture intensity [27].• The analysis performed assumed that all fractures intersected by the synthetic wellswere identified. In practice, there are biases when gathering fracture intensity fromboreholes. For example, besides the orientation bias, it is common that during boreholelogging, not all natural discontinuities are counted, and some artificial fractures (e.g.,mechanical breaks) are counted as natural fractures. Even if the fracture intensity iscalculated using televiewer data, the quality of the data would depend on the qualityof the images and the interpretation performed for each discontinuity.• The methodology proposed to minimize boundary effects may be impracticable whenthe fractures are too big. For big fractures, it would be necessary to generate anextremely big generation box, that may require extensive computation time.Author Contributions: Conceptualization, P.O., methodology, P.O. and D.E.; formal analysis, P.O.;investigation, P.O. and D.E.; writing—original draft preparation, P.O.; writing—review and editing,P.O., D.E., S.R. and A.B.; supervision, D.E. All authors have read and agreed to the published versionof the manuscript.Funding: This research received no external funding.Geosciences 2023, 13, 187 21 of 22Data Availability Statement: Data available upon request pending confidentiality of some of theinput data.Acknowledgments: We would to thank WSP Golder for providing an academic license for FracMan.Conflicts of Interest: The authors declare no conflict of interest.References1. Dershowitz, W.S.; Einstein, H.H. Characterizing rock joint geometry with joint system models. Rock Mech. Rock Engng. 1988, 21,21–51. [CrossRef]2. Elmo, D.; Liu, Y.; Rogers, S. Principles of discrete fracture network modelling for geotechnical applications. In Proceedings of theInternational Discrete Fracture Networking Conference, DFNE 2014, Vancouver, BC, Canada, 19–23 October 2014; p. 238.3. Elmo, D.; Stead, D. An integrated numerical modelling—Discrete fracture network approach applied to the characterisation ofrock mass strength of naturally fractured pillars. Rock Mech. Rock Eng. 2010, 43, 3–19. [CrossRef]4. Pierce, M.; Cundall, P.; Potyondy, P.; Mas Ivars, D. A synthetic rock mass model for jointed rock. In Proceedings of the 1stCanada—US Rock Mechanics Symposium, Vancouver, BC, Canada, 27–31 May 2007; Volume 1, pp. 341–349.5. Rogers, S.; Elmo, D.; Beddoes, R.; Dershowitz, B. Mine scale DFN modelling and rapid upscaling in geomechanical simulations oflarge open pits. In Proceedings of the International Symposium on Rock Slope Stability in Open Pit Mining and Civil Engineering,Santiago, Chile, 9–11 November 2009.6. Rogers, S.F.; Elmo, D.; Webb, G.; Guajardo, C. DFN Modelling of Major Structural Instabilities in a Large Open Pit for Endof Life Planning Purposes. In Proceedings of the 50th U.S. Rock Mechanics/Geomechanics Symposium, Houston, TX, USA,26–29 June 2016.7. Brzovic, A.; Rogers, S.; Webb, G.; Hurtado, J.P.; Marin, N.; Schachter, P.; Alvarez, J.; Baraona, K. Discrete fracture networkmodelling to quantify rock mass pre-conditioning at the El Teniente Mine, Chile. Min. Technol. 2015, 124, 163–177. [CrossRef]8. Dershowitz, W.S.; Finnila, A.; Rogers, S.; Hamdi, P.; Moffitt, K.M. Step Path Rock Bridge Percentage for Analysis of Slope Stability.In Proceedings of the 51st U.S. Rock Mechanics/Geomechanics Symposium, San Francisco, CA, USA, 25–28 June 2017.9. Rogers, S.; Elmo, D.; Webb, G.; Catalan, A. Volumetric Fracture Intensity Measurement for Improved Rock Mass Characterisationand Fragmentation Assessment in Block Caving Operations. Rock Mech. Rock Eng. 2015, 48, 633–649. [CrossRef]10. Dershowitz, W.S.; Herda, H.H. Interpretation of fracture spacing and intensity. In Proceedings of the 33rd U.S. Rock Mechan-ics/Geomechanics Symposium, Santa Fe, NM, USA, 8–10 June 1992.11. Zhang, L.; Einstein, H.H. Estimating the intensity of rock discontinuities. Int. J. Rock Mech. Min. Sci. 2000, 37, 819–837. [CrossRef]12. Wang, X. Stereological Interpretation of Rock Fracture Traces on Borehole Walls and Other Cylindrical Surfaces. DoctorDissertation, Faculty of the Virginia Polytechnic Institute, State University, Blacksburg, VA, USA, 2005.13. Chilès, J.P.; Wackernagel, H.; Beucher, H.; Lantuéjoul, C.; Elion, P. Estimating fracture density from a linear or areal survey.In Proceedings of the Eighth International Geostatistics Congress, Santiago, Chile, 1–5 December 2008; Ortiz, J.M., Emery, X., Eds.;Gecamin Ltda: Santiago, Chile, 2008; pp. 535–544.14. Terzaghi, R. Sources of error in joint surveys. Geotechnique 1965, 15, 287–303. [CrossRef]15. Priest, S.D. Discontinuity Analysis for Rock Engineering; Chapman and Hall: London, UK; Springer: Dordrecht, The Netherlands,1993.16. Rocsience. Dip User Guide. 2022. Available online: https://www.rocscience.com/help/dips/documentation/stereonet-plots/terzaghi-weighting (accessed on 25 April 2023).17. Samaniego, J.A.; Priest, S.D. The prediction of water flows through discontinuity networks into underground excavation. In Pro-ceedings of the Design and Performance of Underground Excavations: ISRM Symposium, Cambridge, UK, 3–6 September 1984;pp. 157–164.18. Junkin, W.R.; Ben-Awuah, E.; Fava, L. DFN Variability Analysis Through Voxelization. In Proceedings of the 53rd U.S. RockMechanics/Geomechanics Symposium, New York, NY, USA, 23–26 June 2019.19. Golder Associates. Golder Associates (UK) Ltd. FracMan 8.0 Manual; Golder Associates: Bourne End, UK.20. Baecher, G.B.; Lanney, N.A.; Einstein, H.H. Statistical Description of Rock Properties and Sampling. In Proceedings of the 18thU.S. Symposium on Rock Mechanics, Golden, CO, USA, 22–24 June 1977; p. ARMA-77-0400.21. Hekmatnejad, A.; Emery, X.; Brzovic, A.; Schachter, P.; Vallejos, J.A. Spatial modeling of discontinuity intensity from boreholeobservations at El Teniente mine, Chile. Eng. Geol. 2017, 228, 97–106. [CrossRef]22. Elmo, D.; Stead, D. The concept of representative elementary length (REL) as an effective tool to study scale effects in rockengineering problems. Boletín Geológico Min. 2021, 131, 355–362. [CrossRef]23. Yang, B.; Mitelman, A.; Elmo, D.; Stead, D. Why the future of rock mass classification systems requires revisiting its empiricalpast. Q. J. Eng. Geol. Hydrogeol. 2022, 55. [CrossRef]24. Deere, D.U.; Hendron, A.J.; Patton, F.D.; Cording, E.J. Design of surface and near surface construction in rock. In Failure andBreakage of Rock, Proceedings of the 8th U.S. Symposium of Rock Mechanics, Minneapolis, MN, USA, 15–17 September 1966; Fairhurst, C.,Ed.; OnePetro: New York, NY, USA, 1967; pp. 237–302.25. Warburton, P.M. A stereological interpretation of joint trace data. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1988, 17, 181–190.[CrossRef]Geosciences 2023, 13, 187 22 of 2226. Zhang, L.; Einstein, H.H. Estimating the mean trace length of rock discontinuities. Rock Mech. Rock Engng. 1998, 31, 217–235.[CrossRef]27. McCaffrey, K.J.W.; Sleight, J.M.; Pugliese, S.; Holdsworth, R.E. Fracture formation and evolution in crystalline rocks: Insightsfrom attribute analysis. Geol. Soc. Lond. Spec. Publ. 2003, 214, 109–124. [CrossRef]Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individualauthor(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury topeople or property resulting from any ideas, methods, instructions or products referred to in the content.
Article
eng
Reviewed
Vancouver : University of British Columbia Library
Multidisciplinary Digital Publishing Institute
CC BY 4.0
Faculty
Researcher
Other
volumetric fracture intensity
discrete fracture network
borehole data
boundary effects
Discrete Fracture Network (DFN) Analysis to Quantify the Reliability of Borehole-Derived Volumetric Fracture Intensity
Text