Correlations between Geometric Properties and Permeability of 2D Fracture Networks

*e 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. *is 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. *e results show that the equivalent permeability has power law relationships with the geometric properties of fracture networks. *e 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. *e equivalent permeability of a fracture network increases with the increment of fracture density and/or fractal dimensions proportionally.


Introduction
e estimation of equivalent permeability of fractured rock masses plays an important role in many environmental and engineering applications, such as CO 2 sequestration, underground nuclear waste repositories, and geothermal and heat storages [1][2][3][4]. e 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][6][7][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][10][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][13][14]. e 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. e 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. e performance of conventional statistical parameters on predicting equivalent permeability of DFNs is compared with that using fractal dimensions.

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 Table 1. 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. us, 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. e 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. ree 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. e basic assumptions are made that the fluid flow only occurs in fractures and obeys the cubic law, and the rock matrix is impermeable. e equivalent permeability is calculated under two boundary conditions (horizontal and vertical) with a constant hydraulic gradient of 10 kPa/m. e equivalent permeability is back-calculated using the following equation: 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.

Fractal Evaluation
e fractal dimensions of the fracture network structures represented by D f and the intersection points represented by D in 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. e process of calculating the fractal dimension of DFNs is described in [10] in detail. e calculated fractal dimensions for three examples are shown in Figure 1. e value of D f is always larger than 1, while D in may have values less than 1 because D in 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.

Results and Analysis
Figures 2(a) ∼ 2(d) show the relationships between K and ρ in , ρ seg , d m , and d in , respectively. Here, ρ in is the number density of intersections, ρ seg is the number density of segments, d m is the mass density of fractures, and d in is the number of intersections per meter at the boundary. e 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 ρ. e flow path of a fracture network is strongly influenced by ρ seg and d m . A larger value of ρ seg or d m 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. e orientation of fractures is depicted by d in . A larger value of d in 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, d in at the horizontal boundaries is larger than 0, whereas d in in the direction perpendicular to horizontal is 0, resulting in the stronger permeability in the horizontal direction. e geometric parameters, i.e., ρ in , ρ seg , d m , and d in , 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 > 10 0 .
is 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. e 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 , d m , and d in ), 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, D f and D in , are utilized to represent the distributions of fractures and intersections, respectively. e geometry of a fracture network can be redepicted based on the images of fracture outcrops in a field [15]. en, D f and D in 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 R 2 > 0.96 as follows:    K � 9.27 × 10 − 9 ρ 0.55 in , R 2 � 0.9885, K � 8.10 × 10 − 9 D 2.16 in , R 2 � 0.9813. Substituting equations (2) and (4) into equations (6) and (7), respectively, yields Equations (8) ∼ (9) imply that D f and D in are also correlated with d m and ρ in following power law functions, with R 2 > 0.97. Figure 3 shows the comparisons between the predicted results of K, D f , and D in using equations (2) ∼ (9) with the calculated results obtained from flow simulations (K) and fractal evaluations (D f and D in ). e results show that the predicted values agree well with the calculated results, and although when D f < 1.15 and D in < 0.4, the predicted values of D f and D in are slightly larger than the calculated results, due to the effect of randomness of fracture locations. erefore, the permeability of a DFN consisting of lognormally distributed fractures can be well estimated using the geometric parameters and fractal dimensions. Figure 4 shows the variations of K + and V K + with ρ varying from ρ 0 to 64ρ 0 , in which K + and V K + are defined as follows: where K + is the dimensionless equivalent permeability, K v and K h are the equivalent permeability of the vertical flow (Figure 1(c)) and horizontal flow (Figure 1(a)), respectively, V K + represents the variation in K + for different DFNs, and K + max and K + min are the maximum and minimum values of K + , respectively.
When ρ is small, i.e., ρ < 16 ρ0 , K + is more scattered, and the value of V K + 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 V K + is smaller. K + ranges from 0.23 to 1.81 when ρ � ρ 0 and varies from 0.94 to 1.12 when ρ � 64 ρ 0 . e effect of randomness of fracture locations can be neglected when ρ > 24 ρ0 , corresponding to D f > 1.55 and D in > 1.07, above which the permeability of a DFN increases with the increment of fracture density and/or fractal dimension.

Conclusions
is 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. e mathematical expressions were proposed to link the equivalent permeability with the geometric parameters and fractal dimensions that describe the geometric characters of fracture networks. e 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.

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

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.