Abstract

Estimating equivalent permeability at grid block scale of numerical models is a critical issue for large-scale fractured porous rocks. However, it is difficult to constrain the permeability distributions for equivalent fracture models as these are strongly influenced by complex fracture properties. This study quantitatively investigated equivalent permeability distributions for fractured porous rocks, considering the impact of the correlated fracture aperture and length model. Two-dimensional discrete fracture models are generated with varied correlation exponent ranges from 0.5 to 1, which indicates different geomechanical properties of fractured porous rock. The equivalent fracture models are built by the multiple boundary upscaling method. Results indicate that the spatial distribution of equivalent permeability varied with the correlation exponent. When the minimum fracture length and the number of fractures increase, the process that the diagonal equivalent permeability tensor components change from a power law like to a lognormal like and to a normal-like distribution slows down as the correlation exponent increases. The average dimensionless equivalent permeability for the equivalent fracture models is well described by an exponential relationship with the correlation exponent. A power law model is built between the equivalent permeability of equivalent fracture models and fracture density of discrete fracture models for the correlated aperture-length models. The results demonstrate that both the fracture density and length-aperture model influence the equivalent permeability of equivalent fracture models interactively.

1. Introduction

Fractures are among the most common structures in the brittle rocks of the Earth’s crust and span a wide range of length scales from millimeter to kilometers. The interconnected fracture networks in the rock matrix are the primary pathway for fluid to flow underground and have significant importance on practical applications, such as aquifer management, groundwater contamination, hydrocarbon and geothermal energy exploitation, geological disposal of nuclear waste, and coal mine water inrush [16]. As a quantitative measure of a rock mass’s ability to conduct fluid, permeability is one of the most important properties affecting the flow in fractured porous rocks.

Due to the high heterogeneity and anisotropy of fractured rocks, the permeability of fractured porous rock generally has a symmetric full tensor form [7] and exhibits scale effects [8, 9]. Permeability of fractured rocks can be estimated by laboratory experiment [10], borehole flowmeter [11], field pumping test [12], discrete fracture characterization data from core samples [13], borehole well logs [14], and outcrops [15]. Such measurements can be roughly classified into continuum approaches and discrete fracture approaches [16]. The continuum approaches are often applied to explore the bulk rock flow regime for a scale of the measurement. In contrast, numerical or analytical calculating equivalent permeability from discrete fracture geometries provides an efficient way for estimating rock permeability and links the small scale fracture geometries and equivalent permeability for field-scale models, which is also known as upscaling [17, 18]. An understanding of the influence of discrete fractures on the equivalent permeability of grid blocks opens the possibility of creating efficient and accurate equivalent fracture models based on the multiple-scale fracture characterization data.

For investigating the correlation between equivalent permeability and fracture geometries, analytical models and numerical models were developed since the 1960s [19] and since 1980s [20], respectively. With the development of fracture characterization, numerical models, and theoretic models, continuous efforts have been made for correlating the equivalent permeability and fracture geometries [2123]. de Dreuzy et al. [24] studied the influence of a power law length distribution and a lognormal aperture distribution on the equivalent permeability for fractured rocks with an impervious matrix. Min et al. [25] developed a numerical procedure to calculate the permeability tensors of fractured rocks, using a stochastic REV concept based on discrete fracture network models. Based on a newly developed correlation equation, Baghbanan and Jing [26] investigated the permeability of fractured rocks, considering the correlation between distributed fracture aperture and length. Klimczak et al. [27] modeled flow through fracture networks for both correlated and uncorrelated fracture length-to-aperture relationships, calculated permeability of fractured rocks, and confirmed the importance of the correlated square root relationship of the aperture to length scaling. Leung and Zimmerman [28] established a methodology for estimating the macroscopic effective hydraulic conductivity based on the geometric parameters of the fracture network rather than solving the flow equations. However, compared to porous rocks (e.g., [2931]), determining the bulk properties of fractured rocks is still challenging regarding the complexity of fracture geometries.

More recently, Miao et al. [32] derived an analytical model for estimating permeability for fractured rocks based on the fractal geometry theory and the cubic law for laminar flow in fractures. Liu et al. [33] proposed a fractal model to link the fractal characteristics of both the fluid flow tortuosity and fracture geometry with the equivalent permeability of the fracture networks. Hyman et al. [34] characterized how different fracture size and aperture relationships influence flow and transport simulations through sparse three-dimensional discrete fracture networks and observe that networks with a correlated relationship have consistently higher effective permeability values. The above study assumed that rock matrix permeability is impervious. Bisdom et al. [35] compared the different aperture models and critical stress criteria and conclude that their impacts on equivalent permeability depend on matrix permeability for naturally fractured reservoirs. Xiong et al. [36] developed a numerical procedure about nonlinear flow in three-dimension discrete fracture networks by solving the Reynolds equation and Forchheimer equation. They found that the required hydraulic head gradient for the onset of nonlinear fluid flows reduces when the fracture surface becomes smooth and the fracture connectivity increases. Yin et al. [37] developed a high-precision apparatus to investigate the influence of shear processes on nonlinear flow within three-dimensional rough-walled fractures experimentally and found that the hydraulic aperture of fracture enlarges increases as the shear displacement increases. The studies are conducted at the scale of the laboratory measurement. Nevertheless, the complexity of the interconnected fractures at the field scale has yet to be fully considered.

While many previous studies propose relationships between fracture length and aperture and illustrated their effects on fluid flow through fractured rocks (e.g., [33, 35, 38, 39]), yet they mainly address the equivalent permeability for a grid block at different scales. The influence of correlated length and aperture on equivalent permeability distributions of grid blocks for large-scale fractured porous media has yet to be clarified. The novelty of this study is that equivalent permeability distribution with correlated fracture aperture and length is analyzed for fractured porous media with state of the art upscaling method [40], which provides a potential link between fracture properties and stochastic equivalent fracture models and is possible for reducing uncertainty for equivalent fracture models.

The paper is organized as follows: first, the discrete fracture models with correlated fracture aperture and length are presented, and the multiple boundary upscaling method is introduced. Then, the equivalent permeability distribution for the discrete fracture models is analyzed stochastically based on histograms. Third, the relationship between the average equivalent permeability and correlated aperture and length model is interpreted. Forth, the effects of fracture geometry on the equivalent permeability are investigated based on the correlated aperture and length model. Finally, the results are discussed, and the conclusions are drawn.

2. Methodology

Two main procedures were conducted for analyzing the effects of fracture length-aperture correlation on equivalent permeability distribution for the two-dimensional fractured porous rocks: generating discrete fracture models correlating the fracture length and aperture and upscaling equivalent permeability for the discrete fracture model by using the multiple boundary method. The related techniques are described below.

2.1. Discrete Fracture Model Generation

Fracture apertures vary due to the mechanical misfit of fracture walls, chemical dissolution, and normal pressure due to the depth of overburden [41]. Fracture apertures cover a wide range of scales. They can be measured by a wide variety of methods, including direct measurements in core or outcrop and deduction from flow data. The aperture for fractured porous rocks can be described by constant (e.g., [24]), statistical models [42], correlation with fracture [27], and mechanical models [38]. The power law relationship between fracture aperture and length (including linear and sublinear) applied in this study is widely observed at the field scale [43] and is derived by linear elastic fracture mechanics [44], which can be expressed by [45]

where is fracture aperture, is fracture length, is the coefficient related to mechanical properties of fractured rock, and the correlation exponent indicates the mechanical interaction between closely spaced fractures (Figure 1). For isolated veins, faults, and shear deformation bands with constant driving stress, the correlation between fracture aperture and length tends to be linear, i.e., =1.0 [45, 46]. For more complex open-mode fractures with a constant fracture toughness, the ability of a rock containing a fracture to resist further fracturing, is around 0.5 [44, 47]. The correlation exponent is greater than 0.5, and less than 1.0 could result from post-jointing relaxation [27]. In this study, the correlation exponent ranges from 0.5 to 1.0, and is assumed to be according to the field data [48].

Fracture length covers a wide scale and can be described by a power law model [41]: where is the fracture length, is the number of fractures with sizes in the range , is a density constant, and is a power law exponent. As the length of fractures is measured upon a specific scale and is limited by the rock size in the Earth’s crust, it has lower and upper bounds for the power law distribution. The power law exponent represents the growth properties of the fractures and varies from 1.3 to 3.5 [41, 49], and in this study, is assumed to be 2.5.

In this study, a series of two-dimensional synthetic fracture networks is generated in a squared domain of size  m (Figure 2). The location and orientation of fractures are assumed purely random, which is nominally homogeneous. The fracture length follows power law distribution as shown in equation (1) and the lower bound and the upper bound are 4 m and 40 m, respectively. The number of fractures is 50. The fracture geometries vary by increasing from 4 m to 6 m, 8 m, and 10 m and by rising from 50 to 75, 100, and 125. Thus, a total of 16 sets of discrete fracture networks with different geometric parameters are created with the Monte Carlo method based mainly on the ADFNE software [50]. In the following study, both the permeability of fractures and rock matrix are taken into account, which constitutes discrete fracture models.

2.2. Upscaling Permeability

The equivalent fracture models are built based on upscaling permeability for the discrete fracture models. The dimension of the Cartesian grid for the equivalent fracture model is . After subdividing the fractures of the whole domain into the Cartesian grid, the multiple boundary upscaling method is applied for calculating equivalent permeability for each grid block. The multiple boundary upscaling method is applicable for both two-dimensional and three-dimensional discrete fracture models, which was tested for a rotated fracture by fitting with analytical solutions as well as for fracture networks [40].

The multiple boundary upscaling method involves mainly three steps. First, the linear boundary conditions, which mimic the flow across fractured porous rocks underground, are applied for the grid block with the pressure gradient of 1 Pa/m along the -axis (Figure 3). Then, the steady flow problem is solved, and the fluxes and from both the fracture and rock matrix are calculated with a multiple boundary expression [40]. Lastly, the equivalent permeability component and is calculated according to Darcy’s law inversely. The equivalent permeability for two-dimensional discrete fracture models is a second-rank symmetric tensor with four components. and can be calculated by changing the direction of the linear boundary condition along the -axis. It should be noted that the upscaled equivalent permeability is no inherently symmetric, and the symmetric permeability tensor is acquired by averaging the off-diagonal components.

The MRST code [51] is applied for solving flow problems in the fractures and rock matrix during the upscaling procedure. Flow equations in the fractures and rock matrix are based on mass conservation and Darcy’s law, which are coupled based on the flow continuity on the fracture-rock matrix interface and are solved with a multipoint flux approximation [52]. The rock matrix and the fracture are represented by two-dimensional triangular grids of 0.2 m and one-dimensional line grids of 0.1 m, respectively. The matrix permeability is assumed as (1 md) according to the actual data of fractured hydrocarbon reservoirs [53], and the fracture permeability is calculated according to the fracture aperture, that is, . The numerical upscaling procedure is validated against analytical solutions. Considering a fully penetrated fracture with different correlation exponent , the equivalent permeability can be calculated analytically and numerically with the multiple boundary method (Figure 4). It is shown that the numerical solutions match well with the analytical solution with varied correlation exponent .

3. Equivalent Permeability Distribution

In this section, the influence of length-aperture correlation parameters on the histograms of equivalent permeability is demonstrated. The aperture-length correlation exponent varied from 0.5 to 1.0 with a step of 0.1, which indicates fracture mechanical state changes from constant fracture toughness to constant fracture driving force. For a specific correlation exponent , 16 discrete fracture networks with varied fracture network geometries, i.e., the minimum fracture length and the number of fractures , are generated. For generalizing the results, ten realizations are created for each group of fracture geometric parameters and thus a total of 960 discrete fracture models are generated for the following analysis.

Figure 5(a) shows the Cartesian grid of the equivalent fracture model for a realization of the discrete fracture models with , and . The spatial distributions of equivalent permeability components , , and are plotted in Figures 5(b)–5(d), respectively. The equivalent permeability ranges from (matrix permeability) to . For the diagonal components and , the high permeable grid blocks are penetrated by long fractures which have large apertures as shown in the discrete fracture model. However, the spatial distribution for and is different from each other due to the random orientation of the fractures. For the off-diagonal component in Figure 5(d), they are either positive or negative, and the grid blocks with high absolute values correspond to those of large and in Figures 5(b) and 5(c).

For investigating the statistical distribution of equivalent permeability tensors, the histograms of , , and are plotted in Figures 5(e)–5(g), respectively. The fitting curves for histograms are also plotted. It is shown that both and exhibit a power law-like distribution whereas exhibits a normal distribution, which is similar to the results in three-dimensional fractured porous rock with poorly fracture connectivity [54].

For the discrete fracture model with different fracture geometric parameters, the fitting curves of the equivalent permeability tensor histograms for , , and are plotted in Figure 6. Rather than a single realization, the curves in Figure 6 represent fitting of all ten realizations for a given set of fracture geometric parameters. For , it shows that the shape and the magnitude of the diagonal components and are similar and tend to a power law distribution, whereas tends to a normal distribution. With the increase in and , and change gradually from a power law to a lognormal to a normal distribution. The lognormal distribution of equivalent permeability is commonly assumed in the realistic reservoir models [55, 56]. Such transformation of the equivalent permeability distribution is partially due to the connectivity of fractures within grid blocks [49]. The shape of expands, which means an increase of the equivalent permeability. Furthermore, for evaluating the strength of the relationship for the fitting curves in Figure 6, the Spearman rank correlation coefficient between equivalent permeability and the frequency is calculated. The correlation coefficients for and are similar, with an average of about -0.75, which indicates a relatively strong negative correlation. For , the coefficient is less than -0.04, which almost denotes no correlation. It is reasonable as the off-diagonal component tends to be evenly distributed around zero due to the random fracture orientations.

With the increase in the correlation exponent , the magnitude of the equivalent permeability tensor components also increases, which is due to that the equivalent permeability of the grid blocks is more controlled by long fractures with high aperture values [26]. However, with the increase in and , the evolution of the histogram fitting curves differs for varied correlation exponent . For a higher (e.g.,), the histograms do not show a clear normal distribution for large and compared to those for a lower , e.g., They look like a transition shape for a lower changing from a power law to a lognormal to a normal distribution. The difference is partially due to that the high correlation exponent increases the heterogeneity of the discrete fracture model. And it is reasonable to speculate that with the increase in and , the diagonal components tend to a lognormal distribution or normal distribution for a high . Accordingly, for a specific fracture geometry, a high correlation exponent means a higher heterogeneity which requires a sufficient number of fractures to reach a lognormal like distribution.

4. Correlating the Equivalent Permeability with the Aperture-Length Model

Before studying the effect of the length-aperture correlation exponent on the equivalent permeability for the equivalent fracture models, the average dimensionless permeability for each discrete fracture model is defined as

where is the total number of grid blocks for the equivalent fracture model, and the rock matrix permeability is equal to 1 md as assumed in the previous section. Under a specific aperture-length correlation exponent , the average dimensionless permeability for all 160 realizations of discrete fracture models with varied fracture geometries is calculated. Figure 7 shows the relationship between the average dimensionless permeability and the aperture-length correlation exponent of the discrete fracture models.

It shows in Figures 7(a) and 7(b) that for the diagonal components of equivalent permeability tensors, log and log, increase linearly with the increase in the correlation exponent , which can be expressed as the following exponential function:

where the dimensionless coefficients and range from 101.4 to 102.4 and from 3.4 to 4.3, respectively. For the off-diagonal component , as the negative values exist, the absolute value is used to analyze the correlation between the average dimensionless permeability and the aperture-length correlation exponent (Figure 7(c)). Likewise, the correlation can also be built in the framework of equation (4), which indicates that the exponential model can describe the relationship between the average dimensionless permeability and the aperture-length correlation exponent . The observations also agree with the conclusions of Bisdom et al. [35] that the aperture-fracture model has an important effect on the magnitude of equivalent permeability.

We should note that although the correlation between of the equivalent fracture models and of the discrete fracture models can be described by equation (4), the dimensionless coefficient for and increases as the fracture geometric parameters and raise (Figures 7(a) and 7(b)). In contrast, the dimensionless coefficient does not show a clear change. It is reasonable as the large fracture length and number of fractures increase the permeability and connectivity of fractures and result in increased equivalent permeability. The coefficients and for are of the same magnitude as those for and . However, for the coefficient , the range is wider than those of and , which may be mainly due to the process of averaging the off-diagonal components. For the coefficient , it is smaller than those of diagonal components. It is reasonable as the off-diagonal components are smaller than the diagonal components for a given discrete fracture model as shown in Figures 5(b)–5(d).

5. Effects of Fracture Geometry on the Equivalent Permeability

The influence of fracture geometries on the average dimensionless equivalent permeability for correlated length-aperture models is analyzed in this section. For each discrete fracture model, the dimensionless fracture density is defined as [28]

where is the domain area, is the length of the -th fracture, and is the total number of fractures of the discrete fracture model.

The correlations between the log() and log are plotted for different length-aperture correlation exponents in Figure 8. It shows that for the diagonal components (Figures 8(a) and 8(b)), log() increases linearly with the increase in log for a given length-aperture correlation exponent , which can be fitted by the following power law relationship:

where and are dimensionless coefficient in the ranges 103-105 and 1.1-1.3, respectively. Whereas for the off-diagonal components (Figure 8(c)), they are more scattered compared to the diagonal components. This characteristic is similar to that shown in Figure 7(c). It is shown in Figure 8 that the average dimensionless equivalent permeability scatters widely around the fitting line when the correlation exponent increases. This is mainly because a higher correlation exponent results in an enlarged difference in equivalent permeability. It indicates further that the higher in the correlated length-aperture model, the higher heterogeneity for the equivalent fracture permeability field. It should be noted that closing to 1 means a linear correlation between and , which has been observed in the results for macroscopic hydraulic conductivity by Leung and Zimmerman [28]. Furthermore, when the correlation exponent increases, and also increase (Figure 9). It shows that the aperture-length model influences the dimensionless coefficients and for the power model described by equation (6).

6. Discussion

In this study, the equivalent permeability distribution for discrete fracture models with correlated fracture aperture and length is investigated. Our results support the results of Klimczak et al. [27], Leung and Zimmerman [28], and Bisdom et al. [35], who highlight the importance of varied aperture models on equivalent permeability. Furthermore, the matrix permeability is considered for the fractured porous rocks as well in this study. In contrast to estimating the equivalent permeability at a given scale as in the above studies, this study investigates the statistical properties of equivalent permeability and the relationship between the averaged equivalent permeability rocks and the discrete fracture properties based on the multiple boundary upscaling method [54]. The methodology in the paper could be applied for efficient estimation of equivalent permeability distributions for fractured porous rocks correlating fracture length and aperture, which can be incorporated into stochastic equivalent fracture models. Thus, the results are meaningful for reducing uncertainty in numerical models for groundwater flow, solute, and heat transport processes (e.g., [56]) and for identifying key fracture geometric properties influencing rock hydraulic properties before fracture sampling at the field scale [2, 57].

The correlation exponent varies from 0.5 to 1.0, which indicates fracture change from constant fracture toughness to constant driving stress [27, 44, 46]. The equivalent permeability for a fractured rock mass could be highly anisotropic, and the directional permeability can be estimated by the rotation of grid blocks [58]. In this study, the orientation and degree of anisotropy of equivalent permeability are analyzed based on coordinate rotations of the permeability tensor ellipse [39]. The directional permeability for grid blocks of equivalent fracture models is plotted in Figure 10. It should be noted that the permeability tensor ellipses are normalized, which only indicates the direction and anisotropy of equivalent permeability for a grid block rather than the magnitude of equivalent permeability for the grid blocks. It is shown that with the correlation exponent changes from 0.5 to 1.0, the ratio of (the major axis of the ellipse) to (the minor axis of the ellipse) generally increased for the grid blocks, which suggest that the anisotropy of equivalent permeability increased.

Compared to the constant aperture model, the correlation of fracture length-aperture increases the heterogeneity of the discrete fracture model as well as that of the upscaled equivalent fracture model [26, 35]. The study emphasizes the importance of aperture-length model for building discrete fracture models and further for the equivalent permeability distributions. Previous studies mainly concentrate on the comparison between constant and varied aperture models on equivalent permeability of one grid block (e.g., [35]); this study investigates the influence of correlated fracture aperture and length model on equivalent permeability distributions, which helps link discrete fractures at the small scales to equivalent permeability of reservoir models at the field scales. However, we should note that the fractures are highly influenced by the mechanical properties and stress field of surrounding rocks [59, 60], and the variation of fracture aperture is mechanical properties also be determined by constitutive models (e.g., [38]). Thus, how different aperture models affect on equivalent permeability distributions for grid blocks of equivalent fracture models should be further studied. Another limitation of the study is that the discrete fracture models are in two dimensions, which is a simplification for three-dimensional models. For investigating three-dimensional fractured porous media, realistic fracture geometric data and high computation platforms for meshing and simulation flow in discrete fracture models (e.g., [61]) are needed.

7. Conclusions

To conclude, a series of equivalent fracture models are built based on discrete fracture models with different aperture-length correlations and fracture geometries by the multiple boundary upscaling method, and the conclusions are as follows: (1)For the correlation exponent , with the increase in the minimum fracture length and the number of fractures , the diagonal components of equivalent permeability tensors change from a power law like to a lognormal like to a normal like distribution(2)With the increase in correlation exponent , and change “slower” from a power law distribution to a normal distribution, which is mainly due to the increase in heterogeneity. For , the spatial distribution keeps as normal distributions with a median of zero regarding the random orientation of fractures(3)The relationship between the average dimensionless equivalent permeability and the correlation exponent follows an exponential function. The relationships for different aperture-length models are similar, which are partially related to the fracture geometries. The dimensionless coefficients for the exponential function are influenced by fracture geometry. That is, they increase as and raise(4)The correlation between the average dimensionless equivalent permeability and dimensionless fracture density model follows a power law for the correlated aperture-fracture models. Likewise, the dimensionless coefficients for the power law model are influenced by the correlation exponent and increase as raises. The analysis demonstrates that both the fracture geometric parameters and the length-aperture model influence the equivalent permeability of equivalent fracture models interactively

Data Availability

The data that support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The author declared that he has no conflicts of interest to this work.

Acknowledgments

This work was funded by the Natural Science Foundation of Shandong Province, China (ZR2019BD028 and ZR2019MD013), the Scientific Research Foundation of Shandong University of Science and Technology for Recruited Talents (2019RCJJ004), and the National Natural Science Foundation of China (41702305). The author would also like to acknowledge the developers for the open-source codes MRST and ADFNE applied in the study.