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.
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 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. 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 to 64ρ0, and for each , 10 sets of different random numbers are utilized. Here, 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.
(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: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 by and the intersection points represented by 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 is always larger than 1, while may have values less than 1 because 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
Figures 2(a) ∼ 2(d) show the relationships between K and , , , and , respectively. Here, is the number density of intersections, is the number density of segments, dm is the mass density of fractures, and 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 and . The flow path of a fracture network is strongly influenced by and . A larger value of or 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 . A larger value of 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, at the horizontal boundaries is larger than 0, whereas in the direction perpendicular to horizontal is 0, resulting in the stronger permeability in the horizontal direction.
(a)
(b)
(c)
(d)
(e)
(f)
The geometric parameters, i.e., , , , and , 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., , 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., . 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., , , , and ), 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, and , 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, and 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:
Substituting equations (2) and (4) into equations (6) and (7), respectively, yields
Equations (8) ∼ (9) imply that and are also correlated with and following power law functions, with R2 > 0.97.
Figure 3 shows the comparisons between the predicted results of K, , and using equations (2) ∼ (9) with the calculated results obtained from flow simulations (K) and fractal evaluations ( and ). The results show that the predicted values agree well with the calculated results, and although when and , the predicted values of and 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.
(a)
(b)
(c)
Figure 4 shows the variations of and with ρ varying from ρ0 to 64ρ0, in which and are defined as follows:where is the dimensionless equivalent permeability, and are the equivalent permeability of the vertical flow (Figure 1(c)) and horizontal flow (Figure 1(a)), respectively, represents the variation in for different DFNs, and and are the maximum and minimum values of , respectively.
When is small, i.e., , is more scattered, and the value of is larger, comparing with the cases with a larger , i.e., , in which converges to a smaller range of variations, and the value of is smaller. ranges from 0.23 to 1.81 when and varies from 0.94 to 1.12 when . The effect of randomness of fracture locations can be neglected when , corresponding to and , 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., , , , , , and . The values of and can be calculated using the box-counting method to represent and with power law functions. The predicted values of K, , and are consistent with those of the calculated results, indicating that the proposed expressions of K with the geometric properties of DFNs are reliable. When and , 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.
Nomenclature
| A: | Cross-sectional area (L2) |
| K: | Equivalent permeability (L2) |
| P: | Hydraulic pressure (ML−1T−2) |
| : | Equivalent permeability of horizontal flow (L2) |
| : | Equivalent permeability of vertical flow (L2) |
| : | Dimensionless equivalent permeability |
| : | Variation of dimensionless equivalent permeability |
| : | Dynamic viscosity (M/(L·T)) |
| : | Number density of fracture lines (L−2) |
| : | Initial number density of fracture lines (L−2) |
| : | Fractal dimension of fracture network structures |
| : | Fractal dimension of fracture intersections |
| : | Number density of fracture intersections (L−2) |
| : | Number density of fracture segments (L−2) |
| : | Length of fractures per square meter (L.L−2) |
| : | Number of fractures intersecting with the boundary per meter (L−1). |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was partially funded by the National Natural Science Foundation of China (Grant no. 51579239).