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

Authors: Jing Zhang; Liyuan Yu; Hongwen Jing; Richeng Liu
Journal: Geofluids
 (2018)

---

## Abstract

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.


---

## Body

## 1. 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–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 stress-displacement 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.


## 2. 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]:
(1)Q=ρge312μΔhLfW,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, Lf 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–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 Df 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–25]; Miao et al. 2015). The total number of fractures, Nt, can be calculated as follows:
(2)NtL≥lmin=lmaxlminDf/2,where Nt is the cumulative number of fractures, lmin is the minimum fracture length, and lmax is the maximum fracture length. In addition, lmin≪lmax is the necessary condition for the fracture length distribution to follow the fractal scaling law. lmin/lmax≪10−3 is used in the present study. The length of the ith fracture can be calculated as:
(3)li=lmin1−Ri2/Df,where li is the length of the ith fracture, i=1,2,3,…,Nt, Nt is the total number which has been obtained from (2), and Ri 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 Ln of DFNs can be expressed as [24, 25]:
(4)Nt'=1.73NtLn2Df7.14∑i=1Ntli.After updating the total number of fractures fromNt to Nt′, the fractal length distribution of fractures in any Ln 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)).


## 3. 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–32] and [13] discrete fracture flow model ([20]; Oda 1985; [28, 33]). The first approach assumes that the combined 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.


### 3.2. Discrete Fracture Network Generation

The Monte Carlo simulation technique is used for the generation of the DFNs. The theory and technique of the issue have been well studied and presented in previous studies [27, 38–40]. The process of generation includes the following steps:
(1)
Several parameters, such aslmin, lmax, Df, and Ln, should be determined before DFNs are generated. lmin and lmax 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 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 (Ln) of 50 m, 25 m, 15 m, 10 m, 8 m, and 4 m are generated when Df=1.3, 1.4, 1.5, 1.6, and 1.7, respectively. Ln decreases with the increase in Df. If a large Ln (i.e., 50 m) is adopted when Df is large (i.e., Df=1.7), there would be many fractures generated, which may cost a longer time for solving fluid flow, comparing with that with a smaller Df (i.e., Df=1.3).
(2)
The value ofNt can be calculated using (1). The length of the ith fracture can be calculated using (3) after generating Nt random numbers. Nt′ can be obtained by calculating (4). After that, the lengths of all fracture with a number of Nt′ can be generated using (3).
(3)
The DFNs are established using a DEM (discrete element method) code based on a 2D open source software (OSS) UDEC (Itasca Consulting Group Inc. 2004), in which the fractures are represented with line segments.
(4)
Rectangular submodels are inserted into UDEC to calculate REV with differentDf. The submodels, whose area is the square of Ln, are extracted from original large models with Ln=5~45m, 1~20 m, 0.5~10 m, 0.25~7 m, and 0.25~4 m, corresponding to Df=1.3, 1.4, 1.5, 1.6, and 1.7, respectively. The aperture for each fracture is a constant (65 μm) in order to study the effect of fracture length on the permeability. Here, we take Df=1.6 as an example to show the details of generation of original DFNs and extraction of submodels as shown in Figure 1. Figure 1(a) shows the process of extracting submodels from an original DFN with Ln=10 m and Lx:Ly=1.5:1, in which Lx and Ly represent the length and width of the submodels, respectively. Figure 1(b) exhibits the extraction of submodels that have different length-width ratios while maintaining the same area. Table 1 shows the two sides of rectangular submodels compared to the length of square submodels with the same area when Df = 1.6.
(5)
Constant hydraulic head boundary conditions are applied, as shown in Figure2, to numerically calculate the permeability (Kp) 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.Figure 1
Schematic view of the generation of submodels.
(a)
Extraction of submodels from an original DFNs
(b)
Submodels that have different length-width ratiosTable 1
Two sides of rectangular submodels compared to the length of square submodels with the same area whenDf = 1.6.
Length-width ratio (rectangular models)
Lx:Ly=1:1
Lx:Ly=1:2
Lx:Ly=1:1.5
Lx:Ly=1.5:1
Lx:Ly=2:1
Lx=Ly (m)
Lx (m)
Ly (m)
Lx (m)
Ly (m)
Lx (m)
Ly (m)
Lx (m)
Ly (m)
0.25
0.177
0.354
0.204
0.306
0.306
0.204
0.354
0.177
0.5
0.354
0.707
0.408
0.612
0.612
0.408
0.707
0.354
1
0.707
1.414
0.816
1.224
1.224
0.816
1.414
0.707
2
1.414
2.828
1.632
2.449
2.449
1.632
2.828
1.414
3
2.121
4.242
2.449
3.674
3.674
2.449
4.242
2.121
4
2.828
5.569
3.265
4.898
4.898
3.265
5.569
2.828
5
3.535
7.071
4.082
6.123
6.123
4.082
7.071
3.535
6
4.243
8.485
4.898
7.348
7.348
4.898
8.485
4.243
7
4.949
9.899
5.715
8.573
8.573
5.715
9.899
4.949Figure 2
Hydraulic boundary conditions.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 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.


### 3.3. 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. Figure3 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 Df follows a power law function that can be expressed by:
(5)nl,0.5=αl−a,where α is the proportional coefficient, a is the power law exponent, and nl,0.5 represents the fracture number with lengths in the range of l−0.5,l+0.5.Figure 3
Correlation between fracture number and fracture length.WhenDf 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.Table 2
Comparison of the values ofa and achievements in previous studies.
Authors
Year
Value ofa
Dverstop and Anderson
1989
1.7
Tsang et al.
1996
3.0
Bour and Davy
1997
1.0–3.0
De Dreuzy
2001
0.0–3.5
Richeng Liu
2016
1.17–3.39
Present study
2018
0.96–1.29


## 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–32] and [13] discrete fracture flow model ([20]; Oda 1985; [28, 33]). The first approach assumes that the combined 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.


## 3.2. Discrete Fracture Network Generation

The Monte Carlo simulation technique is used for the generation of the DFNs. The theory and technique of the issue have been well studied and presented in previous studies [27, 38–40]. The process of generation includes the following steps:
(1)
Several parameters, such aslmin, lmax, Df, and Ln, should be determined before DFNs are generated. lmin and lmax 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 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 (Ln) of 50 m, 25 m, 15 m, 10 m, 8 m, and 4 m are generated when Df=1.3, 1.4, 1.5, 1.6, and 1.7, respectively. Ln decreases with the increase in Df. If a large Ln (i.e., 50 m) is adopted when Df is large (i.e., Df=1.7), there would be many fractures generated, which may cost a longer time for solving fluid flow, comparing with that with a smaller Df (i.e., Df=1.3).
(2)
The value ofNt can be calculated using (1). The length of the ith fracture can be calculated using (3) after generating Nt random numbers. Nt′ can be obtained by calculating (4). After that, the lengths of all fracture with a number of Nt′ can be generated using (3).
(3)
The DFNs are established using a DEM (discrete element method) code based on a 2D open source software (OSS) UDEC (Itasca Consulting Group Inc. 2004), in which the fractures are represented with line segments.
(4)
Rectangular submodels are inserted into UDEC to calculate REV with differentDf. The submodels, whose area is the square of Ln, are extracted from original large models with Ln=5~45m, 1~20 m, 0.5~10 m, 0.25~7 m, and 0.25~4 m, corresponding to Df=1.3, 1.4, 1.5, 1.6, and 1.7, respectively. The aperture for each fracture is a constant (65 μm) in order to study the effect of fracture length on the permeability. Here, we take Df=1.6 as an example to show the details of generation of original DFNs and extraction of submodels as shown in Figure 1. Figure 1(a) shows the process of extracting submodels from an original DFN with Ln=10 m and Lx:Ly=1.5:1, in which Lx and Ly represent the length and width of the submodels, respectively. Figure 1(b) exhibits the extraction of submodels that have different length-width ratios while maintaining the same area. Table 1 shows the two sides of rectangular submodels compared to the length of square submodels with the same area when Df = 1.6.
(5)
Constant hydraulic head boundary conditions are applied, as shown in Figure2, to numerically calculate the permeability (Kp) 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.Figure 1
Schematic view of the generation of submodels.
(a)
Extraction of submodels from an original DFNs
(b)
Submodels that have different length-width ratiosTable 1
Two sides of rectangular submodels compared to the length of square submodels with the same area whenDf = 1.6.
Length-width ratio (rectangular models)
Lx:Ly=1:1
Lx:Ly=1:2
Lx:Ly=1:1.5
Lx:Ly=1.5:1
Lx:Ly=2:1
Lx=Ly (m)
Lx (m)
Ly (m)
Lx (m)
Ly (m)
Lx (m)
Ly (m)
Lx (m)
Ly (m)
0.25
0.177
0.354
0.204
0.306
0.306
0.204
0.354
0.177
0.5
0.354
0.707
0.408
0.612
0.612
0.408
0.707
0.354
1
0.707
1.414
0.816
1.224
1.224
0.816
1.414
0.707
2
1.414
2.828
1.632
2.449
2.449
1.632
2.828
1.414
3
2.121
4.242
2.449
3.674
3.674
2.449
4.242
2.121
4
2.828
5.569
3.265
4.898
4.898
3.265
5.569
2.828
5
3.535
7.071
4.082
6.123
6.123
4.082
7.071
3.535
6
4.243
8.485
4.898
7.348
7.348
4.898
8.485
4.243
7
4.949
9.899
5.715
8.573
8.573
5.715
9.899
4.949Figure 2
Hydraulic boundary conditions.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 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.


## 3.3. 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. Figure3 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 Df follows a power law function that can be expressed by:
(5)nl,0.5=αl−a,where α is the proportional coefficient, a is the power law exponent, and nl,0.5 represents the fracture number with lengths in the range of l−0.5,l+0.5.Figure 3
Correlation between fracture number and fracture length.WhenDf 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.Table 2
Comparison of the values ofa and achievements in previous studies.
Authors
Year
Value ofa
Dverstop and Anderson
1989
1.7
Tsang et al.
1996
3.0
Bour and Davy
1997
1.0–3.0
De Dreuzy
2001
0.0–3.5
Richeng Liu
2016
1.17–3.39
Present study
2018
0.96–1.29


## 4. Results and Analysis

### 4.1. 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 differentDf. Next, the equivalent permeability is back-calculated using Darcy’s law, as follows:
(6)Kf=μQA∇P,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 Section3.2, the variations in Kf and corresponding root mean square (RMS) are plotted in Figures 4–8. For each Df and each length-width ratio, 10 DFNs are generated using 10 sets of random numbers. Then, Kf and its mean value of the 10 DFNs are calculated. Kf 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 Ln increases to a critical value and the change in Kf is sufficiently slight, can be regarded as the REV size. VK and RMS decrease with the increment of Ln. Here, VK 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.Figure 4
Variations inKp and RMS with different Ln when Df=1.3.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1Figure 5
Variations inKp and RMS with different Ln when Df=1.4.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1The RMS is introduced to quantify whether two datasets agree well with each other [3], defined as:
(7)RMS=21n∑1nKavg−Ksim2Kmax+Kmin,where Kmax is the maximum equivalent permeability, Kmin is the minimum equivalent permeability, Kavg is the average equivalent permeability, and Ksim is the simulation result.Taking the case ofDf=1.6 in Figure 6(a) for an example, Kp ranges from 8.5E − 14 to 5.5E − 13 m2 when Ln=0.5 m, while for Ln=10m, Kp ranges from 3.5E − 13 to 3.1E − 13 m2. Besides, VK 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 Kp of submodels with the same area.Figure 6
Variations inKp and RMS with different Ln when Df=1.5.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1Figure 7
Variations inKp and RMS with different Ln when Df=1.6.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1Figure 8
Variations inKp and RMS with different Ln when Df=1.7.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1Figure 9
Relationship betweenDf and LREV for different Lx:Ly.Figure 10
Relationship betweenKREV and Lx:Ly for different Df.Figure 11
Relationship betweenLREV and Df.In this study,RMS=0.1 is utilized as the threshold value, and the corresponding side length of DFN (LREV) 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 LREV gradually decreases with the increment of Df from 1.3 to 1.7 as shown in Figure 9. When Df is small (i.e., Df=1.3), the model length-width ratio can influence the magnitude of LREV significantly (i.e., from 28.12 to 35.16 m). When Df is large (i.e., Df=1.7), the influence of the model length-width ratio on LREV 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 Df is small (i.e., Df≤14). By fitting these datasets, an exponential function is proposed between LREV and Df, written as:
(8)LREV=0.44+1.24E7∗4.88E−5Df.Thus, whenDf 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 Lx:Ly=1:1, showing an agreement with each other. When Df is small (i.e., 1.3), the difference in LREV between the predicted results and those in Liu et al. [25] is relatively larger (i.e., 8.51 m). However, with the increment of Df from 1.4 to 1.6, the difference is very small (from 2.39 m to 0.43 m). This is acceptable because when Df is small, the randomness of fracture distribution is significant and therefore LREV changes in a large extent. With the increment of Df (i.e., Df>1.4), the effect of randomness of fracture distribution decreases, and as a result, LREV varies slightly. This again verifies the validity of (8).


### 4.2. Effects ofDf and Model Length-Width Ratio on KREV

The variation ofKREV is plotted against the length-width ratio in Figure 10, where KREV is the equivalent permeability at the REV size. As the model length-width ratio of submodels increases, KREV gradually decreases. Besides, the rate of permeability reduction also decreases with increasing model length-width ratio. The effect of the model length-width ratio on KREV when Df=1.3 is much smaller than that when Df=1.7. The variation in KREV is plotted against Df in Figure 11. The influence of Df is determinate, as what have been observed for the influence of the model length-width ratio. The reason may be that the DFNs with larger Df have better connectivity and stronger conductivity/permeability. KREV increases with the increase in Df, while the variations in the KREV increment for Df≥1.6 are much larger than those for Df<1.6.


## 4.1. 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 differentDf. Next, the equivalent permeability is back-calculated using Darcy’s law, as follows:
(6)Kf=μQA∇P,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 Section3.2, the variations in Kf and corresponding root mean square (RMS) are plotted in Figures 4–8. For each Df and each length-width ratio, 10 DFNs are generated using 10 sets of random numbers. Then, Kf and its mean value of the 10 DFNs are calculated. Kf 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 Ln increases to a critical value and the change in Kf is sufficiently slight, can be regarded as the REV size. VK and RMS decrease with the increment of Ln. Here, VK 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.Figure 4
Variations inKp and RMS with different Ln when Df=1.3.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1Figure 5
Variations inKp and RMS with different Ln when Df=1.4.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1The RMS is introduced to quantify whether two datasets agree well with each other [3], defined as:
(7)RMS=21n∑1nKavg−Ksim2Kmax+Kmin,where Kmax is the maximum equivalent permeability, Kmin is the minimum equivalent permeability, Kavg is the average equivalent permeability, and Ksim is the simulation result.Taking the case ofDf=1.6 in Figure 6(a) for an example, Kp ranges from 8.5E − 14 to 5.5E − 13 m2 when Ln=0.5 m, while for Ln=10m, Kp ranges from 3.5E − 13 to 3.1E − 13 m2. Besides, VK 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 Kp of submodels with the same area.Figure 6
Variations inKp and RMS with different Ln when Df=1.5.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1Figure 7
Variations inKp and RMS with different Ln when Df=1.6.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1Figure 8
Variations inKp and RMS with different Ln when Df=1.7.
(a)
Lx:Ly=1:2
(b)
Lx:Ly=1:2
(c)
Lx:Ly=1:1.5
(d)
Lx:Ly=1:1.5
(e)
Lx:Ly=1.5:1
(f)
Lx:Ly=1.5:1
(g)
Lx:Ly=2:1
(h)
Lx:Ly=2:1Figure 9
Relationship betweenDf and LREV for different Lx:Ly.Figure 10
Relationship betweenKREV and Lx:Ly for different Df.Figure 11
Relationship betweenLREV and Df.In this study,RMS=0.1 is utilized as the threshold value, and the corresponding side length of DFN (LREV) 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 LREV gradually decreases with the increment of Df from 1.3 to 1.7 as shown in Figure 9. When Df is small (i.e., Df=1.3), the model length-width ratio can influence the magnitude of LREV significantly (i.e., from 28.12 to 35.16 m). When Df is large (i.e., Df=1.7), the influence of the model length-width ratio on LREV 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 Df is small (i.e., Df≤14). By fitting these datasets, an exponential function is proposed between LREV and Df, written as:
(8)LREV=0.44+1.24E7∗4.88E−5Df.Thus, whenDf 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 Lx:Ly=1:1, showing an agreement with each other. When Df is small (i.e., 1.3), the difference in LREV between the predicted results and those in Liu et al. [25] is relatively larger (i.e., 8.51 m). However, with the increment of Df from 1.4 to 1.6, the difference is very small (from 2.39 m to 0.43 m). This is acceptable because when Df is small, the randomness of fracture distribution is significant and therefore LREV changes in a large extent. With the increment of Df (i.e., Df>1.4), the effect of randomness of fracture distribution decreases, and as a result, LREV varies slightly. This again verifies the validity of (8).


## 4.2. Effects ofDf and Model Length-Width Ratio on KREV

The variation ofKREV is plotted against the length-width ratio in Figure 10, where KREV is the equivalent permeability at the REV size. As the model length-width ratio of submodels increases, KREV gradually decreases. Besides, the rate of permeability reduction also decreases with increasing model length-width ratio. The effect of the model length-width ratio on KREV when Df=1.3 is much smaller than that when Df=1.7. The variation in KREV is plotted against Df in Figure 11. The influence of Df is determinate, as what have been observed for the influence of the model length-width ratio. The reason may be that the DFNs with larger Df have better connectivity and stronger conductivity/permeability. KREV increases with the increase in Df, while the variations in the KREV increment for Df≥1.6 are much larger than those for Df<1.6.


## 5. Conclusions

In this study, we calculated the equivalent permeability of rock fracture networks that have different models length-width 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 inDf according to the exponential relationship between the REV size and Df. The effect of the length-width ratio of submodels on the REV size decreases with increasing Df for Df<1.5, but can be negligible when Df≥1.5. KREV decreases with the increments of the model length-width ratio of submodels and/or Df. The influence of the model length-width ratio on KREV when Df is small (i.e., Df=1.3) is much smaller than that when Df is large (i.e., Df=1.7). Besides, the variations in the increment of KREV when Df≥1.6 are much larger than those when Df<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.



---

*Source: 7206074-2018-10-29.xml*