Equivalent Permeability Distribution for Fractured Porous Rocks: Correlating Fracture Aperture and Length

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.


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 [1][2][3][4][5][6]. 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 gener-ally 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 [21][22][23]. 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., [29][30][31]), 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.

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 w is fracture aperture, l is fracture length, γ is the coefficient related to mechanical properties of fractured rock, 2 Geofluids and the correlation exponent D 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., D=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, D is around 0.5 [44,47]. The correlation exponent D is greater than 0.5, and less than 1.0 could result from post-jointing relaxation [27]. In this study, the correlation exponent D ranges from 0.5 to 1.0, and γ is assumed to be 1:2 × 10 −4 according to the field data [48]. Fracture length covers a wide scale and can be described by a power law model [41]: where l is the fracture length, nðlÞ is the number of fractures with sizes in the range ½l, l + dl, A is a density constant, and a 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 a represents the growth properties of the fractures and varies from 1.3 to 3.5 [41,49], and in this study, a 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 L = 20 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 l min and the upper bound l max are 4 m and 40 m, respectively. The number of fractures N is 50. The fracture geometries vary by increasing l min from 4 m to 6 m, 8 m, and 10 m and by rising N 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 per-meability of fractures and rock matrix are taken into account, which constitutes discrete fracture models.

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 2 m × 2m. 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 twodimensional 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 x-axis ( Figure 3). Then, the steady flow problem is solved, and the fluxes q x and q y from both the fracture and rock matrix are calculated with a multiple boundary expression [40]. Lastly, the equivalent permeability component k xx and k yx 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. k xy and k yy can be calculated by changing the direction of the linear boundary condition along the y-axis. It should be noted that the upscaled equivalent permeability is no inherently symmetric, and the symmetric permeability tensor is acquired by averaging the offdiagonal 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 k m = 9:87 × 10 −16 m 2 (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, k f =w 2 /12. The numerical upscaling procedure is validated against analytical solutions. Considering a fully penetrated fracture with different correlation exponent D, 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 D.

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 D varied from 0.5 to 1.0 with a step of 0.1, which indicates  4 Geofluids fracture mechanical state changes from constant fracture toughness to constant fracture driving force. For a specific correlation exponent D, 16 discrete fracture networks with varied fracture network geometries, i.e., the minimum fracture length l min and the number of fractures N, 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 l min = 4,N = 125, and D = 1. The spatial distributions of equivalent permeability components k xx , k yy , and k xy are plotted in Figures 5(b)-5(d), respectively. The equivalent permeability ranges from 9:87 × 10 −16 m 2 (matrix permeability) to 3:5 × 10 −9 m 2 . For the diagonal components k xx and k yy , 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 k xx and k yy is different from each other due to the random orientation of the fractures. For the offdiagonal component k xy in Figure 5(d), they are either positive or negative, and the grid blocks with high absolute values correspond to those of large k xx and k yy in Figures 5(b) and 5(c).
For investigating the statistical distribution of equivalent permeability tensors, the histograms of k xx , k yy , and k xy are plotted in Figures 5(e)-5(g), respectively. The fitting curves for histograms are also plotted. It is shown that both k xx and k yy exhibit a power law-like distribution whereas k xy 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 k xx , k yy , and k xy 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 D = 0:5, it shows that the shape and the magnitude of the diagonal components k xx and k yy are similar and tend to a power law distribution, whereas k xy tends to a normal distribution. With the increase in l min and N, k xx and k yy 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 k xy 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 k xx and k yy are similar, with an average of about -0.75, which indicates a relatively strong negative correlation. For k xy , the coefficient is less than -0.04, which almost denotes no correlation. It is reasonable as the offdiagonal component k xy tends to be evenly distributed around zero due to the random fracture orientations.
With the increase in the correlation exponent D, 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 l min and N, the evolution of the histogram fitting curves differs for varied correlation exponent D. For a higher D (e.g., D = 1), the histograms do not show a clear normal distribution for large l min and N compared to those for a lower D, e.g., D = 0:5: They look like a transition shape for a lower D changing from a power law to a lognormal to a normal distribution. The difference is partially due to that the high correlation exponent D increases the heterogeneity of the discrete fracture model. And it is reasonable to speculate that with the increase in l min and N, the diagonal components tend to a lognormal distribution or normal distribution for a high D. Accordingly, for a specific fracture geometry, a high correlation exponent D means a higher heterogeneity which requires a sufficient number of fractures to reach a lognormal like distribution.

Correlating the Equivalent Permeability with the Aperture-Length Model
Before studying the effect of the length-aperture correlation exponent D on the equivalent permeability for the equivalent fracture models, the average dimensionless permeability k' for each discrete fracture model is defined as where N m is the total number of grid blocks for the equivalent fracture model, and the rock matrix permeability k m is equal to 1 md as assumed in the previous section. Under a specific aperture-length correlation exponent D, 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 k' and the aperture-length correlation exponent D of the discrete fracture models.
It shows in Figures 7(a) and 7(b) that for the diagonal components of equivalent permeability tensors, logð k ' xx Þ and logð k ' yy Þ, increase linearly with the increase in the correlation exponent D, which can be expressed as the following exponential function: where the dimensionless coefficients A and B range from 10 1.4 to 10 2.4 and from 3.4 to 4.3, respectively. For the off-diagonal component k xy , as the negative values exist, the absolute value jk xy j is used to analyze the correlation between the average dimensionless permeability and the aperture-length correlation exponent D (Figure 7(c)).

Geofluids
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 k' and the aperture-length correlation exponent D. 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 k' of the equivalent fracture models and D of the discrete fracture models can be described by equation (4), the dimensionless coefficient B for k xx and k xy increases as the fracture geometric parameters l min and N raise (Figures 7(a) and  7(b)). In contrast, the dimensionless coefficient A 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 A and B for jk xy j are of the same magnitude as those for k xx and k xy . However, for the coefficient A, the range is wider than those of k xx and k xy , which may be mainly due to the process of averaging the off-diagonal components. For the coefficient B, 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).

Effects of Fracture Geometry on the Equivalent Permeability
The influence of fracture geometries on the average dimensionless equivalent permeability for correlated lengthaperture models is analyzed in this section. For each discrete fracture model, the dimensionless fracture density is defined as [28] where A is the domain area, l i is the length of the i-th fracture, and N is the total number of fractures of the discrete fracture model.
The correlations between the log( k') 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( k') increases linearly with the increase in logðρÞ for a given length-aperture correlation exponent D, which can be fitted by the following power law relationship: where β and C are dimensionless coefficient in the ranges 10 3 -10 5 and 1.1-1.3, respectively. Whereas for the offdiagonal 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 D increases. This is mainly because a higher correlation exponent D results in an enlarged difference in equivalent permeability. It indicates further that the higher D in the correlated length-aperture model, the higher heterogeneity for the equivalent fracture permeability field. It should be noted that C closing to 1 means a linear correlation between k' and ρ, which has been observed in the results for macroscopic hydraulic conductivity by Leung and Zimmerman [28]. Furthermore, when the 6 Geofluids correlation exponent D increases, β and C also increase ( Figure 9). It shows that the aperture-length model influences the dimensionless coefficients β and C for the power model described by equation (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 7 Geofluids 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 D 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  Figure 8: The average dimensionless permeability components of the equivalent fracture models versus dimensionless density of the discrete fracture models for varied aperture-length correlation exponent D. 8 Geofluids 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 D changes from 0.5 to 1.0, the ratio of k max (the major axis of the ellipse) to k min (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 frac-ture 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 threedimensional 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.

Conclusions
To conclude, a series of equivalent fracture models are built based on discrete fracture models with different aperturelength correlations and fracture geometries by the multiple boundary upscaling method, and the conclusions are as follows: (1) For the correlation exponent D = 0:5, with the increase in the minimum fracture length l min and the number of fractures N, the diagonal components of equivalent permeability tensors change from a 9 Geofluids power law like to a lognormal like to a normal like distribution (2) With the increase in correlation exponent D, k xx and k yy change "slower" from a power law distribution to a normal distribution, which is mainly due to the increase in heterogeneity. For k xy , 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 k' and the correlation exponent D 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 l min and N raise (4) The correlation between the average dimensionless equivalent permeability k' 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 D and increase as D 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.