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

Authors: Tao Chen
Journal: Geofluids
 (2020)

---

## 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.


---

## Body

## 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 [1–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 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 [21–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–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.


## 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]
(1)w=γlD,wherew is fracture aperture, l is fracture length, γ is the coefficient related to mechanical properties of fractured rock, 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].Figure 1
Theoretical scaling of fracture aperture with length for geologic structural discontinuities [45].Fracture length covers a wide scale and can be described by a power law model [41]:
(2)nl=Al−adl,where l is the fracture length, nl 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 sizeL=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 lmin and the upper bound lmax are 4 m and 40 m, respectively. The number of fractures N is 50. The fracture geometries vary by increasing lmin 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 permeability of fractures and rock matrix are taken into account, which constitutes discrete fracture models.Figure 2
Discrete fracture models with correlated fracture aperture and length.


### 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 is2m×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 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 thex-axis (Figure 3). Then, the steady flow problem is solved, and the fluxes qx and qy from both the fracture and rock matrix are calculated with a multiple boundary expression [40]. Lastly, the equivalent permeability component kxx and kyx 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. kxy and kyy 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 off-diagonal components.Figure 3
Calculating flow rates on multiple boundaries for a grid block under linear boundary conditions along thex-axis (left) and y-axis (right).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 km=9.87×10−16m2(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, kf=w2/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.Figure 4
Meshes for the discrete fracture model, in which the bold lines denote fracture grids (left), and the corresponding equivalent permeability componentkxx for different correlation exponents D (right).


## 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]
(1)w=γlD,wherew is fracture aperture, l is fracture length, γ is the coefficient related to mechanical properties of fractured rock, 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].Figure 1
Theoretical scaling of fracture aperture with length for geologic structural discontinuities [45].Fracture length covers a wide scale and can be described by a power law model [41]:
(2)nl=Al−adl,where l is the fracture length, nl 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 sizeL=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 lmin and the upper bound lmax are 4 m and 40 m, respectively. The number of fractures N is 50. The fracture geometries vary by increasing lmin 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 permeability of fractures and rock matrix are taken into account, which constitutes discrete fracture models.Figure 2
Discrete fracture models with correlated fracture aperture and length.


## 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 is2m×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 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 thex-axis (Figure 3). Then, the steady flow problem is solved, and the fluxes qx and qy from both the fracture and rock matrix are calculated with a multiple boundary expression [40]. Lastly, the equivalent permeability component kxx and kyx 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. kxy and kyy 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 off-diagonal components.Figure 3
Calculating flow rates on multiple boundaries for a grid block under linear boundary conditions along thex-axis (left) and y-axis (right).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 km=9.87×10−16m2(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, kf=w2/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.Figure 4
Meshes for the discrete fracture model, in which the bold lines denote fracture grids (left), and the corresponding equivalent permeability componentkxx for different correlation exponents D (right).


## 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 exponentD 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 D, 16 discrete fracture networks with varied fracture network geometries, i.e., the minimum fracture length lmin 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.Figure5(a) shows the Cartesian grid of the equivalent fracture model for a realization of the discrete fracture models with lmin=4,N=125, and D=1. The spatial distributions of equivalent permeability components kxx, kyy, and kxy are plotted in Figures 5(b)–5(d), respectively. The equivalent permeability ranges from 9.87×10−16m2 (matrix permeability) to 3.5×10−9m2. For the diagonal components kxx and kyy, 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 kxx and kyy is different from each other due to the random orientation of the fractures. For the off-diagonal component kxy in Figure 5(d), they are either positive or negative, and the grid blocks with high absolute values correspond to those of large kxx and kyy in Figures 5(b) and 5(c).Figure 5
(a) Cartesian grid for the discrete fracture model, (b–d) spatial distributions, and (e–g) histograms for equivalent permeability components based on upscaling.For investigating the statistical distribution of equivalent permeability tensors, the histograms ofkxx, kyy, and kxy are plotted in Figures 5(e)–5(g), respectively. The fitting curves for histograms are also plotted. It is shown that both kxx and kyy exhibit a power law-like distribution whereas kxy 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 forkxx, kyy, and kxy 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 kxx and kyy are similar and tend to a power law distribution, whereas kxy tends to a normal distribution. With the increase in lmin and N, kxx and kyy 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 kxy 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 kxx and kyy are similar, with an average of about -0.75, which indicates a relatively strong negative correlation. For kxy, the coefficient is less than -0.04, which almost denotes no correlation. It is reasonable as the off-diagonal component kxy tends to be evenly distributed around zero due to the random fracture orientations.Figure 6
Fitting curves of histograms for equivalent permeability components with different aperture-length correlation exponentD and increasing lmin and N.With the increase in the correlation exponentD, 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 lmin 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 lmin 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 lmin 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.


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

Before studying the effect of the length-aperture correlation exponentD on the equivalent permeability for the equivalent fracture models, the average dimensionless permeability k’¯ for each discrete fracture model is defined as
(3)k’¯=∑n=1Nmki/kmNm,whereNm is the total number of grid blocks for the equivalent fracture model, and the rock matrix permeability km 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.Figure 7
The average dimensionless permeability components of the equivalent fracture model versus aperture-length correlation exponentD of the discrete fracture model for increasing fracture geometric properties lmin and N.
(a)(b)(c)It shows in Figures7(a) and 7(b) that for the diagonal components of equivalent permeability tensors, logkxx’¯ and logkyy’¯, increase linearly with the increase in the correlation exponent D, which can be expressed as the following exponential function:
(4)k’¯=A∙10B∙D,where the dimensionless coefficientsA and B range from 101.4 to 102.4 and from 3.4 to 4.3, respectively. For the off-diagonal component kxy, as the negative values exist, the absolute value kxy is used to analyze the correlation between the average dimensionless permeability and the aperture-length correlation exponent D (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 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 betweenk’¯ of the equivalent fracture models and D of the discrete fracture models can be described by equation (4), the dimensionless coefficient B for kxx and kxy increases as the fracture geometric parameters lmin 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 kxyare of the same magnitude as those for kxx and kxy. However, for the coefficient A, the range is wider than those of kxx and kxy, 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).


## 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]
(5)ρ=1A∑iNli22,whereA is the domain area, li 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:
(6)k’¯=β∙ρC,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 exponentD.whereβ and C 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 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 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).Figure 9
Variations of the dimensionless coefficient (a)C and (b) β regarding the correlation exponent D. The blue color, red color, and green color denote kxx, kyy, and kxy, respectively.
(a)(b)


## 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 exponentD 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 D changes from 0.5 to 1.0, the ratio of kmax (the major axis of the ellipse) to kmin (the minor axis of the ellipse) generally increased for the grid blocks, which suggest that the anisotropy of equivalent permeability increased.Figure 10
Equivalent permeability ellipse for correlated fracture aperture and length models forD =0.5 and D =1.0.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 exponentD=0.5, with the increase in the minimum fracture length lmin and the number of fractures N, 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 exponentD, kxx and kyychange “slower” from a power law distribution to a normal distribution, which is mainly due to the increase in heterogeneity. For kxy, 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 permeabilityk’¯ 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 lmin and N raise(4)
The correlation between the average dimensionless equivalent permeabilityk’¯ 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



---

*Source: 8834666-2020-10-10.xml*