Permeability tensor of three‐dimensional fractured porous rock and a comparison to trace map predictions

The reduction from three‐ to two‐dimensional analysis of the permeability of a fractured rock mass introduces errors in both the magnitude and direction of principal permeabilities. This error is numerically quantified for porous rock by comparing the equivalent permeability of three‐dimensional fracture networks with the values computed on arbitrarily extracted planar trace maps. A method to compute the full permeability tensor of three‐dimensional discrete fracture and matrix models is described. The method is based on the element‐wise averaging of pressure and flux, obtained from a finite element solution to the Laplace problem, and is validated against analytical expressions for periodic anisotropic porous media. For isotropic networks of power law size‐distributed fractures with length‐correlated aperture, two‐dimensional cut planes are shown to underestimate the magnitude of permeability by up to 3 orders of magnitude near the percolation threshold, approaching an average factor of deviation of 3 with increasing fracture density. At low‐fracture densities, percolation may occur in three dimensions but not in any of the two‐dimensional cut planes. Anisotropy of the equivalent permeability tensor varies accordingly and is more pronounced in two‐dimensional extractions. These results confirm that two‐dimensional analysis cannot be directly used as an approximation of three‐dimensional equivalent permeability. However, an alternative expression of the excluded area relates trace map fracture density to an equivalent three‐dimensional fracture density, yielding comparable minimum and maximum permeability. This formulation can be used to approximate three‐dimensional flow properties in cases where only two‐dimensional analysis is available.


Introduction
The hydraulic response of fractured rock to pressure gradients is of vital interest to various Earth science and engineering disciplines, as the majority of host rocks to hydrocarbon reservoirs and nuclear waste disposal sites are heavily fractured [e.g., van Golf-Racht, 1982;Wu et al., 1999]. Explicit representation of fractures in numerical models used for such analysis is of growing popularity, yet due to limitations of computational efficiency, homogenization of small scale features is still required. With respect to fractured rock, equivalent or block permeability [e.g., Renard and de Marsily, 1997;Amaziane and Hontans, 2002] quantifies the ensemble flow response of fractures and rock matrix. In engineering practice, hydraulic analysis of fractured rock is often reduced to a two-dimensional problem, owing to computational limitations, especially if multiphysics effects are incorporated in the model [Nick et al., 2011;Rutqvist et al., 2013]. Three-dimensional flow properties are then deduced from two-dimensional calculations [Pouya and Fouché, 2009;Reeves et al., 2013]. This paper will focus on quantifying the difference between two-and three-dimensional permeability calculations, by means of direct computation of the full equivalent permeability tensor of three-dimensional networks, and extracted planar trace maps.
Overcoming these limitations, Koudina et al. [1998] led the investigation of effective permeability in three-dimensional discrete fracture models based on unstructured grids. Surface representation of fractures and a finite volume discretization of steady state single phase flow equations was used to compute the effective permeability. Bogdanov et al. [2003] further extended the approach by accounting for the permeability of the rock matrix, introducing the use of discrete fracture and matrix (DFM) models to the investigation of effective permeability of fractured rock. This line of analysis assumes representative, spatially periodic fracture network models. Investigating geologically motivated fracture patterns and employing a control volume finite element discretization of the governing equations, Matthäi and Belayneh [2004] derived permeability along the axes of cuboid DFM models, aligned with the dominant fracture orientation, wherein the authors stressed the uncertainty in the directionality of the effective permeability given the assumption that the block axes are aligned with the directions of principal permeability.
A pressure gradient and flux averaging technique is presented for three-dimensional finite element analysis of the steady state flow. It is applied to unstructured grids of DFM models to find the full equivalent, or block, permeability tensor, with no a priori assumptions regarding the orientation of the eigenvectors. This approach is based on developments by Durlofsky [1991], and subsequent applications restricted to structured grids [Wen et al., 2003], and two-dimensional DFM models [Azizmohammadi and Paul, 2012]. The need for periodic boundary conditions to obtain a symmetric tensor is relaxed through the employment of border regions. Validation problems demonstrate that even with asymmetric permeameter constraints on flow experiments, a representative tensor is accurately obtained. Advantages of the present implementation are numerous: full use can be made of flexible meshing techniques and geological realism of representation of heterogeneities; fractures can be represented as surfaces, relaxing the need for local mesh refinement; no assumption on the orientation of principal permeabilities are made; varying permeability in both rock matrix and individual fractures is permitted; and computational efficiency and simplicity is achieved, since the analysis is based on finite element solution of a Laplace equation and subsequent mesh traversal only. Finally, an alternative density formulation is presented to relate the permeability of trace maps to that of three-dimensional models.

Permeability of Fractured Rock
Fluid flow through fractured rock as a response to a hydraulic head gradient will vary heterogeneously throughout the domain. Within the fractures, flux will vary locally depending on orientation, conductivity, and connectivity of the individual joints. Flux in the rock matrix depends on the matrix permeability and whether or not the fracture network allows for pressure gradients to exist across individual blocks [Matthäi and Belayneh, 2004]. Overall, a certain volumetric flow will result as an ensemble fracture-matrix response. The proportionality factor between the block's pressure gradient and this total flow is referred to as block, or equivalent, permeability [e.g., Renard and de Marsily, 1997;Amaziane and Hontans, 2002]. The distinction to effective permeability is based on observations that a representative elementary volume cannot 10.1002/2014JB011027 always be defined in fractured rock due to the fractal nature of joint geometry and the resulting lack of homogenization scale [e.g., Bonnet et al., 2001]. Thus, the permeability of fractured rock is not an intrinsic property but varies strongly with scale and sample. In the context of this study permeability, k, is a symmetric second-rank tensor [e.g., Durlofsky, 1991] that refers to the equivalent permeability of a fractured rock mass at a specific scale of investigation. It follows that where k ij are the components relating the flux and pressure gradients along the x, y, z axes in a three-dimensional Cartesian system. Principal permeabilities are defined in direction by the eigenvectors, and in magnitude by the corresponding eigenvalues, of k. These eigenvalues are referred to here as the minimum, medium, and maximum equivalent permeability, i.e., k equ,min , k equ,med , and k equ,max , respectively. On occasion, e.g., for concise notation in figures, the subscript equ will be omitted. Anisotropy, , is defined as the ratio of the maximum to minimum principal equivalent permeability:

Dimensionless Density
For fractured porous media, the equivalent permeability can be normalized with respect to the matrix permeability k m , where the latter is assumed to be isotropic and homogeneous. It follows that the dimensionless equivalent permeability k ′ is defined as and will herein further be distinguished as k ′ 2-D and k ′ 3-D to discriminate results in two and three dimensions, respectively. Percolation theory for fractured media makes extensive use of dimensionless density to quantify the degree of fracture connectivity, independent of absolute fracture and domain size [Adler and Thovert, 1999]. For two-and three-dimensional media, where ′ is the dimensionless density, is the absolute density, and A ex , V ex are the excluded area and volume, respectively. Absolute density is defined as the ratio of the number of fracture centers, N f , within the domain of area A, or volume V, for two-and three-dimensional networks, respectively, The concept of excluded volume was introduced to the field of percolation theory by Balberg et al. [1984]. V ex defines the minimum volume that two object centers must not reside within, in order for the objects themselves not to intersect, and thus relates ′ to the degree of connectivity of the domain.
In fractured porous media, dimensionless density based on excluded volume is used to identify the onset of fracture percolation, i.e., percolation threshold density ′ c and the evolution of permeability [e.g., Mourzenko et al., 2005;Bogdanov et al., 2007]. Given a distribution of uniform randomly oriented fracture traces of varying length l on a common plane, an average excluded area is defined as [Balberg et al., 1984;Adler and Thovert, 1999] For identical assumptions regarding orientation, disc-shaped fractures of varying radius R in threedimensional space yield an average excluded volume of [Mourzenko et al., 2005] where the average of equation (7) is computed for all combinations of fracture pairs of R 1 , R 2 within the volume of interest. Without a priori assumptions on the fracture size distribution, the dimensionless density for networks is then defined by using the averages ⟨A ex ⟩ and ⟨V ex ⟩ in place of A ex and V ex in equations (4a) and (4b), respectively.
LANG ET AL.

Ellipse and Ellipsoid Visualization of Equivalent Permeability
To visualize orientation and degree of anisotropy of permeability tensors, ellipses and ellipsoids can be used in two and three dimensions, respectively. Since equivalent permeabilities presented in this study span several orders of magnitude, normalization is performed with respect to local k max , such that the geometric aspect ratio resembles k max ∕k min . This normalization applies to all presented plots and implies that ellipses and ellipsoids may be only used to compare orientation and degree of anisotropy of equivalent permeability, but not absolute magnitudes. Figure 1 illustrates isotropic tensors in which ellipses reduce to circles ( Figure 1b) and ellipsoids to spheres ( Figure 1d). Furthermore, matrix elements are not displayed in presented figures, even though the models used in this analysis consist of watertight finite-element meshes, with triangular elements discretizing surfaces, or fractures, and tetrahedral elements discretizing volumetric parts, or rock matrix. Instead, rock matrix will be shown in transparent gray, with no notion of discretization in space.

Methodology
To arrive at full equivalent permeability tensors, a series of incompressible, single-phase filtration problems are solved over a porous medium domain using the finite element method. For each simulation, the resulting scalar pressure field is postprocessed to yield a pressure gradient vector field. Element-wise constant fluxes are computed using local element pressure gradients, and permeability is computed by applying Darcy's law. Contributions to the equivalent permeability tensor components are a function of element volume weighted averaging of pressure gradients and fluxes. The resulting overconstrained problem is solved using a least squares approximation. The asymmetric effect of permeameter-type flow constraints is relaxed by excluding a boundary adjacent peripheral region from the tensor approximation.

Governing Equations and Numerical Solution
The Laplace problem is solved to model incompressible, single-phase flow in porous media, in which conservation of mass reduces to conservation of volume, and volume flux can in turn be expressed using Darcy's law. It follows that where the solution of pressure, p, is approximated at nodes of unstructured finite element grids in three dimensions, employing standard Galerkin method. The viscosity, , is assumed to be constant; for local variation of the intrinsic and isotropic, permeability, k, is permitted. After solution of the pressure field, element-wise constant barycentric velocity is arrived at by applying Darcy's law, to find which results in a vector field of element piecewise constant velocities, u e , that reflect local element pressure gradients p e and permeabilities, k e . The representation of fractures in the mesh is a discrete lower dimensional one. In three dimensions, fractures are modeled using surface elements in between volumetric elements. Two-dimensional models, correspondingly, use line elements in between surface elements. Flux in fracture elements relates to the local pressure gradient through the aperture squared divided by 12, a relationship that results from the assumption of laminar flow between parallel plates obeying a parabolic velocity profile, commonly referred to as the cubic law [e.g., Witherspoon et al., 1980;Zimmerman et al., 1991].

Equivalent Permeability Tensor: 2-D
Darcy's scalar representation of the relationship between flux and pressure gradient can be expressed in vector form for anisotropic media as where u i is the filtration velocity in the i direction, k ij are the components of the permeability tensor, and p x i is the pressure gradient along x i . Note that i, j represent x, y, z in a Cartesian three-dimensional coordinate system. Volume-averaged values for oriented flux, ⟨u j ⟩, and pressure gradient, ⟨ p x i ⟩, are defined as [Durlofsky, 1991] for all elements e within a domain of the volume V t . Note that fractures in the form of lower dimensional elements require a local thickness attribute for volume integration, and the hydraulic aperture is used for this purpose herein. The matrix form in two-dimensional analysis for two flow experiments along the x and y axis, indicated by I and II, respectively [e.g., Wen et al., 2003;Azizmohammadi and Paul, 2012], is where here and in subsequent equations fluid viscosity is assumed to be unity.
LANG ET AL.
©2014. The Authors. (13) matrix is extended to three dimensions and validated for isotropic and anisotropic periodic media. In 3-D, flux and pressure gradient averages given by equations (11) and (12) are computed for directions aligned with the three global coordinate axes, indicated by superscripts I, II, and III. The permeability tensor is hence overdetermined and assumed to be symmetric and positive definite [e.g., Pouya and Fouché, 2009]. Using symmetry conditions as constraints, it follows that the linear problem to be solved of form Ax = b is given by

Equivalent Permeability Tensor: 3-D Equation
where k ij is the equivalent permeability component, with i, j = x, y, z and ⟨u x,y,z ⟩ I,II,III are the averaged flux components for each experiment. In the (12 × 9) matrix on the left-hand side of equation (14), the averaged pressure gradient submatrix ⟨ P j ⟩ i is defined as and the submatrices that enforce symmetry conditions on the permeability tensor, I j , are This matrix problem constitutes an overdetermined system of 12 equations in 9 unknowns, which also must satisfy the constraint k xy = k yx , k xz = k zx , and k yz = k zy . It follows that to arrive at a solution for the tensor components which permits a least squares minimization [Sanford and Dally, 1979]. The obtained tensor components are hence the result of an approximation, which requires the solution to be checked for symmetry. If necessary, k ij = k ji , i ≠ j is enforced by defining the k ij tensor to be the symmetric part of the computed k ij tensor, i.e., The presented matrix problems are based on Darcy's law in tensor form, which in physical terms implies that the diagonal permeability terms relate flux to pressure gradient along the same direction, and the off-diagonal terms relate flux to the perpendicularly oriented pressure gradient. This applies to both the local (element) and the global (fracture matrix) scale, where fluxes and pressure gradients are discrete and averaged, respectively.

Boundary Conditions
The tensor approximation presented in equation (14) requires that three numerical flow experiments are to be run for a three-dimensional problem. The original approach by Durlofsky [1991] works on structured meshes and constrains opposite boundaries to honor periodicity of fluxes. For unstructured grids of discrete fracture networks, however, such conditions are not practicable. Instead, asymmetric permeameter boundary conditions are used in the presented approach. In each of the flow simulations, the pressures on a pair of opposing boundary surfaces are set to yield a prescribed macroscopic pressure difference Δp. Bounding box surfaces oriented parallel to this macroscopic flow direction are assigned to no-flow conditions, i.e., a zero pressure-gradient normal to their plane. As a consequence, calculated equivalent permeability will be a function of the applied boundary conditions [e.g., Renard and de Marsily, 1997;Farmer, 2002] and their asymmetry to some extent. Close to the boundary, element pressure gradients will reflect the prescribed physical flow constraints rather than the effects of local permeability and geometry. The resulting flow response is illustrated in Figure 2, where boundary-adjacent flux vectors honor no-flow boundary conditions-their contribution to the averaging problems in equations (13) and (14) distorts the expected directionality of fluid flow.

Sampling Domain
To retrieve a representative symmetric tensor from flux and pressure fields that, at least near the periphery, are physically influenced by asymmetric constraints, averages in equations (13) and (14) are computed only over a restricted subvolume of the flow region away from the borders. This approach is referred to here as undersampling, since the pressure diffusion problem is solved over the entire model, but equivalent permeability is calculated by sampling the fluxes and pressure gradients in a subdomain only. Other authors refer to similar approaches as "flow jacket" or "oversampling" and find that this procedure reduces boundary effects and better reflects the connectivity of the sampling domain [e.g., Holden and Lia, 1992;Wu et al., 2002]. The sampling subdomain will be chosen such that only those elements are included that keep a minimum distance to all model boundaries; see Figure 3.

Implementation
The finite element analysis is based on an academic, multi-institutional numerical framework CSMP++ for flow [Matthäi et al., 2007], transport [Nick et al., 2011], and geomechanics [Paluszny and Zimmerman, 2011] developed in C++. Ensuing systems of simultaneous linear equations are solved for employing the Frauenhofer Systems Algebraric Multigrid solver [Stüben, 2001]. All model geometries used for validation purposes in this study were prepared using a commercial Non-Uniform Rational Basis Spline 3-D modeling software. All discretizations of geometries to tetrahedral-triangular finite element meshes were performed by a commercial Octree mesher in a noninteractive manner.

Validation
Two-and three-dimensional periodic media models featuring lower dimensional (fractures) and unidimensional (rock matrix) permeability heterogeneities are used to validate this procedure. The sampling subregion is chosen such that it reflects model periodicity. Matrix permeability, k m , is 1.0×10 −15 m 2 , whereas the high permeable features are of k m,h = k m × 10. For all models, fractures are assumed to have a uniform aperture of a f = 0.001 m corresponding to a permeability of k f = k m × 8.33 × 10 7 m 2 according to the cubic law. The relative error of the numerical results with respect to the analytical solution is defined as where is the relative error, k a is the analytical solution for the equivalent permeability, and k n is the numerical equivalent permeability. Periodic models of layered media are used, for which analytical solutions are obtainable by applying weighted and harmonic averages to arrive at an effective permeability tensor in which the eigenvectors are parallel or perpendicular to the layering. The permeability tensor in an arbitrarily oriented coordinate system can be obtained via matrix rotation. Periodic models in this study are homogeneous beyond their scale; hence, equivalent permeability can be considered as effective permeability. Therefore, in this section on model validation, both terms are used interchangeable.

Sugar Cube Models
In a sugar cube model, the effective permeability is isotropic, and the diagonal terms are equal in magnitude to the thickness-weighted mean of the fracture and matrix permeabilities, k avg,w , [e.g., Snow, 1969] where h cube is the cuboid edge length of 10 m and n f = 6 is the number of equidistant fractures normal to each axis (Figure 1). This orientation of fracture planes yields effective permeability eigenvectors that can be taken to be coaxial to the Cartesian coordinate axes, in which case all off-diagonal terms in the equivalent permeability tensor are zero, so that Figure 1 illustrates (Figure 1a) two-and (Figure 1c) three-dimensional sugar cubes. Numerical results yield small relative errors of O(10 −5 ) in two dimensions and O(10 −3 ) in three dimensions, as presented in Table 1.

Anisotropic Periodic Media
Anisotropic periodic media with nonzero off-diagonal terms in the equivalent permeability tensor are used to validate accurate approximation of off-diagonal values in the permeability tensor. An analytical solution is obtained by applying harmonic, k avg,h , and arithmetic, k avg,w , weighted averages to arrive at a permeability LANG ET AL.
©2014. The Authors. tensor, k eff,coax , for a periodic system in which layering of heterogeneities is normal to the y axis.
Calculation of k avg,w is analogous to the procedure used for isotropic models, and k avg,h , defined as consistent with the above notation. Multiplication by a rotation matrix, R, brings the effective permeability tensor to its final form, k eff , where = ∕4 for a system in which layers are rotated by 45 • with respect to the x-y plane (see Figure 4). This approach applies to all periodic verification models presented here and translates to two dimensions by removal of the third row and column from their corresponding submatrices. Figure 4 illustrates periodic matrix only models in (Figure 4a) two and (Figure 4c) three dimensions. Two parallel layers of identical thickness are rotated by 45 • with respect to the horizontal and exhibit a permeability contrast of a factor of 10. Aspect ratios of resulting ellipses and ellipsoids reflect this contrast and are aligned with the orientation of layering.
DFM validation models use lines and surfaces with an aperture attribute to represent fractures, as illustrated in Figures 5a and 5c, respectively.
The permeability contrast between fractures and matrix is O(10 7 ). Therefore, the minimum principal permeability vanishes when visualizing the tensor. In two dimensions the ellipse reduces to a line (Figure 5b), and in three dimensions the ellipsoid reduces to a disc (Figure 5d).
Numerical results are in good agreement with the analytical solution of equation (28) for both matrix-only and fracture-and-matrix models. Relative errors are of O(10 −3 ) and below for all verification models. Part of this error is attributed to the fact that for lower dimensional representation of fractures, no additional degrees of freedom are introduced. This results in a shortcoming in accounting for flow within the fracture perpendicular to its orientation. This brings about an O(10 −4 ) error in matching the harmonic mean, which reflects that fractures contribute to a minor extent to overall conductivity in their normal direction.
LANG ET AL.
©2014. The Authors. Furthermore, fractures may only act additively in terms of conductivity. In geological terms this implies that the effect of sealing fractures cannot be accounted for in the current state of implementation, unlike in the finite-volume approach presented by Bogdanov et al. [2003]. Table 2 lists errors in permeability tensor components and eigenvalues for matrix-only and DFM models in two and three dimensions.

Experimental Setup
A series of cube-shaped DFM models is created using stochastically generated fracture patterns of increasing fracture density. Hydraulic aperture, a f , is sublinearly length correlated [Olson, 2003]. For each of the three-dimensional models, a series of randomly oriented two-dimensional cut planes is extracted. The sampling domain for all models is chosen in a conservative manner such that a distance of 10% of the model edge length is maintained from the nearest model boundaries.

Discrete Fracture and Matrix Models
A set of eleven three-dimensional, stochastic discrete fracture and matrix (DFM) models is generated. Bounding boxes are of cubic shape with an edge length, L, of 100 m. Figure 6 illustrates a selection of three cubes for a (Figure 6a) low, (Figure 6b) medium, and (Figure 6c) high density. Matrix is shown in transparent gray without notion of spatial finite element discretization. Our DFM models idealize fractures as planar discs. In agreement with field observations [e.g., Odling, 1997, Bonnet et al., 2001, the distribution of radii follows a power law where P(r)dr is the probability of observing fractures of radii [r, r + dr] and C is defined through normalization [e.g., Mourzenko et al., 2005] for r ∈ [r m , r M ], with subscripts m and M denoting the minimum and maximum radius. A single exponent of n = −1.5 and r m = L∕10, r M = L∕2 is used for all stochastic realizations. Orientation is uniformly random distributed, with independent Mersenne twister pseudorandom generators [Matsumoto and Nishimura, 1998] used for the fracture normal vector components and the variates used to generate fracture radii. Generated fracture models have between 25 and 900 fractures and range in dimensionless density ′ 3-D from 0.513 to 21.256 or 0.013 to 0.448 in terms of m 2 per m 3 (P32); see Table 3.
Matrix permeability is uniform over the entire model and in magnitude of k m = 1.0 × 10 −15 m 2 . Fracture permeability is obtained from the cubic law using a sublinear length-to-aperture correlation [Olson, 2003].
It follows where a f 1 is the hydraulic aperture of a fracture with 1 m in radius and set to 65 μm [Leung and Zimmerman, 2012]. Obtained apertures are assigned uniformly over the surface of each fracture, our framework, however, allows for arbitrary variations given the element-wise constant discretization of permeability. Figure 6. (a) Low, (b) medium, (c) high fracture density DFM models and a random cut-plane of each. Ellipses and ellipsoids illustrate model equivalent permeability tensor. The orientation of each rectangular cut-plane is random and independent of the global Cartesian coordinate axes. The cut-plane tensors are shown in a local coordinate system ( , ).

Cut Planes
A set of 75 rectangular cut planes is extracted from each three-dimensional DFM model in a random manner. For each coordinate axis, 25 cut planes are defined by two random variables each. These random variables quantify the position of intersection with two opposing boundaries. Figure 7a shows a single cut plane extracted from the original cube DFM model (Figure 7b) where several selected cut planes are highlighted. For the purpose of computing the equivalent permeability tensor, these planes are rotated internally onto the xy plane through Gram-Schmidt orthogonalization. Each node coordinate is found as a linear combination of the plane-spanning vectors. The obtained coefficients are used to transform from the local to the xy plane using the unit vectors e x and e y . This rotation does not change the geometric properties of the plane but only aids to make it conform with the two-dimensional convention. Due to the random nature of plane orientations, the resulting two-dimensional models are of nonuniform aspect ratio and are plotted within their own original two-dimensional coordinate system ( , ). The two end-members in terms of height-to-width ratio can be identified as unity, for cut planes parallel to boundaries and √ 2 for cut planes spanning between two opposing edges of the cuboid. The hydraulic aperture, and hence the permeability, of each fracture trace is equal to that of the three-dimensional disc it intersects. Note that as a consequence the radius of the three-dimensional fracture determines the hydraulic aperture of its traces in  (30). The error in the computed cut plane permeability with respect to that of the three-dimensional network is defined as where k,equ is the error, k equ,3-D and k equ,2-D is the equivalent permeability obtained for the three-dimensional model and its cut-plane, respectively, and log(k) = log 10 (k).

Results
Trace map permeability is considered an approximation which is somewhat representative of the "true" three-dimensional network permeability. This section describes how trace maps consistently underestimate 3-D source permeability, and proposes a density-based correction derived from the results of comparative 2-D and 3-D numerical simulations. The comparison uses two parameters: the magnitude of the principal values k min and k max in terms of k ′ and anisotropy of the full-permeability tensor, . Both are functions of fracture density; however, the evolution in planar models differs from the one of their cubic source models.

Percolation Threshold Densities
Three-dimensional fracture networks exhibit percolation in the range of 2.2 < ′ 3-D,c < 4.2, which becomes apparent by the increase in the minimum and maximum equivalent permeability by an order of magnitude (see Figure 8). By contrast, extracted trace maps do not exhibit percolating fracture pathways until their source networks double in density, 4.5 < ′ 3-D,c,trace maps < 14. A linear relationship exists between cuboid and cut plane density in terms of ′ 3-D and ′ 2-D , with the mean density in two dimensions consistently lower by a factor of 4 (see Figure 9b). Therefore, permeability is examined as a function of ′ 2-D for each trace map. This direct comparison yields percolation within the range of 1.5 < ′ 2-D,c < 3.4; see Figure 9a.

Permeability Magnitude
Over the entire spectrum of fracture density, none of the permeabilities computed from trace maps reach the equivalent permeability obtained from their three-dimensional source (see Figure 8). Moreover, outliers of at least 1 order of magnitude lower are observed over the entire range. At densities below ′ 3-D,c trace map permeabilities are within an order of magnitude of their three-dimensional counterpart. Between ′ 3-D,c and ′ 3-D,c,trace maps the discrepancy is, in average, 2 to 3 orders of magnitude, and the variance in two-dimensional Figure 7. (a, b) Random cut plane sampling concept. Cut planes are extracted from 100 × 100 × 100 m cube DFM models. Figure 7a illustrates the extracted cut planes, while Figure 7b shows the original 3-D model. results is the largest. Above ′ 3-D,c,trace maps , spreading of trace map permeability reduces to 1 order of magnitude and average difference starts to converge toward a factor of 3. The resulting error in trace map permeability is illustrated in Figure 10. The average error reaches a maximum around ′ 3-D,c whereas the spreading in errors is highest around ′ 3-D,c,trace maps . At densities ′ 3-D > ′ 3-D,c,trace maps the error in permeability magnitude approaches a constant value of a factor of 3, accompanied by a convergence of the relative error toward a constant (see Figure 10). These results are indicative of the existence of a correction factor that allows for the direct comparison of 2-D and 3-D permeabilities (see section 7.3).

Density Formulation for Permeability Equivalence
An alternative expression of the average excluded area ⟨A ex ⟩ yields better agreement between 2-D and 3-D permeability and takes the form 2 ⟨l⟩ 2 , for nonpercolating trace maps ⟨l 2 ⟩, for percolating trace maps (32) The postpercolation expression accumulates ⟨l 2 ⟩, as originally proposed by [Balberg et al., 1984], and drops the geometric factor of 2∕ . This acts to increase the 2-D density of the trace maps and conduces to an improved hydraulic equivalence of the 2-D and 3-D data sets. Figure 11 plots cut plane permeability as  a function of dimensionless density, ′ 2-D , evaluated using equation (32), alongside three-dimensional permeability. For presented trace maps, percolation is identified by a threshold k ′ of 10 3 . Above percolation, the average relative error of maximum cut plane permeability is within 0.25. Figure 12 shows the error distribution, which reflects the tendency of underestimating permeability that results from trace map analysis.

Anisotropy
Single, large fractures dictate the orientation of principal permeability at small densities, ′ 3-D < ′ 3-D,c , and anisotropy is moderate due to the lack of fracture percolation; see Figure 13a. A pronounced peak indicates densities close to the percolation threshold, at which single fracture pathways control directionality of flow. Owing to the random orientation of fractures, anisotropy approaches unity with increasing density, ′ 3-D > ′ 3-D,c where percolating pathways exist in numerous directions, rendering the system isotropic in terms of flow response. A resulting unity value of is observed for ′ 3-D > 11. The evolution of trace map anisotropy is similar in nature but differs in magnitude; see Figure 13b. Moderate values of above unity are observed for ′ 3-D < 5. An increase by 2 orders of magnitude is observed over the range 5 < ′ 3-D < 16 where cut planes start to percolate. At higher cube densities, cut plane flow response converges to an isotropic nature, lagging behind their three-dimensional source networks.
Further investigation of two versus three-dimensional anisotropy is illustrated in Figure 14. For each cube model we compare the xy projection of the permeability ellipsoid against ellipses of xy parallel cut planes through them. This sampling concept is shown for a single model in Figure 14l and results of increasing  density presented in Figures 14a to 14k. Low-density models below percolation density of trace maps (Figures 14a to 14c) show two sets of moderate ellipse aspect ratios, owing to two larger isolated fractures present in the cube model. Aspect ratio of the ellipsoid projection is higher, reflecting percolation in three dimensions. At medium densities around percolation threshold of trace maps ( Figures  14d to 14h), few to many to few highly anisotropic ellipses are present, with a maximum at a density that produces the most single percolating pathways. Ellipsoids already reflect increasingly isotropic percolation in volumetric models. At high fracture intensity (Figures 14i to 14k), both cube and cut plane permeability tensors reflect highly isotropic nature of flow paths.

Discussion
The absence of a clear separation of length scales in most fractured rock engineering problems poses a fundamental issue with respect to their homogenization. Rock joints have been shown to extend to maximum lengths between hundred and thousand meters [Bonnet et al., 2001;Olson, 2003]. This range effectively covers typical distances of externally imposed pressure gradients such as interwell communication and excavations in geological disposal facilities to name two examples. The approach presented herein retrieves the full permeability tensor for blocks of fractured porous rock without the assumption that one is working on a representative volume element of the network and thereby separates itself from previous approaches in effective permeability computation [Koudina et al., 1998;Bogdanov et al., 2003]. It has been shown, by example, that the use of border regions in the subsampling procedure allows representative symmetric tensors to be retrieved, even if the physical flow constraints are asymmetric. This contribution aims to improve on the diagonal-tensor characterization of fractured rock block models [e.g., Gonzalez-Garcia et al., 2000] without requiring an extension to periodic models [Malinouskaya et al., 2014].  The presented application on isotropic fracture networks with planar, disc-shaped joints and strict length dependent transmissivity relationship is, of course, an simplification of a nontrivial problem. Geologically realistic fracture networks show high anisotropy of fracture patterns and individual joint permeability [e.g., Olson et al., 2007] as complex functions of displacement history and diagenetic processes. Recent advances in geomechanical modeling allow for the creation of such models in three dimensions [e.g., Paluszny and Matthai, 2010;Paluszny and Zimmerman, 2011]. Thus, we see a need to transition from analyzing a set of statistically similar, stochastically created models toward a more accurate analysis of single realizations, for which the presented method is well suited. A more accurate procedure to arrive at fractured rock permeability could involve the creation of a library of geomechanically grown fracture models, from which suitable candidates are chosen based on field data.
With the inclusion of more physical processes in the numerical analysis of fracture network comes an increase in computational costs. Thermal, hydraulic, mechanical, and chemical models, therefore, still restrain themselves to two dimensions. Based on observations presented herein, reconstruction of three-dimensional equivalent permeability ellipsoids from randomly chosen cut plane ellipses [Pouya and Fouché, 2009] cannot be accurate in neither orientation nor magnitude. Two-dimensional characterization of three-dimensional fracture network permeability in hydrogeologic applications [e.g., Reeves et al., 2013] likely underestimates block conductivity by a factor of 3 or higher. These numerical findings are consistent with previous observations of percolation theory and stereological analysis of 2-D and 3-D networks [e.g., Berkowitz and Adler, 1998]. The corresponding effects on transport problems are likely to be even more severe, given the large variability in flow velocities within the network of fractures [Matthäi and Belayneh, 2004].
With respect to discrepancies in two-dimensional and three-dimensional treatment of hydraulic analysis, isotropic fracture patterns seem most favorable and the difference in results can be taken as best case scenario apart from trivial geometries in which volumetric models are sheer extrapolation of a two-dimensional fracture map. The latter are represented by vertical, layer-bound fracture networks, e.g., as traced in the Bristol channel [Belayneh et al., 2006]. For isotropic models as presented here, results indicate some form of correction factor to exist that translates two-dimensional to three-dimensional conductivities, at least at high fracture densities, where errors approach a linear function in semilog space (see Figure 10). The existence of such a correction for isotropic networks is further based on two concepts. First, with increasing fracture density, both bulk and trace map models are expected to converge toward Snow's prediction for networks of infinitely long fractures [e.g., Snow, 1969;Mourzenko et al., 2011]. Secondly, concepts of stereology suggest that power law length distributions of disc-shaped fractures translate from two-dimensional cut planes to three-dimensional models by the exponent relationship n 3-D = n 2-D + 1 [Marrett, 1996;Piggott, 1997]. An alternative formulation of the dimensionless density allows to find an equivalent relationship between 2-D and 3-D permeability. This relationship could prove useful in the framework of multiphysics simulation, where computational limitations still constrain dimensionality of fracture network models in many cases.

Conclusions
Hydraulic conductivity of fracture networks is an ongoing field of research, and two-dimensional analysis is popular in multiphysics models combining geomechanical, chemical, and temperature effects. A quantification of how the computed fracture network permeability differs as a consequence of this reduction from volumetric to planar models remained elusive so far. This paper presents a method to compute the full equivalent permeability tensor of three-dimensional fracture-matrix models with explicitly represented fractures. The approach is validated for anisotropic models containing volume elements exclusively and for models which contain fractures discretized by surface triangles. The technique only requires the solution to a Laplace equation and involves the subsequent traversing and averaging of local fluxes and pressure gradients over the unstructured finite element mesh. No assumptions regarding the orientation of principal permeabilities are made, allowing for accurate quantification of block permeability with arbitrary anisotropy. The method was subsequently applied to compare the measured effective permeability of volumetric domains and their seemingly equivalent 2-D counterparts. Results show that fracture connectivity in three dimensions is not accurately captured by planar extractions. It follows that values predicted by 2-D trace maps for anisotropy and principal equivalent permeability are not representative of original 3-D data sets. More specifically, fracture percolation occurs at higher densities in two-dimensional cut planes as compared to their volumetric sources. This leads to the well-known density-dependent, but consistent, underestimation of permeability. In fact, maximum and minimum principal permeabilities are underestimated by up to 3 orders of magnitude near the percolation threshold and, in average, by a factor of 3 for higher densities. An approximation of three-dimensional flow behavior based on two-dimensional analysis is hence not directly applicable. Above their percolation threshold, however, trace maps provide an order of magnitude accurate estimation of volumetric permeability. An alternative definition of the excluded area allows the hydraulic response of trace maps to be corrected to reflect the original behavior of the three-dimensional fracture-matrix system.