Estimating the Effect of Fractal Dimension on Representative Elementary Volume of Randomly Distributed Rock Fracture Networks

The effect of fractal dimension (Df ) on the determination of representative elementary volume (REV) was investigated through numerical experimentations, in which a new method was adopted to extract submodels that have different length-width ratios from original discrete facture networks (DFNs). Fluid flow in 1610 DFNs with different geometric characteristics of fractures and length-width ratios was simulated, and the equivalent permeability was calculated. The results show that the average equivalent permeability (KREV) at the REV size for DFNs increases with the increase in Df . The KREV shows a downward trend with increasing length-width ratio of the submodel. A strong exponent functional relationship is found between the REV size and Df . The REV size decreases with increasing Df . With the increment of the length-width ratio of submodels, the REV size shows a decreasing trend. The effects of length-width ratio and Df on the REV size can be negligible when Df ≥ 1 5, but are significant when Df < 1 5.


Introduction
Fractures play a dominating role in the mechanical and hydraulic properties of rock masses and are sources of discontinuity, anisotropy, and heterogeneity.Therefore, to analyze the performance of structure characteristics in fractured rock masses, it is important to accurately select an appropriate model volume to determine the relevant rock properties.In recent years, numerical simulation techniques such as the Monte Carlo method have been developed for modeling fluid flow in models containing complex fracture systems [1][2][3][4][5].It is definite that the permeability of a fractured rock mass changes significantly depending on the size of the model [6].Beyond a certain sample size (area in two dimensions and/or volume in three dimensions), the average permeability tends to attain a critical value and the variation in permeability will be very small and can be negligible.In such a situation, this size can be termed as REV size for a fractured rock mass with respect to fluid flow behavior [1,6].
The previous works have documented the effectiveness of fracture geometry parameters and average strength in determining REV [1,7,8].The REV size decreases with increasing joint density and joint size [1].The permeability change at low stress levels is more sensitive to model size than at high stress levels due to the nonlinear fracture normal stressdisplacement relation [2,7].The determination of REV becomes more difficult for fracture network models subjected to a lower stress environment.The application of the composite element method can greatly facilitate the preprocess and enable a large number of stochastic tests for the fractured rock samples [9].The existence of REV is illustrated by changing the sizes and orientations of the samples.The REV is used to overcome the obstacle that small models contain only a few joints while the large models contain hundreds of joints that lead to impractical computation run time [10].The REV designation embeds a sufficient number of joints in REV-size models.These innovations can reduce computation time by two orders of magnitude from hundreds of hours to a few hours.
The studies of REV in three dimensions have previously been performed.A new method that gives through thorough consideration to the fracture features of rock masses was proposed based on models of 3D fracture networks [11].Subsequently, the REV size was determined by volumetric joint count calculation.The specimens, which are generated by collecting structure data, were introduced into a 3D particle flow code (PFC3D) to create synthetic rock mass (SRM) samples.A series of T-tests and F-tests were carried out to determine the REV size, in which the T-test was used to assess the difference of samples and the F-test was used to describe the calculated variance [12].The basic assumptions were that the rock matrix is impermeable and linearly elastic and that the fluid flows only in fractures [7].
Fractal characterization of complex fracture networks has been proven to be effective through statistical analysis of natural geological rock masses [13,14].The index of fracture density, which is fractal and often scale-invariant, can be predicted through fractal geometry [15].The analytical expression for gas permeability is derived based on dual-porosity media [16].In the dual porosity, the original channel diameter of embedded fractal-like tree networks follows a fractal distribution.A percolation term (ρ ′ − ρ c ′ ) was obtained using a series of geometric properties of fracture networks and had a high correlation coefficient between the actual and the predicted equivalent fracture network permeability.
However, most of the previous studies on REV are focused on square submodels in both 2D and 3D.Few of the studies consider the relationship between fractal dimension and REV.A method is presented to compute the equivalent permeability of submodels extracted from original discrete fracture networks (DFNs) that are generated using the Monte Carlo method.Finally, the relationship between fractal dimension and REV size is established, considering different model length-width ratios.

Theory
Fractures in DFNs are typically treated as parallel plates, in which the flow rate is proportional to the cube of aperture, as follows [17]: where Q is the flow rate, e is the hydraulic aperture of a fracture, μ is the dynamic viscosity, g is the gravitational acceleration, Δh is the hydraulic head difference, L f is the length of a fracture, and W is the aperture of a fracture.
In the deep fractured rock masses, the fluid flow is commonly in the linear regime and obeys the cubic law, which is a kind of Darcy's flow [18][19][20][21].Therefore, we just calculated the laminar flow and did not consider turbulent flow.The permeability of the rock matrix can be negligible when comparing to the permeability of fractures in tight rock masses, i.e., granite and basalt.The fractal dimension D f can be calculated using gauge method and grid method.A fractal permeability model for bi-dispersed porous media is developed based on the fractal characteristics of pores in the media, and a series of functions are derived to establish the fractal permeability model [22].Besides, many studies have successfully proven that the fracture length distribution follows the fractal scaling law ( [23][24][25]; Miao et al. 2015).The total number of fractures, N t , can be calculated as follows: where N t is the cumulative number of fractures, l min is the minimum fracture length, and l max is the maximum fracture length.In addition, l min ≪ l max is the necessary condition for the fracture length distribution to follow the fractal scaling law.l min /l max ≪ 10 −3 is used in the present study.The length of the ith fracture can be calculated as: where l i is the length of the ith fracture, i = 1, 2, 3, … , N t , N t is the total number which has been obtained from (2), and R i is a uniformly distributed random number that varies from 0 to 1. Based on a series of assumptions and derivations, the formula of new total number of fractures corresponding to any side length L n of DFNs can be expressed as [24,25]: After updating the total number of fractures from N t to N t ′ , the fractal length distribution of fractures in any L n of a DFN model can be generated using (3).Thus, a series of stochastic DFNs are generated using the Monte Carlo method and fluid flow is modeled by solving the cubic law (see (1)).

Generation of DFNs
3.1.Fracture System Description.Stochastic models provide validate ways for representing fracture networks based on available structural data.It was noted that fracture system models can be treated as an entity to represent rock mass geometry [26].They further described the Orthogonal, Baecher, Veneziano, Dershowitz, and Mosaic Tessellation models.The paper [27] suggested an incrementally linear elastic, orthotropic constitutive model to represent the equivalent continuum prefailure mechanical behavior of the jointed rock masses.
In general, two fundamental approaches are suggested for modeling fluid flow through fractured rock masses: [28] equivalent continuum flow model [17,[29][30][31][32] and [13] discrete fracture flow model ( [20]; Oda 1985; [28,33]).The first approach assumes that the combined 2 Geofluids hydraulic effect of fractures and rock matrix can be represented by an equivalent continuum model.The second approach treats fractures as separate elements having significantly higher hydraulic conductivity compared to that of the rock matrix.The former model is used to simulate fractures of small size and having a great quantity, while the latter is applicable for large-scale fractures [9].More reasonable models such as continuum-discrete coupling model are also explored.This approach takes the permeability of fractures and capillary in the matrix into account.This model is close to reality, but the workload of numerical simulation is commonly unavailable.In most cases, fracture locations are stochastically distributed and fracture length is specified directly or indirectly [1,12,18,24,34].Many models can accommodate a large quantity of fractures that are intersected and/or terminated in rock masses.A number of fractures can be located in the same plane in 2D.Some sophisticated models can exhibit not only geometric characteristics but also geological structures [35].
The previous studies have shown that the permeability predicted using the 2D model is slightly less than that predicted using the 3D model, i.e., less than one order of magnitude [36,37].Besides, although the 3D model can characterize the spatial distributions of geometric parameters of fractures, these parameters are difficult to obtain because the fractures are buried in rock masses and are nonvisual.In contrast, the 2D model can use the data from outcrops to calculate parameters such as fracture length, aperture, and orientation.Therefore, the 2D model is currently widely accepted by engineers and hydrologists.

Discrete
Fracture Network Generation.The Monte Carlo simulation technique is used for the generation of the DFNs.
Table 1: Two sides of rectangular submodels compared to the length of square submodels with the same area when D f = 1.6.

Geofluids
The theory and technique of the issue have been well studied and presented in previous studies [27,[38][39][40].The process of generation includes the following steps: (1) Several parameters, such as l min , l max , D f , and L n , should be determined before DFNs are generated.l min and l max are set to be 0.5 and 500 m, respectively, which satisfies the necessary condition for the size distribution of fracture length to follow the fractal       (2) The value of N t can be calculated using (1).The length of the ith fracture can be calculated using (3) after generating N t random numbers.N t ′ can be obtained by calculating (4).After that, the lengths of all fracture with a number of N t ′ can be generated using (3). ) 6 Geofluids area.Table 1 shows the two sides of rectangular submodels compared to the length of square submodels with the same area when D f = 1.6.
(5) Constant hydraulic head boundary conditions are applied, as shown in Figure 2, to numerically calculate the permeability (K p ) along the direction of the hydraulic gradient when fluid flow is in the steady state.The upper and bottom sides are impermeable, and the horizontal flow from the left side to the right side is set for the calculation.The hydraulic pressure gradient is 10 kPa/m for each DFN.
When the geometric parameters such as fracture length, aperture, and orientation are available, our model can simulate the fluid flow in real situations.However, characterizing these parameters is not a simple work, especially for the 3D model.Generally, for simplification, some functions are used to generate fracture parameters.For example, fracture length distribution follows a power law function [41,42]; fracture    7 Geofluids aperture is correlated with fracture length [18,19]; and fracture orientation follows the Fisher distribution [2,8].By doing so, the calculated results (i.e., permeability) are still close to the in situ results.

Discrete Fracture Network Validation.
A fundamental concern is whether the generated DFNs are representative of observed field conditions.This is often difficult to quantify due to the limited available field data.This can be overcome by comparing the statistical information of generated fracture length distribution with those reported in literature.The feasibility and effectiveness of fractal length distribution of fractures are illustrated.Figure 3 shows the results of the relationship between the length and numbers of fractures.It is obvious that each set of length-number relationship corresponding to each D f follows a power law function that can be expressed by: n l, 0 5 = αl −a , 5 where α is the proportional coefficient, a is the power law exponent, and n l, 0 5 represents the fracture number with lengths in the range of l − 0 5, l + 0 5 .L n (m) When D f varies from 1.3 to 1.7, a is in the range of 0 96, 1 29 , which fits well with the values reported in previous studies as shown in Table 2.This verifies the validity of ( 2) and (3), and indirectly verifies the validity of the generated DFNs.In practice, the choice of the model depends on how it can be correlated to the available field data and to be engineering needs of the project.Therefore, the validity of the generated models still needs further verification using in situ data collected from different scales and different rock types.

Determinations of Equivalent Permeability and REV Size.
Using the abovementioned methods, the flow rates are obtained under a fixed hydraulic pressure gradient for models that have different D f .Next, the equivalent permeability is back-calculated using Darcy's law, as follows: ) ) where μ is the fluid viscosity, Q is the total flow rate, A is the cross section of percolation, and ∇P is the hydraulic pressure gradient that equals 10 kPa/m.For each DFN, the average width of the cross section is assumed to be 1.0 m because the DFN model is two-dimensional and thus A equals the side length of the inlet boundary.
Based on the extracted DFNs with different length-width ratios as shown in Section 3.2, the variations in K f and corresponding root mean square (RMS) are plotted in Figures 4-8.For each D f and each length-width ratio, 10 DFNs are generated using 10 sets of random numbers.Then, K f and its mean value of the 10 DFNs are calculated.K f tends to stabilize with increasing model size.Before the curve stabilizes, the permeability of the models depends on the connectivity of the fractures.The REV can be defined theoretically as the size beyond which the permeability varies in a sufficiently small range.In practice, the corresponding model scale, in which L n increases to a critical value and the change in K f is sufficiently slight, can be regarded as the REV size.V K and RMS decrease with the increment of L n .Here, V K is defined as the difference in equivalent permeability between the maximum and minimum values.The RMS is determined by the ratio of standard deviation of simulation results to the mean value of the maximum and minimum equivalent permeability.
The RMS is introduced to quantify whether two datasets agree well with each other [3], defined as: where K max is the maximum equivalent permeability, K min is the minimum equivalent permeability, K avg is the average equivalent permeability, and K sim is the simulation result.
Taking the case of D f = 1 6 in Figure 6(a) for an example, K p ranges from 8.5E − 14 to 5.5E − 13 m 2 when L n = 0 5 m, while for L n = 10 m, K p ranges from 3.5E − 13 to 3.1E − 13 m 2 .Besides, V K decreases from 4.95E − 13 to 4E − 14 and RMS decreases from 0.443 to 0.032, indicating that the connectivity of the fractures tends to be homogeneous.However, when the length-width ratio of submodels changes, the number of connected fractures will change, and so does the K p of submodels with the same area.
In this study, RMS = 0 1 is utilized as the threshold value, and the corresponding side length of DFN (L REV ) is treated as the REV size.The variations in permeability change significantly when RMS > 0 1 and therefore RMS = 0 1 is used here rather than other values such as 0.2 [3].The results show that L REV gradually decreases with the increment of D f from 1.3 to 1.7 as shown in Figure 9.When D f is small (i.e., D f = 1 3), the model length-width ratio can influence the magnitude of L REV significantly (i.e., from 28.12 to 35.16 m).When D f is large (i.e., D f = 1 7), the influence of the model length-width ratio on L REV can be negligible (i.e., from 1.32 to 1.36 m).Therefore, when calculating the REV of a fractured rock mass, the shape of the submodels should be considered when D f is small (i.e., D f ≤ 14).By fitting these datasets, an exponential function is proposed between L REV and D f , written as: Thus, when D f of a fractured rock mass is calculated using the geological survey data from the outcrops, the REV size can be approximately estimated using (8).Besides, the calculated results are compared with those in Liu et al. [25], in which the REV size is analyzed using rectangular submodels with L x L y = 1 1, showing an agreement with each other.When D f is small (i.e., 1.3), the difference L REV between the predicted results and those in Liu et al. [25] is relatively larger (i.e., 8.51 m).However, with the increment of D f from 1.4 to 1.6, the difference is very small (from 2.39 m to 0.43 m).This is acceptable because when D f is small, the randomness of fracture distribution is significant and therefore L REV changes in a large extent.With the increment of D f (i.e., D f > 1 4), the effect of randomness of fracture distribution decreases, and as a result, L REV varies slightly.This again verifies the validity of (8).

Effects of D f and Model
Length-Width Ratio on K REV .The variation of K REV is plotted against the length-width ratio in Figure 10, where K REV is the equivalent permeability at the REV size.As the model length-width ratio of submodels increases, K REV gradually decreases.Besides, the rate of permeability reduction also decreases with increasing model length-width ratio.The effect of the model lengthwidth ratio on K REV when D f = 1 3 is much smaller than that when D f = 1 7.The variation in K REV is plotted against D f in Figure 11.The influence of D f is determinate, as what have been observed for the influence of the model lengthwidth ratio.The reason may be that the DFNs with larger D f have better connectivity and stronger conductivity/permeability. K REV increases with the increase in D f , while the variations in the K REV increment for D f ≥ 1 6 are much larger than those for D f < 1 6.  10 Geofluids

Conclusions
In this study, we calculated the equivalent permeability of rock fracture networks that have different models lengthwidth ratios.This study, using a numerical method, shows another possible way to study the permeability characteristics of fractured rock masses.Submodels of different model length-width ratios, instead of square submodels, are extracted from original DFNs.Since this study requires a large quantity of numerical computations, an effective and reliable algorithm is essential.A DEM (discrete element method) code based on a 2D open-source software (OSS) UDEC is written to generate models and calculate flow rates.Generated DFNs with different model length-width ratios are used to simulate Darcy flow in rock fractures.Through systematic and large quantity numerical computations, the permeability characteristics of DFNs and the existence of REV are identified.The agreement between the predicted REV size and those reported in literatures verifies the validity of the proposed method.In conclusion, a new method using fractal dimension of the fracture networks for calculating REV is proposed.The REV size of the fractured rock mass decreases with the increase in D f according to the exponential relationship between the REV size and D f .The effect of the length-width ratio of submodels on the REV size decreases with increasing D f for D f < 1 5, but can be negligible when D f ≥ 1 5. K REV decreases with the increments of the model length-width ratio of submodels and/ or D f .The influence of the model length-width ratio on K REV when D f is small (i.e., D f = 1 3) much smaller than that when D f is large (i.e., D f = 1 7).Besides, the variations in the increment of K REV when D f ≥ 1 6 are much larger than those when D f < 1 6.
The above conclusions are obtained based on the 2D DFNs using the assumption that the fluid flow obeys Darcy's law.Further investigations considering non-Darcy flow in 3D DFNs with variable fracture apertures will be performed in order to enrich the results and strengthen the conclusions.However, this requires a relatively high computation capacity.In addition, laboratory and field tests as well as engineering practices can be used to verify the feasibility and effectiveness of the proposed method.

Figure 1 :
Figure 1: Schematic view of the generation of submodels.

Figure 3 :
Figure 3: Correlation between fracture number and fracture length.
scaling law.The fractures with lengths that are smaller than 0.5 m contribute negligibly to the flow rate with respect to long fractures.The location of fractures can be determined by setting the center point, orientation, and fracture length.The fracture orientation and center point distribution are assumed to be uniformly and randomly distributed in order to focus on the effects of fractal length distribution.Original square DFNs having side lengths (L n ) of 50m, 25 m, 15 m, 10 m, 8 m, and 4 m are generated when D f = 1 3, 1.4, 1.5, 1.6, and 1.7, respectively.L n decreases with the increase in D f .If a large L n (i.e., 50 m) is adopted when D f is large (i.e., D f = 1 7), there would be many fractures generated, which may cost a longer time for solving fluid flow, comparing with that with a smaller D f (i.e., D f = 1 3).

Figure 6 :
Figure 6: Variations in K p and RMS with different L n when D f = 1 5.

Figure 7 :
Figure 7: Variations in K p and RMS with different L n when D f = 1 6.

Figure 8 :
Figure 8: Variations in K p and RMS with different L n when D f = 1 7.

Figure 9 :
Figure 9: Relationship between D f and L REV for different L x :L y .

Figure 11 :
Figure 11: Relationship between L REV and D f .

Figure 10 :
Figure 10: Relationship between K REV and L x L y for different D f .

Table 2 :
Comparison of the values of a and achievements in previous studies.