# Equivalent Permeability Distribution for Fractured Porous Rocks: The Influence of Fracture Network Properties

Authors: Tao Chen
Journal: Geofluids
 (2020)

---

## Abstract

Equivalent fracture models are widely used for simulations of groundwater exploitation, geothermal reservoir production, and solute transport in groundwater systems. Equivalent permeability has a great impact on such processes. In this study, equivalent permeability distributions are investigated based on a state-of-the-art numerical upscaling method (i.e., the multiple boundary method) for fractured porous rocks. An ensemble of discrete fracture models is created based on power law length distributions. The equivalent permeability is upscaled from discrete fracture models based on the multiple boundary method. The results show that the statistical distributions of equivalent permeability tensor components are highly related to fracture geometry and differ from each other. For the histograms of the equivalent permeability, the shapes ofkxx and kyy change from a power law-like distribution to a lognormal-like distribution when the fracture length and the number of fractures increase. For the off-diagonal component kxy, it is a normal-like distribution and its range expands when the fracture length and the number of fractures increase. The mean of diagonal equivalent permeability tensor components increases linearly with the fracture density. The analysis helps in generating stochastic equivalent permeability models in fractured porous rocks and reduces uncertainties when applying equivalent fracture models.


---

## Body

## 1. Introduction

Modeling flow and its coupled processes in fractured porous rocks are important for subsurface environmental and energy problems [1, 2], such as solute transport in groundwater or at radioactive waste disposal sites [3–5], geothermal energy exploitation [6, 7], gas hydrate stability [8], and coal mine water inflows [9]. The numerical models can be roughly divided into two categories: equivalent fracture models and discrete fracture models [10–12]. For discrete fracture models, fracture geometry is depicted either deterministically or stochastically (e.g., [13]); however, the simulations at the field scale are too computationally expensive to applications when the fracture density is large. For equivalent fracture models, permeability is upscaled from discrete fracture geometry to grid blocks, which reduces the complexity of simulations and makes the equivalent fracture models applicable at the field scale (e.g., [14]).For equivalent fracture models, a major part of uncertainty comes from the characterizations of permeability fields (e.g., [15]). Upscaling permeability for heterogeneous porous rocks has been comprehensively reviewed (e.g., [16, 17]). Nevertheless, the permeability of fractured porous rocks is relatively complex as it is highly heterogeneous and anisotropic and varies with the scale of measurement (e.g., [18–21]). Upscaling permeability from multiple-scale discrete fractures to equivalent fracture models was widely studied in the last decades, which contributes to a better knowledge of permeability underground as well as reduces the uncertainty of equivalent fracture models [2, 22, 23]. Bogdanov et al. [24] calculated equivalent permeability numerically by solving flow problems in three-dimensional models consisting of both discrete fractures and rock matrix. Öhman and Niemi [25] compared two upscaling methods for fractured media, a classical fracture network-based approach and a new enhanced stochastic continuum approach. Frampton and Cvetkovic [26] studied reactive tracer transport through discrete fracture networks using a stochastic Lagrangian framework, combined with the methodology for upscaling particle arrival times. Gong et al. [27] developed an upscaling method to build a dual-porosity, dual-permeability model from discrete fracture characterizations. Ebigbo et al. [28] studied the inclusion-based effective medium models to calculate the permeability of three-dimensional fractured porous rocks. Recently, Rajeh et al. [29] developed a modified semianalytical superposition method considering the hydraulic conductivity of the porous matrix and the fractures, as well as the connectivity of the conductive fracture network. Sweeney et al. [30] presented a computational geometry-based upscaling approach to accurately capture the dynamic processes considering the fracture geometric properties and the matrix properties. We should note that flow can also occur between disconnected fractures if the rock mass has some porosity and permeability, which would mean some differences in calculations of equivalent permeability [31, 32]. Furthermore, Laubach et al. [33] reviewed the diagenetic or chemical processes in fractured porous rocks, which tend to disconnect fracture trace patterns.Specifically, the relationship between equivalent permeability and fracture geometry has been extensively investigated based on either analytical or numerical upscaling methods [34–36]. de Dreuzy et al. [37] analyzed the influences of fracture length distribution on equivalent permeability for two-dimensional discrete fracture networks in which the matrix is assumed as impervious. Baghbanan and Jing [38] studied the equivalent permeability and the representative elementary volume (REV) of two-dimensional fractured rocks considering the correlation between fracture aperture and fracture length. Lang et al. [39] developed a numerical method to compute the equivalent permeability for three-dimensional discrete fracture models based on the volume averaging of pressure and flux. For discrete fracture models with power law length distribution and with aperture-length correlation, two-dimensional cut planes could underestimate the equivalent permeability by up to three orders of magnitude compared to that for three-dimensional models. Li et al. [40] established analytical expressions for equivalent permeability and fractal aperture distribution based on a multiple fractal modeling considering both the rock matrix and discrete fracture network (i.e., the discrete fracture model). Maillot et al. [41] calculated equivalent permeability for three-dimensional discrete fracture networks with a self-similar power law fracture length distribution and compared the difference between Poisson models in which the fractures are distributed randomly and kinematic fracture model which mimics dynamic fracturing processes. Hardebol et al. [42] integrated multiple-scale datasets for generating discrete fracture models. Their simulation results show that the equivalent permeability of fractured porous rocks can be two to three orders of magnitude higher than the matrix permeability, which is highly controlled by the connectivity of the fractures. Recently, Chen et al. [43] calculated the equivalent permeability distributions for a well-connected model and a poorly connected discrete fracture model and analyzed equivalent permeability distribution statistically based on histograms. However, only a few studies investigate the fracture geometry to equivalent permeability distributions for an equivalent fracture model, which has important implications for creating stochastic permeability fields for fractured porous rocks.In this study, the influences of fracture network geometry on equivalent permeability distributions in fractured porous media are systematically quantified. The originality of the work is studying the correlation between the equivalent permeability distributions and the fracture geometry by using the novel upscaling method, namely, the multiple boundary method [44] for fractured porous rocks. This paper is structured as follows: first, the geometrical properties of an ensemble of synthetic two-dimensional discrete fracture models are described. Then, the equivalent permeability for the discrete fracture models is upscaled to equivalent fracture models by using the multiple boundary method. Lastly, the statistical distribution of the equivalent permeability and the correlations between equivalent permeability and the fracture density are analyzed, which quantitatively illustrates how equivalent permeability distributions are influenced by the fracture network geometry.


## 2. Generating and Upscaling for Discrete Fracture Models

For investigating the equivalent permeability distributions in fractured porous rocks, two main procedures are involved. First, discrete fracture models are created stochastically based on fracture geometric parameters. Then, the equivalent permeability fields are calculated from the discrete fracture models with the multiple boundary method.


### 2.1. Stochastic Discrete Fracture Models

A discrete fracture model, which includes fractures and rock matrix, depicts fracture geometry explicitly. The fracture geometry mainly includes fracture length, fracture orientation, fracture position, and fracture density. According to multiple-scale fracture characterization data, the fracture length can be described by a power law distribution in natural fracture systems [45]:
(1)nl=Al−a,where l is the fracture length, nldl is the number of fractures with sizes in the range l,l+dl, A is a constant, and a is a power law exponent. As the length of fractures is limited by the rock size in the Earth’s crust, it should have upper and lower bounds in the power law distribution. The power law exponent a represents the growth properties of the fractures and varies from 1.3 and 3.5 [45, 46].In this study, two-dimensional discrete fracture models with varied geometric properties are generated stochastically in a20m×20m domain based on the ADFNE code [47]. The fracture centers are seeded over in the domain. The fracture length follows power law distribution, and the lower bound lmin and the upper bound lmax are 4 m and 40 m, respectively. For the fractures partially within the domain, they are truncated by the boundaries. For the multiscale length fracture network, it is hard to avoid the censor effect in a given domain for a randomly distributed fracture network, which may underestimate the real fracture length for the fracture network. Nevertheless, for a power-law distribution of lengths, it is not usually a serious problem in this study [48].The lower bound of the fracture lengthlmin increases from 4 m to 6 m, 8 m, and 10 m. The number of fractures N increases from 50 to 75, 100, and 125. Accordingly, the varied fracture length and fracture density yield sixteen discrete fracture models (Figure 1). The fracture orientation and the fracture position are generated randomly. Fracture apertures range from μm to cm and depend on the scale of measurement [45]. As investigating the effect of fracture network properties on equivalent permeability distributions is the primary goal of this study, a constant fracture aperture of 1.2×10−4m is assumed to ensure that all fractures influence fluid flow and to get rid of other affecting factors [42, 49]. For the domain of 20m×20m, the aperture value in this study is comparable to those in Baghbanan and Jing [38], Lei et al. [50], and Agheshlui et al. [51]. However, we should note that for naturally fractured porous rocks, the rock stress and rock-fluid interaction may result in varied apertures of fractures [18], which will be elaborated further in the discussion section. The fracture permeability is 1.2×10−9m2 according to the cubic law [52]. The matrix permeability is assumed as km=9.87×10−16m2(1 md).Figure 1
Fracture network geometry with different geometric parameters.


### 2.2. Equivalent Permeability Upscaling

The upscaling methods for fractured porous rocks can be roughly classified as analytical and numerical methods [22]. The analytical method calculates equivalent permeability efficiently by using an analytical function; however, its uncertainty increases as the fracture geometry tends to multiscale and complex. For the numerical method, after dividing the whole fracture network into grid blocks of an equivalent fracture model (Figure 2(a)), it requires to mesh the discrete fracture model (Figure 2(b)) and solve steady flow problems numerically with specified boundary conditions for each grid block (Figures 2(c) and 2(d)). Based on the results of the solution, the equivalent permeability is calculated inversely according to Darcy’s law. The procedures of the numerical method are more complex than that of the analytical method as the meshing code, and the flow simulation code are required, whereas the accuracy will be improved regarding the complex fracture network geometry (e.g., [53]).Figure 2
(a) Unstructured meshes of the discrete fracture model within the grid block. (b) Cartesian grid for the equivalent fracture model withlmin=4m and N=50; the corresponding solutions of the steady-state flow problems with the linear boundary conditions (c) along the x-axis and (d) along the y-axis.In this study, the multiple boundary method is applied [44] for calculating the equivalent permeability. Inspired by the single boundary method (e.g., [54]), the multiple boundary method utilizes flow rate information on multiple boundaries for calculating equivalent permeability. For example, when the linear boundary conditions are applied along the x-axis, the flow rates along the x- and y-axis, qx and qy, are calculated as follows:
(2)qx=∫0lyvr∙nxdy+∫0lxvu∙nxdx+∫0lxvl∙nxdx,qy=∫0lxvu∙nydx+∫0lxvl∙nydx+∫0lyvr∙nydy,where vr, vu, and vl are the specific discharges on the right, upper, and lower boundaries; nx and ny are unit vectors along the x- and y-axes; and lx and ly are dimensions of the grid block in the x- and y-directions. Then, the equivalent permeability components can be calculated based on the matrix form of Darcy’s law:
(3)−1μkxxkxykyxkyy∇Px∇Py=qxqy,where μ is the dynamic viscosity and ∇Px and ∇Py are pressure gradients along the x-axis and y-axis, respectively. In the case of applying linear boundary conditions along thex-axis, ∇Py is equal to zero and kxx and kyx can be calculated. By using the linear boundary condition along the y-axis, ∇Px equals zero and the other equivalent permeability components kxy and kyy can be calculated in the same way. As the calculated equivalent permeability components kxy and kyx are not inherently consistent, the symmetric permeability tensor can be obtained by averaging the off-diagonal components.It should be noted that the flow boundary conditions affect the resulting equivalent permeability. They can be mainly classified into three categories: no-flow, linear, and periodic boundary conditions (e.g., [27, 49, 54, 55]). Regarding the multiple upscaling method, the linear boundary condition is applied as it mimics the realistic flow condition underground to some extent and the resulting equivalent permeability matches well with the analytical solutions [44]. By contrast, applying no-flow boundary conditions or the single boundary method, which calculates the flow rate qx or qy only on the right or the top boundary (see Equation (2)), may underestimate the flow rate and therefore result in a smaller equivalent permeability than that of the analytical solution. In this study, the dimensions of the grid block lx and ly are 2 m and the pressure gradient of the linear boundary conditions is 1 Pa/m. During the upscaling procedure, the MATLAB Reservoir Simulation Toolbox (MRST) is used for meshing and solving flow problems for the discrete fracture models [56, 57].


## 2.1. Stochastic Discrete Fracture Models

A discrete fracture model, which includes fractures and rock matrix, depicts fracture geometry explicitly. The fracture geometry mainly includes fracture length, fracture orientation, fracture position, and fracture density. According to multiple-scale fracture characterization data, the fracture length can be described by a power law distribution in natural fracture systems [45]:
(1)nl=Al−a,where l is the fracture length, nldl is the number of fractures with sizes in the range l,l+dl, A is a constant, and a is a power law exponent. As the length of fractures is limited by the rock size in the Earth’s crust, it should have upper and lower bounds in the power law distribution. The power law exponent a represents the growth properties of the fractures and varies from 1.3 and 3.5 [45, 46].In this study, two-dimensional discrete fracture models with varied geometric properties are generated stochastically in a20m×20m domain based on the ADFNE code [47]. The fracture centers are seeded over in the domain. The fracture length follows power law distribution, and the lower bound lmin and the upper bound lmax are 4 m and 40 m, respectively. For the fractures partially within the domain, they are truncated by the boundaries. For the multiscale length fracture network, it is hard to avoid the censor effect in a given domain for a randomly distributed fracture network, which may underestimate the real fracture length for the fracture network. Nevertheless, for a power-law distribution of lengths, it is not usually a serious problem in this study [48].The lower bound of the fracture lengthlmin increases from 4 m to 6 m, 8 m, and 10 m. The number of fractures N increases from 50 to 75, 100, and 125. Accordingly, the varied fracture length and fracture density yield sixteen discrete fracture models (Figure 1). The fracture orientation and the fracture position are generated randomly. Fracture apertures range from μm to cm and depend on the scale of measurement [45]. As investigating the effect of fracture network properties on equivalent permeability distributions is the primary goal of this study, a constant fracture aperture of 1.2×10−4m is assumed to ensure that all fractures influence fluid flow and to get rid of other affecting factors [42, 49]. For the domain of 20m×20m, the aperture value in this study is comparable to those in Baghbanan and Jing [38], Lei et al. [50], and Agheshlui et al. [51]. However, we should note that for naturally fractured porous rocks, the rock stress and rock-fluid interaction may result in varied apertures of fractures [18], which will be elaborated further in the discussion section. The fracture permeability is 1.2×10−9m2 according to the cubic law [52]. The matrix permeability is assumed as km=9.87×10−16m2(1 md).Figure 1
Fracture network geometry with different geometric parameters.


## 2.2. Equivalent Permeability Upscaling

The upscaling methods for fractured porous rocks can be roughly classified as analytical and numerical methods [22]. The analytical method calculates equivalent permeability efficiently by using an analytical function; however, its uncertainty increases as the fracture geometry tends to multiscale and complex. For the numerical method, after dividing the whole fracture network into grid blocks of an equivalent fracture model (Figure 2(a)), it requires to mesh the discrete fracture model (Figure 2(b)) and solve steady flow problems numerically with specified boundary conditions for each grid block (Figures 2(c) and 2(d)). Based on the results of the solution, the equivalent permeability is calculated inversely according to Darcy’s law. The procedures of the numerical method are more complex than that of the analytical method as the meshing code, and the flow simulation code are required, whereas the accuracy will be improved regarding the complex fracture network geometry (e.g., [53]).Figure 2
(a) Unstructured meshes of the discrete fracture model within the grid block. (b) Cartesian grid for the equivalent fracture model withlmin=4m and N=50; the corresponding solutions of the steady-state flow problems with the linear boundary conditions (c) along the x-axis and (d) along the y-axis.In this study, the multiple boundary method is applied [44] for calculating the equivalent permeability. Inspired by the single boundary method (e.g., [54]), the multiple boundary method utilizes flow rate information on multiple boundaries for calculating equivalent permeability. For example, when the linear boundary conditions are applied along the x-axis, the flow rates along the x- and y-axis, qx and qy, are calculated as follows:
(2)qx=∫0lyvr∙nxdy+∫0lxvu∙nxdx+∫0lxvl∙nxdx,qy=∫0lxvu∙nydx+∫0lxvl∙nydx+∫0lyvr∙nydy,where vr, vu, and vl are the specific discharges on the right, upper, and lower boundaries; nx and ny are unit vectors along the x- and y-axes; and lx and ly are dimensions of the grid block in the x- and y-directions. Then, the equivalent permeability components can be calculated based on the matrix form of Darcy’s law:
(3)−1μkxxkxykyxkyy∇Px∇Py=qxqy,where μ is the dynamic viscosity and ∇Px and ∇Py are pressure gradients along the x-axis and y-axis, respectively. In the case of applying linear boundary conditions along thex-axis, ∇Py is equal to zero and kxx and kyx can be calculated. By using the linear boundary condition along the y-axis, ∇Px equals zero and the other equivalent permeability components kxy and kyy can be calculated in the same way. As the calculated equivalent permeability components kxy and kyx are not inherently consistent, the symmetric permeability tensor can be obtained by averaging the off-diagonal components.It should be noted that the flow boundary conditions affect the resulting equivalent permeability. They can be mainly classified into three categories: no-flow, linear, and periodic boundary conditions (e.g., [27, 49, 54, 55]). Regarding the multiple upscaling method, the linear boundary condition is applied as it mimics the realistic flow condition underground to some extent and the resulting equivalent permeability matches well with the analytical solutions [44]. By contrast, applying no-flow boundary conditions or the single boundary method, which calculates the flow rate qx or qy only on the right or the top boundary (see Equation (2)), may underestimate the flow rate and therefore result in a smaller equivalent permeability than that of the analytical solution. In this study, the dimensions of the grid block lx and ly are 2 m and the pressure gradient of the linear boundary conditions is 1 Pa/m. During the upscaling procedure, the MATLAB Reservoir Simulation Toolbox (MRST) is used for meshing and solving flow problems for the discrete fracture models [56, 57].


## 3. Results

For the discrete fracture models with fractures shown in Figure1, the equivalent permeability is upscaled from the discrete fracture models to the equivalent fracture models. For yielding generalized results for the statistical distributions of the equivalent permeability, ten realizations are created for a given set of fracture geometric parameters. The equivalent permeability tensor consists of four components, kxx, kxy, kyx, and kyy, for two-dimensional models as shown in Equation (3). Herein, the equivalent permeability tensor is assumed to be symmetric, i.e., kxy=kyx. The statistical distributions of the equivalent permeability and the correlation between the equivalent permeability and the fracture geometric parameter are analyzed in the following.


### 3.1. Statistical Distributions of the Equivalent Permeability

The Cartesian grid of the equivalent fracture models and the upscaled equivalent permeability field forkxx, kxy, and kyy and for lmin=4m and N=50 are plotted in Figure 3. The diagonal components of equivalent permeability tensor kxx and kyy are of the same order of magnitude and less than 3×10−13 m2. However, the range of kxy lies between −2×10−13 and 2×10−13 m2. This is because the diagonal components kxx and kyy should always be positive according to the physical meaning of the permeability in Darcy’s law, whereas the off-diagonal term kxy could be positive or negative which depends on the angle between the axes of the equivalent permeability tensor ellipse and the coordinate axes. The range of absolute value of kxy is smaller than those of kxx and kyy as the equivalent permeability tensor should be positive definite matrices [58]. For the grid blocks where the fractures have a large density and their orientations tend to along the x-axis, kxx is bigger. This is also similar for kyy. It should be noted that although the range of kxx and kyy is close, their spatial distributions are different (Figure 3) which is primarily dominated by fracture orientations.Figure 3
The spatial distribution (top) and the histograms (bottom) of the equivalent permeability tensor componentskxx, kyy, and kxy for the discrete fracture model corresponding to that in Figure 2(b).The histograms of equivalent permeability for the fracture networks are plotted in Figure3. The ranges of kxx, kyy, and kxy are divided equally into 10 bins and are gathered on the same frame. It shows that the histograms of kxx, kyy, and kxy exhibit different shapes. For kxy, they tend to be symmetrically distributed with a middle value around 0 m2. For kxx and kyy, they are all positive and exhibit a power law-like distribution. The histograms for kxx, kyy, and kxy are also fitted by the solid lines with the least square method (Figure 3).The fitted lines of the histograms ofkxx, kyy, and kxy for all realizations are gathered on the same frame (Figure 4). The statistical properties of equivalent permeability also change with the fracture geometric parameters. For the discrete fracture model with lmin=4 m and N=50, the modes of kxx and kyy are less than 1×10−13 m2 and kxy is distributed roughly between −1×10−13 and 1×10−13 m2. The shapes for kxx and kyy tend to a power law shape which are similar to that of the upscaled equivalent permeability for a three-dimensional poorly connected discrete fracture model [43]. This is primarily due to the fact that a lot of grid blocks are dominated by the matrix permeability where the fractures inside do not penetrate the grid blocks or below the percolation threshold. When lmin or N increase, the shapes for kxx and kyy change: the modes migrate from a low value to a high value and the range of kxy expands symmetrically. It should be noted that the shapes of the histograms for kxx and kyy tend to lognormal-like distributions whereas they keep as a normal-like distribution for kxy during the increase of lmin and N. The Gaussian distribution or log-normal distribution for equivalent permeability is commonly assumed in equivalent fracture models (e.g., [15]) and is supported by field measurement data (e.g., [59]). The results in Figure 4 support this assumption, which are based on upscaling permeability for fractured porous rocks. In particular, the fitting curves in Figure 4 further show how the shapes of the histograms of equivalent permeability migrated from power law like to lognormal-like distributions with the change of fracture geometries. The results highlight the importance of fracture geometries on equivalent permeability distributions.Figure 4
Fitted curves for the histograms of equivalent permeability regarding varied fracture geometric parameterslmin and N for all realizations. The blue, red, and green colors denote kxx, kyy, and kxy, respectively.


### 3.2. Correlation between Equivalent Permeability Distribution and Fracture Network Geometry

For analyzing the correlations between the equivalent permeability and the fracture network geometry, the dimensionless equivalent permeabilityk’ of the grid block and the dimensionless fracture density ρ for the discrete fracture model are defined [60], respectively, as follows:
(4)k′=kkm,(5)ρ=1A∑iNli22,where k denotes the components of the equivalent permeability, li denotes the fracture length of the i-th fracture, and A is the area of the discrete fracture model which is 400 m2 in this model. For each discrete fracture model, the mean and standard deviation of the dimensionless equivalent permeability components k′xx¯, k′yy¯, and k′xy¯, and σk′xx, σk′yy, and σk′xy as well as the dimensionless fracture density ρ are calculated. The statistical properties of equivalent permeability tensor components for all 160 discrete fracture models are plotted in Figure 5.Figure 5
Plots of the mean and standard deviation of the dimensionless equivalent permeability tensor components with the dimensionless fracture densityρ for varied fracture geometric parameters.
(a)(b)(c)(d)(e)(f)It is shown that a strong correlation exists between the diagonal components of equivalent permeability tensor and the fracture network geometry:k′xx¯ and k′yy¯ increase with the dimensionless fracture density. The standard deviation for k′xx and k′yy also increases with ρ which is reasonable as few grid blocks may be interconnected by several fractures. The correlation between k′xx¯ and k′yy¯ and ρ can be built by the fitted line:
(6)k′diag¯=a∙ρ+b,where ρis in the range of 2-14, k′diag¯ is the diagonal components of the equivalent permeability tensor which is in the range of 50-350, a is the slope about 20, and b is the constant about 30 here. Leung and Zimmerman [60] investigated the equivalent permeability of a single grid block, and the results show that the equivalent permeability increases linearly with the fracture density. Similar results are also observed for the mean equivalent permeability of grid blocks, as shown in Figure 5. The slope coefficient of a in this study is comparable for those in a single grid block. Although the main scope of the paper concentrates on the influence of fracture network geometries on equivalent permeability distribution, the factor of varied fracture aperture, which makes the model complex, needs to be considered in further studies.Forkxy, the mean k′xy¯ shows a small fluctuation around zero with the changes of ρ. Furthermore, for kxy, its range varies when the fracture network properties change (Figure 4). When ρincreases, the ranges for k′xy¯ and σk′xy also increase. This is mainly because for a grid block with an inclined fracture with azimuth θ, the absolute value of kxy increases with kxx and kyy, which is given by [61] the following:
(7)kθ=kxxkxykyxkyy=kx∗cos2θ+ky∗sin2θkx∗−ky∗cosθ∙sinθkx∗−ky∗cosθ∙sinθkx∗sin2θ+ky∗cos2θ,where kx∗ and ky∗ are the components of the equivalent permeability tensor in x- and y-axes when the fracture azimuth is 0°.


## 3.1. Statistical Distributions of the Equivalent Permeability

The Cartesian grid of the equivalent fracture models and the upscaled equivalent permeability field forkxx, kxy, and kyy and for lmin=4m and N=50 are plotted in Figure 3. The diagonal components of equivalent permeability tensor kxx and kyy are of the same order of magnitude and less than 3×10−13 m2. However, the range of kxy lies between −2×10−13 and 2×10−13 m2. This is because the diagonal components kxx and kyy should always be positive according to the physical meaning of the permeability in Darcy’s law, whereas the off-diagonal term kxy could be positive or negative which depends on the angle between the axes of the equivalent permeability tensor ellipse and the coordinate axes. The range of absolute value of kxy is smaller than those of kxx and kyy as the equivalent permeability tensor should be positive definite matrices [58]. For the grid blocks where the fractures have a large density and their orientations tend to along the x-axis, kxx is bigger. This is also similar for kyy. It should be noted that although the range of kxx and kyy is close, their spatial distributions are different (Figure 3) which is primarily dominated by fracture orientations.Figure 3
The spatial distribution (top) and the histograms (bottom) of the equivalent permeability tensor componentskxx, kyy, and kxy for the discrete fracture model corresponding to that in Figure 2(b).The histograms of equivalent permeability for the fracture networks are plotted in Figure3. The ranges of kxx, kyy, and kxy are divided equally into 10 bins and are gathered on the same frame. It shows that the histograms of kxx, kyy, and kxy exhibit different shapes. For kxy, they tend to be symmetrically distributed with a middle value around 0 m2. For kxx and kyy, they are all positive and exhibit a power law-like distribution. The histograms for kxx, kyy, and kxy are also fitted by the solid lines with the least square method (Figure 3).The fitted lines of the histograms ofkxx, kyy, and kxy for all realizations are gathered on the same frame (Figure 4). The statistical properties of equivalent permeability also change with the fracture geometric parameters. For the discrete fracture model with lmin=4 m and N=50, the modes of kxx and kyy are less than 1×10−13 m2 and kxy is distributed roughly between −1×10−13 and 1×10−13 m2. The shapes for kxx and kyy tend to a power law shape which are similar to that of the upscaled equivalent permeability for a three-dimensional poorly connected discrete fracture model [43]. This is primarily due to the fact that a lot of grid blocks are dominated by the matrix permeability where the fractures inside do not penetrate the grid blocks or below the percolation threshold. When lmin or N increase, the shapes for kxx and kyy change: the modes migrate from a low value to a high value and the range of kxy expands symmetrically. It should be noted that the shapes of the histograms for kxx and kyy tend to lognormal-like distributions whereas they keep as a normal-like distribution for kxy during the increase of lmin and N. The Gaussian distribution or log-normal distribution for equivalent permeability is commonly assumed in equivalent fracture models (e.g., [15]) and is supported by field measurement data (e.g., [59]). The results in Figure 4 support this assumption, which are based on upscaling permeability for fractured porous rocks. In particular, the fitting curves in Figure 4 further show how the shapes of the histograms of equivalent permeability migrated from power law like to lognormal-like distributions with the change of fracture geometries. The results highlight the importance of fracture geometries on equivalent permeability distributions.Figure 4
Fitted curves for the histograms of equivalent permeability regarding varied fracture geometric parameterslmin and N for all realizations. The blue, red, and green colors denote kxx, kyy, and kxy, respectively.


## 3.2. Correlation between Equivalent Permeability Distribution and Fracture Network Geometry

For analyzing the correlations between the equivalent permeability and the fracture network geometry, the dimensionless equivalent permeabilityk’ of the grid block and the dimensionless fracture density ρ for the discrete fracture model are defined [60], respectively, as follows:
(4)k′=kkm,(5)ρ=1A∑iNli22,where k denotes the components of the equivalent permeability, li denotes the fracture length of the i-th fracture, and A is the area of the discrete fracture model which is 400 m2 in this model. For each discrete fracture model, the mean and standard deviation of the dimensionless equivalent permeability components k′xx¯, k′yy¯, and k′xy¯, and σk′xx, σk′yy, and σk′xy as well as the dimensionless fracture density ρ are calculated. The statistical properties of equivalent permeability tensor components for all 160 discrete fracture models are plotted in Figure 5.Figure 5
Plots of the mean and standard deviation of the dimensionless equivalent permeability tensor components with the dimensionless fracture densityρ for varied fracture geometric parameters.
(a)(b)(c)(d)(e)(f)It is shown that a strong correlation exists between the diagonal components of equivalent permeability tensor and the fracture network geometry:k′xx¯ and k′yy¯ increase with the dimensionless fracture density. The standard deviation for k′xx and k′yy also increases with ρ which is reasonable as few grid blocks may be interconnected by several fractures. The correlation between k′xx¯ and k′yy¯ and ρ can be built by the fitted line:
(6)k′diag¯=a∙ρ+b,where ρis in the range of 2-14, k′diag¯ is the diagonal components of the equivalent permeability tensor which is in the range of 50-350, a is the slope about 20, and b is the constant about 30 here. Leung and Zimmerman [60] investigated the equivalent permeability of a single grid block, and the results show that the equivalent permeability increases linearly with the fracture density. Similar results are also observed for the mean equivalent permeability of grid blocks, as shown in Figure 5. The slope coefficient of a in this study is comparable for those in a single grid block. Although the main scope of the paper concentrates on the influence of fracture network geometries on equivalent permeability distribution, the factor of varied fracture aperture, which makes the model complex, needs to be considered in further studies.Forkxy, the mean k′xy¯ shows a small fluctuation around zero with the changes of ρ. Furthermore, for kxy, its range varies when the fracture network properties change (Figure 4). When ρincreases, the ranges for k′xy¯ and σk′xy also increase. This is mainly because for a grid block with an inclined fracture with azimuth θ, the absolute value of kxy increases with kxx and kyy, which is given by [61] the following:
(7)kθ=kxxkxykyxkyy=kx∗cos2θ+ky∗sin2θkx∗−ky∗cosθ∙sinθkx∗−ky∗cosθ∙sinθkx∗sin2θ+ky∗cos2θ,where kx∗ and ky∗ are the components of the equivalent permeability tensor in x- and y-axes when the fracture azimuth is 0°.


## 4. Discussion and Conclusions

The results and their implications are mainly twofold: (1) the shapes of histograms for equivalent permeability tensors varied with fracture network properties, which helps generate stochastic permeability fields at the continuum scale for fractures porous rocks, and (2) correlating fracture network properties with the equivalent permeability tensors quantitatively; i.e., it would be possible to infer the mean of equivalent permeability tensors from the fracture density. The connectivity, which can be defined as the intersection between fractures or the percolation parameter related to fracture lengths in a volume [34, 37], is an important property for discrete fractures. Thus, two key parameters of the fracture network are considered: the fracture length and the fracture number, which have a great effect on hydraulic property or connectivity of fracture networks. Although generated synthetically, the fracture network geometry is comparable to those of realistic models. The domain size (20m×20m) is a compromise of the availability of fracture geometric data at the field scale and the computation effort. In this study, the fixed domain size is comparable with related studies (e.g., [60]) and coordinates with the fracture geometries. Specifically, the areal fracture intensity P21, which denotes the length of fracture traces per unit area [62], varies from 0.79 to 4.03 m/m2 in these discrete fracture models and falls within a reasonable range of values from the measurement data in the Äspö Hard Rock laboratory [63]. Regarding the power-law distribution of fracture length, ten realizations for a set of fracture geometric parameters seem a bit few to get a representative sample. Based on the results shown in Figure 4 that all ten realizations under the same fracture geometric parameters yield similar results and the histograms vary obviously with the fracture geometric parameters, these mostly indicated that the realizations yield relatively stable statistical results. Accordingly, the domain size and the number of realizations fulfill the requirement of the study. Whereas the current model only considers the fracture network geometries, when adding other factors (e.g., varied aperture), which will make complex models, the multiple domain size would be considered [42] for yielding comprehensive and statistically stable results.The spatial and statistical distributions of equivalent permeability are influenced by fracture network geometry as well as the matrix permeability (e.g., for the discrete fracture model withlmin=4 m and N=50). For the equivalent permeability tensor components, their spatial distributions differ from each other for a given discrete fracture model, which reflects the high anisotropy of the equivalent permeability in fractured porous rocks. Such differences emphasize the importance of considering a full tensor form of equivalent permeability for modeling flow in naturally fractured porous rocks.The fracture network geometries are broadly regarded as critical for understanding hydraulic properties (e.g., [37]). Considering the main scope of this study is to investigate how fracture network geometries (i.e., the connectivity) affect equivalent permeability distributions and devoid of other relevant factors; the aperture is assumed as constant, which indicates that each fracture contributes equally to the fluid flow in fracture networks with an identical permeability, even though the constant aperture is a common assumption for modeling flow in fractured reservoirs or related studies [37, 42, 64]. It should be noted that the varied apertures seem more realistic for the naturally fractured rocks, which is another key property of naturally fractured rocks [45]. They can be characterized by statistical models [34, 65] and are correlated with fracture length [66, 67] and with rock mechanical stress state [18, 68, 69]. The constant aperture can lead to an increase of the equivalent permeability of up to a factor of six as compared to aperture heterogeneities [70]. When the fracture aperture is correlated with fracture trace length, the equivalent permeability of the model can be controlled by fewer numbers of large fractures with large apertures [71]. Such aperture models open the way to study equivalent permeability distributions for naturally fractured rocks and should be addressed in further work.Another limitation of this research is that the current fracture model is analyzed in two dimensions, which corresponds to layered models and is a simplification for naturally fractured porous rocks based on outcrop data [50, 55, 69]. In the two-dimensional models, the fractures are assumed to extend vertically, whereas when fractures are not vertical and their dip angles scatter randomly, the fracture connectivity and equivalent permeability may vary from two-dimensional to three-dimensional models [39, 72]. So it is not very straightforward for correlating results from two dimensions to three dimensions, which depends on specific geological settings and fracture geometric properties. Although when fracture densities are high, three-dimensional equivalent permeability can be predicted with proper corrections based on two-dimensional results [39]. The methodology applied in this study does not have any limitations on three-dimensional naturally fractured porous rocks. For the upscaling procedure, the flow rate is calculated based on three boundaries of the grid block for the two-dimensional discrete fracture model, as shown in Equation (2). In contrast, for three-dimensional models, the flow rated is calculated based on five boundaries, which may have outflows [43]. The analysis of three-dimensional fractured porous rocks has several limitations: (1) the availability of multiscale three-dimensional fracture geometry data, (2) corrections of sampling bias [73], (3) lack of benchmarking models for naturally fractured rocks, and (4) efficient computational platforms for simulating flow in three-dimensional discrete fracture models.To conclude, the equivalent permeability distributions of two-dimensional stochastic discrete fracture models are studied. The fracture length follows a power law distribution, and the equivalent permeability was upscaled from a discrete fracture model to an equivalent fracture model by using the multiple boundary method. The spatial and statistical distributions of equivalent permeability depend on the fracture network geometry as well as the matrix permeability. The spatial distributions ofkxx, kyy, and kxy differ from each other which is primarily dominated by the fracture geometry. The histograms of equivalent permeability tensor components are also different. Initially, the diagonal components, kxx and kyy, tend to a power law-like distribution. The off-diagonal component, kxy, tends to a normal-like distribution and is symmetrically distributed around zero. When the fracture length and fracture number increase, the shapes of the histograms of kxx and kyy change to a lognormal-like distribution and the ranges of kxy expand. The mean of the diagonal equivalent permeability tensor components kxx and kyy increases linearly with the fracture density with the coefficient about 20. The standard deviation of kxx and kyy also increases with the fracture density. Nevertheless, the mean of kxy distribution keeps as zero whereas its standard deviation increased.



---

*Source: 6751349-2020-12-02.xml*