# Correlations between Geometric Properties and Permeability of 2D Fracture Networks

Authors: Xiaolin Wang; Liyuan Yu; Hanqing Yang
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The equivalent permeability of fractured rock masses plays an important role in understanding the fluid flow and solute transport properties in underground engineering, yet the effective predictive models have not been proposed. This study established mathematical expressions to link permeability of 2D fracture networks to the geometric properties of fractured rock masses, including number density of fracture lines, total length of fractures per square meter, and fractal dimensions of fracture network structures and intersections. The results show that the equivalent permeability has power law relationships with the geometric properties of fracture networks. The fractal dimensions that can be easily obtained from an engineering site can be used to predict the permeability of a rock fracture network. When the fractal dimensions of fracture network structures and intersections exceed the critical values, the effect of randomness of fracture locations is negligible. The equivalent permeability of a fracture network increases with the increment of fracture density and/or fractal dimensions proportionally.


---

## Body

## 1. Introduction

The estimation of equivalent permeability of fractured rock masses plays an important role in many environmental and engineering applications, such as CO2 sequestration, underground nuclear waste repositories, and geothermal and heat storages [1–4]. The previous studies have revealed that the equivalent permeability of fracture networks is strongly correlated with the geometric properties of fractures, such as fracture length, density, aperture, and orientation [2, 5–8]. Recently, some works have shown that the fractal dimension is an effective tool to describe the geometric properties of fracture networks, and it is quantitatively correlated with the equivalent permeability [9–11].A number of predictive models have been proposed to calculate the equivalent permeability of both porous and fractured media. In these models, the length of fractures has different types of distributions, such as random distribution [8, 9], power-law distribution [6, 7], and fractal-like tree distribution [12–14]. The lognormal distribution is another important geometric description of the fracture length in networks, which, however, was not considered in any previous predictive models.In this study, the lognormal distribution of fracture length is incorporated into the discrete fracture network (DFN) models with different number densities of fracture lines, and the equivalent permeability of each model is calculated via fluid flow simulations. The mathematical expressions that link equivalent permeability to the geometric properties, i.e., connectivity, density, and orientation, are proposed, and their validities are verified by comparing the predicted results with the calculated results of DFNs. The performance of conventional statistical parameters on predicting equivalent permeability of DFNs is compared with that using fractal dimensions.


## 2. DFN Generation and Fluid Flow Simulation

To generate a DFN model using the Monte Carlo method, the following three aspects need to be addressed. First, the geometric parameters (i.e., length, aperture, orientation, and location) of each fracture need to be assigned, assuming that the fracture length follows a lognormal distribution and the orientation follows a normal distribution, as shown in Table1. Second, an original model with a side length of 150 m is generated, followed by extractions from the target model with a side length of 100 m from the original model. Thus, the unevenly distributed fractures close to the boundaries are deleted from the model. Because the fluid flow only takes place in the connected fractures, the fracture segments located out of the model, the isolated fracture segments, and the “dead-end” elements of fractures, which are the ends of fractures that are not connected with other fractures and do not contribute to the fluid flow, are deleted from the model. Finally, the mass continuity equations at each fracture intersection are solved utilizing an iteration scheme under a given boundary condition. In the present study, the aperture has a constant value of 0.18 mm for each fracture, and the location of a fracture is randomly and uniformly distributed. Here, a constant value of apertures is taken into account because the present study aims to investigate the relationships between fractal dimensions and hydraulic properties of fracture networks. The number density of fracture lines, ρ, varies from ρ0 to 64ρ0, and for each ρ, 10 sets of different random numbers are utilized. Here, ρ0 is the initial number density of fracture lines, which have values of 0.00075 m−2 and 0.00025 m−2 for the two fracture sets, respectively. Three examples of the established DFNs and distribution of fracture intersections are shown in Figure 1. For each model, both the horizontal and vertical flows are simulated under a constant horizontal hydraulic gradient (see Figure 1(a)) and a vertical hydraulic gradient (see Figure 1(c)), respectively.Table 1
Geometric parameters used for DFN generation.
SetFracture length (m)Orientation (°)Aperture (mm)DFN side length (m)Initial number density,ρ0 (m−2)Avg.Dev.Avg.Dev.153.65107.77343.950.181000.000752688.7279.75286.620.181000.00025Figure 1
Distributions of fractures ((a), (c), and (e)) and intersections ((b), (d), and (f)) in three DFNs with different densities, and the two kinds of boundary conditions applied to DFNs in which the flow directions are perpendicular to each other. (a)ρ=64ρ0 and Df=1.62. (b) ρ=64ρ0 and Din=1.28. (c) ρ=16ρ0 and Df=1.49. (d) ρ=16ρ0 and Din=0.95. (e) ρ=ρ0 and Df=1.16. (f) ρ=ρ0 and Din=0.39.
(a)(b)(c)(d)(e)(f)The basic assumptions are made that the fluid flow only occurs in fractures and obeys the cubic law, and the rock matrix is impermeable. The equivalent permeability is calculated under two boundary conditions (horizontal and vertical) with a constant hydraulic gradient of 10 kPa/m. The equivalent permeability is back-calculated using the following equation:(1)Q=AKμ∂P∂L,where Q is the flow rate, A is the cross-sectional area, K is the permeability, μ is the dynamic viscosity, L is the fracture length, and P is the hydraulic pressure.


## 3. Fractal Evaluation

The fractal dimensions of the fracture network structures represented byDf and the intersection points represented by Din are calculated using the box-counting method. Each image of fracture network structures or intersections is constituted by 300 pixels in length and 300 pixels in width, and each pixel has a value of either 1 or 0 (null). An image is covered with square boxes of different dimensions from 300 × 300 pixels (1 box) to 1 × 1 pixel (300 × 300 boxes). If there is a part of a fracture in a box, this box will be assigned with a value of 1, otherwise 0. A log-log plot of the box count vs. the number of total boxes can be plotted, and the slope represents the fractal dimension. The process of calculating the fractal dimension of DFNs is described in [10] in detail. The calculated fractal dimensions for three examples are shown in Figure 1. The value of Df is always larger than 1, while Din may have values less than 1 because Din describes the fractal properties of dispersedly distributed intersection points in a plane, which may have values larger or smaller than 1 depending on the density of points.


## 4. Results and Analysis

Figures2(a) ∼ 2(d) show the relationships between K and ρin, ρseg, dm, and din, respectively. Here, ρin is the number density of intersections, ρseg is the number density of segments, dm is the mass density of fractures, and din is the number of intersections per meter at the boundary. The connectivity of a DFN, defined as the total number of intersection points divided by the total number of fracture lines [9], can be represented by a function of ρin and ρ. The flow path of a fracture network is strongly influenced by ρseg and dm. A larger value of ρseg or dm indicates a denser fracture distribution and a larger number of intersections, resulting in a larger number of flow paths connecting the inlet boundary to the outlet boundary. The orientation of fractures is depicted by din. A larger value of din at an orientation means a larger number of fractures intersecting with the boundary and a stronger conductivity. For a limit case, when all of the fractures are horizontal, din at the horizontal boundaries is larger than 0, whereas din in the direction perpendicular to horizontal is 0, resulting in the stronger permeability in the horizontal direction.Figure 2
Relationships betweenK and (a) ρin, (b) ρseg, (c) dm, (d) din, (e) Df, and (f) Din, respectively.
(a)(b)(c)(d)(e)(f)The geometric parameters, i.e.,ρin, ρseg, dm, and din, increase with the increment of ρ, and K of DFNs increases with these geometric parameters following power-law functions, as shown in equations (2) ∼ (5). When the values of the geometric parameters are small, i.e., ρin<10−1, the value of K varies within a large range, which gradually converges to the best-fitted curves when the values of the geometric parameters become large, i.e., ρin>100. This is because when the density of fractures and intersections is small, the location of fractures is significantly influenced by the randomly distributed fracture center point that can have a large dispersion. The effect of randomness of center points on the geometric characters of networks diminishes with increasing fracture density [10]. Although the value of K is well estimated utilizing the geometric parameters (i.e., ρin, ρseg, dm, and din), it is still a challenging and time-consuming task to obtain their values for a real-fractured rock mass. To facilitate the process, the fractal dimensions, Df and Din, are utilized to represent the distributions of fractures and intersections, respectively. The geometry of a fracture network can be redepicted based on the images of fracture outcrops in a field [15]. Then, Df and Din can be easily calculated using the box-counting method. Figures 2(e) ∼ 2(f) exhibit that K is also correlated with the fractal dimensions following power law functions, as shown in equations (6) ∼ (7). All of the calculated data fit equations (2) ∼ (7) fairly well with the correlation coefficient R2 > 0.96 as follows:(2)K=9.27×10−9ρin0.55,R2=0.9885,(3)K=6.86×10−9ρseg0.53,R2=0.9902,(4)K=3.64×10−9dm0.55,R2=0.9859,(5)K=1.49×10−8din1.02,R2=0.9668,(6)K=4.36×10−10Df7.14,R2=0.9746,(7)K=8.10×10−9Din2.16,R2=0.9813.Substituting equations (2) and (4) into equations (6) and (7), respectively, yields(8)Df=1.33dm0.079,R2=0.9790,(9)Din=1.06ρin0.25,R2=0.9830.Equations (8) ∼ (9) imply that Df and Din are also correlated with dm and ρin following power law functions, with R2 > 0.97.Figure3 shows the comparisons between the predicted results of K, Df, and Din using equations (2) ∼ (9) with the calculated results obtained from flow simulations (K) and fractal evaluations (Df and Din). The results show that the predicted values agree well with the calculated results, and although when Df<1.15 and Din<0.4, the predicted values of Df and Din are slightly larger than the calculated results, due to the effect of randomness of fracture locations. Therefore, the permeability of a DFN consisting of lognormally distributed fractures can be well estimated using the geometric parameters and fractal dimensions.Figure 3
Comparisons between predicted and calculated results of (a)K, (b) Df, and (c) Din. (a) Predicted K by equations (1)–(6) vs. calculated K. (b) Predicted Df by equation (7) vs. calculated Df. (c) Predicted Din by equation (8) vs. calculated Din.
(a)(b)(c)Figure4 shows the variations of K+ and VK+ with ρ varying from ρ0 to 64ρ0, in which K+ and VK+ are defined as follows:(10)K+=KvKh,(11)VK+=Kmax+−Kmin+,where K+ is the dimensionless equivalent permeability, Kv and Kh are the equivalent permeability of the vertical flow (Figure 1(c)) and horizontal flow (Figure 1(a)), respectively, VK+ represents the variation in K+ for different DFNs, and Kmax+ and Kmin+ are the maximum and minimum values of K+, respectively.Figure 4
Variations inK+ and VK+ with varying ρ.Whenρ is small, i.e., ρ<16ρ0, K+ is more scattered, and the value of VK+ is larger, comparing with the cases with a larger ρ, i.e., ρ>24ρ0, in which K+ converges to a smaller range of variations, and the value of VK+ is smaller. K+ ranges from 0.23 to 1.81 when ρ=ρ0 and varies from 0.94 to 1.12 when ρ=64ρ0. The effect of randomness of fracture locations can be neglected when ρ>24ρ0, corresponding to Df>1.55 and Din>1.07, above which the permeability of a DFN increases with the increment of fracture density and/or fractal dimension.


## 5. Conclusions

This study established a great number of discrete fracture networks (DNFs) with different fracture densities and random numbers for generating the center point and orientation of fractures, in which the length of fractures follows a lognormal distribution. The mathematical expressions were proposed to link the equivalent permeability with the geometric parameters and fractal dimensions that describe the geometric characters of fracture networks. The validity was verified by comparing the predicted results with the calculated results of DFN simulations.The results show that the equivalent permeability has power law relationships with the geometric properties of DFNs, i.e.,ρin, ρseg, dm, din, Df, and Din. The values of Df and Din can be calculated using the box-counting method to represent dm and ρin with power law functions. The predicted values of K, Df, and Din are consistent with those of the calculated results, indicating that the proposed expressions of K with the geometric properties of DFNs are reliable. When Df>1.55 and Din>1.07, the effect of randomness of fracture locations are negligible, and the permeability of a DFN increases with the increment of fracture density and/or fractal dimensions.The future works will focus on the effect of aperture variation on the equivalent permeability of 2D fracture networks. Besides, we will also study the relationship between the equivalent permeability and fractal dimension of 3D rough fracture networks because the 3D fracture network can better characterize the roughness, orientation, connectivity, and permeability tensor of real-fractured rock masses than 2D fracture networks. Equations (2) ∼ (9) are applicable for 2D fracture networks, whose applicability to 3D fracture networks will be investigated in future works.



---

*Source: 6645238-2021-01-27.xml*