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

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 di ﬀ er from each other. For the histograms of the equivalent permeability, the shapes of k xx and k yy change from a power law-like distribution to a lognormal-like distribution when the fracture length and the number of fractures increase. For the o ﬀ diagonal component k xy , 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.


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][4][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][11][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][19][20][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 threedimensional 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][35][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 twodimensional 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, twodimensional 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 threedimensional 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] calcu-lated the equivalent permeability distributions for a wellconnected 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.

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]: where l is the fracture length, nðlÞdl 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 a 20 m × 20 m domain based on the ADFNE code [47]. The fracture centers are seeded over in the domain.

Geofluids
The fracture length follows power law distribution, and the lower bound l min and the upper bound l max 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 length l min 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 −4 m is assumed to ensure that all fractures influence fluid flow and to get rid of other affecting factors [42,49]. For the domain of 20 m × 20 m, 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 −9 m 2 according to the cubic law [52]. The matrix permeability is assumed as k m = 9:87 × 10 −16 m 2 (1 md).

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]).
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 xand y-axis, q x and q y , are calculated as follows: where v r , v u , and v l are the specific discharges on the right, upper, and lower boundaries; n x and n y are unit vectors along the xand y-axes; and l x and l y are dimensions of the grid block in the xand y-directions. Then, the equivalent permeability components can be calculated based on the matrix form of Darcy's law: where μ is the dynamic viscosity and ∇P x and ∇P y are pressure gradients along the x-axis and y-axis, respectively. In the case of applying linear boundary conditions along the x -axis, ∇P y is equal to zero and k xx and k yx can be calculated. By using the linear boundary condition along the y-axis, ∇ P x equals zero and the other equivalent permeability components k xy and k yy can be calculated in the same way. As the calculated equivalent permeability components k xy and k yx 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 q x or q y 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 l x and l y 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].

Results
For the discrete fracture models with fractures shown in Figure 1, 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 3 Geofluids equivalent permeability tensor consists of four components, k xx , k xy , k yx , and k yy , for two-dimensional models as shown in Equation (3). Herein, the equivalent permeability tensor is assumed to be symmetric, i.e., k xy = k yx . The statistical distributions of the equivalent permeability and the correlation between the equivalent permeability and the fracture geometric parameter are analyzed in the following.

Statistical Distributions of the Equivalent Permeability.
The Cartesian grid of the equivalent fracture models and the upscaled equivalent permeability field for k xx , k xy , and k yy and for l min = 4 m and N = 50 are plotted in Figure 3. The diagonal components of equivalent permeability tensor k xx and k yy are of the same order of magnitude and less than 3 × 10 −13 m 2 . However, the range of k xy lies between −2 × 10 −13 and 2 × 10 −13 m 2 . This is because the diagonal components k xx and k yy should always be positive according to the physical meaning of the permeability in Darcy's law, whereas the off-diagonal term k xy 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 k xy is smaller than those of k xx and k yy 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, k xx is bigger. This is also similar for k yy . It should be noted that although the range of k xx and k yy is close, their spatial distributions are different (Figure 3) which is primarily dominated by fracture orientations.
The histograms of equivalent permeability for the fracture networks are plotted in Figure 3. The ranges of k xx , k yy , and k xy are divided equally into 10 bins and are gathered on the same frame. It shows that the histograms of k xx , k yy , and k xy exhibit different shapes. For k xy , they tend to be symmetrically distributed with a middle value around 0 m 2 . For k xx and k yy , they are all positive and exhibit a power lawlike distribution. The histograms for k xx , k yy , and k xy are also fitted by the solid lines with the least square method (Figure 3).
The fitted lines of the histograms of k xx , k yy , and k xy 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 l min = 4 m and N = 50, the modes of k xx and k yy are less than 1 × 10 −13 m 2 and k xy is distributed roughly between −1 × 10 −13 and 1 × 10 −13 m 2 . The shapes for k xx and k yy 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

Geofluids
l min or N increase, the shapes for k xx and k yy change: the modes migrate from a low value to a high value and the range of k xy expands symmetrically. It should be noted that the shapes of the histograms for k xx and k yy tend to lognormallike distributions whereas they keep as a normal-like distribution for k xy during the increase of l min 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.

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 permeabil-ity k ' of the grid block and the dimensionless fracture density ρ for the discrete fracture model are defined [60], respectively, as follows: where k denotes the components of the equivalent permeability, l i denotes the fracture length of the i-th fracture, and A is the area of the discrete fracture model which is 400 m 2 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.

Geofluids
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: 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. For k xy , the mean k ′ xy shows a small fluctuation around zero with the changes of ρ. Furthermore, for k xy , 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 k xy increases with k xx and k yy , which is given by [61] the following: where k * x and k * y are the components of the equivalent permeability tensor in xand y-axes when the fracture azimuth is 0°.

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 8 Geofluids (20 m × 20 m) 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 P 21 , which denotes the length of fracture traces per unit area [62], varies from 0.79 to 4.03 m/m 2 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 with l min = 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, threedimensional 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 threedimensional 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 of k xx , k yy , and k xy 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, k xx and k yy , tend to a power law-like distribution. The off-diagonal component, k xy , 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 k xx and k yy change to a lognormal-like distribution and the ranges of k xy expand. The mean of the diagonal equivalent permeability tensor components k xx and k yy increases linearly with the fracture density with the coefficient about 20. The standard deviation of k xx and k yy also increases with the fracture density. Nevertheless, the mean of k xy distribution keeps as zero whereas its standard deviation increased.

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