# Numerical Investigations on the Effect of Fracture Length Distribution on the Representative Elementary Volume of 3D Discrete Fracture Networks

Authors: Ding Liu; Hai Pu; Shiru Guo; Ziheng Sha; Chong Li
Journal: Geofluids
 (2022)

---

## Abstract

Determination of the representative elementary volume (REV) of fractured rock masses based on equivalent permeability (K) 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 K 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 K of the 3D DFN model has an exponential relationship with the power exponent (a) of fracture length distribution and the evaluation of K 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 (Ln) is small, the K varies significantly due to the influence of random number seeds used to generate fracture length, location, and orientation. The K holds a constant value after Ln 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 VREV. The VREV varies within a narrow range when a≤4.0. When a=4.5, the VREV rapidly increases to more than 3.4 times than that when a=4.0. The fluid flow becomes more inhomogeneous due to the small nonpersistent fractures that dominate the preferential flow paths when a exceeds a certain value (i.e., 4.5). The K at the REV size decreases exponentially with the increment of a. This tendency can be explained by the decrease of the average intersection length (Li) with the increment of a, which is a geometric parameter for reflecting the connectivity of the fracture network.


---

## Body

## 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 (a) 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 (P32) and different a were established to calculate the K 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 a on REV size and K of fracture networks at the REV size was systematically investigated. The flow chart for describing this study is shown in Figure 1.Figure 1
Flow chart.


## 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 (P32) 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
(1)nl=α∙l−a,where α is the coefficient of proportionality, l is the fracture length which distributed between the lower (lmin) and upper (lmax) bounds, and a 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 Figure2(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).Figure 2
Generation and verification of the EPN model.
(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 ofK with different aTable 1
Parameters used for the generation of DFNs in Section2 (after Huang et al. [24]).
ParametersDistributionParametersDomain size (m)—Ln=30 (m)PositionsUniform random—Orientations (°)FisherSet 1:DD=10, DA=90, and κ=200Set 2: DD=120, DA=40, and κ=200Set 3: DD=70, DA=60, and κ=200Set 4: DD=320, DA=80, and κ=200Length (m)Power law2<a<4.5, 5<Lf<30Density (m2/m3)—P32=0.4DA stands for dip angle, DD represents dip direction, and κ means Fisher constant.The 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 (Figure2(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 (AX=B) is set up for every inner node as follows:
(2)Hj=∑i=1nCij·Hi∑i=1nCij,n=2,3,4,where Hj is the pressure head of node j, Cij is the conductance of the pathway between node i and node j, and Hi is the pressure head of node i. Cij can be expressed as
(3)Cij=ga3b12υLp,where g is the gravity acceleration, a is the aperture which is equal to the aperture of fracture in this study, b is third of fracture dimension, v is the kinematic viscosity of the fluid, and Lp is the length of pipe. The expressions for Hj and Cij 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 (AX=B) is fully prepared, a simple linear arithmetic technology (X=A/B) 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
(4)Q=CΔH,where Q is flow rate and ΔH is pressure head difference between two nodes. After the flow rate is obtained, the K can be back-calculated according to Darcy’s law, as follows:
(5)K=μQAρgJ,where A 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 Table1, 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 a, 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.


## 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 (P32) 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
(1)nl=α∙l−a,where α is the coefficient of proportionality, l is the fracture length which distributed between the lower (lmin) and upper (lmax) bounds, and a 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 Figure2(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).Figure 2
Generation and verification of the EPN model.
(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 ofK with different aTable 1
Parameters used for the generation of DFNs in Section2 (after Huang et al. [24]).
ParametersDistributionParametersDomain size (m)—Ln=30 (m)PositionsUniform random—Orientations (°)FisherSet 1:DD=10, DA=90, and κ=200Set 2: DD=120, DA=40, and κ=200Set 3: DD=70, DA=60, and κ=200Set 4: DD=320, DA=80, and κ=200Length (m)Power law2<a<4.5, 5<Lf<30Density (m2/m3)—P32=0.4DA stands for dip angle, DD represents dip direction, and κ means Fisher constant.The 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 (Figure2(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 (AX=B) is set up for every inner node as follows:
(2)Hj=∑i=1nCij·Hi∑i=1nCij,n=2,3,4,where Hj is the pressure head of node j, Cij is the conductance of the pathway between node i and node j, and Hi is the pressure head of node i. Cij can be expressed as
(3)Cij=ga3b12υLp,where g is the gravity acceleration, a is the aperture which is equal to the aperture of fracture in this study, b is third of fracture dimension, v is the kinematic viscosity of the fluid, and Lp is the length of pipe. The expressions for Hj and Cij 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 (AX=B) is fully prepared, a simple linear arithmetic technology (X=A/B) 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
(4)Q=CΔH,where Q is flow rate and ΔH is pressure head difference between two nodes. After the flow rate is obtained, the K can be back-calculated according to Darcy’s law, as follows:
(5)K=μQAρgJ,where A 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 Table1, 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 a, 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 (N) ranging from 6000 to 15000 are generated to evaluate the relationship between the number of fractures and the fracture density (P32) of the 3D DFN model. The expression of P32 can be written as
(6)P32=∑SV,where S is the fracture area in the sampling volume and V 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 N and P32 for different a. The P32 increases with the increment of N in a linear way for different a. The parameters of fitting curves as shown in Figure 3 are listed in Table 3. And the variation of parameter c with a is demonstrated in Figure 3(b). Figure 3(b) shows that the increasing rate of P32 (c) decreases with the increasing a. Note that a smaller a would lead to more larger fractures and the variation of P32 is more sensitive to a smaller a. The N for different a when P32=0.5 m2/m3 is determined by the fitting curve in Figure 3(a) and demonstrated in Figure 3(c). The N with P32=0.5 m2/m3 increases (from 6351 to 12800) with the increasing a.Table 2
Parameters used for the generation of DFNs.
ParametersDistributionParametersDomain size (m)ConstantL=100 (m)PositionsUniform random—Orientations (°)FisherDA ~ [0,π/2)DD ~ [0, 2π)Length (m)Power law2<a<4.5, 5<Lf<30Density (m2/m3)ConstantP32=0.5Figure 3
Relationships betweenP32 and number of fractures for different a.
(a)
Variation ofP32 with number of fractures(b)
Variation ofc with a(c)
Number of fractures for different a withP32=0.5Table 3
Parameters of fitting curves as shown in Figure3.
abcR22.0-1.3200E-037.9028E-050.99992.58.6509E-046.7083E-050.99983.0-7.7872E-045.7719E-050.99993.5-5.2400E-035.0560E-050.99994.05.9300E-034.3726E-050.99974.52.9600E-033.8837E-050.9997The average length of fracture intersections (Li) of 3D DFN models with different N and same a changes in relatively small range as shown in Figure 4(a). Figure 4(b) shows that the average value of Li decreases with the increment of a following an exponential relationship. The equivalent permeability of the 3D models for different N decreases exponentially with increasing a 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 a. 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 a is better than that with larger a when N keeps the same. In addition, the K of DFN models with larger N must be greater than that with smaller N when a 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 a on K at the REV size will be investigated in the next section.Figure 4
Variation ofLi with the number of fractures and a.
(a)
Variation ofLi with number of fractures(b)
Variation of the average value ofLi with aFigure 5
Relationships betweenK and a for different N.Table 4
Parameters of fitting curves as shown in Figure5.
NdfR260004.4555E-04-10.03280.563370002.4670E-11-1.46260.985280001.4264E-10-1.95280.996790005.8037E-11-1.43120.9873100001.0347E-10-1.52970.9889110009.0353E-11-1.32760.9870120008.8187E-11-1.22480.9803130001.3115E-10-1.29080.9980140001.0064E-10-1.09440.9797150001.2774E-10-1.11180.9961


### 3.2. Generation and Extraction of 3D DFN Submodel with Differenta

A series of original 3D DFN models with differenta were generated as shown in Figure 6. For each a, 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 Li decreases with increasing a 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 a on the REV size of fracture networks.Figure 6
Three-dimensional DFN model with differenta.
(a)
a=2.0(b)
a=2.5(c)
a=3.0(d)
a=3.5(e)
a=4.0(f)
a=4.5Figure 7
DFN models generated using ten sets of random number seeds.
(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 10Figure 8
Relationships betweenLi and a.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 (H) 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 H1=1 m and the right side is outside with a free boundary (H2=0 m), while the other boundaries are impermeable. The left side is the yz plane for x=0 m and the right side is the yz plane for x 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 (H=1 m).Figure 9
Examples of extracting submodels from the original DFN model.
(a)
Original DFN model (Ln=100 m)(b)
Submodel (Ln=30 m)(c)
Submodel (Ln=50 m)(d)
Submodel (Ln=80 m)Figure 10
Hydraulic boundary conditions of cube submodels and corresponding distributions of water heads.
(a)
A three-dimensional DFN submodel(b)
Flow direction (left to right)(c)
Water head distribution


## 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 (N) ranging from 6000 to 15000 are generated to evaluate the relationship between the number of fractures and the fracture density (P32) of the 3D DFN model. The expression of P32 can be written as
(6)P32=∑SV,where S is the fracture area in the sampling volume and V 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 N and P32 for different a. The P32 increases with the increment of N in a linear way for different a. The parameters of fitting curves as shown in Figure 3 are listed in Table 3. And the variation of parameter c with a is demonstrated in Figure 3(b). Figure 3(b) shows that the increasing rate of P32 (c) decreases with the increasing a. Note that a smaller a would lead to more larger fractures and the variation of P32 is more sensitive to a smaller a. The N for different a when P32=0.5 m2/m3 is determined by the fitting curve in Figure 3(a) and demonstrated in Figure 3(c). The N with P32=0.5 m2/m3 increases (from 6351 to 12800) with the increasing a.Table 2
Parameters used for the generation of DFNs.
ParametersDistributionParametersDomain size (m)ConstantL=100 (m)PositionsUniform random—Orientations (°)FisherDA ~ [0,π/2)DD ~ [0, 2π)Length (m)Power law2<a<4.5, 5<Lf<30Density (m2/m3)ConstantP32=0.5Figure 3
Relationships betweenP32 and number of fractures for different a.
(a)
Variation ofP32 with number of fractures(b)
Variation ofc with a(c)
Number of fractures for different a withP32=0.5Table 3
Parameters of fitting curves as shown in Figure3.
abcR22.0-1.3200E-037.9028E-050.99992.58.6509E-046.7083E-050.99983.0-7.7872E-045.7719E-050.99993.5-5.2400E-035.0560E-050.99994.05.9300E-034.3726E-050.99974.52.9600E-033.8837E-050.9997The average length of fracture intersections (Li) of 3D DFN models with different N and same a changes in relatively small range as shown in Figure 4(a). Figure 4(b) shows that the average value of Li decreases with the increment of a following an exponential relationship. The equivalent permeability of the 3D models for different N decreases exponentially with increasing a 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 a. 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 a is better than that with larger a when N keeps the same. In addition, the K of DFN models with larger N must be greater than that with smaller N when a 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 a on K at the REV size will be investigated in the next section.Figure 4
Variation ofLi with the number of fractures and a.
(a)
Variation ofLi with number of fractures(b)
Variation of the average value ofLi with aFigure 5
Relationships betweenK and a for different N.Table 4
Parameters of fitting curves as shown in Figure5.
NdfR260004.4555E-04-10.03280.563370002.4670E-11-1.46260.985280001.4264E-10-1.95280.996790005.8037E-11-1.43120.9873100001.0347E-10-1.52970.9889110009.0353E-11-1.32760.9870120008.8187E-11-1.22480.9803130001.3115E-10-1.29080.9980140001.0064E-10-1.09440.9797150001.2774E-10-1.11180.9961


## 3.2. Generation and Extraction of 3D DFN Submodel with Differenta

A series of original 3D DFN models with differenta were generated as shown in Figure 6. For each a, 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 Li decreases with increasing a 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 a on the REV size of fracture networks.Figure 6
Three-dimensional DFN model with differenta.
(a)
a=2.0(b)
a=2.5(c)
a=3.0(d)
a=3.5(e)
a=4.0(f)
a=4.5Figure 7
DFN models generated using ten sets of random number seeds.
(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 10Figure 8
Relationships betweenLi and a.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 (H) 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 H1=1 m and the right side is outside with a free boundary (H2=0 m), while the other boundaries are impermeable. The left side is the yz plane for x=0 m and the right side is the yz plane for x 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 (H=1 m).Figure 9
Examples of extracting submodels from the original DFN model.
(a)
Original DFN model (Ln=100 m)(b)
Submodel (Ln=30 m)(c)
Submodel (Ln=50 m)(d)
Submodel (Ln=80 m)Figure 10
Hydraulic boundary conditions of cube submodels and corresponding distributions of water heads.
(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. Figure11 shows the calculated results from ten realizations as the scattered data points of equivalent permeability, along with the bounding curves of their maximum (Kmax), minimum (Kmin), mean (Kavg), the difference between the maximum and minimum (Vk) 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 (Ln). 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 (P32=0.5 m2/m3), RMS=0.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 VREV.Figure 11
Variations inK and RMS with different a in cube submodels.
(a)
Relationship betweenK and Ln with a=2.0(b)
Relationship between RMS andLn with a=2.0(c)
Relationship betweenK and Ln with a=2.5(d)
Relationship between RMS andLn with a=2.5(e)
Relationship betweenK and Ln with a=3.0(f)
Relationship between RMS andLn with a=3.0(g)
Relationship betweenK and Ln with a=3.5(h)
Relationship between RMS andLn with a=3.5(i)
Relationship betweenK and Ln with a=4.0(j)
Relationship between RMS andLn with a=4.0(k)
Relationship betweenK and Ln with a=4.5(l)
Relationship between RMS andLn with a=4.5Figure11 presents the variations in K of cube submodels with different Ln and corresponding RMS. The K fluctuates dramatically because of the impacts of random number seeds used to create fracture location, orientation, and length when Ln is small. When Ln crosses specific thresholds, K remains constant regardless Ln, allowing the model scale to be used as the REV size and the VREV to be determined. The K equals to 1.092E-12 m2, 8.1E-13 m2, 5.18E-13 m2, 4.7E-13 m2, 2.54E-13 m2, and 1.78E-13 m2 when RMS=0.3 is reached at the approximated VREV size of 1.96E+05 m3, 1.88E+05 m3, 1.19E+05 m3, 1.84E+05 m3, 1.96E+05 m3, and 6.70E+05 m3 for the cases in Figures 11(a), 11(c), 11(e), 11(g), 11(i), and 11(k), respectively. The K becomes less distributed with increasing Ln, with a noticeable trend to attain a REV size. Note that the K of DFN models decreases rapidly with the increment of a. 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 a. The framework is gradually predominantly structured by smaller fractures. The RMS continuously decreases with increasing Ln, especially before the Ln 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 ofLn, 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

Figure12(a) depicts the relationship between VREV and a. The VREV varies in a relatively small range for a=2.0~4.0 and increases significantly when a increases from 4.0 to 4.5. When a varies from 2.0 to 4.0, the VREV varies in a small range (from 1.96E+05 m3 to 1.19E+05 m3). The VREV (6.70E+05 m3) for a=4.5 is approximately 3.4 times larger than that (1.96E+05 m3) for a=4.0. The main reason for the rapidly increasing of VREV is that the flow paths are constructed by small fractures when a=4.5 and the small fractures that are far less than the model size increase the heterogeneity of the fracture network.Figure 12
Variations inVREV and K with increasing a from 2.0 to 4.5 m for P32=0.5 m2/m3.
(a)
Relationship betweenVREV and a(b)
Relationship betweenK and aFigure12(b) indicates that VREV decreases with the increment of a, since a larger a increases the proportion of small fractures and the small fractures contribute less continuous pathways for the flow. The a has a significant influence on the determination of the K at the REV size than the VREV, especially for a≤4.0, which indicates that the K at the REV size is more sensitive than VREV when the P32 of the fracture networks keep the same.


## 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. Figure11 shows the calculated results from ten realizations as the scattered data points of equivalent permeability, along with the bounding curves of their maximum (Kmax), minimum (Kmin), mean (Kavg), the difference between the maximum and minimum (Vk) 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 (Ln). 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 (P32=0.5 m2/m3), RMS=0.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 VREV.Figure 11
Variations inK and RMS with different a in cube submodels.
(a)
Relationship betweenK and Ln with a=2.0(b)
Relationship between RMS andLn with a=2.0(c)
Relationship betweenK and Ln with a=2.5(d)
Relationship between RMS andLn with a=2.5(e)
Relationship betweenK and Ln with a=3.0(f)
Relationship between RMS andLn with a=3.0(g)
Relationship betweenK and Ln with a=3.5(h)
Relationship between RMS andLn with a=3.5(i)
Relationship betweenK and Ln with a=4.0(j)
Relationship between RMS andLn with a=4.0(k)
Relationship betweenK and Ln with a=4.5(l)
Relationship between RMS andLn with a=4.5Figure11 presents the variations in K of cube submodels with different Ln and corresponding RMS. The K fluctuates dramatically because of the impacts of random number seeds used to create fracture location, orientation, and length when Ln is small. When Ln crosses specific thresholds, K remains constant regardless Ln, allowing the model scale to be used as the REV size and the VREV to be determined. The K equals to 1.092E-12 m2, 8.1E-13 m2, 5.18E-13 m2, 4.7E-13 m2, 2.54E-13 m2, and 1.78E-13 m2 when RMS=0.3 is reached at the approximated VREV size of 1.96E+05 m3, 1.88E+05 m3, 1.19E+05 m3, 1.84E+05 m3, 1.96E+05 m3, and 6.70E+05 m3 for the cases in Figures 11(a), 11(c), 11(e), 11(g), 11(i), and 11(k), respectively. The K becomes less distributed with increasing Ln, with a noticeable trend to attain a REV size. Note that the K of DFN models decreases rapidly with the increment of a. 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 a. The framework is gradually predominantly structured by smaller fractures. The RMS continuously decreases with increasing Ln, especially before the Ln 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 ofLn, 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

Figure12(a) depicts the relationship between VREV and a. The VREV varies in a relatively small range for a=2.0~4.0 and increases significantly when a increases from 4.0 to 4.5. When a varies from 2.0 to 4.0, the VREV varies in a small range (from 1.96E+05 m3 to 1.19E+05 m3). The VREV (6.70E+05 m3) for a=4.5 is approximately 3.4 times larger than that (1.96E+05 m3) for a=4.0. The main reason for the rapidly increasing of VREV is that the flow paths are constructed by small fractures when a=4.5 and the small fractures that are far less than the model size increase the heterogeneity of the fracture network.Figure 12
Variations inVREV and K with increasing a from 2.0 to 4.5 m for P32=0.5 m2/m3.
(a)
Relationship betweenVREV and a(b)
Relationship betweenK and aFigure12(b) indicates that VREV decreases with the increment of a, since a larger a increases the proportion of small fractures and the small fractures contribute less continuous pathways for the flow. The a has a significant influence on the determination of the K at the REV size than the VREV, especially for a≤4.0, which indicates that the K at the REV size is more sensitive than VREV when the P32 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 (K) 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 (a) on the representative elementary volume (REV) size and the K 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 (N) and fracture density (P32).The results show that theK calculated by EPN models agrees well with that reported in the literature, which confirms the validity of the EPN model for calculating the K of 3D DFN models. With the increment of Ln, the values of K are less scattered and show a clear tendency of reaching a REV size. When a increases from 2.0 to 4.0, the REV size of 3D DFN models (VREV) varies in a relatively small range (from 1.96E+05 m3 to 1.19E+05 m3), whereas the VREV (a=4.5) is 3.4 times larger than that when a=4.0. The rapid increment of VREV results from the inhomogeneous fluid flow in which the small nonpersistent fractures dominate the preferential flow paths when a exceeds a certain value (i.e., 4.5). The K at the REV size decreases following an exponential distribution with the increment of a 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.



---

*Source: 8073013-2022-06-09.xml*