Discrete Fracture Network (DFN) Analysis to Quantify the Reliability of Borehole-Derived Volumetric Fracture Intensity

Abstract

Volumetric fracture intensity (P₃₂) 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 P₃₂ at a resolution compatible with geotechnical applications is not yet possible. As a result, P₃₂ 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 P₃₂ using information from boreholes (1D data). When calculating P₃₂ 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 P₃₂. Additionally, the size of the interval has a significant impact on the variability of the calculated P₃₂. We propose a methodology to calculate P₃₂ 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.

1. Introduction

A key component of understanding the behaviour of a rock mass is the characterization of 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 than individually [1], using specific properties, namely discontinuity orientation, size, shape, spatial location, and intensity.

The necessity of modelling discontinuities for different engineering and geological applications 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, as discussed in [2]. Note that the term fracture is herein applied in a broad context to refer to a variety of discontinuities at different scales (faults, joints, veins, etc.). DFN modelling has become increasingly popular in recent years amongst geotechnical practitioners and engineers, including the generation of synthetic rock mass properties [3,4], geomechanical simulation of open pits [5,6], quantification of rock mass pre-conditioning [7], estimation of rock bridge percentage for stability analysis [8] and rock mass fragmentation and calculation of fragmentation distribution at cave mine scale [9].

One of the most important inputs for DFN analysis is fracture intensity, which can be expressed in various ways, including (i) the total sum of fracture surface area per unit of volume, P₃₂; (ii) the total sum of fracture length per unit area, P₂₁; and (iii) the number of fractures per unit of length, P₁₀.

P₃₂ is a nondirectional intrinsic measure of fracture intensity, and these characteristics make P₃₂ the preferred measure of fracture intensity in DFN modelling [9]. While P₃₂ cannot be directly measured, it can be obtained from P₁₀ or P₂₁ 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.

Zhang and Einstein [11] proposed a methodology to estimate the volumetric intensity 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 assumptions about 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, it requires exposed rock faces, which are not always available, particularly in the initial stages of a project. In addition, even if surface data are available, this may not represent the underground conditions. Wang [12] proposed a numerical approximation between P₁₀ and P₃₂; the limitation of this methodology is that it assumes that the fracture population in a single set follows a Univariate Fisher hemispherical probability distribution, with the method not being suitable for use with other types of probability distributions. To overcome those limitations Chilès et al. [13] developed a methodology based on Terzaghi Weighting [14] to estimate P₃₂ directly from borehole intensity. Chilès’ method [13] has been widely adopted in the DFN community. This paper specifically addresses the question remaining of using self-assumed practical industry standards to excuse the lack of geological information (e.g., fracture size).

2. Chilès’ Methodology to Calculate P₃₂ 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 observe. 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 underrepresented by a factor of “sinα”; this factor is known as the Terzaghi Weighting [14].

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.

$$d=\frac{N}{L \mathrm{s}\mathrm{i}\mathrm{n} \alpha }$$

(1)

According to Chilès et al. [13], Equation (1) can be rewritten in terms of the predicted volumetric intensity P₃₂.

$${\widehat{P}}_{32}=\frac{1}{L}{\sum }_1^n\frac{1}{\left|\mathrm{c}\mathrm{o}\mathrm{s}({\omega }_i-\beta )\right|}$$

(2)

Chilès et al. [13] presented a generalization of Equation (2) for the case in which fractures are sampled along a line with varying orientation (βi) or several lines with total length L:

$${\widehat{P}}_{32}=\frac{1}{L}{\sum }_1^n\frac{1}{\left|\mathrm{c}\mathrm{o}\mathrm{s}({\omega }_i-{\beta }_i)\right|}$$

(3)

where L is the length of the scanline, βi is the local orientation of the scanline at the location of 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 unit vectors, β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 the fracture i (αi) as

$${\widehat{P}}_{32}=\frac{1}{L}{\sum }_1^n\frac{1}{\mathrm{s}\mathrm{i}\mathrm{n}\left({\alpha }_i\right)}$$

(4)

Since the weighting factor tends to infinity as α approaches zero, applying a maximum limit to α (minimum bias angle) is expected. Priest [15] suggests using a maximum weighting of 10; this weight corresponds to an α value of approximately 6°, while other authors recommend an α angle between 5° and 20° [16]. In practice, a limit of 15° is often used to improve the robustness of the estimator. However, this has the inconvenience of introducing some bias by discarding valid data [13].

3. DFN Models and Boundary Effects

The engineering community often disregards boundary effects when generating DFN models, even though their impact can be significant when the volume used to create the DFN 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 the volume of interest to minimize these boundary effects. Samaniego and Priest [17] indicated that a generation volume of four times the volume of interest would be sufficient to reduce boundary effects. Note that depending on the model size, this increase in the size of the generation volume may significantly increase the computation time. On the other hand, Junkin et al. [18] indicate that it could be possible to mitigate boundary effects by including mapped traces or fracture seeds on the exterior boundaries of a DFN region. They proposed to sample the DFN model using planes to estimate P₂₁ intensity in areas where the boundary effect does not exist and then use this intensity to stochastically seed traces on the boundary of the DFN.

This Section is dedicated to studying the boundary effects that may influence the intensity 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 on DFN models’ intensity.

3.1. DFN Modelling to Investigate Boundary Effects on Intensity

To determine and quantify the impact of boundary effects on volumetric fracture intensity (P₃₂), a set of DFN models was generated within a box of 150 m per side. P₃₂ values were then calculated using progressively smaller boxes (Figure 2). The input P₃₂ values were compared with the P₃₂ obtained for the different volumes. This modelling exercise aimed to quantify the generation box’s effect on the models’ target intensity.

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 P₃₂ 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)

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)

3.2. Results

Figure 4 summarizes the results for a mean fracture radius size of 20 m and a P₃₂ input of 2 m⁻¹. The results show a difference of up to 30% between the input and target intensity, thus demonstrating that boundary effects can affect the intensity. P₃₂ values sampled at a scale corresponding to the generation box yield the input P₃₂. 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 P₃₂ value increases, showing an asymptotic behaviour. Contrary to what may be expected, the models with fractures generated using their centres as reference present a better agreement with the input P₃₂, and the average curve stabilizes faster than the one for the models in which fractures are generated using the surface points option.

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 P₃₂ and sampled P₃₂ 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 P₃₂. This correction factor can be defined as the ratio between the input intensity and the average P₃₂ 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.

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 P₃₂ of 2 m⁻¹. Then average P₃₂ values were calculated using progressively smaller boxes. To define a target intensity P₃₂ of 2 m⁻¹ 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% of the target intensity within the volume of interest. This reduction of 84% corresponds to the factor in which the average curve of P₃₂ stabilizes (2.00/2.38 = 0.84 of the average sampled P₃₂ presented in Figure 5a). When using this methodology, obtaining the target P₃₂ of 2 m⁻¹ 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.

3.2.2. Required Number of Realizations

To investigate the number of realizations required to obtain a meaningful sampled average P₃₂ in the volume of interest, the cumulative average of the sampled P₃₂ was plotted against the number of realizations (Figure 7). This provides a graphical representation that determines whether the number of realizations converges to a constant sampled P₃₂.

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 realizations are needed to obtain a meaningful average P₃₂. 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 P₃₂. Note that the number of realizations should increase if the plots do not converge.

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 P₃₂ and the P₃₂ 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. Percentage variation between input P₃₂ and sampled P₃₂ 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%

shows the results of this exercise for a mean fracture radius of 25 m.

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 properties are maintained constant, the correction factor is independent of the fracture intensity; 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 relationship 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 described in the step-by-step guide below and illustrated in Figure 8.

• Define a volume of interest, which corresponds to the volume in which we want to define the target P₃₂.
• Generate fractures in a region bigger than the volume of interest. The size of the generation volume depends on the mean fracture size and its standard deviation. A sensitivity analysis could be performed to estimate the size of the generation volume.
• Run an adequate number of realizations with an initial P₃₂ close to the target P₃₂. 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 realizations. If the plots do not converge, then the number of realizations should be increased. Conversely, if the curves converge quickly to the average sampled P₃₂, then it is possible to decrease the number of realizations. This will reduce the total computation time, obtaining similar results.
• Calculate P₃₂ in the volume of interest and other smaller volumes to check that the curve of average sampled P₃₂ has stabilized for the volume of interest. If the curve has not stabilized, it is necessary to increase the size of the generation region (go back to step 2).
• If for the volume of interest, the curve is in the asymptotic section, calculate the average P₃₂ in the volume of interest and the ratio between input P₃₂ and the average P₃₂ in the volume of interest. This radio corresponds to the correction factor.
• Rerun the model using the initial P₃₂ multiplied by the correction factor.
• Check that the P₃₂ in the volume of interest corresponds to the target P₃₂ within an acceptable 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 the acceptable margin of error is recommended. For this paper, a difference of ±1% of the target intensity has been adopted.

4. Variability of Linear Fracture Intensity and Its Relationship with Volumetric Intensity

The relationship between the volumetric intensity (P₃₂) and the linear fracture intensity (P₁₀) was investigated using a series of DFN models of known intensity. The DFN models were then sampled well with synthetics to obtain P₁₀.

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. Input parameter used in the models to investigate the reliability of borehole-derived intensity.

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 (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

. 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 P₃₂ was sampled in the volume of interest and only models with a difference smaller than 1% of the target intensity were considered.

The DFN models were sampled using 10 synthetic wells (Figure 9), and the orientations 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 the locus of the poles of joints parallel to the drill hole and less likely to be observed. Note that 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°) are less likely to be observed in that particular well.

When sampling the DFN models with the synthetic wells, a linear relationship between input P₃₂ and sampled P₁₀ was observed (Figure 10). As expected, greater sampled P₁₀ dispersion is observed for greater values of input P₃₂. The linear relationship between P₃₂ and P₁₀ was already demonstrated using DFN modelling by Dershowitz et al. [10].

When the models are generated using the methodology proposed to mitigate the boundary effects, the relationship between P₃₂ and P₁₀ is independent of the fracture size used (Figure 11). While this statement suggests that P₃₂ 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.

The results show that the mean and standard deviation of the sampled P₁₀ follows a linear relationship. Therefore, estimating the P₁₀ variability for any given P₃₂ using a normal probability function for a given fracture orientation is possible. Since P₁₀ and P₃₂ follow a linear relationship, using probability functions and interpolation, it is then possible to estimate the variability of P₃₂ for any given P₁₀ (Figure 12).

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 intensity [13]. This section presents the impact of changing interval length and the minimum bias angle for the calculated P₃₂. A methodology to calculate P₃₂ 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 intervals.

5.1. Effect of the Minimum Bias Angle Considered to Estimate P₃₂ from P₁₀

Using the DFN model presented in Section 4, the volumetric intensity (P₃₂) was calculated from the linear intensity (P₁₀), using Equation (4). Note that P₃₂ was calculated using the whole well length, and the analyses focused on evaluating the effect of the minimum bias angle α in the estimation of P₃₂. For this, three cases were considered: minimum α of 15°, 5° and 1°.

When limiting the alpha angle to a minimum of 15°, the calculated P₃₂ is underestimated (

Table 4. Ratio between input P₃₂ and calculated P₃₂ from borehole intensity.

Well	Minimum α of 15°		Minimum α of 5°		Minimum α of 1°
	Average Ratio	SD of Ratio	Average Ratio	SD of Ratio	Average Ratio	SD of Ratio
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
Well 6	1.11	0.10	1.04	0.10	1.01	0.11
Well 7	1.11	0.10	1.04	0.10	1.02	0.12
Well 8	1.16	0.12	1.05	0.12	1.02	0.13
Well 9	1.15	0.11	1.05	0.12	1.02	0.13
Well 10	1.13	0.10	1.05	0.11	1.03	0.12
Total	1.15	0.12	1.05	0.12	1.02	0.13

), 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 minimum angle decreases, suggesting that the dispersion of the data increases for low angles. Table 4 shows a good agreement between the input P₃₂ and the calculated P₃₂ 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 P₃₂ when several boreholes or scanlines with different orientations are available, but it does not allow for calculating the spatial variability of intensity along the wells.

5.2. Effect of Interval Length on the Calculated P₃₂ from P₁₀

Because of the tendency to calculate linear intensity using run lengths, smaller angles will result in intervals with high artificial variability in the calculated intensity. Note that the practice of using run lengths as reference intervals has no scientific basis and is largely dictated by empirical field practices, which are difficult to change. It is suggested to use linear intensity values determined for geotechnical domains, which may include several adjacent cores runs.

The effects of varying the interval length on the calculated P₃₂ were investigated using simple DFN models with constant fracture orientations. When the fracture orientation and the sampling length are constant, Equation (4) can be expressed as

$${\widehat{P}}_{32}=\frac{1}{L}{\sum }_1^n\frac{1}{\mathrm{s}\mathrm{i}\mathrm{n} \alpha }=\frac{N}{L \mathrm{s}\mathrm{i}\mathrm{n} \alpha }=\frac{P_{10}}{\mathrm{s}\mathrm{i}\mathrm{n} \alpha }$$

(5)

A volume of interest was defined within a box of 100 m per side, and DFN models were generated using a box size according to the recommendations presented in Section 3, and using the properties presented in

Table 5. Input parameters used in the models to investigate the effects of interval length for the calculated P₃₂.

Property	Comments
Orientation	Constant Plunge of 90°-45°-30°-15°-5°-2-1°
Spatial model	Enhanced Baecher with generation locations at fracture centre
Intensity P32 (m−1)	1-2-4-8
Size: Lognormal Distribution	Radius Mean (Xmean): 2-10-20 (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	10 realizations per model (840 realizations in total from the combinations of intensity, size and α)

. DFN models were sampled with a 100 m long vertical well. Under these conditions, the α angle is equal to the fracture’s input plunge, and the results are independent of the fracture trend.

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 the intervals were used to calculate the actual P₃₂ per interval (sampled P₃₂). The actual P₃₂ per interval was then compared with the P₃₂ calculated using Equation (5), based on the P₁₀ values sampled per Well interval. Note that for this exercise, estimating the P₃₂ correction factor was unnecessary, since the actual P₃₂ in each box was compared with P₃₂ calculated in each interval. For the same reason, the actual and sample P₃₂ values are higher than the input P₃₂ values presented in Table 5.

Figure 13 presents the results of comparing actual P₃₂ and calculated P₃₂ for different interval sizes. For small α angles, the dispersion increases significantly for small size intervals, with many values out of the graph scale. On the other hand, for 50 m and 100 m intervals, the values tend to follow the line that defines the linear fit of the data. Note that 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 a small interval is low. Therefore, many of the small intervals present P₃₂ values of zero. 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 P₃₂ values close to the actual intensity (

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
	2 m	5 m	10 m	20 m	50 m	100 m	2 m	5 m	10 m	20 m	50 m	100 m
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

). This effect was already discussed by Hekmatnejad et al. [21], who found that on a composite scale (interval scale), fluctuations in calculated P₃₂ are large and there may be a significant deviation between the calculated and the actual value of P₃₂. However, on average, the values tend toward the actual P_32.

It can be observed that even when the P₃₂ calculated in small intervals present a high dispersion, their mean tent to represent the actual P₃₂ (Table 6). This interesting result suggests 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 calculated P₃₂ was plotted for each interval length. Figure 14 shows the maximum calculated P₃₂ per interval for input intensities of 1 m⁻¹ and 8 m⁻¹. 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.

5.3. Proposed Methodology to Calculate P₃₂ from Borehole Intensity

Based on the analyses performed, a methodology to calculate P₃₂ from borehole intensity is proposed. This methodology can be easily implemented using code, and the idea behind this is to allow the calculation of P₃₂ 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 interpolation values in the block model.

The following methodology is proposed to calculate P₃₂ from borehole intensity:

• Since P₃₂ is an additive variable, it is possible to discretize a borehole in different overlapping intervals and calculate P₃₂ as the addition of the P₃₂ 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 P₃₂ using the correction proposed by Chilès et al. [13] (Equation (4)), using regular intervals whose length depends on the magnitude of the acute angle (α) between the scanline and the fracture.

  Table 7. Recommended interval length to calculate P₃₂.

  α Angle (°)	Interval Length to Calculate P32	Comments
  15 ≤ α	10 m to 20 m	Use a length similar to the size of the block model cells
  2 ≤ α <15	20 m to 50 m	P32 must be calculated within the same Structural Domain
  1≤ α <2	50 m to total borehole length	It can be extended at the full well length; P32 must be calculated within the same Structural Domain

  presents the recommended length intervals for each angle.
• Add the P₃₂ values calculated using the smallest intervals; this addition corresponds to the total P₃₂.

Note that if only the average intensity is required, it is recommended to calculate P₃₂ using the whole length of the borehole or the length within the structural domain of interest. It is still recommended to limit the minimum α angle to 1°, to avoid extreme high-weighting factors or division by zero.

The proposed methodology was tested using the DFN model and synthetic wells presented in Section 4. The interval lengths used, depending on the acute angle between the well and the fracture, are presented in

Table 8. Intervals used to test the methodology proposed to calculate P₃₂ depending on α angle.

α Angle (°)	Interval Length to Calculate P32	Maximum Weighing Factor per Fracture
15 ≤ α	15 m	3.9
2≤ α <15	30 m	28.7
1≤ α <2	90 m (well total length)	57.3

.

Independent of the borehole direction, a good agreement was observed between the input P₃₂ and the median (50% percentile) of the calculated P₃₂. Figure 15 presents a comparison between the proposed methodology and the P₃₂ values calculated in Well 5 using different interval sizes and minimum α angles, while

Table 9. Comparison between the proposed methodology and the P₃₂ 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
		% Variation Median (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
		% Variation Median (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
		% Variation Median (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
		% Variation Median (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%
Proposed Methodology		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
		% Variation Median (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%

⁽¹⁾ Linear fit factor between P₃₂ and calculated P₃₂ (Input P₃₂ = factor x Calculated P₃₂), note that this factor is equivalent to 1/slope of the linear fit lines presented on Figure 15; ⁽²⁾ Percentage variation between input P₃₂ and median (Percentile 50%) of the calculated P₃₂.

summarizes the result obtained for each one of the 10 wells analyzed.

The methodology herein proposed is the one that provides the best result in terms of a good agreement between input P₃₂ and calculated P₃₂, and at the same time maintaining a relatively low variability. Note that when a minimum angle of 15° is considered the mean P₃₂ values calculated are underestimated. This is especially relevant considering it is common practice to use a minimum angle of 15° [13] to restrict the maximum weighting factor.

An interval of 3 m corresponds to the typical length of a drill run. It is common practice (but not necessarily best practice) to calculate linear intensity on a run basis and then estimate P₃₂ using the same run length. The analyses show that using a length of 3 m, even when restricting the minimum angle to 15 °, produces a high variability on the calculated P₃₂, meaning that the calculated P₃₂ will often be under- or overestimated. Note that, in practice, depending on the drilling operation and rock quality, the actual drilling run lengths are variable and may be significantly smaller than 3 m. This means that if P₃₂ is calculated using those shorter intervals, the P₃₂ variability may be much higher than the one calculated using 3 m intervals. These results agree with the work by Elmo and Stead and Yang et al. [22,23] on the problems of using run lengths to calculate RQD (rock quality designation [24], which is an empirical rock mass quality parameter linked to fracture frequency (i.e., linear intensity).

As another validation method, P₃₂ 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, P₃₂ is calculated as the area of fractures within the grid cell divided by the grid cell volume. Figure 16 presents an example of the P₃₂ values calculated in the grid for one realization using an input P₃₂ of 6 m⁻¹ 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 to capture the spatial variation of P₃₂ and simultaneously obtain a good agreement between the actual P₃₂ (P₃₂ values calculated in the grid) and the P₃₂ calculated along the wells. On the other hand, when P₃₂ 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 presents the result of 100 realizations as cumulative frequency curves of actual P₃₂ in the grid cells and the calculated P₃₂ in wells. The proposed methodology is the one that produces the best results, especially when compared with the common practice of limiting the minimum angle to 15° and calculating P₃₂ using the drill run length (3 m).

6. Summary and Conclusions

A main component of rock mass characterization is the characterization of discontinuities, which play a major role in rock mass’s mechanical and hydraulic properties. This paper focused on volumetric fracture intensity (P₃₂), which remains one of the key properties of fracture characterization for DFN modelling, albeit it is not possible to measure it directly.

The problem of estimating P₃₂ from information observed along boreholes was addressed by generating a series of DFN models and then comparing the volumetric intensity of those models with the volumetric intensity calculated from borehole intensity, using 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 P₃₂. Subsequently, the paper has focused on presenting a new method for calculating P₃₂ 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 P₃₂ calculated in borehole intervals. This is very useful when the interval intensity is used as input for block modelling of P₃₂ 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 P₃₂ 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 fracture orientations.
• The analysis showed that when the boundary effects are mitigated, the relationship between P₁₀ and P₃₂ 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 result should not be construed to suggest that fracture size is not an important parameter; indeed, fracture size remains a required input for generating a DFN model.
• The analysis showed that simulations can quantify the variability and distribution of the volumetric intensity for a given linear intensity. Note that more than a unique P₃₂ value can be obtained for a given P₁₀ value, but with this methodology, it is possible to quantify the probability of obtaining a certain P₃₂ for a given P₁₀.
• The general practice of constraining the minimum bias angle to 15° and calculating P₃₂ using small intervals (drill run lengths) results in an underestimation and a high variability in the calculated P₃₂. When using the methodology proposed by Chilès et al. [13] to calculate P₃₂, the interval’s size greatly impacts the variability of the calculated P₃₂. Using small intervals will artificially increase the variability of the calculated P₃₂, even when the average P₃₂ tends to be close to the actual P₃₂. On the other hand, when limiting the minimum angle between the fracture and the borehole to a minimum of 15°, the mean calculate P₃₂ values are underestimated.
• The proposed methodology to calculate P₃₂ using variable lengths, depending on the angle between the fractures and the borehole, allows the capture of the spatial variation in intensity and at the same time avoid increasing or decreasing artificially the intensity of the intervals. This is very useful when the interval intensity is used as input for interpolating P₃₂ values in block models.
• When only average volumetric intensity is required, it is recommended to calculate P₃₂ using the total length of the borehole or the length within the structural domain of interest. In this case, limiting the minimum α angle to 1° is recommended to avoid extremely high-weighting factors or division by zero.

Limitations and Assumptions

The research described previously is subject to the following limitations and assumptions:

• The main limitation of this study is that the analyses presented are purely simulated and have not yet been applied to real data. The concept of validating modelling results with field data does not apply to P₃₂ values since it is impossible to directly compare calculated P₃₂ with the actual P₃₂ of the rock mass (which remains unknown).
• All fractures have been to be assumed planar and extremely thin. Although fractures in nature can be curved and present thickness, it is common practice to assume that discontinuities 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, some investigators indicate that fractures are in reality not equidimensional [26].
• It was assumed in the models that fractures are randomly and independently distributed in space. In reality, fracture intensity may be dependent on geological features; for example, a common geological observation is that the distance to major faults controls fracture intensity [27].
• The analysis performed assumed that all fractures intersected by the synthetic wells were identified. In practice, there are biases when gathering fracture intensity from boreholes. For example, besides the orientation bias, it is common that during borehole logging, not all natural discontinuities are counted, and some artificial fractures (e.g., mechanical breaks) are counted as natural fractures. Even if the fracture intensity is calculated using televiewer data, the quality of the data would depend on the quality of the images and the interpretation performed for each discontinuity.
• The methodology proposed to minimize boundary effects may be impracticable when the fractures are too big. For big fractures, it would be necessary to generate an extremely big generation box, that may require extensive computation time.