# Permeability and percolation of anisotropic three-dimensional fracture networks

M. Khamforoush, K. Shams, Jean-François Thovert, P.-M. Adler

► To cite this version:

M. Khamforoush, K. Shams, Jean-François Thovert, P.-M. Adler. Permeability and percolation of anisotropic three-dimensional fracture networks. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics [2001-2015], 2008, 77 (5), pp.056307-056317. 10.1103/PhysRevE.77.056307. hal-00414154

HAL Id: hal-00414154

https://hal.science/hal-00414154v1

Submitted on 16 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Permeability and percolation of anisotropic three-dimensional fracture networks

M. Khamforoush,¹ K. Shams,²,* J.-F. Thovert,³ and P. M. Adler⁴

¹Faculty of Engineering, University of Kurdistan, Sanadaj, Iran

²Department of Chemical Engineering, Isfahan University of Technology, Isfahan 84156, Iran

³LCD, SP2MI, Bôite Postale 30179, 86962 Futuroscope Cedex, France

⁴UPMC-Sisyphe, 4 Place Jussieu, 75252 Paris Cedex 05, France

(Received 27 May 2007; revised manuscript received 27 March 2008; published 12 May 2008)

The percolation properties and permeability of a group of anisotropic three-dimensional fracture networks are studied numerically. Finite-size scaling is used to extrapolate the percolation thresholds of infinite networks in three spatial directions, i.e., X, Y, and Z directions. The influence of the angular dispersion parameter of fracture orientations on percolation thresholds is analyzed. In this analysis, we considered a family of fractures in a three-dimensional space that are oriented around the Z axis based on the Fisher distribution. We revealed that increased anisotropy leads to decreased percolation thresholds in both X and Y directions, and in these two directions percolation thresholds in anisotropic networks demonstrate a declining trend as anisotropy goes up. However, in the Z direction the trend is the opposite. The fracture networks are triangulated via an advancing front technique and the macroscopic permeability of the networks is determined by solving the two-dimensional Darcy equation in each fracture. We found that the macroscopic permeability in the X and Y directions is higher than the associated permeability of isotropic fracture networks, and this property for anisotropic networks in the Z direction is lower compared with that of the isotropic case. Furthermore, as the anisotropy of networks increases the differences become more remarkable.

DOI: 10.1103/PhysRevE.77.056307

PACS number(s): 47.55.—t, 05.40.—a

## I. INTRODUCTION

Fracture in different disciplines has different meanings. In geology and hydrogeology, a single fracture is commonly considered as a subplanar discontinuity in a rock, soils, or geological formations formed mostly by mechanical stresses. A subplanar discontinuity is composed of two solid surfaces which surround a three-dimensional interstitial space. In the earth's crust the gap between single sub-planar discontinuities of fractures may vary between fractions of a millimeter to several kilometers. At large scales, in rocks and geological formations these fractures usually intersect randomly and form a stochastic structure that is called a fracture network. The macroscopic transport properties of the fracture network are a combination of the individual properties of the fractures and depend on their connectivity and intersecting properties.

Permeability is the most significant property of fractured and porous media, and its estimation and measurement is essential in the treatment of flow problems in such environments. It is commonly defined as the property of a medium that allows fluids to pass through it without changing the structure of the medium or displacement of its parts [1]. Permeability is a key transport parameter with great theoretical and practical significance in membrane science, drug delivery systems, medicine, hydrogeology, and reservoir engineering. For example, in reservoir engineering, accurate knowledge of permeability provides petroleum engineers with a tool for efficient management of the production of oil fields, and design and management of enhanced or improved oil recovery operations. In hydrogeology, this transport property is determined by the fracture networks in fractured geological formations. Various aspects of this issue are discussed by Bear et al. [2], Sahimi [3], and Adler and Thovert [4].

The first systematic study of permeability of three-dimensional (3D) fracture networks was initiated by Koudina et al. [5]. The permeability of a 3D network of polygonal

fractures was determined by triangulating the network and solving the 2D Darcy equation in each fracture. Bogdanov et al. [6,7] extended the results reported by Koudina et al. [5] to fractured porous media for one- and two-phase flows and for arbitrary distributions of permeability in the porous matrix and in the fractures. Moreover, the influences of the fracture shape and of their size were studied and rationalized by means of the excluded volume by Mourzenko et al. [8].

Since all these efforts have been restricted to isotropic fracture networks, the main purpose of the present contribution is to extend these results to anisotropic three-dimensional fracture networks which are closer to physical reality. It is well known that the fracture orientation distributions are not uniform in space and that the distribution function of fracture orientations can be characterized by spherical distribution functions such as the Fisher distribution [4,9–13]. An extensive literature search including a recent review article by Berkowitz et al. [14] indicates that a systematic treatment of the permeability problem in anisotropic 3D fracture networks has not been done yet in a systematic way. However, attempts have been made to model specific sites. For example, Sisavath et al. [15] studied a real fracture network which is located in the Baget watershed basin in the southeast of France. They utilized line survey data of the basin to reconstruct several types of anisotropic fracture networks, and studied the percolation properties and permeability of the site.

In the present study, the percolation threshold and permeability of anisotropic three-dimensional fracture networks are systematically investigated. In Sec. II, the geometrical aspects and general features of the anisotropic fracture network model are described. The numerical approaches to determine the excluded volume, percolation threshold, and macroscopic permeability of anisotropic fracture networks are described in Sec. III. Results are gathered in Sec. IV. A systematic account of analytical and numerical results that were ob-
tained in the course of this study is given and discussed. The excluded volume of anisotropic fracture networks is calculated. Then, the influence of anisotropy on the directional percolation threshold is investigated. Finally, the Darcy equation is solved in fracture networks and the macroscopic network permeability is determined and compared with analytical results and isotropic permeability reported in the literature. Furthermore, the effect of anisotropy on permeability is studied. Finally, some concluding remarks are given in Sec. V.

## II. GENERAL

### A. Fundamentals of the fracture network model

The fracture network generator is described in detail by Huseby et al. [16] and Khamforoush and Shams [17]. Briefly, the fracture networks are composed of monodisperse (fixed size) polygons. Each fracture is approximated by a regular plane hexagon inscribed in a disk of radius R.

Since the fracture network is assumed to be statistically homogeneous, it is replaced by a spatially periodic network and only a unit cell of size L³ is generated. L is generally large with respect to the size of the fractures, i.e., L significantly exceeds the size of the fracture (L/R≥4). Hence, L is a homogenization scale over the fractured medium. Therefore, the medium is assumed to be spatially periodic at large scales. It is good to mention that the assumption of a spatially periodic medium is a standard practice in studying infinite media by analytical [18] and/or numerical approaches [16,19]. The general properties of spatially periodic media are given by Adler [19]. The main characteristic of such a medium is that its geometrical properties are invariant under translations, i.e.,

$$ \mathbf{R}_n = n_1 \mathbf{I}_1 + n_2 \mathbf{I}_2 + n_3 \mathbf{I}_3, \quad (1) $$

where n₁, n₂, and n₃ are integers and **I**₁, **I**₂, and **I**₃ are three-vectors that constitute the unit cell to study the fractured system. The whole network is therefore composed of periodic juxtaposition of unit cells in space, translated by **R**ₙ. In this study, we consider only cubic unit cells with |**I**₁|=|**I**₂| = |**I**₃|=L. We assume that the fracture centers are distributed uniformly within the cell. The orientation of each fracture is specified by the unit normal vector of its plane. In order to generate anisotropic fracture networks, the fracture normal vectors are directed in an anisotropic fashion according to the Fisher distribution. The finite size of the unit cells induces problems close to the percolation threshold as will be seen in Sec. III B.

### B. Fisher distribution

Consider three-dimensional fracture networks where each fracture is perpendicular to its unit normal vector **n**. In standard polar coordinates, this unit normal vector is defined by the two angles θ and φ. Therefore, the orientation vector **n** can be considered as a point on the surface of a unit sphere. Since the fracturing process in natural rocks is usually induced by mechanical stresses [20], the fracture orientations are often nonuniformly distributed in space, and usually families of fractures with a preferred orientation along the

stress direction are formed. The fracture orientations of each set are more or less widely spread around a mean direction [4,12,13].

One of the successful models used to describe the real angular dispersion of fractures is the Fisher distribution [4,11–13]. The Fisher distribution is the analog of the Gaussian distribution on the sphere. If (θ₀, φ₀) are the initial polar coordinates, the Fisher distribution can be deduced as follows [21]:

$$ f(\theta, \phi) = \frac{\kappa}{4\pi \sinh \kappa} \sin \theta \exp\{\kappa[\cos \theta_0 \cos \theta + \sin \theta_0 \cos(\phi - \phi_0)]\}, \quad (2) $$

where 0≤θ≤π, 0≤φ≤2π, and the distribution is rotationally symmetric around the initial direction (mean pole). For the particular case where (θ₀, φ₀) are equal to zero, the Fisher distribution reduces to

$$ f(\theta, \phi) = \frac{\kappa}{4\pi \sinh \kappa} \sin \theta \exp(\kappa \cos \theta). \quad (3) $$

Therefore, the normal vectors are oriented around the Z axis. Here, κ is the dispersion parameter about the mean direction. It is called the Fisher distribution parameter; the larger κ, the more the distribution is concentrated around the direction of the mean pole. This behavior is illustrated in Fig. 1 and a typical anisotropic fracture network is represented in Fig. 2. In this network, the θ angles of fractures are allocated according to Eq. (3) and the φ angles of fractures are distributed uniformly in [0, 2π].

However, in order to homogenize the rotational distribution of polygons around the normal vectors of their planes, a rotation angle ψ about the normal vector of each fracture on the polygon plane is defined. The value of this angle is uniformly chosen in the interval [0, 2π]. Such networks with different values of κ are used to determine the percolation threshold and macroscopic permeability of anisotropic systems in this study.

## III. COMPUTATIONAL PROCEDURE

### A. Excluded volume

In all the subsequent calculations, it is assumed that the fractures are hexagonal. It should be mentioned that it is not necessary to consider the fracture shape effect in this study since it has been shown that, where fracture network densities are rendered dimensionless by means of the excluded volume, the macroscopic network properties are independent of fracture shape [5–8,16,22]. The excluded volume Vex of an object in the context of fractured media was defined for the first time by Balberg et al. [23] as the volume into which the center of another similar object should not enter, if one wants to avoid overlap of two objects. For identical plane polygons in isotropic fracture networks, Vex is defined by Adler and Thovert [4] as

$$ V_{\text{ex}} = \frac{1}{2}AP \quad (4) $$
<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Probability density of $\theta$ according to the Fisher distribution for $\kappa=1, 5, 10, 20$, and 50 (a). Distribution of fracture normal vectors over a unit sphere for $\kappa=10$ (b), 20 (c), and 50 (d). Each dot corresponds to a unit normal. The colors are used to distinguish the different hexagons and carry no further information.</figcaption></figure>

<figure><img src="image_3.png" /></figure>

where A and P are the polygon area and perimeter, respectively.

The density $\rho$ of fracture networks is defined as the number of fractures per unit volume $L^3$. Therefore, if $N_{fr}$ hexagons are distributed in a cubic cell of dimension L, the density $\rho$ will be equal to

$$\rho = \frac{N_{fr}}{L^3}. \qquad (5)$$

The dimensionless density $\rho'$ is defined as the number of fractures per excluded volume or equivalently

$$\rho' = \rho V_{ex}. \qquad (6)$$

Huseby et al. [16] have shown that $\rho'$ is exactly equal to the mean number of intersections per fracture. This parameter represents a direct measure of the network connectivity; therefore

$$\rho' = \rho V_{ex} = \langle n_I \rangle. \qquad (7)$$

It should be emphasized that Eq. (4) is valid only for isotropic fracture orientations. For the anisotropic case, an analytical derivation of the excluded volume $V_{ex}(\kappa)$ seems impossible [4], especially for a complex distribution function such as the Fisher distribution. However, it is expected that the

mean number of intersections per fracture is still related to $V_{ex}(\kappa)$ by a relation analogous to (7) because of the definitions of the various quantities involved. For the sake of completeness, let us write this relation as

$$V_{ex}(\kappa) = \frac{\langle n_I \rangle}{\rho}, \qquad (8)$$

where $\langle n_I \rangle$ is the average number of intersections per fracture and $V_{ex}(\kappa)$ is the excluded volume of the anisotropic fracture networks.

This obvious relation is extremely useful since it will be used in order to estimate $V_{ex}(\kappa)$ by means of Monte Carlo calculations. Random anisotropic networks are generated for a given value of $\kappa$; it is a simple matter to determine the total number of intersections in these networks and to derive $\langle n_I \rangle$ and $V_{ex}(\kappa)$. Note that Eq. (8) can also be used for isotropic networks, because it presents the frequency of intersections, and is independent of fracture orientations.

## B. Percolation threshold

To estimate the percolation threshold of the anisotropic fracture network, the classical finite-size scaling technique described by Stauffer and Aharony [24] is used. The percolation threshold is
<figure><img src="image_4.png" /><figcaption>FIG. 2. (Color online) Typical illustration of anisotropic 3D fracture networks composed of monodisperse hexagons. The volume of size $L^3$ contains 450 hexagons with $\kappa=20$ and $L=12R$. The colors are used to distinguish the different hexagons and carry no further information.</figcaption></figure>

lation threshold is studied independently along the X, Y, and Z directions in spatially periodic unit cells. The probability of percolation (i.e., the existence of continuous clusters which span the system) is a function of relative cell size $L/R$, of the Fisher parameter $\kappa$, and of the fracture network density $\rho$. However, the critical density of fracture networks, $\rho_{c_A}$, at the percolation threshold depends on the shape of the fractures. When the excluded volume is used, the dimensionless percolation threshold $\rho'_{c_A}$, which is the critical number of fractures per excluded volume, is constant for a wide range of polygons [16]. Therefore, the anisotropic dimensionless density is used instead of the simple density.

The relative cell size $L/R$ is changed from 4 to 12, and the dimensionless density is gradually increased. To investigate the effect of the Fisher constant on percolation thresholds the value of $\kappa$ is varied as $\kappa=\{2, 5, 10, 20, 50\}$. For given values of $L$, $\rho'_A$, and $\kappa$, the probability of existence of a percolating cluster, $\Pi(L, \kappa, \rho'_A)$, is determined along each direction from 500 realizations of the system. This value is a trade-off between precision and CPU time recommended in [16]. The critical density $\rho'_{c_{A,L}}$, or percolation threshold of a finite system of size $L$, can be defined as the value for which $\Pi(L, \kappa, \rho'_A)$ is equal to 0.5. It is determined by fitting the data for $\Pi(L, \kappa, \rho'_A)$ with a two-parameter error function of the following form [24,25]:

$$ \Pi(L, \kappa, \rho'_A) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\rho'_A} \frac{1}{\Delta L} \exp\left(-\frac{(\xi - \rho'_{c_{A,L}})^2}{2(\Delta L)^2}\right) d\xi, \quad (9) $$

where $\rho'_{c_{A,L}}$ and the width of the transition region $\Delta L$ are adjustable parameters. Since the finite-size percolation threshold $\rho'_{c_{A,L}}$ is size dependent, it is only in the limit of infinite cell size that one can define the real percolation threshold $\rho'_{c_{A,\infty}}$, whose value for an infinite cell is denoted by the subscript $\infty$.

The value of the infinite percolation threshold $\rho'_{c_{A,\infty}}$ can be extrapolated from the data obtained for finite cells by two scaling relations [26,27]:

$$ \rho'_{c_{A,L}} - \rho'_{c_{A,\infty}} \propto \left(\frac{1}{L}\right)^{1/\nu}, \quad (10) $$

$$ \Delta L \propto \left(\frac{1}{L}\right)^{1/\nu}. \quad (11) $$

Practically speaking, in order to calculate the real percolation threshold, $\rho'_{c_{A,L}}$ is plotted as a function of $\Delta L$, a linear fit is done, and $\rho'_{c_{A,\infty}}$ is determined when $\Delta L$ approaches to zero.

### C. The flow problem

Since nonpercolating systems have zero permeability, flow calculations are limited to densities above the percolation threshold. The fractures are assumed to be embedded in an impermeable solid matrix. On a local scale within a typical fracture of aperture $b$, where the Reynolds number is low, the Newtonian fluid flow is governed by the Stokes equations. According to the classical Poiseuille law, the typical permeability $\sigma$ of a fracture with aperture $b$ is of the order

$$ \sigma = \frac{b^3}{12}, \quad (12) $$

however, $b$ is assumed to be much smaller than the typical lateral extent $2R$ of the fracture, and flow on scales between $b$ and $2R$ is governed by the Darcy equation
$$ \mathbf{q} = -\frac{1}{\mu} \boldsymbol{\sigma} \cdot \overline{\nabla p}, \qquad (13) $$

where **q** is the locally averaged flow rate per unit width (dimensions $L^2 T^{-1}$), $\mu$ is the fluid viscosity, $\boldsymbol{\sigma}$ (dimensions $L^3$) is the fracture permeability tensor, and $\overline{\nabla} p$ is the pressure gradient. In this paper, $\boldsymbol{\sigma}$ is taken to be constant and isotropic over each fracture. The dependence of $\boldsymbol{\sigma}$ upon the fracture characteristics was investigated by Mourzenko et al. [28] and O'Brien et al. [29]. Therefore, the tensor $\boldsymbol{\sigma}$ reduces to $\boldsymbol{\sigma}\mathbf{I}$ where $\mathbf{I}$ is the unit tensor and $\boldsymbol{\sigma}$ the scalar fracture permeability.

The mass conservation equation for an incompressible fluid is

$$ \nabla_s \cdot \mathbf{q} = 0, \qquad (14) $$

where $\nabla_s$ is the two-dimensional gradient operator in the mean fracture plane.

The fractured medium is composed of the periodic juxta-position of identical unit cells along the X, Y, and Z directions and a macroscopic pressure gradient $\overline{\nabla} p$ is applied as boundary condition upon the unbounded medium of fracture networks. Fluid flow is described by Eqs. (13) and (14), together with periodic conditions for $\mathbf{v}$, $\mathbf{q}$, and $\overline{\nabla} p$. The macroscopic pressure gradient $\overline{\nabla} p$ can be expressed as follows:

$$ \overline{\nabla} p = \frac{1}{v} \int_{\partial v} p \, ds, \qquad (15) $$

where $v$ is the volume of the cubic cell.

The local equations are solved according to the methodology described in [5]. When the local fields are known, the macroscopic quantities of interest such as the seepage velocity $\overline{v}$ can be evaluated as

$$ \overline{v} = \frac{1}{v} \int_{v_f} \mathbf{v} \, dv = \frac{1}{v} \int_{s_f} \mathbf{q} \, ds, \qquad (16) $$

where $v_f$ is the interstitial volume of the fractures and $s_f$ is their projection on their mean planes. This flux is related to the macroscopic pressure gradient by a Darcy law [4] valid for the fracture network,

$$ \overline{v} = -\frac{1}{\mu} \mathbf{K} \cdot \overline{\nabla} p, \qquad (17) $$

where **K** is the macroscopic permeability tensor (dimensions $L^2$), to be derived from Eqs. (16) and (17), when Eqs. (13) and (14) are solved. It must be noted that **K** is a symmetric and positive definite tensor when averaged over all realizations of networks,

$$ \mathbf{K} = \mathbf{K}\mathbf{I}. \qquad (18) $$

In this study, $b$ is defined as the characteristic length for the fracture aperture; $R$ is the typical fracture extent, and $p_0$ is a reference pressure. These three parameters can be used to recast the equations in a dimensionless form. We define

$$ \sigma_0 = \frac{b^3}{12}, \quad q_0 = \frac{\sigma_0 p_0}{\mu R}, \quad U_0 = \frac{\sigma_0 p_0}{\mu R^2}, \quad K_0 = \frac{\sigma_0}{R}. \qquad (19) $$

The dimensionless quantities and gradient operator (with primes) are defined as follows:

$$ \begin{align} p &= p_0 p', \quad \nabla &= \frac{1}{R} \nabla', \notag \\ \mathbf{v} &= U_0 \mathbf{v}', \quad \mathbf{q} = q_0 \mathbf{q}', \tag{20} \\ \boldsymbol{\sigma} &= \boldsymbol{\sigma} \boldsymbol{\sigma}', \quad \mathbf{K} = K_0 \mathbf{K}'. \notag \end{align} $$

This dimensionless formulation is used in the following developments. The dimensionless permeability tensor $\mathbf{K}'$ can be simplified because of the statistical symmetry of the generated networks. Since $\theta_0$ is equal to zero, the two axes X and Y play similar roles; therefore,

$$ \mathbf{K}' = \begin{bmatrix} K'_{xx} & 0 & 0 \\ 0 & K'_{yy} & 0 \\ 0 & 0 & K'_{zz} \end{bmatrix} = \begin{bmatrix} K'_{\perp} & 0 & 0 \\ 0 & K'_{\perp} & 0 \\ 0 & 0 & K'_{\parallel} \end{bmatrix}, \qquad (21) $$

where $K'_{\perp}$ is the dimensionless scalar permeability in the X and Y directions and $K'_{\parallel}$ is the dimensionless scalar permeability in the Z direction.

The numerical approach and code used to solve the flow problem are the same as implemented by Koudina et al. [5]. First, the fracture network is triangulated by the advancing front technique [30,31]. This technique is an appropriate method for unstructured triangulation. The mesh is characterized by the prescribed maximum edge length $\delta_M$ which is set equal to $R/4$ [5]. The flow equations are discretized by the classical finite volume method which consists of integrating the equation over elementary volumes. The value of the pressure $p'$ is determined at each point of the triangular mesh via a flux balance condition by the finite volume technique. Finally, the linear algebraic equations of discretized transport equation are solved to specify the unknown pressure at each point.

## IV. RESULTS AND DISCUSSION

In this section, a systematic account of analytical and numerical results that were obtained in the course of this study is given and discussed. First, the excluded volumes and the percolation thresholds are presented. Then, the influence of $\kappa$ (angular dispersion parameter) on analytical and numerical values of the permeability is introduced.

### A. Excluded volume

Excluded volumes of anisotropic fracture networks for various values of the Fisher parameter, as well as the excluded volume of an isotropic fracture network, are calculated for various densities, ranging from slightly below the percolation threshold up to the connected networks. In each case, $\langle n_i \rangle$ is determined and the excluded volume is calculated by Eq. (8). The corresponding excluded volumes are averaged over 1000 realizations for each fracture network.
Since the fracture networks are generated in a spatially periodic medium, the average number of intersections is not affected by the cell size.

The results obtained for the isotropic case and anisotropic cases for different values of the Fisher parameter are given in Table I. As expected, the excluded volumes are independent of density and the statistical fluctuations of the averages influence only the second decimal digit of the quantity. Consequently, we have used the average of these values as the final excluded volume in the subsequent calculations.

The average value of the calculated excluded volume of the isotropic network ($\kappa=0$) over various densities is equal to 7.7926, and the theoretical value of the excluded volume for hexagons by a straightforward application of Eq. (4) is equal to $9\sqrt{3}/2=7.7942$. The comparison between these two values shows that the numerical excluded volume of isotropic fracture networks is approximately equal to the theoretical value. This shows that Eq. (8) is reliable for computing the excluded volume.

The final results for anisotropic excluded volumes for different values of the Fisher parameter are illustrated in Fig. 3, where $V_{ex_A}$ is a decreasing function of $\kappa$. Increasing the anisotropy $\kappa$, as shown in Fig. 3, induces the reduction of the excluded volume of the network. This behavior can be explained based on the directional configuration of fractures. For $\kappa=0$, fracture networks are isotropic and the fracture normal vectors are uniformly distributed over the whole surface of the unit sphere. Conceptually, the average number of intersections should be a maximum in this case. However, for strictly positive values of $\kappa$, the fracture normal vectors are distributed in a narrow region about the mean pole and as the value of $\kappa$ increases this region gradually shrinks (see Fig. 1). Hence the fractures are oriented in a nearly parallel situation, and the intersections between the fractures are drastically decreased. This implies that the excluded volume is reduced.

## B. Percolation threshold

The infinite percolation thresholds (average number of connections per object) along the X, Y, and Z directions are systematically reported in Table II and represented in Fig. 4 for various values of the Fisher parameter $\kappa$. The following comments are in order to explain these data.

For the isotropic case, $\kappa=0$, the percolation thresholds in the X, Y, and Z directions are nearly equal. This is to be expected, since, in a fracture network with elements uniformly oriented in 3D space, the chances of finding the percolating clusters should be the same in all directions. The average of the isotropic percolation thresholds in the three directions is equal to 2.3019 which is approximately equal to the reported value 2.3 of Huseby *et al.* [16].

For nonisotropic cases, $\kappa>0$, when $\kappa$ increases, the values of the percolation thresholds in the X, Y and Z directions are changed. For a constant value of $\kappa$, the percolation thresholds along the X and Y directions are nearly identical. However, the percolation threshold in the Z direction is higher than the corresponding values in the X and Y directions. In other words, because of the special configuration of the fracture normal vectors, the average number of necessary

<table><caption>TABLE I. The excluded volume $V_{ex_A}(\kappa)$ of isotropic and anisotropic fracture networks for various densities versus the Fisher distribution parameter $\kappa \cdot \sigma_{V_{ex_A}}$ is the standard deviation of $V_{ex_A}$ for 1000 realizations.</caption><thead><tr><th>κ</th><th>ρ=N<sub>fr</sub>/L<sup>3</sup></th><th>V<sub>exA</sub>(κ)</th><th>σ<sub>V<sub>exA</sub></sub></th></tr></thead><tbody><tr><td rowspan="10">0.0</td><td>0.2441</td><td>7.7925</td><td>0.0717</td></tr><tr><td>0.2563</td><td>7.7931</td><td>0.0655</td></tr><tr><td>0.2686</td><td>7.7891</td><td>0.0569</td></tr><tr><td>0.2808</td><td>7.7929</td><td>0.0515</td></tr><tr><td>0.2930</td><td>7.7976</td><td>0.0481</td></tr><tr><td>0.3052</td><td>7.7954</td><td>0.0453</td></tr><tr><td>0.3174</td><td>7.7933</td><td>0.0427</td></tr><tr><td>0.3296</td><td>7.7891</td><td>0.0377</td></tr><tr><td>0.3418</td><td>7.7865</td><td>0.0356</td></tr><tr><td>0.3540</td><td>7.7968</td><td>0.0351</td></tr><tr><td rowspan="10">2.0</td><td>0.2441</td><td>7.6141</td><td>0.0672</td></tr><tr><td>0.2563</td><td>7.6064</td><td>0.0567</td></tr><tr><td>0.2686</td><td>7.6091</td><td>0.0556</td></tr><tr><td>0.2808</td><td>7.6143</td><td>0.0540</td></tr><tr><td>0.2930</td><td>7.6085</td><td>0.0451</td></tr><tr><td>0.3052</td><td>7.6025</td><td>0.0431</td></tr><tr><td>0.3174</td><td>7.6058</td><td>0.0379</td></tr><tr><td>0.3296</td><td>7.6071</td><td>0.0380</td></tr><tr><td>0.3418</td><td>7.6076</td><td>0.0370</td></tr><tr><td>0.3540</td><td>7.6153</td><td>0.0325</td></tr><tr><td rowspan="10">5.0</td><td>0.2441</td><td>6.4608</td><td>0.0583</td></tr><tr><td>0.2563</td><td>6.4616</td><td>0.0523</td></tr><tr><td>0.2686</td><td>6.4537</td><td>0.0513</td></tr><tr><td>0.2808</td><td>6.4675</td><td>0.0476</td></tr><tr><td>0.2930</td><td>6.4621</td><td>0.0419</td></tr><tr><td>0.3052</td><td>6.4585</td><td>0.0410</td></tr><tr><td>0.3174</td><td>6.4670</td><td>0.0364</td></tr><tr><td>0.3296</td><td>6.4517</td><td>0.0327</td></tr><tr><td>0.3418</td><td>6.4611</td><td>0.0318</td></tr><tr><td>0.3540</td><td>6.4639</td><td>0.0299</td></tr><tr><td rowspan="10">10.0</td><td>0.2441</td><td>5.1216</td><td>0.0458</td></tr><tr><td>0.2563</td><td>5.1327</td><td>0.0441</td></tr><tr><td>0.2686</td><td>5.1182</td><td>0.0413</td></tr><tr><td>0.2808</td><td>5.1234</td><td>0.0372</td></tr><tr><td>0.2930</td><td>5.1221</td><td>0.0329</td></tr><tr><td>0.3052</td><td>5.1232</td><td>0.0327</td></tr><tr><td>0.3174</td><td>5.1289</td><td>0.0288</td></tr><tr><td>0.3296</td><td>5.1256</td><td>0.0264</td></tr><tr><td>0.3418</td><td>5.1301</td><td>0.0273</td></tr><tr><td>0.3540</td><td>5.1265</td><td>0.0243</td></tr><tr><td rowspan="10">20.0</td><td>0.2441</td><td>3.7672</td><td>0.0346</td></tr><tr><td>0.2563</td><td>3.7685</td><td>0.0310</td></tr><tr><td>0.2686</td><td>3.7575</td><td>0.0299</td></tr><tr><td>0.2808</td><td>3.7608</td><td>0.0271</td></tr><tr><td>0.2930</td><td>3.7629</td><td>0.0234</td></tr><tr><td>0.3052</td><td>3.7521</td><td>0.0223</td></tr><tr><td>0.3174</td><td>3.7625</td><td>0.0217</td></tr><tr><td>0.3296</td><td>3.7634</td><td>0.0198</td></tr><tr><td>0.3418</td><td>3.7692</td><td>0.0185</td></tr><tr><td>0.3540</td><td>3.7613</td><td>0.018</td></tr></tbody></table>
<table><caption>TABLE I. (Continued.)</caption><thead><tr><th>κ</th><th>ρ=N<sub>fr</sub>/L<sup>3</sup></th><th>V<sub>exA</sub>(κ)</th><th>σ<sub>V<sub>exA</sub></sub></th></tr></thead><tbody><tr><td rowspan="8">50.0</td><td>0.2441</td><td>2.4521</td><td>0.0225</td></tr><tr><td>0.2563</td><td>2.4578</td><td>0.0201</td></tr><tr><td>0.2686</td><td>2.4516</td><td>0.0187</td></tr><tr><td>0.2808</td><td>2.4503</td><td>0.0168</td></tr><tr><td>0.2930</td><td>2.4632</td><td>0.0153</td></tr><tr><td>0.3052</td><td>2.4524</td><td>0.0140</td></tr><tr><td>0.3174</td><td>2.4529</td><td>0.0143</td></tr><tr><td>0.3296</td><td>2.4553</td><td>0.0128</td></tr><tr><td>0.3418</td><td>2.4603</td><td>0.0118</td></tr><tr><td>0.3540</td><td>2.4588</td><td>0.0116</td></tr></tbody></table>

connections per object for percolation in the Z direction compared with the X and Y directions occurs at higher densities. This means that the chances of percolation for a given density in the Z direction are lower than the chances of percolation along the X and Y directions. The reason should be the nonisotropic orientations of fractures about the Z axis. That is, the arrangement of fracture orientations as shown in the typical illustration of fracture networks in Fig. 2 has a format that hinders the percolation in the Z direction. This is because all the fractures have been oriented about a preferred direction along the Z axis and it is intuitive that the intersections of fractures with each other in the Z direction are prevented when compared to the X and Y directions. Furthermore, the Fisher distribution is rotationally symmetric about its mean pole; therefore the percolation thresholds in the X and Y directions must be the same. This is confirmed by the results listed in Table II. As the anisotropy of networks increases, the angular dispersion of fracture normal vectors is restricted to a narrow region around their pole. For larger values of $κ$ this region becomes smaller, so the fractures tend to a relatively parallel situation. The direct result of this transformation is the reduction of fracture intersections as is shown in Fig. 3.

In addition, Balberg [32] has established the following universal criterion for percolation of convex objects in 3D systems whether they are isotropic or not:

<figure><img src="image_5.png" /><figcaption>FIG. 3. (Color online) Excluded volume of anisotropic 3D fracture networks $V_{exA}$, versus $κ$ (▲). The straight line is the isotropic excluded volume (■).</figcaption></figure>

<table><caption>TABLE II. The asymptotic values of percolation thresholds for isotropic and anisotropic fracture networks versus the Fisher distribution parameter in the X, Y, and Z directions.</caption><thead><tr><th>κ</th><th>ρ'<sub>c∞x,A</sub></th><th>ρ'<sub>c∞y,A</sub></th><th>ρ'<sub>c∞z,A</sub></th></tr></thead><tbody><tr><td>(Isotropic case)</td><td>2.3083</td><td>2.2939</td><td>2.3036</td></tr><tr><td>2</td><td>2.2795</td><td>2.2704</td><td>2.3227</td></tr><tr><td>5</td><td>2.2295</td><td>2.2335</td><td>2.3720</td></tr><tr><td>10</td><td>2.1821</td><td>2.1773</td><td>2.4022</td></tr><tr><td>20</td><td>2.1363</td><td>2.1402</td><td>2.4207</td></tr><tr><td>50</td><td>2.1063</td><td>2.1094</td><td>2.4381</td></tr></tbody></table>

$$0.7 < \rho'_{c} < 2.8. \qquad (22)$$

The previous results agree with the percolation threshold limits reported by [32].

### C. Permeability of anisotropic fracture networks

#### 1. Analytical (Snow) model

Snow [10] considered networks made of infinite planes with an arbitrary orientation distribution. The general expression of the permeability tensor $K_A$ of an anisotropic network made of infinite fractures is given in [4,10] as

$$K_A = \int \int_{\Omega} S(\mathbf{n}) K_S(\mathbf{n}) (\mathbf{I} - \mathbf{n}\mathbf{n}) d^2\mathbf{n}, \qquad (23)$$

where $S(\mathbf{n})$ is the surface area per unit volume of the family of fractures with normal $\mathbf{n}$ and $K_S(\mathbf{n})$ is the surface permeability of these fractures.

Let us restrict this general expression to a network whose orientations obey the Fisher distribution (2) and whose surface permeability is a constant equal to $\sigma$. Then,

$$K_A = \rho A \sigma \int \int_{\Omega} f(\theta, \varphi) (\mathbf{I} - \mathbf{n}\mathbf{n}) \sin \theta \, d\theta \, d\varphi. \qquad (24)$$

When $\mathbf{n}$ is expressed in spherical coordinates, the integral in $\theta$ is evaluated with the trivial change of variable $X=\cos \theta$. We easily obtain

<figure><img src="image_6.png" /><figcaption>FIG. 4. (Color online) Percolation thresholds for infinite anisotropic 3D fracture networks along the X (●), Y (■), and Z directions (▲) against the Fisher distribution parameter $κ$.</figcaption></figure>
$$K_{\perp} = K_n \frac{3}{2 \sinh \kappa} \left[ \sinh \kappa \left( 1 + \frac{1}{\kappa^2} \right) - \frac{1}{\kappa} \cosh \kappa \right], \quad (25)$$

$$K_{\parallel} = K_n \frac{3}{\kappa^2 \sinh \kappa} (\kappa \cosh \kappa - \sinh \kappa), \quad (26)$$

where $K_n$ is the permeability of an isotropic network ($\kappa=0$). For regular hexagons Koudina et al. [5] have deduced that $K_n=2\rho'/9$.

## 2. Numerical results

The permeability of anisotropic networks of regular monodisperse hexagonal fractures is systematically determined in a cubic cell with periodic boundary conditions for a normalized density $\rho'$ ranging from slightly below the percolation threshold up to $\rho'=18$ for various values of the Fisher parameter, i.e., $\kappa=\{2,5,10,20,50\}$. As for the percolation study, it is not necessary to investigate the effect of fracture shape on permeability, because Koudina et al. [5] have shown that, if the density of fracture networks is rendered dimensionless by the excluded volume, the final results for the permeability will be independent of the fracture shape; also they have shown that in the proximity of the percolation threshold the permeability of polygonal fracture networks is slightly affected by the cell size, but for dimensionless densities $\rho'$ greater than 3.5 the cell size effect is vanished. Therefore, the permeability of anisotropic fracture networks in this study is predicted at a relative cell size $L/R=4$. All the results along the three spatial directions for $\rho' \le 7.0$ and $\rho' > 7.0$ are statistical averages over 400 and 100 random realizations of fracture networks, respectively.

As we mentioned earlier, the flow properties of anisotropic fracture networks have not been systematically studied yet, and we could not find a similar systematic study in the literature to compare our results with. Therefore, we have tried to validate our permeability results by comparing them with three sets of data. The first two are reported by Koudina et al. [5] for isotropic monodisperse hexagonal fracture networks and the third are the analytical results which were derived in the previous section for anisotropic fracture networks with infinite fractures. The first set of data reported in [5] is based on the solution of the two-dimensional Darcy equation in each fracture; we have reproduced them in Figs. 5–10 by thin solid lines passing through solid circles. The second and third sets of data are based on analytic expressions that have been obtained based on the Snow model [10] for isotropic and anisotropic networks, respectively. According to Koudina et al.’s discussion [5], the Snow model is valid only in the limit of very dense networks. In subsequent figures we have shown this model for isotropic and anisotropic networks by thick solid lines.

As mentioned before, the permeability tensor is symmetric, and as expected this tensor is anisotropic. In other words, due to the special arrangement of fracture orientations, the value of permeability in the Z direction is different from the associated values in the X and Y directions. This behavior is illustrated for a typical value of the Fisher parameter (e.g., $\kappa=10$) in Figs. 5 and 6 along with numerical and analytical results. Both figures present the same results. The

<figure><img src="image_7.png" /><figcaption>FIG. 5. (Color online) Dimensionless permeability components in X (+), Y (○), and Z directions (□) for anisotropic 3D fracture networks with $\kappa=10$ versus the normalized density $\rho'_A$ in a cubic cell with size $L=4R$ in Cartesian coordinates. The thick solid lines are the Snow model, and isotropic permeability is indicated by (●).</figcaption></figure>

only difference is that the former is a Cartesian plot and the latter is a log-log plot. As anticipated, the dimensionless permeability in the Z direction, $K'_{zz}$, is much smaller than the corresponding values in the X and Y directions, i.e., $K'_{xx}$ and $K'_{yy}$. This is because fractures are oriented about the Z axis and the fracture connections in the Z direction are drastically lower than the fracture connections in the X and Y directions. As shown in Figs. 5 and 6, there is good agreement between the numerical and analytical results.

Moreover, as revealed in Figs. 5 and 6, $K'_{xx}$ and $K'_{yy}$ are larger than the associated permeability of isotropic fracture networks, while $K'_{zz}$ is lower than the isotropic case, in agreement with physical intuition. One can see from Figs. 5 and 6

<figure><img src="image_8.png" /><figcaption>FIG. 6. (Color online) Dimensionless permeability components in X (+), Y (○), and Z directions (□) for anisotropic 3D fracture networks with $\kappa=10$ versus $\rho'_A-\rho'_{c\infty,A}$ in a cubic cell with size $L=4R$ in log-log coordinates. The thick solid lines are the Snow model, and isotropic permeability is indicated by (●).</figcaption></figure>
<figure><img src="image_9.png" /></figure>

FIG. 7. (Color online) Average dimensionless permeability in X and Y directions, $(K'_x+K'_y)/2$, of anisotropic 3D fracture networks for various values of the Fisher distribution parameters $\kappa=2$ (+), 5 (◇), 10 (□), 20 (×), and 50 (Δ), versus the normalized density $\rho'_A$. The results are for $L=4R$. The thick solid lines are the Snow model and the isotropic case is indicated by (●).

that the values of $K'_{xx}$ and $K'_{yy}$ are approximately equal. So we define the arithmetic average permeability as $(K'_{xx}+K'_{yy})/2$; then to demonstrate the effect of anisotropy on permeability it suffices that we consider the variations of $(K'_{xx}+K'_{yy})/2$ and $K'_{zz}$ with the respect to the Fisher parameter separately. This is shown in Figs. 7–10. In these figures the effect of the Fisher dispersion parameter $\kappa$ has been taken into account in the macroscopic permeability.

When the density is rendered dimensionless according to the anisotropic excluded volume $V_{exA}$, from Figs. 7 and 8 one can see that increased $\kappa$ leads to considerable growth of per-

<figure><img src="image_10.png" /></figure>

FIG. 9. (Color online) Dimensionless permeability component in the Z direction, $K'_z$, of anisotropic 3D fracture networks for various values of Fisher distribution parameters, $\kappa=2$ (+), 5 (◇), 10 (□), 20 (×), and 50 (Δ), versus the normalized density $\rho'_A$. The results are for $L=4R$. The thick solid lines are the Snow model and the isotropic case is indicated by (●).

meability along the X and Y directions, especially when compared with isotropic networks, but the opposite trend is observed for the permeability in the Z direction as is illustrated in Figs. 9 and 10. This behavior can be interpreted as follows: for low values of the Fisher parameter, the fracture normal vectors are distributed in a wide region over the surface of the unit sphere; therefore the fracture networks tend to be isotropic and the difference between $(K'_{xx}+K'_{yy})/2$ and $K'_{zz}$ is not considerable. However, for large values of $\kappa$, as intuition suggests, the difference between $(K'_{xx}+K'_{yy})/2$ and $K'_{zz}$ has become more profound. In other words, when the

<figure><img src="image_11.png" /></figure>

FIG. 8. (Color online) Average dimensionless permeability in X and Y directions, $(K'_x+K'_y)/2$, of anisotropic 3D fracture networks for various values of Fisher distribution parameters $\kappa=2$ (+), 5 (◇), 10 (□), 20 (×), and 50 (Δ), versus $\rho'_A-\rho'_{c,\infty,A}$. The results are for $L=4R$. The thick solid lines are the Snow model and the isotropic case is indicated by (●).

<figure><img src="image_12.png" /></figure>

FIG. 10. (Color online) Dimensionless permeability component in the Z direction, $K'_{zz}$, of anisotropic 3D fracture networks for various values of Fisher distribution parameters, $\kappa=2$ (+), 5 (◇), 10 (□), 20 (×), and 50 (Δ), versus $\rho'_A-\rho'_{c,\infty,A}$. The results are for $L=4R$. The thick solid lines are the Snow model and the isotropic case is indicated by (●).
Fisher parameter increases, fracture normal vectors are more and more inclined toward the Z axis. As a result, the fractures tend to a parallel situation along the X and Y directions. Therefore, the average number of connections in the Z direction is decreased and the chances of finding the connected spanning paths are decreased too. Consequently, since there are not enough connections to transmit fluid, the permeability in the Z direction, $K'_{zz}$, is drastically reduced.

In contrast, the results of Table II reveal that the probability of existence of percolating clusters in the X and Y directions is increased; therefore, the number of connected spanning paths is increased too. As a matter of fact, this behavior is confirmed by the results shown on Figs. 7 and 8. Increased permeability due to larger values of $\kappa$ along the X and Y directions can be attributed to the generation of new connected spanning paths which provide open channels for fluid transmission. These findings comply with the analytical results as illustrated in Figs. 7–10.

## V. CONCLUDING REMARKS

In this study a group of anisotropic fracture networks is constructed according to the Fisher distribution. The anisotropy of fracture networks is controlled by the Fisher dispersion parameter. As the Fisher parameter is increased, the dis-

tribution of fracture normal vectors becomes more dense and inclined toward the Z axis. The excluded volume, percolation properties, and permeability of fracture networks are studied systematically for various values of angular dispersion parameter. For the anisotropic case, at low values of the dispersion parameter, fracture networks tend to be isotropic, and the difference between $(K'_{xx}+K'_{yy})/2$ and $K'_{zz}$ is not appreciable. However, for large values of $\kappa$ the difference between $(K'_{xx}+K'_{yy})/2$ and $K'_{zz}$ is increased. It is found that, as the Fisher parameter is increased the probability of existence of percolating paths along the X and Y directions goes up compared with the probability of existence of percolating paths in the Z direction. Furthermore, the permeability of fracture networks in the X and Y directions is increased when normal vectors of fractures are distributed in a smaller region about the mean pole. Also, the Z component of permeability of fracture networks decreases in comparison with the isotropic case. This result is in good agreement with our physical perception.

## ACKNOWLEDGMENTS

M.K. would like to acknowledge the support of University of Kurdistan and of the Ministry of Science, Research, and Technology of Iran that made his visit to Institute de Physique du Globe de Paris possible.

[1] F. K. North, *Petroleum Geology* (Chapman & Hall, London, 1994).

[2] J. Bear, C.-F. Tsang, and G. de Marsily, *Flow and Contaminant Transport in Fractured Rock* (Academic Press, San Diego, 1993).

[3] M. Sahimi, *Flow and Transport in Porous Media and Fractured Rock* (VCH, New York, 1995).

[4] P. M. Adler and J.-F. Thovert, *Fractures and Fracture Networks* (Kluwer Academic, Dordrecht, 1999).

[5] N. Koudina, R. Gonzalez Garcia, J.-F. Thovert, and P. M. Adler, *Phys. Rev. E* **57**, 4466 (1998).

[6] I. I. Bogdanov, V. V. Mourzenko, J.-F. Thovert, and P. M. Adler, *Water Resour. Res.* **39**, 1023 (2003).

[7] I. I. Bogdanov, V. V. Mourzenko, J.-F. Thovert, and P. M. Adler, *Phys. Rev. E* **68**, 026703 (2003).

[8] V. V. Mourzenko, J.-F. Thovert, and P. M. Adler, *Phys. Rev. E* **69**, 066307 (2004).

[9] R. Fisher, *Proc. R. Soc. London, Ser. A* **217**, 295 (1953).

[10] D. T. Snow, *Water Resour. Res.* **5**, 1273 (1969).

[11] W. S. Dershowitz and H. H. Einstein, *Rock Mech. Rock Eng.* **21**, 21 (1988).

[12] S. D. Priest, *Discontinuity Analysis for Rock Engineering* (Chapman & Hall, London, 1993).

[13] J. Kemeny and P. Randy, *Comput. Geosci.* **29**, 65 (2003).

[14] B. Berkowitz, *Adv. Water Resour.* **25**, 861 (2002).

[15] S. Sisavath, V. V. Mourzenko, P. Genthon, J.-F. Thovert, and P. M. Adler, *Geophys. J. Int.* **157**, 917 (2004).

[16] O. Huseby, J.-F. Thovert, and P. M. Adler, *J. Phys. A* **30**, 1415 (1997).

[17] M. Khamforoush and K. Shams (unpublished).

[18] S. Whitaker, *The Method of Volume Averaging* (Kluwer, Dordrecht, 1998).

[19] P. M. Adler, *Porous Media* (Butterworth-Heinemann, Stoneham, MA, 1992).

[20] L. Valentini, D. Perugini, and G. Poli, *Physica A* **377**, 323 (2007).

[21] K. V. Mardia and P. E. Jupp, *Directional Statistics* (John Wiley & Sons, Chichester, U.K., 2000).

[22] O. Huseby, J.-F. Thovert, and P. M. Adler, *Phys. Fluids* **13**, 594 (2001).

[23] I. Balberg, C. H. Anderson, S. Alexander, and N. Wagner, *Phys. Rev. B* **30**, 3933 (1984).

[24] D. Stauffer and A. Aharony, *Introduction to Percolation Theory* (Taylor & Francis, London, 1994).

[25] P. J. Reynolds, H. E. Stanley, and W. Klein, *Phys. Rev. B* **21**, 1223 (1980).

[26] M. E. Fisher, in *Critical Phenomena*, Proceedings of the International School of Physics “Enrico Fermi,” Course LI, Varenna, 1971, edited by M. S. Green (Academic Press, New York, 1971).

[27] E. Charlaix, *J. Phys. A* **19**, L533 (1986).

[28] V. V. Mourzenko, J. F. Thovert, and P. M. Adler, *J. Phys. II* **5**, 465 (1995).

[29] G. S. O’Brien, C. J. Bean, and F. McDermott, *Earth Planet. Sci. Lett.* **213**, 271 (2003).

[30] P. L. George, *Automatic Mesh Generation* (John Wiley & Sons, New York, 1991).

[31] T. J. Chung, *Computational Fluid Dynamics* (Cambridge University Press, Cambridge, U.K., 2003).

[32] I. Balberg, *Phys. Rev. B* **31**, 4053 (1985).