#### Abstract

Determination of the representative elementary volume (REV) of fractured rock masses based on equivalent permeability () is significantly dependent on the geometric characteristics of fractures. In this work, a series of numerical simulations were performed to analyze the relationship between geometric characteristics of fractures and the REV size, in which fracture length follows a power-law distribution. A method to evaluate the of a three-dimensional (3D) discrete fracture network (DFN) by extracting the equivalent pipe network (EPN) model from the DFN model was utilized and verified. The results show that of the 3D DFN model has an exponential relationship with the power exponent () of fracture length distribution and the evaluation of agrees well with that reported in previous studies, confirming the reliability of the EPN model for calculating seepage properties of complex 3D DFN models. When the side length of submodels () is small, the varies significantly due to the influence of random number seeds used to generate fracture length, location, and orientation. The holds a constant value after exceeds some specific value. The critical model scale is determined as the REV size, and the corresponding volume of the 3D DFN model is represented by . The varies within a narrow range when . When , the rapidly increases to more than 3.4 times than that when . The fluid flow becomes more inhomogeneous due to the small nonpersistent fractures that dominate the preferential flow paths when exceeds a certain value (i.e., 4.5). The at the REV size decreases exponentially with the increment of . This tendency can be explained by the decrease of the average intersection length () with the increment of , which is a geometric parameter for reflecting the connectivity of the fracture network.

#### 1. Introduction

The hydraulic properties of fracture networks which constitute the main flow paths are of significant importance for fluid flow and transport in fractured rock masses with negligibly low permeability of matrix [1–6]. Natural fracture networks usually show a strong hydraulic variability which comes from the overall statistical geometry characteristics [7–13]. The discrete fracture network (DFN) models are conceived for modeling the spatial structures and geometric features of naturally fractured rocks quite well [7–10, 14–22]. Therefore, DFN modeling is applicable for investigating the effect of geometric characteristics on the hydraulic properties in fractured rock masses.

Many efforts have been denoted to investigate the effect of statistical geometry characteristics on the hydraulic properties in fractured rock masses using DFNs. De Dreuzy et al. [7] developed a two-dimensional (2D) fractured model and gave the dependence of the equivalent permeability of fracture network on the power-law exponents of the fracture length. Min et al. [14] determined the equivalent permeability tensors of fractured rocks using the stochastic representative elementary volume (REV) concept that utilizes numerous realizations of DFN models. And the methodology is applied to the field data in which a power law characterizes the fracture trace lengths. The anisotropy of permeability of studied fracture rock masses was not high, and the existence of REV was confirmed. Lei et al. [15] compared the hydraulic behavior of the analogue fracture network (AFN) and DFN when geomechanical changes are present. The results revealed that DFNs could accurately reflect an AFN under stable mechanical conditions; such representations may be unreliable when mechanical changes occur. However, the statistically modeled 2D DFN model is short for capturing the geometric properties (i.e., orientation towards the out-of-plane direction) of 3D fracture networks in fractured rock masses [22]. Some studies have shown that directly estimating 3D fracture network permeability using 2D fracture network permeability causes considerable inaccuracies in both direction and magnitude [10, 23].

Therefore, 3D DFN modeling is applicable for predicting the hydraulic behaviors of the natural fracture network. Bang et al. [9] developed a 3D DFN model to estimate the REV of granitic rock mass at Korea Atomic Energy Research concerning the hydraulic properties. Huang et al. [24] investigated the effect of power exponent on the equivalent permeability of the 3D DFN model by solving the Reynolds equation using the Galerkin method, which is consisted of fractures following power-law distribution. Xu et al. [19] assessed the permeability of the 3D DFN model built with field measurement by calculating the fluid flow through the pipe network equaled from the DFN model. And the size of REV of the studied fractured rock mass was determined. The previous studies clearly show that the power law is a well-known statistical distribution of fracture length. The equivalent permeability can determine the REV size of the DFN model. However, the effects of power exponent of fracture length distribution () on the REV size of 3D DFN and the equivalent permeability at the REV size have not been considered.

In this study, a series of 3D DFN models with the same fracture density () and different were established to calculate the of fracture networks, in which the fracture length distribution follows the power law. In DFN models, the fracture discs are equivalent to the 3D equivalent pipe network (EPN) model to raise the computational efficiency. The validity of the 3D EPN model was verified by comparing the simulation results of the EPN model with the simulation results of DFN models reported in the literature. The REV size is determined by calculating the equivalent permeability of 3D DFN submodels extracted from 3D original DFN models. The influence of the on REV size and of fracture networks at the REV size was systematically investigated. The flow chart for describing this study is shown in Figure 1.

#### 2. Flow Modeling and Verification

Efficient 3D DFN models were generated with the Monte Carlo method to mimic fracture networks based on geometric characteristics such as fracture orientation and length, which follow the statistical geometry characteristics of fracture network. DFN model was equivalent to the EPN model for calculating fluid flow through fracture network. The flow behavior of the EPN model can be characterized by the finite difference method (FDM). The validity of the EPN model was verified by comparing the simulation result with that reported in the literature.

##### 2.1. Generation of 3D DFN Model and Corresponding EPN Model

Natural fractures come at a wide variety of scales, from microscopic to regional [25, 26]. The 3D DFN method has been frequently used to simulate natural fractures, including properties like fracture length and orientation that follow a probability distribution [10, 19]. Circles, ellipses, and polygons with varied aspect ratios are commonly used to model fracture shapes [27, 28]. All fractures in this study are assumed to be circular disks shaped in the impermeable matrix, whose center locations are uniformly and randomly distributed in the interest region. The fracture orientations follow the fisher distribution. Aperture heterogeneity is not considered because the fracture apertures are all fixed to 1 mm. Based on these assumptions, the fracture length distribution and fracture density are the two most critical characters that determine the hydraulic properties of 3D DFN models. Fracture density () of the 3D DFN model is a classical parameter defined as the total area of fracture surface per unit volume [29]. A power-law function is used to model the length distribution of fractures, written as where is the coefficient of proportionality, is the fracture length which distributed between the lower () and upper () bounds, and is the power exponent of fracture length distribution varying between 2.0 and 4.5 [30].

The fracture network is generated by the Monte Carlo method. The fracture network in Figure 2(a), for instance, is generated based on the statistical geometry parameters as listed in Table 1, whose power exponent of length distribution is 4.5. The fracture network comprises four groups of fractures whose fracture lengths follow power-law distribution, as shown in Figure 2(b).

**(a) Three-dimensional DFN model**

**(b) Statistical histogram of the fracture length**

**(c) Intersections**

**(d) Diagram of equivalent pipe model**

**(e) Backbone**

**(f)**Comparison of results of with differentThe 3D equivalent pipe network (EPN) model can be extracted from the 3D DFN model to effectively estimate the fracture network’s connectivity and meet engineering needs. The EPN model is constructed by connecting the center of fracture disk and the center of intersection lines (Figure 2(c)) with other fractures as shown in Figure 2(d) based on an open-source software ADFNE [28]. Figure 2(e) shows the final pathway for the 3D DFN model in Figure 2(a) after eliminating the isolated and practically isolated pipes.

##### 2.2. Verification of the EPN Model

For modeling the fluid flow in the EPN model, a handful and simple technique is established by applying the finite difference method [31], whose basics are pretty simple, well studied, and documented. All flow rates from income pipes and outcome pipes have to fully satisfy each inner node’s mass conservation law. A system of matrix () is set up for every inner node as follows: where is the pressure head of node , is the conductance of the pathway between node and node , and is the pressure head of node . can be expressed as where is the gravity acceleration, is the aperture which is equal to the aperture of fracture in this study, is third of fracture dimension, is the kinematic viscosity of the fluid, and is the length of pipe. The expressions for and are based on the following assumptions: the rock matrix is incompressible and impermeable; the fluid is incompressible; the surface of circle disc fracture is parallel and smooth. That is, the fluid flow obeys the laminar flow of incompressible fluid.

When the equation system () is fully prepared, a simple linear arithmetic technology () helps to find the solution for pressure head values of all inner nodes. After solving the system of equations, a smooth pressure spectrum for pipes is presented as Figure 2(e), in which the inlet pressure head is 1 m (dark area) and the outlet pressure head is 0 m (dark blue), respectively. The total flow rate is the sum of the flow rate of all pipes on the inlet boundary or outlet boundary. The flow rate on these two boundaries is equal. The flow rate through a pipe can be calculated using where is flow rate and is pressure head difference between two nodes. After the flow rate is obtained, the can be back-calculated according to Darcy’s law, as follows: where is the cross-sectional area. Alghalandis [28] described more details about the solving process.

To verify the validity of the EPN model, the equivalent permeability coefficients of six 3D DFN models were calculated by the EPN model. The six DFN models are generated with the statistical geometry parameters as tabulated in Table 1, in which the power exponent increases from 2.0 to 4.5 with an interval of 0.5. As shown in Figure 2(f), the equivalent permeability of fracture network decreases exponentially with the increment of , and the simulation results agree well with the values reported in the literature. The differences between simulated and reported results may be caused by the uncertainty of the number of fractures of each group. Therefore, the EPN model is reliable for simulating the flow characteristics of the 3D DFN model.

#### 3. Flow Simulations in 3D Fracture Networks

A series of 3D DFN models were generated to analyze the effect of power exponent of fracture length distribution on the representative elementary volume (REV) of fractured rock masses and the equivalent permeability at the REV size, in which fracture length follows a power-law distribution. The 3D DFN models are equivalent to the EPN models to avoid a large number of meshes performed on each fracture that may lead to time-consuming work.

##### 3.1. Determination of the Number of Fractures in DFN Models

For each fracture length distribution, a series of 3D DFN models with different numbers of fractures () ranging from 6000 to 15000 are generated to evaluate the relationship between the number of fractures and the fracture density () of the 3D DFN model. The expression of can be written as
where is the fracture area in the sampling volume and is the sampling volume. Table 2 shows the geometric parameters used for generating the 3D DFN models. The side length of the 3D DFN models is fixed to 100 m. All fractures are assumed to be parallel-plane circle disks, in which the diameter of the disc denotes the fracture length. The fracture lengths range from 5 m to 20 m with power-law exponents of 2.0, 2.2, 3.0, 3.5, 4.0, and 4.5. The location of center points of fracture discs is assumed to be uniformly and randomly distributed. Figure 3(a) shows the relationship between and for different . The increases with the increment of in a linear way for different . The parameters of fitting curves as shown in Figure 3 are listed in Table 3. And the variation of parameter with is demonstrated in Figure 3(b). Figure 3(b) shows that the increasing rate of () decreases with the increasing . Note that a smaller would lead to more larger fractures and the variation of is more sensitive to a smaller . The for different when m^{2}/m^{3} is determined by the fitting curve in Figure 3(a) and demonstrated in Figure 3(c). The with m^{2}/m^{3} increases (from 6351 to 12800) with the increasing .

**(a)**Variation of with number of fractures

**(b)**Variation of with

**(c)**Number of fractures for different a withThe average length of fracture intersections () of 3D DFN models with different and same changes in relatively small range as shown in Figure 4(a). Figure 4(b) shows that the average value of decreases with the increment of following an exponential relationship. The equivalent permeability of the 3D models for different decreases exponentially with increasing as shown in Figure 5; the parameters of the fitting lines are tabulated in Table 4. This is because the proportion of small fractures increases with the increment of . The large fractures are more likely to contribute continuous pathways for the fluid flow in the fracture network. The connectivity of DFN models with smaller is better than that with larger when keeps the same. In addition, the of DFN models with larger must be greater than that with smaller when keeps the same. Note that the scale effect, which would influence the magnitude of permeability, is not considered. It is uncertain whether the model size reaches REV size. The influence of on *K* at the REV size will be investigated in the next section.

**(a)**Variation of with number of fractures

**(b)**Variation of the average value of with##### 3.2. Generation and Extraction of 3D DFN Submodel with Different

A series of original 3D DFN models with different were generated as shown in Figure 6. For each , ten original 3D DFN models with the same geometric parameters and different randomness were developed to determine the REV size of fractured rock masses. Figure 7 shows the examples of DFN models generated using ten sets of random number seeds. The statistical analysis indicates that the decreases with increasing as shown in Figure 8. Geometric characteristics in this case cannot evaluate the connectivity of the 3D DFN model directly. Therefore, the permeability of the 3D model would be calculated for predicting the influence of on the REV size of fracture networks.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)****(a) Set 1**

**(b) Set 2**

**(c) Set 3**

**(d) Set 4**

**(e) Set 5**

**(f) Set 6**

**(g) Set 7**

**(h) Set 8**

**(i) Set 9**

**(j) Set 10**

Cube submodels are extracted with the side length ranging from 10 m to 90 m with an interval of 10 m from the original 3D models. For all submodels in this investigation, the sampling domain is chosen conservatively, keeping a distance of at least 10% of the model edge length from the original model boundaries. Other studies refer to similar techniques as “flow jacketing” or “oversampling” and claim that this method decreases boundary effects and better reflects the sampling domain’s connectivity [10, 32, 33]. Therefore, the range of side length of the extracted submodel is chosen as 10-90 m when the side length of the original cube DFN model is 100 m. Figure 9 shows the examples of extracting submodels from the original 3D model. A hydraulic water head () difference of 1 m is applied at the opposite sides of cube submodels as shown in Figure 10(b). The left side is inlet with m and the right side is outside with a free boundary ( m), while the other boundaries are impermeable. The left side is the plane for m and the right side is the plane for is equal to the side length of the submodel. As described in Section 2, the flow through the 3D DFN model is calculated by the EPN model to avoid too much meshes performed on each fracture that may lead to time-consuming work. Figure 10(c) shows the water head distribution of the 3D submodel (Figure 10(a)) under the constant hydraulic water head difference ( m).

**(a)**Original DFN model ( m)

**(b)**Submodel ( m)

**(c)**Submodel ( m)

**(d)**Submodel ( m)**(a) A three-dimensional DFN submodel**

**(b) Flow direction (left to right)**

**(c) Water head distribution**

#### 4. Results and Analysis

##### 4.1. Scale Effect

For all cases, a total of 540 cube submodels were extracted from 60 original 3D DFN models, and the permeability of all submodels was calculated. Figure 11 shows the calculated results from ten realizations as the scattered data points of equivalent permeability, along with the bounding curves of their maximum (), minimum (), mean (), the difference between the maximum and minimum () values (Figures 11(a), 11(c), 11(e), 11(g), 11(i), and 11(k)), and corresponding root mean square (RMS) (Figures 11(b), 11(d), 11(f), 11(h), 11(j), and 11(l)). There are ten equivalent permeability coefficients for each side length of the submodel (). The equivalent permeability coefficients of ten submodels with same side length were extracted from ten original 3D DFN models with same geometric parameters and different random numbers seeds. The RMS is used to determine if two datasets are properly aligned. Since the fracture density is relatively small ( m^{2}/m^{3}), is utilized as the threshold value, which is slightly larger than 0.2 reported in the literature [34], and the corresponding volume of DFN is represented with .

**(a)**Relationship between and with

**(b)**Relationship between RMS and with

**(c)**Relationship between and with

**(d)**Relationship between RMS and with

**(e)**Relationship between and with

**(f)**Relationship between RMS and with

**(g)**Relationship between and with

**(h)**Relationship between RMS and with

**(i)**Relationship between and with

**(j)**Relationship between RMS and with

**(k)**Relationship between and with

**(l)**Relationship between RMS and withFigure 11 presents the variations in of cube submodels with different and corresponding RMS. The fluctuates dramatically because of the impacts of random number seeds used to create fracture location, orientation, and length when is small. When crosses specific thresholds, remains constant regardless , allowing the model scale to be used as the REV size and the to be determined. The equals to 1.092E-12 m^{2}, 8.1E-13 m^{2}, 5.18E-13 m^{2}, 4.7E-13 m^{2}, 2.54E-13 m^{2}, and 1.78E-13 m^{2} when is reached at the approximated size of 1.96E+05 m^{3}, 1.88E+05 m^{3}, 1.19E+05 m^{3}, 1.84E+05 m^{3}, 1.96E+05 m^{3}, and 6.70E+05 m^{3} for the cases in Figures 11(a), 11(c), 11(e), 11(g), 11(i), and 11(k), respectively. The becomes less distributed with increasing , with a noticeable trend to attain a REV size. Note that the of DFN models decreases rapidly with the increment of . Since the DFN models are generated based on stochastically assigned fracture lengths that follow a power-law distribution, the proportion of larger fractures decreases due to the increment of . The framework is gradually predominantly structured by smaller fractures. The RMS continuously decreases with increasing , especially before the archives the side length of the REV size as shown in Figures 11(b), 11(d), 11(f), 11(h), 11(j), and 11(l).

The results show an obvious scale dependency of the permeability in the fracture network. With the increment of , the permeability of DFN submodels gradually varies in a relatively small range, which is similar to that reported in the literature [35].

##### 4.2. Effect of Fracture Length on REV and Equivalent Permeability

Figure 12(a) depicts the relationship between and . The varies in a relatively small range for and increases significantly when increases from 4.0 to 4.5. When varies from 2.0 to 4.0, the varies in a small range (from 1.96E+05 m^{3} to 1.19E+05 m^{3}). The (6.70E+05 m^{3}) for is approximately 3.4 times larger than that (1.96E+05 m^{3}) for . The main reason for the rapidly increasing of is that the flow paths are constructed by small fractures when and the small fractures that are far less than the model size increase the heterogeneity of the fracture network.

**(a)**Relationship between and

**(b)**Relationship between andFigure 12(b) indicates that decreases with the increment of , since a larger increases the proportion of small fractures and the small fractures contribute less continuous pathways for the flow. The has a significant influence on the determination of the at the REV size than the , especially for , which indicates that the at the REV size is more sensitive than when the of the fracture networks keep the same.

#### 5. Conclusion

The present study applies the equivalent pipe network (EPN) modeling method for calculating the fluid flow through the three-dimensional (3D) discrete fracture network (DFN). The validity of the EPN approach is verified by comparing the simulated equivalent permeability () with that reported in the literature. A series of 3D DFN models are generated to investigate the influence of power exponent of fracture length distribution () on the representative elementary volume (REV) size and the at the REV size. The number of fractures used to generate the DFN models is determined by analyzing the relationship between the number of fractures () and fracture density ().

The results show that the calculated by EPN models agrees well with that reported in the literature, which confirms the validity of the EPN model for calculating the of 3D DFN models. With the increment of , the values of are less scattered and show a clear tendency of reaching a REV size. When increases from 2.0 to 4.0, the REV size of 3D DFN models () varies in a relatively small range (from 1.96E+05 m^{3} to 1.19E+05 m^{3}), whereas the () is 3.4 times larger than that when . The rapid increment of results from the inhomogeneous fluid flow in which the small nonpersistent fractures dominate the preferential flow paths when exceeds a certain value (i.e., 4.5). The at the REV size decreases following an exponential distribution with the increment of because the proportion of small fractures increases. The small fractures are less likely to contribute continuous pathways for the fluid flow in the fracture network.

This study assumed that the location and orientation are uniform and random, while the two parameters might follow some other distributions according to the field investigation results. Future studies will investigate the evolution of DFN models’ hydraulic properties with different fracture locations and orientations.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (51974296, 52061135111, and 51874277).